| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Implementation for G4UTrap wrapper class
27 //
28 // 13.09.13 G.Cosmo, CERN
29 // --------------------------------------------------------------------
30
31 #include "G4Trap.hh"
32 #include "G4UTrap.hh"
33
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
35
36 #include "G4AffineTransform.hh"
37 #include "G4VPVParameterisation.hh"
38 #include "G4BoundingEnvelope.hh"
39
40 using namespace CLHEP;
41
42 /////////////////////////////////////////////////////////////////////////
43 //
44 // Constructors
45 //
46 G4UTrap::G4UTrap( const G4String& pName,
47 G4double pdz,
48 G4double pTheta, G4double pPhi,
49 G4double pdy1, G4double pdx1, G4double pdx2,
50 G4double pAlp1,
51 G4double pdy2, G4double pdx3, G4double pdx4,
52 G4double pAlp2 )
53 : Base_t(pName, pdz, pTheta, pPhi, pdy1, pdx1, pdx2,
54 pAlp1, pdy2, pdx3, pdx4, pAlp2)
55 {
56 G4ThreeVector pt[8];
57 CheckParameters();
58 GetVertices(pt);
59 CheckPlanarity(pt);
60 }
61
62 G4UTrap::G4UTrap( const G4String& pName,
63 const G4ThreeVector pt[8] )
64 : Base_t(pName)
65 {
66 // Start with check of centering - the center of gravity trap line
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() ) >= kCarTolerance
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].y()+pt[6].y()) >= kCarTolerance
86 || std::abs(pt[0].x()+pt[1].x()+pt[4].x()+pt[5].x() +
87 pt[2].x()+pt[3].x()+pt[6].x()+pt[7].x()) >= kCarTolerance )
88 {
89 std::ostringstream message;
90 message << "Invalid vertice coordinates for Solid: " << GetName();
91 G4Exception("G4UTrap::G4UTrap()", "GeomSolids0002",
92 FatalException, message);
93 }
94
95 SetPlanes(pt);
96 CheckParameters();
97 CheckPlanarity(pt);
98 }
99
100 // Constructor for Right Angular Wedge from STEP; this is different from
101 // the UnplacedTrapezoid constructor taking four double's and constructing
102 // a Trd1.
103 G4UTrap::G4UTrap( const G4String& pName,
104 G4double pZ,
105 G4double pY,
106 G4double pX, G4double pLTX )
107 : Base_t(pName, 0.5*pZ, /*theta=*/0, /*phi=*/0, /*dy1=*/0.5*pY,
108 /*dx1=*/0.5*pX, /*dx2=*/0.5*pLTX, /*Alpha1=*/0.5*(pLTX - pX)/pY,
109 /*dy2=*/0.5*pY, /*dx3=*/0.5*pX, /*dx4=*/0.5*pLTX,
110 /*Alpha2=*/0.5*(pLTX - pX)/pY)
111 {
112 CheckParameters();
113 }
114
115 G4UTrap::G4UTrap( const G4String& pName,
116 G4double pdx1, G4double pdx2,
117 G4double pdy1, G4double pdy2,
118 G4double pdz )
119 : Base_t(pName, pdx1, pdx2, pdy1, pdy2, pdz)
120 {
121 CheckParameters();
122 }
123
124 G4UTrap::G4UTrap(const G4String& pName,
125 G4double pdx, G4double pdy, G4double pdz,
126 G4double pAlpha, G4double pTheta, G4double pPhi )
127 : Base_t(pName, pdx, pdy, pdz, pAlpha, pTheta, pPhi)
128 {
129 CheckParameters();
130 }
131
132 G4UTrap::G4UTrap( const G4String& pName )
133 : Base_t(pName)
134 {
135 }
136
137 ///////////////////////////////////////////////////////////////////////
138 //
139 // Fake default constructor - sets only member data and allocates memory
140 // for usage restricted to object persistency.
141 //
142 G4UTrap::G4UTrap( __void__& a )
143 : Base_t(a)
144 {
145 }
146
147 //////////////////////////////////////////////////////////////////////////
148 //
149 // Destructor
150 //
151 G4UTrap::~G4UTrap() = default;
152
153 //////////////////////////////////////////////////////////////////////////
154 //
155 // Copy constructor
156 //
157 G4UTrap::G4UTrap(const G4UTrap& rhs)
158 : Base_t(rhs)
159 {
160 }
161
162 //////////////////////////////////////////////////////////////////////////
163 //
164 // Assignment operator
165 //
166 G4UTrap& G4UTrap::operator = (const G4UTrap& rhs)
167 {
168 // Check assignment to self
169 //
170 if (this == &rhs) { return *this; }
171
172 // Copy base class data
173 //
174 Base_t::operator=(rhs);
175
176 return *this;
177 }
178
179 //////////////////////////////////////////////////////////////////////////
180 //
181 // Accessors
182 //
183 G4double G4UTrap::GetZHalfLength() const
184 {
185 return GetDz();
186 }
187 G4double G4UTrap::GetYHalfLength1() const
188 {
189 return GetDy1();
190 }
191 G4double G4UTrap::GetXHalfLength1() const
192 {
193 return GetDx1();
194 }
195 G4double G4UTrap::GetXHalfLength2() const
196 {
197 return GetDx2();
198 }
199 G4double G4UTrap::GetTanAlpha1() const
200 {
201 return Base_t::GetTanAlpha1();
202 }
203 G4double G4UTrap::GetYHalfLength2() const
204 {
205 return GetDy2();
206 }
207 G4double G4UTrap::GetXHalfLength3() const
208 {
209 return GetDx3();
210 }
211 G4double G4UTrap::GetXHalfLength4() const
212 {
213 return GetDx4();
214 }
215 G4double G4UTrap::GetTanAlpha2() 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 {
229 return Base_t::GetAlpha1();
230 }
231 G4double G4UTrap::GetAlpha2() const
232 {
233 return Base_t::GetAlpha2();
234 }
235 TrapSidePlane G4UTrap::GetSidePlane(G4int n) const
236 {
237 TrapSidePlane plane;
238 plane.a = GetStruct().GetPlane(n).fA;
239 plane.b = GetStruct().GetPlane(n).fB;
240 plane.c = GetStruct().GetPlane(n).fC;
241 plane.d = GetStruct().GetPlane(n).fD;
242 return plane;
243 }
244 G4ThreeVector G4UTrap::GetSymAxis() const
245 {
246 G4double tanThetaSphi = GetTanThetaSinPhi();
247 G4double tanThetaCphi = GetTanThetaCosPhi();
248 G4double tan2Theta = tanThetaSphi*tanThetaSphi + tanThetaCphi*tanThetaCphi;
249 G4double cosTheta = 1.0 / std::sqrt(1 + tan2Theta);
250 return {tanThetaCphi*cosTheta, tanThetaSphi*cosTheta, cosTheta};
251 }
252
253 //////////////////////////////////////////////////////////////////////////
254 //
255 // Modifier
256 //
257 void G4UTrap::SetAllParameters(G4double pDz, G4double pTheta, G4double pPhi,
258 G4double pDy1, G4double pDx1, G4double pDx2,
259 G4double pAlp1,
260 G4double pDy2, G4double pDx3, G4double pDx4,
261 G4double pAlp2)
262 {
263 SetDz(pDz);
264 SetDy1(pDy1);
265 SetDy2(pDy2);
266 SetDx1(pDx1);
267 SetDx2(pDx2);
268 SetDx3(pDx3);
269 SetDx4(pDx4);
270 SetTanAlpha1(std::tan(pAlp1));
271 SetTanAlpha2(std::tan(pAlp2));
272 // last two will also reset cached variables
273 SetTheta(pTheta);
274 SetPhi(pPhi);
275 fRebuildPolyhedron = true;
276
277 G4ThreeVector pt[8];
278 CheckParameters();
279 GetVertices(pt);
280 CheckPlanarity(pt);
281 }
282
283 /////////////////////////////////////////////////////////////////////////
284 //
285 // Set parameters using eight vertices
286 //
287 void G4UTrap::SetPlanes(const G4ThreeVector pt[8])
288 {
289 U3Vector upt[8];
290 for (unsigned int i=0; i<8; ++i)
291 {
292 upt[i] = U3Vector(pt[i].x(), pt[i].y(), pt[i].z());
293 }
294 fromCornersToParameters(upt);
295 fRebuildPolyhedron = true;
296 }
297
298 /////////////////////////////////////////////////////////////////////////
299 //
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 Solid: " << GetName()
318 << "\n X - " <<fDx1<<", "<<fDx2<<", "<<fDx3<<", "<<fDx4
319 << "\n Y - " <<fDy1<<", "<<fDy2
320 << "\n Z - " <<fDz;
321 G4Exception("G4UTrap::CheckParameters()", "GeomSolids0002",
322 FatalException, message);
323 }
324 }
325
326 /////////////////////////////////////////////////////////////////////////
327 //
328 // Compute coordinates of vertices
329 //
330 void G4UTrap::GetVertices(G4ThreeVector pt[8]) const
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*GetTanThetaCosPhi();
343 G4double DzTthetaSphi = fDz*GetTanThetaSinPhi();
344 G4double Dy1Talpha1 = fDy1*fTalpha1;
345 G4double Dy2Talpha2 = fDy2*fTalpha2;
346
347 pt[0].set(-DzTthetaCphi-Dy1Talpha1-fDx1,-DzTthetaSphi-fDy1,-fDz);
348 pt[1].set(-DzTthetaCphi-Dy1Talpha1+fDx1,-DzTthetaSphi-fDy1,-fDz);
349 pt[2].set(-DzTthetaCphi+Dy1Talpha1-fDx2,-DzTthetaSphi+fDy1,-fDz);
350 pt[3].set(-DzTthetaCphi+Dy1Talpha1+fDx2,-DzTthetaSphi+fDy1,-fDz);
351 pt[4].set( DzTthetaCphi-Dy2Talpha2-fDx3, DzTthetaSphi-fDy2, fDz);
352 pt[5].set( DzTthetaCphi-Dy2Talpha2+fDx3, DzTthetaSphi-fDy2, fDz);
353 pt[6].set( DzTthetaCphi+Dy2Talpha2-fDx4, DzTthetaSphi+fDy2, fDz);
354 pt[7].set( DzTthetaCphi+Dy2Talpha2+fDx4, DzTthetaSphi+fDy2, fDz);
355 }
356
357 /////////////////////////////////////////////////////////////////////////
358 //
359 // Check planarity of lateral planes
360 //
361 void G4UTrap::CheckPlanarity(const G4ThreeVector pt[8]) const
362 {
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", "~+Y", "~-X", "~+X" };
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*p.y() + plane.c*p.z() + plane.d;
374 if (std::abs(dist) > std::abs(dmax)) dmax = dist;
375 }
376 if (std::abs(dmax) > 1000 * kCarTolerance)
377 {
378 std::ostringstream message;
379 message << "Side face " << side[i] << " is not planar for solid: "
380 << GetName() << "\nDiscrepancy: " << dmax/mm << " mm\n";
381 StreamInfo(message);
382 G4Exception("G4UTrap::CheckPlanarity()", "GeomSolids0002",
383 FatalException, message);
384 }
385 }
386 }
387
388 /////////////////////////////////////////////////////////////////////////
389 //
390 // Dispatch to parameterisation for replication mechanism dimension
391 // computation & modification.
392 //
393 void G4UTrap::ComputeDimensions( G4VPVParameterisation* p,
394 const G4int n,
395 const G4VPhysicalVolume* pRep)
396 {
397 p->ComputeDimensions(*(G4Trap*)this,n,pRep);
398 }
399
400 //////////////////////////////////////////////////////////////////////////
401 //
402 // Make a clone of the object
403 //
404 G4VSolid* G4UTrap::Clone() const
405 {
406 return new G4UTrap(*this);
407 }
408
409 //////////////////////////////////////////////////////////////////////////
410 //
411 // Get bounding box
412
413 void G4UTrap::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
414 {
415 static G4bool checkBBox = true;
416
417 TrapSidePlane planes[4];
418 for (G4int i=0; i<4; ++i) { planes[i] = GetSidePlane(i); }
419
420 G4double xmin = kInfinity, xmax = -kInfinity;
421 G4double ymin = kInfinity, ymax = -kInfinity;
422 G4double dz = GetZHalfLength();
423 for (G4int i=0; i<8; ++i)
424 {
425 G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
426 G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
427 G4double z = (i < 4) ? -dz : dz;
428 G4double y = -(planes[iy].c*z + planes[iy].d)/planes[iy].b;
429 G4double x = -(planes[ix].b*y + planes[ix].c*z + planes[ix].d)/planes[ix].a;
430 if (x < xmin) xmin = x;
431 if (x > xmax) xmax = x;
432 if (y < ymin) ymin = y;
433 if (y > ymax) ymax = y;
434 }
435
436 pMin.set(xmin,ymin,-dz);
437 pMax.set(xmax,ymax, dz);
438
439 // Check correctness of the bounding box
440 //
441 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
442 {
443 std::ostringstream message;
444 message << "Bad bounding box (min >= max) for solid: "
445 << GetName() << " !"
446 << "\npMin = " << pMin
447 << "\npMax = " << pMax;
448 G4Exception("G4UTrap::BoundingLimits()", "GeomMgt0001",
449 JustWarning, message);
450 StreamInfo(G4cout);
451 }
452
453 // Check consistency of bounding boxes
454 //
455 if (checkBBox)
456 {
457 G4double tolerance = kCarTolerance;
458 U3Vector vmin, vmax;
459 Extent(vmin,vmax);
460 if (std::abs(pMin.x()-vmin.x()) > tolerance ||
461 std::abs(pMin.y()-vmin.y()) > tolerance ||
462 std::abs(pMin.z()-vmin.z()) > tolerance ||
463 std::abs(pMax.x()-vmax.x()) > tolerance ||
464 std::abs(pMax.y()-vmax.y()) > tolerance ||
465 std::abs(pMax.z()-vmax.z()) > tolerance)
466 {
467 std::ostringstream message;
468 message << "Inconsistency in bounding boxes for solid: "
469 << GetName() << " !"
470 << "\nBBox min: wrapper = " << pMin << " solid = " << vmin
471 << "\nBBox max: wrapper = " << pMax << " solid = " << vmax;
472 G4Exception("G4UTrap::BoundingLimits()", "GeomMgt0001",
473 JustWarning, message);
474 checkBBox = false;
475 }
476 }
477 }
478
479 //////////////////////////////////////////////////////////////////////////
480 //
481 // Calculate extent under transform and specified limit
482
483 G4bool
484 G4UTrap::CalculateExtent(const EAxis pAxis,
485 const G4VoxelLimits& pVoxelLimit,
486 const G4AffineTransform& pTransform,
487 G4double& pMin, G4double& pMax) const
488 {
489 G4ThreeVector bmin, bmax;
490 G4bool exist;
491
492 // Check bounding box (bbox)
493 //
494 BoundingLimits(bmin,bmax);
495 G4BoundingEnvelope bbox(bmin,bmax);
496 #ifdef G4BBOX_EXTENT
497 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
498 #endif
499 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
500 {
501 return exist = pMin < pMax;
502 }
503
504 // Set bounding envelope (benv) and calculate extent
505 //
506 TrapSidePlane planes[4];
507 for (G4int i=0; i<4; ++i) { planes[i] = GetSidePlane(i); }
508
509 G4ThreeVector pt[8];
510 G4double dz = GetZHalfLength();
511 for (G4int i=0; i<8; ++i)
512 {
513 G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
514 G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
515 G4double z = (i < 4) ? -dz : dz;
516 G4double y = -(planes[iy].c*z + planes[iy].d)/planes[iy].b;
517 G4double x = -(planes[ix].b*y + planes[ix].c*z + planes[ix].d)/planes[ix].a;
518 pt[i].set(x,y,z);
519 }
520
521 G4ThreeVectorList baseA(4), baseB(4);
522 baseA[0] = pt[0];
523 baseA[1] = pt[1];
524 baseA[2] = pt[3];
525 baseA[3] = pt[2];
526
527 baseB[0] = pt[4];
528 baseB[1] = pt[5];
529 baseB[2] = pt[7];
530 baseB[3] = pt[6];
531
532 std::vector<const G4ThreeVectorList *> polygons(2);
533 polygons[0] = &baseA;
534 polygons[1] = &baseB;
535
536 G4BoundingEnvelope benv(bmin,bmax,polygons);
537 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
538 return exist;
539 }
540
541 //////////////////////////////////////////////////////////////////////////
542 //
543 // Create polyhedron for visualization
544 //
545 G4Polyhedron* G4UTrap::CreatePolyhedron() const
546 {
547 return new G4PolyhedronTrap(GetZHalfLength(), GetTheta(), GetPhi(),
548 GetYHalfLength1(),
549 GetXHalfLength1(), GetXHalfLength2(), GetAlpha1(),
550 GetYHalfLength2(),
551 GetXHalfLength3(), GetXHalfLength4(), GetAlpha2());
552 }
553
554 #endif // G4GEOM_USE_USOLIDS
555