ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UTrap.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4UTrap.cc
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/PH
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 
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 }
57 
58 G4UTrap::G4UTrap( const G4String& pName,
59  const G4ThreeVector pt[8] )
60  : Base_t(pName)
61 {
62  SetPlanes(pt);
63 }
64 
65 G4UTrap::G4UTrap( const G4String& pName,
66  G4double pZ,
67  G4double pY,
68  G4double pX, G4double pLTX )
69  : Base_t(pName, pZ, pY, pX, pLTX)
70 {
71 }
72 
73 G4UTrap::G4UTrap( const G4String& pName,
74  G4double pdx1, G4double pdx2,
75  G4double pdy1, G4double pdy2,
76  G4double pdz )
77  : Base_t(pName, pdx1, pdx2, pdy1, pdy2, pdz)
78 {
79 }
80 
81 G4UTrap::G4UTrap(const G4String& pName,
82  G4double pdx, G4double pdy, G4double pdz,
83  G4double pAlpha, G4double pTheta, G4double pPhi )
84  : Base_t(pName, pdx, pdy, pdz, pAlpha, pTheta, pPhi)
85 {
86 }
87 
88 G4UTrap::G4UTrap( const G4String& pName )
89  : Base_t(pName)
90 {
91 }
92 
94 //
95 // Fake default constructor - sets only member data and allocates memory
96 // for usage restricted to object persistency.
97 //
98 G4UTrap::G4UTrap( __void__& a )
99  : Base_t(a)
100 {
101 }
102 
104 //
105 // Destructor
106 //
107 G4UTrap::~G4UTrap()
108 {
109 }
110 
112 //
113 // Copy constructor
114 //
115 G4UTrap::G4UTrap(const G4UTrap& rhs)
116  : Base_t(rhs)
117 {
118 }
119 
121 //
122 // Assignment operator
123 //
124 G4UTrap& G4UTrap::operator = (const G4UTrap& rhs)
125 {
126  // Check assignment to self
127  //
128  if (this == &rhs) { return *this; }
129 
130  // Copy base class data
131  //
132  Base_t::operator=(rhs);
133 
134  return *this;
135 }
136 
138 //
139 // Accessors & modifiers
140 
141 G4double G4UTrap::GetZHalfLength() const
142 {
143  return GetDz();
144 }
145 G4double G4UTrap::GetYHalfLength1() const
146 {
147  return GetDy1();
148 }
149 G4double G4UTrap::GetXHalfLength1() const
150 {
151  return GetDx1();
152 }
153 G4double G4UTrap::GetXHalfLength2() const
154 {
155  return GetDx2();
156 }
157 G4double G4UTrap::GetYHalfLength2() const
158 {
159  return GetDy2();
160 }
161 G4double G4UTrap::GetXHalfLength3() const
162 {
163  return GetDx3();
164 }
165 G4double G4UTrap::GetXHalfLength4() const
166 {
167  return GetDx4();
168 }
169 G4double G4UTrap::GetThetaCphi() const
170 {
171  return GetTanThetaCosPhi();
172 }
173 G4double G4UTrap::GetThetaSphi() const
174 {
175  return GetTanThetaSinPhi();
176 }
177 TrapSidePlane G4UTrap::GetSidePlane(G4int n) const
178 {
179  TrapSidePlane plane;
180  plane.a = GetStruct().GetPlane(n).fA;
181  plane.b = GetStruct().GetPlane(n).fB;
182  plane.c = GetStruct().GetPlane(n).fC;
183  plane.d = GetStruct().GetPlane(n).fD;
184  return plane;
185 }
186 G4ThreeVector G4UTrap::GetSymAxis() const
187 {
188  G4double tanThetaSphi = GetTanThetaSinPhi();
189  G4double tanThetaCphi = GetTanThetaCosPhi();
190  G4double tan2Theta = tanThetaSphi*tanThetaSphi + tanThetaCphi*tanThetaCphi;
191  G4double cosTheta = 1.0 / std::sqrt(1 + tan2Theta);
192  return G4ThreeVector(tanThetaCphi*cosTheta, tanThetaSphi*cosTheta, cosTheta);
193 }
194 
195 void G4UTrap::SetAllParameters(G4double pDz, G4double pTheta, G4double pPhi,
196  G4double pDy1, G4double pDx1, G4double pDx2,
197  G4double pAlp1,
198  G4double pDy2, G4double pDx3, G4double pDx4,
199  G4double pAlp2)
200 {
201  SetDz(pDz);
202  SetDy1(pDy1);
203  SetDy2(pDy2);
204  SetDx1(pDx1);
205  SetDx2(pDx2);
206  SetDx3(pDx3);
207  SetDx4(pDx4);
208  SetTanAlpha1(std::tan(pAlp1));
209  SetTanAlpha1(std::tan(pAlp2));
210  // last two will also reset cached variables
211  SetTheta(pTheta);
212  SetPhi(pPhi);
213  fRebuildPolyhedron = true;
214 }
215 
216 void G4UTrap::SetPlanes(const G4ThreeVector pt[8])
217 {
218  U3Vector upt[8];
219  for (unsigned int i=0; i<8; ++i)
220  {
221  upt[i] = U3Vector(pt[i].x(), pt[i].y(), pt[i].z());
222  }
223  fromCornersToParameters(upt);
224  fRebuildPolyhedron = true;
225 }
226 
228 //
229 // Dispatch to parameterisation for replication mechanism dimension
230 // computation & modification.
231 //
232 void G4UTrap::ComputeDimensions( G4VPVParameterisation* p,
233  const G4int n,
234  const G4VPhysicalVolume* pRep)
235 {
236  p->ComputeDimensions(*(G4Trap*)this,n,pRep);
237 }
238 
240 //
241 // Make a clone of the object
242 //
243 G4VSolid* G4UTrap::Clone() const
244 {
245  return new G4UTrap(*this);
246 }
247 
249 //
250 // Get bounding box
251 
252 void G4UTrap::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
253 {
254  static G4bool checkBBox = true;
255 
256  TrapSidePlane planes[4];
257  for (G4int i=0; i<4; ++i) { planes[i] = GetSidePlane(i); }
258 
259  G4double xmin = kInfinity, xmax = -kInfinity;
260  G4double ymin = kInfinity, ymax = -kInfinity;
261  G4double dz = GetZHalfLength();
262  for (G4int i=0; i<8; ++i)
263  {
264  G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
265  G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
266  G4double z = (i < 4) ? -dz : dz;
267  G4double y = -(planes[iy].c*z + planes[iy].d)/planes[iy].b;
268  G4double x = -(planes[ix].b*y + planes[ix].c*z + planes[ix].d)/planes[ix].a;
269  if (x < xmin) xmin = x;
270  if (x > xmax) xmax = x;
271  if (y < ymin) ymin = y;
272  if (y > ymax) ymax = y;
273  }
274 
275  pMin.set(xmin,ymin,-dz);
276  pMax.set(xmax,ymax, dz);
277 
278  // Check correctness of the bounding box
279  //
280  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
281  {
282  std::ostringstream message;
283  message << "Bad bounding box (min >= max) for solid: "
284  << GetName() << " !"
285  << "\npMin = " << pMin
286  << "\npMax = " << pMax;
287  G4Exception("G4UTrap::BoundingLimits()", "GeomMgt0001",
288  JustWarning, message);
289  StreamInfo(G4cout);
290  }
291 
292  // Check consistency of bounding boxes
293  //
294  if (checkBBox)
295  {
296  G4double tolerance = kCarTolerance;
297  U3Vector vmin, vmax;
298  Extent(vmin,vmax);
299  if (std::abs(pMin.x()-vmin.x()) > tolerance ||
300  std::abs(pMin.y()-vmin.y()) > tolerance ||
301  std::abs(pMin.z()-vmin.z()) > tolerance ||
302  std::abs(pMax.x()-vmax.x()) > tolerance ||
303  std::abs(pMax.y()-vmax.y()) > tolerance ||
304  std::abs(pMax.z()-vmax.z()) > tolerance)
305  {
306  std::ostringstream message;
307  message << "Inconsistency in bounding boxes for solid: "
308  << GetName() << " !"
309  << "\nBBox min: wrapper = " << pMin << " solid = " << vmin
310  << "\nBBox max: wrapper = " << pMax << " solid = " << vmax;
311  G4Exception("G4UTrap::BoundingLimits()", "GeomMgt0001",
312  JustWarning, message);
313  checkBBox = false;
314  }
315  }
316 }
317 
319 //
320 // Calculate extent under transform and specified limit
321 
322 G4bool
323 G4UTrap::CalculateExtent(const EAxis pAxis,
324  const G4VoxelLimits& pVoxelLimit,
325  const G4AffineTransform& pTransform,
326  G4double& pMin, G4double& pMax) const
327 {
328  G4ThreeVector bmin, bmax;
329  G4bool exist;
330 
331  // Check bounding box (bbox)
332  //
333  BoundingLimits(bmin,bmax);
334  G4BoundingEnvelope bbox(bmin,bmax);
335 #ifdef G4BBOX_EXTENT
336  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
337 #endif
338  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
339  {
340  return exist = (pMin < pMax) ? true : false;
341  }
342 
343  // Set bounding envelope (benv) and calculate extent
344  //
345  TrapSidePlane planes[4];
346  for (G4int i=0; i<4; ++i) { planes[i] = GetSidePlane(i); }
347 
348  G4ThreeVector pt[8];
349  G4double dz = GetZHalfLength();
350  for (G4int i=0; i<8; ++i)
351  {
352  G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1;
353  G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3;
354  G4double z = (i < 4) ? -dz : dz;
355  G4double y = -(planes[iy].c*z + planes[iy].d)/planes[iy].b;
356  G4double x = -(planes[ix].b*y + planes[ix].c*z + planes[ix].d)/planes[ix].a;
357  pt[i].set(x,y,z);
358  }
359 
360  G4ThreeVectorList baseA(4), baseB(4);
361  baseA[0] = pt[0];
362  baseA[1] = pt[1];
363  baseA[2] = pt[3];
364  baseA[3] = pt[2];
365 
366  baseB[0] = pt[4];
367  baseB[1] = pt[5];
368  baseB[2] = pt[7];
369  baseB[3] = pt[6];
370 
371  std::vector<const G4ThreeVectorList *> polygons(2);
372  polygons[0] = &baseA;
373  polygons[1] = &baseB;
374 
375  G4BoundingEnvelope benv(bmin,bmax,polygons);
376  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
377  return exist;
378 }
379 
381 //
382 // Create polyhedron for visualization
383 //
384 G4Polyhedron* G4UTrap::CreatePolyhedron() const
385 {
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,
394  GetYHalfLength1(),
395  GetXHalfLength1(), GetXHalfLength2(), alpha1,
396  GetYHalfLength2(),
397  GetXHalfLength3(), GetXHalfLength4(), alpha2);
398 }
399 
400 #endif // G4GEOM_USE_USOLIDS