| 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 G4UTrap wrapper class 26 // Implementation for G4UTrap wrapper class
27 // 27 //
28 // 13.09.13 G.Cosmo, CERN 28 // 13.09.13 G.Cosmo, CERN
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4Trap.hh" 31 #include "G4Trap.hh"
32 #include "G4UTrap.hh" 32 #include "G4UTrap.hh"
33 33
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G 34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
35 35
36 #include "G4AffineTransform.hh" 36 #include "G4AffineTransform.hh"
37 #include "G4VPVParameterisation.hh" 37 #include "G4VPVParameterisation.hh"
38 #include "G4BoundingEnvelope.hh" 38 #include "G4BoundingEnvelope.hh"
39 39
40 using namespace CLHEP; 40 using namespace CLHEP;
41 41
42 ////////////////////////////////////////////// 42 /////////////////////////////////////////////////////////////////////////
43 // 43 //
44 // Constructors 44 // Constructors
45 // 45 //
46 G4UTrap::G4UTrap( const G4String& pName, 46 G4UTrap::G4UTrap( const G4String& pName,
47 G4double pdz, 47 G4double pdz,
48 G4double pTheta, G4dou 48 G4double pTheta, G4double pPhi,
49 G4double pdy1, G4doubl 49 G4double pdy1, G4double pdx1, G4double pdx2,
50 G4double pAlp1, 50 G4double pAlp1,
51 G4double pdy2, G4doubl 51 G4double pdy2, G4double pdx3, G4double pdx4,
52 G4double pAlp2 ) 52 G4double pAlp2 )
53 : Base_t(pName, pdz, pTheta, pPhi, pdy1, pdx 53 : Base_t(pName, pdz, pTheta, pPhi, pdy1, pdx1, pdx2,
54 pAlp1, pdy2, pdx3, pdx4, pAlp2) 54 pAlp1, pdy2, pdx3, pdx4, pAlp2)
55 { 55 {
56 G4ThreeVector pt[8]; 56 G4ThreeVector pt[8];
57 CheckParameters(); 57 CheckParameters();
58 GetVertices(pt); 58 GetVertices(pt);
59 CheckPlanarity(pt); 59 CheckPlanarity(pt);
60 } 60 }
61 61
62 G4UTrap::G4UTrap( const G4String& pName, 62 G4UTrap::G4UTrap( const G4String& pName,
63 const G4ThreeVector pt[8] ) 63 const G4ThreeVector pt[8] )
64 : Base_t(pName) 64 : Base_t(pName)
65 { 65 {
66 // Start with check of centering - the cente 66 // Start with check of centering - the center of gravity trap line
67 // should cross the origin of frame 67 // should cross the origin of frame
68 if ( pt[0].z() >= 0 68 if ( pt[0].z() >= 0
69 || pt[0].z() != pt[1].z() 69 || pt[0].z() != pt[1].z()
70 || pt[0].z() != pt[2].z() 70 || pt[0].z() != pt[2].z()
71 || pt[0].z() != pt[3].z() 71 || pt[0].z() != pt[3].z()
72 72
73 || pt[4].z() <= 0 73 || pt[4].z() <= 0
74 || pt[4].z() != pt[5].z() 74 || pt[4].z() != pt[5].z()
75 || pt[4].z() != pt[6].z() 75 || pt[4].z() != pt[6].z()
76 || pt[4].z() != pt[7].z() 76 || pt[4].z() != pt[7].z()
77 77
78 || std::abs( pt[0].z() + pt[4].z() ) > 78 || std::abs( pt[0].z() + pt[4].z() ) >= kCarTolerance
79 79
80 || pt[0].y() != pt[1].y() 80 || pt[0].y() != pt[1].y()
81 || pt[2].y() != pt[3].y() 81 || pt[2].y() != pt[3].y()
82 || pt[4].y() != pt[5].y() 82 || pt[4].y() != pt[5].y()
83 || pt[6].y() != pt[7].y() 83 || pt[6].y() != pt[7].y()
84 84
85 || std::abs(pt[0].y()+pt[2].y()+pt[4]. 85 || std::abs(pt[0].y()+pt[2].y()+pt[4].y()+pt[6].y()) >= kCarTolerance
86 || std::abs(pt[0].x()+pt[1].x()+pt[4]. 86 || std::abs(pt[0].x()+pt[1].x()+pt[4].x()+pt[5].x() +
87 pt[2].x()+pt[3].x()+pt[6]. 87 pt[2].x()+pt[3].x()+pt[6].x()+pt[7].x()) >= kCarTolerance )
88 { 88 {
89 std::ostringstream message; 89 std::ostringstream message;
90 message << "Invalid vertice coordinates fo 90 message << "Invalid vertice coordinates for Solid: " << GetName();
91 G4Exception("G4UTrap::G4UTrap()", "GeomSol 91 G4Exception("G4UTrap::G4UTrap()", "GeomSolids0002",
92 FatalException, message); 92 FatalException, message);
93 } 93 }
94 94
95 SetPlanes(pt); 95 SetPlanes(pt);
96 CheckParameters(); 96 CheckParameters();
97 CheckPlanarity(pt); 97 CheckPlanarity(pt);
98 } 98 }
99 99
100 // Constructor for Right Angular Wedge from ST 100 // Constructor for Right Angular Wedge from STEP; this is different from
101 // the UnplacedTrapezoid constructor taking fo 101 // the UnplacedTrapezoid constructor taking four double's and constructing
102 // a Trd1. 102 // a Trd1.
103 G4UTrap::G4UTrap( const G4String& pName, 103 G4UTrap::G4UTrap( const G4String& pName,
104 G4double pZ, 104 G4double pZ,
105 G4double pY, 105 G4double pY,
106 G4double pX, G4double 106 G4double pX, G4double pLTX )
107 : Base_t(pName, 0.5*pZ, /*theta=*/0, /*phi=* 107 : Base_t(pName, 0.5*pZ, /*theta=*/0, /*phi=*/0, /*dy1=*/0.5*pY,
108 /*dx1=*/0.5*pX, /*dx2=*/0.5*pLTX, / 108 /*dx1=*/0.5*pX, /*dx2=*/0.5*pLTX, /*Alpha1=*/0.5*(pLTX - pX)/pY,
109 /*dy2=*/0.5*pY, /*dx3=*/0.5*pX, /*d 109 /*dy2=*/0.5*pY, /*dx3=*/0.5*pX, /*dx4=*/0.5*pLTX,
110 /*Alpha2=*/0.5*(pLTX - pX)/pY) 110 /*Alpha2=*/0.5*(pLTX - pX)/pY)
111 { 111 {
112 CheckParameters(); 112 CheckParameters();
113 } 113 }
114 114
115 G4UTrap::G4UTrap( const G4String& pName, 115 G4UTrap::G4UTrap( const G4String& pName,
116 G4double pdx1, G4doub 116 G4double pdx1, G4double pdx2,
117 G4double pdy1, G4doub 117 G4double pdy1, G4double pdy2,
118 G4double pdz ) 118 G4double pdz )
119 : Base_t(pName, pdx1, pdx2, pdy1, pdy2, pdz) 119 : Base_t(pName, pdx1, pdx2, pdy1, pdy2, pdz)
120 { 120 {
121 CheckParameters(); 121 CheckParameters();
122 } 122 }
123 123
124 G4UTrap::G4UTrap(const G4String& pName, 124 G4UTrap::G4UTrap(const G4String& pName,
125 G4double pdx, G4double 125 G4double pdx, G4double pdy, G4double pdz,
126 G4double pAlpha, G4doub 126 G4double pAlpha, G4double pTheta, G4double pPhi )
127 : Base_t(pName, pdx, pdy, pdz, pAlpha, pThet 127 : Base_t(pName, pdx, pdy, pdz, pAlpha, pTheta, pPhi)
128 { 128 {
129 CheckParameters(); 129 CheckParameters();
130 } 130 }
131 131
132 G4UTrap::G4UTrap( const G4String& pName ) 132 G4UTrap::G4UTrap( const G4String& pName )
133 : Base_t(pName) 133 : Base_t(pName)
134 { 134 {
135 } 135 }
136 136
137 ////////////////////////////////////////////// 137 ///////////////////////////////////////////////////////////////////////
138 // 138 //
139 // Fake default constructor - sets only member 139 // Fake default constructor - sets only member data and allocates memory
140 // for usage restri 140 // for usage restricted to object persistency.
141 // 141 //
142 G4UTrap::G4UTrap( __void__& a ) 142 G4UTrap::G4UTrap( __void__& a )
143 : Base_t(a) 143 : Base_t(a)
144 { 144 {
145 } 145 }
146 146
147 ////////////////////////////////////////////// 147 //////////////////////////////////////////////////////////////////////////
148 // 148 //
149 // Destructor 149 // Destructor
150 // 150 //
151 G4UTrap::~G4UTrap() = default; 151 G4UTrap::~G4UTrap() = default;
152 152
153 ////////////////////////////////////////////// 153 //////////////////////////////////////////////////////////////////////////
154 // 154 //
155 // Copy constructor 155 // Copy constructor
156 // 156 //
157 G4UTrap::G4UTrap(const G4UTrap& rhs) 157 G4UTrap::G4UTrap(const G4UTrap& rhs)
158 : Base_t(rhs) 158 : Base_t(rhs)
159 { 159 {
160 } 160 }
161 161
162 ////////////////////////////////////////////// 162 //////////////////////////////////////////////////////////////////////////
163 // 163 //
164 // Assignment operator 164 // Assignment operator
165 // 165 //
166 G4UTrap& G4UTrap::operator = (const G4UTrap& r 166 G4UTrap& G4UTrap::operator = (const G4UTrap& rhs)
167 { 167 {
168 // Check assignment to self 168 // Check assignment to self
169 // 169 //
170 if (this == &rhs) { return *this; } 170 if (this == &rhs) { return *this; }
171 171
172 // Copy base class data 172 // Copy base class data
173 // 173 //
174 Base_t::operator=(rhs); 174 Base_t::operator=(rhs);
175 175
176 return *this; 176 return *this;
177 } 177 }
178 178
179 ////////////////////////////////////////////// 179 //////////////////////////////////////////////////////////////////////////
180 // 180 //
181 // Accessors 181 // Accessors
182 // 182 //
183 G4double G4UTrap::GetZHalfLength() const 183 G4double G4UTrap::GetZHalfLength() const
184 { 184 {
185 return GetDz(); 185 return GetDz();
186 } 186 }
187 G4double G4UTrap::GetYHalfLength1() const 187 G4double G4UTrap::GetYHalfLength1() const
188 { 188 {
189 return GetDy1(); 189 return GetDy1();
190 } 190 }
191 G4double G4UTrap::GetXHalfLength1() const 191 G4double G4UTrap::GetXHalfLength1() const
192 { 192 {
193 return GetDx1(); 193 return GetDx1();
194 } 194 }
195 G4double G4UTrap::GetXHalfLength2() const 195 G4double G4UTrap::GetXHalfLength2() const
196 { 196 {
197 return GetDx2(); 197 return GetDx2();
198 } 198 }
199 G4double G4UTrap::GetTanAlpha1() const 199 G4double G4UTrap::GetTanAlpha1() const
200 { 200 {
201 return Base_t::GetTanAlpha1(); 201 return Base_t::GetTanAlpha1();
202 } 202 }
203 G4double G4UTrap::GetYHalfLength2() const 203 G4double G4UTrap::GetYHalfLength2() const
204 { 204 {
205 return GetDy2(); 205 return GetDy2();
206 } 206 }
207 G4double G4UTrap::GetXHalfLength3() const 207 G4double G4UTrap::GetXHalfLength3() const
208 { 208 {
209 return GetDx3(); 209 return GetDx3();
210 } 210 }
211 G4double G4UTrap::GetXHalfLength4() const 211 G4double G4UTrap::GetXHalfLength4() const
212 { 212 {
213 return GetDx4(); 213 return GetDx4();
214 } 214 }
215 G4double G4UTrap::GetTanAlpha2() const 215 G4double G4UTrap::GetTanAlpha2() const
216 { 216 {
217 return Base_t::GetTanAlpha2(); 217 return Base_t::GetTanAlpha2();
218 } 218 }
219 G4double G4UTrap::GetPhi() const 219 G4double G4UTrap::GetPhi() const
220 { 220 {
221 return Base_t::GetPhi(); 221 return Base_t::GetPhi();
222 } 222 }
223 G4double G4UTrap::GetTheta() const 223 G4double G4UTrap::GetTheta() const
224 { 224 {
225 return Base_t::GetTheta(); 225 return Base_t::GetTheta();
226 } 226 }
227 G4double G4UTrap::GetAlpha1() const 227 G4double G4UTrap::GetAlpha1() const
228 { 228 {
229 return Base_t::GetAlpha1(); 229 return Base_t::GetAlpha1();
230 } 230 }
231 G4double G4UTrap::GetAlpha2() const 231 G4double G4UTrap::GetAlpha2() const
232 { 232 {
233 return Base_t::GetAlpha2(); 233 return Base_t::GetAlpha2();
234 } 234 }
235 TrapSidePlane G4UTrap::GetSidePlane(G4int n) c 235 TrapSidePlane G4UTrap::GetSidePlane(G4int n) const
236 { 236 {
237 TrapSidePlane plane; 237 TrapSidePlane plane;
238 plane.a = GetStruct().GetPlane(n).fA; 238 plane.a = GetStruct().GetPlane(n).fA;
239 plane.b = GetStruct().GetPlane(n).fB; 239 plane.b = GetStruct().GetPlane(n).fB;
240 plane.c = GetStruct().GetPlane(n).fC; 240 plane.c = GetStruct().GetPlane(n).fC;
241 plane.d = GetStruct().GetPlane(n).fD; 241 plane.d = GetStruct().GetPlane(n).fD;
242 return plane; 242 return plane;
243 } 243 }
244 G4ThreeVector G4UTrap::GetSymAxis() const 244 G4ThreeVector G4UTrap::GetSymAxis() const
245 { 245 {
246 G4double tanThetaSphi = GetTanThetaSinPhi(); 246 G4double tanThetaSphi = GetTanThetaSinPhi();
247 G4double tanThetaCphi = GetTanThetaCosPhi(); 247 G4double tanThetaCphi = GetTanThetaCosPhi();
248 G4double tan2Theta = tanThetaSphi*tanThetaSp 248 G4double tan2Theta = tanThetaSphi*tanThetaSphi + tanThetaCphi*tanThetaCphi;
249 G4double cosTheta = 1.0 / std::sqrt(1 + tan2 249 G4double cosTheta = 1.0 / std::sqrt(1 + tan2Theta);
250 return {tanThetaCphi*cosTheta, tanThetaSphi* 250 return {tanThetaCphi*cosTheta, tanThetaSphi*cosTheta, cosTheta};
251 } 251 }
252 252
253 ////////////////////////////////////////////// 253 //////////////////////////////////////////////////////////////////////////
254 // 254 //
255 // Modifier 255 // Modifier
256 // 256 //
257 void G4UTrap::SetAllParameters(G4double pDz, G 257 void G4UTrap::SetAllParameters(G4double pDz, G4double pTheta, G4double pPhi,
258 G4double pDy1, 258 G4double pDy1, G4double pDx1, G4double pDx2,
259 G4double pAlp1, 259 G4double pAlp1,
260 G4double pDy2, 260 G4double pDy2, G4double pDx3, G4double pDx4,
261 G4double pAlp2) 261 G4double pAlp2)
262 { 262 {
263 SetDz(pDz); 263 SetDz(pDz);
264 SetDy1(pDy1); 264 SetDy1(pDy1);
265 SetDy2(pDy2); 265 SetDy2(pDy2);
266 SetDx1(pDx1); 266 SetDx1(pDx1);
267 SetDx2(pDx2); 267 SetDx2(pDx2);
268 SetDx3(pDx3); 268 SetDx3(pDx3);
269 SetDx4(pDx4); 269 SetDx4(pDx4);
270 SetTanAlpha1(std::tan(pAlp1)); 270 SetTanAlpha1(std::tan(pAlp1));
271 SetTanAlpha2(std::tan(pAlp2)); 271 SetTanAlpha2(std::tan(pAlp2));
272 // last two will also reset cached variables 272 // last two will also reset cached variables
273 SetTheta(pTheta); 273 SetTheta(pTheta);
274 SetPhi(pPhi); 274 SetPhi(pPhi);
275 fRebuildPolyhedron = true; 275 fRebuildPolyhedron = true;
276 276
277 G4ThreeVector pt[8]; 277 G4ThreeVector pt[8];
278 CheckParameters(); 278 CheckParameters();
279 GetVertices(pt); 279 GetVertices(pt);
280 CheckPlanarity(pt); 280 CheckPlanarity(pt);
281 } 281 }
282 282
283 ////////////////////////////////////////////// 283 /////////////////////////////////////////////////////////////////////////
284 // 284 //
285 // Set parameters using eight vertices 285 // Set parameters using eight vertices
286 // 286 //
287 void G4UTrap::SetPlanes(const G4ThreeVector pt 287 void G4UTrap::SetPlanes(const G4ThreeVector pt[8])
288 { 288 {
289 U3Vector upt[8]; 289 U3Vector upt[8];
290 for (unsigned int i=0; i<8; ++i) 290 for (unsigned int i=0; i<8; ++i)
291 { 291 {
292 upt[i] = U3Vector(pt[i].x(), pt[i].y(), pt 292 upt[i] = U3Vector(pt[i].x(), pt[i].y(), pt[i].z());
293 } 293 }
294 fromCornersToParameters(upt); 294 fromCornersToParameters(upt);
295 fRebuildPolyhedron = true; 295 fRebuildPolyhedron = true;
296 } 296 }
297 297
298 ////////////////////////////////////////////// 298 /////////////////////////////////////////////////////////////////////////
299 // 299 //
300 // Check dimensions 300 // Check dimensions
301 // 301 //
302 void G4UTrap::CheckParameters() const 302 void G4UTrap::CheckParameters() const
303 { 303 {
304 G4double fDz = GetZHalfLength(); 304 G4double fDz = GetZHalfLength();
305 G4double fDy1 = GetYHalfLength1(); 305 G4double fDy1 = GetYHalfLength1();
306 G4double fDx1 = GetXHalfLength1(); 306 G4double fDx1 = GetXHalfLength1();
307 G4double fDx2 = GetXHalfLength2(); 307 G4double fDx2 = GetXHalfLength2();
308 G4double fDy2 = GetYHalfLength2(); 308 G4double fDy2 = GetYHalfLength2();
309 G4double fDx3 = GetXHalfLength3(); 309 G4double fDx3 = GetXHalfLength3();
310 G4double fDx4 = GetXHalfLength4(); 310 G4double fDx4 = GetXHalfLength4();
311 311
312 if (fDz<=0 || 312 if (fDz<=0 ||
313 fDy1<=0 || fDx1<=0 || fDx2<=0 || 313 fDy1<=0 || fDx1<=0 || fDx2<=0 ||
314 fDy2<=0 || fDx3<=0 || fDx4<=0) 314 fDy2<=0 || fDx3<=0 || fDx4<=0)
315 { 315 {
316 std::ostringstream message; 316 std::ostringstream message;
317 message << "Invalid Length Parameters for 317 message << "Invalid Length Parameters for Solid: " << GetName()
318 << "\n X - " <<fDx1<<", "<<fDx2<< 318 << "\n X - " <<fDx1<<", "<<fDx2<<", "<<fDx3<<", "<<fDx4
319 << "\n Y - " <<fDy1<<", "<<fDy2 319 << "\n Y - " <<fDy1<<", "<<fDy2
320 << "\n Z - " <<fDz; 320 << "\n Z - " <<fDz;
321 G4Exception("G4UTrap::CheckParameters()", 321 G4Exception("G4UTrap::CheckParameters()", "GeomSolids0002",
322 FatalException, message); 322 FatalException, message);
323 } 323 }
324 } 324 }
325 325
326 ////////////////////////////////////////////// 326 /////////////////////////////////////////////////////////////////////////
327 // 327 //
328 // Compute coordinates of vertices 328 // Compute coordinates of vertices
329 // 329 //
330 void G4UTrap::GetVertices(G4ThreeVector pt[8]) 330 void G4UTrap::GetVertices(G4ThreeVector pt[8]) const
331 { 331 {
332 G4double fDz = GetZHalfLength(); 332 G4double fDz = GetZHalfLength();
333 G4double fDy1 = GetYHalfLength1(); 333 G4double fDy1 = GetYHalfLength1();
334 G4double fDx1 = GetXHalfLength1(); 334 G4double fDx1 = GetXHalfLength1();
335 G4double fDx2 = GetXHalfLength2(); 335 G4double fDx2 = GetXHalfLength2();
336 G4double fDy2 = GetYHalfLength2(); 336 G4double fDy2 = GetYHalfLength2();
337 G4double fDx3 = GetXHalfLength3(); 337 G4double fDx3 = GetXHalfLength3();
338 G4double fDx4 = GetXHalfLength4(); 338 G4double fDx4 = GetXHalfLength4();
339 G4double fTalpha1 = GetTanAlpha1(); 339 G4double fTalpha1 = GetTanAlpha1();
340 G4double fTalpha2 = GetTanAlpha2(); 340 G4double fTalpha2 = GetTanAlpha2();
341 341
342 G4double DzTthetaCphi = fDz*GetTanThetaCosPh 342 G4double DzTthetaCphi = fDz*GetTanThetaCosPhi();
343 G4double DzTthetaSphi = fDz*GetTanThetaSinPh 343 G4double DzTthetaSphi = fDz*GetTanThetaSinPhi();
344 G4double Dy1Talpha1 = fDy1*fTalpha1; 344 G4double Dy1Talpha1 = fDy1*fTalpha1;
345 G4double Dy2Talpha2 = fDy2*fTalpha2; 345 G4double Dy2Talpha2 = fDy2*fTalpha2;
346 346
347 pt[0].set(-DzTthetaCphi-Dy1Talpha1-fDx1,-DzT 347 pt[0].set(-DzTthetaCphi-Dy1Talpha1-fDx1,-DzTthetaSphi-fDy1,-fDz);
348 pt[1].set(-DzTthetaCphi-Dy1Talpha1+fDx1,-DzT 348 pt[1].set(-DzTthetaCphi-Dy1Talpha1+fDx1,-DzTthetaSphi-fDy1,-fDz);
349 pt[2].set(-DzTthetaCphi+Dy1Talpha1-fDx2,-DzT 349 pt[2].set(-DzTthetaCphi+Dy1Talpha1-fDx2,-DzTthetaSphi+fDy1,-fDz);
350 pt[3].set(-DzTthetaCphi+Dy1Talpha1+fDx2,-DzT 350 pt[3].set(-DzTthetaCphi+Dy1Talpha1+fDx2,-DzTthetaSphi+fDy1,-fDz);
351 pt[4].set( DzTthetaCphi-Dy2Talpha2-fDx3, DzT 351 pt[4].set( DzTthetaCphi-Dy2Talpha2-fDx3, DzTthetaSphi-fDy2, fDz);
352 pt[5].set( DzTthetaCphi-Dy2Talpha2+fDx3, DzT 352 pt[5].set( DzTthetaCphi-Dy2Talpha2+fDx3, DzTthetaSphi-fDy2, fDz);
353 pt[6].set( DzTthetaCphi+Dy2Talpha2-fDx4, DzT 353 pt[6].set( DzTthetaCphi+Dy2Talpha2-fDx4, DzTthetaSphi+fDy2, fDz);
354 pt[7].set( DzTthetaCphi+Dy2Talpha2+fDx4, DzT 354 pt[7].set( DzTthetaCphi+Dy2Talpha2+fDx4, DzTthetaSphi+fDy2, fDz);
355 } 355 }
356 356
357 ////////////////////////////////////////////// 357 /////////////////////////////////////////////////////////////////////////
358 // 358 //
359 // Check planarity of lateral planes 359 // Check planarity of lateral planes
360 // 360 //
361 void G4UTrap::CheckPlanarity(const G4ThreeVect 361 void G4UTrap::CheckPlanarity(const G4ThreeVector pt[8]) const
362 { 362 {
363 constexpr G4int iface[4][4] = { {0,4,5,1}, { 363 constexpr G4int iface[4][4] = { {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3} };
364 const static G4String side[4] = { "~-Y", "~+ 364 const static G4String side[4] = { "~-Y", "~+Y", "~-X", "~+X" };
365 365
366 for (G4int i=0; i<4; ++i) 366 for (G4int i=0; i<4; ++i)
367 { 367 {
368 TrapSidePlane plane = GetSidePlane(i); 368 TrapSidePlane plane = GetSidePlane(i);
369 G4double dmax = 0; 369 G4double dmax = 0;
370 for (G4int k=0; k<4; ++k) 370 for (G4int k=0; k<4; ++k)
371 { 371 {
372 const G4ThreeVector p = pt[iface[i][k]]; 372 const G4ThreeVector p = pt[iface[i][k]];
373 G4double dist = plane.a*p.x() + plane.b* 373 G4double dist = plane.a*p.x() + plane.b*p.y() + plane.c*p.z() + plane.d;
374 if (std::abs(dist) > std::abs(dmax)) dma 374 if (std::abs(dist) > std::abs(dmax)) dmax = dist;
375 } 375 }
376 if (std::abs(dmax) > 1000 * kCarTolerance) 376 if (std::abs(dmax) > 1000 * kCarTolerance)
377 { 377 {
378 std::ostringstream message; 378 std::ostringstream message;
379 message << "Side face " << side[i] << " 379 message << "Side face " << side[i] << " is not planar for solid: "
380 << GetName() << "\nDiscrepancy: 380 << GetName() << "\nDiscrepancy: " << dmax/mm << " mm\n";
381 StreamInfo(message); 381 StreamInfo(message);
382 G4Exception("G4UTrap::CheckPlanarity()", 382 G4Exception("G4UTrap::CheckPlanarity()", "GeomSolids0002",
383 FatalException, message); 383 FatalException, message);
384 } 384 }
385 } 385 }
386 } 386 }
387 387
388 ////////////////////////////////////////////// 388 /////////////////////////////////////////////////////////////////////////
389 // 389 //
390 // Dispatch to parameterisation for replicatio 390 // Dispatch to parameterisation for replication mechanism dimension
391 // computation & modification. 391 // computation & modification.
392 // 392 //
393 void G4UTrap::ComputeDimensions( G4VPVPar 393 void G4UTrap::ComputeDimensions( G4VPVParameterisation* p,
394 const G4int n, 394 const G4int n,
395 const G4VPhysi 395 const G4VPhysicalVolume* pRep)
396 { 396 {
397 p->ComputeDimensions(*(G4Trap*)this,n,pRep); 397 p->ComputeDimensions(*(G4Trap*)this,n,pRep);
398 } 398 }
399 399
400 ////////////////////////////////////////////// 400 //////////////////////////////////////////////////////////////////////////
401 // 401 //
402 // Make a clone of the object 402 // Make a clone of the object
403 // 403 //
404 G4VSolid* G4UTrap::Clone() const 404 G4VSolid* G4UTrap::Clone() const
405 { 405 {
406 return new G4UTrap(*this); 406 return new G4UTrap(*this);
407 } 407 }
408 408
409 ////////////////////////////////////////////// 409 //////////////////////////////////////////////////////////////////////////
410 // 410 //
411 // Get bounding box 411 // Get bounding box
412 412
413 void G4UTrap::BoundingLimits(G4ThreeVector& pM 413 void G4UTrap::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
414 { 414 {
415 static G4bool checkBBox = true; 415 static G4bool checkBBox = true;
416 416
417 TrapSidePlane planes[4]; 417 TrapSidePlane planes[4];
418 for (G4int i=0; i<4; ++i) { planes[i] = GetS 418 for (G4int i=0; i<4; ++i) { planes[i] = GetSidePlane(i); }
419 419
420 G4double xmin = kInfinity, xmax = -kInfinity 420 G4double xmin = kInfinity, xmax = -kInfinity;
421 G4double ymin = kInfinity, ymax = -kInfinity 421 G4double ymin = kInfinity, ymax = -kInfinity;
422 G4double dz = GetZHalfLength(); 422 G4double dz = GetZHalfLength();
423 for (G4int i=0; i<8; ++i) 423 for (G4int i=0; i<8; ++i)
424 { 424 {
425 G4int iy = (i==0 || i==1 || i==4 || i==5) 425 G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
426 G4int ix = (i==0 || i==2 || i==4 || i==6) 426 G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
427 G4double z = (i < 4) ? -dz : dz; 427 G4double z = (i < 4) ? -dz : dz;
428 G4double y = -(planes[iy].c*z + planes[iy] 428 G4double y = -(planes[iy].c*z + planes[iy].d)/planes[iy].b;
429 G4double x = -(planes[ix].b*y + planes[ix] 429 G4double x = -(planes[ix].b*y + planes[ix].c*z + planes[ix].d)/planes[ix].a;
430 if (x < xmin) xmin = x; 430 if (x < xmin) xmin = x;
431 if (x > xmax) xmax = x; 431 if (x > xmax) xmax = x;
432 if (y < ymin) ymin = y; 432 if (y < ymin) ymin = y;
433 if (y > ymax) ymax = y; 433 if (y > ymax) ymax = y;
434 } 434 }
435 435
436 pMin.set(xmin,ymin,-dz); 436 pMin.set(xmin,ymin,-dz);
437 pMax.set(xmax,ymax, dz); 437 pMax.set(xmax,ymax, dz);
438 438
439 // Check correctness of the bounding box 439 // Check correctness of the bounding box
440 // 440 //
441 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 441 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
442 { 442 {
443 std::ostringstream message; 443 std::ostringstream message;
444 message << "Bad bounding box (min >= max) 444 message << "Bad bounding box (min >= max) for solid: "
445 << GetName() << " !" 445 << GetName() << " !"
446 << "\npMin = " << pMin 446 << "\npMin = " << pMin
447 << "\npMax = " << pMax; 447 << "\npMax = " << pMax;
448 G4Exception("G4UTrap::BoundingLimits()", " 448 G4Exception("G4UTrap::BoundingLimits()", "GeomMgt0001",
449 JustWarning, message); 449 JustWarning, message);
450 StreamInfo(G4cout); 450 StreamInfo(G4cout);
451 } 451 }
452 452
453 // Check consistency of bounding boxes 453 // Check consistency of bounding boxes
454 // 454 //
455 if (checkBBox) 455 if (checkBBox)
456 { 456 {
457 G4double tolerance = kCarTolerance; 457 G4double tolerance = kCarTolerance;
458 U3Vector vmin, vmax; 458 U3Vector vmin, vmax;
459 Extent(vmin,vmax); 459 Extent(vmin,vmax);
460 if (std::abs(pMin.x()-vmin.x()) > toleranc 460 if (std::abs(pMin.x()-vmin.x()) > tolerance ||
461 std::abs(pMin.y()-vmin.y()) > toleranc 461 std::abs(pMin.y()-vmin.y()) > tolerance ||
462 std::abs(pMin.z()-vmin.z()) > toleranc 462 std::abs(pMin.z()-vmin.z()) > tolerance ||
463 std::abs(pMax.x()-vmax.x()) > toleranc 463 std::abs(pMax.x()-vmax.x()) > tolerance ||
464 std::abs(pMax.y()-vmax.y()) > toleranc 464 std::abs(pMax.y()-vmax.y()) > tolerance ||
465 std::abs(pMax.z()-vmax.z()) > toleranc 465 std::abs(pMax.z()-vmax.z()) > tolerance)
466 { 466 {
467 std::ostringstream message; 467 std::ostringstream message;
468 message << "Inconsistency in bounding bo 468 message << "Inconsistency in bounding boxes for solid: "
469 << GetName() << " !" 469 << GetName() << " !"
470 << "\nBBox min: wrapper = " << p 470 << "\nBBox min: wrapper = " << pMin << " solid = " << vmin
471 << "\nBBox max: wrapper = " << p 471 << "\nBBox max: wrapper = " << pMax << " solid = " << vmax;
472 G4Exception("G4UTrap::BoundingLimits()", 472 G4Exception("G4UTrap::BoundingLimits()", "GeomMgt0001",
473 JustWarning, message); 473 JustWarning, message);
474 checkBBox = false; 474 checkBBox = false;
475 } 475 }
476 } 476 }
477 } 477 }
478 478
479 ////////////////////////////////////////////// 479 //////////////////////////////////////////////////////////////////////////
480 // 480 //
481 // Calculate extent under transform and specif 481 // Calculate extent under transform and specified limit
482 482
483 G4bool 483 G4bool
484 G4UTrap::CalculateExtent(const EAxis pAxis, 484 G4UTrap::CalculateExtent(const EAxis pAxis,
485 const G4VoxelLimits& 485 const G4VoxelLimits& pVoxelLimit,
486 const G4AffineTransfo 486 const G4AffineTransform& pTransform,
487 G4double& pMin, 487 G4double& pMin, G4double& pMax) const
488 { 488 {
489 G4ThreeVector bmin, bmax; 489 G4ThreeVector bmin, bmax;
490 G4bool exist; 490 G4bool exist;
491 491
492 // Check bounding box (bbox) 492 // Check bounding box (bbox)
493 // 493 //
494 BoundingLimits(bmin,bmax); 494 BoundingLimits(bmin,bmax);
495 G4BoundingEnvelope bbox(bmin,bmax); 495 G4BoundingEnvelope bbox(bmin,bmax);
496 #ifdef G4BBOX_EXTENT 496 #ifdef G4BBOX_EXTENT
497 if (true) return bbox.CalculateExtent(pAxis, 497 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
498 #endif 498 #endif
499 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 499 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
500 { 500 {
501 return exist = pMin < pMax; 501 return exist = pMin < pMax;
502 } 502 }
503 503
504 // Set bounding envelope (benv) and calculat 504 // Set bounding envelope (benv) and calculate extent
505 // 505 //
506 TrapSidePlane planes[4]; 506 TrapSidePlane planes[4];
507 for (G4int i=0; i<4; ++i) { planes[i] = GetS 507 for (G4int i=0; i<4; ++i) { planes[i] = GetSidePlane(i); }
508 508
509 G4ThreeVector pt[8]; 509 G4ThreeVector pt[8];
510 G4double dz = GetZHalfLength(); 510 G4double dz = GetZHalfLength();
511 for (G4int i=0; i<8; ++i) 511 for (G4int i=0; i<8; ++i)
512 { 512 {
513 G4int iy = (i==0 || i==1 || i==4 || i==5) 513 G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
514 G4int ix = (i==0 || i==2 || i==4 || i==6) 514 G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
515 G4double z = (i < 4) ? -dz : dz; 515 G4double z = (i < 4) ? -dz : dz;
516 G4double y = -(planes[iy].c*z + planes[iy] 516 G4double y = -(planes[iy].c*z + planes[iy].d)/planes[iy].b;
517 G4double x = -(planes[ix].b*y + planes[ix] 517 G4double x = -(planes[ix].b*y + planes[ix].c*z + planes[ix].d)/planes[ix].a;
518 pt[i].set(x,y,z); 518 pt[i].set(x,y,z);
519 } 519 }
520 520
521 G4ThreeVectorList baseA(4), baseB(4); 521 G4ThreeVectorList baseA(4), baseB(4);
522 baseA[0] = pt[0]; 522 baseA[0] = pt[0];
523 baseA[1] = pt[1]; 523 baseA[1] = pt[1];
524 baseA[2] = pt[3]; 524 baseA[2] = pt[3];
525 baseA[3] = pt[2]; 525 baseA[3] = pt[2];
526 526
527 baseB[0] = pt[4]; 527 baseB[0] = pt[4];
528 baseB[1] = pt[5]; 528 baseB[1] = pt[5];
529 baseB[2] = pt[7]; 529 baseB[2] = pt[7];
530 baseB[3] = pt[6]; 530 baseB[3] = pt[6];
531 531
532 std::vector<const G4ThreeVectorList *> polyg 532 std::vector<const G4ThreeVectorList *> polygons(2);
533 polygons[0] = &baseA; 533 polygons[0] = &baseA;
534 polygons[1] = &baseB; 534 polygons[1] = &baseB;
535 535
536 G4BoundingEnvelope benv(bmin,bmax,polygons); 536 G4BoundingEnvelope benv(bmin,bmax,polygons);
537 exist = benv.CalculateExtent(pAxis,pVoxelLim 537 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
538 return exist; 538 return exist;
539 } 539 }
540 540
541 ////////////////////////////////////////////// 541 //////////////////////////////////////////////////////////////////////////
542 // 542 //
543 // Create polyhedron for visualization 543 // Create polyhedron for visualization
544 // 544 //
545 G4Polyhedron* G4UTrap::CreatePolyhedron() cons 545 G4Polyhedron* G4UTrap::CreatePolyhedron() const
546 { 546 {
547 return new G4PolyhedronTrap(GetZHalfLength() 547 return new G4PolyhedronTrap(GetZHalfLength(), GetTheta(), GetPhi(),
548 GetYHalfLength1( 548 GetYHalfLength1(),
549 GetXHalfLength1( 549 GetXHalfLength1(), GetXHalfLength2(), GetAlpha1(),
550 GetYHalfLength2( 550 GetYHalfLength2(),
551 GetXHalfLength3( 551 GetXHalfLength3(), GetXHalfLength4(), GetAlpha2());
552 } 552 }
553 553
554 #endif // G4GEOM_USE_USOLIDS 554 #endif // G4GEOM_USE_USOLIDS
555 555