ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GenericPolycone.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GenericPolycone.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 // G4GenericPolycone implementation
27 //
28 // Authors: T.Nikitina, G.Cosmo - CERN
29 // --------------------------------------------------------------------
30 
31 #include "G4GenericPolycone.hh"
32 
33 #if !defined(G4GEOM_USE_UGENERICPOLYCONE)
34 
35 #include "G4PolyconeSide.hh"
36 #include "G4PolyPhiFace.hh"
37 
38 #include "G4GeomTools.hh"
39 #include "G4VoxelLimits.hh"
40 #include "G4AffineTransform.hh"
41 #include "G4BoundingEnvelope.hh"
42 
43 #include "Randomize.hh"
44 
45 #include "G4Polyhedron.hh"
46 #include "G4EnclosingCylinder.hh"
47 #include "G4ReduciblePolygon.hh"
48 #include "G4VPVParameterisation.hh"
49 
50 using namespace CLHEP;
51 
52 // Constructor (generic parameters)
53 //
55  G4double phiStart,
56  G4double phiTotal,
57  G4int numRZ,
58  const G4double r[],
59  const G4double z[] )
60  : G4VCSGfaceted( name )
61 {
62 
63  G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
64 
65  Create( phiStart, phiTotal, rz );
66 
67  // Set original_parameters struct for consistency
68  //
69  //SetOriginalParameters(rz);
70 
71  delete rz;
72 }
73 
74 // Create
75 //
76 // Generic create routine, called by each constructor after
77 // conversion of arguments
78 //
80  G4double phiTotal,
81  G4ReduciblePolygon *rz )
82 {
83  //
84  // Perform checks of rz values
85  //
86  if (rz->Amin() < 0.0)
87  {
88  std::ostringstream message;
89  message << "Illegal input parameters - " << GetName() << G4endl
90  << " All R values must be >= 0 !";
91  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
92  FatalErrorInArgument, message);
93  }
94 
95  G4double rzArea = rz->Area();
96  if (rzArea < -kCarTolerance)
97  {
98  rz->ReverseOrder();
99  }
100  else if (rzArea < kCarTolerance)
101  {
102  std::ostringstream message;
103  message << "Illegal input parameters - " << GetName() << G4endl
104  << " R/Z cross section is zero or near zero: " << rzArea;
105  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
106  FatalErrorInArgument, message);
107  }
108 
110  || (!rz->RemoveRedundantVertices( kCarTolerance )) )
111  {
112  std::ostringstream message;
113  message << "Illegal input parameters - " << GetName() << G4endl
114  << " Too few unique R/Z values !";
115  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
116  FatalErrorInArgument, message);
117  }
118 
119  if (rz->CrossesItself(1/kInfinity))
120  {
121  std::ostringstream message;
122  message << "Illegal input parameters - " << GetName() << G4endl
123  << " R/Z segments cross !";
124  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
125  FatalErrorInArgument, message);
126  }
127 
128  numCorner = rz->NumVertices();
129 
130  //
131  // Phi opening? Account for some possible roundoff, and interpret
132  // nonsense value as representing no phi opening
133  //
134  if (phiTotal <= 0 || phiTotal > twopi-1E-10)
135  {
136  phiIsOpen = false;
137  startPhi = 0;
138  endPhi = twopi;
139  }
140  else
141  {
142  phiIsOpen = true;
143 
144  //
145  // Convert phi into our convention
146  //
147  startPhi = phiStart;
148  while( startPhi < 0 ) // Loop checking, 13.08.2015, G.Cosmo
149  startPhi += twopi;
150 
151  endPhi = phiStart+phiTotal;
152  while( endPhi < startPhi ) // Loop checking, 13.08.2015, G.Cosmo
153  endPhi += twopi;
154  }
155 
156  //
157  // Allocate corner array.
158  //
160 
161  //
162  // Copy corners
163  //
164  G4ReduciblePolygonIterator iterRZ(rz);
165 
166  G4PolyconeSideRZ* next = corners;
167  iterRZ.Begin();
168  do // Loop checking, 13.08.2015, G.Cosmo
169  {
170  next->r = iterRZ.GetA();
171  next->z = iterRZ.GetB();
172  } while( ++next, iterRZ.Next() );
173 
174  //
175  // Allocate face pointer array
176  //
178  faces = new G4VCSGface*[numFace];
179 
180  //
181  // Construct conical faces
182  //
183  // But! Don't construct a face if both points are at zero radius!
184  //
185  G4PolyconeSideRZ *corner = corners,
186  *prev = corners + numCorner-1,
187  *nextNext;
188  G4VCSGface** face = faces;
189  do // Loop checking, 13.08.2015, G.Cosmo
190  {
191  next = corner+1;
192  if (next >= corners+numCorner) next = corners;
193  nextNext = next+1;
194  if (nextNext >= corners+numCorner) nextNext = corners;
195 
196  if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
197 
198  //
199  // We must decide here if we can dare declare one of our faces
200  // as having a "valid" normal (i.e. allBehind = true). This
201  // is never possible if the face faces "inward" in r.
202  //
203  G4bool allBehind;
204  if (corner->z > next->z)
205  {
206  allBehind = false;
207  }
208  else
209  {
210  //
211  // Otherwise, it is only true if the line passing
212  // through the two points of the segment do not
213  // split the r/z cross section
214  //
215  allBehind = !rz->BisectedBy( corner->r, corner->z,
216  next->r, next->z, kCarTolerance );
217  }
218 
219  *face++ = new G4PolyconeSide( prev, corner, next, nextNext,
220  startPhi, endPhi-startPhi, phiIsOpen, allBehind );
221  } while( prev=corner, corner=next, corner > corners );
222 
223  if (phiIsOpen)
224  {
225  //
226  // Construct phi open edges
227  //
228  *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi );
229  *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi );
230  }
231 
232  //
233  // We might have dropped a face or two: recalculate numFace
234  //
235  numFace = face-faces;
236 
237  //
238  // Make enclosingCylinder
239  //
241  new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
242 }
243 
244 // Fake default constructor - sets only member data and allocates memory
245 // for usage restricted to object persistency.
246 //
248  : G4VCSGfaceted(a), startPhi(0.), endPhi(0.), numCorner(0)
249 {
250 }
251 
252 // Destructor
253 //
255 {
256  delete [] corners;
257  delete enclosingCylinder;
258 }
259 
260 // Copy constructor
261 //
263  : G4VCSGfaceted( source )
264 {
265  CopyStuff( source );
266 }
267 
268 // Assignment operator
269 //
272 {
273  if (this == &source) return *this;
274 
275  G4VCSGfaceted::operator=( source );
276 
277  delete [] corners;
278  // if (original_parameters) delete original_parameters;
279 
280  delete enclosingCylinder;
281 
282  CopyStuff( source );
283 
284  return *this;
285 }
286 
287 // CopyStuff
288 //
290 {
291  //
292  // Simple stuff
293  //
294  startPhi = source.startPhi;
295  endPhi = source.endPhi;
296  phiIsOpen = source.phiIsOpen;
297  numCorner = source.numCorner;
298 
299  //
300  // The corner array
301  //
303 
304  G4PolyconeSideRZ *corn = corners,
305  *sourceCorn = source.corners;
306  do // Loop checking, 13.08.2015, G.Cosmo
307  {
308  *corn = *sourceCorn;
309  } while( ++sourceCorn, ++corn < corners+numCorner );
310 
311  //
312  // Enclosing cylinder
313  //
315 
316  fRebuildPolyhedron = false;
317  fpPolyhedron = nullptr;
318 }
319 
320 // Reset
321 //
323 {
324  std::ostringstream message;
325  message << "Solid " << GetName() << " built using generic construct."
326  << G4endl << "Not applicable to the generic construct !";
327  G4Exception("G4GenericPolycone::Reset()", "GeomSolids1001",
328  JustWarning, message, "Parameters NOT resetted.");
329  return true;
330 }
331 
332 // Inside
333 //
334 // This is an override of G4VCSGfaceted::Inside, created in order
335 // to speed things up by first checking with G4EnclosingCylinder.
336 //
338 {
339  //
340  // Quick test
341  //
342  if (enclosingCylinder->MustBeOutside(p)) return kOutside;
343 
344  //
345  // Long answer
346  //
347  return G4VCSGfaceted::Inside(p);
348 }
349 
350 // DistanceToIn
351 //
352 // This is an override of G4VCSGfaceted::Inside, created in order
353 // to speed things up by first checking with G4EnclosingCylinder.
354 //
356  const G4ThreeVector& v ) const
357 {
358  //
359  // Quick test
360  //
361  if (enclosingCylinder->ShouldMiss(p,v))
362  return kInfinity;
363 
364  //
365  // Long answer
366  //
367  return G4VCSGfaceted::DistanceToIn( p, v );
368 }
369 
370 // DistanceToIn
371 //
373 {
374  return G4VCSGfaceted::DistanceToIn(p);
375 }
376 
377 // BoundingLimits
378 //
379 // Get bounding box
380 //
381 void
383  G4ThreeVector& pMax) const
384 {
385  G4double rmin = kInfinity, rmax = -kInfinity;
386  G4double zmin = kInfinity, zmax = -kInfinity;
387 
388  for (G4int i=0; i<GetNumRZCorner(); ++i)
389  {
390  G4PolyconeSideRZ corner = GetCorner(i);
391  if (corner.r < rmin) rmin = corner.r;
392  if (corner.r > rmax) rmax = corner.r;
393  if (corner.z < zmin) zmin = corner.z;
394  if (corner.z > zmax) zmax = corner.z;
395  }
396 
397  if (IsOpen())
398  {
399  G4TwoVector vmin,vmax;
400  G4GeomTools::DiskExtent(rmin,rmax,
403  vmin,vmax);
404  pMin.set(vmin.x(),vmin.y(),zmin);
405  pMax.set(vmax.x(),vmax.y(),zmax);
406  }
407  else
408  {
409  pMin.set(-rmax,-rmax, zmin);
410  pMax.set( rmax, rmax, zmax);
411  }
412 
413  // Check correctness of the bounding box
414  //
415  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
416  {
417  std::ostringstream message;
418  message << "Bad bounding box (min >= max) for solid: "
419  << GetName() << " !"
420  << "\npMin = " << pMin
421  << "\npMax = " << pMax;
422  G4Exception("GenericG4Polycone::BoundingLimits()", "GeomMgt0001",
423  JustWarning, message);
424  DumpInfo();
425  }
426 }
427 
428 // CalculateExtent
429 //
430 // Calculate extent under transform and specified limit
431 //
432 G4bool
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  G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
460  G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
461 
462  // get RZ contour, ensure anticlockwise order of corners
463  for (G4int i=0; i<GetNumRZCorner(); ++i)
464  {
465  G4PolyconeSideRZ corner = GetCorner(i);
466  contourRZ.push_back(G4TwoVector(corner.r,corner.z));
467  }
468  G4double area = G4GeomTools::PolygonArea(contourRZ);
469  if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
470 
471  // triangulate RZ countour
472  if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
473  {
474  std::ostringstream message;
475  message << "Triangulation of RZ contour has failed for solid: "
476  << GetName() << " !"
477  << "\nExtent has been calculated using boundary box";
478  G4Exception("G4GenericPolycone::CalculateExtent()",
479  "GeomMgt1002", JustWarning, message);
480  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
481  }
482 
483  // set trigonometric values
484  const G4int NSTEPS = 24; // number of steps for whole circle
485  G4double astep = twopi/NSTEPS; // max angle for one step
486 
487  G4double sphi = GetStartPhi();
488  G4double ephi = GetEndPhi();
489  G4double dphi = IsOpen() ? ephi-sphi : twopi;
490  G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
491  G4double ang = dphi/ksteps;
492 
493  G4double sinHalf = std::sin(0.5*ang);
494  G4double cosHalf = std::cos(0.5*ang);
495  G4double sinStep = 2.*sinHalf*cosHalf;
496  G4double cosStep = 1. - 2.*sinHalf*sinHalf;
497 
498  G4double sinStart = GetSinStartPhi();
499  G4double cosStart = GetCosStartPhi();
500  G4double sinEnd = GetSinEndPhi();
501  G4double cosEnd = GetCosEndPhi();
502 
503  // define vectors and arrays
504  std::vector<const G4ThreeVectorList *> polygons;
505  polygons.resize(ksteps+2);
506  G4ThreeVectorList pols[NSTEPS+2];
507  for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(6);
508  for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
509  G4double r0[6],z0[6]; // contour with original edges of triangle
510  G4double r1[6]; // shifted radii of external edges of triangle
511 
512  // main loop along triangles
513  pMin = kInfinity;
514  pMax =-kInfinity;
515  G4int ntria = triangles.size()/3;
516  for (G4int i=0; i<ntria; ++i)
517  {
518  G4int i3 = i*3;
519  for (G4int k=0; k<3; ++k)
520  {
521  G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3;
522  G4int k2 = k*2;
523  // set contour with original edges of triangle
524  r0[k2+0] = triangles[e0].x(); z0[k2+0] = triangles[e0].y();
525  r0[k2+1] = triangles[e1].x(); z0[k2+1] = triangles[e1].y();
526  // set shifted radii
527  r1[k2+0] = r0[k2+0];
528  r1[k2+1] = r0[k2+1];
529  if (z0[k2+1] - z0[k2+0] <= 0) continue;
530  r1[k2+0] /= cosHalf;
531  r1[k2+1] /= cosHalf;
532  }
533 
534  // rotate countour, set sequence of 6-sided polygons
535  G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
536  G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
537  for (G4int j=0; j<6; ++j)
538  {
539  pols[0][j].set(r0[j]*cosStart,r0[j]*sinStart,z0[j]);
540  }
541  for (G4int k=1; k<ksteps+1; ++k)
542  {
543  for (G4int j=0; j<6; ++j)
544  {
545  pols[k][j].set(r1[j]*cosCur,r1[j]*sinCur,z0[j]);
546  }
547  G4double sinTmp = sinCur;
548  sinCur = sinCur*cosStep + cosCur*sinStep;
549  cosCur = cosCur*cosStep - sinTmp*sinStep;
550  }
551  for (G4int j=0; j<6; ++j)
552  {
553  pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j]*sinEnd,z0[j]);
554  }
555 
556  // set sub-envelope and adjust extent
558  G4BoundingEnvelope benv(polygons);
559  if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
560  if (emin < pMin) pMin = emin;
561  if (emax > pMax) pMax = emax;
562  if (eminlim > pMin && emaxlim < pMax) return true; // max possible extent
563  }
564  return (pMin < pMax);
565 }
566 
567 // GetEntityType
568 //
570 {
571  return G4String("G4GenericPolycone");
572 }
573 
574 // Clone
575 //
576 // Make a clone of the object
577 //
579 {
580  return new G4GenericPolycone(*this);
581 }
582 
583 // StreamInfo
584 //
585 // Stream object contents to an output stream
586 //
587 std::ostream& G4GenericPolycone::StreamInfo( std::ostream& os ) const
588 {
589  G4int oldprc = os.precision(16);
590  os << "-----------------------------------------------------------\n"
591  << " *** Dump for solid - " << GetName() << " ***\n"
592  << " ===================================================\n"
593  << " Solid type: G4GenericPolycone\n"
594  << " Parameters: \n"
595  << " starting phi angle : " << startPhi/degree << " degrees \n"
596  << " ending phi angle : " << endPhi/degree << " degrees \n";
597  G4int i=0;
598 
599  os << " number of RZ points: " << numCorner << "\n"
600  << " RZ values (corners): \n";
601  for (i=0; i<numCorner; i++)
602  {
603  os << " "
604  << corners[i].r << ", " << corners[i].z << "\n";
605  }
606  os << "-----------------------------------------------------------\n";
607  os.precision(oldprc);
608 
609  return os;
610 }
611 
612 // GetPointOnSurface
613 //
615 {
616  return GetPointOnSurfaceGeneric();
617 
618 }
619 
620 // CreatePolyhedron
621 //
623 {
624  // The following code prepares for:
625  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
626  // const double xyz[][3],
627  // const int faces_vec[][4])
628  // Here is an extract from the header file HepPolyhedron.h:
646  const G4int numSide =
648  * (endPhi - startPhi) / twopi) + 1;
649  G4int nNodes;
650  G4int nFaces;
651  typedef G4double double3[3];
652  double3* xyz;
653  typedef G4int int4[4];
654  int4* faces_vec;
655  if (phiIsOpen)
656  {
657  // Triangulate open ends. Simple ear-chopping algorithm...
658  // I'm not sure how robust this algorithm is (J.Allison).
659  //
660  std::vector<G4bool> chopped(numCorner, false);
661  std::vector<G4int*> triQuads;
662  G4int remaining = numCorner;
663  G4int iStarter = 0;
664  while (remaining >= 3) // Loop checking, 13.08.2015, G.Cosmo
665  {
666  // Find unchopped corners...
667  //
668  G4int A = -1, B = -1, C = -1;
669  G4int iStepper = iStarter;
670  do // Loop checking, 13.08.2015, G.Cosmo
671  {
672  if (A < 0) { A = iStepper; }
673  else if (B < 0) { B = iStepper; }
674  else if (C < 0) { C = iStepper; }
675  do // Loop checking, 13.08.2015, G.Cosmo
676  {
677  if (++iStepper >= numCorner) { iStepper = 0; }
678  }
679  while (chopped[iStepper]);
680  }
681  while (C < 0 && iStepper != iStarter);
682 
683  // Check triangle at B is pointing outward (an "ear").
684  // Sign of z cross product determines...
685  //
686  G4double BAr = corners[A].r - corners[B].r;
687  G4double BAz = corners[A].z - corners[B].z;
688  G4double BCr = corners[C].r - corners[B].r;
689  G4double BCz = corners[C].z - corners[B].z;
690  if (BAr * BCz - BAz * BCr < kCarTolerance)
691  {
692  G4int* tq = new G4int[3];
693  tq[0] = A + 1;
694  tq[1] = B + 1;
695  tq[2] = C + 1;
696  triQuads.push_back(tq);
697  chopped[B] = true;
698  --remaining;
699  }
700  else
701  {
702  do // Loop checking, 13.08.2015, G.Cosmo
703  {
704  if (++iStarter >= numCorner) { iStarter = 0; }
705  }
706  while (chopped[iStarter]);
707  }
708  }
709  // Transfer to faces...
710  //
711  nNodes = (numSide + 1) * numCorner;
712  nFaces = numSide * numCorner + 2 * triQuads.size();
713  faces_vec = new int4[nFaces];
714  G4int iface = 0;
715  G4int addition = numCorner * numSide;
716  G4int d = numCorner - 1;
717  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
718  {
719  for (size_t i = 0; i < triQuads.size(); ++i)
720  {
721  // Negative for soft/auxiliary/normally invisible edges...
722  //
723  G4int a, b, c;
724  if (iEnd == 0)
725  {
726  a = triQuads[i][0];
727  b = triQuads[i][1];
728  c = triQuads[i][2];
729  }
730  else
731  {
732  a = triQuads[i][0] + addition;
733  b = triQuads[i][2] + addition;
734  c = triQuads[i][1] + addition;
735  }
736  G4int ab = std::abs(b - a);
737  G4int bc = std::abs(c - b);
738  G4int ca = std::abs(a - c);
739  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
740  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
741  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
742  faces_vec[iface][3] = 0;
743  ++iface;
744  }
745  }
746 
747  // Continue with sides...
748 
749  xyz = new double3[nNodes];
750  const G4double dPhi = (endPhi - startPhi) / numSide;
752  G4int ixyz = 0;
753  for (G4int iSide = 0; iSide < numSide; ++iSide)
754  {
755  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
756  {
757  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
758  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
759  xyz[ixyz][2] = corners[iCorner].z;
760  if (iSide == 0) // startPhi
761  {
762  if (iCorner < numCorner - 1)
763  {
764  faces_vec[iface][0] = ixyz + 1;
765  faces_vec[iface][1] = -(ixyz + numCorner + 1);
766  faces_vec[iface][2] = ixyz + numCorner + 2;
767  faces_vec[iface][3] = ixyz + 2;
768  }
769  else
770  {
771  faces_vec[iface][0] = ixyz + 1;
772  faces_vec[iface][1] = -(ixyz + numCorner + 1);
773  faces_vec[iface][2] = ixyz + 2;
774  faces_vec[iface][3] = ixyz - numCorner + 2;
775  }
776  }
777  else if (iSide == numSide - 1) // endPhi
778  {
779  if (iCorner < numCorner - 1)
780  {
781  faces_vec[iface][0] = ixyz + 1;
782  faces_vec[iface][1] = ixyz + numCorner + 1;
783  faces_vec[iface][2] = ixyz + numCorner + 2;
784  faces_vec[iface][3] = -(ixyz + 2);
785  }
786  else
787  {
788  faces_vec[iface][0] = ixyz + 1;
789  faces_vec[iface][1] = ixyz + numCorner + 1;
790  faces_vec[iface][2] = ixyz + 2;
791  faces_vec[iface][3] = -(ixyz - numCorner + 2);
792  }
793  }
794  else
795  {
796  if (iCorner < numCorner - 1)
797  {
798  faces_vec[iface][0] = ixyz + 1;
799  faces_vec[iface][1] = -(ixyz + numCorner + 1);
800  faces_vec[iface][2] = ixyz + numCorner + 2;
801  faces_vec[iface][3] = -(ixyz + 2);
802  }
803  else
804  {
805  faces_vec[iface][0] = ixyz + 1;
806  faces_vec[iface][1] = -(ixyz + numCorner + 1);
807  faces_vec[iface][2] = ixyz + 2;
808  faces_vec[iface][3] = -(ixyz - numCorner + 2);
809  }
810  }
811  ++iface;
812  ++ixyz;
813  }
814  phi += dPhi;
815  }
816 
817  // Last corners...
818 
819  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
820  {
821  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
822  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
823  xyz[ixyz][2] = corners[iCorner].z;
824  ++ixyz;
825  }
826  }
827  else // !phiIsOpen - i.e., a complete 360 degrees.
828  {
829  nNodes = numSide * numCorner;
830  nFaces = numSide * numCorner;;
831  xyz = new double3[nNodes];
832  faces_vec = new int4[nFaces];
833  const G4double dPhi = (endPhi - startPhi) / numSide;
835  G4int ixyz = 0, iface = 0;
836  for (G4int iSide = 0; iSide < numSide; ++iSide)
837  {
838  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
839  {
840  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
841  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
842  xyz[ixyz][2] = corners[iCorner].z;
843 
844  if (iSide < numSide - 1)
845  {
846  if (iCorner < numCorner - 1)
847  {
848  faces_vec[iface][0] = ixyz + 1;
849  faces_vec[iface][1] = -(ixyz + numCorner + 1);
850  faces_vec[iface][2] = ixyz + numCorner + 2;
851  faces_vec[iface][3] = -(ixyz + 2);
852  }
853  else
854  {
855  faces_vec[iface][0] = ixyz + 1;
856  faces_vec[iface][1] = -(ixyz + numCorner + 1);
857  faces_vec[iface][2] = ixyz + 2;
858  faces_vec[iface][3] = -(ixyz - numCorner + 2);
859  }
860  }
861  else // Last side joins ends...
862  {
863  if (iCorner < numCorner - 1)
864  {
865  faces_vec[iface][0] = ixyz + 1;
866  faces_vec[iface][1] = -(ixyz + numCorner - nFaces + 1);
867  faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
868  faces_vec[iface][3] = -(ixyz + 2);
869  }
870  else
871  {
872  faces_vec[iface][0] = ixyz + 1;
873  faces_vec[iface][1] = -(ixyz - nFaces + numCorner + 1);
874  faces_vec[iface][2] = ixyz - nFaces + 2;
875  faces_vec[iface][3] = -(ixyz - numCorner + 2);
876  }
877  }
878  ++ixyz;
879  ++iface;
880  }
881  phi += dPhi;
882  }
883  }
884  G4Polyhedron* polyhedron = new G4Polyhedron;
885  G4int prob = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
886  delete [] faces_vec;
887  delete [] xyz;
888  if (prob)
889  {
890  std::ostringstream message;
891  message << "Problem creating G4Polyhedron for: " << GetName();
892  G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002",
893  JustWarning, message);
894  delete polyhedron;
895  return nullptr;
896  }
897  else
898  {
899  return polyhedron;
900  }
901 }
902 
903 #endif