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