Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4Polyhedra.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 //
27 // $Id: G4Polyhedra.cc 104316 2017-05-24 13:04:23Z 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  G4ThreeVector& pMax) const
578 {
579  G4double rmin = kInfinity, rmax = -kInfinity;
580  G4double zmin = kInfinity, zmax = -kInfinity;
581  for (G4int i=0; i<GetNumRZCorner(); ++i)
582  {
583  G4PolyhedraSideRZ corner = GetCorner(i);
584  if (corner.r < rmin) rmin = corner.r;
585  if (corner.r > rmax) rmax = corner.r;
586  if (corner.z < zmin) zmin = corner.z;
587  if (corner.z > zmax) zmax = corner.z;
588  }
589 
590  G4double sphi = GetStartPhi();
591  G4double ephi = GetEndPhi();
592  G4double dphi = IsOpen() ? ephi-sphi : twopi;
593  G4int ksteps = GetNumSide();
594  G4double astep = dphi/ksteps;
595  G4double sinStep = std::sin(astep);
596  G4double cosStep = std::cos(astep);
597 
598  G4double sinCur = GetSinStartPhi();
599  G4double cosCur = GetCosStartPhi();
600  if (!IsOpen()) rmin = 0;
601  G4double xmin = rmin*cosCur, xmax = xmin;
602  G4double ymin = rmin*sinCur, ymax = ymin;
603  for (G4int k=0; k<ksteps+1; ++k)
604  {
605  G4double x = rmax*cosCur;
606  if (x < xmin) xmin = x;
607  if (x > xmax) xmax = x;
608  G4double y = rmax*sinCur;
609  if (y < ymin) ymin = y;
610  if (y > ymax) ymax = y;
611  if (rmin > 0)
612  {
613  G4double xx = rmin*cosCur;
614  if (xx < xmin) xmin = xx;
615  if (xx > xmax) xmax = xx;
616  G4double yy = rmin*sinCur;
617  if (yy < ymin) ymin = yy;
618  if (yy > ymax) ymax = yy;
619  }
620  G4double sinTmp = sinCur;
621  sinCur = sinCur*cosStep + cosCur*sinStep;
622  cosCur = cosCur*cosStep - sinTmp*sinStep;
623  }
624  pMin.set(xmin,ymin,zmin);
625  pMax.set(xmax,ymax,zmax);
626 
627  // Check correctness of the bounding box
628  //
629  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
630  {
631  std::ostringstream message;
632  message << "Bad bounding box (min >= max) for solid: "
633  << GetName() << " !"
634  << "\npMin = " << pMin
635  << "\npMax = " << pMax;
636  G4Exception("G4Polyhedra::BoundingLimits()", "GeomMgt0001",
637  JustWarning, message);
638  DumpInfo();
639  }
640 }
641 
643 //
644 // Calculate extent under transform and specified limit
645 
647  const G4VoxelLimits& pVoxelLimit,
648  const G4AffineTransform& pTransform,
649  G4double& pMin, G4double& pMax) const
650 {
651  G4ThreeVector bmin, bmax;
652  G4bool exist;
653 
654  // Check bounding box (bbox)
655  //
656  BoundingLimits(bmin,bmax);
657  G4BoundingEnvelope bbox(bmin,bmax);
658 #ifdef G4BBOX_EXTENT
659  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
660 #endif
661  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
662  {
663  return exist = (pMin < pMax) ? true : false;
664  }
665 
666  // To find the extent, RZ contour of the polycone is subdivided
667  // in triangles. The extent is calculated as cumulative extent of
668  // all sub-polycones formed by rotation of triangles around Z
669  //
670  G4TwoVectorList contourRZ;
671  G4TwoVectorList triangles;
672  std::vector<G4int> iout;
673  G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
674  G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
675 
676  // get RZ contour, ensure anticlockwise order of corners
677  for (G4int i=0; i<GetNumRZCorner(); ++i)
678  {
679  G4PolyhedraSideRZ corner = GetCorner(i);
680  contourRZ.push_back(G4TwoVector(corner.r,corner.z));
681  }
683  G4double area = G4GeomTools::PolygonArea(contourRZ);
684  if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
685 
686  // triangulate RZ countour
687  if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
688  {
689  std::ostringstream message;
690  message << "Triangulation of RZ contour has failed for solid: "
691  << GetName() << " !"
692  << "\nExtent has been calculated using boundary box";
693  G4Exception("G4Polyhedra::CalculateExtent()",
694  "GeomMgt1002",JustWarning,message);
695  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
696  }
697 
698  // set trigonometric values
699  G4double sphi = GetStartPhi();
700  G4double ephi = GetEndPhi();
701  G4double dphi = IsOpen() ? ephi-sphi : twopi;
702  G4int ksteps = GetNumSide();
703  G4double astep = dphi/ksteps;
704  G4double sinStep = std::sin(astep);
705  G4double cosStep = std::cos(astep);
706  G4double sinStart = GetSinStartPhi();
707  G4double cosStart = GetCosStartPhi();
708 
709  // allocate vector lists
710  std::vector<const G4ThreeVectorList *> polygons;
711  polygons.resize(ksteps+1);
712  for (G4int k=0; k<ksteps+1; ++k) {
713  polygons[k] = new G4ThreeVectorList(3);
714  }
715 
716  // main loop along triangles
717  pMin = kInfinity;
718  pMax = -kInfinity;
719  G4int ntria = triangles.size()/3;
720  for (G4int i=0; i<ntria; ++i)
721  {
722  G4double sinCur = sinStart;
723  G4double cosCur = cosStart;
724  G4int i3 = i*3;
725  for (G4int k=0; k<ksteps+1; ++k) // rotate triangle
726  {
727  G4ThreeVectorList* ptr = const_cast<G4ThreeVectorList*>(polygons[k]);
728  G4ThreeVectorList::iterator iter = ptr->begin();
729  iter->set(triangles[i3+0].x()*cosCur,
730  triangles[i3+0].x()*sinCur,
731  triangles[i3+0].y());
732  iter++;
733  iter->set(triangles[i3+1].x()*cosCur,
734  triangles[i3+1].x()*sinCur,
735  triangles[i3+1].y());
736  iter++;
737  iter->set(triangles[i3+2].x()*cosCur,
738  triangles[i3+2].x()*sinCur,
739  triangles[i3+2].y());
740 
741  G4double sinTmp = sinCur;
742  sinCur = sinCur*cosStep + cosCur*sinStep;
743  cosCur = cosCur*cosStep - sinTmp*sinStep;
744  }
745 
746  // set sub-envelope and adjust extent
747  G4double emin,emax;
748  G4BoundingEnvelope benv(polygons);
749  if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
750  if (emin < pMin) pMin = emin;
751  if (emax > pMax) pMax = emax;
752  if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
753  }
754  // free memory
755  for (G4int k=0; k<ksteps+1; ++k) { delete polygons[k]; polygons[k]=0;}
756  return (pMin < pMax);
757 }
758 
759 //
760 // ComputeDimensions
761 //
763  const G4int n,
764  const G4VPhysicalVolume* pRep )
765 {
766  p->ComputeDimensions(*this,n,pRep);
767 }
768 
769 
770 //
771 // GetEntityType
772 //
774 {
775  return G4String("G4Polyhedra");
776 }
777 
778 
779 //
780 // Make a clone of the object
781 //
783 {
784  return new G4Polyhedra(*this);
785 }
786 
787 
788 //
789 // Stream object contents to an output stream
790 //
791 std::ostream& G4Polyhedra::StreamInfo( std::ostream& os ) const
792 {
793  G4int oldprc = os.precision(16);
794  os << "-----------------------------------------------------------\n"
795  << " *** Dump for solid - " << GetName() << " ***\n"
796  << " ===================================================\n"
797  << " Solid type: G4Polyhedra\n"
798  << " Parameters: \n"
799  << " starting phi angle : " << startPhi/degree << " degrees \n"
800  << " ending phi angle : " << endPhi/degree << " degrees \n"
801  << " number of sides : " << numSide << " \n";
802  G4int i=0;
803  if (!genericPgon)
804  {
806  os << " number of Z planes: " << numPlanes << "\n"
807  << " Z values: \n";
808  for (i=0; i<numPlanes; i++)
809  {
810  os << " Z plane " << i << ": "
811  << original_parameters->Z_values[i] << "\n";
812  }
813  os << " Tangent distances to inner surface (Rmin): \n";
814  for (i=0; i<numPlanes; i++)
815  {
816  os << " Z plane " << i << ": "
817  << original_parameters->Rmin[i] << "\n";
818  }
819  os << " Tangent distances to outer surface (Rmax): \n";
820  for (i=0; i<numPlanes; i++)
821  {
822  os << " Z plane " << i << ": "
823  << original_parameters->Rmax[i] << "\n";
824  }
825  }
826  os << " number of RZ points: " << numCorner << "\n"
827  << " RZ values (corners): \n";
828  for (i=0; i<numCorner; i++)
829  {
830  os << " "
831  << corners[i].r << ", " << corners[i].z << "\n";
832  }
833  os << "-----------------------------------------------------------\n";
834  os.precision(oldprc);
835 
836  return os;
837 }
838 
839 
840 //
841 // GetPointOnPlane
842 //
843 // Auxiliary method for get point on surface
844 //
847  G4ThreeVector p2, G4ThreeVector p3) const
848 {
849  G4double lambda1, lambda2, chose,aOne,aTwo;
850  G4ThreeVector t, u, v, w, Area, normal;
851  aOne = 1.;
852  aTwo = 1.;
853 
854  t = p1 - p0;
855  u = p2 - p1;
856  v = p3 - p2;
857  w = p0 - p3;
858 
859  chose = G4RandFlat::shoot(0.,aOne+aTwo);
860  if( (chose>=0.) && (chose < aOne) )
861  {
862  lambda1 = G4RandFlat::shoot(0.,1.);
863  lambda2 = G4RandFlat::shoot(0.,lambda1);
864  return (p2+lambda1*v+lambda2*w);
865  }
866 
867  lambda1 = G4RandFlat::shoot(0.,1.);
868  lambda2 = G4RandFlat::shoot(0.,lambda1);
869  return (p0+lambda1*t+lambda2*u);
870 }
871 
872 
873 //
874 // GetPointOnTriangle
875 //
876 // Auxiliary method for get point on surface
877 //
879  G4ThreeVector p2,
880  G4ThreeVector p3) const
881 {
882  G4double lambda1,lambda2;
883  G4ThreeVector v=p3-p1, w=p1-p2;
884 
885  lambda1 = G4RandFlat::shoot(0.,1.);
886  lambda2 = G4RandFlat::shoot(0.,lambda1);
887 
888  return (p2 + lambda1*w + lambda2*v);
889 }
890 
891 
892 //
893 // GetPointOnSurface
894 //
896 {
897  if( !genericPgon ) // Polyhedra by faces
898  {
899  G4int j, numPlanes = original_parameters->Num_z_planes, Flag=0;
900  G4double chose, totArea=0., Achose1, Achose2,
901  rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
902  G4double a, b, l2, rang, totalPhi, ksi,
903  area, aTop=0., aBottom=0., zVal=0.;
904 
905  G4ThreeVector p0, p1, p2, p3;
906  std::vector<G4double> aVector1;
907  std::vector<G4double> aVector2;
908  std::vector<G4double> aVector3;
909 
910  totalPhi= (phiIsOpen) ? (endPhi-startPhi) : twopi;
911  ksi = totalPhi/numSide;
912  G4double cosksi = std::cos(ksi/2.);
913 
914  // Below we generate the areas relevant to our solid
915  //
916  for(j=0; j<numPlanes-1; j++)
917  {
918  a = original_parameters->Rmax[j+1];
919  b = original_parameters->Rmax[j];
921  -original_parameters->Z_values[j+1]) + sqr(b-a);
922  area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
923  aVector1.push_back(area);
924  }
925 
926  for(j=0; j<numPlanes-1; j++)
927  {
928  a = original_parameters->Rmin[j+1];//*cosksi;
929  b = original_parameters->Rmin[j];//*cosksi;
931  -original_parameters->Z_values[j+1]) + sqr(b-a);
932  area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
933  aVector2.push_back(area);
934  }
935 
936  for(j=0; j<numPlanes-1; j++)
937  {
938  if(phiIsOpen == true)
939  {
940  aVector3.push_back(0.5*(original_parameters->Rmax[j]
943  -original_parameters->Rmin[j+1])
944  *std::fabs(original_parameters->Z_values[j+1]
946  }
947  else { aVector3.push_back(0.); }
948  }
949 
950  for(j=0; j<numPlanes-1; j++)
951  {
952  totArea += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
953  }
954 
955  // Must include top and bottom areas
956  //
957  if(original_parameters->Rmax[numPlanes-1] != 0.)
958  {
959  a = original_parameters->Rmax[numPlanes-1];
960  b = original_parameters->Rmin[numPlanes-1];
961  l2 = sqr(a-b);
962  aTop = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
963  }
964 
965  if(original_parameters->Rmax[0] != 0.)
966  {
967  a = original_parameters->Rmax[0];
968  b = original_parameters->Rmin[0];
969  l2 = sqr(a-b);
970  aBottom = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
971  }
972 
973  Achose1 = 0.;
974  Achose2 = numSide*(aVector1[0]+aVector2[0])+2.*aVector3[0];
975 
976  chose = G4RandFlat::shoot(0.,totArea+aTop+aBottom);
977  if( (chose >= 0.) && (chose < aTop + aBottom) )
978  {
979  chose = G4RandFlat::shoot(startPhi,startPhi+totalPhi);
980  rang = std::floor((chose-startPhi)/ksi-0.01);
981  if(rang<0) { rang=0; }
982  rang = std::fabs(rang);
983  sinphi1 = std::sin(startPhi+rang*ksi);
984  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
985  cosphi1 = std::cos(startPhi+rang*ksi);
986  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
987  chose = G4RandFlat::shoot(0., aTop + aBottom);
988  if(chose>=0. && chose<aTop)
989  {
990  rad1 = original_parameters->Rmin[numPlanes-1];
991  rad2 = original_parameters->Rmax[numPlanes-1];
992  zVal = original_parameters->Z_values[numPlanes-1];
993  }
994  else
995  {
996  rad1 = original_parameters->Rmin[0];
997  rad2 = original_parameters->Rmax[0];
998  zVal = original_parameters->Z_values[0];
999  }
1000  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
1001  p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
1002  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
1003  p3 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
1004  return GetPointOnPlane(p0,p1,p2,p3);
1005  }
1006  else
1007  {
1008  for (j=0; j<numPlanes-1; j++)
1009  {
1010  if( ((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes-2) )
1011  {
1012  Flag = j; break;
1013  }
1014  Achose1 += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
1015  Achose2 = Achose1 + numSide*(aVector1[j+1]+aVector2[j+1])
1016  + 2.*aVector3[j+1];
1017  }
1018  }
1019 
1020  // At this point we have chosen a subsection
1021  // between to adjacent plane cuts...
1022 
1023  j = Flag;
1024 
1025  totArea = numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
1026  chose = G4RandFlat::shoot(0.,totArea);
1027 
1028  if( (chose>=0.) && (chose<numSide*aVector1[j]) )
1029  {
1030  chose = G4RandFlat::shoot(startPhi,startPhi+totalPhi);
1031  rang = std::floor((chose-startPhi)/ksi-0.01);
1032  if(rang<0) { rang=0; }
1033  rang = std::fabs(rang);
1034  rad1 = original_parameters->Rmax[j];
1035  rad2 = original_parameters->Rmax[j+1];
1036  sinphi1 = std::sin(startPhi+rang*ksi);
1037  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
1038  cosphi1 = std::cos(startPhi+rang*ksi);
1039  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
1040  zVal = original_parameters->Z_values[j];
1041 
1042  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
1043  p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
1044 
1045  zVal = original_parameters->Z_values[j+1];
1046 
1047  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
1048  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
1049  return GetPointOnPlane(p0,p1,p2,p3);
1050  }
1051  else if ( (chose >= numSide*aVector1[j])
1052  && (chose <= numSide*(aVector1[j]+aVector2[j])) )
1053  {
1054  chose = G4RandFlat::shoot(startPhi,startPhi+totalPhi);
1055  rang = std::floor((chose-startPhi)/ksi-0.01);
1056  if(rang<0) { rang=0; }
1057  rang = std::fabs(rang);
1058  rad1 = original_parameters->Rmin[j];
1059  rad2 = original_parameters->Rmin[j+1];
1060  sinphi1 = std::sin(startPhi+rang*ksi);
1061  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
1062  cosphi1 = std::cos(startPhi+rang*ksi);
1063  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
1064  zVal = original_parameters->Z_values[j];
1065 
1066  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
1067  p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
1068 
1069  zVal = original_parameters->Z_values[j+1];
1070 
1071  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
1072  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
1073  return GetPointOnPlane(p0,p1,p2,p3);
1074  }
1075 
1076  chose = G4RandFlat::shoot(0.,2.2);
1077  if( (chose>=0.) && (chose < 1.) )
1078  {
1079  rang = startPhi;
1080  }
1081  else
1082  {
1083  rang = endPhi;
1084  }
1085 
1086  cosphi1 = std::cos(rang); rad1 = original_parameters->Rmin[j];
1087  sinphi1 = std::sin(rang); rad2 = original_parameters->Rmax[j];
1088 
1089  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
1091  p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
1093 
1094  rad1 = original_parameters->Rmax[j+1];
1095  rad2 = original_parameters->Rmin[j+1];
1096 
1097  p2 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
1099  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
1101  return GetPointOnPlane(p0,p1,p2,p3);
1102  }
1103  else // Generic polyhedra
1104  {
1105  return GetPointOnSurfaceGeneric();
1106  }
1107 }
1108 
1109 //
1110 // CreatePolyhedron
1111 //
1113 {
1114  if (!genericPgon)
1115  {
1123  }
1124  else
1125  {
1126  // The following code prepares for:
1127  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
1128  // const double xyz[][3],
1129  // const int faces_vec[][4])
1130  // Here is an extract from the header file HepPolyhedron.h:
1148  G4int nNodes;
1149  G4int nFaces;
1150  typedef G4double double3[3];
1151  double3* xyz;
1152  typedef G4int int4[4];
1153  int4* faces_vec;
1154  if (phiIsOpen)
1155  {
1156  // Triangulate open ends. Simple ear-chopping algorithm...
1157  // I'm not sure how robust this algorithm is (J.Allison).
1158  //
1159  std::vector<G4bool> chopped(numCorner, false);
1160  std::vector<G4int*> triQuads;
1161  G4int remaining = numCorner;
1162  G4int iStarter = 0;
1163  while (remaining >= 3) // Loop checking, 13.08.2015, G.Cosmo
1164  {
1165  // Find unchopped corners...
1166  //
1167  G4int A = -1, B = -1, C = -1;
1168  G4int iStepper = iStarter;
1169  do // Loop checking, 13.08.2015, G.Cosmo
1170  {
1171  if (A < 0) { A = iStepper; }
1172  else if (B < 0) { B = iStepper; }
1173  else if (C < 0) { C = iStepper; }
1174  do // Loop checking, 13.08.2015, G.Cosmo
1175  {
1176  if (++iStepper >= numCorner) iStepper = 0;
1177  }
1178  while (chopped[iStepper]);
1179  }
1180  while (C < 0 && iStepper != iStarter);
1181 
1182  // Check triangle at B is pointing outward (an "ear").
1183  // Sign of z cross product determines...
1184 
1185  G4double BAr = corners[A].r - corners[B].r;
1186  G4double BAz = corners[A].z - corners[B].z;
1187  G4double BCr = corners[C].r - corners[B].r;
1188  G4double BCz = corners[C].z - corners[B].z;
1189  if (BAr * BCz - BAz * BCr < kCarTolerance)
1190  {
1191  G4int* tq = new G4int[3];
1192  tq[0] = A + 1;
1193  tq[1] = B + 1;
1194  tq[2] = C + 1;
1195  triQuads.push_back(tq);
1196  chopped[B] = true;
1197  --remaining;
1198  }
1199  else
1200  {
1201  do // Loop checking, 13.08.2015, G.Cosmo
1202  {
1203  if (++iStarter >= numCorner) { iStarter = 0; }
1204  }
1205  while (chopped[iStarter]);
1206  }
1207  }
1208 
1209  // Transfer to faces...
1210 
1211  nNodes = (numSide + 1) * numCorner;
1212  nFaces = numSide * numCorner + 2 * triQuads.size();
1213  faces_vec = new int4[nFaces];
1214  G4int iface = 0;
1215  G4int addition = numCorner * numSide;
1216  G4int d = numCorner - 1;
1217  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
1218  {
1219  for (size_t i = 0; i < triQuads.size(); ++i)
1220  {
1221  // Negative for soft/auxiliary/normally invisible edges...
1222  //
1223  G4int a, b, c;
1224  if (iEnd == 0)
1225  {
1226  a = triQuads[i][0];
1227  b = triQuads[i][1];
1228  c = triQuads[i][2];
1229  }
1230  else
1231  {
1232  a = triQuads[i][0] + addition;
1233  b = triQuads[i][2] + addition;
1234  c = triQuads[i][1] + addition;
1235  }
1236  G4int ab = std::abs(b - a);
1237  G4int bc = std::abs(c - b);
1238  G4int ca = std::abs(a - c);
1239  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
1240  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
1241  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
1242  faces_vec[iface][3] = 0;
1243  ++iface;
1244  }
1245  }
1246 
1247  // Continue with sides...
1248 
1249  xyz = new double3[nNodes];
1250  const G4double dPhi = (endPhi - startPhi) / numSide;
1251  G4double phi = startPhi;
1252  G4int ixyz = 0;
1253  for (G4int iSide = 0; iSide < numSide; ++iSide)
1254  {
1255  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1256  {
1257  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1258  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1259  xyz[ixyz][2] = corners[iCorner].z;
1260  if (iCorner < numCorner - 1)
1261  {
1262  faces_vec[iface][0] = ixyz + 1;
1263  faces_vec[iface][1] = ixyz + numCorner + 1;
1264  faces_vec[iface][2] = ixyz + numCorner + 2;
1265  faces_vec[iface][3] = ixyz + 2;
1266  }
1267  else
1268  {
1269  faces_vec[iface][0] = ixyz + 1;
1270  faces_vec[iface][1] = ixyz + numCorner + 1;
1271  faces_vec[iface][2] = ixyz + 2;
1272  faces_vec[iface][3] = ixyz - numCorner + 2;
1273  }
1274  ++iface;
1275  ++ixyz;
1276  }
1277  phi += dPhi;
1278  }
1279 
1280  // Last corners...
1281 
1282  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1283  {
1284  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1285  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1286  xyz[ixyz][2] = corners[iCorner].z;
1287  ++ixyz;
1288  }
1289  }
1290  else // !phiIsOpen - i.e., a complete 360 degrees.
1291  {
1292  nNodes = numSide * numCorner;
1293  nFaces = numSide * numCorner;;
1294  xyz = new double3[nNodes];
1295  faces_vec = new int4[nFaces];
1296  // const G4double dPhi = (endPhi - startPhi) / numSide;
1297  const G4double dPhi = twopi / numSide;
1298  // !phiIsOpen endPhi-startPhi = 360 degrees.
1299  G4double phi = startPhi;
1300  G4int ixyz = 0, iface = 0;
1301  for (G4int iSide = 0; iSide < numSide; ++iSide)
1302  {
1303  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1304  {
1305  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1306  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1307  xyz[ixyz][2] = corners[iCorner].z;
1308  if (iSide < numSide - 1)
1309  {
1310  if (iCorner < numCorner - 1)
1311  {
1312  faces_vec[iface][0] = ixyz + 1;
1313  faces_vec[iface][1] = ixyz + numCorner + 1;
1314  faces_vec[iface][2] = ixyz + numCorner + 2;
1315  faces_vec[iface][3] = ixyz + 2;
1316  }
1317  else
1318  {
1319  faces_vec[iface][0] = ixyz + 1;
1320  faces_vec[iface][1] = ixyz + numCorner + 1;
1321  faces_vec[iface][2] = ixyz + 2;
1322  faces_vec[iface][3] = ixyz - numCorner + 2;
1323  }
1324  }
1325  else // Last side joins ends...
1326  {
1327  if (iCorner < numCorner - 1)
1328  {
1329  faces_vec[iface][0] = ixyz + 1;
1330  faces_vec[iface][1] = ixyz + numCorner - nFaces + 1;
1331  faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
1332  faces_vec[iface][3] = ixyz + 2;
1333  }
1334  else
1335  {
1336  faces_vec[iface][0] = ixyz + 1;
1337  faces_vec[iface][1] = ixyz - nFaces + numCorner + 1;
1338  faces_vec[iface][2] = ixyz - nFaces + 2;
1339  faces_vec[iface][3] = ixyz - numCorner + 2;
1340  }
1341  }
1342  ++ixyz;
1343  ++iface;
1344  }
1345  phi += dPhi;
1346  }
1347  }
1348  G4Polyhedron* polyhedron = new G4Polyhedron;
1349  G4int problem = polyhedron->createPolyhedron(nNodes,nFaces,xyz,faces_vec);
1350  delete [] faces_vec;
1351  delete [] xyz;
1352  if (problem)
1353  {
1354  std::ostringstream message;
1355  message << "Problem creating G4Polyhedron for: " << GetName();
1356  G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
1357  JustWarning, message);
1358  delete polyhedron;
1359  return 0;
1360  }
1361  else
1362  {
1363  return polyhedron;
1364  }
1365  }
1366 }
1367 
1368 
1370 {
1371  G4int numPlanes = (G4int)numCorner;
1372  G4bool isConvertible=true;
1373  G4double Zmax=rz->Bmax();
1374  rz->StartWithZMin();
1375 
1376  // Prepare vectors for storage
1377  //
1378  std::vector<G4double> Z;
1379  std::vector<G4double> Rmin;
1380  std::vector<G4double> Rmax;
1381 
1382  G4int countPlanes=1;
1383  G4int icurr=0;
1384  G4int icurl=0;
1385 
1386  // first plane Z=Z[0]
1387  //
1388  Z.push_back(corners[0].z);
1389  G4double Zprev=Z[0];
1390  if (Zprev == corners[1].z)
1391  {
1392  Rmin.push_back(corners[0].r);
1393  Rmax.push_back (corners[1].r);icurr=1;
1394  }
1395  else if (Zprev == corners[numPlanes-1].z)
1396  {
1397  Rmin.push_back(corners[numPlanes-1].r);
1398  Rmax.push_back (corners[0].r);
1399  icurl=numPlanes-1;
1400  }
1401  else
1402  {
1403  Rmin.push_back(corners[0].r);
1404  Rmax.push_back (corners[0].r);
1405  }
1406 
1407  // next planes until last
1408  //
1409  G4int inextr=0, inextl=0;
1410  for (G4int i=0; i < numPlanes-2; i++)
1411  {
1412  inextr=1+icurr;
1413  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1414 
1415  if((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax)) { break; }
1416 
1417  G4double Zleft = corners[inextl].z;
1418  G4double Zright = corners[inextr].z;
1419  if(Zright>Zleft)
1420  {
1421  Z.push_back(Zleft);
1422  countPlanes++;
1423  G4double difZr=corners[inextr].z - corners[icurr].z;
1424  G4double difZl=corners[inextl].z - corners[icurl].z;
1425 
1426  if(std::fabs(difZl) < kCarTolerance)
1427  {
1428  if(std::fabs(difZr) < kCarTolerance)
1429  {
1430  Rmin.push_back(corners[inextl].r);
1431  Rmax.push_back(corners[icurr].r);
1432  }
1433  else
1434  {
1435  Rmin.push_back(corners[inextl].r);
1436  Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1437  *(corners[inextr].r - corners[icurr].r));
1438  }
1439  }
1440  else if (difZl >= kCarTolerance)
1441  {
1442  if(std::fabs(difZr) < kCarTolerance)
1443  {
1444  Rmin.push_back(corners[icurl].r);
1445  Rmax.push_back(corners[icurr].r);
1446  }
1447  else
1448  {
1449  Rmin.push_back(corners[icurl].r);
1450  Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1451  *(corners[inextr].r - corners[icurr].r));
1452  }
1453  }
1454  else
1455  {
1456  isConvertible=false; break;
1457  }
1458  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1459  }
1460  else if(std::fabs(Zright-Zleft)<kCarTolerance) // Zright=Zleft
1461  {
1462  Z.push_back(Zleft);
1463  countPlanes++;
1464  icurr++;
1465 
1466  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1467 
1468  Rmin.push_back(corners[inextl].r);
1469  Rmax.push_back (corners[inextr].r);
1470  }
1471  else // Zright<Zleft
1472  {
1473  Z.push_back(Zright);
1474  countPlanes++;
1475 
1476  G4double difZr=corners[inextr].z - corners[icurr].z;
1477  G4double difZl=corners[inextl].z - corners[icurl].z;
1478  if(std::fabs(difZr) < kCarTolerance)
1479  {
1480  if(std::fabs(difZl) < kCarTolerance)
1481  {
1482  Rmax.push_back(corners[inextr].r);
1483  Rmin.push_back(corners[icurr].r);
1484  }
1485  else
1486  {
1487  Rmin.push_back(corners[icurl].r + (Zright-corners[icurl].z)/difZl
1488  * (corners[inextl].r - corners[icurl].r));
1489  Rmax.push_back(corners[inextr].r);
1490  }
1491  icurr++;
1492  } // plate
1493  else if (difZr >= kCarTolerance)
1494  {
1495  if(std::fabs(difZl) < kCarTolerance)
1496  {
1497  Rmax.push_back(corners[inextr].r);
1498  Rmin.push_back (corners[icurr].r);
1499  }
1500  else
1501  {
1502  Rmax.push_back(corners[inextr].r);
1503  Rmin.push_back (corners[icurl].r+(Zright-corners[icurl].z)/difZl
1504  * (corners[inextl].r - corners[icurl].r));
1505  }
1506  icurr++;
1507  }
1508  else
1509  {
1510  isConvertible=false; break;
1511  }
1512  }
1513  } // end for loop
1514 
1515  // last plane Z=Zmax
1516  //
1517  Z.push_back(Zmax);
1518  countPlanes++;
1519  inextr=1+icurr;
1520  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1521 
1522  Rmax.push_back(corners[inextr].r);
1523  Rmin.push_back(corners[inextl].r);
1524 
1525  // Set original parameters Rmin,Rmax,Z
1526  //
1527  if(isConvertible)
1528  {
1531  original_parameters->Z_values = new G4double[countPlanes];
1532  original_parameters->Rmin = new G4double[countPlanes];
1533  original_parameters->Rmax = new G4double[countPlanes];
1534 
1535  for(G4int j=0; j < countPlanes; j++)
1536  {
1537  original_parameters->Z_values[j] = Z[j];
1538  original_parameters->Rmax[j] = Rmax[j];
1539  original_parameters->Rmin[j] = Rmin[j];
1540  }
1543  original_parameters->Num_z_planes = countPlanes;
1544 
1545  }
1546  else // Set parameters(r,z) with Rmin==0 as convention
1547  {
1548 #ifdef G4SPECSDEBUG
1549  std::ostringstream message;
1550  message << "Polyhedra " << GetName() << G4endl
1551  << "cannot be converted to Polyhedra with (Rmin,Rmaz,Z) parameters!";
1552  G4Exception("G4Polyhedra::SetOriginalParameters()",
1553  "GeomSolids0002", JustWarning, message);
1554 #endif
1557  original_parameters->Z_values = new G4double[numPlanes];
1558  original_parameters->Rmin = new G4double[numPlanes];
1559  original_parameters->Rmax = new G4double[numPlanes];
1560 
1561  for(G4int j=0; j < numPlanes; j++)
1562  {
1564  original_parameters->Rmax[j] = corners[j].r;
1565  original_parameters->Rmin[j] = 0.0;
1566  }
1569  original_parameters->Num_z_planes = numPlanes;
1570  }
1571  //return isConvertible;
1572 }
1573 
1574 #endif
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
Float_t x
Definition: compare.C:6
void set(double x, double y, double z)
const XML_Char * name
Definition: expat.h:151
static G4bool TriangulatePolygon(const G4TwoVectorList &polygon, G4TwoVectorList &result)
Definition: G4GeomTools.cc:199
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
CLHEP::Hep3Vector G4ThreeVector
G4bool Reset()
Definition: G4Polyhedra.cc:486
Double_t xx
G4Polyhedron * CreatePolyhedron() const
G4double GetMaxExtent(const EAxis pAxis) const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Polyhedra.cc:548
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4GeometryType GetEntityType() const
Definition: G4Polyhedra.cc:773
void CopyStuff(const G4Polyhedra &source)
Definition: G4Polyhedra.cc:437
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
const char * p
Definition: xmltok.h:285
CLHEP::Hep2Vector G4TwoVector
Definition: G4TwoVector.hh:42
void ScaleA(G4double scale)
Double_t z
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
G4bool phiIsOpen
Definition: G4Polyhedra.hh:193
G4bool IsOpen() const
double z() const
G4double GetSinStartPhi() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
static const G4double emax
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Polyhedra.cc:646
G4bool ShouldMiss(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool MustBeOutside(const G4ThreeVector &p) const
virtual EInside Inside(const G4ThreeVector &p) const
void SetOriginalParameters(G4PolyhedraHistorical *pars)
int Zmax
Float_t Z
G4Polyhedra & operator=(const G4Polyhedra &source)
Definition: G4Polyhedra.cc:417
static const G4double ab
EInside Inside(const G4ThreeVector &p) const
Definition: G4Polyhedra.cc:528
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
std::vector< G4ThreeVector > G4ThreeVectorList
G4double endPhi
Definition: G4Polyhedra.hh:192
G4double GetStartPhi() const
G4double GetMinExtent(const EAxis pAxis) const
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
G4bool RemoveDuplicateVertices(G4double tolerance)
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4PolyhedraHistorical * original_parameters
Definition: G4Polyhedra.hh:197
static constexpr double twopi
Definition: G4SIunits.hh:76
G4int numCorner
Definition: G4Polyhedra.hh:195
G4double GetEndPhi() const
double A(double temperature)
#define DBL_EPSILON
Definition: templates.hh:87
Float_t d
G4double Amin() const
G4VSolid * Clone() const
Definition: G4Polyhedra.cc:782
G4bool CrossesItself(G4double tolerance)
ThreeVector shoot(const G4int Ap, const G4int Af)
G4VCSGface ** faces
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Polyhedra.cc:762
G4PolyhedraSideRZ GetCorner(const G4int index) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4bool genericPgon
Definition: G4Polyhedra.hh:194
int G4int
Definition: G4Types.hh:78
EInside
Definition: geomdefs.hh:58
G4int NumVertices() const
std::vector< G4TwoVector > G4TwoVectorList
Definition: G4GeomTools.hh:50
G4bool RemoveRedundantVertices(G4double tolerance)
G4double startPhi
Definition: G4Polyhedra.hh:191
EAxis
Definition: geomdefs.hh:54
G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
G4String GetName() const
void DumpInfo() const
double C(double temp)
G4double GetCosStartPhi() const
G4int GetNumRZCorner() const
G4ThreeVector GetPointOnSurfaceGeneric() const
static constexpr double degree
Definition: G4SIunits.hh:144
G4ThreeVector GetPointOnTriangle(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2) const
Definition: G4Polyhedra.cc:878
double x() const
static void RemoveRedundantVertices(G4TwoVectorList &polygon, std::vector< G4int > &iout, G4double tolerance=0)
Definition: G4GeomTools.cc:311
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
virtual ~G4Polyhedra()
Definition: G4Polyhedra.cc:395
Char_t n[5]
G4PolyhedraSideRZ * corners
Definition: G4Polyhedra.hh:196
T sqr(const T &x)
Definition: templates.hh:145
double y() const
G4ThreeVector GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3) const
Definition: G4Polyhedra.cc:846
G4double Bmax() const
void Create(G4double phiStart, G4double phiTotal, G4int numSide, G4ReduciblePolygon *rz)
Definition: G4Polyhedra.cc:206
double B(double temperature)
G4bool fRebuildPolyhedron
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Polyhedra.cc:791
G4EnclosingCylinder * enclosingCylinder
Definition: G4Polyhedra.hh:199
G4ThreeVector GetPointOnSurface() const
Definition: G4Polyhedra.cc:895
G4Polyhedron * fpPolyhedron
static G4double PolygonArea(const G4TwoVectorList &polygon)
Definition: G4GeomTools.cc:82
G4int GetNumSide() const
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4Polyhedra.cc:576