Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Polyhedra.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: G4Polyhedra.cc 101819 2016-12-01 08:13:36Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4Polyhedra.cc
35 //
36 // Implementation of a CSG polyhedra, as an inherited class of G4VCSGfaceted.
37 //
38 // To be done:
39 // * Cracks: there are probably small cracks in the seams between the
40 // phi face (G4PolyPhiFace) and sides (G4PolyhedraSide) that are not
41 // entirely leakproof. Also, I am not sure all vertices are leak proof.
42 // * Many optimizations are possible, but not implemented.
43 // * Visualization needs to be updated outside of this routine.
44 //
45 // Utility classes:
46 // * G4EnclosingCylinder: I decided a quick check of geometry would be a
47 // good idea (for CPU speed). If the quick check fails, the regular
48 // full-blown G4VCSGfaceted version is invoked.
49 // * G4ReduciblePolygon: Really meant as a check of input parameters,
50 // this utility class also "converts" the GEANT3-like PGON/PCON
51 // arguments into the newer ones.
52 // Both these classes are implemented outside this file because they are
53 // shared with G4Polycone.
54 //
55 // --------------------------------------------------------------------
56 
57 #include "G4Polyhedra.hh"
58 
59 #if !defined(G4GEOM_USE_UPOLYHEDRA)
60 
61 #include "G4PolyhedraSide.hh"
62 #include "G4PolyPhiFace.hh"
63 
64 #include "G4GeomTools.hh"
65 #include "G4VoxelLimits.hh"
66 #include "G4AffineTransform.hh"
67 #include "G4BoundingEnvelope.hh"
68 
69 #include "Randomize.hh"
70 
71 #include "G4EnclosingCylinder.hh"
72 #include "G4ReduciblePolygon.hh"
73 #include "G4VPVParameterisation.hh"
74 
75 #include <sstream>
76 
77 using namespace CLHEP;
78 
79 //
80 // Constructor (GEANT3 style parameters)
81 //
82 // GEANT3 PGON radii are specified in the distance to the norm of each face.
83 //
85  G4double phiStart,
86  G4double thePhiTotal,
87  G4int theNumSide,
88  G4int numZPlanes,
89  const G4double zPlane[],
90  const G4double rInner[],
91  const G4double rOuter[] )
92  : G4VCSGfaceted( name ), genericPgon(false)
93 {
94  if (theNumSide <= 0)
95  {
96  std::ostringstream message;
97  message << "Solid must have at least one side - " << GetName() << G4endl
98  << " No sides specified !";
99  G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
100  FatalErrorInArgument, message);
101  }
102 
103  //
104  // Calculate conversion factor from G3 radius to G4 radius
105  //
106  G4double phiTotal = thePhiTotal;
107  if ( (phiTotal <=0) || (phiTotal >= twopi*(1-DBL_EPSILON)) )
108  { phiTotal = twopi; }
109  G4double convertRad = std::cos(0.5*phiTotal/theNumSide);
110 
111  //
112  // Some historical stuff
113  //
115 
116  original_parameters->numSide = theNumSide;
117  original_parameters->Start_angle = phiStart;
119  original_parameters->Num_z_planes = numZPlanes;
120  original_parameters->Z_values = new G4double[numZPlanes];
121  original_parameters->Rmin = new G4double[numZPlanes];
122  original_parameters->Rmax = new G4double[numZPlanes];
123 
124  G4int i;
125  for (i=0; i<numZPlanes; i++)
126  {
127  if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
128  {
129  if( (rInner[i] > rOuter[i+1])
130  ||(rInner[i+1] > rOuter[i]) )
131  {
132  DumpInfo();
133  std::ostringstream message;
134  message << "Cannot create a Polyhedra with no contiguous segments."
135  << G4endl
136  << " Segments are not contiguous !" << G4endl
137  << " rMin[" << i << "] = " << rInner[i]
138  << " -- rMax[" << i+1 << "] = " << rOuter[i+1] << G4endl
139  << " rMin[" << i+1 << "] = " << rInner[i+1]
140  << " -- rMax[" << i << "] = " << rOuter[i];
141  G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
142  FatalErrorInArgument, message);
143  }
144  }
145  original_parameters->Z_values[i] = zPlane[i];
146  original_parameters->Rmin[i] = rInner[i]/convertRad;
147  original_parameters->Rmax[i] = rOuter[i]/convertRad;
148  }
149 
150 
151  //
152  // Build RZ polygon using special PCON/PGON GEANT3 constructor
153  //
154  G4ReduciblePolygon *rz =
155  new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
156  rz->ScaleA( 1/convertRad );
157 
158  //
159  // Do the real work
160  //
161  Create( phiStart, phiTotal, theNumSide, rz );
162 
163  delete rz;
164 }
165 
166 
167 //
168 // Constructor (generic parameters)
169 //
171  G4double phiStart,
172  G4double phiTotal,
173  G4int theNumSide,
174  G4int numRZ,
175  const G4double r[],
176  const G4double z[] )
177  : G4VCSGfaceted( name ), genericPgon(true)
178 {
179  if (theNumSide <= 0)
180  {
181  std::ostringstream message;
182  message << "Solid must have at least one side - " << GetName() << G4endl
183  << " No sides specified !";
184  G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
185  FatalErrorInArgument, message);
186  }
187 
188  G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
189 
190  Create( phiStart, phiTotal, theNumSide, rz );
191 
192  // Set original_parameters struct for consistency
193  //
195 
196  delete rz;
197 }
198 
199 
200 //
201 // Create
202 //
203 // Generic create routine, called by each constructor
204 // after conversion of arguments
205 //
207  G4double phiTotal,
208  G4int theNumSide,
209  G4ReduciblePolygon *rz )
210 {
211  //
212  // Perform checks of rz values
213  //
214  if (rz->Amin() < 0.0)
215  {
216  std::ostringstream message;
217  message << "Illegal input parameters - " << GetName() << G4endl
218  << " All R values must be >= 0 !";
219  G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
220  FatalErrorInArgument, message);
221  }
222 
223  G4double rzArea = rz->Area();
224  if (rzArea < -kCarTolerance)
225  {
226  rz->ReverseOrder();
227  }
228  else if (rzArea < kCarTolerance)
229  {
230  std::ostringstream message;
231  message << "Illegal input parameters - " << GetName() << G4endl
232  << " R/Z cross section is zero or near zero: " << rzArea;
233  G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
234  FatalErrorInArgument, message);
235  }
236 
238  || (!rz->RemoveRedundantVertices( kCarTolerance )) )
239  {
240  std::ostringstream message;
241  message << "Illegal input parameters - " << GetName() << G4endl
242  << " Too few unique R/Z values !";
243  G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
244  FatalErrorInArgument, message);
245  }
246 
247  if (rz->CrossesItself( 1/kInfinity ))
248  {
249  std::ostringstream message;
250  message << "Illegal input parameters - " << GetName() << G4endl
251  << " R/Z segments cross !";
252  G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
253  FatalErrorInArgument, message);
254  }
255 
256  numCorner = rz->NumVertices();
257 
258 
259  startPhi = phiStart;
260  while( startPhi < 0 ) // Loop checking, 13.08.2015, G.Cosmo
261  startPhi += twopi;
262  //
263  // Phi opening? Account for some possible roundoff, and interpret
264  // nonsense value as representing no phi opening
265  //
266  if ( (phiTotal <= 0) || (phiTotal > twopi*(1-DBL_EPSILON)) )
267  {
268  phiIsOpen = false;
269  endPhi = phiStart+twopi;
270  }
271  else
272  {
273  phiIsOpen = true;
274 
275  //
276  // Convert phi into our convention
277  //
278  endPhi = phiStart+phiTotal;
279  while( endPhi < startPhi ) // Loop checking, 13.08.2015, G.Cosmo
280  endPhi += twopi;
281  }
282 
283  //
284  // Save number sides
285  //
286  numSide = theNumSide;
287 
288  //
289  // Allocate corner array.
290  //
292 
293  //
294  // Copy corners
295  //
296  G4ReduciblePolygonIterator iterRZ(rz);
297 
298  G4PolyhedraSideRZ *next = corners;
299  iterRZ.Begin();
300  do // Loop checking, 13.08.2015, G.Cosmo
301  {
302  next->r = iterRZ.GetA();
303  next->z = iterRZ.GetB();
304  } while( ++next, iterRZ.Next() );
305 
306  //
307  // Allocate face pointer array
308  //
310  faces = new G4VCSGface*[numFace];
311 
312  //
313  // Construct side faces
314  //
315  // To do so properly, we need to keep track of four successive RZ
316  // corners.
317  //
318  // But! Don't construct a face if both points are at zero radius!
319  //
320  G4PolyhedraSideRZ *corner = corners,
321  *prev = corners + numCorner-1,
322  *nextNext;
323  G4VCSGface **face = faces;
324  do // Loop checking, 13.08.2015, G.Cosmo
325  {
326  next = corner+1;
327  if (next >= corners+numCorner) next = corners;
328  nextNext = next+1;
329  if (nextNext >= corners+numCorner) nextNext = corners;
330 
331  if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
332 /*
333  // We must decide here if we can dare declare one of our faces
334  // as having a "valid" normal (i.e. allBehind = true). This
335  // is never possible if the face faces "inward" in r *unless*
336  // we have only one side
337  //
338  G4bool allBehind;
339  if ((corner->z > next->z) && (numSide > 1))
340  {
341  allBehind = false;
342  }
343  else
344  {
345  //
346  // Otherwise, it is only true if the line passing
347  // through the two points of the segment do not
348  // split the r/z cross section
349  //
350  allBehind = !rz->BisectedBy( corner->r, corner->z,
351  next->r, next->z, kCarTolerance );
352  }
353 */
354  *face++ = new G4PolyhedraSide( prev, corner, next, nextNext,
356  } while( prev=corner, corner=next, corner > corners );
357 
358  if (phiIsOpen)
359  {
360  //
361  // Construct phi open edges
362  //
363  *face++ = new G4PolyPhiFace( rz, startPhi, phiTotal/numSide, endPhi );
364  *face++ = new G4PolyPhiFace( rz, endPhi, phiTotal/numSide, startPhi );
365  }
366 
367  //
368  // We might have dropped a face or two: recalculate numFace
369  //
370  numFace = face-faces;
371 
372  //
373  // Make enclosingCylinder
374  //
376  new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
377 }
378 
379 
380 //
381 // Fake default constructor - sets only member data and allocates memory
382 // for usage restricted to object persistency.
383 //
385  : G4VCSGfaceted(a), numSide(0), startPhi(0.), endPhi(0.),
386  phiIsOpen(false), genericPgon(false), numCorner(0), corners(0),
387  original_parameters(0), enclosingCylinder(0)
388 {
389 }
390 
391 
392 //
393 // Destructor
394 //
396 {
397  delete [] corners;
399 
400  delete enclosingCylinder;
401 }
402 
403 
404 //
405 // Copy constructor
406 //
408  : G4VCSGfaceted( source )
409 {
410  CopyStuff( source );
411 }
412 
413 
414 //
415 // Assignment operator
416 //
418 {
419  if (this == &source) return *this;
420 
421  G4VCSGfaceted::operator=( source );
422 
423  delete [] corners;
425 
426  delete enclosingCylinder;
427 
428  CopyStuff( source );
429 
430  return *this;
431 }
432 
433 
434 //
435 // CopyStuff
436 //
437 void G4Polyhedra::CopyStuff( const G4Polyhedra &source )
438 {
439  //
440  // Simple stuff
441  //
442  numSide = source.numSide;
443  startPhi = source.startPhi;
444  endPhi = source.endPhi;
445  phiIsOpen = source.phiIsOpen;
446  numCorner = source.numCorner;
447  genericPgon= source.genericPgon;
448 
449  //
450  // The corner array
451  //
453 
454  G4PolyhedraSideRZ *corn = corners,
455  *sourceCorn = source.corners;
456  do // Loop checking, 13.08.2015, G.Cosmo
457  {
458  *corn = *sourceCorn;
459  } while( ++sourceCorn, ++corn < corners+numCorner );
460 
461  //
462  // Original parameters
463  //
464  if (source.original_parameters)
465  {
468  }
469 
470  //
471  // Enclosing cylinder
472  //
474 
475  fRebuildPolyhedron = false;
476  fpPolyhedron = 0;
477 }
478 
479 
480 //
481 // Reset
482 //
483 // Recalculates and reshapes the solid, given pre-assigned scaled
484 // original_parameters.
485 //
487 {
488  if (genericPgon)
489  {
490  std::ostringstream message;
491  message << "Solid " << GetName() << " built using generic construct."
492  << G4endl << "Not applicable to the generic construct !";
493  G4Exception("G4Polyhedra::Reset()", "GeomSolids1001",
494  JustWarning, message, "Parameters NOT resetted.");
495  return 1;
496  }
497 
498  //
499  // Clear old setup
500  //
502  delete [] corners;
503  delete enclosingCylinder;
504 
505  //
506  // Rebuild polyhedra
507  //
508  G4ReduciblePolygon *rz =
516  delete rz;
517 
518  return 0;
519 }
520 
521 
522 //
523 // Inside
524 //
525 // This is an override of G4VCSGfaceted::Inside, created in order
526 // to speed things up by first checking with G4EnclosingCylinder.
527 //
529 {
530  //
531  // Quick test
532  //
533  if (enclosingCylinder->MustBeOutside(p)) return kOutside;
534 
535  //
536  // Long answer
537  //
538  return G4VCSGfaceted::Inside(p);
539 }
540 
541 
542 //
543 // DistanceToIn
544 //
545 // This is an override of G4VCSGfaceted::Inside, created in order
546 // to speed things up by first checking with G4EnclosingCylinder.
547 //
549  const G4ThreeVector &v ) const
550 {
551  //
552  // Quick test
553  //
554  if (enclosingCylinder->ShouldMiss(p,v))
555  return kInfinity;
556 
557  //
558  // Long answer
559  //
560  return G4VCSGfaceted::DistanceToIn( p, v );
561 }
562 
563 
564 //
565 // DistanceToIn
566 //
568 {
569  return G4VCSGfaceted::DistanceToIn(p);
570 }
571 
573 //
574 // Get bounding box
575 
577 {
578  G4double rmin = kInfinity, rmax = -kInfinity;
579  G4double zmin = kInfinity, zmax = -kInfinity;
580  for (G4int i=0; i<GetNumRZCorner(); ++i)
581  {
582  G4PolyhedraSideRZ corner = GetCorner(i);
583  if (corner.r < rmin) rmin = corner.r;
584  if (corner.r > rmax) rmax = corner.r;
585  if (corner.z < zmin) zmin = corner.z;
586  if (corner.z > zmax) zmax = corner.z;
587  }
588 
589  G4double sphi = GetStartPhi();
590  G4double ephi = GetEndPhi();
591  G4double dphi = IsOpen() ? ephi-sphi : twopi;
592  G4int ksteps = GetNumSide();
593  G4double astep = dphi/ksteps;
594  G4double sinStep = std::sin(astep);
595  G4double cosStep = std::cos(astep);
596 
597  G4double sinCur = GetSinStartPhi();
598  G4double cosCur = GetCosStartPhi();
599  if (!IsOpen()) rmin = 0;
600  G4double xmin = rmin*cosCur, xmax = xmin;
601  G4double ymin = rmin*sinCur, ymax = ymin;
602  for (G4int k=0; k<ksteps+1; ++k)
603  {
604  G4double x = rmax*cosCur;
605  if (x < xmin) xmin = x;
606  if (x > xmax) xmax = x;
607  G4double y = rmax*sinCur;
608  if (y < ymin) ymin = y;
609  if (y > ymax) ymax = y;
610  if (rmin > 0)
611  {
612  G4double xx = rmin*cosCur;
613  if (xx < xmin) xmin = xx;
614  if (xx > xmax) xmax = xx;
615  G4double yy = rmin*sinCur;
616  if (yy < ymin) ymin = yy;
617  if (yy > ymax) ymax = yy;
618  }
619  G4double sinTmp = sinCur;
620  sinCur = sinCur*cosStep + cosCur*sinStep;
621  cosCur = cosCur*cosStep - sinTmp*sinStep;
622  }
623  pMin.set(xmin,ymin,zmin);
624  pMax.set(xmax,ymax,zmax);
625 
626  // Check correctness of the bounding box
627  //
628  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
629  {
630  std::ostringstream message;
631  message << "Bad bounding box (min >= max) for solid: "
632  << GetName() << " !"
633  << "\npMin = " << pMin
634  << "\npMax = " << pMax;
635  G4Exception("G4Polyhedra::Extent()", "GeomMgt0001", JustWarning, message);
636  DumpInfo();
637  }
638 }
639 
641 //
642 // Calculate extent under transform and specified limit
643 
645  const G4VoxelLimits& pVoxelLimit,
646  const G4AffineTransform& pTransform,
647  G4double& pMin, G4double& pMax) const
648 {
649  G4ThreeVector bmin, bmax;
650  G4bool exist;
651 
652  // Check bounding box (bbox)
653  //
654  Extent(bmin,bmax);
655  G4BoundingEnvelope bbox(bmin,bmax);
656 #ifdef G4BBOX_EXTENT
657  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
658 #endif
659  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
660  {
661  return exist = (pMin < pMax) ? true : false;
662  }
663 
664  // To find the extent, RZ contour of the polycone is subdivided
665  // in triangles. The extent is calculated as cumulative extent of
666  // all sub-polycones formed by rotation of triangles around Z
667  //
668  G4TwoVectorList contourRZ;
669  G4TwoVectorList triangles;
670  std::vector<G4int> iout;
671  G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
672  G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
673 
674  // get RZ contour, ensure anticlockwise order of corners
675  for (G4int i=0; i<GetNumRZCorner(); ++i)
676  {
677  G4PolyhedraSideRZ corner = GetCorner(i);
678  contourRZ.push_back(G4TwoVector(corner.r,corner.z));
679  }
681  G4double area = G4GeomTools::PolygonArea(contourRZ);
682  if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
683 
684  // triangulate RZ countour
685  if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
686  {
687  std::ostringstream message;
688  message << "Triangulation of RZ contour has failed for solid: "
689  << GetName() << " !"
690  << "\nExtent has been calculated using boundary box";
691  G4Exception("G4Polyhedra::CalculateExtent()",
692  "GeomMgt1002",JustWarning,message);
693  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
694  }
695 
696  // set trigonometric values
697  G4double sphi = GetStartPhi();
698  G4double ephi = GetEndPhi();
699  G4double dphi = IsOpen() ? ephi-sphi : twopi;
700  G4int ksteps = GetNumSide();
701  G4double astep = dphi/ksteps;
702  G4double sinStep = std::sin(astep);
703  G4double cosStep = std::cos(astep);
704  G4double sinStart = GetSinStartPhi();
705  G4double cosStart = GetCosStartPhi();
706 
707  // allocate vector lists
708  std::vector<const G4ThreeVectorList *> polygons;
709  polygons.resize(ksteps+1);
710  for (G4int k=0; k<ksteps+1; ++k) {
711  polygons[k] = new G4ThreeVectorList(3);
712  }
713 
714  // main loop along triangles
715  pMin = kInfinity;
716  pMax = -kInfinity;
717  G4int ntria = triangles.size()/3;
718  for (G4int i=0; i<ntria; ++i)
719  {
720  G4double sinCur = sinStart;
721  G4double cosCur = cosStart;
722  G4int i3 = i*3;
723  for (G4int k=0; k<ksteps+1; ++k) // rotate triangle
724  {
725  G4ThreeVectorList* ptr = const_cast<G4ThreeVectorList*>(polygons[k]);
726  G4ThreeVectorList::iterator iter = ptr->begin();
727  iter->set(triangles[i3+0].x()*cosCur,
728  triangles[i3+0].x()*sinCur,
729  triangles[i3+0].y());
730  iter++;
731  iter->set(triangles[i3+1].x()*cosCur,
732  triangles[i3+1].x()*sinCur,
733  triangles[i3+1].y());
734  iter++;
735  iter->set(triangles[i3+2].x()*cosCur,
736  triangles[i3+2].x()*sinCur,
737  triangles[i3+2].y());
738 
739  G4double sinTmp = sinCur;
740  sinCur = sinCur*cosStep + cosCur*sinStep;
741  cosCur = cosCur*cosStep - sinTmp*sinStep;
742  }
743 
744  // set sub-envelope and adjust extent
745  G4double emin,emax;
746  G4BoundingEnvelope benv(polygons);
747  if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
748  if (emin < pMin) pMin = emin;
749  if (emax > pMax) pMax = emax;
750  if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
751  }
752  // free memory
753  for (G4int k=0; k<ksteps+1; ++k) { delete polygons[k]; polygons[k]=0;}
754  return (pMin < pMax);
755 }
756 
757 //
758 // ComputeDimensions
759 //
761  const G4int n,
762  const G4VPhysicalVolume* pRep )
763 {
764  p->ComputeDimensions(*this,n,pRep);
765 }
766 
767 
768 //
769 // GetEntityType
770 //
772 {
773  return G4String("G4Polyhedra");
774 }
775 
776 
777 //
778 // Make a clone of the object
779 //
781 {
782  return new G4Polyhedra(*this);
783 }
784 
785 
786 //
787 // Stream object contents to an output stream
788 //
789 std::ostream& G4Polyhedra::StreamInfo( std::ostream& os ) const
790 {
791  G4int oldprc = os.precision(16);
792  os << "-----------------------------------------------------------\n"
793  << " *** Dump for solid - " << GetName() << " ***\n"
794  << " ===================================================\n"
795  << " Solid type: G4Polyhedra\n"
796  << " Parameters: \n"
797  << " starting phi angle : " << startPhi/degree << " degrees \n"
798  << " ending phi angle : " << endPhi/degree << " degrees \n"
799  << " number of sides : " << numSide << " \n";
800  G4int i=0;
801  if (!genericPgon)
802  {
804  os << " number of Z planes: " << numPlanes << "\n"
805  << " Z values: \n";
806  for (i=0; i<numPlanes; i++)
807  {
808  os << " Z plane " << i << ": "
809  << original_parameters->Z_values[i] << "\n";
810  }
811  os << " Tangent distances to inner surface (Rmin): \n";
812  for (i=0; i<numPlanes; i++)
813  {
814  os << " Z plane " << i << ": "
815  << original_parameters->Rmin[i] << "\n";
816  }
817  os << " Tangent distances to outer surface (Rmax): \n";
818  for (i=0; i<numPlanes; i++)
819  {
820  os << " Z plane " << i << ": "
821  << original_parameters->Rmax[i] << "\n";
822  }
823  }
824  os << " number of RZ points: " << numCorner << "\n"
825  << " RZ values (corners): \n";
826  for (i=0; i<numCorner; i++)
827  {
828  os << " "
829  << corners[i].r << ", " << corners[i].z << "\n";
830  }
831  os << "-----------------------------------------------------------\n";
832  os.precision(oldprc);
833 
834  return os;
835 }
836 
837 
838 //
839 // GetPointOnPlane
840 //
841 // Auxiliary method for get point on surface
842 //
845  G4ThreeVector p2, G4ThreeVector p3) const
846 {
847  G4double lambda1, lambda2, chose,aOne,aTwo;
848  G4ThreeVector t, u, v, w, Area, normal;
849  aOne = 1.;
850  aTwo = 1.;
851 
852  t = p1 - p0;
853  u = p2 - p1;
854  v = p3 - p2;
855  w = p0 - p3;
856 
857  chose = G4RandFlat::shoot(0.,aOne+aTwo);
858  if( (chose>=0.) && (chose < aOne) )
859  {
860  lambda1 = G4RandFlat::shoot(0.,1.);
861  lambda2 = G4RandFlat::shoot(0.,lambda1);
862  return (p2+lambda1*v+lambda2*w);
863  }
864 
865  lambda1 = G4RandFlat::shoot(0.,1.);
866  lambda2 = G4RandFlat::shoot(0.,lambda1);
867  return (p0+lambda1*t+lambda2*u);
868 }
869 
870 
871 //
872 // GetPointOnTriangle
873 //
874 // Auxiliary method for get point on surface
875 //
877  G4ThreeVector p2,
878  G4ThreeVector p3) const
879 {
880  G4double lambda1,lambda2;
881  G4ThreeVector v=p3-p1, w=p1-p2;
882 
883  lambda1 = G4RandFlat::shoot(0.,1.);
884  lambda2 = G4RandFlat::shoot(0.,lambda1);
885 
886  return (p2 + lambda1*w + lambda2*v);
887 }
888 
889 
890 //
891 // GetPointOnSurface
892 //
894 {
895  if( !genericPgon ) // Polyhedra by faces
896  {
897  G4int j, numPlanes = original_parameters->Num_z_planes, Flag=0;
898  G4double chose, totArea=0., Achose1, Achose2,
899  rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
900  G4double a, b, l2, rang, totalPhi, ksi,
901  area, aTop=0., aBottom=0., zVal=0.;
902 
903  G4ThreeVector p0, p1, p2, p3;
904  std::vector<G4double> aVector1;
905  std::vector<G4double> aVector2;
906  std::vector<G4double> aVector3;
907 
908  totalPhi= (phiIsOpen) ? (endPhi-startPhi) : twopi;
909  ksi = totalPhi/numSide;
910  G4double cosksi = std::cos(ksi/2.);
911 
912  // Below we generate the areas relevant to our solid
913  //
914  for(j=0; j<numPlanes-1; j++)
915  {
916  a = original_parameters->Rmax[j+1];
917  b = original_parameters->Rmax[j];
919  -original_parameters->Z_values[j+1]) + sqr(b-a);
920  area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
921  aVector1.push_back(area);
922  }
923 
924  for(j=0; j<numPlanes-1; j++)
925  {
926  a = original_parameters->Rmin[j+1];//*cosksi;
927  b = original_parameters->Rmin[j];//*cosksi;
929  -original_parameters->Z_values[j+1]) + sqr(b-a);
930  area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
931  aVector2.push_back(area);
932  }
933 
934  for(j=0; j<numPlanes-1; j++)
935  {
936  if(phiIsOpen == true)
937  {
938  aVector3.push_back(0.5*(original_parameters->Rmax[j]
941  -original_parameters->Rmin[j+1])
942  *std::fabs(original_parameters->Z_values[j+1]
944  }
945  else { aVector3.push_back(0.); }
946  }
947 
948  for(j=0; j<numPlanes-1; j++)
949  {
950  totArea += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
951  }
952 
953  // Must include top and bottom areas
954  //
955  if(original_parameters->Rmax[numPlanes-1] != 0.)
956  {
957  a = original_parameters->Rmax[numPlanes-1];
958  b = original_parameters->Rmin[numPlanes-1];
959  l2 = sqr(a-b);
960  aTop = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
961  }
962 
963  if(original_parameters->Rmax[0] != 0.)
964  {
965  a = original_parameters->Rmax[0];
966  b = original_parameters->Rmin[0];
967  l2 = sqr(a-b);
968  aBottom = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
969  }
970 
971  Achose1 = 0.;
972  Achose2 = numSide*(aVector1[0]+aVector2[0])+2.*aVector3[0];
973 
974  chose = G4RandFlat::shoot(0.,totArea+aTop+aBottom);
975  if( (chose >= 0.) && (chose < aTop + aBottom) )
976  {
977  chose = G4RandFlat::shoot(startPhi,startPhi+totalPhi);
978  rang = std::floor((chose-startPhi)/ksi-0.01);
979  if(rang<0) { rang=0; }
980  rang = std::fabs(rang);
981  sinphi1 = std::sin(startPhi+rang*ksi);
982  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
983  cosphi1 = std::cos(startPhi+rang*ksi);
984  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
985  chose = G4RandFlat::shoot(0., aTop + aBottom);
986  if(chose>=0. && chose<aTop)
987  {
988  rad1 = original_parameters->Rmin[numPlanes-1];
989  rad2 = original_parameters->Rmax[numPlanes-1];
990  zVal = original_parameters->Z_values[numPlanes-1];
991  }
992  else
993  {
994  rad1 = original_parameters->Rmin[0];
995  rad2 = original_parameters->Rmax[0];
996  zVal = original_parameters->Z_values[0];
997  }
998  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
999  p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
1000  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
1001  p3 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
1002  return GetPointOnPlane(p0,p1,p2,p3);
1003  }
1004  else
1005  {
1006  for (j=0; j<numPlanes-1; j++)
1007  {
1008  if( ((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes-2) )
1009  {
1010  Flag = j; break;
1011  }
1012  Achose1 += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
1013  Achose2 = Achose1 + numSide*(aVector1[j+1]+aVector2[j+1])
1014  + 2.*aVector3[j+1];
1015  }
1016  }
1017 
1018  // At this point we have chosen a subsection
1019  // between to adjacent plane cuts...
1020 
1021  j = Flag;
1022 
1023  totArea = numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
1024  chose = G4RandFlat::shoot(0.,totArea);
1025 
1026  if( (chose>=0.) && (chose<numSide*aVector1[j]) )
1027  {
1028  chose = G4RandFlat::shoot(startPhi,startPhi+totalPhi);
1029  rang = std::floor((chose-startPhi)/ksi-0.01);
1030  if(rang<0) { rang=0; }
1031  rang = std::fabs(rang);
1032  rad1 = original_parameters->Rmax[j];
1033  rad2 = original_parameters->Rmax[j+1];
1034  sinphi1 = std::sin(startPhi+rang*ksi);
1035  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
1036  cosphi1 = std::cos(startPhi+rang*ksi);
1037  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
1038  zVal = original_parameters->Z_values[j];
1039 
1040  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
1041  p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
1042 
1043  zVal = original_parameters->Z_values[j+1];
1044 
1045  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
1046  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
1047  return GetPointOnPlane(p0,p1,p2,p3);
1048  }
1049  else if ( (chose >= numSide*aVector1[j])
1050  && (chose <= numSide*(aVector1[j]+aVector2[j])) )
1051  {
1052  chose = G4RandFlat::shoot(startPhi,startPhi+totalPhi);
1053  rang = std::floor((chose-startPhi)/ksi-0.01);
1054  if(rang<0) { rang=0; }
1055  rang = std::fabs(rang);
1056  rad1 = original_parameters->Rmin[j];
1057  rad2 = original_parameters->Rmin[j+1];
1058  sinphi1 = std::sin(startPhi+rang*ksi);
1059  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
1060  cosphi1 = std::cos(startPhi+rang*ksi);
1061  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
1062  zVal = original_parameters->Z_values[j];
1063 
1064  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
1065  p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
1066 
1067  zVal = original_parameters->Z_values[j+1];
1068 
1069  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
1070  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
1071  return GetPointOnPlane(p0,p1,p2,p3);
1072  }
1073 
1074  chose = G4RandFlat::shoot(0.,2.2);
1075  if( (chose>=0.) && (chose < 1.) )
1076  {
1077  rang = startPhi;
1078  }
1079  else
1080  {
1081  rang = endPhi;
1082  }
1083 
1084  cosphi1 = std::cos(rang); rad1 = original_parameters->Rmin[j];
1085  sinphi1 = std::sin(rang); rad2 = original_parameters->Rmax[j];
1086 
1087  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
1089  p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
1091 
1092  rad1 = original_parameters->Rmax[j+1];
1093  rad2 = original_parameters->Rmin[j+1];
1094 
1095  p2 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
1097  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
1099  return GetPointOnPlane(p0,p1,p2,p3);
1100  }
1101  else // Generic polyhedra
1102  {
1103  return GetPointOnSurfaceGeneric();
1104  }
1105 }
1106 
1107 //
1108 // CreatePolyhedron
1109 //
1111 {
1112  if (!genericPgon)
1113  {
1121  }
1122  else
1123  {
1124  // The following code prepares for:
1125  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
1126  // const double xyz[][3],
1127  // const int faces_vec[][4])
1128  // Here is an extract from the header file HepPolyhedron.h:
1146  G4int nNodes;
1147  G4int nFaces;
1148  typedef G4double double3[3];
1149  double3* xyz;
1150  typedef G4int int4[4];
1151  int4* faces_vec;
1152  if (phiIsOpen)
1153  {
1154  // Triangulate open ends. Simple ear-chopping algorithm...
1155  // I'm not sure how robust this algorithm is (J.Allison).
1156  //
1157  std::vector<G4bool> chopped(numCorner, false);
1158  std::vector<G4int*> triQuads;
1159  G4int remaining = numCorner;
1160  G4int iStarter = 0;
1161  while (remaining >= 3) // Loop checking, 13.08.2015, G.Cosmo
1162  {
1163  // Find unchopped corners...
1164  //
1165  G4int A = -1, B = -1, C = -1;
1166  G4int iStepper = iStarter;
1167  do // Loop checking, 13.08.2015, G.Cosmo
1168  {
1169  if (A < 0) { A = iStepper; }
1170  else if (B < 0) { B = iStepper; }
1171  else if (C < 0) { C = iStepper; }
1172  do // Loop checking, 13.08.2015, G.Cosmo
1173  {
1174  if (++iStepper >= numCorner) iStepper = 0;
1175  }
1176  while (chopped[iStepper]);
1177  }
1178  while (C < 0 && iStepper != iStarter);
1179 
1180  // Check triangle at B is pointing outward (an "ear").
1181  // Sign of z cross product determines...
1182 
1183  G4double BAr = corners[A].r - corners[B].r;
1184  G4double BAz = corners[A].z - corners[B].z;
1185  G4double BCr = corners[C].r - corners[B].r;
1186  G4double BCz = corners[C].z - corners[B].z;
1187  if (BAr * BCz - BAz * BCr < kCarTolerance)
1188  {
1189  G4int* tq = new G4int[3];
1190  tq[0] = A + 1;
1191  tq[1] = B + 1;
1192  tq[2] = C + 1;
1193  triQuads.push_back(tq);
1194  chopped[B] = true;
1195  --remaining;
1196  }
1197  else
1198  {
1199  do // Loop checking, 13.08.2015, G.Cosmo
1200  {
1201  if (++iStarter >= numCorner) { iStarter = 0; }
1202  }
1203  while (chopped[iStarter]);
1204  }
1205  }
1206 
1207  // Transfer to faces...
1208 
1209  nNodes = (numSide + 1) * numCorner;
1210  nFaces = numSide * numCorner + 2 * triQuads.size();
1211  faces_vec = new int4[nFaces];
1212  G4int iface = 0;
1213  G4int addition = numCorner * numSide;
1214  G4int d = numCorner - 1;
1215  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
1216  {
1217  for (size_t i = 0; i < triQuads.size(); ++i)
1218  {
1219  // Negative for soft/auxiliary/normally invisible edges...
1220  //
1221  G4int a, b, c;
1222  if (iEnd == 0)
1223  {
1224  a = triQuads[i][0];
1225  b = triQuads[i][1];
1226  c = triQuads[i][2];
1227  }
1228  else
1229  {
1230  a = triQuads[i][0] + addition;
1231  b = triQuads[i][2] + addition;
1232  c = triQuads[i][1] + addition;
1233  }
1234  G4int ab = std::abs(b - a);
1235  G4int bc = std::abs(c - b);
1236  G4int ca = std::abs(a - c);
1237  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
1238  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
1239  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
1240  faces_vec[iface][3] = 0;
1241  ++iface;
1242  }
1243  }
1244 
1245  // Continue with sides...
1246 
1247  xyz = new double3[nNodes];
1248  const G4double dPhi = (endPhi - startPhi) / numSide;
1249  G4double phi = startPhi;
1250  G4int ixyz = 0;
1251  for (G4int iSide = 0; iSide < numSide; ++iSide)
1252  {
1253  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1254  {
1255  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1256  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1257  xyz[ixyz][2] = corners[iCorner].z;
1258  if (iCorner < numCorner - 1)
1259  {
1260  faces_vec[iface][0] = ixyz + 1;
1261  faces_vec[iface][1] = ixyz + numCorner + 1;
1262  faces_vec[iface][2] = ixyz + numCorner + 2;
1263  faces_vec[iface][3] = ixyz + 2;
1264  }
1265  else
1266  {
1267  faces_vec[iface][0] = ixyz + 1;
1268  faces_vec[iface][1] = ixyz + numCorner + 1;
1269  faces_vec[iface][2] = ixyz + 2;
1270  faces_vec[iface][3] = ixyz - numCorner + 2;
1271  }
1272  ++iface;
1273  ++ixyz;
1274  }
1275  phi += dPhi;
1276  }
1277 
1278  // Last corners...
1279 
1280  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1281  {
1282  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1283  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1284  xyz[ixyz][2] = corners[iCorner].z;
1285  ++ixyz;
1286  }
1287  }
1288  else // !phiIsOpen - i.e., a complete 360 degrees.
1289  {
1290  nNodes = numSide * numCorner;
1291  nFaces = numSide * numCorner;;
1292  xyz = new double3[nNodes];
1293  faces_vec = new int4[nFaces];
1294  // const G4double dPhi = (endPhi - startPhi) / numSide;
1295  const G4double dPhi = twopi / numSide;
1296  // !phiIsOpen endPhi-startPhi = 360 degrees.
1297  G4double phi = startPhi;
1298  G4int ixyz = 0, iface = 0;
1299  for (G4int iSide = 0; iSide < numSide; ++iSide)
1300  {
1301  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1302  {
1303  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1304  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1305  xyz[ixyz][2] = corners[iCorner].z;
1306  if (iSide < numSide - 1)
1307  {
1308  if (iCorner < numCorner - 1)
1309  {
1310  faces_vec[iface][0] = ixyz + 1;
1311  faces_vec[iface][1] = ixyz + numCorner + 1;
1312  faces_vec[iface][2] = ixyz + numCorner + 2;
1313  faces_vec[iface][3] = ixyz + 2;
1314  }
1315  else
1316  {
1317  faces_vec[iface][0] = ixyz + 1;
1318  faces_vec[iface][1] = ixyz + numCorner + 1;
1319  faces_vec[iface][2] = ixyz + 2;
1320  faces_vec[iface][3] = ixyz - numCorner + 2;
1321  }
1322  }
1323  else // Last side joins ends...
1324  {
1325  if (iCorner < numCorner - 1)
1326  {
1327  faces_vec[iface][0] = ixyz + 1;
1328  faces_vec[iface][1] = ixyz + numCorner - nFaces + 1;
1329  faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
1330  faces_vec[iface][3] = ixyz + 2;
1331  }
1332  else
1333  {
1334  faces_vec[iface][0] = ixyz + 1;
1335  faces_vec[iface][1] = ixyz - nFaces + numCorner + 1;
1336  faces_vec[iface][2] = ixyz - nFaces + 2;
1337  faces_vec[iface][3] = ixyz - numCorner + 2;
1338  }
1339  }
1340  ++ixyz;
1341  ++iface;
1342  }
1343  phi += dPhi;
1344  }
1345  }
1346  G4Polyhedron* polyhedron = new G4Polyhedron;
1347  G4int problem = polyhedron->createPolyhedron(nNodes,nFaces,xyz,faces_vec);
1348  delete [] faces_vec;
1349  delete [] xyz;
1350  if (problem)
1351  {
1352  std::ostringstream message;
1353  message << "Problem creating G4Polyhedron for: " << GetName();
1354  G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
1355  JustWarning, message);
1356  delete polyhedron;
1357  return 0;
1358  }
1359  else
1360  {
1361  return polyhedron;
1362  }
1363  }
1364 }
1365 
1366 
1368 {
1369  G4int numPlanes = (G4int)numCorner;
1370  G4bool isConvertible=true;
1371  G4double Zmax=rz->Bmax();
1372  rz->StartWithZMin();
1373 
1374  // Prepare vectors for storage
1375  //
1376  std::vector<G4double> Z;
1377  std::vector<G4double> Rmin;
1378  std::vector<G4double> Rmax;
1379 
1380  G4int countPlanes=1;
1381  G4int icurr=0;
1382  G4int icurl=0;
1383 
1384  // first plane Z=Z[0]
1385  //
1386  Z.push_back(corners[0].z);
1387  G4double Zprev=Z[0];
1388  if (Zprev == corners[1].z)
1389  {
1390  Rmin.push_back(corners[0].r);
1391  Rmax.push_back (corners[1].r);icurr=1;
1392  }
1393  else if (Zprev == corners[numPlanes-1].z)
1394  {
1395  Rmin.push_back(corners[numPlanes-1].r);
1396  Rmax.push_back (corners[0].r);
1397  icurl=numPlanes-1;
1398  }
1399  else
1400  {
1401  Rmin.push_back(corners[0].r);
1402  Rmax.push_back (corners[0].r);
1403  }
1404 
1405  // next planes until last
1406  //
1407  G4int inextr=0, inextl=0;
1408  for (G4int i=0; i < numPlanes-2; i++)
1409  {
1410  inextr=1+icurr;
1411  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1412 
1413  if((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax)) { break; }
1414 
1415  G4double Zleft = corners[inextl].z;
1416  G4double Zright = corners[inextr].z;
1417  if(Zright>Zleft)
1418  {
1419  Z.push_back(Zleft);
1420  countPlanes++;
1421  G4double difZr=corners[inextr].z - corners[icurr].z;
1422  G4double difZl=corners[inextl].z - corners[icurl].z;
1423 
1424  if(std::fabs(difZl) < kCarTolerance)
1425  {
1426  if(std::fabs(difZr) < kCarTolerance)
1427  {
1428  Rmin.push_back(corners[inextl].r);
1429  Rmax.push_back(corners[icurr].r);
1430  }
1431  else
1432  {
1433  Rmin.push_back(corners[inextl].r);
1434  Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1435  *(corners[inextr].r - corners[icurr].r));
1436  }
1437  }
1438  else if (difZl >= kCarTolerance)
1439  {
1440  if(std::fabs(difZr) < kCarTolerance)
1441  {
1442  Rmin.push_back(corners[icurl].r);
1443  Rmax.push_back(corners[icurr].r);
1444  }
1445  else
1446  {
1447  Rmin.push_back(corners[icurl].r);
1448  Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1449  *(corners[inextr].r - corners[icurr].r));
1450  }
1451  }
1452  else
1453  {
1454  isConvertible=false; break;
1455  }
1456  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1457  }
1458  else if(std::fabs(Zright-Zleft)<kCarTolerance) // Zright=Zleft
1459  {
1460  Z.push_back(Zleft);
1461  countPlanes++;
1462  icurr++;
1463 
1464  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1465 
1466  Rmin.push_back(corners[inextl].r);
1467  Rmax.push_back (corners[inextr].r);
1468  }
1469  else // Zright<Zleft
1470  {
1471  Z.push_back(Zright);
1472  countPlanes++;
1473 
1474  G4double difZr=corners[inextr].z - corners[icurr].z;
1475  G4double difZl=corners[inextl].z - corners[icurl].z;
1476  if(std::fabs(difZr) < kCarTolerance)
1477  {
1478  if(std::fabs(difZl) < kCarTolerance)
1479  {
1480  Rmax.push_back(corners[inextr].r);
1481  Rmin.push_back(corners[icurr].r);
1482  }
1483  else
1484  {
1485  Rmin.push_back(corners[icurl].r + (Zright-corners[icurl].z)/difZl
1486  * (corners[inextl].r - corners[icurl].r));
1487  Rmax.push_back(corners[inextr].r);
1488  }
1489  icurr++;
1490  } // plate
1491  else if (difZr >= kCarTolerance)
1492  {
1493  if(std::fabs(difZl) < kCarTolerance)
1494  {
1495  Rmax.push_back(corners[inextr].r);
1496  Rmin.push_back (corners[icurr].r);
1497  }
1498  else
1499  {
1500  Rmax.push_back(corners[inextr].r);
1501  Rmin.push_back (corners[icurl].r+(Zright-corners[icurl].z)/difZl
1502  * (corners[inextl].r - corners[icurl].r));
1503  }
1504  icurr++;
1505  }
1506  else
1507  {
1508  isConvertible=false; break;
1509  }
1510  }
1511  } // end for loop
1512 
1513  // last plane Z=Zmax
1514  //
1515  Z.push_back(Zmax);
1516  countPlanes++;
1517  inextr=1+icurr;
1518  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1519 
1520  Rmax.push_back(corners[inextr].r);
1521  Rmin.push_back(corners[inextl].r);
1522 
1523  // Set original parameters Rmin,Rmax,Z
1524  //
1525  if(isConvertible)
1526  {
1529  original_parameters->Z_values = new G4double[countPlanes];
1530  original_parameters->Rmin = new G4double[countPlanes];
1531  original_parameters->Rmax = new G4double[countPlanes];
1532 
1533  for(G4int j=0; j < countPlanes; j++)
1534  {
1535  original_parameters->Z_values[j] = Z[j];
1536  original_parameters->Rmax[j] = Rmax[j];
1537  original_parameters->Rmin[j] = Rmin[j];
1538  }
1541  original_parameters->Num_z_planes = countPlanes;
1542 
1543  }
1544  else // Set parameters(r,z) with Rmin==0 as convention
1545  {
1546 #ifdef G4SPECSDEBUG
1547  std::ostringstream message;
1548  message << "Polyhedra " << GetName() << G4endl
1549  << "cannot be converted to Polyhedra with (Rmin,Rmaz,Z) parameters!";
1550  G4Exception("G4Polyhedra::SetOriginalParameters()",
1551  "GeomSolids0002", JustWarning, message);
1552 #endif
1555  original_parameters->Z_values = new G4double[numPlanes];
1556  original_parameters->Rmin = new G4double[numPlanes];
1557  original_parameters->Rmax = new G4double[numPlanes];
1558 
1559  for(G4int j=0; j < numPlanes; j++)
1560  {
1562  original_parameters->Rmax[j] = corners[j].r;
1563  original_parameters->Rmin[j] = 0.0;
1564  }
1567  original_parameters->Num_z_planes = numPlanes;
1568  }
1569  //return isConvertible;
1570 }
1571 
1572 #endif
void set(double x, double y, double z)
void CopyStuff(const G4Polyhedra &source)
Definition: G4Polyhedra.cc:437
G4String GetName() const
G4bool CrossesItself(G4double tolerance)
const XML_Char * name
Definition: expat.h:151
ThreeVector shoot(const G4int Ap, const G4int Af)
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4Polyhedra.cc:576
G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
virtual ~G4Polyhedra()
Definition: G4Polyhedra.cc:395
G4bool fRebuildPolyhedron
G4Polyhedron * fpPolyhedron
static const G4double kInfinity
Definition: geomdefs.hh:42
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Polyhedra.cc:644
CLHEP::Hep3Vector G4ThreeVector
G4double GetCosStartPhi() const
G4GeometryType GetEntityType() const
Definition: G4Polyhedra.cc:771
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])
G4PolyhedraHistorical * original_parameters
Definition: G4Polyhedra.hh:197
std::vector< G4TwoVector > G4TwoVectorList
Definition: G4GeomTools.hh:50
const char * p
Definition: xmltok.h:285
G4ThreeVector GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3) const
Definition: G4Polyhedra.cc:844
G4bool genericPgon
Definition: G4Polyhedra.hh:194
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
static G4double PolygonArea(const G4TwoVectorList &polygon)
Definition: G4GeomTools.cc:82
G4double GetSinStartPhi() const
G4PolyhedraSideRZ * corners
Definition: G4Polyhedra.hh:196
tuple x
Definition: test.py:50
double C(double temp)
G4bool MustBeOutside(const G4ThreeVector &p) const
double B(double temperature)
G4bool Reset()
Definition: G4Polyhedra.cc:486
G4double GetEndPhi() const
int G4int
Definition: G4Types.hh:78
void Create(G4double phiStart, G4double phiTotal, G4int numSide, G4ReduciblePolygon *rz)
Definition: G4Polyhedra.cc:206
G4ThreeVector GetPointOnSurface() const
Definition: G4Polyhedra.cc:893
double z() const
void DumpInfo() const
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double endPhi
Definition: G4Polyhedra.hh:192
G4int GetNumRZCorner() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4bool RemoveDuplicateVertices(G4double tolerance)
G4EnclosingCylinder * enclosingCylinder
Definition: G4Polyhedra.hh:199
tuple b
Definition: test.py:12
G4bool RemoveRedundantVertices(G4double tolerance)
G4int GetNumSide() const
G4ThreeVector GetPointOnSurfaceGeneric() const
double A(double temperature)
G4int numCorner
Definition: G4Polyhedra.hh:195
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Polyhedra.cc:760
G4int NumVertices() const
static constexpr double degree
Definition: G4SIunits.hh:144
G4double Bmax() const
void ScaleA(G4double scale)
bool G4bool
Definition: G4Types.hh:79
G4double startPhi
Definition: G4Polyhedra.hh:191
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
#define DBL_EPSILON
Definition: templates.hh:87
const G4int n
static G4bool TriangulatePolygon(const G4TwoVectorList &polygon, G4TwoVectorList &result)
Definition: G4GeomTools.cc:178
std::vector< G4ThreeVector > G4ThreeVectorList
tuple v
Definition: test.py:18
G4bool IsOpen() const
G4VCSGface ** faces
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4double emax
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Polyhedra.cc:548
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
EInside Inside(const G4ThreeVector &p) const
Definition: G4Polyhedra.cc:528
static const G4double ab
G4Polyhedron * CreatePolyhedron() const
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
double y() const
tuple z
Definition: test.py:28
G4Polyhedra & operator=(const G4Polyhedra &source)
Definition: G4Polyhedra.cc:417
CLHEP::Hep2Vector G4TwoVector
Definition: G4TwoVector.hh:42
G4double GetStartPhi() const
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4PolyhedraSideRZ GetCorner(const G4int index) const
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
void SetOriginalParameters(G4PolyhedraHistorical *pars)
tuple c
Definition: test.py:13
G4double GetMaxExtent(const EAxis pAxis) const
G4ThreeVector GetPointOnTriangle(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2) const
Definition: G4Polyhedra.cc:876
virtual EInside Inside(const G4ThreeVector &p) const
G4bool phiIsOpen
Definition: G4Polyhedra.hh:193
G4bool ShouldMiss(const G4ThreeVector &p, const G4ThreeVector &v) const
tuple astep
Definition: test1.py:13
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Polyhedra.cc:789
G4double GetMinExtent(const EAxis pAxis) const
G4VSolid * Clone() const
Definition: G4Polyhedra.cc:780
G4Polyhedra(const G4String &name, G4double phiStart, G4double phiTotal, G4int numSide, G4int numZPlanes, const G4double zPlane[], const G4double rInner[], const G4double rOuter[])
Definition: G4Polyhedra.cc:84
static void RemoveRedundantVertices(G4TwoVectorList &polygon, std::vector< G4int > &iout, G4double tolerance=0)
Definition: G4GeomTools.cc:293