| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Implementation for G4Trap class 26 // Implementation for G4Trap class
27 // 27 //
28 // 21.03.95 P.Kent: Modified for `tolerant' ge 28 // 21.03.95 P.Kent: Modified for `tolerant' geometry
29 // 09.09.96 V.Grichine: Final modifications be 29 // 09.09.96 V.Grichine: Final modifications before to commit
30 // 08.12.97 J.Allison: Added "nominal" constru 30 // 08.12.97 J.Allison: Added "nominal" constructor and method SetAllParameters
31 // 28.04.05 V.Grichine: new SurfaceNormal acco << 31 // 28.04.05 V.Grichine: new SurfaceNormal according to J.Apostolakis proposal
32 // 18.04.17 E.Tcherniaev: complete revision, s 32 // 18.04.17 E.Tcherniaev: complete revision, speed-up
33 // ------------------------------------------- 33 // --------------------------------------------------------------------
34 34
35 #include "G4Trap.hh" 35 #include "G4Trap.hh"
36 36
37 #if !defined(G4GEOM_USE_UTRAP) 37 #if !defined(G4GEOM_USE_UTRAP)
38 38
39 #include "globals.hh" 39 #include "globals.hh"
40 #include "G4GeomTools.hh" 40 #include "G4GeomTools.hh"
41 41
42 #include "G4VoxelLimits.hh" 42 #include "G4VoxelLimits.hh"
43 #include "G4AffineTransform.hh" 43 #include "G4AffineTransform.hh"
44 #include "G4BoundingEnvelope.hh" 44 #include "G4BoundingEnvelope.hh"
45 45
46 #include "G4VPVParameterisation.hh" 46 #include "G4VPVParameterisation.hh"
47 47
48 #include "G4QuickRand.hh" << 48 #include "Randomize.hh"
49 49
50 #include "G4VGraphicsScene.hh" 50 #include "G4VGraphicsScene.hh"
51 #include "G4Polyhedron.hh" 51 #include "G4Polyhedron.hh"
52 52
53 using namespace CLHEP; 53 using namespace CLHEP;
54 54
55 ////////////////////////////////////////////// 55 //////////////////////////////////////////////////////////////////////////
56 // 56 //
57 // Constructor - check and set half-widths as << 57 // Constructor - check and set half-widths as well as angles:
58 // final check of coplanarity 58 // final check of coplanarity
59 59
60 G4Trap::G4Trap( const G4String& pName, 60 G4Trap::G4Trap( const G4String& pName,
61 G4double pDz, 61 G4double pDz,
62 G4double pTheta, G4doubl 62 G4double pTheta, G4double pPhi,
63 G4double pDy1, G4double 63 G4double pDy1, G4double pDx1, G4double pDx2,
64 G4double pAlp1, 64 G4double pAlp1,
65 G4double pDy2, G4double 65 G4double pDy2, G4double pDx3, G4double pDx4,
66 G4double pAlp2 ) << 66 G4double pAlp2)
67 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 67 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
68 { 68 {
69 fDz = pDz; 69 fDz = pDz;
70 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi 70 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
71 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi 71 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
72 72
73 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp 73 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalpha1 = std::tan(pAlp1);
74 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp 74 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalpha2 = std::tan(pAlp2);
75 75
76 CheckParameters(); 76 CheckParameters();
77 MakePlanes(); 77 MakePlanes();
78 } 78 }
79 79
80 ////////////////////////////////////////////// 80 //////////////////////////////////////////////////////////////////////////
81 // 81 //
82 // Constructor - Design of trapezoid based on << 82 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
83 // which are its vertices. Checking of planari << 83 // which are its vertices. Checking of planarity with preparation of
84 // fPlanes[] and than calculation of other mem 84 // fPlanes[] and than calculation of other members
85 85
86 G4Trap::G4Trap( const G4String& pName, 86 G4Trap::G4Trap( const G4String& pName,
87 const G4ThreeVector pt[8] ) 87 const G4ThreeVector pt[8] )
88 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 88 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
89 { 89 {
90 // Start with check of centering - the cente 90 // Start with check of centering - the center of gravity trap line
91 // should cross the origin of frame 91 // should cross the origin of frame
92 // 92 //
93 if ( pt[0].z() >= 0 << 93 if (!( pt[0].z() < 0
94 || pt[0].z() != pt[1].z() << 94 && pt[0].z() == pt[1].z()
95 || pt[0].z() != pt[2].z() << 95 && pt[0].z() == pt[2].z()
96 || pt[0].z() != pt[3].z() << 96 && pt[0].z() == pt[3].z()
97 << 97
98 || pt[4].z() <= 0 << 98 && pt[4].z() > 0
99 || pt[4].z() != pt[5].z() << 99 && pt[4].z() == pt[5].z()
100 || pt[4].z() != pt[6].z() << 100 && pt[4].z() == pt[6].z()
101 || pt[4].z() != pt[7].z() << 101 && pt[4].z() == pt[7].z()
102 << 102
103 || std::fabs( pt[0].z() + pt[4].z() ) << 103 && std::fabs( pt[0].z() + pt[4].z() ) < kCarTolerance
104 << 104
105 || pt[0].y() != pt[1].y() << 105 && pt[0].y() == pt[1].y()
106 || pt[2].y() != pt[3].y() << 106 && pt[2].y() == pt[3].y()
107 || pt[4].y() != pt[5].y() << 107 && pt[4].y() == pt[5].y()
108 || pt[6].y() != pt[7].y() << 108 && pt[6].y() == pt[7].y()
109 << 109
110 || std::fabs(pt[0].y()+pt[2].y()+pt[4] << 110 && std::fabs(pt[0].y()+pt[2].y()+pt[4].y()+pt[6].y()) < kCarTolerance
111 || std::fabs(pt[0].x()+pt[1].x()+pt[4] << 111 && std::fabs(pt[0].x()+pt[1].x()+pt[4].x()+pt[5].x() +
112 pt[2].x()+pt[3].x()+pt[6] << 112 pt[2].x()+pt[3].x()+pt[6].x()+pt[7].x()) < kCarTolerance ))
113 { 113 {
114 std::ostringstream message; 114 std::ostringstream message;
115 message << "Invalid vertice coordinates fo 115 message << "Invalid vertice coordinates for Solid: " << GetName();
116 G4Exception("G4Trap::G4Trap()", "GeomSolid 116 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
117 FatalException, message); 117 FatalException, message);
118 } 118 }
119 << 119
120 // Set parameters 120 // Set parameters
121 // 121 //
122 fDz = (pt[7]).z(); 122 fDz = (pt[7]).z();
123 << 123
124 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5; 124 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5;
125 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5; 125 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5;
126 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5; 126 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5;
127 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]). 127 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1;
128 128
129 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5; 129 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5;
130 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5; 130 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5;
131 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5; 131 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5;
132 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]). 132 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2;
133 133
134 fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx 134 fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz;
135 fTthetaSphi = ((pt[4]).y()+fDy2)/fDz; 135 fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;
136 136
137 CheckParameters(); 137 CheckParameters();
138 MakePlanes(pt); 138 MakePlanes(pt);
139 } 139 }
140 140
141 ////////////////////////////////////////////// 141 //////////////////////////////////////////////////////////////////////////
142 // 142 //
143 // Constructor for Right Angular Wedge from ST 143 // Constructor for Right Angular Wedge from STEP
144 144
145 G4Trap::G4Trap( const G4String& pName, 145 G4Trap::G4Trap( const G4String& pName,
146 G4double pZ, 146 G4double pZ,
147 G4double pY, 147 G4double pY,
148 G4double pX, G4double pL 148 G4double pX, G4double pLTX )
149 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 149 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
150 { 150 {
151 fDz = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi 151 fDz = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi = 0;
152 fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLT 152 fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLTX; fTalpha1 = 0.5*(pLTX - pX)/pY;
153 fDy2 = fDy1; fDx3 = fDx1; fDx4 = fDx2; 153 fDy2 = fDy1; fDx3 = fDx1; fDx4 = fDx2; fTalpha2 = fTalpha1;
154 154
155 CheckParameters(); 155 CheckParameters();
156 MakePlanes(); 156 MakePlanes();
157 } 157 }
158 158
159 ////////////////////////////////////////////// 159 //////////////////////////////////////////////////////////////////////////
160 // 160 //
161 // Constructor for G4Trd 161 // Constructor for G4Trd
162 162
163 G4Trap::G4Trap( const G4String& pName, 163 G4Trap::G4Trap( const G4String& pName,
164 G4double pDx1, G4double 164 G4double pDx1, G4double pDx2,
165 G4double pDy1, G4double 165 G4double pDy1, G4double pDy2,
166 G4double pDz ) 166 G4double pDz )
167 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 167 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance), fTrapType(0)
168 { 168 {
169 fDz = pDz; fTthetaCphi = 0; fTthetaSphi = 169 fDz = pDz; fTthetaCphi = 0; fTthetaSphi = 0;
170 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalp 170 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalpha1 = 0;
171 fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalp 171 fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalpha2 = 0;
172 172
173 CheckParameters(); 173 CheckParameters();
174 MakePlanes(); 174 MakePlanes();
175 } 175 }
176 176
177 ////////////////////////////////////////////// 177 //////////////////////////////////////////////////////////////////////////
178 // 178 //
179 // Constructor for G4Para 179 // Constructor for G4Para
180 180
181 G4Trap::G4Trap( const G4String& pName, 181 G4Trap::G4Trap( const G4String& pName,
182 G4double pDx, G4double p 182 G4double pDx, G4double pDy,
183 G4double pDz, 183 G4double pDz,
184 G4double pAlpha, 184 G4double pAlpha,
185 G4double pTheta, G4doubl 185 G4double pTheta, G4double pPhi )
186 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 186 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
187 { 187 {
188 fDz = pDz; 188 fDz = pDz;
189 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi 189 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
190 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi 190 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
191 191
192 fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1 192 fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1 = std::tan(pAlpha);
193 fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2 193 fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2 = fTalpha1;
194 194
195 CheckParameters(); 195 CheckParameters();
196 MakePlanes(); 196 MakePlanes();
197 } 197 }
198 198
199 ////////////////////////////////////////////// 199 //////////////////////////////////////////////////////////////////////////
200 // 200 //
201 // Nominal constructor for G4Trap whose parame 201 // Nominal constructor for G4Trap whose parameters are to be set by
202 // a G4VParamaterisation later. Check and set 202 // a G4VParamaterisation later. Check and set half-widths as well as
203 // angles: final check of coplanarity 203 // angles: final check of coplanarity
204 204
205 G4Trap::G4Trap( const G4String& pName ) 205 G4Trap::G4Trap( const G4String& pName )
206 : G4CSGSolid (pName), halfCarTolerance(0.5*k 206 : G4CSGSolid (pName), halfCarTolerance(0.5*kCarTolerance),
207 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), 207 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
208 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.) 208 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
209 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.) 209 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
210 { 210 {
211 MakePlanes(); 211 MakePlanes();
212 } 212 }
213 213
214 ////////////////////////////////////////////// 214 //////////////////////////////////////////////////////////////////////////
215 // 215 //
216 // Fake default constructor - sets only member 216 // Fake default constructor - sets only member data and allocates memory
217 // for usage restri 217 // for usage restricted to object persistency.
218 // 218 //
219 G4Trap::G4Trap( __void__& a ) 219 G4Trap::G4Trap( __void__& a )
220 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo 220 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTolerance),
221 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), 221 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
222 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.) 222 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
223 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.) 223 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
224 { 224 {
225 MakePlanes(); 225 MakePlanes();
226 } 226 }
227 227
228 ////////////////////////////////////////////// 228 //////////////////////////////////////////////////////////////////////////
229 // 229 //
230 // Destructor 230 // Destructor
231 231
232 G4Trap::~G4Trap() = default; << 232 G4Trap::~G4Trap()
>> 233 {
>> 234 }
233 235
234 ////////////////////////////////////////////// 236 //////////////////////////////////////////////////////////////////////////
235 // 237 //
236 // Copy constructor 238 // Copy constructor
237 239
238 G4Trap::G4Trap(const G4Trap& rhs) 240 G4Trap::G4Trap(const G4Trap& rhs)
239 : G4CSGSolid(rhs), halfCarTolerance(rhs.half 241 : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance),
240 fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi) 242 fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
241 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.f 243 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
242 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.f 244 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
243 { 245 {
244 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs 246 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
245 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs. <<
246 fTrapType = rhs.fTrapType; 247 fTrapType = rhs.fTrapType;
247 } 248 }
248 249
249 ////////////////////////////////////////////// 250 //////////////////////////////////////////////////////////////////////////
250 // 251 //
251 // Assignment operator 252 // Assignment operator
252 253
253 G4Trap& G4Trap::operator = (const G4Trap& rhs) << 254 G4Trap& G4Trap::operator = (const G4Trap& rhs)
254 { 255 {
255 // Check assignment to self 256 // Check assignment to self
256 // 257 //
257 if (this == &rhs) { return *this; } 258 if (this == &rhs) { return *this; }
258 259
259 // Copy base class data 260 // Copy base class data
260 // 261 //
261 G4CSGSolid::operator=(rhs); 262 G4CSGSolid::operator=(rhs);
262 263
263 // Copy data 264 // Copy data
264 // 265 //
265 halfCarTolerance = rhs.halfCarTolerance; 266 halfCarTolerance = rhs.halfCarTolerance;
266 fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi 267 fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi; fTthetaSphi = rhs.fTthetaSphi;
267 fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs 268 fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
268 fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs 269 fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
269 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs 270 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
270 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs. <<
271 fTrapType = rhs.fTrapType; 271 fTrapType = rhs.fTrapType;
272 return *this; 272 return *this;
273 } 273 }
274 274
275 ////////////////////////////////////////////// 275 //////////////////////////////////////////////////////////////////////////
276 // 276 //
277 // Set all parameters, as for constructor - ch 277 // Set all parameters, as for constructor - check and set half-widths
278 // as well as angles: final check of coplanari 278 // as well as angles: final check of coplanarity
279 279
280 void G4Trap::SetAllParameters ( G4double pDz, 280 void G4Trap::SetAllParameters ( G4double pDz,
281 G4double pThet 281 G4double pTheta,
282 G4double pPhi, 282 G4double pPhi,
283 G4double pDy1, 283 G4double pDy1,
284 G4double pDx1, 284 G4double pDx1,
285 G4double pDx2, 285 G4double pDx2,
286 G4double pAlp1 286 G4double pAlp1,
287 G4double pDy2, 287 G4double pDy2,
288 G4double pDx3, 288 G4double pDx3,
289 G4double pDx4, 289 G4double pDx4,
290 G4double pAlp2 290 G4double pAlp2 )
291 { 291 {
292 // Reset data of the base class 292 // Reset data of the base class
293 fCubicVolume = 0; 293 fCubicVolume = 0;
294 fSurfaceArea = 0; 294 fSurfaceArea = 0;
295 fRebuildPolyhedron = true; 295 fRebuildPolyhedron = true;
296 296
297 // Set parameters 297 // Set parameters
298 fDz = pDz; 298 fDz = pDz;
299 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi 299 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
300 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi 300 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
301 301
302 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp 302 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalpha1 = std::tan(pAlp1);
303 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp 303 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalpha2 = std::tan(pAlp2);
304 304
305 CheckParameters(); 305 CheckParameters();
306 MakePlanes(); 306 MakePlanes();
307 } 307 }
308 308
309 ////////////////////////////////////////////// 309 //////////////////////////////////////////////////////////////////////////
310 // 310 //
311 // Check length parameters 311 // Check length parameters
312 312
313 void G4Trap::CheckParameters() 313 void G4Trap::CheckParameters()
314 { 314 {
315 if (fDz<=0 || 315 if (fDz<=0 ||
316 fDy1<=0 || fDx1<=0 || fDx2<=0 || 316 fDy1<=0 || fDx1<=0 || fDx2<=0 ||
317 fDy2<=0 || fDx3<=0 || fDx4<=0) 317 fDy2<=0 || fDx3<=0 || fDx4<=0)
318 { 318 {
319 std::ostringstream message; 319 std::ostringstream message;
320 message << "Invalid Length Parameters for 320 message << "Invalid Length Parameters for Solid: " << GetName()
321 << "\n X - " <<fDx1<<", "<<fDx2<< 321 << "\n X - " <<fDx1<<", "<<fDx2<<", "<<fDx3<<", "<<fDx4
322 << "\n Y - " <<fDy1<<", "<<fDy2 322 << "\n Y - " <<fDy1<<", "<<fDy2
323 << "\n Z - " <<fDz; 323 << "\n Z - " <<fDz;
324 G4Exception("G4Trap::CheckParameters()", " 324 G4Exception("G4Trap::CheckParameters()", "GeomSolids0002",
325 FatalException, message); 325 FatalException, message);
326 } 326 }
327 } 327 }
328 328
329 ////////////////////////////////////////////// 329 //////////////////////////////////////////////////////////////////////////
330 // 330 //
331 // Compute vertices and set side planes 331 // Compute vertices and set side planes
332 332
333 void G4Trap::MakePlanes() 333 void G4Trap::MakePlanes()
334 { 334 {
335 G4double DzTthetaCphi = fDz*fTthetaCphi; 335 G4double DzTthetaCphi = fDz*fTthetaCphi;
336 G4double DzTthetaSphi = fDz*fTthetaSphi; 336 G4double DzTthetaSphi = fDz*fTthetaSphi;
337 G4double Dy1Talpha1 = fDy1*fTalpha1; 337 G4double Dy1Talpha1 = fDy1*fTalpha1;
338 G4double Dy2Talpha2 = fDy2*fTalpha2; 338 G4double Dy2Talpha2 = fDy2*fTalpha2;
339 339
340 G4ThreeVector pt[8] = 340 G4ThreeVector pt[8] =
341 { 341 {
342 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx 342 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx1,-DzTthetaSphi-fDy1,-fDz),
343 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx 343 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx1,-DzTthetaSphi-fDy1,-fDz),
344 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx 344 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx2,-DzTthetaSphi+fDy1,-fDz),
345 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx 345 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx2,-DzTthetaSphi+fDy1,-fDz),
346 G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx 346 G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx3, DzTthetaSphi-fDy2, fDz),
347 G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx 347 G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx3, DzTthetaSphi-fDy2, fDz),
348 G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx 348 G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx4, DzTthetaSphi+fDy2, fDz),
349 G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx 349 G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx4, DzTthetaSphi+fDy2, fDz)
350 }; 350 };
351 351
352 MakePlanes(pt); 352 MakePlanes(pt);
353 } 353 }
354 354
355 ////////////////////////////////////////////// 355 //////////////////////////////////////////////////////////////////////////
356 // 356 //
357 // Set side planes, check planarity 357 // Set side planes, check planarity
358 358
359 void G4Trap::MakePlanes(const G4ThreeVector pt 359 void G4Trap::MakePlanes(const G4ThreeVector pt[8])
360 { 360 {
361 constexpr G4int iface[4][4] = { {0,4,5,1}, { << 361 G4int iface[4][4] = { {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3} };
362 const static G4String side[4] = { "~-Y", "~+ << 362 G4String side[4] = { "~-Y", "~+Y", "~-X", "~+X" };
363 363
364 for (G4int i=0; i<4; ++i) 364 for (G4int i=0; i<4; ++i)
365 { 365 {
366 if (MakePlane(pt[iface[i][0]], 366 if (MakePlane(pt[iface[i][0]],
367 pt[iface[i][1]], 367 pt[iface[i][1]],
368 pt[iface[i][2]], 368 pt[iface[i][2]],
369 pt[iface[i][3]], 369 pt[iface[i][3]],
370 fPlanes[i])) continue; 370 fPlanes[i])) continue;
371 371
372 // Non planar side face 372 // Non planar side face
373 G4ThreeVector normal(fPlanes[i].a,fPlanes[ 373 G4ThreeVector normal(fPlanes[i].a,fPlanes[i].b,fPlanes[i].c);
374 G4double dmax = 0; 374 G4double dmax = 0;
375 for (G4int k=0; k<4; ++k) 375 for (G4int k=0; k<4; ++k)
376 { 376 {
377 G4double dist = normal.dot(pt[iface[i][k 377 G4double dist = normal.dot(pt[iface[i][k]]) + fPlanes[i].d;
378 if (std::abs(dist) > std::abs(dmax)) dma 378 if (std::abs(dist) > std::abs(dmax)) dmax = dist;
379 } 379 }
380 std::ostringstream message; 380 std::ostringstream message;
381 message << "Side face " << side[i] << " is 381 message << "Side face " << side[i] << " is not planar for solid: "
382 << GetName() << "\nDiscrepancy: " 382 << GetName() << "\nDiscrepancy: " << dmax/mm << " mm\n";
383 StreamInfo(message); 383 StreamInfo(message);
384 G4Exception("G4Trap::MakePlanes()", "GeomS 384 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
385 FatalException, message); 385 FatalException, message);
386 } 386 }
387 387
388 // Re-compute parameters << 388 // Define type of trapezoid
389 SetCachedValues(); << 389 fTrapType = 0;
>> 390 if (fPlanes[0].b == -1 && fPlanes[1].b == 1 &&
>> 391 std::abs(fPlanes[0].a) < DBL_EPSILON &&
>> 392 std::abs(fPlanes[0].c) < DBL_EPSILON &&
>> 393 std::abs(fPlanes[1].a) < DBL_EPSILON &&
>> 394 std::abs(fPlanes[1].c) < DBL_EPSILON)
>> 395 {
>> 396 fTrapType = 1; // YZ section is a rectangle ...
>> 397 if (std::abs(fPlanes[2].a + fPlanes[3].a) < DBL_EPSILON &&
>> 398 std::abs(fPlanes[2].c - fPlanes[3].c) < DBL_EPSILON &&
>> 399 fPlanes[2].b == 0 &&
>> 400 fPlanes[3].b == 0)
>> 401 {
>> 402 fTrapType = 2; // ... and XZ section is a isosceles trapezoid
>> 403 fPlanes[2].a = -fPlanes[3].a;
>> 404 fPlanes[2].c = fPlanes[3].c;
>> 405 }
>> 406 if (std::abs(fPlanes[2].a + fPlanes[3].a) < DBL_EPSILON &&
>> 407 std::abs(fPlanes[2].b - fPlanes[3].b) < DBL_EPSILON &&
>> 408 fPlanes[2].c == 0 &&
>> 409 fPlanes[3].c == 0)
>> 410 {
>> 411 fTrapType = 3; // ... and XY section is a isosceles trapezoid
>> 412 fPlanes[2].a = -fPlanes[3].a;
>> 413 fPlanes[2].b = fPlanes[3].b;
>> 414 }
>> 415 }
390 } 416 }
391 417
392 ////////////////////////////////////////////// << 418 ///////////////////////////////////////////////////////////////////////
393 // 419 //
394 // Calculate the coef's of the plane p1->p2->p 420 // Calculate the coef's of the plane p1->p2->p3->p4->p1
395 // where the ThreeVectors 1-4 are in anti-cloc 421 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed
396 // from infront of the plane (i.e. from normal 422 // from infront of the plane (i.e. from normal direction).
397 // 423 //
398 // Return true if the points are coplanar, fal 424 // Return true if the points are coplanar, false otherwise
399 425
400 G4bool G4Trap::MakePlane( const G4ThreeVector& 426 G4bool G4Trap::MakePlane( const G4ThreeVector& p1,
401 const G4ThreeVector& 427 const G4ThreeVector& p2,
402 const G4ThreeVector& 428 const G4ThreeVector& p3,
403 const G4ThreeVector& 429 const G4ThreeVector& p4,
404 TrapSidePlane& 430 TrapSidePlane& plane )
405 { 431 {
406 G4ThreeVector normal = ((p4 - p2).cross(p3 - 432 G4ThreeVector normal = ((p4 - p2).cross(p3 - p1)).unit();
407 if (std::abs(normal.x()) < DBL_EPSILON) norm << 433 if (std::abs(normal.x()) < DBL_EPSILON) normal.setX(0);
408 if (std::abs(normal.y()) < DBL_EPSILON) norm << 434 if (std::abs(normal.y()) < DBL_EPSILON) normal.setY(0);
409 if (std::abs(normal.z()) < DBL_EPSILON) norm << 435 if (std::abs(normal.z()) < DBL_EPSILON) normal.setZ(0);
410 normal = normal.unit(); 436 normal = normal.unit();
411 437
412 G4ThreeVector centre = (p1 + p2 + p3 + p4)*0 438 G4ThreeVector centre = (p1 + p2 + p3 + p4)*0.25;
413 plane.a = normal.x(); 439 plane.a = normal.x();
414 plane.b = normal.y(); 440 plane.b = normal.y();
415 plane.c = normal.z(); 441 plane.c = normal.z();
416 plane.d = -normal.dot(centre); 442 plane.d = -normal.dot(centre);
417 443
418 // compute distances and check planarity 444 // compute distances and check planarity
419 G4double d1 = std::abs(normal.dot(p1) + plan 445 G4double d1 = std::abs(normal.dot(p1) + plane.d);
420 G4double d2 = std::abs(normal.dot(p2) + plan 446 G4double d2 = std::abs(normal.dot(p2) + plane.d);
421 G4double d3 = std::abs(normal.dot(p3) + plan 447 G4double d3 = std::abs(normal.dot(p3) + plane.d);
422 G4double d4 = std::abs(normal.dot(p4) + plan 448 G4double d4 = std::abs(normal.dot(p4) + plane.d);
423 G4double dmax = std::max(std::max(std::max(d 449 G4double dmax = std::max(std::max(std::max(d1,d2),d3),d4);
424 << 450
425 return dmax <= 1000 * kCarTolerance; << 451 return (dmax > 1000 * kCarTolerance) ? false : true;
426 } <<
427 <<
428 ////////////////////////////////////////////// <<
429 // <<
430 // Recompute parameters using planes <<
431 <<
432 void G4Trap::SetCachedValues() <<
433 { <<
434 // Set indeces <<
435 constexpr G4int iface[6][4] = <<
436 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2, <<
437 <<
438 // Get vertices <<
439 G4ThreeVector pt[8]; <<
440 GetVertices(pt); <<
441 <<
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 } 452 }
481 453
482 ////////////////////////////////////////////// << 454 ///////////////////////////////////////////////////////////////////////
483 // 455 //
484 // Get volume 456 // Get volume
485 457
486 G4double G4Trap::GetCubicVolume() 458 G4double G4Trap::GetCubicVolume()
487 { 459 {
488 if (fCubicVolume == 0) 460 if (fCubicVolume == 0)
489 { 461 {
490 G4ThreeVector pt[8]; 462 G4ThreeVector pt[8];
491 GetVertices(pt); 463 GetVertices(pt);
492 << 464
493 G4double dz = pt[4].z() - pt[0].z(); 465 G4double dz = pt[4].z() - pt[0].z();
494 G4double dy1 = pt[2].y() - pt[0].y(); 466 G4double dy1 = pt[2].y() - pt[0].y();
495 G4double dx1 = pt[1].x() - pt[0].x(); 467 G4double dx1 = pt[1].x() - pt[0].x();
496 G4double dx2 = pt[3].x() - pt[2].x(); 468 G4double dx2 = pt[3].x() - pt[2].x();
497 G4double dy2 = pt[6].y() - pt[4].y(); 469 G4double dy2 = pt[6].y() - pt[4].y();
498 G4double dx3 = pt[5].x() - pt[4].x(); 470 G4double dx3 = pt[5].x() - pt[4].x();
499 G4double dx4 = pt[7].x() - pt[6].x(); 471 G4double dx4 = pt[7].x() - pt[6].x();
500 472
501 fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(d 473 fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(dy1 + dy2) +
502 (dx4 + dx3 - dx2 - dx1)*(d 474 (dx4 + dx3 - dx2 - dx1)*(dy2 - dy1)/3)*dz*0.125;
503 } 475 }
504 return fCubicVolume; 476 return fCubicVolume;
505 } 477 }
506 478
507 ////////////////////////////////////////////// << 479 ///////////////////////////////////////////////////////////////////////
508 // 480 //
509 // Get surface area 481 // Get surface area
510 482
511 G4double G4Trap::GetSurfaceArea() 483 G4double G4Trap::GetSurfaceArea()
512 { 484 {
513 if (fSurfaceArea == 0) 485 if (fSurfaceArea == 0)
514 { 486 {
515 G4ThreeVector pt[8]; 487 G4ThreeVector pt[8];
516 G4int iface [6][4] = 488 G4int iface [6][4] =
517 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2, 489 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
518 490
519 GetVertices(pt); 491 GetVertices(pt);
520 for (const auto & i : iface) << 492 for (G4int i=0; i<6; ++i)
521 { 493 {
522 fSurfaceArea += G4GeomTools::QuadAreaNor << 494 fSurfaceArea += G4GeomTools::QuadAreaNormal(pt[iface[i][0]],
523 << 495 pt[iface[i][1]],
524 << 496 pt[iface[i][2]],
525 << 497 pt[iface[i][3]]).mag();
526 } 498 }
527 } 499 }
528 return fSurfaceArea; 500 return fSurfaceArea;
529 } 501 }
530 502
531 ////////////////////////////////////////////// << 503 ///////////////////////////////////////////////////////////////////////
532 // 504 //
533 // Dispatch to parameterisation for replicatio 505 // Dispatch to parameterisation for replication mechanism dimension
534 // computation & modification. 506 // computation & modification.
535 507
536 void G4Trap::ComputeDimensions( G4VPVPar 508 void G4Trap::ComputeDimensions( G4VPVParameterisation* p,
537 const G4int n, 509 const G4int n,
538 const G4VPhysi 510 const G4VPhysicalVolume* pRep )
539 { 511 {
540 p->ComputeDimensions(*this,n,pRep); 512 p->ComputeDimensions(*this,n,pRep);
541 } 513 }
542 514
543 ////////////////////////////////////////////// << 515 ///////////////////////////////////////////////////////////////////////
544 // 516 //
545 // Get bounding box 517 // Get bounding box
546 518
547 void G4Trap::BoundingLimits(G4ThreeVector& pMi 519 void G4Trap::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
548 { 520 {
549 G4ThreeVector pt[8]; 521 G4ThreeVector pt[8];
550 GetVertices(pt); 522 GetVertices(pt);
551 523
552 G4double xmin = kInfinity, xmax = -kInfinity 524 G4double xmin = kInfinity, xmax = -kInfinity;
553 G4double ymin = kInfinity, ymax = -kInfinity 525 G4double ymin = kInfinity, ymax = -kInfinity;
554 for (const auto & i : pt) << 526 for (G4int i=0; i<8; ++i)
555 { 527 {
556 G4double x = i.x(); << 528 G4double x = pt[i].x();
557 if (x < xmin) xmin = x; 529 if (x < xmin) xmin = x;
558 if (x > xmax) xmax = x; 530 if (x > xmax) xmax = x;
559 G4double y = i.y(); << 531 G4double y = pt[i].y();
560 if (y < ymin) ymin = y; 532 if (y < ymin) ymin = y;
561 if (y > ymax) ymax = y; 533 if (y > ymax) ymax = y;
562 } 534 }
563 535
564 G4double dz = GetZHalfLength(); 536 G4double dz = GetZHalfLength();
565 pMin.set(xmin,ymin,-dz); 537 pMin.set(xmin,ymin,-dz);
566 pMax.set(xmax,ymax, dz); 538 pMax.set(xmax,ymax, dz);
567 539
568 // Check correctness of the bounding box 540 // Check correctness of the bounding box
569 // 541 //
570 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 542 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
571 { 543 {
572 std::ostringstream message; 544 std::ostringstream message;
573 message << "Bad bounding box (min >= max) 545 message << "Bad bounding box (min >= max) for solid: "
574 << GetName() << " !" 546 << GetName() << " !"
575 << "\npMin = " << pMin 547 << "\npMin = " << pMin
576 << "\npMax = " << pMax; 548 << "\npMax = " << pMax;
577 G4Exception("G4Trap::BoundingLimits()", "G 549 G4Exception("G4Trap::BoundingLimits()", "GeomMgt0001",
578 JustWarning, message); 550 JustWarning, message);
579 DumpInfo(); 551 DumpInfo();
580 } 552 }
581 } 553 }
582 554
583 ////////////////////////////////////////////// << 555 ///////////////////////////////////////////////////////////////////////
584 // 556 //
585 // Calculate extent under transform and specif 557 // Calculate extent under transform and specified limit
586 558
587 G4bool G4Trap::CalculateExtent( const EAxis pA 559 G4bool G4Trap::CalculateExtent( const EAxis pAxis,
588 const G4VoxelL 560 const G4VoxelLimits& pVoxelLimit,
589 const G4Affine 561 const G4AffineTransform& pTransform,
590 G4double 562 G4double& pMin, G4double& pMax) const
591 { 563 {
592 G4ThreeVector bmin, bmax; 564 G4ThreeVector bmin, bmax;
593 G4bool exist; 565 G4bool exist;
594 566
595 // Check bounding box (bbox) 567 // Check bounding box (bbox)
596 // 568 //
597 BoundingLimits(bmin,bmax); 569 BoundingLimits(bmin,bmax);
598 G4BoundingEnvelope bbox(bmin,bmax); 570 G4BoundingEnvelope bbox(bmin,bmax);
599 #ifdef G4BBOX_EXTENT 571 #ifdef G4BBOX_EXTENT
600 return bbox.CalculateExtent(pAxis,pVoxelLimi 572 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
601 #endif 573 #endif
602 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 574 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
603 { 575 {
604 return exist = pMin < pMax; << 576 return exist = (pMin < pMax) ? true : false;
605 } 577 }
606 578
607 // Set bounding envelope (benv) and calculat 579 // Set bounding envelope (benv) and calculate extent
608 // 580 //
609 G4ThreeVector pt[8]; 581 G4ThreeVector pt[8];
610 GetVertices(pt); 582 GetVertices(pt);
611 583
612 G4ThreeVectorList baseA(4), baseB(4); 584 G4ThreeVectorList baseA(4), baseB(4);
613 baseA[0] = pt[0]; 585 baseA[0] = pt[0];
614 baseA[1] = pt[1]; 586 baseA[1] = pt[1];
615 baseA[2] = pt[3]; 587 baseA[2] = pt[3];
616 baseA[3] = pt[2]; 588 baseA[3] = pt[2];
617 589
618 baseB[0] = pt[4]; 590 baseB[0] = pt[4];
619 baseB[1] = pt[5]; 591 baseB[1] = pt[5];
620 baseB[2] = pt[7]; 592 baseB[2] = pt[7];
621 baseB[3] = pt[6]; 593 baseB[3] = pt[6];
622 594
623 std::vector<const G4ThreeVectorList *> polyg 595 std::vector<const G4ThreeVectorList *> polygons(2);
624 polygons[0] = &baseA; 596 polygons[0] = &baseA;
625 polygons[1] = &baseB; 597 polygons[1] = &baseB;
626 598
627 G4BoundingEnvelope benv(bmin,bmax,polygons); 599 G4BoundingEnvelope benv(bmin,bmax,polygons);
628 exist = benv.CalculateExtent(pAxis,pVoxelLim 600 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
629 return exist; 601 return exist;
630 } 602 }
631 603
632 ////////////////////////////////////////////// << 604 ///////////////////////////////////////////////////////////////////////
633 // 605 //
634 // Return whether point is inside/outside/on_s 606 // Return whether point is inside/outside/on_surface
635 607
636 EInside G4Trap::Inside( const G4ThreeVector& p 608 EInside G4Trap::Inside( const G4ThreeVector& p ) const
637 { 609 {
>> 610 G4double dz = std::abs(p.z())-fDz;
>> 611 if (dz > halfCarTolerance) return kOutside;
>> 612
638 switch (fTrapType) 613 switch (fTrapType)
639 { 614 {
640 case 0: // General case 615 case 0: // General case
641 { 616 {
642 G4double dz = std::abs(p.z())-fDz; <<
643 G4double dy1 = fPlanes[0].b*p.y()+fPlane 617 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
644 G4double dy2 = fPlanes[1].b*p.y()+fPlane 618 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
645 G4double dy = std::max(dz,std::max(dy1,d 619 G4double dy = std::max(dz,std::max(dy1,dy2));
646 620
647 G4double dx1 = fPlanes[2].a*p.x()+fPlane 621 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
648 + fPlanes[2].c*p.z()+fPlane 622 + fPlanes[2].c*p.z()+fPlanes[2].d;
649 G4double dx2 = fPlanes[3].a*p.x()+fPlane 623 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
650 + fPlanes[3].c*p.z()+fPlane 624 + fPlanes[3].c*p.z()+fPlanes[3].d;
651 G4double dist = std::max(dy,std::max(dx1 625 G4double dist = std::max(dy,std::max(dx1,dx2));
652 626
653 return (dist > halfCarTolerance) ? kOuts << 627 if (dist > halfCarTolerance) return kOutside;
654 ((dist > -halfCarTolerance) ? kSurface << 628 return (dist > -halfCarTolerance) ? kSurface : kInside;
655 } 629 }
656 case 1: // YZ section is a rectangle 630 case 1: // YZ section is a rectangle
657 { 631 {
658 G4double dz = std::abs(p.z())-fDz; <<
659 G4double dy = std::max(dz,std::abs(p.y() 632 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
660 G4double dx1 = fPlanes[2].a*p.x()+fPlane 633 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
661 + fPlanes[2].c*p.z()+fPlane 634 + fPlanes[2].c*p.z()+fPlanes[2].d;
662 G4double dx2 = fPlanes[3].a*p.x()+fPlane 635 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
663 + fPlanes[3].c*p.z()+fPlane 636 + fPlanes[3].c*p.z()+fPlanes[3].d;
664 G4double dist = std::max(dy,std::max(dx1 637 G4double dist = std::max(dy,std::max(dx1,dx2));
665 638
666 return (dist > halfCarTolerance) ? kOuts << 639 if (dist > halfCarTolerance) return kOutside;
667 ((dist > -halfCarTolerance) ? kSurface << 640 return (dist > -halfCarTolerance) ? kSurface : kInside;
668 } 641 }
669 case 2: // YZ section is a rectangle and 642 case 2: // YZ section is a rectangle and
670 { // XZ section is an isosceles trap 643 { // XZ section is an isosceles trapezoid
671 G4double dz = std::abs(p.z())-fDz; <<
672 G4double dy = std::max(dz,std::abs(p.y() 644 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
673 G4double dx = fPlanes[3].a*std::abs(p.x( 645 G4double dx = fPlanes[3].a*std::abs(p.x())
674 + fPlanes[3].c*p.z()+fPlanes 646 + fPlanes[3].c*p.z()+fPlanes[3].d;
675 G4double dist = std::max(dy,dx); 647 G4double dist = std::max(dy,dx);
676 648
677 return (dist > halfCarTolerance) ? kOuts << 649 if (dist > halfCarTolerance) return kOutside;
678 ((dist > -halfCarTolerance) ? kSurface << 650 return (dist > -halfCarTolerance) ? kSurface : kInside;
679 } 651 }
680 case 3: // YZ section is a rectangle and 652 case 3: // YZ section is a rectangle and
681 { // XY section is an isosceles trap 653 { // XY section is an isosceles trapezoid
682 G4double dz = std::abs(p.z())-fDz; <<
683 G4double dy = std::max(dz,std::abs(p.y() 654 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
684 G4double dx = fPlanes[3].a*std::abs(p.x( 655 G4double dx = fPlanes[3].a*std::abs(p.x())
685 + fPlanes[3].b*p.y()+fPlanes 656 + fPlanes[3].b*p.y()+fPlanes[3].d;
686 G4double dist = std::max(dy,dx); 657 G4double dist = std::max(dy,dx);
687 658
688 return (dist > halfCarTolerance) ? kOuts << 659 if (dist > halfCarTolerance) return kOutside;
689 ((dist > -halfCarTolerance) ? kSurface << 660 return (dist > -halfCarTolerance) ? kSurface : kInside;
690 } 661 }
691 } 662 }
692 return kOutside; << 663 return kOutside;
693 } 664 }
694 665
695 ////////////////////////////////////////////// << 666 ///////////////////////////////////////////////////////////////////////
696 // 667 //
697 // Determine side, and return corresponding no 668 // Determine side, and return corresponding normal
698 669
699 G4ThreeVector G4Trap::SurfaceNormal( const G4T 670 G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const
700 { 671 {
>> 672 G4int nsurf = 0; // number of surfaces where p is placed
701 G4double nx = 0, ny = 0, nz = 0; 673 G4double nx = 0, ny = 0, nz = 0;
702 G4double dz = std::abs(p.z()) - fDz; 674 G4double dz = std::abs(p.z()) - fDz;
703 nz = std::copysign(G4double(std::abs(dz) <= << 675 if (std::abs(dz) <= halfCarTolerance)
>> 676 {
>> 677 nz = (p.z() < 0) ? -1 : 1;
>> 678 ++nsurf;
>> 679 }
704 680
705 switch (fTrapType) 681 switch (fTrapType)
706 { 682 {
707 case 0: // General case 683 case 0: // General case
708 { 684 {
709 for (G4int i=0; i<2; ++i) 685 for (G4int i=0; i<2; ++i)
710 { 686 {
711 G4double dy = fPlanes[i].b*p.y() + fPl 687 G4double dy = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
712 if (std::abs(dy) > halfCarTolerance) c 688 if (std::abs(dy) > halfCarTolerance) continue;
713 ny = fPlanes[i].b; 689 ny = fPlanes[i].b;
714 nz += fPlanes[i].c; 690 nz += fPlanes[i].c;
>> 691 ++nsurf;
715 break; 692 break;
716 } 693 }
717 for (G4int i=2; i<4; ++i) 694 for (G4int i=2; i<4; ++i)
718 { 695 {
719 G4double dx = fPlanes[i].a*p.x() + 696 G4double dx = fPlanes[i].a*p.x() +
720 fPlanes[i].b*p.y() + fPl 697 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
721 if (std::abs(dx) > halfCarTolerance) c 698 if (std::abs(dx) > halfCarTolerance) continue;
722 nx = fPlanes[i].a; 699 nx = fPlanes[i].a;
723 ny += fPlanes[i].b; 700 ny += fPlanes[i].b;
724 nz += fPlanes[i].c; 701 nz += fPlanes[i].c;
>> 702 ++nsurf;
725 break; 703 break;
726 } 704 }
727 break; 705 break;
728 } 706 }
729 case 1: // YZ section - rectangle << 707 case 1: // YZ section is a rectangle
730 { 708 {
731 G4double dy = std::abs(p.y()) + fPlanes[ 709 G4double dy = std::abs(p.y()) + fPlanes[1].d;
732 ny = std::copysign(G4double(std::abs(dy) << 710 if (std::abs(dy) <= halfCarTolerance) ny = (p.y() < 0) ? -1 : 1;
733 for (G4int i=2; i<4; ++i) 711 for (G4int i=2; i<4; ++i)
734 { 712 {
735 G4double dx = fPlanes[i].a*p.x() + 713 G4double dx = fPlanes[i].a*p.x() +
736 fPlanes[i].b*p.y() + fPl 714 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
737 if (std::abs(dx) > halfCarTolerance) c 715 if (std::abs(dx) > halfCarTolerance) continue;
738 nx = fPlanes[i].a; 716 nx = fPlanes[i].a;
739 ny += fPlanes[i].b; 717 ny += fPlanes[i].b;
740 nz += fPlanes[i].c; 718 nz += fPlanes[i].c;
>> 719 ++nsurf;
741 break; 720 break;
742 } 721 }
743 break; 722 break;
744 } 723 }
745 case 2: // YZ section - rectangle, XZ sect << 724 case 2: // YZ section is a rectangle and
746 { << 725 { // XZ section is an isosceles trapezoid
747 G4double dy = std::abs(p.y()) + fPlanes[ 726 G4double dy = std::abs(p.y()) + fPlanes[1].d;
748 ny = std::copysign(G4double(std::abs(dy) << 727 if (std::abs(dy) <= halfCarTolerance) ny = (p.y() < 0) ? -1 : 1;
749 G4double dx = fPlanes[3].a*std::abs(p.x( 728 G4double dx = fPlanes[3].a*std::abs(p.x()) +
750 fPlanes[3].c*p.z() + fPlan 729 fPlanes[3].c*p.z() + fPlanes[3].d;
751 G4double k = std::abs(dx) <= halfCarTole << 730 if (std::abs(dx) <= halfCarTolerance)
752 nx = std::copysign(k, p.x())*fPlanes[3] << 731 {
753 nz += k*fPlanes[3].c; << 732 nx = (p.x() < 0) ? -fPlanes[3].a : fPlanes[3].a;
>> 733 nz += fPlanes[3].c;
>> 734 ++nsurf;
>> 735 }
754 break; 736 break;
755 } 737 }
756 case 3: // YZ section - rectangle, XY sect << 738 case 3: // YZ section is a rectangle and
757 { << 739 { // XY section is an isosceles trapezoid
758 G4double dy = std::abs(p.y()) + fPlanes[ 740 G4double dy = std::abs(p.y()) + fPlanes[1].d;
759 ny = std::copysign(G4double(std::abs(dy) << 741 if (std::abs(dy) <= halfCarTolerance) ny = (p.y() < 0) ? -1 : 1;
760 G4double dx = fPlanes[3].a*std::abs(p.x( 742 G4double dx = fPlanes[3].a*std::abs(p.x()) +
761 fPlanes[3].b*p.y() + fPlan 743 fPlanes[3].b*p.y() + fPlanes[3].d;
762 G4double k = std::abs(dx) <= halfCarTole << 744 if (std::abs(dx) <= halfCarTolerance)
763 nx = std::copysign(k, p.x())*fPlanes[3] << 745 {
764 ny += k*fPlanes[3].b; << 746 nx = (p.x() < 0) ? -fPlanes[3].a : fPlanes[3].a;
>> 747 ny += fPlanes[3].b;
>> 748 ++nsurf;
>> 749 }
765 break; 750 break;
766 } 751 }
767 } 752 }
768 753
769 // Return normal 754 // Return normal
770 // 755 //
771 G4double mag2 = nx*nx + ny*ny + nz*nz; << 756 if (nsurf == 1) return G4ThreeVector(nx,ny,nz);
772 if (mag2 == 1) return { nx,ny,nz }; << 757 else if (nsurf != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner
773 else if (mag2 != 0) return G4ThreeVector(nx, <<
774 else 758 else
775 { 759 {
776 // Point is not on the surface 760 // Point is not on the surface
777 // 761 //
778 #ifdef G4CSGDEBUG 762 #ifdef G4CSGDEBUG
779 std::ostringstream message; 763 std::ostringstream message;
780 G4long oldprc = message.precision(16); << 764 G4int oldprc = message.precision(16);
781 message << "Point p is not on surface (!?) 765 message << "Point p is not on surface (!?) of solid: "
782 << GetName() << G4endl; 766 << GetName() << G4endl;
783 message << "Position:\n"; 767 message << "Position:\n";
784 message << " p.x() = " << p.x()/mm << " 768 message << " p.x() = " << p.x()/mm << " mm\n";
785 message << " p.y() = " << p.y()/mm << " 769 message << " p.y() = " << p.y()/mm << " mm\n";
786 message << " p.z() = " << p.z()/mm << " 770 message << " p.z() = " << p.z()/mm << " mm";
787 G4cout.precision(oldprc) ; 771 G4cout.precision(oldprc) ;
788 G4Exception("G4Trap::SurfaceNormal(p)", "G 772 G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
789 JustWarning, message ); 773 JustWarning, message );
790 DumpInfo(); 774 DumpInfo();
791 #endif 775 #endif
792 return ApproxSurfaceNormal(p); 776 return ApproxSurfaceNormal(p);
793 } 777 }
794 } 778 }
795 779
796 ////////////////////////////////////////////// << 780 ///////////////////////////////////////////////////////////////////////
797 // 781 //
798 // Algorithm for SurfaceNormal() following the 782 // Algorithm for SurfaceNormal() following the original specification
799 // for points not on the surface 783 // for points not on the surface
800 784
801 G4ThreeVector G4Trap::ApproxSurfaceNormal( con 785 G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
802 { 786 {
803 G4double dist = -DBL_MAX; 787 G4double dist = -DBL_MAX;
804 G4int iside = 0; 788 G4int iside = 0;
805 for (G4int i=0; i<4; ++i) 789 for (G4int i=0; i<4; ++i)
806 { 790 {
807 G4double d = fPlanes[i].a*p.x() + 791 G4double d = fPlanes[i].a*p.x() +
808 fPlanes[i].b*p.y() + 792 fPlanes[i].b*p.y() +
809 fPlanes[i].c*p.z() + fPlanes[ 793 fPlanes[i].c*p.z() + fPlanes[i].d;
810 if (d > dist) { dist = d; iside = i; } 794 if (d > dist) { dist = d; iside = i; }
811 } 795 }
812 796
813 G4double distz = std::abs(p.z()) - fDz; 797 G4double distz = std::abs(p.z()) - fDz;
814 if (dist > distz) 798 if (dist > distz)
815 return { fPlanes[iside].a, fPlanes[iside]. << 799 return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
816 else 800 else
817 return { 0, 0, (G4double)((p.z() < 0) ? -1 << 801 return G4ThreeVector(0, 0, (p.z() < 0) ? -1 : 1);
818 } 802 }
819 803
820 ////////////////////////////////////////////// << 804 ///////////////////////////////////////////////////////////////////////
821 // 805 //
822 // Calculate distance to shape from outside 806 // Calculate distance to shape from outside
823 // - return kInfinity if no intersection 807 // - return kInfinity if no intersection
824 808
825 G4double G4Trap::DistanceToIn(const G4ThreeVec 809 G4double G4Trap::DistanceToIn(const G4ThreeVector& p,
826 const G4ThreeVec 810 const G4ThreeVector& v ) const
827 { 811 {
828 // Z intersections 812 // Z intersections
829 // 813 //
830 if ((std::abs(p.z()) - fDz) >= -halfCarToler 814 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0)
831 return kInfinity; 815 return kInfinity;
832 G4double invz = (-v.z() == 0) ? DBL_MAX : -1 816 G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z();
833 G4double dz = (invz < 0) ? fDz : -fDz; << 817 G4double dz = (invz < 0) ? fDz : -fDz;
834 G4double tzmin = (p.z() + dz)*invz; 818 G4double tzmin = (p.z() + dz)*invz;
835 G4double tzmax = (p.z() - dz)*invz; 819 G4double tzmax = (p.z() - dz)*invz;
836 820
837 // Y intersections 821 // Y intersections
838 // 822 //
839 G4double tymin = 0, tymax = DBL_MAX; 823 G4double tymin = 0, tymax = DBL_MAX;
840 G4int i = 0; 824 G4int i = 0;
841 for ( ; i<2; ++i) 825 for ( ; i<2; ++i)
842 { << 826 {
843 G4double cosa = fPlanes[i].b*v.y() + fPlan 827 G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
844 G4double dist = fPlanes[i].b*p.y() + fPlan 828 G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
845 if (dist >= -halfCarTolerance) 829 if (dist >= -halfCarTolerance)
846 { 830 {
847 if (cosa >= 0) return kInfinity; 831 if (cosa >= 0) return kInfinity;
848 G4double tmp = -dist/cosa; 832 G4double tmp = -dist/cosa;
849 if (tymin < tmp) tymin = tmp; 833 if (tymin < tmp) tymin = tmp;
850 } 834 }
851 else if (cosa > 0) 835 else if (cosa > 0)
852 { 836 {
853 G4double tmp = -dist/cosa; 837 G4double tmp = -dist/cosa;
854 if (tymax > tmp) tymax = tmp; 838 if (tymax > tmp) tymax = tmp;
855 } << 839 }
856 } 840 }
857 841
858 // Z intersections 842 // Z intersections
859 // 843 //
860 G4double txmin = 0, txmax = DBL_MAX; 844 G4double txmin = 0, txmax = DBL_MAX;
861 for ( ; i<4; ++i) 845 for ( ; i<4; ++i)
862 { << 846 {
863 G4double cosa = fPlanes[i].a*v.x()+fPlanes 847 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
864 G4double dist = fPlanes[i].a*p.x()+fPlanes 848 G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].c*p.z() +
865 fPlanes[i].d; 849 fPlanes[i].d;
866 if (dist >= -halfCarTolerance) 850 if (dist >= -halfCarTolerance)
867 { 851 {
868 if (cosa >= 0) return kInfinity; 852 if (cosa >= 0) return kInfinity;
869 G4double tmp = -dist/cosa; 853 G4double tmp = -dist/cosa;
870 if (txmin < tmp) txmin = tmp; 854 if (txmin < tmp) txmin = tmp;
871 } 855 }
872 else if (cosa > 0) 856 else if (cosa > 0)
873 { 857 {
874 G4double tmp = -dist/cosa; 858 G4double tmp = -dist/cosa;
875 if (txmax > tmp) txmax = tmp; 859 if (txmax > tmp) txmax = tmp;
876 } << 860 }
877 } 861 }
878 862
879 // Find distance 863 // Find distance
880 // 864 //
881 G4double tmin = std::max(std::max(txmin,tymi 865 G4double tmin = std::max(std::max(txmin,tymin),tzmin);
882 G4double tmax = std::min(std::min(txmax,tyma 866 G4double tmax = std::min(std::min(txmax,tymax),tzmax);
883 << 867
884 if (tmax <= tmin + halfCarTolerance) return 868 if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit
885 return (tmin < halfCarTolerance ) ? 0. : tmi 869 return (tmin < halfCarTolerance ) ? 0. : tmin;
886 } 870 }
887 871
888 ////////////////////////////////////////////// << 872 ////////////////////////////////////////////////////////////////////////////
889 // 873 //
890 // Calculate exact shortest distance to any bo 874 // Calculate exact shortest distance to any boundary from outside
891 // This is the best fast estimation of the sho 875 // This is the best fast estimation of the shortest distance to trap
892 // - return 0 if point is inside 876 // - return 0 if point is inside
893 877
894 G4double G4Trap::DistanceToIn( const G4ThreeVe 878 G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const
895 { 879 {
896 switch (fTrapType) 880 switch (fTrapType)
897 { 881 {
898 case 0: // General case 882 case 0: // General case
899 { 883 {
900 G4double dz = std::abs(p.z())-fDz; 884 G4double dz = std::abs(p.z())-fDz;
901 G4double dy1 = fPlanes[0].b*p.y()+fPlane 885 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
902 G4double dy2 = fPlanes[1].b*p.y()+fPlane 886 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
903 G4double dy = std::max(dz,std::max(dy1,d 887 G4double dy = std::max(dz,std::max(dy1,dy2));
904 888
905 G4double dx1 = fPlanes[2].a*p.x()+fPlane 889 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
906 + fPlanes[2].c*p.z()+fPlane 890 + fPlanes[2].c*p.z()+fPlanes[2].d;
907 G4double dx2 = fPlanes[3].a*p.x()+fPlane 891 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
908 + fPlanes[3].c*p.z()+fPlane 892 + fPlanes[3].c*p.z()+fPlanes[3].d;
909 G4double dist = std::max(dy,std::max(dx1 893 G4double dist = std::max(dy,std::max(dx1,dx2));
910 return (dist > 0) ? dist : 0.; 894 return (dist > 0) ? dist : 0.;
911 } 895 }
912 case 1: // YZ section is a rectangle 896 case 1: // YZ section is a rectangle
913 { 897 {
914 G4double dz = std::abs(p.z())-fDz; 898 G4double dz = std::abs(p.z())-fDz;
915 G4double dy = std::max(dz,std::abs(p.y() 899 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
916 G4double dx1 = fPlanes[2].a*p.x()+fPlane 900 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
917 + fPlanes[2].c*p.z()+fPlane 901 + fPlanes[2].c*p.z()+fPlanes[2].d;
918 G4double dx2 = fPlanes[3].a*p.x()+fPlane 902 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
919 + fPlanes[3].c*p.z()+fPlane 903 + fPlanes[3].c*p.z()+fPlanes[3].d;
920 G4double dist = std::max(dy,std::max(dx1 904 G4double dist = std::max(dy,std::max(dx1,dx2));
921 return (dist > 0) ? dist : 0.; 905 return (dist > 0) ? dist : 0.;
922 } 906 }
923 case 2: // YZ section is a rectangle and 907 case 2: // YZ section is a rectangle and
924 { // XZ section is an isosceles trap 908 { // XZ section is an isosceles trapezoid
925 G4double dz = std::abs(p.z())-fDz; 909 G4double dz = std::abs(p.z())-fDz;
926 G4double dy = std::max(dz,std::abs(p.y() 910 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
927 G4double dx = fPlanes[3].a*std::abs(p.x( 911 G4double dx = fPlanes[3].a*std::abs(p.x())
928 + fPlanes[3].c*p.z()+fPlanes 912 + fPlanes[3].c*p.z()+fPlanes[3].d;
929 G4double dist = std::max(dy,dx); 913 G4double dist = std::max(dy,dx);
930 return (dist > 0) ? dist : 0.; 914 return (dist > 0) ? dist : 0.;
931 } 915 }
932 case 3: // YZ section is a rectangle and 916 case 3: // YZ section is a rectangle and
933 { // XY section is an isosceles trap 917 { // XY section is an isosceles trapezoid
934 G4double dz = std::abs(p.z())-fDz; 918 G4double dz = std::abs(p.z())-fDz;
935 G4double dy = std::max(dz,std::abs(p.y() 919 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
936 G4double dx = fPlanes[3].a*std::abs(p.x( 920 G4double dx = fPlanes[3].a*std::abs(p.x())
937 + fPlanes[3].b*p.y()+fPlanes 921 + fPlanes[3].b*p.y()+fPlanes[3].d;
938 G4double dist = std::max(dy,dx); 922 G4double dist = std::max(dy,dx);
939 return (dist > 0) ? dist : 0.; 923 return (dist > 0) ? dist : 0.;
940 } 924 }
941 } 925 }
942 return 0.; 926 return 0.;
943 } 927 }
944 928
945 ////////////////////////////////////////////// << 929 ////////////////////////////////////////////////////////////////////////////
946 // 930 //
947 // Calculate distance to surface of shape from 931 // Calculate distance to surface of shape from inside and
948 // find normal at exit point, if required 932 // find normal at exit point, if required
949 // - when leaving the surface, return 0 933 // - when leaving the surface, return 0
950 934
951 G4double G4Trap::DistanceToOut(const G4ThreeVe 935 G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
952 const G4bool ca 936 const G4bool calcNorm,
953 G4bool* v 937 G4bool* validNorm, G4ThreeVector* n) const
954 { 938 {
955 // Z intersections 939 // Z intersections
956 // 940 //
957 if ((std::abs(p.z()) - fDz) >= -halfCarToler 941 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0)
958 { 942 {
959 if (calcNorm) 943 if (calcNorm)
960 { 944 {
961 *validNorm = true; 945 *validNorm = true;
962 n->set(0, 0, (p.z() < 0) ? -1 : 1); 946 n->set(0, 0, (p.z() < 0) ? -1 : 1);
963 } 947 }
964 return 0; 948 return 0;
965 } 949 }
966 G4double vz = v.z(); 950 G4double vz = v.z();
967 G4double tmax = (vz == 0) ? DBL_MAX : (std:: 951 G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz;
968 G4int iside = (vz < 0) ? -4 : -2; // little 952 G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
969 953
970 // Y intersections 954 // Y intersections
971 // 955 //
972 G4int i = 0; 956 G4int i = 0;
973 for ( ; i<2; ++i) 957 for ( ; i<2; ++i)
974 { 958 {
975 G4double cosa = fPlanes[i].b*v.y() + fPlan 959 G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z();
976 if (cosa > 0) 960 if (cosa > 0)
977 { 961 {
978 G4double dist = fPlanes[i].b*p.y() + fPl 962 G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
979 if (dist >= -halfCarTolerance) 963 if (dist >= -halfCarTolerance)
980 { 964 {
981 if (calcNorm) 965 if (calcNorm)
982 { 966 {
983 *validNorm = true; 967 *validNorm = true;
984 n->set(0, fPlanes[i].b, fPlanes[i].c 968 n->set(0, fPlanes[i].b, fPlanes[i].c);
985 } 969 }
986 return 0; 970 return 0;
987 } 971 }
988 G4double tmp = -dist/cosa; 972 G4double tmp = -dist/cosa;
989 if (tmax > tmp) { tmax = tmp; iside = i; 973 if (tmax > tmp) { tmax = tmp; iside = i; }
990 } 974 }
991 } 975 }
992 976
993 // X intersections 977 // X intersections
994 // 978 //
995 for ( ; i<4; ++i) 979 for ( ; i<4; ++i)
996 { 980 {
997 G4double cosa = fPlanes[i].a*v.x()+fPlanes 981 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
998 if (cosa > 0) 982 if (cosa > 0)
999 { 983 {
1000 G4double dist = fPlanes[i].a*p.x() + 984 G4double dist = fPlanes[i].a*p.x() +
1001 fPlanes[i].b*p.y() + fP 985 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d;
1002 if (dist >= -halfCarTolerance) 986 if (dist >= -halfCarTolerance)
1003 { 987 {
1004 if (calcNorm) 988 if (calcNorm)
1005 { 989 {
1006 *validNorm = true; 990 *validNorm = true;
1007 n->set(fPlanes[i].a, fPlanes[i].b, 991 n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c);
1008 } 992 }
1009 return 0; 993 return 0;
1010 } 994 }
1011 G4double tmp = -dist/cosa; 995 G4double tmp = -dist/cosa;
1012 if (tmax > tmp) { tmax = tmp; iside = i 996 if (tmax > tmp) { tmax = tmp; iside = i; }
1013 } 997 }
1014 } 998 }
1015 999
1016 // Set normal, if required, and return dist 1000 // Set normal, if required, and return distance
1017 // 1001 //
1018 if (calcNorm) << 1002 if (calcNorm)
1019 { 1003 {
1020 *validNorm = true; 1004 *validNorm = true;
1021 if (iside < 0) 1005 if (iside < 0)
1022 n->set(0, 0, iside + 3); // (-4+3)=-1, 1006 n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1
1023 else 1007 else
1024 n->set(fPlanes[iside].a, fPlanes[iside] 1008 n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
1025 } 1009 }
1026 return tmax; 1010 return tmax;
1027 } 1011 }
1028 1012
1029 ///////////////////////////////////////////// << 1013 ////////////////////////////////////////////////////////////////////////////
1030 // 1014 //
1031 // Calculate exact shortest distance to any b 1015 // Calculate exact shortest distance to any boundary from inside
1032 // - Returns 0 is ThreeVector outside 1016 // - Returns 0 is ThreeVector outside
1033 1017
1034 G4double G4Trap::DistanceToOut( const G4Three 1018 G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const
1035 { 1019 {
1036 #ifdef G4CSGDEBUG 1020 #ifdef G4CSGDEBUG
1037 if( Inside(p) == kOutside ) 1021 if( Inside(p) == kOutside )
1038 { 1022 {
1039 std::ostringstream message; 1023 std::ostringstream message;
1040 G4long oldprc = message.precision(16); << 1024 G4int oldprc = message.precision(16);
1041 message << "Point p is outside (!?) of so 1025 message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
1042 message << "Position:\n"; 1026 message << "Position:\n";
1043 message << " p.x() = " << p.x()/mm << " 1027 message << " p.x() = " << p.x()/mm << " mm\n";
1044 message << " p.y() = " << p.y()/mm << " 1028 message << " p.y() = " << p.y()/mm << " mm\n";
1045 message << " p.z() = " << p.z()/mm << " 1029 message << " p.z() = " << p.z()/mm << " mm";
1046 G4cout.precision(oldprc); << 1030 G4cout.precision(oldprc) ;
1047 G4Exception("G4Trap::DistanceToOut(p)", " 1031 G4Exception("G4Trap::DistanceToOut(p)", "GeomSolids1002",
1048 JustWarning, message ); 1032 JustWarning, message );
1049 DumpInfo(); 1033 DumpInfo();
1050 } 1034 }
1051 #endif 1035 #endif
1052 switch (fTrapType) 1036 switch (fTrapType)
1053 { 1037 {
1054 case 0: // General case 1038 case 0: // General case
1055 { 1039 {
1056 G4double dz = std::abs(p.z())-fDz; 1040 G4double dz = std::abs(p.z())-fDz;
1057 G4double dy1 = fPlanes[0].b*p.y()+fPlan 1041 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
1058 G4double dy2 = fPlanes[1].b*p.y()+fPlan 1042 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
1059 G4double dy = std::max(dz,std::max(dy1, 1043 G4double dy = std::max(dz,std::max(dy1,dy2));
1060 1044
1061 G4double dx1 = fPlanes[2].a*p.x()+fPlan 1045 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
1062 + fPlanes[2].c*p.z()+fPlan 1046 + fPlanes[2].c*p.z()+fPlanes[2].d;
1063 G4double dx2 = fPlanes[3].a*p.x()+fPlan 1047 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
1064 + fPlanes[3].c*p.z()+fPlan 1048 + fPlanes[3].c*p.z()+fPlanes[3].d;
1065 G4double dist = std::max(dy,std::max(dx 1049 G4double dist = std::max(dy,std::max(dx1,dx2));
1066 return (dist < 0) ? -dist : 0.; 1050 return (dist < 0) ? -dist : 0.;
1067 } 1051 }
1068 case 1: // YZ section is a rectangle 1052 case 1: // YZ section is a rectangle
1069 { 1053 {
1070 G4double dz = std::abs(p.z())-fDz; 1054 G4double dz = std::abs(p.z())-fDz;
1071 G4double dy = std::max(dz,std::abs(p.y( 1055 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1072 G4double dx1 = fPlanes[2].a*p.x()+fPlan 1056 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()
1073 + fPlanes[2].c*p.z()+fPlan 1057 + fPlanes[2].c*p.z()+fPlanes[2].d;
1074 G4double dx2 = fPlanes[3].a*p.x()+fPlan 1058 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()
1075 + fPlanes[3].c*p.z()+fPlan 1059 + fPlanes[3].c*p.z()+fPlanes[3].d;
1076 G4double dist = std::max(dy,std::max(dx 1060 G4double dist = std::max(dy,std::max(dx1,dx2));
1077 return (dist < 0) ? -dist : 0.; 1061 return (dist < 0) ? -dist : 0.;
1078 } 1062 }
1079 case 2: // YZ section is a rectangle and 1063 case 2: // YZ section is a rectangle and
1080 { // XZ section is an isosceles tra 1064 { // XZ section is an isosceles trapezoid
1081 G4double dz = std::abs(p.z())-fDz; 1065 G4double dz = std::abs(p.z())-fDz;
1082 G4double dy = std::max(dz,std::abs(p.y( 1066 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1083 G4double dx = fPlanes[3].a*std::abs(p.x 1067 G4double dx = fPlanes[3].a*std::abs(p.x())
1084 + fPlanes[3].c*p.z()+fPlane 1068 + fPlanes[3].c*p.z()+fPlanes[3].d;
1085 G4double dist = std::max(dy,dx); 1069 G4double dist = std::max(dy,dx);
1086 return (dist < 0) ? -dist : 0.; 1070 return (dist < 0) ? -dist : 0.;
1087 } 1071 }
1088 case 3: // YZ section is a rectangle and 1072 case 3: // YZ section is a rectangle and
1089 { // XY section is an isosceles tra 1073 { // XY section is an isosceles trapezoid
1090 G4double dz = std::abs(p.z())-fDz; 1074 G4double dz = std::abs(p.z())-fDz;
1091 G4double dy = std::max(dz,std::abs(p.y( 1075 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d);
1092 G4double dx = fPlanes[3].a*std::abs(p.x 1076 G4double dx = fPlanes[3].a*std::abs(p.x())
1093 + fPlanes[3].b*p.y()+fPlane 1077 + fPlanes[3].b*p.y()+fPlanes[3].d;
1094 G4double dist = std::max(dy,dx); 1078 G4double dist = std::max(dy,dx);
1095 return (dist < 0) ? -dist : 0.; 1079 return (dist < 0) ? -dist : 0.;
1096 } 1080 }
1097 } 1081 }
1098 return 0.; 1082 return 0.;
1099 } 1083 }
1100 1084
1101 ///////////////////////////////////////////// << 1085 ////////////////////////////////////////////////////////////////////////////
1102 // 1086 //
1103 // GetEntityType 1087 // GetEntityType
1104 1088
1105 G4GeometryType G4Trap::GetEntityType() const 1089 G4GeometryType G4Trap::GetEntityType() const
1106 { 1090 {
1107 return {"G4Trap"}; << 1091 return G4String("G4Trap");
1108 } <<
1109 <<
1110 ///////////////////////////////////////////// <<
1111 // <<
1112 // IsFaceted <<
1113 <<
1114 G4bool G4Trap::IsFaceted() const <<
1115 { <<
1116 return true; <<
1117 } 1092 }
1118 1093
1119 ///////////////////////////////////////////// 1094 //////////////////////////////////////////////////////////////////////////
1120 // 1095 //
1121 // Make a clone of the object 1096 // Make a clone of the object
1122 // 1097 //
1123 G4VSolid* G4Trap::Clone() const 1098 G4VSolid* G4Trap::Clone() const
1124 { 1099 {
1125 return new G4Trap(*this); 1100 return new G4Trap(*this);
1126 } 1101 }
1127 1102
1128 ///////////////////////////////////////////// 1103 //////////////////////////////////////////////////////////////////////////
1129 // 1104 //
1130 // Stream object contents to an output stream 1105 // Stream object contents to an output stream
1131 1106
1132 std::ostream& G4Trap::StreamInfo( std::ostrea 1107 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1133 { 1108 {
1134 G4double phi = GetPhi(); << 1109 G4double phi = std::atan2(fTthetaSphi,fTthetaCphi);
1135 G4double theta = GetTheta(); << 1110 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
1136 G4double alpha1 = GetAlpha1(); << 1111 +fTthetaSphi*fTthetaSphi));
1137 G4double alpha2 = GetAlpha2(); << 1112 G4double alpha1 = std::atan(fTalpha1);
>> 1113 G4double alpha2 = std::atan(fTalpha2);
>> 1114 G4String signDegree = "\u00B0";
1138 1115
1139 G4long oldprc = os.precision(16); << 1116 G4int oldprc = os.precision(16);
1140 os << "------------------------------------ 1117 os << "-----------------------------------------------------------\n"
1141 << " *** Dump for solid: " << GetName 1118 << " *** Dump for solid: " << GetName() << " ***\n"
1142 << " ================================ 1119 << " ===================================================\n"
1143 << " Solid type: G4Trap\n" 1120 << " Solid type: G4Trap\n"
1144 << " Parameters:\n" 1121 << " Parameters:\n"
1145 << " half length Z: " << fDz/mm << " 1122 << " half length Z: " << fDz/mm << " mm\n"
1146 << " half length Y, face -Dz: " << fD 1123 << " half length Y, face -Dz: " << fDy1/mm << " mm\n"
1147 << " half length X, face -Dz, side -D 1124 << " half length X, face -Dz, side -Dy1: " << fDx1/mm << " mm\n"
1148 << " half length X, face -Dz, side +D 1125 << " half length X, face -Dz, side +Dy1: " << fDx2/mm << " mm\n"
1149 << " half length Y, face +Dz: " << fD 1126 << " half length Y, face +Dz: " << fDy2/mm << " mm\n"
1150 << " half length X, face +Dz, side -D 1127 << " half length X, face +Dz, side -Dy2: " << fDx3/mm << " mm\n"
1151 << " half length X, face +Dz, side +D 1128 << " half length X, face +Dz, side +Dy2: " << fDx4/mm << " mm\n"
1152 << " theta: " << theta/degree << " de << 1129 << " theta: " << theta/degree << signDegree << "\n"
1153 << " phi: " << phi/degree << " degr << 1130 << " phi: " << phi/degree << signDegree << "\n"
1154 << " alpha, face -Dz: " << alpha1/deg << 1131 << " alpha, face -Dz: " << alpha1/degree << signDegree << "\n"
1155 << " alpha, face +Dz: " << alpha2/deg << 1132 << " alpha, face +Dz: " << alpha2/degree << signDegree << "\n"
1156 << "------------------------------------ 1133 << "-----------------------------------------------------------\n";
1157 os.precision(oldprc); 1134 os.precision(oldprc);
1158 1135
1159 return os; 1136 return os;
1160 } 1137 }
1161 1138
1162 ///////////////////////////////////////////// 1139 //////////////////////////////////////////////////////////////////////////
1163 // 1140 //
1164 // Compute vertices from planes 1141 // Compute vertices from planes
1165 1142
1166 void G4Trap::GetVertices(G4ThreeVector pt[8]) 1143 void G4Trap::GetVertices(G4ThreeVector pt[8]) const
1167 { 1144 {
1168 for (G4int i=0; i<8; ++i) 1145 for (G4int i=0; i<8; ++i)
1169 { 1146 {
1170 G4int iy = (i==0 || i==1 || i==4 || i==5) 1147 G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
1171 G4int ix = (i==0 || i==2 || i==4 || i==6) 1148 G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
1172 G4double z = (i < 4) ? -fDz : fDz; 1149 G4double z = (i < 4) ? -fDz : fDz;
1173 G4double y = -(fPlanes[iy].c*z + fPlanes[ 1150 G4double y = -(fPlanes[iy].c*z + fPlanes[iy].d)/fPlanes[iy].b;
1174 G4double x = -(fPlanes[ix].b*y + fPlanes[ 1151 G4double x = -(fPlanes[ix].b*y + fPlanes[ix].c*z
1175 + fPlanes[ix].d)/fPlanes[i 1152 + fPlanes[ix].d)/fPlanes[ix].a;
1176 pt[i].set(x,y,z); 1153 pt[i].set(x,y,z);
1177 } 1154 }
1178 } 1155 }
1179 1156
1180 ///////////////////////////////////////////// << 1157 /////////////////////////////////////////////////////////////////////////
1181 // 1158 //
1182 // Generate random point on the surface 1159 // Generate random point on the surface
1183 1160
1184 G4ThreeVector G4Trap::GetPointOnSurface() con 1161 G4ThreeVector G4Trap::GetPointOnSurface() const
1185 { 1162 {
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]; 1163 G4ThreeVector pt[8];
>> 1164 G4int iface [6][4] =
>> 1165 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} };
>> 1166 G4double sface[6];
>> 1167
1192 GetVertices(pt); 1168 GetVertices(pt);
>> 1169 G4double stotal = 0;
>> 1170 for (G4int i=0; i<6; ++i)
>> 1171 {
>> 1172 G4double ss = G4GeomTools::QuadAreaNormal(pt[iface[i][0]],
>> 1173 pt[iface[i][1]],
>> 1174 pt[iface[i][2]],
>> 1175 pt[iface[i][3]]).mag();
>> 1176 stotal += ss;
>> 1177 sface[i] = stotal;
>> 1178 }
1193 1179
1194 // Select face 1180 // Select face
1195 // 1181 //
1196 G4double select = fAreas[5]*G4QuickRand(); << 1182 G4double select = stotal*G4UniformRand();
1197 G4int k = 5; 1183 G4int k = 5;
1198 k -= (G4int)(select <= fAreas[4]); << 1184 if (select <= sface[4]) k = 4;
1199 k -= (G4int)(select <= fAreas[3]); << 1185 if (select <= sface[3]) k = 3;
1200 k -= (G4int)(select <= fAreas[2]); << 1186 if (select <= sface[2]) k = 2;
1201 k -= (G4int)(select <= fAreas[1]); << 1187 if (select <= sface[1]) k = 1;
1202 k -= (G4int)(select <= fAreas[0]); << 1188 if (select <= sface[0]) k = 0;
1203 1189
1204 // Select sub-triangle 1190 // Select sub-triangle
1205 // 1191 //
1206 G4int i0 = iface[k][0]; 1192 G4int i0 = iface[k][0];
1207 G4int i1 = iface[k][1]; 1193 G4int i1 = iface[k][1];
1208 G4int i2 = iface[k][2]; 1194 G4int i2 = iface[k][2];
1209 G4int i3 = iface[k][3]; 1195 G4int i3 = iface[k][3];
>> 1196 G4double s1 = G4GeomTools::TriangleAreaNormal(pt[i0],pt[i1],pt[i3]).mag();
1210 G4double s2 = G4GeomTools::TriangleAreaNorm 1197 G4double s2 = G4GeomTools::TriangleAreaNormal(pt[i2],pt[i1],pt[i3]).mag();
1211 if (select > fAreas[k] - s2) i0 = i2; << 1198 if ((s1+s2)*G4UniformRand() > s1) i0 = i2;
1212 1199
1213 // Generate point 1200 // Generate point
1214 // 1201 //
1215 G4double u = G4QuickRand(); << 1202 G4double u = G4UniformRand();
1216 G4double v = G4QuickRand(); << 1203 G4double v = G4UniformRand();
1217 if (u + v > 1.) { u = 1. - u; v = 1. - v; } 1204 if (u + v > 1.) { u = 1. - u; v = 1. - v; }
1218 return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3 1205 return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3];
1219 } 1206 }
1220 1207
1221 ///////////////////////////////////////////// 1208 //////////////////////////////////////////////////////////////////////////
1222 // 1209 //
1223 // Methods for visualisation 1210 // Methods for visualisation
1224 1211
1225 void G4Trap::DescribeYourselfTo ( G4VGraphics 1212 void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1226 { 1213 {
1227 scene.AddSolid (*this); 1214 scene.AddSolid (*this);
1228 } 1215 }
1229 1216
1230 G4Polyhedron* G4Trap::CreatePolyhedron () con 1217 G4Polyhedron* G4Trap::CreatePolyhedron () const
1231 { 1218 {
1232 G4double phi = std::atan2(fTthetaSphi, fTth 1219 G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1233 G4double alpha1 = std::atan(fTalpha1); 1220 G4double alpha1 = std::atan(fTalpha1);
1234 G4double alpha2 = std::atan(fTalpha2); 1221 G4double alpha2 = std::atan(fTalpha2);
1235 G4double theta = std::atan(std::sqrt(fTthet 1222 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
1236 +fTthet 1223 +fTthetaSphi*fTthetaSphi));
1237 1224
1238 return new G4PolyhedronTrap(fDz, theta, phi 1225 return new G4PolyhedronTrap(fDz, theta, phi,
1239 fDy1, fDx1, fDx 1226 fDy1, fDx1, fDx2, alpha1,
1240 fDy2, fDx3, fDx 1227 fDy2, fDx3, fDx4, alpha2);
1241 } 1228 }
1242 1229
1243 #endif 1230 #endif
1244 1231