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