ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UPolycone.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4UPolycone.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 of G4UPolycone wrapper class
27 //
28 // 31.10.13 G.Cosmo, CERN
29 // --------------------------------------------------------------------
30 
31 #include "G4Polycone.hh"
32 #include "G4UPolycone.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 (GEANT3 style parameters)
46 //
47 G4UPolycone::G4UPolycone( const G4String& name,
48  G4double phiStart,
49  G4double phiTotal,
50  G4int numZPlanes,
51  const G4double zPlane[],
52  const G4double rInner[],
53  const G4double rOuter[] )
54  : Base_t(name, phiStart, phiTotal, numZPlanes, zPlane, rInner, rOuter)
55 {
56  fGenericPcon = false;
57  SetOriginalParameters();
58  wrStart = phiStart;
59  while (wrStart < 0)
60  {
61  wrStart += twopi;
62  }
63  wrDelta = phiTotal;
64  if (wrDelta <= 0 || wrDelta >= twopi*(1-DBL_EPSILON))
65  {
66  wrStart = 0;
67  wrDelta = twopi;
68  }
69  rzcorners.resize(0);
70  for (G4int i=0; i<numZPlanes; ++i)
71  {
72  G4double z = zPlane[i];
73  G4double r = rOuter[i];
74  rzcorners.push_back(G4TwoVector(r,z));
75  }
76  for (G4int i=numZPlanes-1; i>=0; --i)
77  {
78  G4double z = zPlane[i];
79  G4double r = rInner[i];
80  rzcorners.push_back(G4TwoVector(r,z));
81  }
82  std::vector<G4int> iout;
83  G4GeomTools::RemoveRedundantVertices(rzcorners,iout,2*kCarTolerance);
84 }
85 
86 
88 //
89 // Constructor (generic parameters)
90 //
91 G4UPolycone::G4UPolycone(const G4String& name,
92  G4double phiStart,
93  G4double phiTotal,
94  G4int numRZ,
95  const G4double r[],
96  const G4double z[] )
97  : Base_t(name, phiStart, phiTotal, numRZ, r, z)
98 {
99  fGenericPcon = true;
100  SetOriginalParameters();
101  wrStart = phiStart; while (wrStart < 0) wrStart += twopi;
102  wrDelta = phiTotal;
103  if (wrDelta <= 0 || wrDelta >= twopi*(1-DBL_EPSILON))
104  {
105  wrStart = 0;
106  wrDelta = twopi;
107  }
108  rzcorners.resize(0);
109  for (G4int i=0; i<numRZ; ++i)
110  {
111  rzcorners.push_back(G4TwoVector(r[i],z[i]));
112  }
113  std::vector<G4int> iout;
114  G4GeomTools::RemoveRedundantVertices(rzcorners,iout,2*kCarTolerance);
115 }
116 
117 
119 //
120 // Fake default constructor - sets only member data and allocates memory
121 // for usage restricted to object persistency.
122 //
123 G4UPolycone::G4UPolycone( __void__& a )
124  : Base_t(a)
125 {
126 }
127 
128 
130 //
131 // Destructor
132 //
133 G4UPolycone::~G4UPolycone()
134 {
135 }
136 
137 
139 //
140 // Copy constructor
141 //
142 G4UPolycone::G4UPolycone( const G4UPolycone& source )
143  : Base_t( source )
144 {
145  fGenericPcon = source.fGenericPcon;
146  fOriginalParameters = source.fOriginalParameters;
147  wrStart = source.wrStart;
148  wrDelta = source.wrDelta;
149  rzcorners = source.rzcorners;
150 }
151 
152 
154 //
155 // Assignment operator
156 //
157 G4UPolycone& G4UPolycone::operator=( const G4UPolycone& source )
158 {
159  if (this == &source) return *this;
160 
161  Base_t::operator=( source );
162  fGenericPcon = source.fGenericPcon;
163  fOriginalParameters = source.fOriginalParameters;
164  wrStart = source.wrStart;
165  wrDelta = source.wrDelta;
166  rzcorners = source.rzcorners;
167 
168  return *this;
169 }
170 
171 
173 //
174 // Accessors & modifiers
175 //
176 G4double G4UPolycone::GetStartPhi() const
177 {
178  return wrStart;
179 }
180 G4double G4UPolycone::GetDeltaPhi() const
181 {
182  return wrDelta;
183 }
184 G4double G4UPolycone::GetEndPhi() const
185 {
186  return (wrStart + wrDelta);
187 }
188 G4double G4UPolycone::GetSinStartPhi() const
189 {
190  if (!IsOpen()) return 0.;
191  G4double phi = GetStartPhi();
192  return std::sin(phi);
193 }
194 G4double G4UPolycone::GetCosStartPhi() const
195 {
196  if (!IsOpen()) return 1.;
197  G4double phi = GetStartPhi();
198  return std::cos(phi);
199 }
200 G4double G4UPolycone::GetSinEndPhi() const
201 {
202  if (!IsOpen()) return 0.;
203  G4double phi = GetEndPhi();
204  return std::sin(phi);
205 }
206 G4double G4UPolycone::GetCosEndPhi() const
207 {
208  if (!IsOpen()) return 1.;
209  G4double phi = GetEndPhi();
210  return std::cos(phi);
211 }
212 G4bool G4UPolycone::IsOpen() const
213 {
214  return (wrDelta < twopi);
215 }
216 G4int G4UPolycone::GetNumRZCorner() const
217 {
218  return rzcorners.size();
219 }
220 G4PolyconeSideRZ G4UPolycone::GetCorner(G4int index) const
221 {
222  G4TwoVector rz = rzcorners.at(index);
223  G4PolyconeSideRZ psiderz = { rz.x(), rz.y() };
224 
225  return psiderz;
226 }
227 G4PolyconeHistorical* G4UPolycone::GetOriginalParameters() const
228 {
229  return new G4PolyconeHistorical(fOriginalParameters);
230 }
231 void G4UPolycone::SetOriginalParameters()
232 {
233  vecgeom::PolyconeHistorical* original_parameters = Base_t::GetOriginalParameters();
234 
235  fOriginalParameters.Start_angle = original_parameters->fHStart_angle;
236  fOriginalParameters.Opening_angle = original_parameters->fHOpening_angle;
237  fOriginalParameters.Num_z_planes = original_parameters->fHNum_z_planes;
238 
239  delete [] fOriginalParameters.Z_values;
240  delete [] fOriginalParameters.Rmin;
241  delete [] fOriginalParameters.Rmax;
242 
243  G4int numPlanes = fOriginalParameters.Num_z_planes;
244  fOriginalParameters.Z_values = new G4double[numPlanes];
245  fOriginalParameters.Rmin = new G4double[numPlanes];
246  fOriginalParameters.Rmax = new G4double[numPlanes];
247  for (G4int i=0; i<numPlanes; ++i)
248  {
249  fOriginalParameters.Z_values[i] = original_parameters->fHZ_values[i];
250  fOriginalParameters.Rmin[i] = original_parameters->fHRmin[i];
251  fOriginalParameters.Rmax[i] = original_parameters->fHRmax[i];
252  }
253 }
254 void G4UPolycone::SetOriginalParameters(G4PolyconeHistorical* pars)
255 {
256  fOriginalParameters = *pars;
257  fRebuildPolyhedron = true;
258  Reset();
259 }
260 
261 G4bool G4UPolycone::Reset()
262 {
263  if (fGenericPcon)
264  {
265  std::ostringstream message;
266  message << "Solid " << GetName() << " built using generic construct."
267  << G4endl << "Not applicable to the generic construct !";
268  G4Exception("G4UPolycone::Reset()", "GeomSolids1001",
269  JustWarning, message, "Parameters NOT resetted.");
270  return true; // error code set
271  }
272 
273  //
274  // Rebuild polycone based on original parameters
275  //
276  wrStart = fOriginalParameters.Start_angle;
277  while (wrStart < 0.)
278  {
279  wrStart += twopi;
280  }
281  wrDelta = fOriginalParameters.Opening_angle;
282  if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL_EPSILON))
283  {
284  wrStart = 0.;
285  wrDelta = twopi;
286  }
287  rzcorners.resize(0);
288  for (G4int i=0; i<fOriginalParameters.Num_z_planes; ++i)
289  {
290  G4double z = fOriginalParameters.Z_values[i];
291  G4double r = fOriginalParameters.Rmax[i];
292  rzcorners.push_back(G4TwoVector(r,z));
293  }
294  for (G4int i=fOriginalParameters.Num_z_planes-1; i>=0; --i)
295  {
296  G4double z = fOriginalParameters.Z_values[i];
297  G4double r = fOriginalParameters.Rmin[i];
298  rzcorners.push_back(G4TwoVector(r,z));
299  }
300  std::vector<G4int> iout;
301  G4GeomTools::RemoveRedundantVertices(rzcorners,iout,2*kCarTolerance);
302 
303  return false; // error code unset
304 }
305 
307 //
308 // Dispatch to parameterisation for replication mechanism dimension
309 // computation & modification.
310 //
311 void G4UPolycone::ComputeDimensions(G4VPVParameterisation* p,
312  const G4int n,
313  const G4VPhysicalVolume* pRep)
314 {
315  p->ComputeDimensions(*(G4Polycone*)this,n,pRep);
316 }
317 
318 
320 //
321 // Make a clone of the object
322 
323 G4VSolid* G4UPolycone::Clone() const
324 {
325  return new G4UPolycone(*this);
326 }
327 
329 //
330 // Get bounding box
331 
332 void G4UPolycone::BoundingLimits(G4ThreeVector& pMin,
333  G4ThreeVector& pMax) const
334 {
335  static G4bool checkBBox = true;
336  static G4bool checkPhi = true;
337 
338  G4double rmin = kInfinity, rmax = -kInfinity;
339  G4double zmin = kInfinity, zmax = -kInfinity;
340 
341  for (G4int i=0; i<GetNumRZCorner(); ++i)
342  {
343  G4PolyconeSideRZ corner = GetCorner(i);
344  if (corner.r < rmin) rmin = corner.r;
345  if (corner.r > rmax) rmax = corner.r;
346  if (corner.z < zmin) zmin = corner.z;
347  if (corner.z > zmax) zmax = corner.z;
348  }
349 
350  if (IsOpen())
351  {
352  G4TwoVector vmin,vmax;
353  G4GeomTools::DiskExtent(rmin,rmax,
354  GetSinStartPhi(),GetCosStartPhi(),
355  GetSinEndPhi(),GetCosEndPhi(),
356  vmin,vmax);
357  pMin.set(vmin.x(),vmin.y(),zmin);
358  pMax.set(vmax.x(),vmax.y(),zmax);
359  }
360  else
361  {
362  pMin.set(-rmax,-rmax, zmin);
363  pMax.set( rmax, rmax, zmax);
364  }
365 
366  // Check correctness of the bounding box
367  //
368  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
369  {
370  std::ostringstream message;
371  message << "Bad bounding box (min >= max) for solid: "
372  << GetName() << " !"
373  << "\npMin = " << pMin
374  << "\npMax = " << pMax;
375  G4Exception("G4UPolycone::BoundingLimits()", "GeomMgt0001",
376  JustWarning, message);
377  StreamInfo(G4cout);
378  }
379 
380  // Check consistency of bounding boxes
381  //
382  if (checkBBox)
383  {
384  U3Vector vmin, vmax;
385  Extent(vmin,vmax);
386  if (std::abs(pMin.x()-vmin.x()) > kCarTolerance ||
387  std::abs(pMin.y()-vmin.y()) > kCarTolerance ||
388  std::abs(pMin.z()-vmin.z()) > kCarTolerance ||
389  std::abs(pMax.x()-vmax.x()) > kCarTolerance ||
390  std::abs(pMax.y()-vmax.y()) > kCarTolerance ||
391  std::abs(pMax.z()-vmax.z()) > kCarTolerance)
392  {
393  std::ostringstream message;
394  message << "Inconsistency in bounding boxes for solid: "
395  << GetName() << " !"
396  << "\nBBox min: wrapper = " << pMin << " solid = " << vmin
397  << "\nBBox max: wrapper = " << pMax << " solid = " << vmax;
398  G4Exception("G4UPolycone::BoundingLimits()", "GeomMgt0001",
399  JustWarning, message);
400  checkBBox = false;
401  }
402  }
403 
404  // Check consistency of angles
405  //
406  if (checkPhi)
407  {
408  if (GetStartPhi() != Base_t::GetStartPhi() ||
409  GetEndPhi() != Base_t::GetEndPhi() ||
410  IsOpen() != (Base_t::GetDeltaPhi() < twopi))
411  {
412  std::ostringstream message;
413  message << "Inconsistency in Phi angles or # of sides for solid: "
414  << GetName() << " !"
415  << "\nPhi start : wrapper = " << GetStartPhi()
416  << " solid = " << Base_t::GetStartPhi()
417  << "\nPhi end : wrapper = " << GetEndPhi()
418  << " solid = " << Base_t::GetEndPhi()
419  << "\nPhi is open: wrapper = " << (IsOpen() ? "true" : "false")
420  << " solid = "
421  << ((Base_t::GetDeltaPhi() < twopi) ? "true" : "false");
422  G4Exception("G4UPolycone::BoundingLimits()", "GeomMgt0001",
423  JustWarning, message);
424  checkPhi = false;
425  }
426  }
427 }
428 
430 //
431 // Calculate extent under transform and specified limit
432 
433 G4bool G4UPolycone::CalculateExtent(const EAxis pAxis,
434  const G4VoxelLimits& pVoxelLimit,
435  const G4AffineTransform& pTransform,
436  G4double& pMin, G4double& pMax) const
437 {
438  G4ThreeVector bmin, bmax;
439  G4bool exist;
440 
441  // Check bounding box (bbox)
442  //
443  BoundingLimits(bmin,bmax);
444  G4BoundingEnvelope bbox(bmin,bmax);
445 #ifdef G4BBOX_EXTENT
446  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
447 #endif
448  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
449  {
450  return exist = (pMin < pMax) ? true : false;
451  }
452 
453  // To find the extent, RZ contour of the polycone is subdivided
454  // in triangles. The extent is calculated as cumulative extent of
455  // all sub-polycones formed by rotation of triangles around Z
456  //
457  G4TwoVectorList contourRZ;
458  G4TwoVectorList triangles;
459  std::vector<G4int> iout;
460  G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
461  G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
462 
463  // get RZ contour, ensure anticlockwise order of corners
464  for (G4int i=0; i<GetNumRZCorner(); ++i)
465  {
466  G4PolyconeSideRZ corner = GetCorner(i);
467  contourRZ.push_back(G4TwoVector(corner.r,corner.z));
468  }
469  G4GeomTools::RemoveRedundantVertices(contourRZ,iout,2*kCarTolerance);
470  G4double area = G4GeomTools::PolygonArea(contourRZ);
471  if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
472 
473  // triangulate RZ countour
474  if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
475  {
476  std::ostringstream message;
477  message << "Triangulation of RZ contour has failed for solid: "
478  << GetName() << " !"
479  << "\nExtent has been calculated using boundary box";
480  G4Exception("G4UPolycone::CalculateExtent()",
481  "GeomMgt1002", JustWarning, message);
482  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
483  }
484 
485  // set trigonometric values
486  const G4int NSTEPS = 24; // number of steps for whole circle
487  G4double astep = twopi/NSTEPS; // max angle for one step
488 
489  G4double sphi = GetStartPhi();
490  G4double ephi = GetEndPhi();
491  G4double dphi = IsOpen() ? ephi-sphi : twopi;
492  G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
493  G4double ang = dphi/ksteps;
494 
495  G4double sinHalf = std::sin(0.5*ang);
496  G4double cosHalf = std::cos(0.5*ang);
497  G4double sinStep = 2.*sinHalf*cosHalf;
498  G4double cosStep = 1. - 2.*sinHalf*sinHalf;
499 
500  G4double sinStart = GetSinStartPhi();
501  G4double cosStart = GetCosStartPhi();
502  G4double sinEnd = GetSinEndPhi();
503  G4double cosEnd = GetCosEndPhi();
504 
505  // define vectors and arrays
506  std::vector<const G4ThreeVectorList *> polygons;
507  polygons.resize(ksteps+2);
508  G4ThreeVectorList pols[NSTEPS+2];
509  for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(6);
510  for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
511  G4double r0[6],z0[6]; // contour with original edges of triangle
512  G4double r1[6]; // shifted radii of external edges of triangle
513 
514  // main loop along triangles
515  pMin = kInfinity;
516  pMax =-kInfinity;
517  G4int ntria = triangles.size()/3;
518  for (G4int i=0; i<ntria; ++i)
519  {
520  G4int i3 = i*3;
521  for (G4int k=0; k<3; ++k)
522  {
523  G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3;
524  G4int k2 = k*2;
525  // set contour with original edges of triangle
526  r0[k2+0] = triangles[e0].x(); z0[k2+0] = triangles[e0].y();
527  r0[k2+1] = triangles[e1].x(); z0[k2+1] = triangles[e1].y();
528  // set shifted radii
529  r1[k2+0] = r0[k2+0];
530  r1[k2+1] = r0[k2+1];
531  if (z0[k2+1] - z0[k2+0] <= 0) continue;
532  r1[k2+0] /= cosHalf;
533  r1[k2+1] /= cosHalf;
534  }
535 
536  // rotate countour, set sequence of 6-sided polygons
537  G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
538  G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
539  for (G4int j=0; j<6; ++j) pols[0][j].set(r0[j]*cosStart,r0[j]*sinStart,z0[j]);
540  for (G4int k=1; k<ksteps+1; ++k)
541  {
542  for (G4int j=0; j<6; ++j) pols[k][j].set(r1[j]*cosCur,r1[j]*sinCur,z0[j]);
543  G4double sinTmp = sinCur;
544  sinCur = sinCur*cosStep + cosCur*sinStep;
545  cosCur = cosCur*cosStep - sinTmp*sinStep;
546  }
547  for (G4int j=0; j<6; ++j) pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j]*sinEnd,z0[j]);
548 
549  // set sub-envelope and adjust extent
550  G4double emin,emax;
551  G4BoundingEnvelope benv(polygons);
552  if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
553  if (emin < pMin) pMin = emin;
554  if (emax > pMax) pMax = emax;
555  if (eminlim > pMin && emaxlim < pMax) return true; // max possible extent
556  }
557  return (pMin < pMax);
558 }
559 
561 //
562 // CreatePolyhedron
563 //
564 G4Polyhedron* G4UPolycone::CreatePolyhedron() const
565 {
566  G4PolyhedronPcon*
567  polyhedron = new G4PolyhedronPcon( fOriginalParameters.Start_angle,
568  fOriginalParameters.Opening_angle,
569  fOriginalParameters.Num_z_planes,
570  fOriginalParameters.Z_values,
571  fOriginalParameters.Rmin,
572  fOriginalParameters.Rmax );
573  return polyhedron;
574 }
575 
576 #endif // G4GEOM_USE_USOLIDS