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