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