Geant4  10.00.p01
UPolycone.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 // UPolycone
12 //
13 // 19.04.13 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include "UUtils.hh"
18 #include <string>
19 #include <cmath>
20 #include <sstream>
21 #include "UPolycone.hh"
22 
23 #include "UEnclosingCylinder.hh"
24 #include "UReduciblePolygon.hh"
25 
26 #include "UTubs.hh"
27 #include "UCons.hh"
28 #include "UTransform3D.hh"
29 
30 using namespace std;
31 
32 UPolycone::UPolycone(const std::string& name,
33  double phiStart,
34  double phiTotal,
35  int numZPlanes,
36  const double zPlane[],
37  const double rInner[],
38  const double rOuter[])
39  : VUSolid(name) //, fNumSides(0)
40 {
41  fCubicVolume = 0;
42  fSurfaceArea = 0;
43  Init(phiStart, phiTotal, numZPlanes, zPlane, rInner, rOuter);
44 
45 }
46 
47 UPolycone::UPolycone(const std::string& name,
48  double phiStart,
49  double phiTotal,
50  int numRZ,
51  const double r[],
52  const double z[])
53  : VUSolid(name)
54 {
55  UReduciblePolygon* rz = new UReduciblePolygon(r, z, numRZ);
56 
57  // Create( phiStart, phiTotal, rz );
58 
59  // Set original_parameters struct for consistency
60  //
61 
62  bool convertible = SetOriginalParameters(rz);
63 
64  if (!convertible)
65  {
66  std::ostringstream message;
67  message << "Polycone " << GetName() << "cannot be converted" << std::endl
68  << "to Polycone with (Rmin,Rmaz,Z) parameters! Use GenericPolycone" ;
69  UUtils::Exception("UPolycone::UPolycone()", "GeomSolids0002",
70  FatalError, 1, message.str().c_str());
71  // JustWarning, message, "Use G4GenericPolycone instead!");
72 
73  }
74  else
75  {
76  std::cout << "INFO: Converting polycone " << GetName() << std::endl
77  << "to optimized polycone with (Rmin,Rmaz,Z) parameters !"
78  << std::endl;
79  double* Z, *R1, *R2;
81  Z = new double[num];
82  R1 = new double[num];
83  R2 = new double[num];
84  for (int i = 0; i < num; i++)
85  {
86  Z[i] = fOriginalParameters->fZValues[i];
87  R1[i] = fOriginalParameters->Rmin[i];
88  R2[i] = fOriginalParameters->Rmax[i];
89  }
90 
91  delete(fOriginalParameters);
92  Init(phiStart, phiTotal, num, Z, R1, R2);
93  delete [] R1;
94  delete [] Z;
95  delete [] R2;
96  }
97 
98  delete rz;
99 }
100 
101 
102 //
103 // Constructor (GEANT3 style parameters)
104 //
105 void UPolycone::Init(double phiStart,
106  double phiTotal,
107  int numZPlanes,
108  const double zPlane[],
109  const double rInner[],
110  const double rOuter[])
111 {
112  //Convertion for angles
113 
114  if (phiTotal <= 0 || phiTotal > UUtils::kTwoPi-1E-10)
115  {
116  phiIsOpen=false;
117  startPhi = 0;
119  }
120  else
121  {
122  //
123  // Convert phi into our convention
124  //
125  phiIsOpen=true;
126  startPhi = phiStart;
127  while( startPhi < 0 ) startPhi += UUtils::kTwoPi;
128 
129  endPhi = phiStart+phiTotal;
130  while( endPhi < startPhi ) endPhi += UUtils::kTwoPi;
131  }
132  // Set Parameters
136  fOriginalParameters->fNumZPlanes = numZPlanes;
137  fOriginalParameters->fZValues.resize(numZPlanes);
138  fOriginalParameters->Rmin.resize(numZPlanes);
139  fOriginalParameters->Rmax.resize(numZPlanes);
140 
141  // Calculate RMax of Polycone in order to determine convexity of sections
142  //
143  double RMaxextent=rOuter[0];
144  for (int j=1; j < numZPlanes; j++)
145  {
146  if (rOuter[j] > RMaxextent) RMaxextent=rOuter[j];
147  if (rInner[j]>rOuter[j])
148  {
149  std::ostringstream message;
150  message << "Cannot create Polycone with rInner > rOuter for the same Z"
151  << std::endl
152  << " rInner > rOuter for the same Z !" << std::endl
153  << " rMin[" << j << "] = " << rInner[j]
154  << " -- rMax[" << j << "] = " << rOuter[j];
155  UUtils::Exception("UPolycone::UPolycone()", "GeomSolids0002",
156  FatalErrorInArguments, 1, message.str().c_str());
157  }
158  }
159  //
160  double prevZ = 0, prevRmax = 0, prevRmin = 0;
161  int dirZ = 1;
162  if (zPlane[1] < zPlane[0])dirZ = -1;
163 // int curSolid = 0;
164 
165  int i;
166  for (i = 0; i < numZPlanes; i++)
167  {
168  if ((i < numZPlanes - 1) && (zPlane[i] == zPlane[i + 1]))
169  {
170  if ((rInner[i] > rOuter[i + 1])
171  || (rInner[i + 1] > rOuter[i]))
172  {
173 
174  std::ostringstream message;
175  message << "Cannot create a Polycone with no contiguous segments."
176  << std::endl
177  << " Segments are not contiguous !" << std::endl
178  << " rMin[" << i << "] = " << rInner[i]
179  << " -- rMax[" << i + 1 << "] = " << rOuter[i + 1] << std::endl
180  << " rMin[" << i + 1 << "] = " << rInner[i + 1]
181  << " -- rMax[" << i << "] = " << rOuter[i];
182  UUtils::Exception("UPolycone::UPolycone()", "GeomSolids0002",
183  FatalErrorInArguments, 1, message.str().c_str());
184  }
185  }
186 
187 
188 
189  double rMin = rInner[i];
190  double rMax = rOuter[i];
191  double z = zPlane[i];
192 
193  if (i > 0)
194  {
195  if (z > prevZ)
196  {
197  if (dirZ < 0)
198  {
199  std::ostringstream message;
200  message << "Cannot create a Polycone with different Z directions.Use GenericPolycone."
201  << std::endl
202  << " ZPlane is changing direction !" << std::endl
203  << " zPlane[0] = " << zPlane[0]
204  << " -- zPlane[1] = " << zPlane[1] << std::endl
205  << " zPlane[" << i - 1 << "] = " << zPlane[i - 1]
206  << " -- rPlane[" << i << "] = " << zPlane[i];
207  UUtils::Exception("UPolycone::UPolycone()", "GeomSolids0002",
208  FatalErrorInArguments, 1, message.str().c_str());
209 
210 
211 
212  }
213  VUSolid* solid;
214  double dz = (z - prevZ) / 2;
215 
216  bool tubular = (rMin == prevRmin && prevRmax == rMax);
217 
218 // if (fNumSides == 0)
219  {
220  if (tubular)
221  {
222  solid = new UTubs("", rMin, rMax, dz, phiStart, phiTotal);
223  }
224  else
225  {
226  solid = new UCons("", prevRmin, prevRmax, rMin, rMax, dz, phiStart, phiTotal);
227  }
228  }
229 // else
230 // {
231 // solid = new UHedra("", prevRmin, prevRmax, rMin, rMax, dz, phiStart, phiTotal, fNumSides);
232 // }
233 
234  fZs.push_back(z);
235 
236  int zi = fZs.size() - 1;
237  double shift = fZs[zi - 1] + 0.5 * (fZs[zi] - fZs[zi - 1]);
238 
239  UPolyconeSection section;
240  section.shift = shift;
241  section.tubular = tubular;
242  section.solid = solid;
243  if(tubular)
244  {
245  if (rMax < RMaxextent) { section.convex = false;}
246  else { section.convex = true;}
247  }
248  else
249  {
250  if ((rMax<prevRmax)||(rMax < RMaxextent)||(prevRmax < RMaxextent))
251  { section.convex = false;}
252  else
253  { section.convex = true;}
254  }
255  fSections.push_back(section);
256  }
257  else
258  {
259  ;// i = i;
260  }
261  }
262  else fZs.push_back(z);
263 
264  fOriginalParameters->fZValues[i] = zPlane[i];
265  fOriginalParameters->Rmin[i] = rInner[i];
266  fOriginalParameters->Rmax[i] = rOuter[i];
267 
268  prevZ = z;
269  prevRmin = rMin;
270  prevRmax = rMax;
271  }
272 
273  fMaxSection = fZs.size() - 2;
274 
275  //
276  // Build RZ polygon using special PCON/PGON GEANT3 constructor
277  //
278  UReduciblePolygon* rz = new UReduciblePolygon(rInner, rOuter, zPlane, numZPlanes);
279 
280  double mxy = rz->Amax();
281 // double alfa = UUtils::kPi / fNumSides;
282 
283  double r = rz->Amax();
284 //
285 // Perform checks of rz values
286 //
287  if (rz->Amin() < 0.0)
288  {
289  std::ostringstream message;
290  message << "Illegal input parameters - " << GetName() << std::endl
291  << " All R values must be >= 0 !";
292  UUtils::Exception("UPolycone::Init()", "GeomSolids0002",
293  FatalErrorInArguments,1, message.str().c_str());
294  }
295 
296 
297  /*
298  if (fNumSides != 0)
299  {
300  // mxy *= std::sqrt(2.0); // this is old and wrong, works only for n = 4
301  // double k = std::tan(alfa) * mxy;
302  double l = mxy / std::cos(alfa);
303  mxy = l;
304  r = l;
305  }
306  */
307 
308  mxy += fgTolerance;
309 
310  fBox.Set(mxy, mxy, (rz->Bmax() - rz->Bmin()) / 2);
311 
312  //
313  // Make enclosingCylinder
314  //
315 
316  enclosingCylinder = new UEnclosingCylinder(r, rz->Bmax(), rz->Bmin(), phiIsOpen, phiStart, phiTotal);
317 
318  delete rz;
319 }
320 
321 
322 /*
323 //
324 // Constructor (generic parameters)
325 //
326 UPolycone3::UPolycone3( const std::string& name,
327  double phiStart,
328  double phiTotal,
329  int numRZ,
330  const double r[],
331  const double z[] )
332  : VUSolid( name )
333 {
334  UReduciblePolygon *rz = new UReduciblePolygon( r, z, numRZ );
335 
336  box.Set(rz->Amax(), rz->Amax(), (rz->Bmax() - rz->Bmin()) /2);
337 
338  // Set fOriginalParameters struct for consistency
339  //
340  SetOriginalParameters();
341 
342  delete rz;
343 }
344 */
345 
346 
347 //
348 // Destructor
349 //
351 {
352  //delete [] corners;
353  //delete fOriginalParameters;
354 }
355 
356 
357 
358 //
359 // Stream object contents to an output stream
360 //
361 std::ostream& UPolycone::StreamInfo(std::ostream& os) const
362 {
363  int oldprc = os.precision(16);
364  os << "-----------------------------------------------------------\n"
365  << " *** Dump for solid - " << GetName() << " ***\n"
366  << " ===================================================\n"
367  << " Solid type: UPolycone3\n"
368  << " Parameters: \n"
369  << " starting phi angle : " << startPhi / (UUtils::kPi / 180.0) << " degrees \n"
370  << " ending phi angle : " << endPhi / (UUtils::kPi / 180.0) << " degrees \n";
371  int i = 0;
372  int numPlanes = fOriginalParameters->fNumZPlanes;
373  os << " number of Z planes: " << numPlanes << "\n"
374  << " Z values: \n";
375  for (i = 0; i < numPlanes; i++)
376  {
377  os << " Z plane " << i << ": "
378  << fOriginalParameters->fZValues[i] << "\n";
379  }
380  os << " Tangent distances to inner surface (Rmin): \n";
381  for (i = 0; i < numPlanes; i++)
382  {
383  os << " Z plane " << i << ": "
384  << fOriginalParameters->Rmin[i] << "\n";
385  }
386  os << " Tangent distances to outer surface (Rmax): \n";
387  for (i = 0; i < numPlanes; i++)
388  {
389  os << " Z plane " << i << ": "
390  << fOriginalParameters->Rmax[i] << "\n";
391  }
392  os << "-----------------------------------------------------------\n";
393  os.precision(oldprc);
394 
395  return os;
396 }
397 
398 
400 {
401  const UPolyconeSection& section = fSections[index];
402  UVector3 ps(p.x, p.y, p.z - section.shift);
403 
404 // if (fNumSides) return section.solid->Inside(ps);
405 
406  double rMinPlus, rMaxPlus; //, rMinMinus, rMaxMinus;
407  double dz;
408  static double halfTolerance = fgTolerance * 0.5;
409 
410  if (section.tubular)
411  {
412  UTubs* tubs = (UTubs*) section.solid;
413  rMinPlus = tubs->GetRMin() + halfTolerance;
414  rMaxPlus = tubs->GetRMax() + halfTolerance;
415  dz = tubs->GetDz();//GetZHalfLength();
416  }
417  else
418  {
419  UCons* cons = (UCons*) section.solid;
420 
421  double rMax1 = cons->GetRmax1();
422  double rMax2 = cons->GetRmax2();
423  double rMin1 = cons->GetRmin1();
424  double rMin2 = cons->GetRmin2();
425 
426  dz = cons->GetDz();
427  double ratio = (ps.z + dz) / (2 * dz);
428  rMinPlus = rMin1 + (rMin2 - rMin1) * ratio + halfTolerance;
429  rMaxPlus = rMax1 + (rMax2 - rMax1) * ratio + halfTolerance;
430  }
431 
432  double rMinMinus = rMinPlus - fgTolerance;
433  double rMaxMinus = rMaxPlus - fgTolerance;
434 
435  double r2 = p.x * p.x + p.y * p.y;
436 
437  if (r2 < rMinMinus * rMinMinus || r2 > rMaxPlus * rMaxPlus) return eOutside;
438  if (r2 < rMinPlus * rMinPlus || r2 > rMaxMinus * rMaxMinus) return eSurface;
439 
440  if (! phiIsOpen )
441  {
442  if (ps.z < -dz + halfTolerance || ps.z > dz - halfTolerance)
443  return eSurface;
444  return eInside;
445  }
446 
447  if (r2 < 1e-10) return eInside;
448 
449  double phi = std::atan2(p.y, p.x); // * UUtils::kTwoPi;
450  if ((phi < 0)||(endPhi > UUtils::kTwoPi)) phi += UUtils::kTwoPi;
451 
452  double ddp = phi - startPhi;
453  if (ddp < 0) ddp += UUtils::kTwoPi;
454  if ((phi <= endPhi + frTolerance)&&(phi>= startPhi-frTolerance))
455  {
456  if (ps.z < -dz + halfTolerance || ps.z > dz - halfTolerance)
457  return eSurface;
458 
459  if (std::fabs(endPhi - phi) < frTolerance)
460  return eSurface;
461  if (std::fabs(startPhi - phi) < frTolerance)
462  return eSurface;
463 
464  return eInside;
465  }
466  return eOutside;
467 }
468 
469 
471 {
472  double shift = fZs[0] + fBox.GetZHalfLength();
473  UVector3 pb(p.x, p.y, p.z - shift);
474  if (fBox.Inside(pb) == eOutside)
475  return eOutside;
476 
477  static const double htolerance = 0.5 * fgTolerance;
478  int index = GetSection(p.z);
479 
480  EnumInside pos = InsideSection(index, p);
481  if (pos == eInside) return eInside;
482 
483  int nextSection;
484  EnumInside nextPos;
485 
486  if (index > 0 && p.z - fZs[index] < htolerance)
487  {
488  nextSection = index - 1;
489  nextPos = InsideSection(nextSection, p);
490  }
491  else if (index < fMaxSection && fZs[index + 1] - p.z < htolerance)
492  {
493  nextSection = index + 1;
494  nextPos = InsideSection(nextSection, p);
495  }
496  else
497  return pos;
498 
499  if (nextPos == eInside) return eInside;
500 
501  if (pos == eSurface && nextPos == eSurface)
502  {
503  UVector3 n, n2;
504  NormalSection(index, p, n);
505  NormalSection(nextSection, p, n2);
506  if ((n + n2).Mag2() < 1000 * frTolerance)
507  return eInside;
508  }
509 
510  return (nextPos == eSurface || pos == eSurface) ? eSurface : eOutside;
511 
512 // return (res == VUSolid::eOutside) ? nextPos : res;
513 }
514 
515 /*
516 if (p.z < fZs.front() - htolerance || p.z > fZs.back() + htolerance) return VUSolid::eOutside;
517 */
518 
520  const UVector3& v, double) const
521 {
522  double shift = fZs[0] + fBox.GetZHalfLength();
523  UVector3 pb(p.x, p.y, p.z - shift);
524 
525  double idistance;
526 
527  idistance = fBox.DistanceToIn(pb, v);
528  // using only box, this appears
529  // to be faster than: idistance = enclosingCylinder->DistanceTo(pb, v);
530  if (idistance >= UUtils::kInfinity) return idistance;
531 
532  // this line can be here or not. not a big difference in performance
533  // TODO: fix enclosingCylinder for polyhedra!!! - the current radius
534  // appears to be too small
535  // if (enclosingCylinder->ShouldMiss(p, v)) return UUtils::kInfinity;
536 
537  // this just takes too much time
538  // idistance = enclosingCylinder->DistanceTo(pb, v);
539  // if (idistance == UUtils::kInfinity) return idistance;
540 
541  pb = p + idistance * v;
542  int index = GetSection(pb.z);
543  pb = p;
544  int increment = (v.z > 0) ? 1 : -1;
545  if (std::fabs(v.z) < fgTolerance) increment = 0;
546 
547  double distance = UUtils::kInfinity;
548  do
549  {
550  const UPolyconeSection& section = fSections[index];
551  pb.z -= section.shift;
552  distance = section.solid->DistanceToIn(pb, v);
553  if (distance < UUtils::kInfinity || !increment)
554  break;
555  index += increment;
556  pb.z += section.shift;
557  }
558  while (index >= 0 && index <= fMaxSection);
559  //if(Inside(p)==eInside)return UUtils::kInfinity;
560  return distance;
561 }
562 
563 double UPolycone::DistanceToOut(const UVector3& p, const UVector3& v,
564  UVector3& n, bool& convex, double /*aPstep*/) const
565 {
566  UVector3 pn(p);
568  {
569  const UPolyconeSection& section = fSections[0];
570  pn.z -= section.shift;
571  return section.solid->DistanceToOut(pn, v, n, convex);
572  }
573  int index = GetSection(p.z);
574  double totalDistance = 0;
575  int increment = (v.z > 0) ? 1 : -1;
576  bool convexloc=true;
577 
578  do
579  {
580  const UPolyconeSection& section = fSections[index];
581  if (totalDistance != 0)
582  {
583  pn = p + (totalDistance /*+ 0 * 1e-8*/) * v; // point must be shifted,
584  // so it could eventually
585  pn.z -= section.shift; // get into another solid
586  if (section.solid->Inside(pn) == eOutside)
587  {
588  break;
589  }
590  }
591  else pn.z -= section.shift;
592 
593  double distance = section.solid->DistanceToOut(pn, v, n, convex);
594  if ((convexloc) && (section.convex)) { convexloc = true; }
595  else { convexloc = false; }
596 
597  totalDistance += distance;
598  index += increment;
599  }
600  while (index >= 0 && index <= fMaxSection);
601 
602  convex=convexloc;
603  if (convex)
604  {
605  // Check final convexity for dz
606  //
607  pn = p + (totalDistance) * v;
608  double dz1 = std::fabs(pn.z-fZs[index]);
609  double dz2 = std::fabs(pn.z-fZs[index+1]);
610  double halfTolerance=0.5*fgTolerance;
611  if((dz1 < halfTolerance) && (index>0)) convex = false;
612  if((dz2 < halfTolerance) && (index<fMaxSection)) convex = false;
613  }
614 
615  return totalDistance;
616 }
617 
618 double UPolycone::SafetyFromInside(const UVector3& p, bool) const
619 {
620  int index = UVoxelizer::BinarySearch(fZs, p.z);
621  if (index < 0 || index > fMaxSection) return 0;
622 
623  double minSafety = SafetyFromInsideSection(index, p);
624  if (minSafety == UUtils::kInfinity) return 0;
625  if (minSafety < 1e-6) return 0;
626 
627  double zbase = fZs[index + 1];
628  for (int i = index + 1; i <= fMaxSection; ++i)
629  {
630  double dz = fZs[i] - zbase;
631  if (dz >= minSafety) break;
632  double safety = SafetyFromOutsideSection(i, p); // safety from Inside cannot be called in this context, because point is not inside, we have to call SafetyFromOutside for given section
633  if (safety < minSafety) minSafety = safety;
634  }
635 
636  if (index > 0)
637  {
638  zbase = fZs[index - 1];
639  for (int i = index - 1; i >= 0; --i)
640  {
641  double dz = zbase - fZs[i];
642  if (dz >= minSafety) break;
643  double safety = SafetyFromOutsideSection(i, p);
644  if (safety < minSafety) minSafety = safety;
645  }
646  }
647  return minSafety;
648 }
649 
650 
651 double UPolycone::SafetyFromOutside(const UVector3& p, bool aAccurate) const
652 {
653  if (!aAccurate)
655 
656  int index = GetSection(p.z);
657  double minSafety = SafetyFromOutsideSection(index, p);
658  if (minSafety < 1e-6) return minSafety;
659 
660  double zbase = fZs[index + 1];
661  for (int i = index + 1; i <= fMaxSection; ++i)
662  {
663  double dz = fZs[i] - zbase;
664  if (dz >= minSafety) break;
665  double safety = SafetyFromOutsideSection(i, p);
666  if (safety < minSafety) minSafety = safety;
667  }
668 
669  zbase = fZs[index - 1];
670  for (int i = index - 1; i >= 0; --i)
671  {
672  double dz = zbase - fZs[i];
673  if (dz >= minSafety) break;
674  double safety = SafetyFromOutsideSection(i, p);
675  if (safety < minSafety) minSafety = safety;
676  }
677  return minSafety;
678 }
679 
680 bool UPolycone::Normal(const UVector3& p, UVector3& n) const
681 {
682  double htolerance = 0.5 * fgTolerance;
683  int index = GetSection(p.z);
684 
685  EnumInside nextPos;
686  int nextSection;
687 
688  if (index > 0 && p.z - fZs[index] < htolerance)
689  {
690  nextSection = index - 1;
691  nextPos = InsideSection(nextSection, p);
692  }
693  else if (index < fMaxSection && fZs[index + 1] - p.z < htolerance)
694  {
695  nextSection = index + 1;
696  nextPos = InsideSection(nextSection, p);
697  }
698  else
699  {
700  const UPolyconeSection& section = fSections[index];
701  UVector3 ps(p.x, p.y, p.z - section.shift);
702  bool res = section.solid->Normal(ps, n);
703 
704  return res;
705 
706  // the code bellow is not used can be deleted
707 
708  nextPos = section.solid->Inside(ps);
709  if (nextPos == eSurface)
710  {
711  return res;
712  }
713  else
714  {
715  //TODO: here should be implementation for case point was not on surface. We would have to look also at other sections. It is not clear if it is possible to solve this problem at all, since we would need precise safety... If it is outside, than it might be OK, but if it is inside..., than I beleive we do not have precise safety
716 
717  // ... or we should at least warn that this case is not supported. actually,
718  // i do not see any algorithm which would obtain right normal of point closest to surface.;
719 
720  return false;
721  }
722  }
723 
724  // even if it says we are on the surface, actually it do not have to be
725 
726 // "TODO special case when point is on the border of two z-sections",
727 // "we should implement this after safety's";
728 
729  EnumInside pos = InsideSection(index, p);
730 
731  if (nextPos == eInside)
732  {
733  //UVector3 n;
734  NormalSection(index, p, n);
735  return false;
736  }
737 
738  if (pos == eSurface && nextPos == eSurface)
739  {
740  //UVector3 n, n2;
741  UVector3 n2;
742  NormalSection(index, p, n);
743  NormalSection(nextSection, p, n2);
744  if ((n + n2).Mag2() < 1000 * frTolerance)
745  {
746  // "we are inside. see TODO above";
747  NormalSection(index, p, n);
748  return false;
749  }
750  }
751 
752  if (nextPos == eSurface || pos == eSurface)
753  {
754  if (pos != eSurface) index = nextSection;
755  bool res = NormalSection(index, p, n);
756  return res;
757  }
758 
759  NormalSection(index, p, n);
760  // "we are outside. see TODO above";
761  return false;
762 }
763 
764 
765 void UPolycone::Extent(UVector3& aMin, UVector3& aMax) const
766 {
767  double r = enclosingCylinder->radius;
768  aMin.Set(-r, -r, fZs.front());
769  aMax.Set(r, r, fZs.back());
770 }
771 
773 {
774  if (fCubicVolume != 0.)
775  {
776  ;
777  }
778  else
779  {
780  for (int i = 0; i < fMaxSection; i++)
781  {
782  UPolyconeSection& section = fSections[i];
783  fCubicVolume += section.solid->Capacity();
784  }
785  }
786  return fCubicVolume;
787 }
788 
790 {
791  if (fSurfaceArea != 0)
792  {
793  ;
794  }
795  else
796  {
797  double Area = 0, totArea = 0;
798  int i = 0;
799  int numPlanes = fOriginalParameters->fNumZPlanes;
800 
801 
802  std::vector<double> areas; // (numPlanes+1);
803  std::vector<UVector3> points; // (numPlanes-1);
804 
805  areas.push_back(UUtils::kPi * (UUtils::sqr(fOriginalParameters->Rmax[0])
807 
808  for (i = 0; i < numPlanes - 1; i++)
809  {
810  Area = (fOriginalParameters->Rmin[i] + fOriginalParameters->Rmin[i + 1])
811  * std::sqrt(UUtils::sqr(fOriginalParameters->Rmin[i]
812  - fOriginalParameters->Rmin[i + 1]) +
815 
816  Area += (fOriginalParameters->Rmax[i] + fOriginalParameters->Rmax[i + 1])
817  * std::sqrt(UUtils::sqr(fOriginalParameters->Rmax[i]
818  - fOriginalParameters->Rmax[i + 1]) +
821 
822  Area *= 0.5 * (endPhi - startPhi);
823 
824  if (startPhi == 0. && endPhi == 2 * UUtils::kPi)
825  {
826  Area += std::fabs(fOriginalParameters->fZValues[i + 1]
829  + fOriginalParameters->Rmax[i + 1]
831  - fOriginalParameters->Rmin[i + 1]);
832  }
833  areas.push_back(Area);
834  totArea += Area;
835  }
836 
837  areas.push_back(UUtils::kPi * (UUtils::sqr(fOriginalParameters->Rmax[numPlanes - 1]) -
838  UUtils::sqr(fOriginalParameters->Rmin[numPlanes - 1])));
839 
840  totArea += (areas[0] + areas[numPlanes]);
841  fSurfaceArea = totArea;
842 
843  }
844 
845  return fSurfaceArea;
846 }
847 
849 //
850 // GetPointOnSurface
851 //
852 // GetPointOnCone
853 //
854 // Auxiliary method for Get Point On Surface
855 //
856 UVector3 UPolycone::GetPointOnCone(double fRmin1, double fRmax1,
857  double fRmin2, double fRmax2,
858  double zOne, double zTwo,
859  double& totArea) const
860 {
861  // declare working variables
862  //
863  double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
864  double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
865  double fDz = (zTwo - zOne) / 2., afDz = std::fabs(fDz);
866  UVector3 point, offset = UVector3(0., 0., 0.5 * (zTwo + zOne));
867  fDPhi = endPhi - startPhi;
868  rone = (fRmax1 - fRmax2) / (2.*fDz);
869  rtwo = (fRmin1 - fRmin2) / (2.*fDz);
870  if (fRmax1 == fRmax2)
871  {
872  qone = 0.;
873  }
874  else
875  {
876  qone = fDz * (fRmax1 + fRmax2) / (fRmax1 - fRmax2);
877  }
878  if (fRmin1 == fRmin2)
879  {
880  qtwo = 0.;
881  }
882  else
883  {
884  qtwo = fDz * (fRmin1 + fRmin2) / (fRmin1 - fRmin2);
885  }
886  Aone = 0.5 * fDPhi * (fRmax2 + fRmax1) * (UUtils::sqr(fRmin1 - fRmin2) + UUtils::sqr(zTwo - zOne));
887  Atwo = 0.5 * fDPhi * (fRmin2 + fRmin1) * (UUtils::sqr(fRmax1 - fRmax2) + UUtils::sqr(zTwo - zOne));
888  Afive = fDz * (fRmax1 - fRmin1 + fRmax2 - fRmin2);
889  totArea = Aone + Atwo + 2.*Afive;
890 
891  phi = UUtils::Random(startPhi, endPhi);
892  cosu = std::cos(phi);
893  sinu = std::sin(phi);
894 
895 
896  if ((startPhi == 0) && (endPhi == 2 * UUtils::kPi))
897  {
898  Afive = 0;
899  }
900  chose = UUtils::Random(0., Aone + Atwo + 2.*Afive);
901  if ((chose >= 0) && (chose < Aone))
902  {
903  if (fRmax1 != fRmax2)
904  {
905  zRand = UUtils::Random(-1.*afDz, afDz);
906  point = UVector3(rone * cosu * (qone - zRand),
907  rone * sinu * (qone - zRand), zRand);
908  }
909  else
910  {
911  point = UVector3(fRmax1 * cosu, fRmax1 * sinu,
912  UUtils::Random(-1.*afDz, afDz));
913 
914  }
915  }
916  else if (chose >= Aone && chose < Aone + Atwo)
917  {
918  if (fRmin1 != fRmin2)
919  {
920  zRand = UUtils::Random(-1.*afDz, afDz);
921  point = UVector3(rtwo * cosu * (qtwo - zRand),
922  rtwo * sinu * (qtwo - zRand), zRand);
923 
924  }
925  else
926  {
927  point = UVector3(fRmin1 * cosu, fRmin1 * sinu,
928  UUtils::Random(-1.*afDz, afDz));
929  }
930  }
931  else if ((chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive))
932  {
933  zRand = UUtils::Random(-1.*afDz, afDz);
934  rmin = fRmin2 - ((zRand - fDz) / (2.*fDz)) * (fRmin1 - fRmin2);
935  rmax = fRmax2 - ((zRand - fDz) / (2.*fDz)) * (fRmax1 - fRmax2);
936  rRand1 = std::sqrt(UUtils::Random() * (UUtils::sqr(rmax) - UUtils::sqr(rmin)) + UUtils::sqr(rmin));
937  point = UVector3(rRand1 * std::cos(startPhi),
938  rRand1 * std::sin(startPhi), zRand);
939  }
940  else
941  {
942  zRand = UUtils::Random(-1.*afDz, afDz);
943  rmin = fRmin2 - ((zRand - fDz) / (2.*fDz)) * (fRmin1 - fRmin2);
944  rmax = fRmax2 - ((zRand - fDz) / (2.*fDz)) * (fRmax1 - fRmax2);
945  rRand1 = std::sqrt(UUtils::Random() * (UUtils::sqr(rmax) - UUtils::sqr(rmin)) + UUtils::sqr(rmin));
946  point = UVector3(rRand1 * std::cos(endPhi),
947  rRand1 * std::sin(endPhi), zRand);
948 
949  }
950 
951  return point + offset;
952 }
953 
954 
955 //
956 // GetPointOnTubs
957 //
958 // Auxiliary method for GetPoint On Surface
959 //
960 UVector3 UPolycone::GetPointOnTubs(double fRMin, double fRMax,
961  double zOne, double zTwo,
962  double& totArea) const
963 {
964  double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
965  aOne, aTwo, aFou, rRand, fDz, fSPhi, fDPhi;
966  fDz = std::fabs(0.5 * (zTwo - zOne));
967  fSPhi = startPhi;
968  fDPhi = endPhi - startPhi;
969 
970  aOne = 2.*fDz * fDPhi * fRMax;
971  aTwo = 2.*fDz * fDPhi * fRMin;
972  aFou = 2.*fDz * (fRMax - fRMin);
973  totArea = aOne + aTwo + 2.*aFou;
974  phi = UUtils::Random(startPhi, endPhi);
975  cosphi = std::cos(phi);
976  sinphi = std::sin(phi);
977  rRand = fRMin + (fRMax - fRMin) * std::sqrt(UUtils::Random());
978 
979  if (startPhi == 0 && endPhi == 2 * UUtils::kPi)
980  aFou = 0;
981 
982  chose = UUtils::Random(0., aOne + aTwo + 2.*aFou);
983  if ((chose >= 0) && (chose < aOne))
984  {
985  xRand = fRMax * cosphi;
986  yRand = fRMax * sinphi;
987  zRand = UUtils::Random(-1.*fDz, fDz);
988  return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
989  }
990  else if ((chose >= aOne) && (chose < aOne + aTwo))
991  {
992  xRand = fRMin * cosphi;
993  yRand = fRMin * sinphi;
994  zRand = UUtils::Random(-1.*fDz, fDz);
995  return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
996  }
997  else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aFou))
998  {
999  xRand = rRand * std::cos(fSPhi + fDPhi);
1000  yRand = rRand * std::sin(fSPhi + fDPhi);
1001  zRand = UUtils::Random(-1.*fDz, fDz);
1002  return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1003  }
1004 
1005  // else
1006 
1007  xRand = rRand * std::cos(fSPhi + fDPhi);
1008  yRand = rRand * std::sin(fSPhi + fDPhi);
1009  zRand = UUtils::Random(-1.*fDz, fDz);
1010  return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1011 }
1012 
1013 
1014 //
1015 // GetPointOnRing
1016 //
1017 // Auxiliary method for GetPoint On Surface
1018 //
1019 UVector3 UPolycone::GetPointOnRing(double fRMin1, double fRMax1,
1020  double fRMin2, double fRMax2,
1021  double zOne) const
1022 {
1023  double xRand, yRand, phi, cosphi, sinphi, rRand1, rRand2, A1, Atot, rCh;
1024  phi = UUtils::Random(startPhi, endPhi);
1025  cosphi = std::cos(phi);
1026  sinphi = std::sin(phi);
1027 
1028  if (fRMin1 == fRMin2)
1029  {
1030  rRand1 = fRMin1;
1031  A1 = 0.;
1032  }
1033  else
1034  {
1035  rRand1 = UUtils::Random(fRMin1, fRMin2);
1036  A1 = std::fabs(fRMin2 * fRMin2 - fRMin1 * fRMin1);
1037  }
1038  if (fRMax1 == fRMax2)
1039  {
1040  rRand2 = fRMax1;
1041  Atot = A1;
1042  }
1043  else
1044  {
1045  rRand2 = UUtils::Random(fRMax1, fRMax2);
1046  Atot = A1 + std::fabs(fRMax2 * fRMax2 - fRMax1 * fRMax1);
1047  }
1048  rCh = UUtils::Random(0., Atot);
1049 
1050  if (rCh > A1)
1051  {
1052  rRand1 = rRand2;
1053  }
1054 
1055  xRand = rRand1 * cosphi;
1056  yRand = rRand1 * sinphi;
1057 
1058  return UVector3(xRand, yRand, zOne);
1059 }
1060 
1061 
1062 //
1063 // GetPointOnCut
1064 //
1065 // Auxiliary method for Get Point On Surface
1066 //
1067 UVector3 UPolycone::GetPointOnCut(double fRMin1, double fRMax1,
1068  double fRMin2, double fRMax2,
1069  double zOne, double zTwo,
1070  double& totArea) const
1071 {
1072  if (zOne == zTwo)
1073  {
1074  return GetPointOnRing(fRMin1, fRMax1, fRMin2, fRMax2, zOne);
1075  }
1076  if ((fRMin1 == fRMin2) && (fRMax1 == fRMax2))
1077  {
1078  return GetPointOnTubs(fRMin1, fRMax1, zOne, zTwo, totArea);
1079  }
1080  return GetPointOnCone(fRMin1, fRMax1, fRMin2, fRMax2, zOne, zTwo, totArea);
1081 }
1082 
1083 
1084 //
1085 // GetPointOnSurface
1086 //
1088 {
1089  double Area = 0, totArea = 0, Achose1 = 0, Achose2 = 0, phi, cosphi, sinphi, rRand;
1090  int i = 0;
1091  int numPlanes = fOriginalParameters->fNumZPlanes;
1092 
1093  phi = UUtils::Random(startPhi, endPhi);
1094  cosphi = std::cos(phi);
1095  sinphi = std::sin(phi);
1096 
1097  rRand = fOriginalParameters->Rmin[0] +
1099  * std::sqrt(UUtils::Random()));
1100 
1101  std::vector<double> areas; // (numPlanes+1);
1102  std::vector<UVector3> points; // (numPlanes-1);
1103 
1104  areas.push_back(UUtils::kPi * (UUtils::sqr(fOriginalParameters->Rmax[0])
1106 
1107  for (i = 0; i < numPlanes - 1; i++)
1108  {
1109  Area = (fOriginalParameters->Rmin[i] + fOriginalParameters->Rmin[i + 1])
1110  * std::sqrt(UUtils::sqr(fOriginalParameters->Rmin[i]
1111  - fOriginalParameters->Rmin[i + 1]) +
1113  - fOriginalParameters->fZValues[i]));
1114 
1115  Area += (fOriginalParameters->Rmax[i] + fOriginalParameters->Rmax[i + 1])
1116  * std::sqrt(UUtils::sqr(fOriginalParameters->Rmax[i]
1117  - fOriginalParameters->Rmax[i + 1]) +
1119  - fOriginalParameters->fZValues[i]));
1120 
1121  Area *= 0.5 * (endPhi - startPhi);
1122 
1123  if (startPhi == 0. && endPhi == 2 * UUtils::kPi)
1124  {
1125  Area += std::fabs(fOriginalParameters->fZValues[i + 1]
1126  - fOriginalParameters->fZValues[i]) *
1128  + fOriginalParameters->Rmax[i + 1]
1130  - fOriginalParameters->Rmin[i + 1]);
1131  }
1132  areas.push_back(Area);
1133  totArea += Area;
1134  }
1135 
1136  areas.push_back(UUtils::kPi * (UUtils::sqr(fOriginalParameters->Rmax[numPlanes - 1]) -
1137  UUtils::sqr(fOriginalParameters->Rmin[numPlanes - 1])));
1138 
1139  totArea += (areas[0] + areas[numPlanes]);
1140  double chose = UUtils::Random(0., totArea);
1141 
1142  if ((chose >= 0.) && (chose < areas[0]))
1143  {
1144  return UVector3(rRand * cosphi, rRand * sinphi,
1146  }
1147 
1148  for (i = 0; i < numPlanes - 1; i++)
1149  {
1150  Achose1 += areas[i];
1151  Achose2 = (Achose1 + areas[i + 1]);
1152  if (chose >= Achose1 && chose < Achose2)
1153  {
1156  fOriginalParameters->Rmin[i + 1],
1157  fOriginalParameters->Rmax[i + 1],
1159  fOriginalParameters->fZValues[i + 1], Area);
1160  }
1161  }
1162 
1163  rRand = fOriginalParameters->Rmin[numPlanes - 1] +
1164  ((fOriginalParameters->Rmax[numPlanes - 1] - fOriginalParameters->Rmin[numPlanes - 1])
1165  * std::sqrt(UUtils::Random()));
1166 
1167  return UVector3(rRand * cosphi, rRand * sinphi,
1168  fOriginalParameters->fZValues[numPlanes - 1]);
1169 
1170 }
1171 
1172 //
1173 // UPolyconeHistorical stuff
1174 //
1175 
1177  : fStartAngle(0.), fOpeningAngle(0.), fNumZPlanes(0),
1178  fZValues(0), Rmin(0), Rmax(0)
1179 {
1180 }
1181 
1183 {
1184 }
1185 
1188 {
1189  fStartAngle = source.fStartAngle;
1190  fOpeningAngle = source.fOpeningAngle;
1191  fNumZPlanes = source.fNumZPlanes;
1192 
1193  fZValues.resize(fNumZPlanes);
1194  Rmin.resize(fNumZPlanes);
1195  Rmax.resize(fNumZPlanes);
1196 
1197  for (int i = 0; i < fNumZPlanes; i++)
1198  {
1199  fZValues[i] = source.fZValues[i];
1200  Rmin[i] = source.Rmin[i];
1201  Rmax[i] = source.Rmax[i];
1202  }
1203 }
1204 
1207 {
1208  if (&right == this) return *this;
1209 
1210  if (&right)
1211  {
1212  fStartAngle = right.fStartAngle;
1213  fOpeningAngle = right.fOpeningAngle;
1214  fNumZPlanes = right.fNumZPlanes;
1215 
1216  fZValues.resize(fNumZPlanes);
1217  Rmin.resize(fNumZPlanes);
1218  Rmax.resize(fNumZPlanes);
1219 
1220  for (int i = 0; i < fNumZPlanes; i++)
1221  {
1222  fZValues[i] = right.fZValues[i];
1223  Rmin[i] = right.Rmin[i];
1224  Rmax[i] = right.Rmax[i];
1225  }
1226  }
1227  return *this;
1228 }
1229 
1231 {
1232  return new UPolycone(*this);
1233 }
1234 //
1235 // Copy constructor
1236 //
1237 UPolycone::UPolycone(const UPolycone& source): VUSolid(source)
1238 {
1239  CopyStuff(source);
1240 }
1241 
1242 
1243 //
1244 // Assignment operator
1245 //
1247 {
1248  if (this == &source) return *this;
1249 
1250  //VUSolid::operator=( source );
1251 
1252  //delete [] corners;
1253 
1254  delete enclosingCylinder;
1255 
1256  CopyStuff(source);
1257 
1258  return *this;
1259 }
1260 //
1261 // CopyStuff
1262 //
1263 void UPolycone::CopyStuff(const UPolycone& source)
1264 {
1265  //
1266  // Simple stuff
1267  //
1268 
1269  startPhi = source.startPhi;
1270  endPhi = source.endPhi;
1271  phiIsOpen = source.phiIsOpen;
1272  fCubicVolume = source.fCubicVolume;
1273  fSurfaceArea = source.fSurfaceArea;
1274  fBox = source.fBox;
1275  //
1276  // The array of planes
1277  //
1279  //
1280  // Enclosing cylinder
1281  //
1283 }
1284 //
1285 // Get Entity Type
1286 //
1288 {
1289  return "Polycone";
1290 }
1291 //
1292 // Set Original Parameters
1293 //
1295 {
1296  int numPlanes = (int)numCorner;
1297  bool isConvertible = true;
1298  double Zmax = rz->Bmax();
1299  rz->StartWithZMin();
1300 
1301  // Prepare vectors for storage
1302  //
1303  std::vector<double> Z;
1304  std::vector<double> Rmin;
1305  std::vector<double> Rmax;
1306 
1307  int countPlanes = 1;
1308  int icurr = 0;
1309  int icurl = 0;
1310 
1311  // first plane Z=Z[0]
1312  //
1313  Z.push_back(corners[0].z);
1314  double Zprev = Z[0];
1315  if (Zprev == corners[1].z)
1316  {
1317  Rmin.push_back(corners[0].r);
1318  Rmax.push_back(corners[1].r);
1319  icurr = 1;
1320  }
1321  else if (Zprev == corners[numPlanes - 1].z)
1322  {
1323  Rmin.push_back(corners[numPlanes - 1].r);
1324  Rmax.push_back(corners[0].r);
1325  icurl = numPlanes - 1;
1326  }
1327  else
1328  {
1329  Rmin.push_back(corners[0].r);
1330  Rmax.push_back(corners[0].r);
1331  }
1332 
1333  // next planes until last
1334  //
1335  int inextr = 0, inextl = 0;
1336  for (int i = 0; i < numPlanes - 2; i++)
1337  {
1338  inextr = 1 + icurr;
1339  inextl = (icurl <= 0) ? numPlanes - 1 : icurl - 1;
1340 
1341  if ((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax))
1342  {
1343  break;
1344  }
1345 
1346  double Zleft = corners[inextl].z;
1347  double Zright = corners[inextr].z;
1348  if (Zright > Zleft)
1349  {
1350  Z.push_back(Zleft);
1351  countPlanes++;
1352  double difZr = corners[inextr].z - corners[icurr].z;
1353  double difZl = corners[inextl].z - corners[icurl].z;
1354 
1355  if (std::fabs(difZl) < frTolerance)
1356  {
1357  if (corners[inextl].r >= corners[icurl].r)
1358  {
1359  Rmin.push_back(corners[icurl].r);
1360  Rmax.push_back(Rmax[countPlanes - 2]);
1361  Rmax[countPlanes - 2] = corners[icurl].r;
1362  }
1363  else
1364  {
1365  Rmin.push_back(corners[inextl].r);
1366  Rmax.push_back(corners[icurl].r);
1367  }
1368  }
1369  else if (difZl >= frTolerance)
1370  {
1371  Rmin.push_back(corners[inextl].r);
1372  Rmax.push_back(corners[icurr].r + (Zleft - corners[icurr].z) / difZr
1373  * (corners[inextr].r - corners[icurr].r));
1374  }
1375  else
1376  {
1377  isConvertible = false;
1378  break;
1379  }
1380  icurl = (icurl == 0) ? numPlanes - 1 : icurl - 1;
1381  }
1382  else if (std::fabs(Zright - Zleft) < frTolerance) // Zright=Zleft
1383  {
1384  Z.push_back(Zleft);
1385  countPlanes++;
1386  icurr++;
1387 
1388  icurl = (icurl == 0) ? numPlanes - 1 : icurl - 1;
1389 
1390  Rmin.push_back(corners[inextl].r);
1391  Rmax.push_back(corners[inextr].r);
1392  }
1393  else // Zright<Zleft
1394  {
1395  Z.push_back(Zright);
1396  countPlanes++;
1397 
1398  double difZr = corners[inextr].z - corners[icurr].z;
1399  double difZl = corners[inextl].z - corners[icurl].z;
1400  if (std::fabs(difZr) < frTolerance)
1401  {
1402  if (corners[inextr].r >= corners[icurr].r)
1403  {
1404  Rmin.push_back(corners[icurr].r);
1405  Rmax.push_back(corners[inextr].r);
1406  }
1407  else
1408  {
1409  Rmin.push_back(corners[inextr].r);
1410  Rmax.push_back(corners[icurr].r);
1411  Rmax[countPlanes - 2] = corners[inextr].r;
1412  }
1413  icurr++;
1414  } // plate
1415  else if (difZr >= frTolerance)
1416  {
1417  if (std::fabs(difZl) < frTolerance)
1418  {
1419  Rmax.push_back(corners[inextr].r);
1420  Rmin.push_back(corners[icurr].r);
1421  }
1422  else
1423  {
1424  Rmax.push_back(corners[inextr].r);
1425  Rmin.push_back(corners[icurl].r + (Zright - corners[icurl].z) / difZl
1426  * (corners[inextl].r - corners[icurl].r));
1427  }
1428  icurr++;
1429  }
1430  else
1431  {
1432  isConvertible = false;
1433  break;
1434  }
1435  }
1436  } // end for loop
1437 
1438  // last plane Z=Zmax
1439  //
1440  Z.push_back(Zmax);
1441  countPlanes++;
1442  inextr = 1 + icurr;
1443  inextl = (icurl <= 0) ? numPlanes - 1 : icurl - 1;
1444 
1445  if (corners[inextr].z == corners[inextl].z)
1446  {
1447  Rmax.push_back(corners[inextr].r);
1448  Rmin.push_back(corners[inextl].r);
1449  }
1450  else
1451  {
1452  Rmax.push_back(corners[inextr].r);
1453  Rmin.push_back(corners[inextl].r);
1454  }
1455 
1456  // Set original parameters Rmin,Rmax,Z
1457  //
1458  if (isConvertible)
1459  {
1461  fOriginalParameters->fZValues.resize(numPlanes);
1462  fOriginalParameters->Rmin.resize(numPlanes);
1463  fOriginalParameters->Rmax.resize(numPlanes);
1464 
1465  for (int j = 0; j < countPlanes; j++)
1466  {
1467  fOriginalParameters->fZValues[j] = Z[j];
1468  fOriginalParameters->Rmax[j] = Rmax[j];
1469  fOriginalParameters->Rmin[j] = Rmin[j];
1470  }
1473  fOriginalParameters->fNumZPlanes = countPlanes;
1474 
1475  }
1476  else // Set parameters(r,z) with Rmin==0 as convention
1477  {
1478  std::ostringstream message;
1479  message << "Polycone " << GetName() << std::endl
1480  << "cannot be converted to Polycone with (Rmin,Rmaz,Z) parameters!";
1481  UUtils::Exception("G4Polycone::SetOriginalParameters()", "GeomSolids0002",
1482  Warning, 1, "can not convert");
1483 
1485 
1486  fOriginalParameters->fZValues.resize(numPlanes);
1487  fOriginalParameters->Rmin.resize(numPlanes);
1488  fOriginalParameters->Rmax.resize(numPlanes);
1489 
1490  for (int j = 0; j < numPlanes; j++)
1491  {
1493  fOriginalParameters->Rmax[j] = corners[j].r;
1494  fOriginalParameters->Rmin[j] = 0.0;
1495  }
1498  fOriginalParameters->fNumZPlanes = numPlanes;
1499  }
1500  return isConvertible;
1501 }
1502 //
1503 // Reset
1504 //
1506 {
1507  //
1508  // Clear old setup
1509  //
1510  delete enclosingCylinder;
1511 
1512  fCubicVolume = 0;
1513  fSurfaceArea = 0;
1514  double phiStart=fOriginalParameters->fStartAngle;
1515  double* Z, *R1, *R2;
1516  int num = fOriginalParameters->fNumZPlanes;
1517  Z = new double[num];
1518  R1 = new double[num];
1519  R2 = new double[num];
1520  for (int i = 0; i < num; i++)
1521  {
1522  Z[i] = fOriginalParameters->fZValues[i];
1523  R1[i] = fOriginalParameters->Rmin[i];
1524  R2[i] = fOriginalParameters->Rmax[i];
1525  }
1526 
1527  Init(phiStart, phiStart+ fOriginalParameters->fOpeningAngle, num, Z, R1, R2);
1528  delete [] R1;
1529  delete [] Z;
1530  delete [] R2;
1531 
1532 }
Definition: UTubs.hh:47
double SafetyFromOutside(const UVector3 &p) const
double GetRMin() const
double Amin() const
double GetRmin1() const
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
Definition: UPolycone.cc:519
std::vector< double > Rmin
Definition: UPolycone.hh:65
static double frTolerance
Definition: VUSolid.hh:31
double GetRmax2() const
std::vector< double > fZValues
Definition: UPolycone.hh:64
UVector3 GetPointOnTubs(double fRMin, double fRMax, double zOne, double zTwo, double &totArea) const
Definition: UPolycone.cc:960
double SafetyFromInsideSection(int index, const UVector3 &p) const
Definition: UPolycone.hh:265
bool NormalSection(int index, const UVector3 &p, UVector3 &n) const
Definition: UPolycone.hh:281
int GetSection(double z) const
Definition: UPolycone.hh:289
G4double z
Definition: TRTMaterials.hh:39
const std::string & GetName() const
Definition: VUSolid.hh:103
double fCubicVolume
Definition: UPolycone.hh:222
static double fgTolerance
Definition: VUSolid.hh:30
G4String name
Definition: TRTMaterials.hh:40
void CopyStuff(const UPolycone &source)
Definition: UPolycone.cc:1263
std::vector< UPolyconeSection > fSections
Definition: UPolycone.hh:260
UVector3 GetPointOnCut(double fRMin1, double fRMax1, double fRMin2, double fRMax2, double zOne, double zTwo, double &totArea) const
Definition: UPolycone.cc:1067
std::ostream & StreamInfo(std::ostream &os) const
Definition: UPolycone.cc:361
UPolyconeHistorical & operator=(const UPolyconeHistorical &right)
Definition: UPolycone.cc:1206
virtual bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const =0
double GetRMax() const
virtual EnumInside Inside(const UVector3 &aPoint) const =0
double x
Definition: UVector3.hh:136
Definition: UCons.hh:49
int numCorner
Definition: UPolycone.hh:219
VUSolid::EnumInside Inside(const UVector3 &p) const
Definition: UPolycone.cc:470
std::vector< double > Rmax
Definition: UPolycone.hh:66
int fMaxSection
Definition: UPolycone.hh:261
bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const
Definition: UPolycone.cc:680
static int BinarySearch(const std::vector< T > &vec, T value)
Definition: UVoxelizer.hh:58
static const double kInfinity
Definition: UUtils.hh:54
double SafetyFromOutsideSection(int index, const UVector3 &p) const
Definition: UPolycone.hh:273
double GetDz() const
bool phiIsOpen
Definition: UPolycone.hh:218
virtual double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const =0
static const double kTwoPi
Definition: UUtils.hh:50
UVector3 GetPointOnSurface() const
Definition: UPolycone.cc:1087
void Init(double phiStart, double phiTotal, int numZPlanes, const double zPlane[], const double rInner[], const double rOuter[])
Definition: UPolycone.cc:105
double endPhi
Definition: UPolycone.hh:217
UPolyconeHistorical * fOriginalParameters
Definition: UPolycone.hh:221
double GetDz() const
UPolyconeSideRZ * corners
Definition: UPolycone.hh:220
UGeometryType GetEntityType() const
Definition: UPolycone.cc:1287
void Reset()
Definition: UPolycone.cc:1505
double Bmin() const
EnumInside
Definition: VUSolid.hh:23
const G4int n
double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const
Definition: UBox.cc:117
double SurfaceArea()
Definition: UPolycone.cc:789
T sqr(const T &x)
Definition: UUtils.hh:142
double Capacity()
Definition: UPolycone.cc:772
double GetZHalfLength() const
Definition: UBox.hh:121
double Amax() const
void SetOriginalParameters()
Definition: UPolycone.hh:226
VUSolid::EnumInside InsideSection(int index, const UVector3 &p) const
Definition: UPolycone.cc:399
void Set(double xx, double yy, double zz)
Definition: UVector3.hh:245
double GetRmax1() const
static const double kPi
Definition: UUtils.hh:49
double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UPolycone.cc:651
double fSurfaceArea
Definition: UPolycone.hh:223
VUSolid * Clone() const
Definition: UPolycone.cc:1230
double SafetyFromInside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UPolycone.cc:618
void Set(double dx, double dy, double dz)
Definition: UBox.cc:43
std::string UGeometryType
Definition: UTypes.hh:70
EnumInside Inside(const UVector3 &aPoint) const
Definition: UBox.cc:97
UVector3 GetPointOnRing(double fRMin, double fRMax, double fRMin2, double fRMax2, double zOne) const
Definition: UPolycone.cc:1019
double Bmax() const
UVector3 GetPointOnCone(double fRmin1, double fRmax1, double fRmin2, double fRmax2, double zOne, double zTwo, double &totArea) const
Definition: UPolycone.cc:856
virtual double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const =0
virtual ~UPolycone()
Definition: UPolycone.cc:350
UPolycone & operator=(const UPolycone &source)
Definition: UPolycone.cc:1246
void Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:177
double z
Definition: UVector3.hh:138
double startPhi
Definition: UPolycone.hh:216
std::vector< double > fZs
Definition: UPolycone.hh:259
UEnclosingCylinder * enclosingCylinder
Definition: UPolycone.hh:248
void Extent(UVector3 &aMin, UVector3 &aMax) const
Definition: UPolycone.cc:765
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
Definition: UPolycone.cc:563
double y
Definition: UVector3.hh:137
UPolycone(const std::string &name)
Definition: UPolycone.hh:82
static const G4double pos
UBox fBox
Definition: UPolycone.hh:224
virtual double Capacity()=0
Definition: VUSolid.cc:199
double GetRmin2() const