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