Geant4  10.00.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 84624 2014-10-17 09:56:00Z 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 ) startPhi += twopi;
255  //
256  // Phi opening? Account for some possible roundoff, and interpret
257  // nonsense value as representing no phi opening
258  //
259  if ( (phiTotal <= 0) || (phiTotal > twopi*(1-DBL_EPSILON)) )
260  {
261  phiIsOpen = false;
262  endPhi = phiStart+twopi;
263  }
264  else
265  {
266  phiIsOpen = true;
267 
268  //
269  // Convert phi into our convention
270  //
271  endPhi = phiStart+phiTotal;
272  while( endPhi < startPhi ) endPhi += twopi;
273  }
274 
275  //
276  // Save number sides
277  //
278  numSide = theNumSide;
279 
280  //
281  // Allocate corner array.
282  //
284 
285  //
286  // Copy corners
287  //
288  G4ReduciblePolygonIterator iterRZ(rz);
289 
290  G4PolyhedraSideRZ *next = corners;
291  iterRZ.Begin();
292  do
293  {
294  next->r = iterRZ.GetA();
295  next->z = iterRZ.GetB();
296  } while( ++next, iterRZ.Next() );
297 
298  //
299  // Allocate face pointer array
300  //
302  faces = new G4VCSGface*[numFace];
303 
304  //
305  // Construct side faces
306  //
307  // To do so properly, we need to keep track of four successive RZ
308  // corners.
309  //
310  // But! Don't construct a face if both points are at zero radius!
311  //
312  G4PolyhedraSideRZ *corner = corners,
313  *prev = corners + numCorner-1,
314  *nextNext;
315  G4VCSGface **face = faces;
316  do
317  {
318  next = corner+1;
319  if (next >= corners+numCorner) next = corners;
320  nextNext = next+1;
321  if (nextNext >= corners+numCorner) nextNext = corners;
322 
323  if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
324 /*
325  // We must decide here if we can dare declare one of our faces
326  // as having a "valid" normal (i.e. allBehind = true). This
327  // is never possible if the face faces "inward" in r *unless*
328  // we have only one side
329  //
330  G4bool allBehind;
331  if ((corner->z > next->z) && (numSide > 1))
332  {
333  allBehind = false;
334  }
335  else
336  {
337  //
338  // Otherwise, it is only true if the line passing
339  // through the two points of the segment do not
340  // split the r/z cross section
341  //
342  allBehind = !rz->BisectedBy( corner->r, corner->z,
343  next->r, next->z, kCarTolerance );
344  }
345 */
346  *face++ = new G4PolyhedraSide( prev, corner, next, nextNext,
348  } while( prev=corner, corner=next, corner > corners );
349 
350  if (phiIsOpen)
351  {
352  //
353  // Construct phi open edges
354  //
355  *face++ = new G4PolyPhiFace( rz, startPhi, phiTotal/numSide, endPhi );
356  *face++ = new G4PolyPhiFace( rz, endPhi, phiTotal/numSide, startPhi );
357  }
358 
359  //
360  // We might have dropped a face or two: recalculate numFace
361  //
362  numFace = face-faces;
363 
364  //
365  // Make enclosingCylinder
366  //
368  new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
369 }
370 
371 
372 //
373 // Fake default constructor - sets only member data and allocates memory
374 // for usage restricted to object persistency.
375 //
377  : G4VCSGfaceted(a), numSide(0), startPhi(0.), endPhi(0.),
378  phiIsOpen(false), genericPgon(false), numCorner(0), corners(0),
379  original_parameters(0), enclosingCylinder(0)
380 {
381 }
382 
383 
384 //
385 // Destructor
386 //
388 {
389  delete [] corners;
391 
392  delete enclosingCylinder;
393 }
394 
395 
396 //
397 // Copy constructor
398 //
400  : G4VCSGfaceted( source )
401 {
402  CopyStuff( source );
403 }
404 
405 
406 //
407 // Assignment operator
408 //
410 {
411  if (this == &source) return *this;
412 
413  G4VCSGfaceted::operator=( source );
414 
415  delete [] corners;
417 
418  delete enclosingCylinder;
419 
420  CopyStuff( source );
421 
422  return *this;
423 }
424 
425 
426 //
427 // CopyStuff
428 //
429 void G4Polyhedra::CopyStuff( const G4Polyhedra &source )
430 {
431  //
432  // Simple stuff
433  //
434  numSide = source.numSide;
435  startPhi = source.startPhi;
436  endPhi = source.endPhi;
437  phiIsOpen = source.phiIsOpen;
438  numCorner = source.numCorner;
439  genericPgon= source.genericPgon;
440 
441  //
442  // The corner array
443  //
445 
446  G4PolyhedraSideRZ *corn = corners,
447  *sourceCorn = source.corners;
448  do
449  {
450  *corn = *sourceCorn;
451  } while( ++sourceCorn, ++corn < corners+numCorner );
452 
453  //
454  // Original parameters
455  //
456  if (source.original_parameters)
457  {
460  }
461 
462  //
463  // Enclosing cylinder
464  //
466 
467  fRebuildPolyhedron = false;
468  fpPolyhedron = 0;
469 }
470 
471 
472 //
473 // Reset
474 //
475 // Recalculates and reshapes the solid, given pre-assigned scaled
476 // original_parameters.
477 //
479 {
480  if (genericPgon)
481  {
482  std::ostringstream message;
483  message << "Solid " << GetName() << " built using generic construct."
484  << G4endl << "Not applicable to the generic construct !";
485  G4Exception("G4Polyhedra::Reset()", "GeomSolids1001",
486  JustWarning, message, "Parameters NOT resetted.");
487  return 1;
488  }
489 
490  //
491  // Clear old setup
492  //
494  delete [] corners;
495  delete enclosingCylinder;
496 
497  //
498  // Rebuild polyhedra
499  //
500  G4ReduciblePolygon *rz =
508  delete rz;
509 
510  return 0;
511 }
512 
513 
514 //
515 // Inside
516 //
517 // This is an override of G4VCSGfaceted::Inside, created in order
518 // to speed things up by first checking with G4EnclosingCylinder.
519 //
521 {
522  //
523  // Quick test
524  //
525  if (enclosingCylinder->MustBeOutside(p)) return kOutside;
526 
527  //
528  // Long answer
529  //
530  return G4VCSGfaceted::Inside(p);
531 }
532 
533 
534 //
535 // DistanceToIn
536 //
537 // This is an override of G4VCSGfaceted::Inside, created in order
538 // to speed things up by first checking with G4EnclosingCylinder.
539 //
541  const G4ThreeVector &v ) const
542 {
543  //
544  // Quick test
545  //
546  if (enclosingCylinder->ShouldMiss(p,v))
547  return kInfinity;
548 
549  //
550  // Long answer
551  //
552  return G4VCSGfaceted::DistanceToIn( p, v );
553 }
554 
555 
556 //
557 // DistanceToIn
558 //
560 {
561  return G4VCSGfaceted::DistanceToIn(p);
562 }
563 
564 
565 //
566 // ComputeDimensions
567 //
569  const G4int n,
570  const G4VPhysicalVolume* pRep )
571 {
572  p->ComputeDimensions(*this,n,pRep);
573 }
574 
575 
576 //
577 // GetEntityType
578 //
580 {
581  return G4String("G4Polyhedra");
582 }
583 
584 
585 //
586 // Make a clone of the object
587 //
589 {
590  return new G4Polyhedra(*this);
591 }
592 
593 
594 //
595 // Stream object contents to an output stream
596 //
597 std::ostream& G4Polyhedra::StreamInfo( std::ostream& os ) const
598 {
599  G4int oldprc = os.precision(16);
600  os << "-----------------------------------------------------------\n"
601  << " *** Dump for solid - " << GetName() << " ***\n"
602  << " ===================================================\n"
603  << " Solid type: G4Polyhedra\n"
604  << " Parameters: \n"
605  << " starting phi angle : " << startPhi/degree << " degrees \n"
606  << " ending phi angle : " << endPhi/degree << " degrees \n";
607  G4int i=0;
608  if (!genericPgon)
609  {
611  os << " number of Z planes: " << numPlanes << "\n"
612  << " Z values: \n";
613  for (i=0; i<numPlanes; i++)
614  {
615  os << " Z plane " << i << ": "
616  << original_parameters->Z_values[i] << "\n";
617  }
618  os << " Tangent distances to inner surface (Rmin): \n";
619  for (i=0; i<numPlanes; i++)
620  {
621  os << " Z plane " << i << ": "
622  << original_parameters->Rmin[i] << "\n";
623  }
624  os << " Tangent distances to outer surface (Rmax): \n";
625  for (i=0; i<numPlanes; i++)
626  {
627  os << " Z plane " << i << ": "
628  << original_parameters->Rmax[i] << "\n";
629  }
630  }
631  os << " number of RZ points: " << numCorner << "\n"
632  << " RZ values (corners): \n";
633  for (i=0; i<numCorner; i++)
634  {
635  os << " "
636  << corners[i].r << ", " << corners[i].z << "\n";
637  }
638  os << "-----------------------------------------------------------\n";
639  os.precision(oldprc);
640 
641  return os;
642 }
643 
644 
645 //
646 // GetPointOnPlane
647 //
648 // Auxiliary method for get point on surface
649 //
651  G4ThreeVector p2, G4ThreeVector p3) const
652 {
653  G4double lambda1, lambda2, chose,aOne,aTwo;
654  G4ThreeVector t, u, v, w, Area, normal;
655  aOne = 1.;
656  aTwo = 1.;
657 
658  t = p1 - p0;
659  u = p2 - p1;
660  v = p3 - p2;
661  w = p0 - p3;
662 
663  chose = RandFlat::shoot(0.,aOne+aTwo);
664  if( (chose>=0.) && (chose < aOne) )
665  {
666  lambda1 = RandFlat::shoot(0.,1.);
667  lambda2 = RandFlat::shoot(0.,lambda1);
668  return (p2+lambda1*v+lambda2*w);
669  }
670 
671  lambda1 = RandFlat::shoot(0.,1.);
672  lambda2 = RandFlat::shoot(0.,lambda1);
673  return (p0+lambda1*t+lambda2*u);
674 }
675 
676 
677 //
678 // GetPointOnTriangle
679 //
680 // Auxiliary method for get point on surface
681 //
683  G4ThreeVector p2,
684  G4ThreeVector p3) const
685 {
686  G4double lambda1,lambda2;
687  G4ThreeVector v=p3-p1, w=p1-p2;
688 
689  lambda1 = RandFlat::shoot(0.,1.);
690  lambda2 = RandFlat::shoot(0.,lambda1);
691 
692  return (p2 + lambda1*w + lambda2*v);
693 }
694 
695 
696 //
697 // GetPointOnSurface
698 //
700 {
701  if( !genericPgon ) // Polyhedra by faces
702  {
703  G4int j, numPlanes = original_parameters->Num_z_planes, Flag=0;
704  G4double chose, totArea=0., Achose1, Achose2,
705  rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
706  G4double a, b, l2, rang, totalPhi, ksi,
707  area, aTop=0., aBottom=0., zVal=0.;
708 
709  G4ThreeVector p0, p1, p2, p3;
710  std::vector<G4double> aVector1;
711  std::vector<G4double> aVector2;
712  std::vector<G4double> aVector3;
713 
714  totalPhi= (phiIsOpen) ? (endPhi-startPhi) : twopi;
715  ksi = totalPhi/numSide;
716  G4double cosksi = std::cos(ksi/2.);
717 
718  // Below we generate the areas relevant to our solid
719  //
720  for(j=0; j<numPlanes-1; j++)
721  {
722  a = original_parameters->Rmax[j+1];
723  b = original_parameters->Rmax[j];
725  -original_parameters->Z_values[j+1]) + sqr(b-a);
726  area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
727  aVector1.push_back(area);
728  }
729 
730  for(j=0; j<numPlanes-1; j++)
731  {
732  a = original_parameters->Rmin[j+1];//*cosksi;
733  b = original_parameters->Rmin[j];//*cosksi;
735  -original_parameters->Z_values[j+1]) + sqr(b-a);
736  area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
737  aVector2.push_back(area);
738  }
739 
740  for(j=0; j<numPlanes-1; j++)
741  {
742  if(phiIsOpen == true)
743  {
744  aVector3.push_back(0.5*(original_parameters->Rmax[j]
747  -original_parameters->Rmin[j+1])
748  *std::fabs(original_parameters->Z_values[j+1]
750  }
751  else { aVector3.push_back(0.); }
752  }
753 
754  for(j=0; j<numPlanes-1; j++)
755  {
756  totArea += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
757  }
758 
759  // Must include top and bottom areas
760  //
761  if(original_parameters->Rmax[numPlanes-1] != 0.)
762  {
763  a = original_parameters->Rmax[numPlanes-1];
764  b = original_parameters->Rmin[numPlanes-1];
765  l2 = sqr(a-b);
766  aTop = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
767  }
768 
769  if(original_parameters->Rmax[0] != 0.)
770  {
771  a = original_parameters->Rmax[0];
772  b = original_parameters->Rmin[0];
773  l2 = sqr(a-b);
774  aBottom = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
775  }
776 
777  Achose1 = 0.;
778  Achose2 = numSide*(aVector1[0]+aVector2[0])+2.*aVector3[0];
779 
780  chose = RandFlat::shoot(0.,totArea+aTop+aBottom);
781  if( (chose >= 0.) && (chose < aTop + aBottom) )
782  {
783  chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
784  rang = std::floor((chose-startPhi)/ksi-0.01);
785  if(rang<0) { rang=0; }
786  rang = std::fabs(rang);
787  sinphi1 = std::sin(startPhi+rang*ksi);
788  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
789  cosphi1 = std::cos(startPhi+rang*ksi);
790  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
791  chose = RandFlat::shoot(0., aTop + aBottom);
792  if(chose>=0. && chose<aTop)
793  {
794  rad1 = original_parameters->Rmin[numPlanes-1];
795  rad2 = original_parameters->Rmax[numPlanes-1];
796  zVal = original_parameters->Z_values[numPlanes-1];
797  }
798  else
799  {
800  rad1 = original_parameters->Rmin[0];
801  rad2 = original_parameters->Rmax[0];
802  zVal = original_parameters->Z_values[0];
803  }
804  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
805  p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
806  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
807  p3 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
808  return GetPointOnPlane(p0,p1,p2,p3);
809  }
810  else
811  {
812  for (j=0; j<numPlanes-1; j++)
813  {
814  if( ((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes-2) )
815  {
816  Flag = j; break;
817  }
818  Achose1 += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
819  Achose2 = Achose1 + numSide*(aVector1[j+1]+aVector2[j+1])
820  + 2.*aVector3[j+1];
821  }
822  }
823 
824  // At this point we have chosen a subsection
825  // between to adjacent plane cuts...
826 
827  j = Flag;
828 
829  totArea = numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
830  chose = RandFlat::shoot(0.,totArea);
831 
832  if( (chose>=0.) && (chose<numSide*aVector1[j]) )
833  {
834  chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
835  rang = std::floor((chose-startPhi)/ksi-0.01);
836  if(rang<0) { rang=0; }
837  rang = std::fabs(rang);
838  rad1 = original_parameters->Rmax[j];
839  rad2 = original_parameters->Rmax[j+1];
840  sinphi1 = std::sin(startPhi+rang*ksi);
841  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
842  cosphi1 = std::cos(startPhi+rang*ksi);
843  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
844  zVal = original_parameters->Z_values[j];
845 
846  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
847  p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
848 
849  zVal = original_parameters->Z_values[j+1];
850 
851  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
852  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
853  return GetPointOnPlane(p0,p1,p2,p3);
854  }
855  else if ( (chose >= numSide*aVector1[j])
856  && (chose <= numSide*(aVector1[j]+aVector2[j])) )
857  {
858  chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
859  rang = std::floor((chose-startPhi)/ksi-0.01);
860  if(rang<0) { rang=0; }
861  rang = std::fabs(rang);
862  rad1 = original_parameters->Rmin[j];
863  rad2 = original_parameters->Rmin[j+1];
864  sinphi1 = std::sin(startPhi+rang*ksi);
865  sinphi2 = std::sin(startPhi+(rang+1)*ksi);
866  cosphi1 = std::cos(startPhi+rang*ksi);
867  cosphi2 = std::cos(startPhi+(rang+1)*ksi);
868  zVal = original_parameters->Z_values[j];
869 
870  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
871  p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
872 
873  zVal = original_parameters->Z_values[j+1];
874 
875  p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
876  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
877  return GetPointOnPlane(p0,p1,p2,p3);
878  }
879 
880  chose = RandFlat::shoot(0.,2.2);
881  if( (chose>=0.) && (chose < 1.) )
882  {
883  rang = startPhi;
884  }
885  else
886  {
887  rang = endPhi;
888  }
889 
890  cosphi1 = std::cos(rang); rad1 = original_parameters->Rmin[j];
891  sinphi1 = std::sin(rang); rad2 = original_parameters->Rmax[j];
892 
893  p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
895  p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
897 
898  rad1 = original_parameters->Rmax[j+1];
899  rad2 = original_parameters->Rmin[j+1];
900 
901  p2 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
903  p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
905  return GetPointOnPlane(p0,p1,p2,p3);
906  }
907  else // Generic polyhedra
908  {
909  return GetPointOnSurfaceGeneric();
910  }
911 }
912 
913 //
914 // CreatePolyhedron
915 //
917 {
918  if (!genericPgon)
919  {
927  }
928  else
929  {
930  // The following code prepares for:
931  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
932  // const double xyz[][3],
933  // const int faces_vec[][4])
934  // Here is an extract from the header file HepPolyhedron.h:
952  G4int nNodes;
953  G4int nFaces;
954  typedef G4double double3[3];
955  double3* xyz;
956  typedef G4int int4[4];
957  int4* faces_vec;
958  if (phiIsOpen)
959  {
960  // Triangulate open ends. Simple ear-chopping algorithm...
961  // I'm not sure how robust this algorithm is (J.Allison).
962  //
963  std::vector<G4bool> chopped(numCorner, false);
964  std::vector<G4int*> triQuads;
965  G4int remaining = numCorner;
966  G4int iStarter = 0;
967  while (remaining >= 3)
968  {
969  // Find unchopped corners...
970  //
971  G4int A = -1, B = -1, C = -1;
972  G4int iStepper = iStarter;
973  do
974  {
975  if (A < 0) { A = iStepper; }
976  else if (B < 0) { B = iStepper; }
977  else if (C < 0) { C = iStepper; }
978  do
979  {
980  if (++iStepper >= numCorner) iStepper = 0;
981  }
982  while (chopped[iStepper]);
983  }
984  while (C < 0 && iStepper != iStarter);
985 
986  // Check triangle at B is pointing outward (an "ear").
987  // Sign of z cross product determines...
988 
989  G4double BAr = corners[A].r - corners[B].r;
990  G4double BAz = corners[A].z - corners[B].z;
991  G4double BCr = corners[C].r - corners[B].r;
992  G4double BCz = corners[C].z - corners[B].z;
993  if (BAr * BCz - BAz * BCr < kCarTolerance)
994  {
995  G4int* tq = new G4int[3];
996  tq[0] = A + 1;
997  tq[1] = B + 1;
998  tq[2] = C + 1;
999  triQuads.push_back(tq);
1000  chopped[B] = true;
1001  --remaining;
1002  }
1003  else
1004  {
1005  do
1006  {
1007  if (++iStarter >= numCorner) { iStarter = 0; }
1008  }
1009  while (chopped[iStarter]);
1010  }
1011  }
1012 
1013  // Transfer to faces...
1014 
1015  nNodes = (numSide + 1) * numCorner;
1016  nFaces = numSide * numCorner + 2 * triQuads.size();
1017  faces_vec = new int4[nFaces];
1018  G4int iface = 0;
1019  G4int addition = numCorner * numSide;
1020  G4int d = numCorner - 1;
1021  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
1022  {
1023  for (size_t i = 0; i < triQuads.size(); ++i)
1024  {
1025  // Negative for soft/auxiliary/normally invisible edges...
1026  //
1027  G4int a, b, c;
1028  if (iEnd == 0)
1029  {
1030  a = triQuads[i][0];
1031  b = triQuads[i][1];
1032  c = triQuads[i][2];
1033  }
1034  else
1035  {
1036  a = triQuads[i][0] + addition;
1037  b = triQuads[i][2] + addition;
1038  c = triQuads[i][1] + addition;
1039  }
1040  G4int ab = std::abs(b - a);
1041  G4int bc = std::abs(c - b);
1042  G4int ca = std::abs(a - c);
1043  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
1044  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
1045  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
1046  faces_vec[iface][3] = 0;
1047  ++iface;
1048  }
1049  }
1050 
1051  // Continue with sides...
1052 
1053  xyz = new double3[nNodes];
1054  const G4double dPhi = (endPhi - startPhi) / numSide;
1055  G4double phi = startPhi;
1056  G4int ixyz = 0;
1057  for (G4int iSide = 0; iSide < numSide; ++iSide)
1058  {
1059  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1060  {
1061  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1062  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1063  xyz[ixyz][2] = corners[iCorner].z;
1064  if (iCorner < numCorner - 1)
1065  {
1066  faces_vec[iface][0] = ixyz + 1;
1067  faces_vec[iface][1] = ixyz + numCorner + 1;
1068  faces_vec[iface][2] = ixyz + numCorner + 2;
1069  faces_vec[iface][3] = ixyz + 2;
1070  }
1071  else
1072  {
1073  faces_vec[iface][0] = ixyz + 1;
1074  faces_vec[iface][1] = ixyz + numCorner + 1;
1075  faces_vec[iface][2] = ixyz + 2;
1076  faces_vec[iface][3] = ixyz - numCorner + 2;
1077  }
1078  ++iface;
1079  ++ixyz;
1080  }
1081  phi += dPhi;
1082  }
1083 
1084  // Last corners...
1085 
1086  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1087  {
1088  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1089  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1090  xyz[ixyz][2] = corners[iCorner].z;
1091  ++ixyz;
1092  }
1093  }
1094  else // !phiIsOpen - i.e., a complete 360 degrees.
1095  {
1096  nNodes = numSide * numCorner;
1097  nFaces = numSide * numCorner;;
1098  xyz = new double3[nNodes];
1099  faces_vec = new int4[nFaces];
1100  // const G4double dPhi = (endPhi - startPhi) / numSide;
1101  const G4double dPhi = twopi / numSide; // !phiIsOpen endPhi-startPhi = 360 degrees.
1102  G4double phi = startPhi;
1103  G4int ixyz = 0, iface = 0;
1104  for (G4int iSide = 0; iSide < numSide; ++iSide)
1105  {
1106  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1107  {
1108  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1109  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1110  xyz[ixyz][2] = corners[iCorner].z;
1111  if (iSide < numSide - 1)
1112  {
1113  if (iCorner < numCorner - 1)
1114  {
1115  faces_vec[iface][0] = ixyz + 1;
1116  faces_vec[iface][1] = ixyz + numCorner + 1;
1117  faces_vec[iface][2] = ixyz + numCorner + 2;
1118  faces_vec[iface][3] = ixyz + 2;
1119  }
1120  else
1121  {
1122  faces_vec[iface][0] = ixyz + 1;
1123  faces_vec[iface][1] = ixyz + numCorner + 1;
1124  faces_vec[iface][2] = ixyz + 2;
1125  faces_vec[iface][3] = ixyz - numCorner + 2;
1126  }
1127  }
1128  else // Last side joins ends...
1129  {
1130  if (iCorner < numCorner - 1)
1131  {
1132  faces_vec[iface][0] = ixyz + 1;
1133  faces_vec[iface][1] = ixyz + numCorner - nFaces + 1;
1134  faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
1135  faces_vec[iface][3] = ixyz + 2;
1136  }
1137  else
1138  {
1139  faces_vec[iface][0] = ixyz + 1;
1140  faces_vec[iface][1] = ixyz - nFaces + numCorner + 1;
1141  faces_vec[iface][2] = ixyz - nFaces + 2;
1142  faces_vec[iface][3] = ixyz - numCorner + 2;
1143  }
1144  }
1145  ++ixyz;
1146  ++iface;
1147  }
1148  phi += dPhi;
1149  }
1150  }
1151  G4Polyhedron* polyhedron = new G4Polyhedron;
1152  G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
1153  delete [] faces_vec;
1154  delete [] xyz;
1155  if (problem)
1156  {
1157  std::ostringstream message;
1158  message << "Problem creating G4Polyhedron for: " << GetName();
1159  G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
1160  JustWarning, message);
1161  delete polyhedron;
1162  return 0;
1163  }
1164  else
1165  {
1166  return polyhedron;
1167  }
1168  }
1169 }
1170 
1171 
1173 {
1174  G4int numPlanes = (G4int)numCorner;
1175  G4bool isConvertible=true;
1176  G4double Zmax=rz->Bmax();
1177  rz->StartWithZMin();
1178 
1179  // Prepare vectors for storage
1180  //
1181  std::vector<G4double> Z;
1182  std::vector<G4double> Rmin;
1183  std::vector<G4double> Rmax;
1184 
1185  G4int countPlanes=1;
1186  G4int icurr=0;
1187  G4int icurl=0;
1188 
1189  // first plane Z=Z[0]
1190  //
1191  Z.push_back(corners[0].z);
1192  G4double Zprev=Z[0];
1193  if (Zprev == corners[1].z)
1194  {
1195  Rmin.push_back(corners[0].r);
1196  Rmax.push_back (corners[1].r);icurr=1;
1197  }
1198  else if (Zprev == corners[numPlanes-1].z)
1199  {
1200  Rmin.push_back(corners[numPlanes-1].r);
1201  Rmax.push_back (corners[0].r);
1202  icurl=numPlanes-1;
1203  }
1204  else
1205  {
1206  Rmin.push_back(corners[0].r);
1207  Rmax.push_back (corners[0].r);
1208  }
1209 
1210  // next planes until last
1211  //
1212  G4int inextr=0, inextl=0;
1213  for (G4int i=0; i < numPlanes-2; i++)
1214  {
1215  inextr=1+icurr;
1216  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1217 
1218  if((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax)) { break; }
1219 
1220  G4double Zleft = corners[inextl].z;
1221  G4double Zright = corners[inextr].z;
1222  if(Zright>Zleft)
1223  {
1224  Z.push_back(Zleft);
1225  countPlanes++;
1226  G4double difZr=corners[inextr].z - corners[icurr].z;
1227  G4double difZl=corners[inextl].z - corners[icurl].z;
1228 
1229  if(std::fabs(difZl) < kCarTolerance)
1230  {
1231  if(corners[inextl].r >= corners[icurl].r)
1232  {
1233  Rmin.push_back(corners[icurl].r);
1234  Rmax.push_back(Rmax[countPlanes-2]);
1235  Rmax[countPlanes-2]=corners[icurl].r;
1236  }
1237  else
1238  {
1239  Rmin.push_back(corners[inextl].r);
1240  Rmax.push_back(corners[icurl].r);
1241  }
1242  }
1243  else if (difZl >= kCarTolerance)
1244  {
1245  Rmin.push_back(corners[inextl].r);
1246  Rmax.push_back (corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1247  *(corners[inextr].r - corners[icurr].r));
1248  }
1249  else
1250  {
1251  isConvertible=false; break;
1252  }
1253  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1254  }
1255  else if(std::fabs(Zright-Zleft)<kCarTolerance) // Zright=Zleft
1256  {
1257  Z.push_back(Zleft);
1258  countPlanes++;
1259  icurr++;
1260 
1261  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1262 
1263  Rmin.push_back(corners[inextl].r);
1264  Rmax.push_back (corners[inextr].r);
1265  }
1266  else // Zright<Zleft
1267  {
1268  Z.push_back(Zright);
1269  countPlanes++;
1270 
1271  G4double difZr=corners[inextr].z - corners[icurr].z;
1272  G4double difZl=corners[inextl].z - corners[icurl].z;
1273  if(std::fabs(difZr) < kCarTolerance)
1274  {
1275  if(corners[inextr].r >= corners[icurr].r)
1276  {
1277  Rmin.push_back(corners[icurr].r);
1278  Rmax.push_back(corners[inextr].r);
1279  }
1280  else
1281  {
1282  Rmin.push_back(corners[inextr].r);
1283  Rmax.push_back(corners[icurr].r);
1284  Rmax[countPlanes-2]=corners[inextr].r;
1285  }
1286  icurr++;
1287  } // plate
1288  else if (difZr >= kCarTolerance)
1289  {
1290  if(std::fabs(difZl)<kCarTolerance)
1291  {
1292  Rmax.push_back(corners[inextr].r);
1293  Rmin.push_back (corners[icurr].r);
1294  }
1295  else
1296  {
1297  Rmax.push_back(corners[inextr].r);
1298  Rmin.push_back (corners[icurl].r+(Zright-corners[icurl].z)/difZl
1299  * (corners[inextl].r - corners[icurl].r));
1300  }
1301  icurr++;
1302  }
1303  else
1304  {
1305  isConvertible=false; break;
1306  }
1307  }
1308  } // end for loop
1309 
1310  // last plane Z=Zmax
1311  //
1312  Z.push_back(Zmax);
1313  countPlanes++;
1314  inextr=1+icurr;
1315  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1316 
1317  if(corners[inextr].z==corners[inextl].z)
1318  {
1319  Rmax.push_back(corners[inextr].r);
1320  Rmin.push_back(corners[inextl].r);
1321  }
1322  else
1323  {
1324  Rmax.push_back(corners[inextr].r);
1325  Rmin.push_back(corners[inextl].r);
1326  }
1327 
1328  // Set original parameters Rmin,Rmax,Z
1329  //
1330  if(isConvertible)
1331  {
1334  original_parameters->Z_values = new G4double[countPlanes];
1335  original_parameters->Rmin = new G4double[countPlanes];
1336  original_parameters->Rmax = new G4double[countPlanes];
1337 
1338  for(G4int j=0; j < countPlanes; j++)
1339  {
1340  original_parameters->Z_values[j] = Z[j];
1341  original_parameters->Rmax[j] = Rmax[j];
1342  original_parameters->Rmin[j] = Rmin[j];
1343  }
1346  original_parameters->Num_z_planes = countPlanes;
1347 
1348  }
1349  else // Set parameters(r,z) with Rmin==0 as convention
1350  {
1351 #ifdef G4SPECSDEBUG
1352  std::ostringstream message;
1353  message << "Polyhedra " << GetName() << G4endl
1354  << "cannot be converted to Polyhedra with (Rmin,Rmaz,Z) parameters!";
1355  G4Exception("G4Polyhedra::SetOriginalParameters()",
1356  "GeomSolids0002", JustWarning, message);
1357 #endif
1360  original_parameters->Z_values = new G4double[numPlanes];
1361  original_parameters->Rmin = new G4double[numPlanes];
1362  original_parameters->Rmax = new G4double[numPlanes];
1363 
1364  for(G4int j=0; j < numPlanes; j++)
1365  {
1367  original_parameters->Rmax[j] = corners[j].r;
1368  original_parameters->Rmin[j] = 0.0;
1369  }
1372  original_parameters->Num_z_planes = numPlanes;
1373  }
1374  //return isConvertible;
1375 }
1376 
1377 #endif
void CopyStuff(const G4Polyhedra &source)
Definition: G4Polyhedra.cc:429
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:387
G4bool fRebuildPolyhedron
G4Polyhedron * fpPolyhedron
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
G4GeometryType GetEntityType() const
Definition: G4Polyhedra.cc:579
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:650
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:478
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:699
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:568
G4int NumVertices() const
G4double Bmax() const
void ScaleA(G4double scale)
bool G4bool
Definition: G4Types.hh:79
G4double startPhi
Definition: G4Polyhedra.hh:181
#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:540
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
EInside Inside(const G4ThreeVector &p) const
Definition: G4Polyhedra.cc:520
static const G4double ab
G4Polyhedron * CreatePolyhedron() const
Definition: G4Polyhedra.cc:916
EInside
Definition: geomdefs.hh:58
static const double degree
Definition: G4SIunits.hh:125
G4Polyhedra & operator=(const G4Polyhedra &source)
Definition: G4Polyhedra.cc:409
#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:682
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:597
G4VSolid * Clone() const
Definition: G4Polyhedra.cc:588
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