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