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