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