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