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