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