Geant4  10.01.p03
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 95301 2016-02-04 11:25:33Z 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 "Randomize.hh"
65 
66 #include "G4EnclosingCylinder.hh"
67 #include "G4ReduciblePolygon.hh"
68 #include "G4VPVParameterisation.hh"
69 
70 #include <sstream>
71 
72 using namespace CLHEP;
73 
74 //
75 // Constructor (GEANT3 style parameters)
76 //
77 // GEANT3 PGON radii are specified in the distance to the norm of each face.
78 //
80  G4double phiStart,
81  G4double thePhiTotal,
82  G4int theNumSide,
83  G4int numZPlanes,
84  const G4double zPlane[],
85  const G4double rInner[],
86  const G4double rOuter[] )
87  : G4VCSGfaceted( name ), genericPgon(false)
88 {
89  if (theNumSide <= 0)
90  {
91  std::ostringstream message;
92  message << "Solid must have at least one side - " << GetName() << G4endl
93  << " No sides specified !";
94  G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
95  FatalErrorInArgument, message);
96  }
97 
98  //
99  // Calculate conversion factor from G3 radius to G4 radius
100  //
101  G4double phiTotal = thePhiTotal;
102  if ( (phiTotal <=0) || (phiTotal >= twopi*(1-DBL_EPSILON)) )
103  { phiTotal = twopi; }
104  G4double convertRad = std::cos(0.5*phiTotal/theNumSide);
105 
106  //
107  // Some historical stuff
108  //
110 
111  original_parameters->numSide = theNumSide;
112  original_parameters->Start_angle = phiStart;
114  original_parameters->Num_z_planes = numZPlanes;
115  original_parameters->Z_values = new G4double[numZPlanes];
116  original_parameters->Rmin = new G4double[numZPlanes];
117  original_parameters->Rmax = new G4double[numZPlanes];
118 
119  G4int i;
120  for (i=0; i<numZPlanes; i++)
121  {
122  if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
123  {
124  if( (rInner[i] > rOuter[i+1])
125  ||(rInner[i+1] > rOuter[i]) )
126  {
127  DumpInfo();
128  std::ostringstream message;
129  message << "Cannot create a Polyhedra with no contiguous segments."
130  << G4endl
131  << " Segments are not contiguous !" << G4endl
132  << " rMin[" << i << "] = " << rInner[i]
133  << " -- rMax[" << i+1 << "] = " << rOuter[i+1] << G4endl
134  << " rMin[" << i+1 << "] = " << rInner[i+1]
135  << " -- rMax[" << i << "] = " << rOuter[i];
136  G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
137  FatalErrorInArgument, message);
138  }
139  }
140  original_parameters->Z_values[i] = zPlane[i];
141  original_parameters->Rmin[i] = rInner[i]/convertRad;
142  original_parameters->Rmax[i] = rOuter[i]/convertRad;
143  }
144 
145 
146  //
147  // Build RZ polygon using special PCON/PGON GEANT3 constructor
148  //
149  G4ReduciblePolygon *rz =
150  new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
151  rz->ScaleA( 1/convertRad );
152 
153  //
154  // Do the real work
155  //
156  Create( phiStart, phiTotal, theNumSide, rz );
157 
158  delete rz;
159 }
160 
161 
162 //
163 // Constructor (generic parameters)
164 //
166  G4double phiStart,
167  G4double phiTotal,
168  G4int theNumSide,
169  G4int numRZ,
170  const G4double r[],
171  const G4double z[] )
172  : G4VCSGfaceted( name ), genericPgon(true)
173 {
174  if (theNumSide <= 0)
175  {
176  std::ostringstream message;
177  message << "Solid must have at least one side - " << GetName() << G4endl
178  << " No sides specified !";
179  G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
180  FatalErrorInArgument, message);
181  }
182 
183  G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
184 
185  Create( phiStart, phiTotal, theNumSide, rz );
186 
187  // Set original_parameters struct for consistency
188  //
190 
191  delete rz;
192 }
193 
194 
195 //
196 // Create
197 //
198 // Generic create routine, called by each constructor
199 // after conversion of arguments
200 //
202  G4double phiTotal,
203  G4int theNumSide,
204  G4ReduciblePolygon *rz )
205 {
206  //
207  // Perform checks of rz values
208  //
209  if (rz->Amin() < 0.0)
210  {
211  std::ostringstream message;
212  message << "Illegal input parameters - " << GetName() << G4endl
213  << " All R values must be >= 0 !";
214  G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
215  FatalErrorInArgument, message);
216  }
217 
218  G4double rzArea = rz->Area();
219  if (rzArea < -kCarTolerance)
220  rz->ReverseOrder();
221 
222  else if (rzArea < -kCarTolerance)
223  {
224  std::ostringstream message;
225  message << "Illegal input parameters - " << GetName() << G4endl
226  << " R/Z cross section is zero or near zero: " << rzArea;
227  G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
228  FatalErrorInArgument, message);
229  }
230 
232  || (!rz->RemoveRedundantVertices( kCarTolerance )) )
233  {
234  std::ostringstream message;
235  message << "Illegal input parameters - " << GetName() << G4endl
236  << " Too few unique R/Z values !";
237  G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
238  FatalErrorInArgument, message);
239  }
240 
241  if (rz->CrossesItself( 1/kInfinity ))
242  {
243  std::ostringstream message;
244  message << "Illegal input parameters - " << GetName() << G4endl
245  << " R/Z segments cross !";
246  G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
247  FatalErrorInArgument, message);
248  }
249 
250  numCorner = rz->NumVertices();
251 
252 
253  startPhi = phiStart;
254  while( startPhi < 0 ) // Loop checking, 13.08.2015, G.Cosmo
255  startPhi += twopi;
256  //
257  // Phi opening? Account for some possible roundoff, and interpret
258  // nonsense value as representing no phi opening
259  //
260  if ( (phiTotal <= 0) || (phiTotal > twopi*(1-DBL_EPSILON)) )
261  {
262  phiIsOpen = false;
263  endPhi = phiStart+twopi;
264  }
265  else
266  {
267  phiIsOpen = true;
268 
269  //
270  // Convert phi into our convention
271  //
272  endPhi = phiStart+phiTotal;
273  while( endPhi < startPhi ) // Loop checking, 13.08.2015, G.Cosmo
274  endPhi += twopi;
275  }
276 
277  //
278  // Save number sides
279  //
280  numSide = theNumSide;
281 
282  //
283  // Allocate corner array.
284  //
286 
287  //
288  // Copy corners
289  //
290  G4ReduciblePolygonIterator iterRZ(rz);
291 
292  G4PolyhedraSideRZ *next = corners;
293  iterRZ.Begin();
294  do // Loop checking, 13.08.2015, G.Cosmo
295  {
296  next->r = iterRZ.GetA();
297  next->z = iterRZ.GetB();
298  } while( ++next, iterRZ.Next() );
299 
300  //
301  // Allocate face pointer array
302  //
304  faces = new G4VCSGface*[numFace];
305 
306  //
307  // Construct side faces
308  //
309  // To do so properly, we need to keep track of four successive RZ
310  // corners.
311  //
312  // But! Don't construct a face if both points are at zero radius!
313  //
314  G4PolyhedraSideRZ *corner = corners,
315  *prev = corners + numCorner-1,
316  *nextNext;
317  G4VCSGface **face = faces;
318  do // Loop checking, 13.08.2015, G.Cosmo
319  {
320  next = corner+1;
321  if (next >= corners+numCorner) next = corners;
322  nextNext = next+1;
323  if (nextNext >= corners+numCorner) nextNext = corners;
324 
325  if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
326 /*
327  // We must decide here if we can dare declare one of our faces
328  // as having a "valid" normal (i.e. allBehind = true). This
329  // is never possible if the face faces "inward" in r *unless*
330  // we have only one side
331  //
332  G4bool allBehind;
333  if ((corner->z > next->z) && (numSide > 1))
334  {
335  allBehind = false;
336  }
337  else
338  {
339  //
340  // Otherwise, it is only true if the line passing
341  // through the two points of the segment do not
342  // split the r/z cross section
343  //
344  allBehind = !rz->BisectedBy( corner->r, corner->z,
345  next->r, next->z, kCarTolerance );
346  }
347 */
348  *face++ = new G4PolyhedraSide( prev, corner, next, nextNext,
350  } while( prev=corner, corner=next, corner > corners );
351 
352  if (phiIsOpen)
353  {
354  //
355  // Construct phi open edges
356  //
357  *face++ = new G4PolyPhiFace( rz, startPhi, phiTotal/numSide, endPhi );
358  *face++ = new G4PolyPhiFace( rz, endPhi, phiTotal/numSide, startPhi );
359  }
360 
361  //
362  // We might have dropped a face or two: recalculate numFace
363  //
364  numFace = face-faces;
365 
366  //
367  // Make enclosingCylinder
368  //
370  new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
371 }
372 
373 
374 //
375 // Fake default constructor - sets only member data and allocates memory
376 // for usage restricted to object persistency.
377 //
379  : G4VCSGfaceted(a), numSide(0), startPhi(0.), endPhi(0.),
380  phiIsOpen(false), genericPgon(false), numCorner(0), corners(0),
381  original_parameters(0), enclosingCylinder(0)
382 {
383 }
384 
385 
386 //
387 // Destructor
388 //
390 {
391  delete [] corners;
393 
394  delete enclosingCylinder;
395 }
396 
397 
398 //
399 // Copy constructor
400 //
402  : G4VCSGfaceted( source )
403 {
404  CopyStuff( source );
405 }
406 
407 
408 //
409 // Assignment operator
410 //
412 {
413  if (this == &source) return *this;
414 
415  G4VCSGfaceted::operator=( source );
416 
417  delete [] corners;
419 
420  delete enclosingCylinder;
421 
422  CopyStuff( source );
423 
424  return *this;
425 }
426 
427 
428 //
429 // CopyStuff
430 //
431 void G4Polyhedra::CopyStuff( const G4Polyhedra &source )
432 {
433  //
434  // Simple stuff
435  //
436  numSide = source.numSide;
437  startPhi = source.startPhi;
438  endPhi = source.endPhi;
439  phiIsOpen = source.phiIsOpen;
440  numCorner = source.numCorner;
441  genericPgon= source.genericPgon;
442 
443  //
444  // The corner array
445  //
447 
448  G4PolyhedraSideRZ *corn = corners,
449  *sourceCorn = source.corners;
450  do // Loop checking, 13.08.2015, G.Cosmo
451  {
452  *corn = *sourceCorn;
453  } while( ++sourceCorn, ++corn < corners+numCorner );
454 
455  //
456  // Original parameters
457  //
458  if (source.original_parameters)
459  {
462  }
463 
464  //
465  // Enclosing cylinder
466  //
468 
469  fRebuildPolyhedron = false;
470  fpPolyhedron = 0;
471 }
472 
473 
474 //
475 // Reset
476 //
477 // Recalculates and reshapes the solid, given pre-assigned scaled
478 // original_parameters.
479 //
481 {
482  if (genericPgon)
483  {
484  std::ostringstream message;
485  message << "Solid " << GetName() << " built using generic construct."
486  << G4endl << "Not applicable to the generic construct !";
487  G4Exception("G4Polyhedra::Reset()", "GeomSolids1001",
488  JustWarning, message, "Parameters NOT resetted.");
489  return 1;
490  }
491 
492  //
493  // Clear old setup
494  //
496  delete [] corners;
497  delete enclosingCylinder;
498 
499  //
500  // Rebuild polyhedra
501  //
502  G4ReduciblePolygon *rz =
510  delete rz;
511 
512  return 0;
513 }
514 
515 
516 //
517 // Inside
518 //
519 // This is an override of G4VCSGfaceted::Inside, created in order
520 // to speed things up by first checking with G4EnclosingCylinder.
521 //
523 {
524  //
525  // Quick test
526  //
527  if (enclosingCylinder->MustBeOutside(p)) return kOutside;
528 
529  //
530  // Long answer
531  //
532  return G4VCSGfaceted::Inside(p);
533 }
534 
535 
536 //
537 // DistanceToIn
538 //
539 // This is an override of G4VCSGfaceted::Inside, created in order
540 // to speed things up by first checking with G4EnclosingCylinder.
541 //
543  const G4ThreeVector &v ) const
544 {
545  //
546  // Quick test
547  //
548  if (enclosingCylinder->ShouldMiss(p,v))
549  return kInfinity;
550 
551  //
552  // Long answer
553  //
554  return G4VCSGfaceted::DistanceToIn( p, v );
555 }
556 
557 
558 //
559 // DistanceToIn
560 //
562 {
563  return G4VCSGfaceted::DistanceToIn(p);
564 }
565 
566 
567 //
568 // ComputeDimensions
569 //
571  const G4int n,
572  const G4VPhysicalVolume* pRep )
573 {
574  p->ComputeDimensions(*this,n,pRep);
575 }
576 
577 
578 //
579 // GetEntityType
580 //
582 {
583  return G4String("G4Polyhedra");
584 }
585 
586 
587 //
588 // Make a clone of the object
589 //
591 {
592  return new G4Polyhedra(*this);
593 }
594 
595 
596 //
597 // Stream object contents to an output stream
598 //
599 std::ostream& G4Polyhedra::StreamInfo( std::ostream& os ) const
600 {
601  G4int oldprc = os.precision(16);
602  os << "-----------------------------------------------------------\n"
603  << " *** Dump for solid - " << GetName() << " ***\n"
604  << " ===================================================\n"
605  << " Solid type: G4Polyhedra\n"
606  << " Parameters: \n"
607  << " starting phi angle : " << startPhi/degree << " degrees \n"
608  << " ending phi angle : " << endPhi/degree << " degrees \n"
609  << " number of sides : " << numSide << " \n";
610  G4int i=0;
611  if (!genericPgon)
612  {
614  os << " number of Z planes: " << numPlanes << "\n"
615  << " Z values: \n";
616  for (i=0; i<numPlanes; i++)
617  {
618  os << " Z plane " << i << ": "
619  << original_parameters->Z_values[i] << "\n";
620  }
621  os << " Tangent distances to inner surface (Rmin): \n";
622  for (i=0; i<numPlanes; i++)
623  {
624  os << " Z plane " << i << ": "
625  << original_parameters->Rmin[i] << "\n";
626  }
627  os << " Tangent distances to outer surface (Rmax): \n";
628  for (i=0; i<numPlanes; i++)
629  {
630  os << " Z plane " << i << ": "
631  << original_parameters->Rmax[i] << "\n";
632  }
633  }
634  os << " number of RZ points: " << numCorner << "\n"
635  << " RZ values (corners): \n";
636  for (i=0; i<numCorner; i++)
637  {
638  os << " "
639  << corners[i].r << ", " << corners[i].z << "\n";
640  }
641  os << "-----------------------------------------------------------\n";
642  os.precision(oldprc);
643 
644  return os;
645 }
646 
647 
648 //
649 // GetPointOnPlane
650 //
651 // Auxiliary method for get point on surface
652 //
654  G4ThreeVector p2, G4ThreeVector p3) const
655 {
656  G4double lambda1, lambda2, chose,aOne,aTwo;
657  G4ThreeVector t, u, v, w, Area, normal;
658  aOne = 1.;
659  aTwo = 1.;
660 
661  t = p1 - p0;
662  u = p2 - p1;
663  v = p3 - p2;
664  w = p0 - p3;
665 
666  chose = RandFlat::shoot(0.,aOne+aTwo);
667  if( (chose>=0.) && (chose < aOne) )
668  {
669  lambda1 = RandFlat::shoot(0.,1.);
670  lambda2 = RandFlat::shoot(0.,lambda1);
671  return (p2+lambda1*v+lambda2*w);
672  }
673 
674  lambda1 = RandFlat::shoot(0.,1.);
675  lambda2 = RandFlat::shoot(0.,lambda1);
676  return (p0+lambda1*t+lambda2*u);
677 }
678 
679 
680 //
681 // GetPointOnTriangle
682 //
683 // Auxiliary method for get point on surface
684 //
687  G4ThreeVector p3) const
688 {
689  G4double lambda1,lambda2;
690  G4ThreeVector v=p3-p1, w=p1-p2;
691 
692  lambda1 = RandFlat::shoot(0.,1.);
693  lambda2 = RandFlat::shoot(0.,lambda1);
694 
695  return (p2 + lambda1*w + lambda2*v);
696 }
697 
698 
699 //
700 // GetPointOnSurface
701 //
703 {
704  if( !genericPgon ) // Polyhedra by faces
705  {
706  G4int j, numPlanes = original_parameters->Num_z_planes, Flag=0;
707  G4double chose, totArea=0., Achose1, Achose2,
708  rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
709  G4double a, b, l2, rang, totalPhi, ksi,
710  area, aTop=0., aBottom=0., zVal=0.;
711 
712  G4ThreeVector p0, p1, p2, p3;
713  std::vector<G4double> aVector1;
714  std::vector<G4double> aVector2;
715  std::vector<G4double> aVector3;
716 
717  totalPhi= (phiIsOpen) ? (endPhi-startPhi) : twopi;
718  ksi = totalPhi/numSide;
719  G4double cosksi = std::cos(ksi/2.);
720 
721  // Below we generate the areas relevant to our solid
722  //
723  for(j=0; j<numPlanes-1; j++)
724  {
725  a = original_parameters->Rmax[j+1];
726  b = original_parameters->Rmax[j];
728  -original_parameters->Z_values[j+1]) + sqr(b-a);
729  area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
730  aVector1.push_back(area);
731  }
732 
733  for(j=0; j<numPlanes-1; j++)
734  {
735  a = original_parameters->Rmin[j+1];//*cosksi;
736  b = original_parameters->Rmin[j];//*cosksi;
738  -original_parameters->Z_values[j+1]) + sqr(b-a);
739  area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
740  aVector2.push_back(area);
741  }
742 
743  for(j=0; j<numPlanes-1; j++)
744  {
745  if(phiIsOpen == true)
746  {
747  aVector3.push_back(0.5*(original_parameters->Rmax[j]
750  -original_parameters->Rmin[j+1])
751  *std::fabs(original_parameters->Z_values[j+1]
753  }
754  else { aVector3.push_back(0.); }
755  }
756 
757  for(j=0; j<numPlanes-1; j++)
758  {
759  totArea += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
760  }
761 
762  // Must include top and bottom areas
763  //
764  if(original_parameters->Rmax[numPlanes-1] != 0.)
765  {
766  a = original_parameters->Rmax[numPlanes-1];
767  b = original_parameters->Rmin[numPlanes-1];
768  l2 = sqr(a-b);
769  aTop = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
770  }
771 
772  if(original_parameters->Rmax[0] != 0.)
773  {
774  a = original_parameters->Rmax[0];
775  b = original_parameters->Rmin[0];
776  l2 = sqr(a-b);
777  aBottom = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
778  }
779 
780  Achose1 = 0.;
781  Achose2 = numSide*(aVector1[0]+aVector2[0])+2.*aVector3[0];
782 
783  chose = RandFlat::shoot(0.,totArea+aTop+aBottom);
784  if( (chose >= 0.) && (chose < aTop + aBottom) )
785  {
786  chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
787  rang = std::floor((chose-startPhi)/ksi-0.01);
788  if(rang<0) { rang=0; }
789  rang = std::fabs(rang);
790  sinphi1 = std::sin(startPhi+rang*ksi);
791  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
792  cosphi1 = std::cos(startPhi+rang*ksi);
793  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
794  chose = RandFlat::shoot(0., aTop + aBottom);
795  if(chose>=0. && chose<aTop)
796  {
797  rad1 = original_parameters->Rmin[numPlanes-1];
798  rad2 = original_parameters->Rmax[numPlanes-1];
799  zVal = original_parameters->Z_values[numPlanes-1];
800  }
801  else
802  {
803  rad1 = original_parameters->Rmin[0];
804  rad2 = original_parameters->Rmax[0];
805  zVal = original_parameters->Z_values[0];
806  }
807  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
808  p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
809  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
810  p3 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
811  return GetPointOnPlane(p0,p1,p2,p3);
812  }
813  else
814  {
815  for (j=0; j<numPlanes-1; j++)
816  {
817  if( ((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes-2) )
818  {
819  Flag = j; break;
820  }
821  Achose1 += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
822  Achose2 = Achose1 + numSide*(aVector1[j+1]+aVector2[j+1])
823  + 2.*aVector3[j+1];
824  }
825  }
826 
827  // At this point we have chosen a subsection
828  // between to adjacent plane cuts...
829 
830  j = Flag;
831 
832  totArea = numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
833  chose = RandFlat::shoot(0.,totArea);
834 
835  if( (chose>=0.) && (chose<numSide*aVector1[j]) )
836  {
837  chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
838  rang = std::floor((chose-startPhi)/ksi-0.01);
839  if(rang<0) { rang=0; }
840  rang = std::fabs(rang);
841  rad1 = original_parameters->Rmax[j];
842  rad2 = original_parameters->Rmax[j+1];
843  sinphi1 = std::sin(startPhi+rang*ksi);
844  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
845  cosphi1 = std::cos(startPhi+rang*ksi);
846  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
847  zVal = original_parameters->Z_values[j];
848 
849  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
850  p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
851 
852  zVal = original_parameters->Z_values[j+1];
853 
854  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
855  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
856  return GetPointOnPlane(p0,p1,p2,p3);
857  }
858  else if ( (chose >= numSide*aVector1[j])
859  && (chose <= numSide*(aVector1[j]+aVector2[j])) )
860  {
861  chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
862  rang = std::floor((chose-startPhi)/ksi-0.01);
863  if(rang<0) { rang=0; }
864  rang = std::fabs(rang);
865  rad1 = original_parameters->Rmin[j];
866  rad2 = original_parameters->Rmin[j+1];
867  sinphi1 = std::sin(startPhi+rang*ksi);
868  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
869  cosphi1 = std::cos(startPhi+rang*ksi);
870  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
871  zVal = original_parameters->Z_values[j];
872 
873  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
874  p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
875 
876  zVal = original_parameters->Z_values[j+1];
877 
878  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
879  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
880  return GetPointOnPlane(p0,p1,p2,p3);
881  }
882 
883  chose = RandFlat::shoot(0.,2.2);
884  if( (chose>=0.) && (chose < 1.) )
885  {
886  rang = startPhi;
887  }
888  else
889  {
890  rang = endPhi;
891  }
892 
893  cosphi1 = std::cos(rang); rad1 = original_parameters->Rmin[j];
894  sinphi1 = std::sin(rang); rad2 = original_parameters->Rmax[j];
895 
896  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
898  p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
900 
901  rad1 = original_parameters->Rmax[j+1];
902  rad2 = original_parameters->Rmin[j+1];
903 
904  p2 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
906  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
908  return GetPointOnPlane(p0,p1,p2,p3);
909  }
910  else // Generic polyhedra
911  {
912  return GetPointOnSurfaceGeneric();
913  }
914 }
915 
916 //
917 // CreatePolyhedron
918 //
920 {
921  if (!genericPgon)
922  {
930  }
931  else
932  {
933  // The following code prepares for:
934  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
935  // const double xyz[][3],
936  // const int faces_vec[][4])
937  // Here is an extract from the header file HepPolyhedron.h:
955  G4int nNodes;
956  G4int nFaces;
957  typedef G4double double3[3];
958  double3* xyz;
959  typedef G4int int4[4];
960  int4* faces_vec;
961  if (phiIsOpen)
962  {
963  // Triangulate open ends. Simple ear-chopping algorithm...
964  // I'm not sure how robust this algorithm is (J.Allison).
965  //
966  std::vector<G4bool> chopped(numCorner, false);
967  std::vector<G4int*> triQuads;
968  G4int remaining = numCorner;
969  G4int iStarter = 0;
970  while (remaining >= 3) // Loop checking, 13.08.2015, G.Cosmo
971  {
972  // Find unchopped corners...
973  //
974  G4int A = -1, B = -1, C = -1;
975  G4int iStepper = iStarter;
976  do // Loop checking, 13.08.2015, G.Cosmo
977  {
978  if (A < 0) { A = iStepper; }
979  else if (B < 0) { B = iStepper; }
980  else if (C < 0) { C = iStepper; }
981  do // Loop checking, 13.08.2015, G.Cosmo
982  {
983  if (++iStepper >= numCorner) iStepper = 0;
984  }
985  while (chopped[iStepper]);
986  }
987  while (C < 0 && iStepper != iStarter);
988 
989  // Check triangle at B is pointing outward (an "ear").
990  // Sign of z cross product determines...
991 
992  G4double BAr = corners[A].r - corners[B].r;
993  G4double BAz = corners[A].z - corners[B].z;
994  G4double BCr = corners[C].r - corners[B].r;
995  G4double BCz = corners[C].z - corners[B].z;
996  if (BAr * BCz - BAz * BCr < kCarTolerance)
997  {
998  G4int* tq = new G4int[3];
999  tq[0] = A + 1;
1000  tq[1] = B + 1;
1001  tq[2] = C + 1;
1002  triQuads.push_back(tq);
1003  chopped[B] = true;
1004  --remaining;
1005  }
1006  else
1007  {
1008  do // Loop checking, 13.08.2015, G.Cosmo
1009  {
1010  if (++iStarter >= numCorner) { iStarter = 0; }
1011  }
1012  while (chopped[iStarter]);
1013  }
1014  }
1015 
1016  // Transfer to faces...
1017 
1018  nNodes = (numSide + 1) * numCorner;
1019  nFaces = numSide * numCorner + 2 * triQuads.size();
1020  faces_vec = new int4[nFaces];
1021  G4int iface = 0;
1022  G4int addition = numCorner * numSide;
1023  G4int d = numCorner - 1;
1024  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
1025  {
1026  for (size_t i = 0; i < triQuads.size(); ++i)
1027  {
1028  // Negative for soft/auxiliary/normally invisible edges...
1029  //
1030  G4int a, b, c;
1031  if (iEnd == 0)
1032  {
1033  a = triQuads[i][0];
1034  b = triQuads[i][1];
1035  c = triQuads[i][2];
1036  }
1037  else
1038  {
1039  a = triQuads[i][0] + addition;
1040  b = triQuads[i][2] + addition;
1041  c = triQuads[i][1] + addition;
1042  }
1043  G4int ab = std::abs(b - a);
1044  G4int bc = std::abs(c - b);
1045  G4int ca = std::abs(a - c);
1046  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
1047  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
1048  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
1049  faces_vec[iface][3] = 0;
1050  ++iface;
1051  }
1052  }
1053 
1054  // Continue with sides...
1055 
1056  xyz = new double3[nNodes];
1057  const G4double dPhi = (endPhi - startPhi) / numSide;
1058  G4double phi = startPhi;
1059  G4int ixyz = 0;
1060  for (G4int iSide = 0; iSide < numSide; ++iSide)
1061  {
1062  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1063  {
1064  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1065  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1066  xyz[ixyz][2] = corners[iCorner].z;
1067  if (iCorner < numCorner - 1)
1068  {
1069  faces_vec[iface][0] = ixyz + 1;
1070  faces_vec[iface][1] = ixyz + numCorner + 1;
1071  faces_vec[iface][2] = ixyz + numCorner + 2;
1072  faces_vec[iface][3] = ixyz + 2;
1073  }
1074  else
1075  {
1076  faces_vec[iface][0] = ixyz + 1;
1077  faces_vec[iface][1] = ixyz + numCorner + 1;
1078  faces_vec[iface][2] = ixyz + 2;
1079  faces_vec[iface][3] = ixyz - numCorner + 2;
1080  }
1081  ++iface;
1082  ++ixyz;
1083  }
1084  phi += dPhi;
1085  }
1086 
1087  // Last corners...
1088 
1089  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1090  {
1091  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1092  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1093  xyz[ixyz][2] = corners[iCorner].z;
1094  ++ixyz;
1095  }
1096  }
1097  else // !phiIsOpen - i.e., a complete 360 degrees.
1098  {
1099  nNodes = numSide * numCorner;
1100  nFaces = numSide * numCorner;;
1101  xyz = new double3[nNodes];
1102  faces_vec = new int4[nFaces];
1103  // const G4double dPhi = (endPhi - startPhi) / numSide;
1104  const G4double dPhi = twopi / numSide; // !phiIsOpen endPhi-startPhi = 360 degrees.
1105  G4double phi = startPhi;
1106  G4int ixyz = 0, iface = 0;
1107  for (G4int iSide = 0; iSide < numSide; ++iSide)
1108  {
1109  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1110  {
1111  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1112  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1113  xyz[ixyz][2] = corners[iCorner].z;
1114  if (iSide < numSide - 1)
1115  {
1116  if (iCorner < numCorner - 1)
1117  {
1118  faces_vec[iface][0] = ixyz + 1;
1119  faces_vec[iface][1] = ixyz + numCorner + 1;
1120  faces_vec[iface][2] = ixyz + numCorner + 2;
1121  faces_vec[iface][3] = ixyz + 2;
1122  }
1123  else
1124  {
1125  faces_vec[iface][0] = ixyz + 1;
1126  faces_vec[iface][1] = ixyz + numCorner + 1;
1127  faces_vec[iface][2] = ixyz + 2;
1128  faces_vec[iface][3] = ixyz - numCorner + 2;
1129  }
1130  }
1131  else // Last side joins ends...
1132  {
1133  if (iCorner < numCorner - 1)
1134  {
1135  faces_vec[iface][0] = ixyz + 1;
1136  faces_vec[iface][1] = ixyz + numCorner - nFaces + 1;
1137  faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
1138  faces_vec[iface][3] = ixyz + 2;
1139  }
1140  else
1141  {
1142  faces_vec[iface][0] = ixyz + 1;
1143  faces_vec[iface][1] = ixyz - nFaces + numCorner + 1;
1144  faces_vec[iface][2] = ixyz - nFaces + 2;
1145  faces_vec[iface][3] = ixyz - numCorner + 2;
1146  }
1147  }
1148  ++ixyz;
1149  ++iface;
1150  }
1151  phi += dPhi;
1152  }
1153  }
1154  G4Polyhedron* polyhedron = new G4Polyhedron;
1155  G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
1156  delete [] faces_vec;
1157  delete [] xyz;
1158  if (problem)
1159  {
1160  std::ostringstream message;
1161  message << "Problem creating G4Polyhedron for: " << GetName();
1162  G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
1163  JustWarning, message);
1164  delete polyhedron;
1165  return 0;
1166  }
1167  else
1168  {
1169  return polyhedron;
1170  }
1171  }
1172 }
1173 
1174 
1176 {
1177  G4int numPlanes = (G4int)numCorner;
1178  G4bool isConvertible=true;
1179  G4double Zmax=rz->Bmax();
1180  rz->StartWithZMin();
1181 
1182  // Prepare vectors for storage
1183  //
1184  std::vector<G4double> Z;
1185  std::vector<G4double> Rmin;
1186  std::vector<G4double> Rmax;
1187 
1188  G4int countPlanes=1;
1189  G4int icurr=0;
1190  G4int icurl=0;
1191 
1192  // first plane Z=Z[0]
1193  //
1194  Z.push_back(corners[0].z);
1195  G4double Zprev=Z[0];
1196  if (Zprev == corners[1].z)
1197  {
1198  Rmin.push_back(corners[0].r);
1199  Rmax.push_back (corners[1].r);icurr=1;
1200  }
1201  else if (Zprev == corners[numPlanes-1].z)
1202  {
1203  Rmin.push_back(corners[numPlanes-1].r);
1204  Rmax.push_back (corners[0].r);
1205  icurl=numPlanes-1;
1206  }
1207  else
1208  {
1209  Rmin.push_back(corners[0].r);
1210  Rmax.push_back (corners[0].r);
1211  }
1212 
1213  // next planes until last
1214  //
1215  G4int inextr=0, inextl=0;
1216  for (G4int i=0; i < numPlanes-2; i++)
1217  {
1218  inextr=1+icurr;
1219  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1220 
1221  if((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax)) { break; }
1222 
1223  G4double Zleft = corners[inextl].z;
1224  G4double Zright = corners[inextr].z;
1225  if(Zright>Zleft)
1226  {
1227  Z.push_back(Zleft);
1228  countPlanes++;
1229  G4double difZr=corners[inextr].z - corners[icurr].z;
1230  G4double difZl=corners[inextl].z - corners[icurl].z;
1231 
1232  if(std::fabs(difZl) < kCarTolerance)
1233  {
1234  if(std::fabs(difZr) < kCarTolerance)
1235  {
1236  Rmin.push_back(corners[inextl].r);
1237  Rmax.push_back(corners[icurr].r);
1238  }
1239  else
1240  {
1241  Rmin.push_back(corners[inextl].r);
1242  Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1243  *(corners[inextr].r - corners[icurr].r));
1244  }
1245  }
1246  else if (difZl >= kCarTolerance)
1247  {
1248  if(std::fabs(difZr) < kCarTolerance)
1249  {
1250  Rmin.push_back(corners[icurl].r);
1251  Rmax.push_back(corners[icurr].r);
1252  }
1253  else
1254  {
1255  Rmin.push_back(corners[icurl].r);
1256  Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1257  *(corners[inextr].r - corners[icurr].r));
1258  }
1259  }
1260  else
1261  {
1262  isConvertible=false; break;
1263  }
1264  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1265  }
1266  else if(std::fabs(Zright-Zleft)<kCarTolerance) // Zright=Zleft
1267  {
1268  Z.push_back(Zleft);
1269  countPlanes++;
1270  icurr++;
1271 
1272  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1273 
1274  Rmin.push_back(corners[inextl].r);
1275  Rmax.push_back (corners[inextr].r);
1276  }
1277  else // Zright<Zleft
1278  {
1279  Z.push_back(Zright);
1280  countPlanes++;
1281 
1282  G4double difZr=corners[inextr].z - corners[icurr].z;
1283  G4double difZl=corners[inextl].z - corners[icurl].z;
1284  if(std::fabs(difZr) < kCarTolerance)
1285  {
1286  if(std::fabs(difZl) < kCarTolerance)
1287  {
1288  Rmax.push_back(corners[inextr].r);
1289  Rmin.push_back(corners[icurr].r);
1290  }
1291  else
1292  {
1293  Rmin.push_back(corners[icurl].r + (Zright-corners[icurl].z)/difZl
1294  * (corners[inextl].r - corners[icurl].r));
1295  Rmax.push_back(corners[inextr].r);
1296  }
1297  icurr++;
1298  } // plate
1299  else if (difZr >= kCarTolerance)
1300  {
1301  if(std::fabs(difZl) < kCarTolerance)
1302  {
1303  Rmax.push_back(corners[inextr].r);
1304  Rmin.push_back (corners[icurr].r);
1305  }
1306  else
1307  {
1308  Rmax.push_back(corners[inextr].r);
1309  Rmin.push_back (corners[icurl].r+(Zright-corners[icurl].z)/difZl
1310  * (corners[inextl].r - corners[icurl].r));
1311  }
1312  icurr++;
1313  }
1314  else
1315  {
1316  isConvertible=false; break;
1317  }
1318  }
1319  } // end for loop
1320 
1321  // last plane Z=Zmax
1322  //
1323  Z.push_back(Zmax);
1324  countPlanes++;
1325  inextr=1+icurr;
1326  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1327 
1328  Rmax.push_back(corners[inextr].r);
1329  Rmin.push_back(corners[inextl].r);
1330 
1331  // Set original parameters Rmin,Rmax,Z
1332  //
1333  if(isConvertible)
1334  {
1337  original_parameters->Z_values = new G4double[countPlanes];
1338  original_parameters->Rmin = new G4double[countPlanes];
1339  original_parameters->Rmax = new G4double[countPlanes];
1340 
1341  for(G4int j=0; j < countPlanes; j++)
1342  {
1343  original_parameters->Z_values[j] = Z[j];
1344  original_parameters->Rmax[j] = Rmax[j];
1345  original_parameters->Rmin[j] = Rmin[j];
1346  }
1349  original_parameters->Num_z_planes = countPlanes;
1350 
1351  }
1352  else // Set parameters(r,z) with Rmin==0 as convention
1353  {
1354 #ifdef G4SPECSDEBUG
1355  std::ostringstream message;
1356  message << "Polyhedra " << GetName() << G4endl
1357  << "cannot be converted to Polyhedra with (Rmin,Rmaz,Z) parameters!";
1358  G4Exception("G4Polyhedra::SetOriginalParameters()",
1359  "GeomSolids0002", JustWarning, message);
1360 #endif
1363  original_parameters->Z_values = new G4double[numPlanes];
1364  original_parameters->Rmin = new G4double[numPlanes];
1365  original_parameters->Rmax = new G4double[numPlanes];
1366 
1367  for(G4int j=0; j < numPlanes; j++)
1368  {
1370  original_parameters->Rmax[j] = corners[j].r;
1371  original_parameters->Rmin[j] = 0.0;
1372  }
1375  original_parameters->Num_z_planes = numPlanes;
1376  }
1377  //return isConvertible;
1378 }
1379 
1380 #endif
void CopyStuff(const G4Polyhedra &source)
Definition: G4Polyhedra.cc:431
G4String GetName() const
G4bool CrossesItself(G4double tolerance)
ThreeVector shoot(const G4int Ap, const G4int Af)
G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
virtual ~G4Polyhedra()
Definition: G4Polyhedra.cc:389
G4bool fRebuildPolyhedron
G4Polyhedron * fpPolyhedron
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
G4GeometryType GetEntityType() const
Definition: G4Polyhedra.cc:581
G4double Amin() const
G4double z
Definition: TRTMaterials.hh:39
G4PolyhedraHistorical * original_parameters
Definition: G4Polyhedra.hh:187
G4String name
Definition: TRTMaterials.hh:40
G4ThreeVector GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3) const
Definition: G4Polyhedra.cc:653
G4bool genericPgon
Definition: G4Polyhedra.hh:184
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4PolyhedraSideRZ * corners
Definition: G4Polyhedra.hh:186
G4double a
Definition: TRTMaterials.hh:39
G4bool MustBeOutside(const G4ThreeVector &p) const
G4bool Reset()
Definition: G4Polyhedra.cc:480
int G4int
Definition: G4Types.hh:78
void Create(G4double phiStart, G4double phiTotal, G4int numSide, G4ReduciblePolygon *rz)
Definition: G4Polyhedra.cc:201
G4ThreeVector GetPointOnSurface() const
Definition: G4Polyhedra.cc:702
void DumpInfo() const
G4double endPhi
Definition: G4Polyhedra.hh:182
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4bool RemoveDuplicateVertices(G4double tolerance)
G4EnclosingCylinder * enclosingCylinder
Definition: G4Polyhedra.hh:189
G4bool RemoveRedundantVertices(G4double tolerance)
G4ThreeVector GetPointOnSurfaceGeneric() const
G4int numCorner
Definition: G4Polyhedra.hh:185
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Polyhedra.cc:570
G4int NumVertices() const
G4double Bmax() const
void ScaleA(G4double scale)
bool G4bool
Definition: G4Types.hh:79
G4double startPhi
Definition: G4Polyhedra.hh:181
const G4double p2
const G4double p1
#define DBL_EPSILON
Definition: templates.hh:87
const G4int n
static const G4double A[nN]
G4VCSGface ** faces
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Polyhedra.cc:542
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
EInside Inside(const G4ThreeVector &p) const
Definition: G4Polyhedra.cc:522
static const G4double ab
const G4double p0
G4Polyhedron * CreatePolyhedron() const
Definition: G4Polyhedra.cc:919
EInside
Definition: geomdefs.hh:58
static const double degree
Definition: G4SIunits.hh:125
G4Polyhedra & operator=(const G4Polyhedra &source)
Definition: G4Polyhedra.cc:411
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:305
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
void SetOriginalParameters(G4PolyhedraHistorical *pars)
G4ThreeVector GetPointOnTriangle(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2) const
Definition: G4Polyhedra.cc:685
virtual EInside Inside(const G4ThreeVector &p) const
G4bool phiIsOpen
Definition: G4Polyhedra.hh:183
G4bool ShouldMiss(const G4ThreeVector &p, const G4ThreeVector &v) const
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Polyhedra.cc:599
G4VSolid * Clone() const
Definition: G4Polyhedra.cc:590
G4Polyhedra(const G4String &name, G4double phiStart, G4double phiTotal, G4int numSide, G4int numZPlanes, const G4double zPlane[], const G4double rInner[], const G4double rOuter[])
Definition: G4Polyhedra.cc:79