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