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