ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UTubs.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4UTubs.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 G4UTubs wrapper class
27 //
28 // 30.10.13 G.Cosmo, CERN/PH
29 // --------------------------------------------------------------------
30 
31 #include "G4Tubs.hh"
32 #include "G4UTubs.hh"
33 
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
35 
36 #include "G4GeomTools.hh"
37 #include "G4AffineTransform.hh"
38 #include "G4VPVParameterisation.hh"
39 #include "G4BoundingEnvelope.hh"
40 
41 using namespace CLHEP;
42 
44 //
45 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
46 // - note if pdphi>2PI then reset to 2PI
47 
48 G4UTubs::G4UTubs( const G4String& pName,
49  G4double pRMin, G4double pRMax,
50  G4double pDz,
51  G4double pSPhi, G4double pDPhi )
52  : Base_t(pName, pRMin, pRMax, pDz, pSPhi, pDPhi)
53 {
54 }
55 
57 //
58 // Fake default constructor - sets only member data and allocates memory
59 // for usage restricted to object persistency.
60 //
61 G4UTubs::G4UTubs( __void__& a )
62  : Base_t(a)
63 {
64 }
65 
67 //
68 // Destructor
69 
70 G4UTubs::~G4UTubs()
71 {
72 }
73 
75 //
76 // Copy constructor
77 
78 G4UTubs::G4UTubs(const G4UTubs& rhs)
79  : Base_t(rhs)
80 {
81 }
82 
84 //
85 // Assignment operator
86 
87 G4UTubs& G4UTubs::operator = (const G4UTubs& rhs)
88 {
89  // Check assignment to self
90  //
91  if (this == &rhs) { return *this; }
92 
93  // Copy base class data
94  //
95  Base_t::operator=(rhs);
96 
97  return *this;
98 }
99 
101 //
102 // Accessors and modifiers
103 
104 G4double G4UTubs::GetInnerRadius() const
105 {
106  return rmin();
107 }
108 G4double G4UTubs::GetOuterRadius() const
109 {
110  return rmax();
111 }
112 G4double G4UTubs::GetZHalfLength() const
113 {
114  return z();
115 }
116 G4double G4UTubs::GetStartPhiAngle() const
117 {
118  return sphi();
119 }
120 G4double G4UTubs::GetDeltaPhiAngle() const
121 {
122  return dphi();
123 }
124 G4double G4UTubs::GetSinStartPhi() const
125 {
126  return std::sin(GetStartPhiAngle());
127 }
128 G4double G4UTubs::GetCosStartPhi() const
129 {
130  return std::cos(GetStartPhiAngle());
131 }
132 G4double G4UTubs::GetSinEndPhi() const
133 {
134  return std::sin(GetStartPhiAngle()+GetDeltaPhiAngle());
135 }
136 G4double G4UTubs::GetCosEndPhi() const
137 {
138  return std::cos(GetStartPhiAngle()+GetDeltaPhiAngle());
139 }
140 
141 void G4UTubs::SetInnerRadius(G4double newRMin)
142 {
143  SetRMin(newRMin);
144  fRebuildPolyhedron = true;
145 }
146 void G4UTubs::SetOuterRadius(G4double newRMax)
147 {
148  SetRMax(newRMax);
149  fRebuildPolyhedron = true;
150 }
151 void G4UTubs::SetZHalfLength(G4double newDz)
152 {
153  SetDz(newDz);
154  fRebuildPolyhedron = true;
155 }
156 void G4UTubs::SetStartPhiAngle(G4double newSPhi, G4bool)
157 {
158  SetSPhi(newSPhi);
159  fRebuildPolyhedron = true;
160 }
161 void G4UTubs::SetDeltaPhiAngle(G4double newDPhi)
162 {
163  SetDPhi(newDPhi);
164  fRebuildPolyhedron = true;
165 }
166 
168 //
169 // Dispatch to parameterisation for replication mechanism dimension
170 // computation & modification.
171 
172 void G4UTubs::ComputeDimensions( G4VPVParameterisation* p,
173  const G4int n,
174  const G4VPhysicalVolume* pRep )
175 {
176  p->ComputeDimensions(*(G4Tubs*)this,n,pRep) ;
177 }
178 
180 //
181 // Make a clone of the object
182 
183 G4VSolid* G4UTubs::Clone() const
184 {
185  return new G4UTubs(*this);
186 }
187 
189 //
190 // Get bounding box
191 
192 void G4UTubs::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
193 {
194  static G4bool checkBBox = true;
195 
196  G4double rmin = GetInnerRadius();
197  G4double rmax = GetOuterRadius();
198  G4double dz = GetZHalfLength();
199 
200  // Find bounding box
201  //
202  if (GetDeltaPhiAngle() < twopi)
203  {
204  G4TwoVector vmin,vmax;
205  G4GeomTools::DiskExtent(rmin,rmax,
206  GetSinStartPhi(),GetCosStartPhi(),
207  GetSinEndPhi(),GetCosEndPhi(),
208  vmin,vmax);
209  pMin.set(vmin.x(),vmin.y(),-dz);
210  pMax.set(vmax.x(),vmax.y(), dz);
211  }
212  else
213  {
214  pMin.set(-rmax,-rmax,-dz);
215  pMax.set( rmax, rmax, dz);
216  }
217 
218  // Check correctness of the bounding box
219  //
220  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
221  {
222  std::ostringstream message;
223  message << "Bad bounding box (min >= max) for solid: "
224  << GetName() << " !"
225  << "\npMin = " << pMin
226  << "\npMax = " << pMax;
227  G4Exception("G4UTubs::BoundingLimits()", "GeomMgt0001",
228  JustWarning, message);
229  StreamInfo(G4cout);
230  }
231 
232  // Check consistency of bounding boxes
233  //
234  if (checkBBox)
235  {
236  U3Vector vmin, vmax;
237  Extent(vmin,vmax);
238  if (std::abs(pMin.x()-vmin.x()) > kCarTolerance ||
239  std::abs(pMin.y()-vmin.y()) > kCarTolerance ||
240  std::abs(pMin.z()-vmin.z()) > kCarTolerance ||
241  std::abs(pMax.x()-vmax.x()) > kCarTolerance ||
242  std::abs(pMax.y()-vmax.y()) > kCarTolerance ||
243  std::abs(pMax.z()-vmax.z()) > kCarTolerance)
244  {
245  std::ostringstream message;
246  message << "Inconsistency in bounding boxes for solid: "
247  << GetName() << " !"
248  << "\nBBox min: wrapper = " << pMin << " solid = " << vmin
249  << "\nBBox max: wrapper = " << pMax << " solid = " << vmax;
250  G4Exception("G4UTubs::BoundingLimits()", "GeomMgt0001",
251  JustWarning, message);
252  checkBBox = false;
253  }
254  }
255 }
256 
258 //
259 // Calculate extent under transform and specified limit
260 
261 G4bool
262 G4UTubs::CalculateExtent(const EAxis pAxis,
263  const G4VoxelLimits& pVoxelLimit,
264  const G4AffineTransform& pTransform,
265  G4double& pMin, G4double& pMax) const
266 {
267  G4ThreeVector bmin, bmax;
268  G4bool exist;
269 
270  // Get bounding box
271  BoundingLimits(bmin,bmax);
272 
273  // Check bounding box
274  G4BoundingEnvelope bbox(bmin,bmax);
275 #ifdef G4BBOX_EXTENT
276  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
277 #endif
278  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
279  {
280  return exist = (pMin < pMax) ? true : false;
281  }
282 
283  // Get parameters of the solid
284  G4double rmin = GetInnerRadius();
285  G4double rmax = GetOuterRadius();
286  G4double dz = GetZHalfLength();
287  G4double dphi = GetDeltaPhiAngle();
288 
289  // Find bounding envelope and calculate extent
290  //
291  const G4int NSTEPS = 24; // number of steps for whole circle
292  G4double astep = twopi/NSTEPS; // max angle for one step
293  G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
294  G4double ang = dphi/ksteps;
295 
296  G4double sinHalf = std::sin(0.5*ang);
297  G4double cosHalf = std::cos(0.5*ang);
298  G4double sinStep = 2.*sinHalf*cosHalf;
299  G4double cosStep = 1. - 2.*sinHalf*sinHalf;
300  G4double rext = rmax/cosHalf;
301 
302  // bounding envelope for full cylinder consists of two polygons,
303  // in other cases it is a sequence of quadrilaterals
304  if (rmin == 0 && dphi == twopi)
305  {
306  G4double sinCur = sinHalf;
307  G4double cosCur = cosHalf;
308 
309  G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
310  for (G4int k=0; k<NSTEPS; ++k)
311  {
312  baseA[k].set(rext*cosCur,rext*sinCur,-dz);
313  baseB[k].set(rext*cosCur,rext*sinCur, dz);
314 
315  G4double sinTmp = sinCur;
316  sinCur = sinCur*cosStep + cosCur*sinStep;
317  cosCur = cosCur*cosStep - sinTmp*sinStep;
318  }
319  std::vector<const G4ThreeVectorList *> polygons(2);
320  polygons[0] = &baseA;
321  polygons[1] = &baseB;
322  G4BoundingEnvelope benv(bmin,bmax,polygons);
323  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
324  }
325  else
326  {
327  G4double sinStart = GetSinStartPhi();
328  G4double cosStart = GetCosStartPhi();
329  G4double sinEnd = GetSinEndPhi();
330  G4double cosEnd = GetCosEndPhi();
331  G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
332  G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
333 
334  // set quadrilaterals
335  G4ThreeVectorList pols[NSTEPS+2];
336  for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
337  pols[0][0].set(rmin*cosStart,rmin*sinStart, dz);
338  pols[0][1].set(rmin*cosStart,rmin*sinStart,-dz);
339  pols[0][2].set(rmax*cosStart,rmax*sinStart,-dz);
340  pols[0][3].set(rmax*cosStart,rmax*sinStart, dz);
341  for (G4int k=1; k<ksteps+1; ++k)
342  {
343  pols[k][0].set(rmin*cosCur,rmin*sinCur, dz);
344  pols[k][1].set(rmin*cosCur,rmin*sinCur,-dz);
345  pols[k][2].set(rext*cosCur,rext*sinCur,-dz);
346  pols[k][3].set(rext*cosCur,rext*sinCur, dz);
347 
348  G4double sinTmp = sinCur;
349  sinCur = sinCur*cosStep + cosCur*sinStep;
350  cosCur = cosCur*cosStep - sinTmp*sinStep;
351  }
352  pols[ksteps+1][0].set(rmin*cosEnd,rmin*sinEnd, dz);
353  pols[ksteps+1][1].set(rmin*cosEnd,rmin*sinEnd,-dz);
354  pols[ksteps+1][2].set(rmax*cosEnd,rmax*sinEnd,-dz);
355  pols[ksteps+1][3].set(rmax*cosEnd,rmax*sinEnd, dz);
356 
357  // set envelope and calculate extent
358  std::vector<const G4ThreeVectorList *> polygons;
359  polygons.resize(ksteps+2);
360  for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
361  G4BoundingEnvelope benv(bmin,bmax,polygons);
362  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
363  }
364  return exist;
365 }
366 
368 //
369 // Create polyhedron for visualization
370 //
371 G4Polyhedron* G4UTubs::CreatePolyhedron() const
372 {
373  return new G4PolyhedronTubs(GetInnerRadius(),
374  GetOuterRadius(),
375  GetZHalfLength(),
376  GetStartPhiAngle(),
377  GetDeltaPhiAngle());
378 }
379 
380 #endif // G4GEOM_USE_USOLIDS