| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 //
>> 27 // $Id: G4Trap.cc,v 1.41 2006/06/29 18:45:41 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-08-01-patch-01 $
>> 29 //
>> 30 // class G4Trap
>> 31 //
26 // Implementation for G4Trap class 32 // Implementation for G4Trap class
27 // 33 //
28 // 21.03.95 P.Kent: Modified for `tolerant' ge << 34 // History:
>> 35 //
>> 36 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal
>> 37 // 26.04.05 V.Grichine: new SurfaceNormal is default
>> 38 // 19.04.05 V.Grichine: bug fixed in G4Trap("name",G4ThreeVector[8] vp)
>> 39 // 12.12.04 V.Grichine: SurfaceNormal with edges/vertices
>> 40 // 15.11.04 V.Grichine: bug fixed in G4Trap("name",G4ThreeVector[8] vp)
>> 41 // 13.12.99 V.Grichine: bug fixed in DistanceToIn(p,v)
>> 42 // 19.11.99 V.Grichine: kUndef was added to Eside enum
>> 43 // 04.06.99 S.Giani: Fixed CalculateExtent in rotated case.
>> 44 // 08.12.97 J.Allison: Added "nominal" constructor and method SetAllParameters.
>> 45 // 01.11.96 V.Grichine: Costructor for Right Angular Wedge from STEP, G4Trd/Para
29 // 09.09.96 V.Grichine: Final modifications be 46 // 09.09.96 V.Grichine: Final modifications before to commit
30 // 08.12.97 J.Allison: Added "nominal" constru << 47 // 21.03.95 P.Kent: Modified for `tolerant' geometry
31 // 28.04.05 V.Grichine: new SurfaceNormal acco << 48 //
32 // 18.04.17 E.Tcherniaev: complete revision, s << 49 ////////////////////////////////////////////////////////////////////////////////////
33 // ------------------------------------------- <<
34 50
35 #include "G4Trap.hh" 51 #include "G4Trap.hh"
36 <<
37 #if !defined(G4GEOM_USE_UTRAP) <<
38 <<
39 #include "globals.hh" 52 #include "globals.hh"
40 #include "G4GeomTools.hh" <<
41 53
42 #include "G4VoxelLimits.hh" 54 #include "G4VoxelLimits.hh"
43 #include "G4AffineTransform.hh" 55 #include "G4AffineTransform.hh"
44 #include "G4BoundingEnvelope.hh" <<
45 56
46 #include "G4VPVParameterisation.hh" 57 #include "G4VPVParameterisation.hh"
47 58
48 #include "G4QuickRand.hh" << 59 #include "Randomize.hh"
49 60
50 #include "G4VGraphicsScene.hh" 61 #include "G4VGraphicsScene.hh"
51 #include "G4Polyhedron.hh" 62 #include "G4Polyhedron.hh"
>> 63 #include "G4NURBS.hh"
>> 64 #include "G4NURBSbox.hh"
52 65
53 using namespace CLHEP; 66 using namespace CLHEP;
54 67
>> 68 ////////////////////////////////////////////////////////////////////////
>> 69 //
>> 70 // Accuracy of coplanarity
>> 71
>> 72 const G4double kCoplanar_Tolerance = 1E-4 ;
>> 73
>> 74 //////////////////////////////////////////////////////////////////////////
>> 75 //
>> 76 // Private enum: Not for external use
>> 77
>> 78 enum Eside {kUndef,ks0,ks1,ks2,ks3,kPZ,kMZ};
>> 79
55 ////////////////////////////////////////////// 80 //////////////////////////////////////////////////////////////////////////
56 // 81 //
57 // Constructor - check and set half-widths as << 82 // Constructor - check and set half-widths as well as angles:
58 // final check of coplanarity 83 // final check of coplanarity
59 84
60 G4Trap::G4Trap( const G4String& pName, 85 G4Trap::G4Trap( const G4String& pName,
61 G4double pDz, 86 G4double pDz,
62 G4double pTheta, G4doubl 87 G4double pTheta, G4double pPhi,
63 G4double pDy1, G4double 88 G4double pDy1, G4double pDx1, G4double pDx2,
64 G4double pAlp1, 89 G4double pAlp1,
65 G4double pDy2, G4double 90 G4double pDy2, G4double pDx3, G4double pDx4,
66 G4double pAlp2 ) << 91 G4double pAlp2)
67 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 92 : G4CSGSolid(pName)
68 { 93 {
69 fDz = pDz; << 94 if ( pDz > 0 && pDy1 > 0 && pDx1 > 0 &&
70 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi << 95 pDx2 > 0 && pDy2 > 0 && pDx3 > 0 && pDx4 > 0 )
71 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi << 96 {
72 << 97 fDz=pDz;
73 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp << 98 fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
74 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp << 99 fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
>> 100
>> 101 fDy1=pDy1;
>> 102 fDx1=pDx1;
>> 103 fDx2=pDx2;
>> 104 fTalpha1=std::tan(pAlp1);
>> 105
>> 106 fDy2=pDy2;
>> 107 fDx3=pDx3;
>> 108 fDx4=pDx4;
>> 109 fTalpha2=std::tan(pAlp2);
75 110
76 CheckParameters(); << 111 MakePlanes();
77 MakePlanes(); << 112 }
>> 113 else
>> 114 {
>> 115 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl
>> 116 << " Invalid dimensions !" << G4endl
>> 117 << " X - "
>> 118 << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
>> 119 << " Y - " << pDy1 << ", " << pDy2 << G4endl
>> 120 << " Z - " << pDz << G4endl;
>> 121 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 122 "Invalid length G4Trap parameters.");
>> 123 }
78 } 124 }
79 125
80 ////////////////////////////////////////////// << 126 ////////////////////////////////////////////////////////////////////////////
81 // 127 //
82 // Constructor - Design of trapezoid based on << 128 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
83 // which are its vertices. Checking of planari << 129 // which are its vertices. Checking of planarity with preparation of
84 // fPlanes[] and than calculation of other mem 130 // fPlanes[] and than calculation of other members
85 131
86 G4Trap::G4Trap( const G4String& pName, 132 G4Trap::G4Trap( const G4String& pName,
87 const G4ThreeVector pt[8] ) 133 const G4ThreeVector pt[8] )
88 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 134 : G4CSGSolid(pName)
89 { 135 {
90 // Start with check of centering - the cente 136 // Start with check of centering - the center of gravity trap line
91 // should cross the origin of frame 137 // should cross the origin of frame
92 // <<
93 if ( pt[0].z() >= 0 <<
94 || pt[0].z() != pt[1].z() <<
95 || pt[0].z() != pt[2].z() <<
96 || pt[0].z() != pt[3].z() <<
97 <<
98 || pt[4].z() <= 0 <<
99 || pt[4].z() != pt[5].z() <<
100 || pt[4].z() != pt[6].z() <<
101 || pt[4].z() != pt[7].z() <<
102 <<
103 || std::fabs( pt[0].z() + pt[4].z() ) <<
104 <<
105 || pt[0].y() != pt[1].y() <<
106 || pt[2].y() != pt[3].y() <<
107 || pt[4].y() != pt[5].y() <<
108 || pt[6].y() != pt[7].y() <<
109 <<
110 || std::fabs(pt[0].y()+pt[2].y()+pt[4] <<
111 || std::fabs(pt[0].x()+pt[1].x()+pt[4] <<
112 pt[2].x()+pt[3].x()+pt[6] <<
113 { <<
114 std::ostringstream message; <<
115 message << "Invalid vertice coordinates fo <<
116 G4Exception("G4Trap::G4Trap()", "GeomSolid <<
117 FatalException, message); <<
118 } <<
119 <<
120 // Set parameters <<
121 // <<
122 fDz = (pt[7]).z(); <<
123 <<
124 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5; <<
125 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5; <<
126 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5; <<
127 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]). <<
128 <<
129 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5; <<
130 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5; <<
131 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5; <<
132 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]). <<
133 138
134 fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx << 139 if ( pt[0].z() < 0
135 fTthetaSphi = ((pt[4]).y()+fDy2)/fDz; << 140 && pt[0].z() == pt[1].z() && pt[0].z() == pt[2].z() && pt[0].z() == pt[3].z()
>> 141 && pt[4].z() > 0
>> 142 && pt[4].z() == pt[5].z() && pt[4].z() == pt[6].z() && pt[4].z() == pt[7].z()
>> 143 && ( pt[0].z() + pt[4].z() ) == 0
>> 144 && pt[0].y() == pt[1].y() && pt[2].y() == pt[3].y()
>> 145 && pt[4].y() == pt[5].y() && pt[6].y() == pt[7].y()
>> 146 && ( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) == 0
>> 147 && ( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() +
>> 148 pt[2].x() + pt[3].x() + pt[6].x() + pt[7].x() ) == 0 )
>> 149 {
>> 150 G4bool good;
>> 151
>> 152 // Bottom side with normal approx. -Y
>> 153
>> 154 good = MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
>> 155
>> 156 if (!good)
>> 157 {
>> 158 DumpInfo();
>> 159 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 160 "Face at ~-Y not planar.");
>> 161 }
>> 162
>> 163 // Top side with normal approx. +Y
>> 164
>> 165 good = MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
>> 166
>> 167 if (!good)
>> 168 {
>> 169 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 170 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 171 "Face at ~+Y not planar.");
>> 172 }
>> 173
>> 174 // Front side with normal approx. -X
>> 175
>> 176 good = MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
>> 177
>> 178 if (!good)
>> 179 {
>> 180 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 181 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 182 "Face at ~-X not planar.");
>> 183 }
>> 184
>> 185 // Back side iwth normal approx. +X
>> 186
>> 187 good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
>> 188 if (!good)
>> 189 {
>> 190 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 191 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 192 "Face at ~+X not planar.");
>> 193 }
>> 194 fDz = (pt[7]).z() ;
>> 195
>> 196 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5;
>> 197 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5;
>> 198 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5;
>> 199 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1;
>> 200
>> 201 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5;
>> 202 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5;
>> 203 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5;
>> 204 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2;
136 205
137 CheckParameters(); << 206 fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz;
138 MakePlanes(pt); << 207 fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;
>> 208 }
>> 209 else
>> 210 {
>> 211 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 212 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 213 "Invalid vertice coordinates.");
>> 214 }
139 } 215 }
140 216
141 ////////////////////////////////////////////// << 217 //////////////////////////////////////////////////////////////////////////////
142 // 218 //
143 // Constructor for Right Angular Wedge from ST 219 // Constructor for Right Angular Wedge from STEP
144 220
145 G4Trap::G4Trap( const G4String& pName, 221 G4Trap::G4Trap( const G4String& pName,
146 G4double pZ, 222 G4double pZ,
147 G4double pY, 223 G4double pY,
148 G4double pX, G4double pL 224 G4double pX, G4double pLTX )
149 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 225 : G4CSGSolid(pName)
>> 226
150 { 227 {
151 fDz = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi << 228 G4bool good;
152 fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLT << 229
153 fDy2 = fDy1; fDx3 = fDx1; fDx4 = fDx2; << 230 if ( pZ>0 && pY>0 && pX>0 && pLTX>0 && pLTX<=pX )
>> 231 {
>> 232 fDz = 0.5*pZ ;
>> 233 fTthetaCphi = 0 ;
>> 234 fTthetaSphi = 0 ;
>> 235
>> 236 fDy1 = 0.5*pY;
>> 237 fDx1 = 0.5*pX ;
>> 238 fDx2 = 0.5*pLTX;
>> 239 fTalpha1 = 0.5*(pLTX - pX)/pY;
>> 240
>> 241 fDy2 = fDy1 ;
>> 242 fDx3 = fDx1;
>> 243 fDx4 = fDx2 ;
>> 244 fTalpha2 = fTalpha1 ;
>> 245
>> 246 G4ThreeVector pt[8] ;
>> 247
>> 248 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
>> 249 -fDz*fTthetaSphi-fDy1,-fDz);
>> 250 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
>> 251 -fDz*fTthetaSphi-fDy1,-fDz);
>> 252 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
>> 253 -fDz*fTthetaSphi+fDy1,-fDz);
>> 254 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
>> 255 -fDz*fTthetaSphi+fDy1,-fDz);
>> 256 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
>> 257 +fDz*fTthetaSphi-fDy2,+fDz);
>> 258 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
>> 259 +fDz*fTthetaSphi-fDy2,+fDz);
>> 260 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
>> 261 +fDz*fTthetaSphi+fDy2,+fDz);
>> 262 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
>> 263 +fDz*fTthetaSphi+fDy2,+fDz);
>> 264
>> 265 // Bottom side with normal approx. -Y
>> 266 //
>> 267 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
>> 268 if (!good)
>> 269 {
>> 270 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 271 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 272 "Face at ~-Y not planar.");
>> 273 }
>> 274
>> 275 // Top side with normal approx. +Y
>> 276 //
>> 277 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
>> 278 if (!good)
>> 279 {
>> 280 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 281 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 282 "Face at ~+Y not planar.");
>> 283 }
154 284
155 CheckParameters(); << 285 // Front side with normal approx. -X
156 MakePlanes(); << 286 //
>> 287 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
>> 288 if (!good)
>> 289 {
>> 290 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 291 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 292 "Face at ~-X not planar.");
>> 293 }
>> 294
>> 295 // Back side iwth normal approx. +X
>> 296 //
>> 297 good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
>> 298 if (!good)
>> 299 {
>> 300 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 301 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 302 "Face at ~+X not planar.");
>> 303 }
>> 304 }
>> 305 else
>> 306 {
>> 307 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 308 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 309 "Invalid length G4Trap parameters.");
>> 310 }
157 } 311 }
158 312
159 ////////////////////////////////////////////// << 313 ///////////////////////////////////////////////////////////////////////////////
160 // 314 //
161 // Constructor for G4Trd 315 // Constructor for G4Trd
162 316
163 G4Trap::G4Trap( const G4String& pName, 317 G4Trap::G4Trap( const G4String& pName,
164 G4double pDx1, G4double 318 G4double pDx1, G4double pDx2,
165 G4double pDy1, G4double 319 G4double pDy1, G4double pDy2,
166 G4double pDz ) 320 G4double pDz )
167 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 321 : G4CSGSolid(pName)
168 { 322 {
169 fDz = pDz; fTthetaCphi = 0; fTthetaSphi = << 323 G4bool good;
170 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalp << 324
171 fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalp << 325 if ( pDz>0 && pDy1>0 && pDx1>0 && pDx2>0 && pDy2>0 )
>> 326 {
>> 327 fDz = pDz;
>> 328 fTthetaCphi = 0 ;
>> 329 fTthetaSphi = 0 ;
>> 330
>> 331 fDy1 = pDy1 ;
>> 332 fDx1 = pDx1 ;
>> 333 fDx2 = pDx1 ;
>> 334 fTalpha1 = 0 ;
>> 335
>> 336 fDy2 = pDy2 ;
>> 337 fDx3 = pDx2 ;
>> 338 fDx4 = pDx2 ;
>> 339 fTalpha2 = 0 ;
>> 340
>> 341 G4ThreeVector pt[8] ;
>> 342
>> 343 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
>> 344 -fDz*fTthetaSphi-fDy1,-fDz);
>> 345 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
>> 346 -fDz*fTthetaSphi-fDy1,-fDz);
>> 347 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
>> 348 -fDz*fTthetaSphi+fDy1,-fDz);
>> 349 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
>> 350 -fDz*fTthetaSphi+fDy1,-fDz);
>> 351 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
>> 352 +fDz*fTthetaSphi-fDy2,+fDz);
>> 353 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
>> 354 +fDz*fTthetaSphi-fDy2,+fDz);
>> 355 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
>> 356 +fDz*fTthetaSphi+fDy2,+fDz);
>> 357 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
>> 358 +fDz*fTthetaSphi+fDy2,+fDz);
>> 359
>> 360 // Bottom side with normal approx. -Y
>> 361 //
>> 362 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
>> 363 if (!good)
>> 364 {
>> 365 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 366 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 367 "Face at ~-Y not planar.");
>> 368 }
>> 369
>> 370 // Top side with normal approx. +Y
>> 371 //
>> 372 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
>> 373 if (!good)
>> 374 {
>> 375 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 376 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 377 "Face at ~+Y not planar.");
>> 378 }
172 379
173 CheckParameters(); << 380 // Front side with normal approx. -X
174 MakePlanes(); << 381 //
>> 382 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
>> 383 if (!good)
>> 384 {
>> 385 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 386 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 387 "Face at ~-X not planar.");
>> 388 }
>> 389
>> 390 // Back side iwth normal approx. +X
>> 391 //
>> 392 good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
>> 393 if (!good)
>> 394 {
>> 395 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 396 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 397 "Face at ~+X not planar.");
>> 398 }
>> 399 }
>> 400 else
>> 401 {
>> 402 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 403 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 404 "Invalid length G4Trap parameters.");
>> 405 }
175 } 406 }
176 407
177 ////////////////////////////////////////////// << 408 ////////////////////////////////////////////////////////////////////////////
178 // 409 //
179 // Constructor for G4Para 410 // Constructor for G4Para
180 411
181 G4Trap::G4Trap( const G4String& pName, 412 G4Trap::G4Trap( const G4String& pName,
182 G4double pDx, G4double p 413 G4double pDx, G4double pDy,
183 G4double pDz, 414 G4double pDz,
184 G4double pAlpha, 415 G4double pAlpha,
185 G4double pTheta, G4doubl << 416 G4double pTheta, G4double pPhi)
186 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 417 : G4CSGSolid(pName)
187 { 418 {
188 fDz = pDz; << 419 G4bool good;
189 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi << 420
190 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi << 421 if ( pDz>0 && pDy>0 && pDx>0 )
>> 422 {
>> 423 fDz = pDz ;
>> 424 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi) ;
>> 425 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi) ;
>> 426
>> 427 fDy1 = pDy ;
>> 428 fDx1 = pDx ;
>> 429 fDx2 = pDx ;
>> 430 fTalpha1 = std::tan(pAlpha) ;
>> 431
>> 432 fDy2 = pDy ;
>> 433 fDx3 = pDx ;
>> 434 fDx4 = pDx ;
>> 435 fTalpha2 = fTalpha1 ;
>> 436
>> 437 G4ThreeVector pt[8] ;
>> 438
>> 439 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
>> 440 -fDz*fTthetaSphi-fDy1,-fDz);
>> 441 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
>> 442 -fDz*fTthetaSphi-fDy1,-fDz);
>> 443 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
>> 444 -fDz*fTthetaSphi+fDy1,-fDz);
>> 445 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
>> 446 -fDz*fTthetaSphi+fDy1,-fDz);
>> 447 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
>> 448 +fDz*fTthetaSphi-fDy2,+fDz);
>> 449 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
>> 450 +fDz*fTthetaSphi-fDy2,+fDz);
>> 451 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
>> 452 +fDz*fTthetaSphi+fDy2,+fDz);
>> 453 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
>> 454 +fDz*fTthetaSphi+fDy2,+fDz);
191 455
192 fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1 << 456 // Bottom side with normal approx. -Y
193 fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2 << 457 //
>> 458 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
>> 459 if (!good)
>> 460 {
>> 461 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 462 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 463 "Face at ~-Y not planar.");
>> 464 }
194 465
195 CheckParameters(); << 466 // Top side with normal approx. +Y
196 MakePlanes(); << 467 //
>> 468 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
>> 469 if (!good)
>> 470 {
>> 471 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 472 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 473 "Face at ~+Y not planar.");
>> 474 }
>> 475
>> 476 // Front side with normal approx. -X
>> 477 //
>> 478 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
>> 479 if (!good)
>> 480 {
>> 481 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 482 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 483 "Face at ~-X not planar.");
>> 484 }
>> 485
>> 486 // Back side iwth normal approx. +X
>> 487 //
>> 488 good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
>> 489 if (!good)
>> 490 {
>> 491 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 492 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 493 "Face at ~+X not planar.");
>> 494 }
>> 495 }
>> 496 else
>> 497 {
>> 498 G4cerr << "ERROR - G4Trap()::G4Trap(): " << GetName() << G4endl;
>> 499 G4Exception("G4Trap::G4Trap()", "InvalidSetup", FatalException,
>> 500 "Invalid length G4Trap parameters.");
>> 501 }
197 } 502 }
198 503
199 ////////////////////////////////////////////// << 504 ///////////////////////////////////////////////////////////////////////////
200 // 505 //
201 // Nominal constructor for G4Trap whose parame 506 // Nominal constructor for G4Trap whose parameters are to be set by
202 // a G4VParamaterisation later. Check and set 507 // a G4VParamaterisation later. Check and set half-widths as well as
203 // angles: final check of coplanarity 508 // angles: final check of coplanarity
204 509
205 G4Trap::G4Trap( const G4String& pName ) 510 G4Trap::G4Trap( const G4String& pName )
206 : G4CSGSolid (pName), halfCarTolerance(0.5*k << 511 : G4CSGSolid (pName),
207 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), << 512 fDz (1.),
208 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.) << 513 fTthetaCphi (0.),
209 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.) << 514 fTthetaSphi (0.),
>> 515 fDy1 (1.),
>> 516 fDx1 (1.),
>> 517 fDx2 (1.),
>> 518 fTalpha1 (0.),
>> 519 fDy2 (1.),
>> 520 fDx3 (1.),
>> 521 fDx4 (1.),
>> 522 fTalpha2 (0.)
210 { 523 {
211 MakePlanes(); << 524 MakePlanes();
212 } 525 }
213 526
214 ////////////////////////////////////////////// << 527 ///////////////////////////////////////////////////////////////////////
215 // 528 //
216 // Fake default constructor - sets only member 529 // Fake default constructor - sets only member data and allocates memory
217 // for usage restri 530 // for usage restricted to object persistency.
218 // 531 //
219 G4Trap::G4Trap( __void__& a ) 532 G4Trap::G4Trap( __void__& a )
220 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo << 533 : G4CSGSolid(a)
221 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), <<
222 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.) <<
223 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.) <<
224 { 534 {
225 MakePlanes(); <<
226 } 535 }
227 536
228 ////////////////////////////////////////////// << 537 ////////////////////////////////////////////////////////////////////////
229 // 538 //
230 // Destructor 539 // Destructor
231 540
232 G4Trap::~G4Trap() = default; << 541 G4Trap::~G4Trap()
233 <<
234 ////////////////////////////////////////////// <<
235 // <<
236 // Copy constructor <<
237 <<
238 G4Trap::G4Trap(const G4Trap& rhs) <<
239 : G4CSGSolid(rhs), halfCarTolerance(rhs.half <<
240 fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi) <<
241 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.f <<
242 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.f <<
243 { <<
244 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs <<
245 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs. <<
246 fTrapType = rhs.fTrapType; <<
247 } <<
248 <<
249 ////////////////////////////////////////////// <<
250 // <<
251 // Assignment operator <<
252 <<
253 G4Trap& G4Trap::operator = (const G4Trap& rhs) <<
254 { 542 {
255 // Check assignment to self <<
256 // <<
257 if (this == &rhs) { return *this; } <<
258 <<
259 // Copy base class data <<
260 // <<
261 G4CSGSolid::operator=(rhs); <<
262 <<
263 // Copy data <<
264 // <<
265 halfCarTolerance = rhs.halfCarTolerance; <<
266 fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi <<
267 fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs <<
268 fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs <<
269 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs <<
270 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs. <<
271 fTrapType = rhs.fTrapType; <<
272 return *this; <<
273 } 543 }
274 544
275 ////////////////////////////////////////////// << 545 ///////////////////////////////////////////////////////////////////////
276 // 546 //
277 // Set all parameters, as for constructor - ch 547 // Set all parameters, as for constructor - check and set half-widths
278 // as well as angles: final check of coplanari 548 // as well as angles: final check of coplanarity
279 549
280 void G4Trap::SetAllParameters ( G4double pDz, 550 void G4Trap::SetAllParameters ( G4double pDz,
281 G4double pThet 551 G4double pTheta,
282 G4double pPhi, 552 G4double pPhi,
283 G4double pDy1, 553 G4double pDy1,
284 G4double pDx1, 554 G4double pDx1,
285 G4double pDx2, 555 G4double pDx2,
286 G4double pAlp1 556 G4double pAlp1,
287 G4double pDy2, 557 G4double pDy2,
288 G4double pDx3, 558 G4double pDx3,
289 G4double pDx4, 559 G4double pDx4,
290 G4double pAlp2 560 G4double pAlp2 )
291 { 561 {
292 // Reset data of the base class << 562 fCubicVolume= 0.;
293 fCubicVolume = 0; << 563 fpPolyhedron = 0;
294 fSurfaceArea = 0; << 564 if ( pDz>0 && pDy1>0 && pDx1>0 && pDx2>0 && pDy2>0 && pDx3>0 && pDx4>0 )
295 fRebuildPolyhedron = true; << 565 {
296 << 566 fDz=pDz;
297 // Set parameters << 567 fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
298 fDz = pDz; << 568 fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
299 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi << 569
300 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi << 570 fDy1=pDy1;
301 << 571 fDx1=pDx1;
302 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp << 572 fDx2=pDx2;
303 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp << 573 fTalpha1=std::tan(pAlp1);
>> 574
>> 575 fDy2=pDy2;
>> 576 fDx3=pDx3;
>> 577 fDx4=pDx4;
>> 578 fTalpha2=std::tan(pAlp2);
304 579
305 CheckParameters(); << 580 MakePlanes();
306 MakePlanes(); << 581 }
307 } << 582 else
308 << 583 {
309 ////////////////////////////////////////////// << 584 G4cerr << "ERROR - G4Trap()::SetAllParameters(): " << GetName() << G4endl
310 // << 585 << " Invalid dimensions !" << G4endl
311 // Check length parameters << 586 << " X - "
312 << 587 << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
313 void G4Trap::CheckParameters() << 588 << " Y - " << pDy1 << ", " << pDy2 << G4endl
314 { << 589 << " Z - " << pDz << G4endl;
315 if (fDz<=0 || << 590 G4Exception("G4Trap::SetAllParameters()", "InvalidSetup",
316 fDy1<=0 || fDx1<=0 || fDx2<=0 || << 591 FatalException, "Invalid Length Parameters.");
317 fDy2<=0 || fDx3<=0 || fDx4<=0) <<
318 { <<
319 std::ostringstream message; <<
320 message << "Invalid Length Parameters for <<
321 << "\n X - " <<fDx1<<", "<<fDx2<< <<
322 << "\n Y - " <<fDy1<<", "<<fDy2 <<
323 << "\n Z - " <<fDz; <<
324 G4Exception("G4Trap::CheckParameters()", " <<
325 FatalException, message); <<
326 } 592 }
327 } 593 }
328 594
329 ////////////////////////////////////////////// 595 //////////////////////////////////////////////////////////////////////////
330 // 596 //
331 // Compute vertices and set side planes << 597 // Checking of coplanarity
332 << 598
333 void G4Trap::MakePlanes() << 599 G4bool G4Trap::MakePlanes()
334 { << 600 {
335 G4double DzTthetaCphi = fDz*fTthetaCphi; << 601 G4bool good = true;
336 G4double DzTthetaSphi = fDz*fTthetaSphi; << 602
337 G4double Dy1Talpha1 = fDy1*fTalpha1; << 603 G4ThreeVector pt[8] ;
338 G4double Dy2Talpha2 = fDy2*fTalpha2; << 604
339 << 605 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
340 G4ThreeVector pt[8] = << 606 -fDz*fTthetaSphi-fDy1,-fDz);
341 { << 607 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
342 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx << 608 -fDz*fTthetaSphi-fDy1,-fDz);
343 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx << 609 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
344 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx << 610 -fDz*fTthetaSphi+fDy1,-fDz);
345 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx << 611 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
346 G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx << 612 -fDz*fTthetaSphi+fDy1,-fDz);
347 G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx << 613 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
348 G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx << 614 +fDz*fTthetaSphi-fDy2,+fDz);
349 G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx << 615 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
350 }; << 616 +fDz*fTthetaSphi-fDy2,+fDz);
351 << 617 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
352 MakePlanes(pt); << 618 +fDz*fTthetaSphi+fDy2,+fDz);
353 } << 619 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
354 << 620 +fDz*fTthetaSphi+fDy2,+fDz);
355 ////////////////////////////////////////////// <<
356 // <<
357 // Set side planes, check planarity <<
358 621
359 void G4Trap::MakePlanes(const G4ThreeVector pt << 622 // Bottom side with normal approx. -Y
360 { << 623 //
361 constexpr G4int iface[4][4] = { {0,4,5,1}, { << 624 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]) ;
362 const static G4String side[4] = { "~-Y", "~+ << 625 if (!good)
>> 626 {
>> 627 G4cerr << "ERROR - G4Trap()::MakePlanes(): " << GetName() << G4endl;
>> 628 G4Exception("G4Trap::MakePlanes()", "InvalidSetup", FatalException,
>> 629 "Face at ~-Y not planar.");
>> 630 }
363 631
364 for (G4int i=0; i<4; ++i) << 632 // Top side with normal approx. +Y
>> 633 //
>> 634 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
>> 635 if (!good)
365 { 636 {
366 if (MakePlane(pt[iface[i][0]], << 637 G4cerr << "ERROR - G4Trap()::MakePlanes(): " << GetName() << G4endl;
367 pt[iface[i][1]], << 638 G4Exception("G4Trap::MakePlanes()", "InvalidSetup", FatalException,
368 pt[iface[i][2]], << 639 "Face at ~+Y not planar.");
369 pt[iface[i][3]], << 640 }
370 fPlanes[i])) continue; <<
371 641
372 // Non planar side face << 642 // Front side with normal approx. -X
373 G4ThreeVector normal(fPlanes[i].a,fPlanes[ << 643 //
374 G4double dmax = 0; << 644 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
375 for (G4int k=0; k<4; ++k) << 645 if (!good)
376 { << 646 {
377 G4double dist = normal.dot(pt[iface[i][k << 647 G4cerr << "ERROR - G4Trap()::MakePlanes(): " << GetName() << G4endl;
378 if (std::abs(dist) > std::abs(dmax)) dma << 648 G4Exception("G4Trap::MakePlanes()", "InvalidSetup", FatalException,
379 } << 649 "Face at ~-X not planar.");
380 std::ostringstream message; << 650 }
381 message << "Side face " << side[i] << " is << 651
382 << GetName() << "\nDiscrepancy: " << 652 // Back side iwth normal approx. +X
383 StreamInfo(message); << 653 //
384 G4Exception("G4Trap::MakePlanes()", "GeomS << 654 good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
385 FatalException, message); << 655 if ( !good )
>> 656 {
>> 657 G4cerr << "ERROR - G4Trap()::MakePlanes(): " << GetName() << G4endl;
>> 658 G4Exception("G4Trap::MakePlanes()", "InvalidSetup", FatalException,
>> 659 "Face at ~+X not planar");
386 } 660 }
387 661
388 // Re-compute parameters << 662 return good;
389 SetCachedValues(); <<
390 } 663 }
391 664
392 ////////////////////////////////////////////// << 665 //////////////////////////////////////////////////////////////////////////////
393 // 666 //
394 // Calculate the coef's of the plane p1->p2->p 667 // Calculate the coef's of the plane p1->p2->p3->p4->p1
395 // where the ThreeVectors 1-4 are in anti-cloc << 668 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed from
396 // from infront of the plane (i.e. from normal << 669 // infront of the plane (i.e. from normal direction).
397 // 670 //
398 // Return true if the points are coplanar, fal << 671 // Return true if the ThreeVectors are coplanar + set coef;s
>> 672 // false if ThreeVectors are not coplanar
399 673
400 G4bool G4Trap::MakePlane( const G4ThreeVector& 674 G4bool G4Trap::MakePlane( const G4ThreeVector& p1,
401 const G4ThreeVector& 675 const G4ThreeVector& p2,
402 const G4ThreeVector& 676 const G4ThreeVector& p3,
403 const G4ThreeVector& 677 const G4ThreeVector& p4,
404 TrapSidePlane& 678 TrapSidePlane& plane )
405 { 679 {
406 G4ThreeVector normal = ((p4 - p2).cross(p3 - << 680 G4double a, b, c, s;
407 if (std::abs(normal.x()) < DBL_EPSILON) norm << 681 G4ThreeVector v12, v13, v14, Vcross;
408 if (std::abs(normal.y()) < DBL_EPSILON) norm <<
409 if (std::abs(normal.z()) < DBL_EPSILON) norm <<
410 normal = normal.unit(); <<
411 <<
412 G4ThreeVector centre = (p1 + p2 + p3 + p4)*0 <<
413 plane.a = normal.x(); <<
414 plane.b = normal.y(); <<
415 plane.c = normal.z(); <<
416 plane.d = -normal.dot(centre); <<
417 <<
418 // compute distances and check planarity <<
419 G4double d1 = std::abs(normal.dot(p1) + plan <<
420 G4double d2 = std::abs(normal.dot(p2) + plan <<
421 G4double d3 = std::abs(normal.dot(p3) + plan <<
422 G4double d4 = std::abs(normal.dot(p4) + plan <<
423 G4double dmax = std::max(std::max(std::max(d <<
424 682
425 return dmax <= 1000 * kCarTolerance; << 683 G4bool good;
426 } <<
427 684
428 ////////////////////////////////////////////// << 685 v12 = p2 - p1;
429 // << 686 v13 = p3 - p1;
430 // Recompute parameters using planes << 687 v14 = p4 - p1;
>> 688 Vcross = v12.cross(v13);
431 689
432 void G4Trap::SetCachedValues() << 690 if (std::fabs(Vcross.dot(v14)/(Vcross.mag()*v14.mag())) > kCoplanar_Tolerance)
433 { <<
434 // Set indeces <<
435 constexpr G4int iface[6][4] = <<
436 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2, <<
437 <<
438 // Get vertices <<
439 G4ThreeVector pt[8]; <<
440 GetVertices(pt); <<
441 <<
442 // Set face areas <<
443 for (G4int i=0; i<6; ++i) <<
444 { <<
445 fAreas[i] = G4GeomTools::QuadAreaNormal(pt <<
446 pt <<
447 pt <<
448 pt <<
449 } <<
450 for (G4int i=1; i<6; ++i) { fAreas[i] += fAr <<
451 <<
452 // Define type of trapezoid <<
453 fTrapType = 0; <<
454 if (fPlanes[0].b == -1 && fPlanes[1].b == 1 <<
455 std::abs(fPlanes[0].a) < DBL_EPSILON && <<
456 std::abs(fPlanes[0].c) < DBL_EPSILON && <<
457 std::abs(fPlanes[1].a) < DBL_EPSILON && <<
458 std::abs(fPlanes[1].c) < DBL_EPSILON) <<
459 { <<
460 fTrapType = 1; // YZ section is a rectangl <<
461 if (std::abs(fPlanes[2].a + fPlanes[3].a) <<
462 std::abs(fPlanes[2].c - fPlanes[3].c) <<
463 fPlanes[2].b == 0 && <<
464 fPlanes[3].b == 0) <<
465 { <<
466 fTrapType = 2; // ... and XZ section is <<
467 fPlanes[2].a = -fPlanes[3].a; <<
468 fPlanes[2].c = fPlanes[3].c; <<
469 } <<
470 if (std::abs(fPlanes[2].a + fPlanes[3].a) <<
471 std::abs(fPlanes[2].b - fPlanes[3].b) <<
472 fPlanes[2].c == 0 && <<
473 fPlanes[3].c == 0) <<
474 { <<
475 fTrapType = 3; // ... and XY section is <<
476 fPlanes[2].a = -fPlanes[3].a; <<
477 fPlanes[2].b = fPlanes[3].b; <<
478 } <<
479 } <<
480 } <<
481 <<
482 ////////////////////////////////////////////// <<
483 // <<
484 // Get volume <<
485 <<
486 G4double G4Trap::GetCubicVolume() <<
487 { <<
488 if (fCubicVolume == 0) <<
489 { 691 {
490 G4ThreeVector pt[8]; << 692 good = false;
491 GetVertices(pt); <<
492 <<
493 G4double dz = pt[4].z() - pt[0].z(); <<
494 G4double dy1 = pt[2].y() - pt[0].y(); <<
495 G4double dx1 = pt[1].x() - pt[0].x(); <<
496 G4double dx2 = pt[3].x() - pt[2].x(); <<
497 G4double dy2 = pt[6].y() - pt[4].y(); <<
498 G4double dx3 = pt[5].x() - pt[4].x(); <<
499 G4double dx4 = pt[7].x() - pt[6].x(); <<
500 <<
501 fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(d <<
502 (dx4 + dx3 - dx2 - dx1)*(d <<
503 } 693 }
504 return fCubicVolume; << 694 else
505 } <<
506 <<
507 ////////////////////////////////////////////// <<
508 // <<
509 // Get surface area <<
510 <<
511 G4double G4Trap::GetSurfaceArea() <<
512 { <<
513 if (fSurfaceArea == 0) <<
514 { 695 {
515 G4ThreeVector pt[8]; << 696 // a,b,c correspond to the x/y/z components of the
516 G4int iface [6][4] = << 697 // normal vector to the plane
517 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2, << 698
518 << 699 // a = (p2.y()-p1.y())*(p1.z()+p2.z())+(p3.y()-p2.y())*(p2.z()+p3.z());
519 GetVertices(pt); << 700 // a += (p4.y()-p3.y())*(p3.z()+p4.z())+(p1.y()-p4.y())*(p4.z()+p1.z()); // ?
520 for (const auto & i : iface) << 701 // b = (p2.z()-p1.z())*(p1.x()+p2.x())+(p3.z()-p2.z())*(p2.x()+p3.x());
>> 702 // b += (p4.z()-p3.z())*(p3.x()+p4.x())+(p1.z()-p4.z())*(p4.x()+p1.x()); // ?
>> 703 // c = (p2.x()-p1.x())*(p1.y()+p2.y())+(p3.x()-p2.x())*(p2.y()+p3.y());
>> 704 // c += (p4.x()-p3.x())*(p3.y()+p4.y())+(p1.x()-p4.x())*(p4.y()+p1.y()); // ?
>> 705
>> 706 // Let create diagonals 4-2 and 3-1 than (4-2)x(3-1) provides
>> 707 // vector perpendicular to the plane directed to outside !!!
>> 708 // and a,b,c, = f(1,2,3,4) external relative to trap normal
>> 709
>> 710 a = +(p4.y() - p2.y())*(p3.z() - p1.z())
>> 711 - (p3.y() - p1.y())*(p4.z() - p2.z());
>> 712
>> 713 b = -(p4.x() - p2.x())*(p3.z() - p1.z())
>> 714 + (p3.x() - p1.x())*(p4.z() - p2.z());
>> 715
>> 716 c = +(p4.x() - p2.x())*(p3.y() - p1.y())
>> 717 - (p3.x() - p1.x())*(p4.y() - p2.y());
>> 718
>> 719 s = std::sqrt( a*a + b*b + c*c ); // so now vector plane.(a,b,c) is unit
>> 720
>> 721 if( s > 0 )
>> 722 {
>> 723 plane.a = a/s;
>> 724 plane.b = b/s;
>> 725 plane.c = c/s;
>> 726 }
>> 727 else
521 { 728 {
522 fSurfaceArea += G4GeomTools::QuadAreaNor << 729 G4cerr << "ERROR - G4Trap()::MakePlane(): " << GetName() << G4endl;
523 << 730 G4Exception("G4Trap::MakePlanes()", "InvalidSetup", FatalException,
524 << 731 "Invalid parameters: norm.mod() <= 0") ;
525 <<
526 } 732 }
>> 733 // Calculate D: p1 in in plane so D=-n.p1.Vect()
>> 734
>> 735 plane.d = -( plane.a*p1.x() + plane.b*p1.y() + plane.c*p1.z() );
>> 736
>> 737 good = true;
527 } 738 }
528 return fSurfaceArea; << 739 return good;
529 } 740 }
530 741
531 ////////////////////////////////////////////// << 742 //////////////////////////////////////////////////////////////////////////////
532 // 743 //
533 // Dispatch to parameterisation for replicatio 744 // Dispatch to parameterisation for replication mechanism dimension
534 // computation & modification. 745 // computation & modification.
535 746
536 void G4Trap::ComputeDimensions( G4VPVPar 747 void G4Trap::ComputeDimensions( G4VPVParameterisation* p,
537 const G4int n, 748 const G4int n,
538 const G4VPhysi 749 const G4VPhysicalVolume* pRep )
539 { 750 {
540 p->ComputeDimensions(*this,n,pRep); 751 p->ComputeDimensions(*this,n,pRep);
541 } 752 }
542 753
543 ////////////////////////////////////////////// <<
544 // <<
545 // Get bounding box <<
546 <<
547 void G4Trap::BoundingLimits(G4ThreeVector& pMi <<
548 { <<
549 G4ThreeVector pt[8]; <<
550 GetVertices(pt); <<
551 754
552 G4double xmin = kInfinity, xmax = -kInfinity << 755 ////////////////////////////////////////////////////////////////////////
553 G4double ymin = kInfinity, ymax = -kInfinity <<
554 for (const auto & i : pt) <<
555 { <<
556 G4double x = i.x(); <<
557 if (x < xmin) xmin = x; <<
558 if (x > xmax) xmax = x; <<
559 G4double y = i.y(); <<
560 if (y < ymin) ymin = y; <<
561 if (y > ymax) ymax = y; <<
562 } <<
563 <<
564 G4double dz = GetZHalfLength(); <<
565 pMin.set(xmin,ymin,-dz); <<
566 pMax.set(xmax,ymax, dz); <<
567 <<
568 // Check correctness of the bounding box <<
569 // <<
570 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
571 { <<
572 std::ostringstream message; <<
573 message << "Bad bounding box (min >= max) <<
574 << GetName() << " !" <<
575 << "\npMin = " << pMin <<
576 << "\npMax = " << pMax; <<
577 G4Exception("G4Trap::BoundingLimits()", "G <<
578 JustWarning, message); <<
579 DumpInfo(); <<
580 } <<
581 } <<
582 <<
583 ////////////////////////////////////////////// <<
584 // 756 //
585 // Calculate extent under transform and specif 757 // Calculate extent under transform and specified limit
586 758
587 G4bool G4Trap::CalculateExtent( const EAxis pA 759 G4bool G4Trap::CalculateExtent( const EAxis pAxis,
588 const G4VoxelL 760 const G4VoxelLimits& pVoxelLimit,
589 const G4Affine 761 const G4AffineTransform& pTransform,
590 G4double 762 G4double& pMin, G4double& pMax) const
591 { 763 {
592 G4ThreeVector bmin, bmax; << 764 G4double xMin, xMax, yMin, yMax, zMin, zMax;
593 G4bool exist; << 765 G4bool flag;
594 766
595 // Check bounding box (bbox) << 767 if (!pTransform.IsRotated())
596 // << 768 {
597 BoundingLimits(bmin,bmax); << 769 // Special case handling for unrotated trapezoids
598 G4BoundingEnvelope bbox(bmin,bmax); << 770 // Compute z/x/y/ mins and maxs respecting limits, with early returns
599 #ifdef G4BBOX_EXTENT << 771 // if outside limits. Then switch() on pAxis
600 return bbox.CalculateExtent(pAxis,pVoxelLimi << 772
601 #endif << 773 G4int i ;
602 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 774 G4double xoffset;
603 { << 775 G4double yoffset;
604 return exist = pMin < pMax; << 776 G4double zoffset;
>> 777 G4double temp[8] ; // some points for intersection with zMin/zMax
>> 778 G4ThreeVector pt[8]; // vertices after translation
>> 779
>> 780 xoffset=pTransform.NetTranslation().x();
>> 781 yoffset=pTransform.NetTranslation().y();
>> 782 zoffset=pTransform.NetTranslation().z();
>> 783
>> 784 pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
>> 785 yoffset-fDz*fTthetaSphi-fDy1,zoffset-fDz);
>> 786 pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
>> 787 yoffset-fDz*fTthetaSphi-fDy1,zoffset-fDz);
>> 788 pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
>> 789 yoffset-fDz*fTthetaSphi+fDy1,zoffset-fDz);
>> 790 pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
>> 791 yoffset-fDz*fTthetaSphi+fDy1,zoffset-fDz);
>> 792 pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
>> 793 yoffset+fDz*fTthetaSphi-fDy2,zoffset+fDz);
>> 794 pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
>> 795 yoffset+fDz*fTthetaSphi-fDy2,zoffset+fDz);
>> 796 pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
>> 797 yoffset+fDz*fTthetaSphi+fDy2,zoffset+fDz);
>> 798 pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
>> 799 yoffset+fDz*fTthetaSphi+fDy2,zoffset+fDz);
>> 800 zMin=zoffset-fDz;
>> 801 zMax=zoffset+fDz;
>> 802
>> 803 if ( pVoxelLimit.IsZLimited() )
>> 804 {
>> 805 if ( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance)
>> 806 || (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
>> 807 {
>> 808 return false;
>> 809 }
>> 810 else
>> 811 {
>> 812 if ( zMin < pVoxelLimit.GetMinZExtent() )
>> 813 {
>> 814 zMin = pVoxelLimit.GetMinZExtent() ;
>> 815 }
>> 816 if ( zMax > pVoxelLimit.GetMaxZExtent() )
>> 817 {
>> 818 zMax = pVoxelLimit.GetMaxZExtent() ;
>> 819 }
>> 820 }
>> 821 }
>> 822 temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMin-pt[0].z())
>> 823 /(pt[4].z()-pt[0].z()) ;
>> 824 temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMax-pt[0].z())
>> 825 /(pt[4].z()-pt[0].z()) ;
>> 826 temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMin-pt[2].z())
>> 827 /(pt[6].z()-pt[2].z()) ;
>> 828 temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMax-pt[2].z())
>> 829 /(pt[6].z()-pt[2].z()) ;
>> 830
>> 831 yMax = yoffset - std::fabs(fDz*fTthetaSphi) - fDy1 - fDy2 ;
>> 832 yMin = -yMax ;
>> 833
>> 834 for( i = 0 ; i < 4 ; i++ )
>> 835 {
>> 836 if( temp[i] > yMax ) yMax = temp[i] ;
>> 837 if( temp[i] < yMin ) yMin = temp[i] ;
>> 838 }
>> 839 if ( pVoxelLimit.IsYLimited() )
>> 840 {
>> 841 if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance)
>> 842 || (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
>> 843 {
>> 844 return false;
>> 845 }
>> 846 else
>> 847 {
>> 848 if ( yMin < pVoxelLimit.GetMinYExtent() )
>> 849 {
>> 850 yMin = pVoxelLimit.GetMinYExtent() ;
>> 851 }
>> 852 if ( yMax > pVoxelLimit.GetMaxYExtent() )
>> 853 {
>> 854 yMax = pVoxelLimit.GetMaxYExtent() ;
>> 855 }
>> 856 }
>> 857 }
>> 858 temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())
>> 859 *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
>> 860 temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())
>> 861 *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
>> 862 temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())
>> 863 *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
>> 864 temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())
>> 865 *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
>> 866 temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())
>> 867 *(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
>> 868 temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())
>> 869 *(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
>> 870 temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())
>> 871 *(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
>> 872 temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())
>> 873 *(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
>> 874
>> 875 xMax = xoffset - std::fabs(fDz*fTthetaCphi) - fDx1 - fDx2 -fDx3 - fDx4 ;
>> 876 xMin = -xMax ;
>> 877
>> 878 for( i = 0 ; i < 8 ; i++ )
>> 879 {
>> 880 if( temp[i] > xMax) xMax = temp[i] ;
>> 881 if( temp[i] < xMin) xMin = temp[i] ;
>> 882 }
>> 883 if (pVoxelLimit.IsXLimited()) // xMax/Min = f(yMax/Min) ?
>> 884 {
>> 885 if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance)
>> 886 || (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
>> 887 {
>> 888 return false;
>> 889 }
>> 890 else
>> 891 {
>> 892 if ( xMin < pVoxelLimit.GetMinXExtent() )
>> 893 {
>> 894 xMin = pVoxelLimit.GetMinXExtent() ;
>> 895 }
>> 896 if ( xMax > pVoxelLimit.GetMaxXExtent() )
>> 897 {
>> 898 xMax = pVoxelLimit.GetMaxXExtent() ;
>> 899 }
>> 900 }
>> 901 }
>> 902 switch (pAxis)
>> 903 {
>> 904 case kXAxis:
>> 905 pMin=xMin;
>> 906 pMax=xMax;
>> 907 break;
>> 908
>> 909 case kYAxis:
>> 910 pMin=yMin;
>> 911 pMax=yMax;
>> 912 break;
>> 913
>> 914 case kZAxis:
>> 915 pMin=zMin;
>> 916 pMax=zMax;
>> 917 break;
>> 918
>> 919 default:
>> 920 break;
>> 921 }
>> 922 pMin -= kCarTolerance;
>> 923 pMax += kCarTolerance;
>> 924
>> 925 flag = true;
605 } 926 }
>> 927 else // General rotated case -
>> 928 {
>> 929 G4bool existsAfterClip = false ;
>> 930 G4ThreeVectorList* vertices;
>> 931 pMin = +kInfinity;
>> 932 pMax = -kInfinity;
>> 933
>> 934 // Calculate rotated vertex coordinates. Operator 'new' is called
>> 935
>> 936 vertices = CreateRotatedVertices(pTransform);
>> 937
>> 938 xMin = +kInfinity; yMin = +kInfinity; zMin = +kInfinity;
>> 939 xMax = -kInfinity; yMax = -kInfinity; zMax = -kInfinity;
>> 940
>> 941 for( G4int nv = 0 ; nv < 8 ; nv++ )
>> 942 {
>> 943 if( (*vertices)[nv].x() > xMax ) xMax = (*vertices)[nv].x();
>> 944 if( (*vertices)[nv].y() > yMax ) yMax = (*vertices)[nv].y();
>> 945 if( (*vertices)[nv].z() > zMax ) zMax = (*vertices)[nv].z();
>> 946
>> 947 if( (*vertices)[nv].x() < xMin ) xMin = (*vertices)[nv].x();
>> 948 if( (*vertices)[nv].y() < yMin ) yMin = (*vertices)[nv].y();
>> 949 if( (*vertices)[nv].z() < zMin ) zMin = (*vertices)[nv].z();
>> 950 }
>> 951 if ( pVoxelLimit.IsZLimited() )
>> 952 {
>> 953 if ( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance)
>> 954 || (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
>> 955 {
>> 956 delete vertices ; // 'new' in the function called
>> 957 return false;
>> 958 }
>> 959 else
>> 960 {
>> 961 if ( zMin < pVoxelLimit.GetMinZExtent() )
>> 962 {
>> 963 zMin = pVoxelLimit.GetMinZExtent() ;
>> 964 }
>> 965 if ( zMax > pVoxelLimit.GetMaxZExtent() )
>> 966 {
>> 967 zMax = pVoxelLimit.GetMaxZExtent() ;
>> 968 }
>> 969 }
>> 970 }
>> 971 if ( pVoxelLimit.IsYLimited() )
>> 972 {
>> 973 if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance)
>> 974 || (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
>> 975 {
>> 976 delete vertices ; // 'new' in the function called
>> 977 return false;
>> 978 }
>> 979 else
>> 980 {
>> 981 if ( yMin < pVoxelLimit.GetMinYExtent() )
>> 982 {
>> 983 yMin = pVoxelLimit.GetMinYExtent() ;
>> 984 }
>> 985 if ( yMax > pVoxelLimit.GetMaxYExtent() )
>> 986 {
>> 987 yMax = pVoxelLimit.GetMaxYExtent() ;
>> 988 }
>> 989 }
>> 990 }
>> 991 if ( pVoxelLimit.IsXLimited() )
>> 992 {
>> 993 if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance)
>> 994 || (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
>> 995 {
>> 996 delete vertices ; // 'new' in the function called
>> 997 return false ;
>> 998 }
>> 999 else
>> 1000 {
>> 1001 if ( xMin < pVoxelLimit.GetMinXExtent() )
>> 1002 {
>> 1003 xMin = pVoxelLimit.GetMinXExtent() ;
>> 1004 }
>> 1005 if ( xMax > pVoxelLimit.GetMaxXExtent() )
>> 1006 {
>> 1007 xMax = pVoxelLimit.GetMaxXExtent() ;
>> 1008 }
>> 1009 }
>> 1010 }
>> 1011 switch (pAxis)
>> 1012 {
>> 1013 case kXAxis:
>> 1014 pMin=xMin;
>> 1015 pMax=xMax;
>> 1016 break;
606 1017
607 // Set bounding envelope (benv) and calculat << 1018 case kYAxis:
608 // << 1019 pMin=yMin;
609 G4ThreeVector pt[8]; << 1020 pMax=yMax;
610 GetVertices(pt); << 1021 break;
611 1022
612 G4ThreeVectorList baseA(4), baseB(4); << 1023 case kZAxis:
613 baseA[0] = pt[0]; << 1024 pMin=zMin;
614 baseA[1] = pt[1]; << 1025 pMax=zMax;
615 baseA[2] = pt[3]; << 1026 break;
616 baseA[3] = pt[2]; << 1027
617 << 1028 default:
618 baseB[0] = pt[4]; << 1029 break;
619 baseB[1] = pt[5]; << 1030 }
620 baseB[2] = pt[7]; << 1031 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
621 baseB[3] = pt[6]; << 1032 {
622 << 1033 existsAfterClip=true;
623 std::vector<const G4ThreeVectorList *> polyg << 1034
624 polygons[0] = &baseA; << 1035 // Add tolerance to avoid precision troubles
625 polygons[1] = &baseB; << 1036 //
626 << 1037 pMin -= kCarTolerance ;
627 G4BoundingEnvelope benv(bmin,bmax,polygons); << 1038 pMax += kCarTolerance ;
628 exist = benv.CalculateExtent(pAxis,pVoxelLim << 1039 }
629 return exist; << 1040 delete vertices ; // 'new' in the function called
>> 1041 flag = existsAfterClip ;
>> 1042 }
>> 1043 return flag;
630 } 1044 }
631 1045
632 ////////////////////////////////////////////// << 1046
>> 1047 ////////////////////////////////////////////////////////////////////////
633 // 1048 //
634 // Return whether point is inside/outside/on_s << 1049 // Return whether point inside/outside/on surface, using tolerance
635 1050
636 EInside G4Trap::Inside( const G4ThreeVector& p 1051 EInside G4Trap::Inside( const G4ThreeVector& p ) const
637 { 1052 {
638 switch (fTrapType) << 1053 EInside in;
>> 1054 G4double Dist;
>> 1055 G4int i;
>> 1056 if ( std::fabs(p.z()) <= fDz-kCarTolerance*0.5)
639 { 1057 {
640 case 0: // General case << 1058 in = kInside;
>> 1059
>> 1060 for ( i = 0;i < 4;i++ )
>> 1061 {
>> 1062 Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
>> 1063 +fPlanes[i].c*p.z() + fPlanes[i].d;
>> 1064
>> 1065 if (Dist > kCarTolerance*0.5) return in = kOutside;
>> 1066 else if (Dist > -kCarTolerance*0.5) in = kSurface;
>> 1067
>> 1068 }
>> 1069 }
>> 1070 else if (std::fabs(p.z()) <= fDz+kCarTolerance*0.5)
>> 1071 {
>> 1072 in = kSurface;
>> 1073
>> 1074 for ( i = 0; i < 4; i++ )
641 { 1075 {
642 G4double dz = std::abs(p.z())-fDz; << 1076 Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
643 G4double dy1 = fPlanes[0].b*p.y()+fPlane << 1077 +fPlanes[i].c*p.z() + fPlanes[i].d;
644 G4double dy2 = fPlanes[1].b*p.y()+fPlane <<
645 G4double dy = std::max(dz,std::max(dy1,d <<
646 <<
647 G4double dx1 = fPlanes[2].a*p.x()+fPlane <<
648 + fPlanes[2].c*p.z()+fPlane <<
649 G4double dx2 = fPlanes[3].a*p.x()+fPlane <<
650 + fPlanes[3].c*p.z()+fPlane <<
651 G4double dist = std::max(dy,std::max(dx1 <<
652 <<
653 return (dist > halfCarTolerance) ? kOuts <<
654 ((dist > -halfCarTolerance) ? kSurface <<
655 } <<
656 case 1: // YZ section is a rectangle <<
657 { <<
658 G4double dz = std::abs(p.z())-fDz; <<
659 G4double dy = std::max(dz,std::abs(p.y() <<
660 G4double dx1 = fPlanes[2].a*p.x()+fPlane <<
661 + fPlanes[2].c*p.z()+fPlane <<
662 G4double dx2 = fPlanes[3].a*p.x()+fPlane <<
663 + fPlanes[3].c*p.z()+fPlane <<
664 G4double dist = std::max(dy,std::max(dx1 <<
665 <<
666 return (dist > halfCarTolerance) ? kOuts <<
667 ((dist > -halfCarTolerance) ? kSurface <<
668 } <<
669 case 2: // YZ section is a rectangle and <<
670 { // XZ section is an isosceles trap <<
671 G4double dz = std::abs(p.z())-fDz; <<
672 G4double dy = std::max(dz,std::abs(p.y() <<
673 G4double dx = fPlanes[3].a*std::abs(p.x( <<
674 + fPlanes[3].c*p.z()+fPlanes <<
675 G4double dist = std::max(dy,dx); <<
676 <<
677 return (dist > halfCarTolerance) ? kOuts <<
678 ((dist > -halfCarTolerance) ? kSurface <<
679 } <<
680 case 3: // YZ section is a rectangle and <<
681 { // XY section is an isosceles trap <<
682 G4double dz = std::abs(p.z())-fDz; <<
683 G4double dy = std::max(dz,std::abs(p.y() <<
684 G4double dx = fPlanes[3].a*std::abs(p.x( <<
685 + fPlanes[3].b*p.y()+fPlanes <<
686 G4double dist = std::max(dy,dx); <<
687 1078
688 return (dist > halfCarTolerance) ? kOuts << 1079 if (Dist > kCarTolerance*0.5) return in = kOutside;
689 ((dist > -halfCarTolerance) ? kSurface <<
690 } 1080 }
691 } 1081 }
692 return kOutside; << 1082 else in = kOutside;
>> 1083
>> 1084 return in;
693 } 1085 }
694 1086
695 ////////////////////////////////////////////// << 1087 /////////////////////////////////////////////////////////////////////////////
696 // 1088 //
697 // Determine side, and return corresponding no << 1089 // Calculate side nearest to p, and return normal
>> 1090 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
698 1091
699 G4ThreeVector G4Trap::SurfaceNormal( const G4T 1092 G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const
700 { 1093 {
701 G4double nx = 0, ny = 0, nz = 0; << 1094 G4int i, imin = 0, noSurfaces = 0;
702 G4double dz = std::abs(p.z()) - fDz; << 1095 G4double dist, distz, distx, disty, distmx, distmy, safe = kInfinity;
703 nz = std::copysign(G4double(std::abs(dz) <= << 1096 G4double delta = 0.5*kCarTolerance;
>> 1097 G4ThreeVector norm, sumnorm(0.,0.,0.);
704 1098
705 switch (fTrapType) << 1099 for (i = 0; i < 4; i++)
706 { 1100 {
707 case 0: // General case << 1101 dist = std::fabs(fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
708 { << 1102 + fPlanes[i].c*p.z() + fPlanes[i].d);
709 for (G4int i=0; i<2; ++i) << 1103 if ( dist < safe )
710 { <<
711 G4double dy = fPlanes[i].b*p.y() + fPl <<
712 if (std::abs(dy) > halfCarTolerance) c <<
713 ny = fPlanes[i].b; <<
714 nz += fPlanes[i].c; <<
715 break; <<
716 } <<
717 for (G4int i=2; i<4; ++i) <<
718 { <<
719 G4double dx = fPlanes[i].a*p.x() + <<
720 fPlanes[i].b*p.y() + fPl <<
721 if (std::abs(dx) > halfCarTolerance) c <<
722 nx = fPlanes[i].a; <<
723 ny += fPlanes[i].b; <<
724 nz += fPlanes[i].c; <<
725 break; <<
726 } <<
727 break; <<
728 } <<
729 case 1: // YZ section - rectangle <<
730 { 1104 {
731 G4double dy = std::abs(p.y()) + fPlanes[ << 1105 safe = dist;
732 ny = std::copysign(G4double(std::abs(dy) << 1106 imin = i;
733 for (G4int i=2; i<4; ++i) <<
734 { <<
735 G4double dx = fPlanes[i].a*p.x() + <<
736 fPlanes[i].b*p.y() + fPl <<
737 if (std::abs(dx) > halfCarTolerance) c <<
738 nx = fPlanes[i].a; <<
739 ny += fPlanes[i].b; <<
740 nz += fPlanes[i].c; <<
741 break; <<
742 } <<
743 break; <<
744 } <<
745 case 2: // YZ section - rectangle, XZ sect <<
746 { <<
747 G4double dy = std::abs(p.y()) + fPlanes[ <<
748 ny = std::copysign(G4double(std::abs(dy) <<
749 G4double dx = fPlanes[3].a*std::abs(p.x( <<
750 fPlanes[3].c*p.z() + fPlan <<
751 G4double k = std::abs(dx) <= halfCarTole <<
752 nx = std::copysign(k, p.x())*fPlanes[3] <<
753 nz += k*fPlanes[3].c; <<
754 break; <<
755 } <<
756 case 3: // YZ section - rectangle, XY sect <<
757 { <<
758 G4double dy = std::abs(p.y()) + fPlanes[ <<
759 ny = std::copysign(G4double(std::abs(dy) <<
760 G4double dx = fPlanes[3].a*std::abs(p.x( <<
761 fPlanes[3].b*p.y() + fPlan <<
762 G4double k = std::abs(dx) <= halfCarTole <<
763 nx = std::copysign(k, p.x())*fPlanes[3] <<
764 ny += k*fPlanes[3].b; <<
765 break; <<
766 } 1107 }
767 } 1108 }
>> 1109 distz = std::fabs( std::fabs( p.z() ) - fDz );
768 1110
769 // Return normal << 1111 distmy = std::fabs( fPlanes[0].a*p.x() + fPlanes[0].b*p.y()
770 // << 1112 + fPlanes[0].c*p.z() + fPlanes[0].d );
771 G4double mag2 = nx*nx + ny*ny + nz*nz; << 1113
772 if (mag2 == 1) return { nx,ny,nz }; << 1114 disty = std::fabs( fPlanes[1].a*p.x() + fPlanes[1].b*p.y()
773 else if (mag2 != 0) return G4ThreeVector(nx, << 1115 + fPlanes[1].c*p.z() + fPlanes[1].d );
774 else << 1116
>> 1117 distmx = std::fabs( fPlanes[2].a*p.x() + fPlanes[2].b*p.y()
>> 1118 + fPlanes[2].c*p.z() + fPlanes[2].d );
>> 1119
>> 1120 distx = std::fabs( fPlanes[3].a*p.x() + fPlanes[3].b*p.y()
>> 1121 + fPlanes[3].c*p.z() + fPlanes[3].d );
>> 1122
>> 1123 G4ThreeVector nX = G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
>> 1124 G4ThreeVector nmX = G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
>> 1125 G4ThreeVector nY = G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
>> 1126 G4ThreeVector nmY = G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
>> 1127 G4ThreeVector nZ = G4ThreeVector(0.,0.,1.0);
>> 1128
>> 1129 if (distx <= delta)
>> 1130 {
>> 1131 noSurfaces ++;
>> 1132 sumnorm += nX;
>> 1133 }
>> 1134 if (distmx <= delta)
>> 1135 {
>> 1136 noSurfaces ++;
>> 1137 sumnorm += nmX;
>> 1138 }
>> 1139 if (disty <= delta)
>> 1140 {
>> 1141 noSurfaces ++;
>> 1142 sumnorm += nY;
>> 1143 }
>> 1144 if (distmy <= delta)
>> 1145 {
>> 1146 noSurfaces ++;
>> 1147 sumnorm += nmY;
>> 1148 }
>> 1149 if (distz <= delta)
>> 1150 {
>> 1151 noSurfaces ++;
>> 1152 if ( p.z() >= 0.) sumnorm += nZ;
>> 1153 else sumnorm -= nZ;
>> 1154 }
>> 1155 if ( noSurfaces == 0 )
775 { 1156 {
776 // Point is not on the surface <<
777 // <<
778 #ifdef G4CSGDEBUG 1157 #ifdef G4CSGDEBUG
779 std::ostringstream message; << 1158 G4Exception("G4Trap::SurfaceNormal(p)", "Notification", JustWarning,
780 G4long oldprc = message.precision(16); << 1159 "Point p is not on surface !?" );
781 message << "Point p is not on surface (!?) << 1160 #endif
782 << GetName() << G4endl; << 1161 norm = ApproxSurfaceNormal(p);
783 message << "Position:\n"; <<
784 message << " p.x() = " << p.x()/mm << " <<
785 message << " p.y() = " << p.y()/mm << " <<
786 message << " p.z() = " << p.z()/mm << " <<
787 G4cout.precision(oldprc) ; <<
788 G4Exception("G4Trap::SurfaceNormal(p)", "G <<
789 JustWarning, message ); <<
790 DumpInfo(); <<
791 #endif <<
792 return ApproxSurfaceNormal(p); <<
793 } 1162 }
>> 1163 else if ( noSurfaces == 1 ) norm = sumnorm;
>> 1164 else norm = sumnorm.unit();
>> 1165 return norm;
794 } 1166 }
795 1167
796 ////////////////////////////////////////////// << 1168 ////////////////////////////////////////////////////////////////////////////////////
797 // 1169 //
798 // Algorithm for SurfaceNormal() following the 1170 // Algorithm for SurfaceNormal() following the original specification
799 // for points not on the surface 1171 // for points not on the surface
800 1172
801 G4ThreeVector G4Trap::ApproxSurfaceNormal( con 1173 G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
802 { 1174 {
803 G4double dist = -DBL_MAX; << 1175 G4double safe=kInfinity,Dist,safez;
804 G4int iside = 0; << 1176 G4int i,imin=0;
805 for (G4int i=0; i<4; ++i) << 1177 for (i=0;i<4;i++)
806 { << 1178 {
807 G4double d = fPlanes[i].a*p.x() + << 1179 Dist=std::fabs(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
808 fPlanes[i].b*p.y() + << 1180 +fPlanes[i].c*p.z()+fPlanes[i].d);
809 fPlanes[i].c*p.z() + fPlanes[ << 1181 if (Dist<safe)
810 if (d > dist) { dist = d; iside = i; } << 1182 {
>> 1183 safe=Dist;
>> 1184 imin=i;
>> 1185 }
>> 1186 }
>> 1187 safez=std::fabs(std::fabs(p.z())-fDz);
>> 1188 if (safe<safez)
>> 1189 {
>> 1190 return G4ThreeVector(fPlanes[imin].a,fPlanes[imin].b,fPlanes[imin].c);
811 } 1191 }
812 <<
813 G4double distz = std::abs(p.z()) - fDz; <<
814 if (dist > distz) <<
815 return { fPlanes[iside].a, fPlanes[iside]. <<
816 else 1192 else
817 return { 0, 0, (G4double)((p.z() < 0) ? -1 << 1193 {
>> 1194 if (p.z()>0)
>> 1195 {
>> 1196 return G4ThreeVector(0,0,1);
>> 1197 }
>> 1198 else
>> 1199 {
>> 1200 return G4ThreeVector(0,0,-1);
>> 1201 }
>> 1202 }
818 } 1203 }
819 1204
820 ////////////////////////////////////////////// << 1205 ////////////////////////////////////////////////////////////////////////////
>> 1206 //
>> 1207 // Calculate distance to shape from outside - return kInfinity if no intersection
821 // 1208 //
822 // Calculate distance to shape from outside << 1209 // ALGORITHM:
823 // - return kInfinity if no intersection << 1210 // For each component, calculate pair of minimum and maximum intersection
>> 1211 // values for which the particle is in the extent of the shape
>> 1212 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
824 1213
825 G4double G4Trap::DistanceToIn(const G4ThreeVec << 1214 G4double G4Trap::DistanceToIn( const G4ThreeVector& p,
826 const G4ThreeVec << 1215 const G4ThreeVector& v ) const
827 { 1216 {
828 // Z intersections <<
829 // <<
830 if ((std::abs(p.z()) - fDz) >= -halfCarToler <<
831 return kInfinity; <<
832 G4double invz = (-v.z() == 0) ? DBL_MAX : -1 <<
833 G4double dz = (invz < 0) ? fDz : -fDz; <<
834 G4double tzmin = (p.z() + dz)*invz; <<
835 G4double tzmax = (p.z() - dz)*invz; <<
836 1217
837 // Y intersections << 1218 G4double snxt; // snxt = default return value
>> 1219 G4double max,smax,smin;
>> 1220 G4double pdist,Comp,vdist;
>> 1221 G4int i;
>> 1222 //
>> 1223 // Z Intersection range
838 // 1224 //
839 G4double tymin = 0, tymax = DBL_MAX; << 1225 if ( v.z() > 0 )
840 G4int i = 0; <<
841 for ( ; i<2; ++i) <<
842 { 1226 {
843 G4double cosa = fPlanes[i].b*v.y() + fPlan << 1227 max = fDz - p.z() ;
844 G4double dist = fPlanes[i].b*p.y() + fPlan << 1228 if (max > 0.5*kCarTolerance)
845 if (dist >= -halfCarTolerance) <<
846 { 1229 {
847 if (cosa >= 0) return kInfinity; << 1230 smax = max/v.z();
848 G4double tmp = -dist/cosa; << 1231 smin = (-fDz-p.z())/v.z();
849 if (tymin < tmp) tymin = tmp; <<
850 } 1232 }
851 else if (cosa > 0) << 1233 else
852 { 1234 {
853 G4double tmp = -dist/cosa; << 1235 return snxt=kInfinity;
854 if (tymax > tmp) tymax = tmp; <<
855 } 1236 }
856 } 1237 }
857 << 1238 else if (v.z() < 0 )
858 // Z intersections <<
859 // <<
860 G4double txmin = 0, txmax = DBL_MAX; <<
861 for ( ; i<4; ++i) <<
862 { 1239 {
863 G4double cosa = fPlanes[i].a*v.x()+fPlanes << 1240 max = - fDz - p.z() ;
864 G4double dist = fPlanes[i].a*p.x()+fPlanes << 1241 if (max < -0.5*kCarTolerance )
865 fPlanes[i].d; <<
866 if (dist >= -halfCarTolerance) <<
867 { 1242 {
868 if (cosa >= 0) return kInfinity; << 1243 smax=max/v.z();
869 G4double tmp = -dist/cosa; << 1244 smin=(fDz-p.z())/v.z();
870 if (txmin < tmp) txmin = tmp; <<
871 } 1245 }
872 else if (cosa > 0) << 1246 else
873 { 1247 {
874 G4double tmp = -dist/cosa; << 1248 return snxt=kInfinity;
875 if (txmax > tmp) txmax = tmp; << 1249 }
>> 1250 }
>> 1251 else
>> 1252 {
>> 1253 if (std::fabs(p.z())<fDz - 0.5*kCarTolerance) // Inside was <=fDz
>> 1254 {
>> 1255 smin=0;
>> 1256 smax=kInfinity;
>> 1257 }
>> 1258 else
>> 1259 {
>> 1260 return snxt=kInfinity;
876 } 1261 }
877 } 1262 }
878 1263
879 // Find distance << 1264 for (i=0;i<4;i++)
>> 1265 {
>> 1266 pdist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
>> 1267 +fPlanes[i].c*p.z()+fPlanes[i].d;
>> 1268 Comp=fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
>> 1269 if ( pdist >= -0.5*kCarTolerance ) // was >0
>> 1270 {
>> 1271 //
>> 1272 // Outside the plane -> this is an extent entry distance
>> 1273 //
>> 1274 if (Comp >= 0) // was >0
>> 1275 {
>> 1276 return snxt=kInfinity ;
>> 1277 }
>> 1278 else
>> 1279 {
>> 1280 vdist=-pdist/Comp;
>> 1281 if (vdist>smin)
>> 1282 {
>> 1283 if (vdist<smax)
>> 1284 {
>> 1285 smin = vdist;
>> 1286 }
>> 1287 else
>> 1288 {
>> 1289 return snxt=kInfinity;
>> 1290 }
>> 1291 }
>> 1292 }
>> 1293 }
>> 1294 else
>> 1295 {
>> 1296 //
>> 1297 // Inside the plane -> couble be an extent exit distance (smax)
>> 1298 //
>> 1299 if (Comp>0) // Will leave extent
>> 1300 {
>> 1301 vdist=-pdist/Comp;
>> 1302 if (vdist<smax)
>> 1303 {
>> 1304 if (vdist>smin)
>> 1305 {
>> 1306 smax=vdist;
>> 1307 }
>> 1308 else
>> 1309 {
>> 1310 return snxt=kInfinity;
>> 1311 }
>> 1312 }
>> 1313 }
>> 1314 }
>> 1315 }
880 // 1316 //
881 G4double tmin = std::max(std::max(txmin,tymi << 1317 // Checks in non z plane intersections ensure smin<smax
882 G4double tmax = std::min(std::min(txmax,tyma << 1318 //
883 << 1319 if (smin >=0 )
884 if (tmax <= tmin + halfCarTolerance) return << 1320 {
885 return (tmin < halfCarTolerance ) ? 0. : tmi << 1321 snxt = smin ;
>> 1322 }
>> 1323 else
>> 1324 {
>> 1325 snxt = 0 ;
>> 1326 }
>> 1327 return snxt;
886 } 1328 }
887 1329
888 ////////////////////////////////////////////// << 1330 ///////////////////////////////////////////////////////////////////////////
889 // 1331 //
890 // Calculate exact shortest distance to any bo 1332 // Calculate exact shortest distance to any boundary from outside
891 // This is the best fast estimation of the sho 1333 // This is the best fast estimation of the shortest distance to trap
892 // - return 0 if point is inside << 1334 // - Returns 0 is ThreeVector inside
893 1335
894 G4double G4Trap::DistanceToIn( const G4ThreeVe 1336 G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const
895 { 1337 {
896 switch (fTrapType) << 1338 G4double safe=0.0,Dist;
>> 1339 G4int i;
>> 1340 safe=std::fabs(p.z())-fDz;
>> 1341 for (i=0;i<4;i++)
897 { 1342 {
898 case 0: // General case << 1343 Dist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
899 { << 1344 +fPlanes[i].c*p.z()+fPlanes[i].d;
900 G4double dz = std::abs(p.z())-fDz; << 1345 if (Dist > safe) safe=Dist;
901 G4double dy1 = fPlanes[0].b*p.y()+fPlane <<
902 G4double dy2 = fPlanes[1].b*p.y()+fPlane <<
903 G4double dy = std::max(dz,std::max(dy1,d <<
904 <<
905 G4double dx1 = fPlanes[2].a*p.x()+fPlane <<
906 + fPlanes[2].c*p.z()+fPlane <<
907 G4double dx2 = fPlanes[3].a*p.x()+fPlane <<
908 + fPlanes[3].c*p.z()+fPlane <<
909 G4double dist = std::max(dy,std::max(dx1 <<
910 return (dist > 0) ? dist : 0.; <<
911 } <<
912 case 1: // YZ section is a rectangle <<
913 { <<
914 G4double dz = std::abs(p.z())-fDz; <<
915 G4double dy = std::max(dz,std::abs(p.y() <<
916 G4double dx1 = fPlanes[2].a*p.x()+fPlane <<
917 + fPlanes[2].c*p.z()+fPlane <<
918 G4double dx2 = fPlanes[3].a*p.x()+fPlane <<
919 + fPlanes[3].c*p.z()+fPlane <<
920 G4double dist = std::max(dy,std::max(dx1 <<
921 return (dist > 0) ? dist : 0.; <<
922 } <<
923 case 2: // YZ section is a rectangle and <<
924 { // XZ section is an isosceles trap <<
925 G4double dz = std::abs(p.z())-fDz; <<
926 G4double dy = std::max(dz,std::abs(p.y() <<
927 G4double dx = fPlanes[3].a*std::abs(p.x( <<
928 + fPlanes[3].c*p.z()+fPlanes <<
929 G4double dist = std::max(dy,dx); <<
930 return (dist > 0) ? dist : 0.; <<
931 } <<
932 case 3: // YZ section is a rectangle and <<
933 { // XY section is an isosceles trap <<
934 G4double dz = std::abs(p.z())-fDz; <<
935 G4double dy = std::max(dz,std::abs(p.y() <<
936 G4double dx = fPlanes[3].a*std::abs(p.x( <<
937 + fPlanes[3].b*p.y()+fPlanes <<
938 G4double dist = std::max(dy,dx); <<
939 return (dist > 0) ? dist : 0.; <<
940 } <<
941 } 1346 }
942 return 0.; << 1347 if (safe<0) safe=0;
>> 1348 return safe;
943 } 1349 }
944 1350
945 ////////////////////////////////////////////// << 1351 /////////////////////////////////////////////////////////////////////////////////
946 // 1352 //
947 // Calculate distance to surface of shape from << 1353 // Calculate distance to surface of shape from inside
948 // find normal at exit point, if required << 1354 // Calculate distance to x/y/z planes - smallest is exiting distance
949 // - when leaving the surface, return 0 <<
950 1355
951 G4double G4Trap::DistanceToOut(const G4ThreeVe 1356 G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
952 const G4bool ca 1357 const G4bool calcNorm,
953 G4bool* v << 1358 G4bool *validNorm, G4ThreeVector *n) const
954 { 1359 {
955 // Z intersections << 1360 Eside side = kUndef;
>> 1361 G4double snxt; // snxt = return value
>> 1362 G4double pdist,Comp,vdist,max;
>> 1363 //
>> 1364 // Z Intersections
956 // 1365 //
957 if ((std::abs(p.z()) - fDz) >= -halfCarToler << 1366 if (v.z()>0)
>> 1367 {
>> 1368 max=fDz-p.z();
>> 1369 if (max>kCarTolerance/2)
>> 1370 {
>> 1371 snxt=max/v.z();
>> 1372 side=kPZ;
>> 1373 }
>> 1374 else
>> 1375 {
>> 1376 if (calcNorm)
>> 1377 {
>> 1378 *validNorm=true;
>> 1379 *n=G4ThreeVector(0,0,1);
>> 1380 }
>> 1381 return snxt=0;
>> 1382 }
>> 1383 }
>> 1384 else if (v.z()<0)
958 { 1385 {
959 if (calcNorm) << 1386 max=-fDz-p.z();
>> 1387 if (max<-kCarTolerance/2)
>> 1388 {
>> 1389 snxt=max/v.z();
>> 1390 side=kMZ;
>> 1391 }
>> 1392 else
960 { 1393 {
961 *validNorm = true; << 1394 if (calcNorm)
962 n->set(0, 0, (p.z() < 0) ? -1 : 1); << 1395 {
>> 1396 *validNorm=true;
>> 1397 *n=G4ThreeVector(0,0,-1);
>> 1398 }
>> 1399 return snxt=0;
963 } 1400 }
964 return 0; <<
965 } 1401 }
966 G4double vz = v.z(); << 1402 else
967 G4double tmax = (vz == 0) ? DBL_MAX : (std:: << 1403 {
968 G4int iside = (vz < 0) ? -4 : -2; // little << 1404 snxt=kInfinity;
>> 1405 }
969 1406
970 // Y intersections <<
971 // 1407 //
972 G4int i = 0; << 1408 // Intersections with planes[0] (expanded because of setting enum)
973 for ( ; i<2; ++i) << 1409 //
>> 1410 pdist=fPlanes[0].a*p.x()+fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
>> 1411 Comp=fPlanes[0].a*v.x()+fPlanes[0].b*v.y()+fPlanes[0].c*v.z();
>> 1412 if (pdist>0)
974 { 1413 {
975 G4double cosa = fPlanes[i].b*v.y() + fPlan << 1414 // Outside the plane
976 if (cosa > 0) << 1415 if (Comp>0)
977 { 1416 {
978 G4double dist = fPlanes[i].b*p.y() + fPl << 1417 // Leaving immediately
979 if (dist >= -halfCarTolerance) << 1418 if (calcNorm)
980 { 1419 {
981 if (calcNorm) << 1420 *validNorm=true;
982 { << 1421 *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
983 *validNorm = true; <<
984 n->set(0, fPlanes[i].b, fPlanes[i].c <<
985 } <<
986 return 0; <<
987 } 1422 }
988 G4double tmp = -dist/cosa; << 1423 return snxt=0;
989 if (tmax > tmp) { tmax = tmp; iside = i; << 1424 }
>> 1425 }
>> 1426 else if (pdist<-kCarTolerance/2)
>> 1427 {
>> 1428 // Inside the plane
>> 1429 if (Comp>0)
>> 1430 {
>> 1431 // Will leave extent
>> 1432 vdist=-pdist/Comp;
>> 1433 if (vdist<snxt)
>> 1434 {
>> 1435 snxt=vdist;
>> 1436 side=ks0;
>> 1437 }
>> 1438 }
>> 1439 }
>> 1440 else
>> 1441 {
>> 1442 // On surface
>> 1443 if (Comp>0)
>> 1444 {
>> 1445 if (calcNorm)
>> 1446 {
>> 1447 *validNorm=true;
>> 1448 *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
>> 1449 }
>> 1450 return snxt=0;
990 } 1451 }
991 } 1452 }
992 1453
993 // X intersections <<
994 // 1454 //
995 for ( ; i<4; ++i) << 1455 // Intersections with planes[1] (expanded because of setting enum)
>> 1456 //
>> 1457 pdist=fPlanes[1].a*p.x()+fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
>> 1458 Comp=fPlanes[1].a*v.x()+fPlanes[1].b*v.y()+fPlanes[1].c*v.z();
>> 1459 if (pdist>0)
996 { 1460 {
997 G4double cosa = fPlanes[i].a*v.x()+fPlanes << 1461 // Outside the plane
998 if (cosa > 0) << 1462 if (Comp>0)
999 { 1463 {
1000 G4double dist = fPlanes[i].a*p.x() + << 1464 // Leaving immediately
1001 fPlanes[i].b*p.y() + fP << 1465 if (calcNorm)
1002 if (dist >= -halfCarTolerance) <<
1003 { 1466 {
1004 if (calcNorm) << 1467 *validNorm=true;
1005 { << 1468 *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1006 *validNorm = true; << 1469 }
1007 n->set(fPlanes[i].a, fPlanes[i].b, << 1470 return snxt=0;
1008 } << 1471 }
1009 return 0; << 1472 }
>> 1473 else if (pdist<-kCarTolerance/2)
>> 1474 {
>> 1475 // Inside the plane
>> 1476 if (Comp>0)
>> 1477 {
>> 1478 // Will leave extent
>> 1479 vdist=-pdist/Comp;
>> 1480 if (vdist<snxt)
>> 1481 {
>> 1482 snxt=vdist;
>> 1483 side=ks1;
1010 } 1484 }
1011 G4double tmp = -dist/cosa; << 1485 }
1012 if (tmax > tmp) { tmax = tmp; iside = i << 1486 }
>> 1487 else
>> 1488 {
>> 1489 // On surface
>> 1490 if (Comp>0)
>> 1491 {
>> 1492 if (calcNorm)
>> 1493 {
>> 1494 *validNorm=true;
>> 1495 *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
>> 1496 }
>> 1497 return snxt=0;
>> 1498 }
>> 1499 }
>> 1500
>> 1501 //
>> 1502 // Intersections with planes[2] (expanded because of setting enum)
>> 1503 //
>> 1504 pdist=fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d;
>> 1505 Comp=fPlanes[2].a*v.x()+fPlanes[2].b*v.y()+fPlanes[2].c*v.z();
>> 1506 if (pdist>0)
>> 1507 {
>> 1508 // Outside the plane
>> 1509 if (Comp>0)
>> 1510 {
>> 1511 // Leaving immediately
>> 1512 if (calcNorm)
>> 1513 {
>> 1514 *validNorm=true;
>> 1515 *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
>> 1516 }
>> 1517 return snxt=0;
>> 1518 }
>> 1519 }
>> 1520 else if (pdist<-kCarTolerance/2)
>> 1521 {
>> 1522 // Inside the plane
>> 1523 if (Comp>0)
>> 1524 {
>> 1525 // Will leave extent
>> 1526 vdist=-pdist/Comp;
>> 1527 if (vdist<snxt)
>> 1528 {
>> 1529 snxt=vdist;
>> 1530 side=ks2;
>> 1531 }
>> 1532 }
>> 1533 }
>> 1534 else
>> 1535 {
>> 1536 // On surface
>> 1537 if (Comp>0)
>> 1538 {
>> 1539 if (calcNorm)
>> 1540 {
>> 1541 *validNorm=true;
>> 1542 *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
>> 1543 }
>> 1544 return snxt=0;
1013 } 1545 }
1014 } 1546 }
1015 1547
1016 // Set normal, if required, and return dist <<
1017 // 1548 //
>> 1549 // Intersections with planes[3] (expanded because of setting enum)
>> 1550 //
>> 1551 pdist=fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d;
>> 1552 Comp=fPlanes[3].a*v.x()+fPlanes[3].b*v.y()+fPlanes[3].c*v.z();
>> 1553 if (pdist>0)
>> 1554 {
>> 1555 // Outside the plane
>> 1556 if (Comp>0)
>> 1557 {
>> 1558 // Leaving immediately
>> 1559 if (calcNorm)
>> 1560 {
>> 1561 *validNorm=true;
>> 1562 *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
>> 1563 }
>> 1564 return snxt=0;
>> 1565 }
>> 1566 }
>> 1567 else if (pdist<-kCarTolerance/2)
>> 1568 {
>> 1569 // Inside the plane
>> 1570 if (Comp>0)
>> 1571 {
>> 1572 // Will leave extent
>> 1573 vdist=-pdist/Comp;
>> 1574 if (vdist<snxt)
>> 1575 {
>> 1576 snxt=vdist;
>> 1577 side=ks3;
>> 1578 }
>> 1579 }
>> 1580 }
>> 1581 else
>> 1582 {
>> 1583 // On surface
>> 1584 if (Comp>0)
>> 1585 {
>> 1586 if (calcNorm)
>> 1587 {
>> 1588 *validNorm=true;
>> 1589 *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
>> 1590 }
>> 1591 return snxt=0;
>> 1592 }
>> 1593 }
>> 1594
>> 1595 // set normal
1018 if (calcNorm) 1596 if (calcNorm)
1019 { 1597 {
1020 *validNorm = true; << 1598 *validNorm=true;
1021 if (iside < 0) << 1599 switch(side)
1022 n->set(0, 0, iside + 3); // (-4+3)=-1, << 1600 {
1023 else << 1601 case ks0:
1024 n->set(fPlanes[iside].a, fPlanes[iside] << 1602 *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
>> 1603 break;
>> 1604 case ks1:
>> 1605 *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
>> 1606 break;
>> 1607 case ks2:
>> 1608 *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
>> 1609 break;
>> 1610 case ks3:
>> 1611 *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
>> 1612 break;
>> 1613 case kMZ:
>> 1614 *n=G4ThreeVector(0,0,-1);
>> 1615 break;
>> 1616 case kPZ:
>> 1617 *n=G4ThreeVector(0,0,1);
>> 1618 break;
>> 1619 default:
>> 1620 G4cout.precision(16);
>> 1621 G4cout << G4endl;
>> 1622 DumpInfo();
>> 1623 G4cout << "Position:" << G4endl << G4endl;
>> 1624 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
>> 1625 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
>> 1626 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
>> 1627 G4cout << "Direction:" << G4endl << G4endl;
>> 1628 G4cout << "v.x() = " << v.x() << G4endl;
>> 1629 G4cout << "v.y() = " << v.y() << G4endl;
>> 1630 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
>> 1631 G4cout << "Proposed distance :" << G4endl << G4endl;
>> 1632 G4cout << "snxt = " << snxt/mm << " mm" << G4endl << G4endl;
>> 1633 G4Exception("G4Trap::DistanceToOut(p,v,..)","Notification",JustWarning,
>> 1634 "Undefined side for valid surface normal to solid.");
>> 1635 break;
>> 1636 }
1025 } 1637 }
1026 return tmax; << 1638 return snxt;
1027 } 1639 }
1028 1640
1029 ///////////////////////////////////////////// << 1641 //////////////////////////////////////////////////////////////////////////////
1030 // 1642 //
1031 // Calculate exact shortest distance to any b 1643 // Calculate exact shortest distance to any boundary from inside
1032 // - Returns 0 is ThreeVector outside 1644 // - Returns 0 is ThreeVector outside
1033 1645
1034 G4double G4Trap::DistanceToOut( const G4Three 1646 G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const
1035 { 1647 {
>> 1648 G4double safe=0.0,Dist;
>> 1649 G4int i;
>> 1650
1036 #ifdef G4CSGDEBUG 1651 #ifdef G4CSGDEBUG
1037 if( Inside(p) == kOutside ) 1652 if( Inside(p) == kOutside )
1038 { 1653 {
1039 std::ostringstream message; << 1654 G4cout.precision(16) ;
1040 G4long oldprc = message.precision(16); << 1655 G4cout << G4endl ;
1041 message << "Point p is outside (!?) of so << 1656 DumpInfo();
1042 message << "Position:\n"; << 1657 G4cout << "Position:" << G4endl << G4endl ;
1043 message << " p.x() = " << p.x()/mm << " << 1658 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1044 message << " p.y() = " << p.y()/mm << " << 1659 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1045 message << " p.z() = " << p.z()/mm << " << 1660 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1046 G4cout.precision(oldprc); << 1661 G4Exception("G4Trap::DistanceToOut(p)",
1047 G4Exception("G4Trap::DistanceToOut(p)", " << 1662 "Notification", JustWarning, "Point p is outside !?" );
1048 JustWarning, message ); <<
1049 DumpInfo(); <<
1050 } 1663 }
1051 #endif 1664 #endif
1052 switch (fTrapType) << 1665
>> 1666 safe=fDz-std::fabs(p.z());
>> 1667 if (safe<0) safe=0;
>> 1668 else
1053 { 1669 {
1054 case 0: // General case << 1670 for (i=0;i<4;i++)
1055 { 1671 {
1056 G4double dz = std::abs(p.z())-fDz; << 1672 Dist=-(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1057 G4double dy1 = fPlanes[0].b*p.y()+fPlan << 1673 +fPlanes[i].c*p.z()+fPlanes[i].d);
1058 G4double dy2 = fPlanes[1].b*p.y()+fPlan << 1674 if (Dist<safe) safe=Dist;
1059 G4double dy = std::max(dz,std::max(dy1, << 1675 }
1060 << 1676 if (safe<0) safe=0;
1061 G4double dx1 = fPlanes[2].a*p.x()+fPlan << 1677 }
1062 + fPlanes[2].c*p.z()+fPlan << 1678 return safe;
1063 G4double dx2 = fPlanes[3].a*p.x()+fPlan <<
1064 + fPlanes[3].c*p.z()+fPlan <<
1065 G4double dist = std::max(dy,std::max(dx <<
1066 return (dist < 0) ? -dist : 0.; <<
1067 } <<
1068 case 1: // YZ section is a rectangle <<
1069 { <<
1070 G4double dz = std::abs(p.z())-fDz; <<
1071 G4double dy = std::max(dz,std::abs(p.y( <<
1072 G4double dx1 = fPlanes[2].a*p.x()+fPlan <<
1073 + fPlanes[2].c*p.z()+fPlan <<
1074 G4double dx2 = fPlanes[3].a*p.x()+fPlan <<
1075 + fPlanes[3].c*p.z()+fPlan <<
1076 G4double dist = std::max(dy,std::max(dx <<
1077 return (dist < 0) ? -dist : 0.; <<
1078 } <<
1079 case 2: // YZ section is a rectangle and <<
1080 { // XZ section is an isosceles tra <<
1081 G4double dz = std::abs(p.z())-fDz; <<
1082 G4double dy = std::max(dz,std::abs(p.y( <<
1083 G4double dx = fPlanes[3].a*std::abs(p.x <<
1084 + fPlanes[3].c*p.z()+fPlane <<
1085 G4double dist = std::max(dy,dx); <<
1086 return (dist < 0) ? -dist : 0.; <<
1087 } <<
1088 case 3: // YZ section is a rectangle and <<
1089 { // XY section is an isosceles tra <<
1090 G4double dz = std::abs(p.z())-fDz; <<
1091 G4double dy = std::max(dz,std::abs(p.y( <<
1092 G4double dx = fPlanes[3].a*std::abs(p.x <<
1093 + fPlanes[3].b*p.y()+fPlane <<
1094 G4double dist = std::max(dy,dx); <<
1095 return (dist < 0) ? -dist : 0.; <<
1096 } <<
1097 } <<
1098 return 0.; <<
1099 } 1679 }
1100 1680
1101 ///////////////////////////////////////////// 1681 //////////////////////////////////////////////////////////////////////////
1102 // 1682 //
1103 // GetEntityType << 1683 // Create a List containing the transformed vertices
1104 << 1684 // Ordering [0-3] -fDz cross section
1105 G4GeometryType G4Trap::GetEntityType() const << 1685 // [4-7] +fDz cross section such that [0] is below [4],
1106 { << 1686 // [1] below [5] etc.
1107 return {"G4Trap"}; << 1687 // Note:
>> 1688 // Caller has deletion resposibility
>> 1689
>> 1690 G4ThreeVectorList*
>> 1691 G4Trap::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
>> 1692 {
>> 1693 G4ThreeVectorList *vertices;
>> 1694 vertices=new G4ThreeVectorList();
>> 1695 vertices->reserve(8);
>> 1696 if (vertices)
>> 1697 {
>> 1698 G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
>> 1699 -fDz*fTthetaSphi-fDy1,-fDz);
>> 1700 G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
>> 1701 -fDz*fTthetaSphi-fDy1,-fDz);
>> 1702 G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
>> 1703 -fDz*fTthetaSphi+fDy1,-fDz);
>> 1704 G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
>> 1705 -fDz*fTthetaSphi+fDy1,-fDz);
>> 1706 G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
>> 1707 +fDz*fTthetaSphi-fDy2,+fDz);
>> 1708 G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
>> 1709 +fDz*fTthetaSphi-fDy2,+fDz);
>> 1710 G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
>> 1711 +fDz*fTthetaSphi+fDy2,+fDz);
>> 1712 G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
>> 1713 +fDz*fTthetaSphi+fDy2,+fDz);
>> 1714
>> 1715 vertices->push_back(pTransform.TransformPoint(vertex0));
>> 1716 vertices->push_back(pTransform.TransformPoint(vertex1));
>> 1717 vertices->push_back(pTransform.TransformPoint(vertex2));
>> 1718 vertices->push_back(pTransform.TransformPoint(vertex3));
>> 1719 vertices->push_back(pTransform.TransformPoint(vertex4));
>> 1720 vertices->push_back(pTransform.TransformPoint(vertex5));
>> 1721 vertices->push_back(pTransform.TransformPoint(vertex6));
>> 1722 vertices->push_back(pTransform.TransformPoint(vertex7));
>> 1723 }
>> 1724 else
>> 1725 {
>> 1726 DumpInfo();
>> 1727 G4Exception("G4Trap::CreateRotatedVertices()",
>> 1728 "FatalError", FatalException,
>> 1729 "Error in allocation of vertices. Out of memory !");
>> 1730 }
>> 1731 return vertices;
1108 } 1732 }
1109 1733
1110 ///////////////////////////////////////////// 1734 //////////////////////////////////////////////////////////////////////////
1111 // 1735 //
1112 // IsFaceted << 1736 // GetEntityType
1113 <<
1114 G4bool G4Trap::IsFaceted() const <<
1115 { <<
1116 return true; <<
1117 } <<
1118 1737
1119 ///////////////////////////////////////////// << 1738 G4GeometryType G4Trap::GetEntityType() const
1120 // <<
1121 // Make a clone of the object <<
1122 // <<
1123 G4VSolid* G4Trap::Clone() const <<
1124 { 1739 {
1125 return new G4Trap(*this); << 1740 return G4String("G4Trap");
1126 } 1741 }
1127 1742
1128 ///////////////////////////////////////////// 1743 //////////////////////////////////////////////////////////////////////////
1129 // 1744 //
1130 // Stream object contents to an output stream 1745 // Stream object contents to an output stream
1131 1746
1132 std::ostream& G4Trap::StreamInfo( std::ostrea 1747 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1133 { 1748 {
1134 G4double phi = GetPhi(); <<
1135 G4double theta = GetTheta(); <<
1136 G4double alpha1 = GetAlpha1(); <<
1137 G4double alpha2 = GetAlpha2(); <<
1138 <<
1139 G4long oldprc = os.precision(16); <<
1140 os << "------------------------------------ 1749 os << "-----------------------------------------------------------\n"
1141 << " *** Dump for solid: " << GetName << 1750 << " *** Dump for solid - " << GetName() << " ***\n"
1142 << " ================================ 1751 << " ===================================================\n"
1143 << " Solid type: G4Trap\n" 1752 << " Solid type: G4Trap\n"
1144 << " Parameters:\n" << 1753 << " Parameters: \n"
1145 << " half length Z: " << fDz/mm << " << 1754 << " half length Z: " << fDz/mm << " mm \n"
1146 << " half length Y, face -Dz: " << fD << 1755 << " half length Y of face -fDz: " << fDy1/mm << " mm \n"
1147 << " half length X, face -Dz, side -D << 1756 << " half length X of side -fDy1, face -fDz: " << fDx1/mm << " mm \n"
1148 << " half length X, face -Dz, side +D << 1757 << " half length X of side +fDy1, face -fDz: " << fDx2/mm << " mm \n"
1149 << " half length Y, face +Dz: " << fD << 1758 << " half length Y of face +fDz: " << fDy2/mm << " mm \n"
1150 << " half length X, face +Dz, side -D << 1759 << " half length X of side -fDy2, face +fDz: " << fDx3/mm << " mm \n"
1151 << " half length X, face +Dz, side +D << 1760 << " half length X of side +fDy2, face +fDz: " << fDx4/mm << " mm \n"
1152 << " theta: " << theta/degree << " de << 1761 << " std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree << " degrees \n"
1153 << " phi: " << phi/degree << " degr << 1762 << " std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree << " degrees \n"
1154 << " alpha, face -Dz: " << alpha1/deg << 1763 << " std::tan(alpha), -fDz: " << fTalpha1/degree << " degrees \n"
1155 << " alpha, face +Dz: " << alpha2/deg << 1764 << " std::tan(alpha), +fDz: " << fTalpha2/degree << " degrees \n"
>> 1765 << " trap side plane equations:\n"
>> 1766 << " " << fPlanes[0].a << " X + " << fPlanes[0].b << " Y + "
>> 1767 << fPlanes[0].c << " Z + " << fPlanes[0].d << " = 0\n"
>> 1768 << " " << fPlanes[1].a << " X + " << fPlanes[1].b << " Y + "
>> 1769 << fPlanes[1].c << " Z + " << fPlanes[1].d << " = 0\n"
>> 1770 << " " << fPlanes[2].a << " X + " << fPlanes[2].b << " Y + "
>> 1771 << fPlanes[2].c << " Z + " << fPlanes[2].d << " = 0\n"
>> 1772 << " " << fPlanes[3].a << " X + " << fPlanes[3].b << " Y + "
>> 1773 << fPlanes[3].c << " Z + " << fPlanes[3].d << " = 0\n"
1156 << "------------------------------------ 1774 << "-----------------------------------------------------------\n";
1157 os.precision(oldprc); <<
1158 1775
1159 return os; 1776 return os;
1160 } 1777 }
1161 1778
1162 ///////////////////////////////////////////// << 1779 /////////////////////////////////////////////////////////////////////////
1163 // 1780 //
1164 // Compute vertices from planes << 1781 // GetPointOnPlane
>> 1782 //
>> 1783 // Auxiliary method for Get Point on Surface
1165 1784
1166 void G4Trap::GetVertices(G4ThreeVector pt[8]) << 1785 G4ThreeVector G4Trap::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1,
1167 { << 1786 G4ThreeVector p2, G4ThreeVector p3,
1168 for (G4int i=0; i<8; ++i) << 1787 G4double& area) const
1169 { << 1788 {
1170 G4int iy = (i==0 || i==1 || i==4 || i==5) << 1789 G4double lambda1, lambda2, chose, aOne, aTwo;
1171 G4int ix = (i==0 || i==2 || i==4 || i==6) << 1790 G4ThreeVector t, u, v, w, Area, normal;
1172 G4double z = (i < 4) ? -fDz : fDz; << 1791
1173 G4double y = -(fPlanes[iy].c*z + fPlanes[ << 1792 t = p1 - p0;
1174 G4double x = -(fPlanes[ix].b*y + fPlanes[ << 1793 u = p2 - p1;
1175 + fPlanes[ix].d)/fPlanes[i << 1794 v = p3 - p2;
1176 pt[i].set(x,y,z); << 1795 w = p0 - p3;
1177 } << 1796
>> 1797 Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
>> 1798 w.z()*v.x() - w.x()*v.z(),
>> 1799 w.x()*v.y() - w.y()*v.x());
>> 1800
>> 1801 aOne = 0.5*Area.mag();
>> 1802
>> 1803 Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
>> 1804 t.z()*u.x() - t.x()*u.z(),
>> 1805 t.x()*u.y() - t.y()*u.x());
>> 1806
>> 1807 aTwo = 0.5*Area.mag();
>> 1808
>> 1809 area = aOne + aTwo;
>> 1810
>> 1811 chose = RandFlat::shoot(0.,aOne+aTwo);
>> 1812
>> 1813 if( (chose>=0.) && (chose < aOne) )
>> 1814 {
>> 1815 lambda1 = RandFlat::shoot(0.,1.);
>> 1816 lambda2 = RandFlat::shoot(0.,lambda1);
>> 1817 return (p2+lambda1*v+lambda2*w);
>> 1818 }
>> 1819
>> 1820 // else
>> 1821
>> 1822 lambda1 = RandFlat::shoot(0.,1.);
>> 1823 lambda2 = RandFlat::shoot(0.,lambda1);
>> 1824
>> 1825 return (p0+lambda1*t+lambda2*u);
1178 } 1826 }
1179 1827
1180 ///////////////////////////////////////////// << 1828 ///////////////////////////////////////////////////////////////
1181 // 1829 //
1182 // Generate random point on the surface << 1830 // GetPointOnSurface
1183 1831
1184 G4ThreeVector G4Trap::GetPointOnSurface() con 1832 G4ThreeVector G4Trap::GetPointOnSurface() const
1185 { 1833 {
1186 // Set indeces << 1834 G4double aOne, aTwo, aThree, aFour, aFive, aSix, chose;
1187 constexpr G4int iface [6][4] = << 1835 G4ThreeVector One, Two, Three, Four, Five, Six, test;
1188 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6 <<
1189 <<
1190 // Set vertices <<
1191 G4ThreeVector pt[8]; 1836 G4ThreeVector pt[8];
1192 GetVertices(pt); << 1837
1193 << 1838 pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
1194 // Select face << 1839 -fDz*fTthetaSphi-fDy1,-fDz);
1195 // << 1840 pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
1196 G4double select = fAreas[5]*G4QuickRand(); << 1841 -fDz*fTthetaSphi-fDy1,-fDz);
1197 G4int k = 5; << 1842 pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
1198 k -= (G4int)(select <= fAreas[4]); << 1843 -fDz*fTthetaSphi+fDy1,-fDz);
1199 k -= (G4int)(select <= fAreas[3]); << 1844 pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
1200 k -= (G4int)(select <= fAreas[2]); << 1845 -fDz*fTthetaSphi+fDy1,-fDz);
1201 k -= (G4int)(select <= fAreas[1]); << 1846 pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
1202 k -= (G4int)(select <= fAreas[0]); << 1847 +fDz*fTthetaSphi-fDy2,+fDz);
1203 << 1848 pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
1204 // Select sub-triangle << 1849 +fDz*fTthetaSphi-fDy2,+fDz);
1205 // << 1850 pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
1206 G4int i0 = iface[k][0]; << 1851 +fDz*fTthetaSphi+fDy2,+fDz);
1207 G4int i1 = iface[k][1]; << 1852 pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
1208 G4int i2 = iface[k][2]; << 1853 +fDz*fTthetaSphi+fDy2,+fDz);
1209 G4int i3 = iface[k][3]; << 1854
1210 G4double s2 = G4GeomTools::TriangleAreaNorm << 1855 // make sure we provide the points in a clockwise fashion
1211 if (select > fAreas[k] - s2) i0 = i2; << 1856
1212 << 1857 One = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
1213 // Generate point << 1858 Two = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
1214 // << 1859 Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
1215 G4double u = G4QuickRand(); << 1860 Four = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour);
1216 G4double v = G4QuickRand(); << 1861 Five = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
1217 if (u + v > 1.) { u = 1. - u; v = 1. - v; } << 1862 Six = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
1218 return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3 << 1863
>> 1864 chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
>> 1865 if( (chose>=0.) && (chose<aOne) )
>> 1866 { return One; }
>> 1867 else if( (chose>=aOne) && (chose<aOne+aTwo) )
>> 1868 { return Two; }
>> 1869 else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) )
>> 1870 { return Three; }
>> 1871 else if( (chose>=aOne+aTwo+aThree) && (chose<aOne+aTwo+aThree+aFour) )
>> 1872 { return Four; }
>> 1873 else if( (chose>=aOne+aTwo+aThree+aFour)
>> 1874 && (chose<aOne+aTwo+aThree+aFour+aFive) )
>> 1875 { return Five; }
>> 1876 return Six;
1219 } 1877 }
1220 1878
1221 ///////////////////////////////////////////// 1879 //////////////////////////////////////////////////////////////////////////
1222 // 1880 //
1223 // Methods for visualisation 1881 // Methods for visualisation
1224 1882
1225 void G4Trap::DescribeYourselfTo ( G4VGraphics 1883 void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1226 { 1884 {
1227 scene.AddSolid (*this); 1885 scene.AddSolid (*this);
1228 } 1886 }
1229 1887
1230 G4Polyhedron* G4Trap::CreatePolyhedron () con 1888 G4Polyhedron* G4Trap::CreatePolyhedron () const
1231 { 1889 {
1232 G4double phi = std::atan2(fTthetaSphi, fTth 1890 G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1233 G4double alpha1 = std::atan(fTalpha1); 1891 G4double alpha1 = std::atan(fTalpha1);
1234 G4double alpha2 = std::atan(fTalpha2); 1892 G4double alpha2 = std::atan(fTalpha2);
1235 G4double theta = std::atan(std::sqrt(fTthet << 1893 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi+fTthetaSphi*fTthetaSphi));
1236 +fTthet <<
1237 1894
1238 return new G4PolyhedronTrap(fDz, theta, phi 1895 return new G4PolyhedronTrap(fDz, theta, phi,
1239 fDy1, fDx1, fDx 1896 fDy1, fDx1, fDx2, alpha1,
1240 fDy2, fDx3, fDx 1897 fDy2, fDx3, fDx4, alpha2);
1241 } 1898 }
1242 1899
1243 #endif << 1900 G4NURBS* G4Trap::CreateNURBS () const
>> 1901 {
>> 1902 // return new G4NURBSbox (fDx, fDy, fDz);
>> 1903 return 0 ;
>> 1904 }
1244 1905