Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4UTrap.cc
Go to the documentation of this file.
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
40using namespace CLHEP;
41
42/////////////////////////////////////////////////////////////////////////
43//
44// Constructors
45//
46G4UTrap::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
62G4UTrap::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.
103G4UTrap::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
115G4UTrap::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
124G4UTrap::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
132G4UTrap::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//
142G4UTrap::G4UTrap( __void__& a )
143 : Base_t(a)
144{
145}
146
147//////////////////////////////////////////////////////////////////////////
148//
149// Destructor
150//
151G4UTrap::~G4UTrap() = default;
152
153//////////////////////////////////////////////////////////////////////////
154//
155// Copy constructor
156//
157G4UTrap::G4UTrap(const G4UTrap& rhs)
158 : Base_t(rhs)
159{
160}
161
162//////////////////////////////////////////////////////////////////////////
163//
164// Assignment operator
165//
166G4UTrap& 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//
183G4double G4UTrap::GetZHalfLength() const
184{
185 return GetDz();
186}
187G4double G4UTrap::GetYHalfLength1() const
188{
189 return GetDy1();
190}
191G4double G4UTrap::GetXHalfLength1() const
192{
193 return GetDx1();
194}
195G4double G4UTrap::GetXHalfLength2() const
196{
197 return GetDx2();
198}
199G4double G4UTrap::GetTanAlpha1() const
200{
201 return Base_t::GetTanAlpha1();
202}
203G4double G4UTrap::GetYHalfLength2() const
204{
205 return GetDy2();
206}
207G4double G4UTrap::GetXHalfLength3() const
208{
209 return GetDx3();
210}
211G4double G4UTrap::GetXHalfLength4() const
212{
213 return GetDx4();
214}
215G4double G4UTrap::GetTanAlpha2() const
216{
217 return Base_t::GetTanAlpha2();
218}
219G4double G4UTrap::GetPhi() const
220{
221 return Base_t::GetPhi();
222}
223G4double G4UTrap::GetTheta() const
224{
225 return Base_t::GetTheta();
226}
227G4double G4UTrap::GetAlpha1() const
228{
229 return Base_t::GetAlpha1();
230}
231G4double G4UTrap::GetAlpha2() const
232{
233 return Base_t::GetAlpha2();
234}
235TrapSidePlane 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}
244G4ThreeVector 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//
257void 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//
287void 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//
302void 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//
330void 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//
361void 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//
393void 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//
404G4VSolid* G4UTrap::Clone() const
405{
406 return new G4UTrap(*this);
407}
408
409//////////////////////////////////////////////////////////////////////////
410//
411// Get bounding box
412
413void 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
483G4bool
484G4UTrap::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//
545G4Polyhedron* 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
kCarTolerance
const G4double kCarTolerance
Definition G4BoundingEnvelope.cc:38
G4ThreeVectorList
std::vector< G4ThreeVector > G4ThreeVectorList
Definition G4BoundingEnvelope.hh:54
JustWarning
@ JustWarning
Definition G4ExceptionSeverity.hh:71
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4VPVParameterisation::ComputeDimensions
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition G4VPVParameterisation.hh:89
EAxis
EAxis
Definition geomdefs.hh:54
TrapSidePlane::b
G4double b
Definition G4Trap.hh:92
TrapSidePlane::c
G4double c
Definition G4Trap.hh:92
TrapSidePlane::d
G4double d
Definition G4Trap.hh:92
TrapSidePlane::a
G4double a
Definition G4Trap.hh:92