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