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