Geant4  10.01.p02
UPolyhedra.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UPolyhedra
12 //
13 // --------------------------------------------------------------------
14 //
15 // To be done:
16 // * Cracks: there are probably small cracks in the seams between the
17 // phi face (UPolyPhiFace) and sides (UPolyhedraSide) that are not
18 // entirely leakproof. Also, I am not sure all vertices are leak proof.
19 // * Many optimizations are possible, but not implemented.
20 // * Visualization needs to be updated outside of this routine.
21 //
22 // Utility classes:
23 // * UEnclosingCylinder: I decided a quick check of geometry would be a
24 // good idea (for CPU speed). If the quick check fails, the regular
25 // full-blown UVCSGfaceted version is invoked.
26 // * UReduciblePolygon: Really meant as a check of input parameters,
27 // this utility class also "converts" the GEANT3-like PGON/PCON
28 // arguments into the newer ones.
29 // Both these classes are implemented outside this file because they are
30 // shared with UPolycone.
31 //
32 // 19.09.13 Marek Gayer
33 // Created from original implementation in Geant4
34 // --------------------------------------------------------------------
35 
36 #include "UUtils.hh"
37 #include <string>
38 #include <cmath>
39 #include <sstream>
40 
41 #include "UPolyhedra.hh"
42 #include "UPolyhedraSide.hh"
43 #include "UPolyPhiFace.hh"
44 #include "UEnclosingCylinder.hh"
45 #include "UReduciblePolygon.hh"
46 
47 using namespace std;
48 
49 UPolyhedra::UPolyhedra(const std::string& name,
50  double phiStart,
51  double thePhiTotal,
52  int thefNumSide,
53  int numZPlanes,
54  const double zPlane[],
55  const double rInner[],
56  const double rOuter[])
57  : UVCSGfaceted(name)
58 {
59  Init(phiStart, thePhiTotal, thefNumSide, numZPlanes, zPlane, rInner, rOuter);
60 }
61 
62 //
63 // Constructor (GEANT3 style parameters)
64 //
65 // GEANT3 PGON radii are specified in the distance to the norm of each face.
66 //
68  double phiStart,
69  double thePhiTotal,
70  int thefNumSide,
71  int numZPlanes,
72  const double zPlane[],
73  const double rInner[],
74  const double rOuter[])
75 {
76  fGenericPgon = false;
77 
78  if (thefNumSide <= 0)
79  {
80  std::ostringstream message;
81  message << "Solid must have at least one side - " << GetName() << std::endl
82  << " No sides specified !";
83  UUtils::Exception("UPolyhedra::UPolyhedra()", "GeomSolids0002",
84  UFatalErrorInArguments, 1, message.str().c_str());
85  }
86 
87  //
88  // Calculate conversion factor from G3 radius to U radius
89  //
90  double phiTotal = thePhiTotal;
91  if ((phiTotal <= 0) || (phiTotal >= 2 * UUtils::kPi * (1 - DBL_EPSILON)))
92  {
93  phiTotal = 2 * UUtils::kPi;
94  }
95  double convertRad = std::cos(0.5 * phiTotal / thefNumSide);
96 
97  //
98  // Some historical stuff
99  //
100 // fOriginalParameters = new UPolyhedraHistorical;
101 
102  fOriginalParameters.fNumSide = thefNumSide;
103  fOriginalParameters.fStartAngle = phiStart;
105  fOriginalParameters.fNumZPlanes = numZPlanes;
106  fOriginalParameters.fZValues.resize(numZPlanes);
107  fOriginalParameters.Rmin.resize(numZPlanes);
108  fOriginalParameters.Rmax.resize(numZPlanes);
109 
110  int i;
111  for (i = 0; i < numZPlanes; i++)
112  {
113  if ((i < numZPlanes - 1) && (zPlane[i] == zPlane[i + 1]))
114  {
115  if ((rInner[i] > rOuter[i + 1])
116  || (rInner[i + 1] > rOuter[i]))
117  {
118 
119  std::ostringstream message;
120  message << "Cannot create a Polyhedra with no contiguous segments."
121  << std::endl
122  << " Segments are not contiguous !" << std::endl
123  << " rMin[" << i << "] = " << rInner[i]
124  << " -- rMax[" << i + 1 << "] = " << rOuter[i + 1] << std::endl
125  << " rMin[" << i + 1 << "] = " << rInner[i + 1]
126  << " -- rMax[" << i << "] = " << rOuter[i];
127  UUtils::Exception("UPolyhedra::UPolyhedra()", "GeomSolids0002",
128  UFatalErrorInArguments, 1, message.str().c_str());
129  }
130  }
131  fOriginalParameters.fZValues[i] = zPlane[i];
132  fOriginalParameters.Rmin[i] = rInner[i] / convertRad;
133  fOriginalParameters.Rmax[i] = rOuter[i] / convertRad;
134  }
135 
136 
137  //
138  // Build RZ polygon using special PCON/PGON GEANT3 constructor
139  //
140  UReduciblePolygon* rz =
141  new UReduciblePolygon(rInner, rOuter, zPlane, numZPlanes);
142  rz->ScaleA(1 / convertRad);
143 
144  //
145  // Do the real work
146  //
147  Create(phiStart, phiTotal, thefNumSide, rz);
148 
149  delete rz;
150 }
151 
152 
153 //
154 // Constructor (generic parameters)
155 //
156 UPolyhedra::UPolyhedra(const std::string& name,
157  double phiStart,
158  double phiTotal,
159  int thefNumSide,
160  int numRZ,
161  const double r[],
162  const double z[])
163  : UVCSGfaceted(name), fGenericPgon(true)
164 {
165  UReduciblePolygon* rz = new UReduciblePolygon(r, z, numRZ);
166 
167  Create(phiStart, phiTotal, thefNumSide, rz);
168 
169  // Set fOriginalParameters struct for consistency
170  //
172 
173  delete rz;
174 }
175 
176 
177 //
178 // Create
179 //
180 // Generic create routine, called by each constructor
181 // after conversion of arguments
182 //
183 void UPolyhedra::Create(double phiStart,
184  double phiTotal,
185  int thefNumSide,
186  UReduciblePolygon* rz)
187 {
188  //
189  // Perform checks of rz values
190  //
191  if (rz->Amin() < 0.0)
192  {
193  std::ostringstream message;
194  message << "Illegal input parameters - " << GetName() << std::endl
195  << " All R values must be >= 0 !";
196  UUtils::Exception("UPolyhedra::Create()", "GeomSolids0002",
197  UFatalErrorInArguments, 1, message.str().c_str());
198  }
199 
200  double rzArea = rz->Area();
201  if (rzArea < -VUSolid::Tolerance())
202  rz->ReverseOrder();
203 
204  else if (rzArea < -VUSolid::Tolerance())
205  {
206  std::ostringstream message;
207  message << "Illegal input parameters - " << GetName() << std::endl
208  << " R/Z Cross section is zero or near zero: " << rzArea;
209  UUtils::Exception("UPolyhedra::Create()", "GeomSolids0002",
210  UFatalErrorInArguments, 1, message.str().c_str());
211  }
212 
215  {
216  std::ostringstream message;
217  message << "Illegal input parameters - " << GetName() << std::endl
218  << " Too few unique R/Z values !";
219  UUtils::Exception("UPolyhedra::Create()", "GeomSolids0002",
220  UFatalErrorInArguments, 1, message.str().c_str());
221  }
222 
223  if (rz->CrossesItself(1 / UUtils::kInfinity))
224  {
225  std::ostringstream message;
226  message << "Illegal input parameters - " << GetName() << std::endl
227  << " R/Z segments Cross !";
228  UUtils::Exception("UPolyhedra::Create()", "GeomSolids0002",
229  UFatalErrorInArguments, 1, message.str().c_str());
230  }
231 
232  fNumCorner = rz->NumVertices();
233 
234  fStartPhi = phiStart;
235  while (fStartPhi < 0) fStartPhi += 2 * UUtils::kPi;
236  //
237  // Phi opening? Account for some possible roundoff, and interpret
238  // nonsense value as representing no phi opening
239  //
240  if ((phiTotal <= 0) || (phiTotal > 2 * UUtils::kPi * (1 - DBL_EPSILON)))
241  {
242  fPhiIsOpen = false;
243  fEndPhi = phiStart + 2 * UUtils::kPi;
244  }
245  else
246  {
247  fPhiIsOpen = true;
248 
249  //
250  // Convert phi into our convention
251  //
252  fEndPhi = phiStart + phiTotal;
253  while (fEndPhi < fStartPhi) fEndPhi += 2 * UUtils::kPi;
254  }
255 
256  //
257  // Save number sides
258  //
259  fNumSides = thefNumSide;
260 
261  //
262  // Allocate corner array.
263  //
265 
266  //
267  // Copy fCorners
268  //
269  UReduciblePolygonIterator iterRZ(rz);
270 
271  UPolyhedraSideRZ* next = fCorners;
272  iterRZ.Begin();
273  do
274  {
275  next->r = iterRZ.GetA();
276  next->z = iterRZ.GetB();
277  }
278  while (++next, iterRZ.Next());
279 
280  //
281  // Allocate face pointer array
282  //
284  faces = new UVCSGface*[numFace];
285 
286  //
287  // Construct side faces
288  //
289  // To do so properly, we need to keep track of four successive RZ
290  // fCorners.
291  //
292  // But! Don't construct a face if both points are at zero radius!
293 
294  //
295  UPolyhedraSideRZ* corner = fCorners,
296  *prev = fCorners + fNumCorner - 1,
297  *nextNext;
298  UVCSGface** face = faces;
299  do
300  {
301  next = corner + 1;
302  if (next >= fCorners + fNumCorner) next = fCorners;
303  nextNext = next + 1;
304  if (nextNext >= fCorners + fNumCorner) nextNext = fCorners;
305 
306  if (corner->r < 1 / UUtils::kInfinity && next->r < 1 / UUtils::kInfinity) continue;
307  /*
308  // We must decide here if we can dare declare one of our faces
309  // as having a "valid" normal (i.e. allBehind = true). This
310  // is never possible if the face faces "inward" in r *unless*
311  // we have only one side
312  //
313  bool allBehind;
314  if ((corner->z > next->z) && (fNumSides > 1))
315  {
316  allBehind = false;
317  }
318  else
319  {
320  //
321  // Otherwise, it is only true if the line passing
322  // through the two points of the segment do not
323  // split the r/z Cross section
324  //
325  allBehind = !rz->BisectedBy( corner->r, corner->z,
326  next->r, next->z, VUSolid::Tolerance() );
327  }
328  */
329  *face++ = new UPolyhedraSide(prev, corner, next, nextNext,
331  }
332  while (prev = corner, corner = next, corner > fCorners);
333 
334  if (fPhiIsOpen)
335  {
336  //
337  // Construct phi open edges
338  //
339  *face++ = new UPolyPhiFace(rz, fStartPhi, phiTotal / fNumSides, fEndPhi);
340  *face++ = new UPolyPhiFace(rz, fEndPhi, phiTotal / fNumSides, fStartPhi);
341  }
342 
343  //
344  // We might have dropped a face or two: recalculate numFace
345  //
346  numFace = face - faces;
347 
348  //
349  // Make fEnclosingCylinder
350  //
351 
352  /*
353  double mxy = rz->Amax();
354  double alfa = UUtils::kPi / fNumSides;
355 
356  double r= rz->Amax();
357 
358  if (fNumSides != 0)
359  {
360  // mxy *= std::sqrt(2.0); // this is old and wrong, works only for n = 4
361  double k = std::tan(alfa) * mxy;
362  double l = mxy / std::cos(alfa);
363  mxy = l;
364  r = l;
365  }
366  mxy += fgTolerance;
367  */
368 
370  new UEnclosingCylinder(rz->Amax(), rz->Bmax(), rz->Bmin(), fPhiIsOpen, phiStart, phiTotal);
371 
373 
374  fNoVoxels = fMaxSection < 2; // minimally, sections with at least numbers 0,1,2 values required, this corresponds to fMaxSection == 2
375 }
376 
377 
378 
379 //
380 // Destructor
381 //
383 {
384  delete [] fCorners;
385 // if (fOriginalParameters) delete fOriginalParameters;
386 
387  delete fEnclosingCylinder;
388 }
389 
390 
391 //
392 // Copy constructor
393 //
395  : UVCSGfaceted(source)
396 {
397  CopyStuff(source);
398 }
399 
400 
401 //
402 // Assignment operator
403 //
405 {
406  if (this == &source) return *this;
407 
408  UVCSGfaceted::operator=(source);
409 
410  delete [] fCorners;
411 // if (fOriginalParameters) delete fOriginalParameters;
412 
413  delete fEnclosingCylinder;
414 
415  CopyStuff(source);
416 
417  return *this;
418 }
419 
420 
421 //
422 // CopyStuff
423 //
425 {
426  //
427  // Simple stuff
428  //
429  fNumSides = source.fNumSides;
430  fStartPhi = source.fStartPhi;
431  fEndPhi = source.fEndPhi;
432  fPhiIsOpen = source.fPhiIsOpen;
433  fNumCorner = source.fNumCorner;
434  fGenericPgon = source.fGenericPgon;
435 
436  //
437  // The corner array
438  //
440 
441  UPolyhedraSideRZ* corn = fCorners,
442  *sourceCorn = source.fCorners;
443  do
444  {
445  *corn = *sourceCorn;
446  }
447  while (++sourceCorn, ++corn < fCorners + fNumCorner);
448 
450 
451  //
452  // Enclosing cylinder
453  //
455 }
456 
457 
458 //
459 // Reset
460 //
461 // Recalculates and reshapes the solid, given pre-assigned scaled
462 // fOriginalParameters.
463 //
465 {
466  if (fGenericPgon)
467  {
468  std::ostringstream message;
469  message << "Solid " << GetName() << " built using generic construct."
470  << std::endl << "Not applicable to the generic construct !";
471  UUtils::Exception("UPolyhedra::Reset(,,)", "GeomSolids1001",
472  UWarning, 1, message.str().c_str());
473  return 1;
474  }
475 
476  //
477  // Clear old setup
478  //
480  delete [] fCorners;
481  delete fEnclosingCylinder;
482 
483  //
484  // Rebuild polyhedra
485  //
486  UReduciblePolygon* rz =
494  delete rz;
495 
496  return 0;
497 }
498 
499 
500 //
501 // Inside
502 //
503 // This is an override of UVCSGfaceted::Inside, created in order
504 // to speed things up by first checking with UEnclosingCylinder.
505 //
507 {
508  //
509  // Quick test
510  //
511  if (fEnclosingCylinder->MustBeOutside(p)) return eOutside;
512 
513  //
514  // Long answer
515  //
516  return UVCSGfaceted::Inside(p);
517 }
518 
519 
520 //
521 // DistanceToIn
522 //
523 double UPolyhedra::SafetyFromOutside(const UVector3& aPoint, bool aAccurate) const
524 {
525  return UVCSGfaceted::SafetyFromOutside(aPoint, aAccurate);
526 }
527 
528 
529 //
530 // GetEntityType
531 //
533 {
534  return std::string("Polyhedra");
535 }
536 
537 
538 //
539 // Make a clone of the object
540 //
542 {
543  return new UPolyhedra(*this);
544 }
545 
546 
547 //
548 // Stream object contents to an output stream
549 //
550 std::ostream& UPolyhedra::StreamInfo(std::ostream& os) const
551 {
552  int oldprc = os.precision(16);
553  os << "-----------------------------------------------------------\n"
554  << " *** Dump for solid - " << GetName() << " ***\n"
555  << " ===================================================\n"
556  << " Solid type: UPolyhedra\n"
557  << " Parameters: \n"
558  << " starting phi angle : " << fStartPhi / (UUtils::kPi / 180.0) << " degrees \n"
559  << " ending phi angle : " << fEndPhi / (UUtils::kPi / 180.0) << " degrees \n";
560  int i = 0;
561  if (!fGenericPgon)
562  {
563  int numPlanes = fOriginalParameters.fNumZPlanes;
564  os << " number of Z planes: " << numPlanes << "\n"
565  << " Z values: \n";
566  for (i = 0; i < numPlanes; i++)
567  {
568  os << " Z plane " << i << ": "
569  << fOriginalParameters.fZValues[i] << "\n";
570  }
571  os << " Tangent distances to inner surface (Rmin): \n";
572  for (i = 0; i < numPlanes; i++)
573  {
574  os << " Z plane " << i << ": "
575  << fOriginalParameters.Rmin[i] << "\n";
576  }
577  os << " Tangent distances to outer surface (Rmax): \n";
578  for (i = 0; i < numPlanes; i++)
579  {
580  os << " Z plane " << i << ": "
581  << fOriginalParameters.Rmax[i] << "\n";
582  }
583  }
584  os << " number of RZ points: " << fNumCorner << "\n"
585  << " RZ values (fCorners): \n";
586  for (i = 0; i < fNumCorner; i++)
587  {
588  os << " "
589  << fCorners[i].r << ", " << fCorners[i].z << "\n";
590  }
591  os << "-----------------------------------------------------------\n";
592  os.precision(oldprc);
593 
594  return os;
595 }
596 
597 
598 //
599 // GetPointOnPlane
600 //
601 // Auxiliary method for get point on surface
602 //
604  UVector3 p2, UVector3 p3) const
605 {
606  double lambda1, lambda2, chose, aOne, aTwo;
607  UVector3 t, u, v, w, Area, normal;
608  aOne = 1.;
609  aTwo = 1.;
610 
611  t = p1 - p0;
612  u = p2 - p1;
613  v = p3 - p2;
614  w = p0 - p3;
615 
616  chose = UUtils::Random(0., aOne + aTwo);
617  if ((chose >= 0.) && (chose < aOne))
618  {
619  lambda1 = UUtils::Random(0., 1.);
620  lambda2 = UUtils::Random(0., lambda1);
621  return (p2 + lambda1 * v + lambda2 * w);
622  }
623 
624  lambda1 = UUtils::Random(0., 1.);
625  lambda2 = UUtils::Random(0., lambda1);
626  return (p0 + lambda1 * t + lambda2 * u);
627 }
628 
629 
630 //
631 // GetPointOnTriangle
632 //
633 // Auxiliary method for get point on surface
634 //
636  UVector3 p2,
637  UVector3 p3) const
638 {
639  double lambda1, lambda2;
640  UVector3 v = p3 - p1, w = p1 - p2;
641 
642  lambda1 = UUtils::Random(0., 1.);
643  lambda2 = UUtils::Random(0., lambda1);
644 
645  return (p2 + lambda1 * w + lambda2 * v);
646 }
647 
648 
649 //
650 // GetPointOnSurface
651 //
653 {
654  if (!fGenericPgon) // Polyhedra by faces
655  {
656  int j, numPlanes = fOriginalParameters.fNumZPlanes, Flag = 0;
657  double chose, totArea = 0., Achose1, Achose2,
658  rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
659  double a, b, l2, rang, totalPhi, ksi,
660  area, aTop = 0., aBottom = 0., zVal = 0.;
661 
662  UVector3 p0, p1, p2, p3;
663  std::vector<double> aVector1;
664  std::vector<double> aVector2;
665  std::vector<double> aVector3;
666 
667  totalPhi = (fPhiIsOpen) ? (fEndPhi - fStartPhi) : 2 * UUtils::kPi;
668  ksi = totalPhi / fNumSides;
669  double cosksi = std::cos(ksi / 2.);
670 
671  // Below we generate the areas relevant to our solid
672  //
673  for (j = 0; j < numPlanes - 1; j++)
674  {
675  a = fOriginalParameters.Rmax[j + 1];
676  b = fOriginalParameters.Rmax[j];
678  - fOriginalParameters.fZValues[j + 1]) + UUtils::sqr(b - a);
679  area = std::sqrt(l2 - UUtils::sqr((a - b) * cosksi)) * (a + b) * cosksi;
680  aVector1.push_back(area);
681  }
682 
683  for (j = 0; j < numPlanes - 1; j++)
684  {
685  a = fOriginalParameters.Rmin[j + 1]; //*cosksi;
686  b = fOriginalParameters.Rmin[j];//*cosksi;
688  - fOriginalParameters.fZValues[j + 1]) + UUtils::sqr(b - a);
689  area = std::sqrt(l2 - UUtils::sqr((a - b) * cosksi)) * (a + b) * cosksi;
690  aVector2.push_back(area);
691  }
692 
693  for (j = 0; j < numPlanes - 1; j++)
694  {
695  if (fPhiIsOpen == true)
696  {
697  aVector3.push_back(0.5 * (fOriginalParameters.Rmax[j]
699  + fOriginalParameters.Rmax[j + 1]
700  - fOriginalParameters.Rmin[j + 1])
701  *std::fabs(fOriginalParameters.fZValues[j + 1]
703  }
704  else
705  {
706  aVector3.push_back(0.);
707  }
708  }
709 
710  for (j = 0; j < numPlanes - 1; j++)
711  {
712  totArea += fNumSides * (aVector1[j] + aVector2[j]) + 2.*aVector3[j];
713  }
714 
715  // Must include top and bottom areas
716  //
717  if (fOriginalParameters.Rmax[numPlanes - 1] != 0.)
718  {
719  a = fOriginalParameters.Rmax[numPlanes - 1];
720  b = fOriginalParameters.Rmin[numPlanes - 1];
721  l2 = UUtils::sqr(a - b);
722  aTop = std::sqrt(l2 - UUtils::sqr((a - b) * cosksi)) * (a + b) * cosksi;
723  }
724 
725  if (fOriginalParameters.Rmax[0] != 0.)
726  {
727  a = fOriginalParameters.Rmax[0];
728  b = fOriginalParameters.Rmin[0];
729  l2 = UUtils::sqr(a - b);
730  aBottom = std::sqrt(l2 - UUtils::sqr((a - b) * cosksi)) * (a + b) * cosksi;
731  }
732 
733  Achose1 = 0.;
734  Achose2 = fNumSides * (aVector1[0] + aVector2[0]) + 2.*aVector3[0];
735 
736  chose = UUtils::Random(0., totArea + aTop + aBottom);
737  if ((chose >= 0.) && (chose < aTop + aBottom))
738  {
739  chose = UUtils::Random(fStartPhi, fStartPhi + totalPhi);
740  rang = std::floor((chose - fStartPhi) / ksi - 0.01);
741  if (rang < 0)
742  {
743  rang = 0;
744  }
745  rang = std::fabs(rang);
746  sinphi1 = std::sin(fStartPhi + rang * ksi);
747  sinphi2 = std::sin(fStartPhi + (rang + 1) * ksi);
748  cosphi1 = std::cos(fStartPhi + rang * ksi);
749  cosphi2 = std::cos(fStartPhi + (rang + 1) * ksi);
750  chose = UUtils::Random(0., aTop + aBottom);
751  if (chose >= 0. && chose < aTop)
752  {
753  rad1 = fOriginalParameters.Rmin[numPlanes - 1];
754  rad2 = fOriginalParameters.Rmax[numPlanes - 1];
755  zVal = fOriginalParameters.fZValues[numPlanes - 1];
756  }
757  else
758  {
759  rad1 = fOriginalParameters.Rmin[0];
760  rad2 = fOriginalParameters.Rmax[0];
761  zVal = fOriginalParameters.fZValues[0];
762  }
763  p0 = UVector3(rad1 * cosphi1, rad1 * sinphi1, zVal);
764  p1 = UVector3(rad2 * cosphi1, rad2 * sinphi1, zVal);
765  p2 = UVector3(rad2 * cosphi2, rad2 * sinphi2, zVal);
766  p3 = UVector3(rad1 * cosphi2, rad1 * sinphi2, zVal);
767  return GetPointOnPlane(p0, p1, p2, p3);
768  }
769  else
770  {
771  for (j = 0; j < numPlanes - 1; j++)
772  {
773  if (((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes - 2))
774  {
775  Flag = j;
776  break;
777  }
778  Achose1 += fNumSides * (aVector1[j] + aVector2[j]) + 2.*aVector3[j];
779  Achose2 = Achose1 + fNumSides * (aVector1[j + 1] + aVector2[j + 1])
780  + 2.*aVector3[j + 1];
781  }
782  }
783 
784  // At this point we have chosen a subsection
785  // between to adjacent plane cuts...
786 
787  j = Flag;
788 
789  totArea = fNumSides * (aVector1[j] + aVector2[j]) + 2.*aVector3[j];
790  chose = UUtils::Random(0., totArea);
791 
792  if ((chose >= 0.) && (chose < fNumSides * aVector1[j]))
793  {
794  chose = UUtils::Random(fStartPhi, fStartPhi + totalPhi);
795  rang = std::floor((chose - fStartPhi) / ksi - 0.01);
796  if (rang < 0)
797  {
798  rang = 0;
799  }
800  rang = std::fabs(rang);
801  rad1 = fOriginalParameters.Rmax[j];
802  rad2 = fOriginalParameters.Rmax[j + 1];
803  sinphi1 = std::sin(fStartPhi + rang * ksi);
804  sinphi2 = std::sin(fStartPhi + (rang + 1) * ksi);
805  cosphi1 = std::cos(fStartPhi + rang * ksi);
806  cosphi2 = std::cos(fStartPhi + (rang + 1) * ksi);
807  zVal = fOriginalParameters.fZValues[j];
808 
809  p0 = UVector3(rad1 * cosphi1, rad1 * sinphi1, zVal);
810  p1 = UVector3(rad1 * cosphi2, rad1 * sinphi2, zVal);
811 
812  zVal = fOriginalParameters.fZValues[j + 1];
813 
814  p2 = UVector3(rad2 * cosphi2, rad2 * sinphi2, zVal);
815  p3 = UVector3(rad2 * cosphi1, rad2 * sinphi1, zVal);
816  return GetPointOnPlane(p0, p1, p2, p3);
817  }
818  else if ((chose >= fNumSides * aVector1[j])
819  && (chose <= fNumSides * (aVector1[j] + aVector2[j])))
820  {
821  chose = UUtils::Random(fStartPhi, fStartPhi + totalPhi);
822  rang = std::floor((chose - fStartPhi) / ksi - 0.01);
823  if (rang < 0)
824  {
825  rang = 0;
826  }
827  rang = std::fabs(rang);
828  rad1 = fOriginalParameters.Rmin[j];
829  rad2 = fOriginalParameters.Rmin[j + 1];
830  sinphi1 = std::sin(fStartPhi + rang * ksi);
831  sinphi2 = std::sin(fStartPhi + (rang + 1) * ksi);
832  cosphi1 = std::cos(fStartPhi + rang * ksi);
833  cosphi2 = std::cos(fStartPhi + (rang + 1) * ksi);
834  zVal = fOriginalParameters.fZValues[j];
835 
836  p0 = UVector3(rad1 * cosphi1, rad1 * sinphi1, zVal);
837  p1 = UVector3(rad1 * cosphi2, rad1 * sinphi2, zVal);
838 
839  zVal = fOriginalParameters.fZValues[j + 1];
840 
841  p2 = UVector3(rad2 * cosphi2, rad2 * sinphi2, zVal);
842  p3 = UVector3(rad2 * cosphi1, rad2 * sinphi1, zVal);
843  return GetPointOnPlane(p0, p1, p2, p3);
844  }
845 
846  chose = UUtils::Random(0., 2.2);
847  if ((chose >= 0.) && (chose < 1.))
848  {
849  rang = fStartPhi;
850  }
851  else
852  {
853  rang = fEndPhi;
854  }
855 
856  cosphi1 = std::cos(rang);
857  rad1 = fOriginalParameters.Rmin[j];
858  sinphi1 = std::sin(rang);
859  rad2 = fOriginalParameters.Rmax[j];
860 
861  p0 = UVector3(rad1 * cosphi1, rad1 * sinphi1,
863  p1 = UVector3(rad2 * cosphi1, rad2 * sinphi1,
865 
866  rad1 = fOriginalParameters.Rmax[j + 1];
867  rad2 = fOriginalParameters.Rmin[j + 1];
868 
869  p2 = UVector3(rad1 * cosphi1, rad1 * sinphi1,
871  p3 = UVector3(rad2 * cosphi1, rad2 * sinphi1,
873  return GetPointOnPlane(p0, p1, p2, p3);
874  }
875  else // Generic polyhedra
876  {
877  return GetPointOnSurfaceGeneric();
878  }
879 }
880 
881 //
882 // UPolyhedraHistorical stuff
883 //
885  : fStartAngle(0.), fOpeningAngle(0.), fNumSide(0), fNumZPlanes(0),
886  fZValues(0), Rmin(0), Rmax(0)
887 {
888 }
889 
891 {
892 }
893 
896 {
897  fStartAngle = source.fStartAngle;
898  fOpeningAngle = source.fOpeningAngle;
899  fNumSide = source.fNumSide;
900  fNumZPlanes = source.fNumZPlanes;
901 
902  fZValues = source.fZValues;
903  Rmin = source.Rmin;
904  Rmax = source.Rmax;
905 }
906 
909 {
910  if (&right == this) return *this;
911 
912  fStartAngle = right.fStartAngle;
914  fNumSide = right.fNumSide;
915  fNumZPlanes = right.fNumZPlanes;
916 
917  fZValues = right.fZValues;
918  Rmin = right.Rmin;
919  Rmax = right.Rmax;
920  return *this;
921 }
922 
923 void UPolyhedra::Extent(UVector3& aMin, UVector3& aMax) const
924 {
925  fEnclosingCylinder->Extent(aMin, aMax);
926 }
927 
928 
929 //
930 // DistanceToIn
931 //
932 // This is an override of G4VCSGfaceted::Inside, created in order
933 // to speed things up by first checking with G4EnclosingCylinder.
934 //
936  const UVector3& v, double aPstep) const
937 {
938  //
939  // Quick test
940  //
942  return UUtils::kInfinity;
943 
944  //
945  // Long answer
946  //
947  return UVCSGfaceted::DistanceToIn(p, v, aPstep);
948 }
double fEndPhi
Definition: UPolyhedra.hh:183
bool Reset()
Definition: UPolyhedra.cc:464
double Amin() const
std::ostream & StreamInfo(std::ostream &os) const
Definition: UPolyhedra.cc:550
int fNumCorner
Definition: UPolyhedra.hh:186
bool RemoveRedundantVertices(double tolerance)
bool fPhiIsOpen
Definition: UPolyhedra.hh:184
UVector3 GetPointOnPlane(UVector3 p0, UVector3 p1, UVector3 p2, UVector3 p3) const
Definition: UPolyhedra.cc:603
bool MustBeOutside(const UVector3 &p) const
G4double z
Definition: TRTMaterials.hh:39
const std::string & GetName() const
Definition: VUSolid.hh:103
UPolyhedraHistorical & operator=(const UPolyhedraHistorical &right)
Definition: UPolyhedra.cc:908
G4String name
Definition: TRTMaterials.hh:40
double fStartPhi
Definition: UPolyhedra.hh:182
std::vector< double > fZValues
Definition: UPolyhedra.hh:59
void DeleteStuff()
G4double a
Definition: TRTMaterials.hh:39
bool RemoveDuplicateVertices(double tolerance)
void Extent(UVector3 &aMin, UVector3 &aMax) const
Definition: UPolyhedra.cc:923
static double Tolerance()
Definition: VUSolid.hh:127
void CopyStuff(const UPolyhedra &source)
Definition: UPolyhedra.cc:424
VUSolid * Clone() const
Definition: UPolyhedra.cc:541
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
UPolyhedra & operator=(const UPolyhedra &source)
Definition: UPolyhedra.cc:404
UVector3 GetPointOnTriangle(UVector3 p0, UVector3 p1, UVector3 p2) const
Definition: UPolyhedra.cc:635
static const double kInfinity
Definition: UUtils.hh:54
void InitVoxels(UReduciblePolygon &z, double radius)
void ScaleA(double scale)
UVCSGfaceted & operator=(const UVCSGfaceted &source)
Definition: UVCSGfaceted.cc:63
int NumVertices() const
UVector3 GetPointOnSurface() const
Definition: UPolyhedra.cc:652
const G4double p2
const G4double p1
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
Definition: UPolyhedra.cc:935
#define DBL_EPSILON
Definition: templates.hh:87
std::vector< double > Rmin
Definition: UPolyhedra.hh:60
double Bmin() const
EnumInside
Definition: VUSolid.hh:23
double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UPolyhedra.cc:523
UEnclosingCylinder * fEnclosingCylinder
Definition: UPolyhedra.hh:189
UGeometryType GetEntityType() const
Definition: UPolyhedra.cc:532
UPolyhedraSideRZ * fCorners
Definition: UPolyhedra.hh:187
void SetOriginalParameters()
T sqr(const T &x)
Definition: UUtils.hh:138
bool ShouldMiss(const UVector3 &p, const UVector3 &v) const
bool CrossesItself(double tolerance)
virtual ~UPolyhedra()
Definition: UPolyhedra.cc:382
void Exception(const char *originOfException, const char *exceptionCode, UExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
double Amax() const
const G4double p0
void Init(double phiStart, double phiTotal, int numSide, int numZPlanes, const double zPlane[], const double rInner[], const double rOuter[])
Definition: UPolyhedra.cc:67
UPolyhedra(const std::string &name)
Definition: UPolyhedra.hh:68
static const double kPi
Definition: UUtils.hh:49
void Create(double phiStart, double phiTotal, int numSide, UReduciblePolygon *rz)
Definition: UPolyhedra.cc:183
std::string UGeometryType
Definition: UTypes.hh:39
std::vector< double > Rmax
Definition: UPolyhedra.hh:61
double Bmax() const
UVCSGface ** faces
UPolyhedraHistorical fOriginalParameters
Definition: UPolyhedra.hh:188
void Extent(UVector3 &aMin, UVector3 &aMax) const
UVector3 GetPointOnSurfaceGeneric() const
virtual double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
virtual double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
virtual VUSolid::EnumInside Inside(const UVector3 &p) const
bool fGenericPgon
Definition: UPolyhedra.hh:185
VUSolid::EnumInside Inside(const UVector3 &p) const
Definition: UPolyhedra.cc:506