ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UGenericPolycone.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4UGenericPolycone.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 G4UGenericPolycone wrapper class
27 //
28 // 30.10.13 G.Cosmo, CERN
29 // --------------------------------------------------------------------
30 
31 #include "G4GenericPolycone.hh"
32 #include "G4UGenericPolycone.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 #include "G4Polyhedron.hh"
42 
43 using namespace CLHEP;
44 
46 //
47 // Constructor (generic parameters)
48 //
49 G4UGenericPolycone::G4UGenericPolycone(const G4String& name,
50  G4double phiStart,
51  G4double phiTotal,
52  G4int numRZ,
53  const G4double r[],
54  const G4double z[] )
55  : Base_t(name, phiStart, phiTotal, numRZ, r, z)
56 {
57  wrStart = phiStart; while (wrStart < 0) wrStart += twopi;
58  wrDelta = phiTotal;
59  if (wrDelta <= 0 || wrDelta >= twopi*(1-DBL_EPSILON))
60  {
61  wrStart = 0;
62  wrDelta = twopi;
63  }
64  rzcorners.resize(0);
65  for (G4int i=0; i<numRZ; ++i)
66  {
67  rzcorners.push_back(G4TwoVector(r[i],z[i]));
68  }
69  std::vector<G4int> iout;
70  G4GeomTools::RemoveRedundantVertices(rzcorners,iout,2*kCarTolerance);
71 }
72 
73 
75 //
76 // Fake default constructor - sets only member data and allocates memory
77 // for usage restricted to object persistency.
78 //
79 G4UGenericPolycone::G4UGenericPolycone(__void__& a)
80  : Base_t(a)
81 {
82 }
83 
84 
86 //
87 // Destructor
88 //
89 G4UGenericPolycone::~G4UGenericPolycone()
90 {
91 }
92 
93 
95 //
96 // Copy constructor
97 //
98 G4UGenericPolycone::G4UGenericPolycone(const G4UGenericPolycone& source)
99  : Base_t(source)
100 {
101  wrStart = source.wrStart;
102  wrDelta = source.wrDelta;
103  rzcorners = source.rzcorners;
104 }
105 
106 
108 //
109 // Assignment operator
110 //
111 G4UGenericPolycone&
112 G4UGenericPolycone::operator=(const G4UGenericPolycone& source)
113 {
114  if (this == &source) return *this;
115 
116  Base_t::operator=( source );
117  wrStart = source.wrStart;
118  wrDelta = source.wrDelta;
119  rzcorners = source.rzcorners;
120 
121  return *this;
122 }
123 
124 G4double G4UGenericPolycone::GetStartPhi() const
125 {
126  return wrStart;
127 }
128 G4double G4UGenericPolycone::GetEndPhi() const
129 {
130  return (wrStart + wrDelta);
131 }
132 G4double G4UGenericPolycone::GetSinStartPhi() const
133 {
134  if (IsOpen()) return 0.;
135  G4double phi = GetStartPhi();
136  return std::sin(phi);
137 }
138 G4double G4UGenericPolycone::GetCosStartPhi() const
139 {
140  if (IsOpen()) return 1.;
141  G4double phi = GetStartPhi();
142  return std::cos(phi);
143 }
144 G4double G4UGenericPolycone::GetSinEndPhi() const
145 {
146  if (IsOpen()) return 0.;
147  G4double phi = GetEndPhi();
148  return std::sin(phi);
149 }
150 G4double G4UGenericPolycone::GetCosEndPhi() const
151 {
152  if (IsOpen()) return 1.;
153  G4double phi = GetEndPhi();
154  return std::cos(phi);
155 }
156 G4bool G4UGenericPolycone::IsOpen() const
157 {
158  return (wrDelta < twopi);
159 }
160 G4int G4UGenericPolycone::GetNumRZCorner() const
161 {
162  return rzcorners.size();
163 }
164 G4PolyconeSideRZ G4UGenericPolycone::GetCorner(G4int index) const
165 {
166  G4TwoVector rz = rzcorners.at(index);
167  G4PolyconeSideRZ psiderz = { rz.x(), rz.y() };
168 
169  return psiderz;
170 }
171 
173 //
174 // Make a clone of the object
175 
176 G4VSolid* G4UGenericPolycone::Clone() const
177 {
178  return new G4UGenericPolycone(*this);
179 }
180 
182 //
183 // Get bounding box
184 
185 void
186 G4UGenericPolycone::BoundingLimits(G4ThreeVector& pMin,
187  G4ThreeVector& pMax) const
188 {
189  G4double rmin = kInfinity, rmax = -kInfinity;
190  G4double zmin = kInfinity, zmax = -kInfinity;
191 
192  for (G4int i=0; i<GetNumRZCorner(); ++i)
193  {
194  G4PolyconeSideRZ corner = GetCorner(i);
195  if (corner.r < rmin) rmin = corner.r;
196  if (corner.r > rmax) rmax = corner.r;
197  if (corner.z < zmin) zmin = corner.z;
198  if (corner.z > zmax) zmax = corner.z;
199  }
200 
201  if (IsOpen())
202  {
203  G4TwoVector vmin,vmax;
204  G4GeomTools::DiskExtent(rmin,rmax,
205  GetSinStartPhi(),GetCosStartPhi(),
206  GetSinEndPhi(),GetCosEndPhi(),
207  vmin,vmax);
208  pMin.set(vmin.x(),vmin.y(),zmin);
209  pMax.set(vmax.x(),vmax.y(),zmax);
210  }
211  else
212  {
213  pMin.set(-rmax,-rmax, zmin);
214  pMax.set( rmax, rmax, zmax);
215  }
216 
217  // Check correctness of the bounding box
218  //
219  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
220  {
221  std::ostringstream message;
222  message << "Bad bounding box (min >= max) for solid: "
223  << GetName() << " !"
224  << "\npMin = " << pMin
225  << "\npMax = " << pMax;
226  G4Exception("G4UGenericPolycone::BoundingLimits()", "GeomMgt0001",
227  JustWarning, message);
228  StreamInfo(G4cout);
229  }
230 }
231 
233 //
234 // Calculate extent under transform and specified limit
235 
236 G4bool
237 G4UGenericPolycone::CalculateExtent(const EAxis pAxis,
238  const G4VoxelLimits& pVoxelLimit,
239  const G4AffineTransform& pTransform,
240  G4double& pMin, G4double& pMax) const
241 {
242  G4ThreeVector bmin, bmax;
243  G4bool exist;
244 
245  // Check bounding box (bbox)
246  //
247  BoundingLimits(bmin,bmax);
248  G4BoundingEnvelope bbox(bmin,bmax);
249 #ifdef G4BBOX_EXTENT
250  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
251 #endif
252  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
253  {
254  return exist = (pMin < pMax) ? true : false;
255  }
256 
257  // To find the extent, RZ contour of the polycone is subdivided
258  // in triangles. The extent is calculated as cumulative extent of
259  // all sub-polycones formed by rotation of triangles around Z
260  //
261  G4TwoVectorList contourRZ;
262  G4TwoVectorList triangles;
263  G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
264  G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
265 
266  // get RZ contour, ensure anticlockwise order of corners
267  for (G4int i=0; i<GetNumRZCorner(); ++i)
268  {
269  G4PolyconeSideRZ corner = GetCorner(i);
270  contourRZ.push_back(G4TwoVector(corner.r,corner.z));
271  }
272  G4double area = G4GeomTools::PolygonArea(contourRZ);
273  if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
274 
275  // triangulate RZ countour
276  if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
277  {
278  std::ostringstream message;
279  message << "Triangulation of RZ contour has failed for solid: "
280  << GetName() << " !"
281  << "\nExtent has been calculated using boundary box";
282  G4Exception("G4UGenericPolycone::CalculateExtent()",
283  "GeomMgt1002", JustWarning, message);
284  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
285  }
286 
287  // set trigonometric values
288  const G4int NSTEPS = 24; // number of steps for whole circle
289  G4double astep = twopi/NSTEPS; // max angle for one step
290 
291  G4double sphi = GetStartPhi();
292  G4double ephi = GetEndPhi();
293  G4double dphi = IsOpen() ? ephi-sphi : twopi;
294  G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
295  G4double ang = dphi/ksteps;
296 
297  G4double sinHalf = std::sin(0.5*ang);
298  G4double cosHalf = std::cos(0.5*ang);
299  G4double sinStep = 2.*sinHalf*cosHalf;
300  G4double cosStep = 1. - 2.*sinHalf*sinHalf;
301 
302  G4double sinStart = GetSinStartPhi();
303  G4double cosStart = GetCosStartPhi();
304  G4double sinEnd = GetSinEndPhi();
305  G4double cosEnd = GetCosEndPhi();
306 
307  // define vectors and arrays
308  std::vector<const G4ThreeVectorList *> polygons;
309  polygons.resize(ksteps+2);
310  G4ThreeVectorList pols[NSTEPS+2];
311  for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(6);
312  for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
313  G4double r0[6],z0[6]; // contour with original edges of triangle
314  G4double r1[6]; // shifted radii of external edges of triangle
315 
316  // main loop along triangles
317  pMin = kInfinity;
318  pMax =-kInfinity;
319  G4int ntria = triangles.size()/3;
320  for (G4int i=0; i<ntria; ++i)
321  {
322  G4int i3 = i*3;
323  for (G4int k=0; k<3; ++k)
324  {
325  G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3;
326  G4int k2 = k*2;
327  // set contour with original edges of triangle
328  r0[k2+0] = triangles[e0].x(); z0[k2+0] = triangles[e0].y();
329  r0[k2+1] = triangles[e1].x(); z0[k2+1] = triangles[e1].y();
330  // set shifted radii
331  r1[k2+0] = r0[k2+0];
332  r1[k2+1] = r0[k2+1];
333  if (z0[k2+1] - z0[k2+0] <= 0) continue;
334  r1[k2+0] /= cosHalf;
335  r1[k2+1] /= cosHalf;
336  }
337 
338  // rotate countour, set sequence of 6-sided polygons
339  G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
340  G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
341  for (G4int j=0; j<6; ++j)
342  {
343  pols[0][j].set(r0[j]*cosStart,r0[j]*sinStart,z0[j]);
344  }
345  for (G4int k=1; k<ksteps+1; ++k)
346  {
347  for (G4int j=0; j<6; ++j)
348  {
349  pols[k][j].set(r1[j]*cosCur,r1[j]*sinCur,z0[j]);
350  }
351  G4double sinTmp = sinCur;
352  sinCur = sinCur*cosStep + cosCur*sinStep;
353  cosCur = cosCur*cosStep - sinTmp*sinStep;
354  }
355  for (G4int j=0; j<6; ++j)
356  {
357  pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j]*sinEnd,z0[j]);
358  }
359 
360  // set sub-envelope and adjust extent
361  G4double emin,emax;
362  G4BoundingEnvelope benv(polygons);
363  if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
364  if (emin < pMin) pMin = emin;
365  if (emax > pMax) pMax = emax;
366  if (eminlim > pMin && emaxlim < pMax) return true; // max possible extent
367  }
368  return (pMin < pMax);
369 }
370 
372 //
373 // CreatePolyhedron
374 
375 G4Polyhedron* G4UGenericPolycone::CreatePolyhedron() const
376 {
377 
378 
379  // The following code prepares for:
380  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
381  // const double xyz[][3],
382  // const int faces_vec[][4])
383  // Here is an extract from the header file HepPolyhedron.h:
401  const G4int numSide =
402  G4int(G4Polyhedron::GetNumberOfRotationSteps()
403  * (GetEndPhi() - GetStartPhi()) / twopi) + 1;
404  G4int nNodes;
405  G4int nFaces;
406  typedef G4double double3[3];
407  double3* xyz;
408  typedef G4int int4[4];
409  int4* faces_vec;
410  if (IsOpen())
411  {
412  // Triangulate open ends. Simple ear-chopping algorithm...
413  // I'm not sure how robust this algorithm is (J.Allison).
414  //
415  std::vector<G4bool> chopped(GetNumRZCorner(), false);
416  std::vector<G4int*> triQuads;
417  G4int remaining = GetNumRZCorner();
418  G4int iStarter = 0;
419  while (remaining >= 3) // Loop checking, 13.08.2015, G.Cosmo
420  {
421  // Find unchopped corners...
422  //
423  G4int A = -1, B = -1, C = -1;
424  G4int iStepper = iStarter;
425  do // Loop checking, 13.08.2015, G.Cosmo
426  {
427  if (A < 0) { A = iStepper; }
428  else if (B < 0) { B = iStepper; }
429  else if (C < 0) { C = iStepper; }
430  do // Loop checking, 13.08.2015, G.Cosmo
431  {
432  if (++iStepper >= GetNumRZCorner()) { iStepper = 0; }
433  }
434  while (chopped[iStepper]);
435  }
436  while (C < 0 && iStepper != iStarter);
437 
438  // Check triangle at B is pointing outward (an "ear").
439  // Sign of z cross product determines...
440  //
441  G4double BAr = GetCorner(A).r - GetCorner(B).r;
442  G4double BAz = GetCorner(A).z - GetCorner(B).z;
443  G4double BCr = GetCorner(C).r - GetCorner(B).r;
444  G4double BCz = GetCorner(C).z - GetCorner(B).z;
445  if (BAr * BCz - BAz * BCr < kCarTolerance)
446  {
447  G4int* tq = new G4int[3];
448  tq[0] = A + 1;
449  tq[1] = B + 1;
450  tq[2] = C + 1;
451  triQuads.push_back(tq);
452  chopped[B] = true;
453  --remaining;
454  }
455  else
456  {
457  do // Loop checking, 13.08.2015, G.Cosmo
458  {
459  if (++iStarter >= GetNumRZCorner()) { iStarter = 0; }
460  }
461  while (chopped[iStarter]);
462  }
463  }
464  // Transfer to faces...
465  //
466  nNodes = (numSide + 1) * GetNumRZCorner();
467  nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size();
468  faces_vec = new int4[nFaces];
469  G4int iface = 0;
470  G4int addition = GetNumRZCorner() * numSide;
471  G4int d = GetNumRZCorner() - 1;
472  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
473  {
474  for (size_t i = 0; i < triQuads.size(); ++i)
475  {
476  // Negative for soft/auxiliary/normally invisible edges...
477  //
478  G4int a, b, c;
479  if (iEnd == 0)
480  {
481  a = triQuads[i][0];
482  b = triQuads[i][1];
483  c = triQuads[i][2];
484  }
485  else
486  {
487  a = triQuads[i][0] + addition;
488  b = triQuads[i][2] + addition;
489  c = triQuads[i][1] + addition;
490  }
491  G4int ab = std::abs(b - a);
492  G4int bc = std::abs(c - b);
493  G4int ca = std::abs(a - c);
494  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
495  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
496  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
497  faces_vec[iface][3] = 0;
498  ++iface;
499  }
500  }
501 
502  // Continue with sides...
503 
504  xyz = new double3[nNodes];
505  const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
506  G4double phi = GetStartPhi();
507  G4int ixyz = 0;
508  for (G4int iSide = 0; iSide < numSide; ++iSide)
509  {
510  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
511  {
512  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
513  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
514  xyz[ixyz][2] = GetCorner(iCorner).z;
515  if (iSide == 0) // startPhi
516  {
517  if (iCorner < GetNumRZCorner() - 1)
518  {
519  faces_vec[iface][0] = ixyz + 1;
520  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
521  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
522  faces_vec[iface][3] = ixyz + 2;
523  }
524  else
525  {
526  faces_vec[iface][0] = ixyz + 1;
527  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
528  faces_vec[iface][2] = ixyz + 2;
529  faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
530  }
531  }
532  else if (iSide == numSide - 1) // endPhi
533  {
534  if (iCorner < GetNumRZCorner() - 1)
535  {
536  faces_vec[iface][0] = ixyz + 1;
537  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
538  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
539  faces_vec[iface][3] = -(ixyz + 2);
540  }
541  else
542  {
543  faces_vec[iface][0] = ixyz + 1;
544  faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
545  faces_vec[iface][2] = ixyz + 2;
546  faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
547  }
548  }
549  else
550  {
551  if (iCorner < GetNumRZCorner() - 1)
552  {
553  faces_vec[iface][0] = ixyz + 1;
554  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
555  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
556  faces_vec[iface][3] = -(ixyz + 2);
557  }
558  else
559  {
560  faces_vec[iface][0] = ixyz + 1;
561  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
562  faces_vec[iface][2] = ixyz + 2;
563  faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
564  }
565  }
566  ++iface;
567  ++ixyz;
568  }
569  phi += dPhi;
570  }
571 
572  // Last corners...
573 
574  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
575  {
576  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
577  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
578  xyz[ixyz][2] = GetCorner(iCorner).z;
579  ++ixyz;
580  }
581  }
582  else // !phiIsOpen - i.e., a complete 360 degrees.
583  {
584  nNodes = numSide * GetNumRZCorner();
585  nFaces = numSide * GetNumRZCorner();;
586  xyz = new double3[nNodes];
587  faces_vec = new int4[nFaces];
588  const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
589  G4double phi = GetStartPhi();
590  G4int ixyz = 0, iface = 0;
591  for (G4int iSide = 0; iSide < numSide; ++iSide)
592  {
593  for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
594  {
595  xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
596  xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
597  xyz[ixyz][2] = GetCorner(iCorner).z;
598 
599  if (iSide < numSide - 1)
600  {
601  if (iCorner < GetNumRZCorner() - 1)
602  {
603  faces_vec[iface][0] = ixyz + 1;
604  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
605  faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
606  faces_vec[iface][3] = -(ixyz + 2);
607  }
608  else
609  {
610  faces_vec[iface][0] = ixyz + 1;
611  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() + 1);
612  faces_vec[iface][2] = ixyz + 2;
613  faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
614  }
615  }
616  else // Last side joins ends...
617  {
618  if (iCorner < GetNumRZCorner() - 1)
619  {
620  faces_vec[iface][0] = ixyz + 1;
621  faces_vec[iface][1] = -(ixyz + GetNumRZCorner() - nFaces + 1);
622  faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2;
623  faces_vec[iface][3] = -(ixyz + 2);
624  }
625  else
626  {
627  faces_vec[iface][0] = ixyz + 1;
628  faces_vec[iface][1] = -(ixyz - nFaces + GetNumRZCorner() + 1);
629  faces_vec[iface][2] = ixyz - nFaces + 2;
630  faces_vec[iface][3] = -(ixyz - GetNumRZCorner() + 2);
631  }
632  }
633  ++ixyz;
634  ++iface;
635  }
636  phi += dPhi;
637  }
638  }
639  G4Polyhedron* polyhedron = new G4Polyhedron;
640  G4int prob = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
641  delete [] faces_vec;
642  delete [] xyz;
643  if (prob)
644  {
645  std::ostringstream message;
646  message << "Problem creating G4Polyhedron for: " << GetName();
647  G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002",
648  JustWarning, message);
649  delete polyhedron;
650  return nullptr;
651  }
652  else
653  {
654  return polyhedron;
655  }
656 }
657 
658 #endif // G4GEOM_USE_USOLIDS