Geant4  10.01.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  }
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  FatalErrorInArguments, 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  FatalErrorInArguments, 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  FatalErrorInArguments, 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  FatalErrorInArguments,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  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  //if(Inside(p)==eInside)return UUtils::kInfinity;
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  if(n.Mag()<0.9)std::cout<<"UPcon Norm 1 ="<<n<<" res="<<res<<" p="<<p<<" name="<<GetName()<<std::endl;
819  return res;
820 
821  // the code bellow is not used can be deleted
822 
823  nextPos = section.solid->Inside(ps);
824  if (nextPos == eSurface)
825  {
826  if(n.Mag()<0.9)std::cout<<"UPcon Norm 2 ="<<n<<" res="<<res<<" p="<<p<<" name="<<GetName()<<std::endl;
827  return res;
828  }
829  else
830  {
831  //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
832 
833  // ... or we should at least warn that this case is not supported. actually,
834  // i do not see any algorithm which would obtain right normal of point closest to surface.;
835  if(n.Mag()<0.9)std::cout<<"UPcon Norm 3 ="<<n<<std::endl;
836  return false;
837  }
838  }
839 
840  // even if it says we are on the surface, actually it do not have to be
841 
842 // "TODO special case when point is on the border of two z-sections",
843 // "we should implement this after safety's";
844 
845  EnumInside pos = InsideSection(index, p);
846 
847  if (nextPos == eInside)
848  {
849  //UVector3 n;
850  NormalSection(index, p, n);
851  if(n.Mag()<0.9)std::cout<<"UPcon Norm 4 ="<<n<<std::endl;
852  return false;
853  }
854 
855  if (pos == eSurface && nextPos == eSurface)
856  {
857  //UVector3 n, n2;
858  UVector3 n2;
859  NormalSection(index, p, n);
860  NormalSection(nextSection, p, n2);
861  if ((n + n2).Mag2() < 1000 * frTolerance)
862  {
863  // "we are inside. see TODO above";
864 
865  NormalSection(index, p, n);
866  if(n.Mag()<0.9)std::cout<<"UPcon Norm 5 ="<<n<<std::endl;
867  return false;
868  }
869  }
870 
871  if (nextPos == eSurface || pos == eSurface)
872  {
873  if (pos != eSurface) index = nextSection;
874  bool res = NormalSection(index, p, n);
875  if(n.Mag()<0.9)std::cout<<"UPcon Norm 6 ="<<n<<std::endl;
876  return res;
877  }
878 
879  NormalSection(index, p, n);
880  if(n.Mag()<0.9)std::cout<<"UPcon Norm 7 ="<<n<<std::endl;
881  // "we are outside. see TODO above";
882  return false;
883 }
884 
885 
886 void UPolycone::Extent(UVector3& aMin, UVector3& aMax) const
887 {
888  double r = enclosingCylinder->radius;
889  aMin.Set(-r, -r, fZs.front());
890  aMax.Set(r, r, fZs.back());
891  //std::cout<<"UPOLy Extent "<<aMin<<" 2="<<" name="<<aMax<<GetName()<<std::endl;
892 }
893 
895 {
896  if (fCubicVolume != 0.)
897  {
898  ;
899  }
900  else
901  {
902  for (int i = 0; i < fMaxSection+1; i++)
903  {
904  UPolyconeSection& section = fSections[i];
905  fCubicVolume += section.solid->Capacity();
906  }
907  }
908  return fCubicVolume;
909 }
910 
912 {
913  if (fSurfaceArea != 0)
914  {
915  ;
916  }
917  else
918  {
919  double Area = 0, totArea = 0;
920  int i = 0;
921  int numPlanes = fOriginalParameters->fNumZPlanes;
922 
923 
924  std::vector<double> areas; // (numPlanes+1);
925  std::vector<UVector3> points; // (numPlanes-1);
926 
927  areas.push_back(UUtils::kPi * (UUtils::sqr(fOriginalParameters->Rmax[0])
929 
930  for (i = 0; i < numPlanes - 1; i++)
931  {
932  Area = (fOriginalParameters->Rmin[i] + fOriginalParameters->Rmin[i + 1])
933  * std::sqrt(UUtils::sqr(fOriginalParameters->Rmin[i]
934  - fOriginalParameters->Rmin[i + 1]) +
937 
938  Area += (fOriginalParameters->Rmax[i] + fOriginalParameters->Rmax[i + 1])
939  * std::sqrt(UUtils::sqr(fOriginalParameters->Rmax[i]
940  - fOriginalParameters->Rmax[i + 1]) +
943 
944  Area *= 0.5 * (endPhi - startPhi);
945 
947  {
948  Area += std::fabs(fOriginalParameters->fZValues[i + 1]
951  + fOriginalParameters->Rmax[i + 1]
953  - fOriginalParameters->Rmin[i + 1]);
954  }
955  areas.push_back(Area);
956  totArea += Area;
957  }
958 
959  areas.push_back(UUtils::kPi * (UUtils::sqr(fOriginalParameters->Rmax[numPlanes - 1]) -
960  UUtils::sqr(fOriginalParameters->Rmin[numPlanes - 1])));
961 
962  totArea += (areas[0] + areas[numPlanes]);
963  fSurfaceArea = totArea;
964 
965  }
966 
967  return fSurfaceArea;
968 }
969 
971 //
972 // GetPointOnSurface
973 //
974 // GetPointOnCone
975 //
976 // Auxiliary method for Get Point On Surface
977 //
978 UVector3 UPolycone::GetPointOnCone(double fRmin1, double fRmax1,
979  double fRmin2, double fRmax2,
980  double zOne, double zTwo,
981  double& totArea) const
982 {
983  // declare working variables
984  //
985  double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
986  double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
987  double fDz = (zTwo - zOne) / 2., afDz = std::fabs(fDz);
988  UVector3 point, offset = UVector3(0., 0., 0.5 * (zTwo + zOne));
989  fDPhi = endPhi - startPhi;
990  rone = (fRmax1 - fRmax2) / (2.*fDz);
991  rtwo = (fRmin1 - fRmin2) / (2.*fDz);
992  if (fRmax1 == fRmax2)
993  {
994  qone = 0.;
995  }
996  else
997  {
998  qone = fDz * (fRmax1 + fRmax2) / (fRmax1 - fRmax2);
999  }
1000  if (fRmin1 == fRmin2)
1001  {
1002  qtwo = 0.;
1003  }
1004  else
1005  {
1006  qtwo = fDz * (fRmin1 + fRmin2) / (fRmin1 - fRmin2);
1007  }
1008  Aone = 0.5 * fDPhi * (fRmax2 + fRmax1) * (UUtils::sqr(fRmin1 - fRmin2) + UUtils::sqr(zTwo - zOne));
1009  Atwo = 0.5 * fDPhi * (fRmin2 + fRmin1) * (UUtils::sqr(fRmax1 - fRmax2) + UUtils::sqr(zTwo - zOne));
1010  Afive = fDz * (fRmax1 - fRmin1 + fRmax2 - fRmin2);
1011  totArea = Aone + Atwo + 2.*Afive;
1012 
1013  phi = UUtils::Random(startPhi, endPhi);
1014  cosu = std::cos(phi);
1015  sinu = std::sin(phi);
1016 
1017 
1018  if ((startPhi == 0) && (endPhi == 2 * UUtils::kPi))
1019  {
1020  Afive = 0;
1021  }
1022  chose = UUtils::Random(0., Aone + Atwo + 2.*Afive);
1023  if ((chose >= 0) && (chose < Aone))
1024  {
1025  if (fRmax1 != fRmax2)
1026  {
1027  zRand = UUtils::Random(-1.*afDz, afDz);
1028  point = UVector3(rone * cosu * (qone - zRand),
1029  rone * sinu * (qone - zRand), zRand);
1030  }
1031  else
1032  {
1033  point = UVector3(fRmax1 * cosu, fRmax1 * sinu,
1034  UUtils::Random(-1.*afDz, afDz));
1035 
1036  }
1037  }
1038  else if (chose >= Aone && chose < Aone + Atwo)
1039  {
1040  if (fRmin1 != fRmin2)
1041  {
1042  zRand = UUtils::Random(-1.*afDz, afDz);
1043  point = UVector3(rtwo * cosu * (qtwo - zRand),
1044  rtwo * sinu * (qtwo - zRand), zRand);
1045 
1046  }
1047  else
1048  {
1049  point = UVector3(fRmin1 * cosu, fRmin1 * sinu,
1050  UUtils::Random(-1.*afDz, afDz));
1051  }
1052  }
1053  else if ((chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive))
1054  {
1055  zRand = UUtils::Random(-1.*afDz, afDz);
1056  rmin = fRmin2 - ((zRand - fDz) / (2.*fDz)) * (fRmin1 - fRmin2);
1057  rmax = fRmax2 - ((zRand - fDz) / (2.*fDz)) * (fRmax1 - fRmax2);
1058  rRand1 = std::sqrt(UUtils::Random() * (UUtils::sqr(rmax) - UUtils::sqr(rmin)) + UUtils::sqr(rmin));
1059  point = UVector3(rRand1 * std::cos(startPhi),
1060  rRand1 * std::sin(startPhi), zRand);
1061  }
1062  else
1063  {
1064  zRand = UUtils::Random(-1.*afDz, afDz);
1065  rmin = fRmin2 - ((zRand - fDz) / (2.*fDz)) * (fRmin1 - fRmin2);
1066  rmax = fRmax2 - ((zRand - fDz) / (2.*fDz)) * (fRmax1 - fRmax2);
1067  rRand1 = std::sqrt(UUtils::Random() * (UUtils::sqr(rmax) - UUtils::sqr(rmin)) + UUtils::sqr(rmin));
1068  point = UVector3(rRand1 * std::cos(endPhi),
1069  rRand1 * std::sin(endPhi), zRand);
1070 
1071  }
1072 
1073  return point + offset;
1074 }
1075 
1076 
1077 //
1078 // GetPointOnTubs
1079 //
1080 // Auxiliary method for GetPoint On Surface
1081 //
1082 UVector3 UPolycone::GetPointOnTubs(double fRMin, double fRMax,
1083  double zOne, double zTwo,
1084  double& totArea) const
1085 {
1086  double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1087  aOne, aTwo, aFou, rRand, fDz, fSPhi, fDPhi;
1088  fDz = std::fabs(0.5 * (zTwo - zOne));
1089  fSPhi = startPhi;
1090  fDPhi = endPhi - startPhi;
1091 
1092  aOne = 2.*fDz * fDPhi * fRMax;
1093  aTwo = 2.*fDz * fDPhi * fRMin;
1094  aFou = 2.*fDz * (fRMax - fRMin);
1095  totArea = aOne + aTwo + 2.*aFou;
1096  phi = UUtils::Random(startPhi, endPhi);
1097  cosphi = std::cos(phi);
1098  sinphi = std::sin(phi);
1099  rRand = fRMin + (fRMax - fRMin) * std::sqrt(UUtils::Random());
1100 
1101  if (startPhi == 0 && endPhi == 2 * UUtils::kPi)
1102  aFou = 0;
1103 
1104  chose = UUtils::Random(0., aOne + aTwo + 2.*aFou);
1105  if ((chose >= 0) && (chose < aOne))
1106  {
1107  xRand = fRMax * cosphi;
1108  yRand = fRMax * sinphi;
1109  zRand = UUtils::Random(-1.*fDz, fDz);
1110  return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1111  }
1112  else if ((chose >= aOne) && (chose < aOne + aTwo))
1113  {
1114  xRand = fRMin * cosphi;
1115  yRand = fRMin * sinphi;
1116  zRand = UUtils::Random(-1.*fDz, fDz);
1117  return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1118  }
1119  else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aFou))
1120  {
1121  xRand = rRand * std::cos(fSPhi + fDPhi);
1122  yRand = rRand * std::sin(fSPhi + fDPhi);
1123  zRand = UUtils::Random(-1.*fDz, fDz);
1124  return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1125  }
1126 
1127  // else
1128 
1129  xRand = rRand * std::cos(fSPhi + fDPhi);
1130  yRand = rRand * std::sin(fSPhi + fDPhi);
1131  zRand = UUtils::Random(-1.*fDz, fDz);
1132  return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1133 }
1134 
1135 
1136 //
1137 // GetPointOnRing
1138 //
1139 // Auxiliary method for GetPoint On Surface
1140 //
1141 UVector3 UPolycone::GetPointOnRing(double fRMin1, double fRMax1,
1142  double fRMin2, double fRMax2,
1143  double zOne) const
1144 {
1145  double xRand, yRand, phi, cosphi, sinphi, rRand1, rRand2, A1, Atot, rCh;
1146  phi = UUtils::Random(startPhi, endPhi);
1147  cosphi = std::cos(phi);
1148  sinphi = std::sin(phi);
1149 
1150  if (fRMin1 == fRMin2)
1151  {
1152  rRand1 = fRMin1;
1153  A1 = 0.;
1154  }
1155  else
1156  {
1157  rRand1 = UUtils::Random(fRMin1, fRMin2);
1158  A1 = std::fabs(fRMin2 * fRMin2 - fRMin1 * fRMin1);
1159  }
1160  if (fRMax1 == fRMax2)
1161  {
1162  rRand2 = fRMax1;
1163  Atot = A1;
1164  }
1165  else
1166  {
1167  rRand2 = UUtils::Random(fRMax1, fRMax2);
1168  Atot = A1 + std::fabs(fRMax2 * fRMax2 - fRMax1 * fRMax1);
1169  }
1170  rCh = UUtils::Random(0., Atot);
1171 
1172  if (rCh > A1)
1173  {
1174  rRand1 = rRand2;
1175  }
1176 
1177  xRand = rRand1 * cosphi;
1178  yRand = rRand1 * sinphi;
1179 
1180  return UVector3(xRand, yRand, zOne);
1181 }
1182 
1183 
1184 //
1185 // GetPointOnCut
1186 //
1187 // Auxiliary method for Get Point On Surface
1188 //
1189 UVector3 UPolycone::GetPointOnCut(double fRMin1, double fRMax1,
1190  double fRMin2, double fRMax2,
1191  double zOne, double zTwo,
1192  double& totArea) const
1193 {
1194  if (zOne == zTwo)
1195  {
1196  return GetPointOnRing(fRMin1, fRMax1, fRMin2, fRMax2, zOne);
1197  }
1198  if ((fRMin1 == fRMin2) && (fRMax1 == fRMax2))
1199  {
1200  return GetPointOnTubs(fRMin1, fRMax1, zOne, zTwo, totArea);
1201  }
1202  return GetPointOnCone(fRMin1, fRMax1, fRMin2, fRMax2, zOne, zTwo, totArea);
1203 }
1204 
1205 
1206 //
1207 // GetPointOnSurface
1208 //
1210 {
1211  double Area = 0, totArea = 0, Achose1 = 0, Achose2 = 0, phi, cosphi, sinphi, rRand;
1212  int i = 0;
1213  int numPlanes = fOriginalParameters->fNumZPlanes;
1214 
1215  phi = UUtils::Random(startPhi, endPhi);
1216  cosphi = std::cos(phi);
1217  sinphi = std::sin(phi);
1218 
1219  rRand = fOriginalParameters->Rmin[0] +
1221  * std::sqrt(UUtils::Random()));
1222 
1223  std::vector<double> areas; // (numPlanes+1);
1224  std::vector<UVector3> points; // (numPlanes-1);
1225 
1226  areas.push_back(UUtils::kPi * (UUtils::sqr(fOriginalParameters->Rmax[0])
1228 
1229  for (i = 0; i < numPlanes - 1; i++)
1230  {
1231  Area = (fOriginalParameters->Rmin[i] + fOriginalParameters->Rmin[i + 1])
1232  * std::sqrt(UUtils::sqr(fOriginalParameters->Rmin[i]
1233  - fOriginalParameters->Rmin[i + 1]) +
1235  - fOriginalParameters->fZValues[i]));
1236 
1237  Area += (fOriginalParameters->Rmax[i] + fOriginalParameters->Rmax[i + 1])
1238  * std::sqrt(UUtils::sqr(fOriginalParameters->Rmax[i]
1239  - fOriginalParameters->Rmax[i + 1]) +
1241  - fOriginalParameters->fZValues[i]));
1242 
1243  Area *= 0.5 * (endPhi - startPhi);
1244 
1245  if (startPhi == 0. && endPhi == 2 * UUtils::kPi)
1246  {
1247  Area += std::fabs(fOriginalParameters->fZValues[i + 1]
1248  - fOriginalParameters->fZValues[i]) *
1250  + fOriginalParameters->Rmax[i + 1]
1252  - fOriginalParameters->Rmin[i + 1]);
1253  }
1254  areas.push_back(Area);
1255  totArea += Area;
1256  }
1257 
1258  areas.push_back(UUtils::kPi * (UUtils::sqr(fOriginalParameters->Rmax[numPlanes - 1]) -
1259  UUtils::sqr(fOriginalParameters->Rmin[numPlanes - 1])));
1260 
1261  totArea += (areas[0] + areas[numPlanes]);
1262  double chose = UUtils::Random(0., totArea);
1263 
1264  if ((chose >= 0.) && (chose < areas[0]))
1265  {
1266  return UVector3(rRand * cosphi, rRand * sinphi,
1268  }
1269 
1270  for (i = 0; i < numPlanes - 1; i++)
1271  {
1272  Achose1 += areas[i];
1273  Achose2 = (Achose1 + areas[i + 1]);
1274  if (chose >= Achose1 && chose < Achose2)
1275  {
1278  fOriginalParameters->Rmin[i + 1],
1279  fOriginalParameters->Rmax[i + 1],
1281  fOriginalParameters->fZValues[i + 1], Area);
1282  }
1283  }
1284 
1285  rRand = fOriginalParameters->Rmin[numPlanes - 1] +
1286  ((fOriginalParameters->Rmax[numPlanes - 1] - fOriginalParameters->Rmin[numPlanes - 1])
1287  * std::sqrt(UUtils::Random()));
1288 
1289  return UVector3(rRand * cosphi, rRand * sinphi,
1290  fOriginalParameters->fZValues[numPlanes - 1]);
1291 
1292 }
1293 
1294 //
1295 // UPolyconeHistorical stuff
1296 //
1297 
1299  : fStartAngle(0.), fOpeningAngle(0.), fNumZPlanes(0),
1300  fZValues(0), Rmin(0), Rmax(0)
1301 {
1302 }
1303 
1305 {
1306 }
1307 
1310 {
1311  fStartAngle = source.fStartAngle;
1312  fOpeningAngle = source.fOpeningAngle;
1313  fNumZPlanes = source.fNumZPlanes;
1314 
1315  fZValues.resize(fNumZPlanes);
1316  Rmin.resize(fNumZPlanes);
1317  Rmax.resize(fNumZPlanes);
1318 
1319  for (int i = 0; i < fNumZPlanes; i++)
1320  {
1321  fZValues[i] = source.fZValues[i];
1322  Rmin[i] = source.Rmin[i];
1323  Rmax[i] = source.Rmax[i];
1324  }
1325 }
1326 
1329 {
1330  if (&right == this) return *this;
1331 
1332  fStartAngle = right.fStartAngle;
1333  fOpeningAngle = right.fOpeningAngle;
1334  fNumZPlanes = right.fNumZPlanes;
1335 
1336  fZValues.resize(fNumZPlanes);
1337  Rmin.resize(fNumZPlanes);
1338  Rmax.resize(fNumZPlanes);
1339 
1340  for (int i = 0; i < fNumZPlanes; i++)
1341  {
1342  fZValues[i] = right.fZValues[i];
1343  Rmin[i] = right.Rmin[i];
1344  Rmax[i] = right.Rmax[i];
1345  }
1346 
1347  return *this;
1348 }
1349 
1351 {
1352  return new UPolycone(*this);
1353 }
1354 //
1355 // Copy constructor
1356 //
1357 UPolycone::UPolycone(const UPolycone& source): VUSolid(source)
1358 {
1359  CopyStuff(source);
1360 }
1361 
1362 
1363 //
1364 // Assignment operator
1365 //
1367 {
1368  if (this == &source) return *this;
1369 
1370  //VUSolid::operator=( source );
1371 
1372  //delete [] corners;
1373 
1374  delete enclosingCylinder;
1375 
1376  CopyStuff(source);
1377 
1378  return *this;
1379 }
1380 //
1381 // CopyStuff
1382 //
1383 void UPolycone::CopyStuff(const UPolycone& source)
1384 {
1385  //
1386  // Simple stuff
1387  //
1388 
1389  startPhi = source.startPhi;
1390  endPhi = source.endPhi;
1391  phiIsOpen = source.phiIsOpen;
1392  fCubicVolume = source.fCubicVolume;
1393  fSurfaceArea = source.fSurfaceArea;
1394  fBox = source.fBox;
1395  //
1396  // The array of planes
1397  //
1399  //
1400  // Enclosing cylinder
1401  //
1403 }
1404 //
1405 // Get Entity Type
1406 //
1408 {
1409  return "Polycone";
1410 }
1411 //
1412 // Set Original Parameters
1413 //
1415 {
1416  int numPlanes = (int)numCorner;
1417  bool isConvertible = true;
1418  double Zmax = rz->Bmax();
1419  rz->StartWithZMin();
1420 
1421  // Prepare vectors for storage
1422  //
1423  std::vector<double> Z;
1424  std::vector<double> Rmin;
1425  std::vector<double> Rmax;
1426 
1427  int countPlanes = 1;
1428  int icurr = 0;
1429  int icurl = 0;
1430 
1431  // first plane Z=Z[0]
1432  //
1433  Z.push_back(corners[0].z);
1434  double Zprev = Z[0];
1435  if (Zprev == corners[1].z)
1436  {
1437  Rmin.push_back(corners[0].r);
1438  Rmax.push_back(corners[1].r);
1439  icurr = 1;
1440  }
1441  else if (Zprev == corners[numPlanes - 1].z)
1442  {
1443  Rmin.push_back(corners[numPlanes - 1].r);
1444  Rmax.push_back(corners[0].r);
1445  icurl = numPlanes - 1;
1446  }
1447  else
1448  {
1449  Rmin.push_back(corners[0].r);
1450  Rmax.push_back(corners[0].r);
1451  }
1452 
1453  // next planes until last
1454  //
1455  int inextr = 0, inextl = 0;
1456  for (int i = 0; i < numPlanes - 2; i++)
1457  {
1458  inextr = 1 + icurr;
1459  inextl = (icurl <= 0) ? numPlanes - 1 : icurl - 1;
1460 
1461  if ((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax))
1462  {
1463  break;
1464  }
1465 
1466  double Zleft = corners[inextl].z;
1467  double Zright = corners[inextr].z;
1468  if (Zright > Zleft)
1469  {
1470  Z.push_back(Zleft);
1471  countPlanes++;
1472  double difZr = corners[inextr].z - corners[icurr].z;
1473  double difZl = corners[inextl].z - corners[icurl].z;
1474 
1475  if (std::fabs(difZl) < frTolerance)
1476  {
1477  if (std::fabs(difZr) < frTolerance)
1478  {
1479  Rmin.push_back(corners[inextl].r);
1480  Rmax.push_back(corners[icurr].r);
1481  }
1482  else
1483  {
1484  Rmin.push_back(corners[inextl].r);
1485  Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1486  *(corners[inextr].r - corners[icurr].r));
1487  }
1488  }
1489  else if (difZl >= frTolerance)
1490  {
1491  if (std::fabs(difZr) < frTolerance)
1492  {
1493  Rmin.push_back(corners[icurl].r);
1494  Rmax.push_back(corners[icurr].r);
1495  }
1496  else
1497  {
1498  Rmin.push_back(corners[icurl].r);
1499  Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
1500  *(corners[inextr].r - corners[icurr].r));
1501  }
1502  }
1503  else
1504  {
1505  isConvertible = false;
1506  break;
1507  }
1508  icurl = (icurl == 0) ? numPlanes - 1 : icurl - 1;
1509  }
1510  else if (std::fabs(Zright - Zleft) < frTolerance) // Zright=Zleft
1511  {
1512  Z.push_back(Zleft);
1513  countPlanes++;
1514  icurr++;
1515 
1516  icurl = (icurl == 0) ? numPlanes - 1 : icurl - 1;
1517 
1518  Rmin.push_back(corners[inextl].r);
1519  Rmax.push_back(corners[inextr].r);
1520  }
1521  else // Zright<Zleft
1522  {
1523  Z.push_back(Zright);
1524  countPlanes++;
1525 
1526  double difZr = corners[inextr].z - corners[icurr].z;
1527  double difZl = corners[inextl].z - corners[icurl].z;
1528  if (std::fabs(difZr) < frTolerance)
1529  {
1530  if (std::fabs(difZl) < frTolerance)
1531  {
1532  Rmax.push_back(corners[inextr].r);
1533  Rmin.push_back (corners[icurr].r);
1534  }
1535  else
1536  {
1537  Rmin.push_back (corners[icurl].r + (Zright-corners[icurl].z)/difZl
1538  * (corners[inextl].r - corners[icurl].r));
1539  Rmax.push_back(corners[inextr].r);
1540  }
1541  icurr++;
1542  } // plate
1543  else if (difZr >= frTolerance)
1544  {
1545  if (std::fabs(difZl) < frTolerance)
1546  {
1547  Rmax.push_back(corners[inextr].r);
1548  Rmin.push_back(corners[icurr].r);
1549  }
1550  else
1551  {
1552  Rmax.push_back(corners[inextr].r);
1553  Rmin.push_back(corners[icurl].r + (Zright - corners[icurl].z) / difZl
1554  * (corners[inextl].r - corners[icurl].r));
1555  }
1556  icurr++;
1557  }
1558  else
1559  {
1560  isConvertible = false;
1561  break;
1562  }
1563  }
1564  } // end for loop
1565 
1566  // last plane Z=Zmax
1567  //
1568  Z.push_back(Zmax);
1569  countPlanes++;
1570  inextr = 1 + icurr;
1571  inextl = (icurl <= 0) ? numPlanes - 1 : icurl - 1;
1572 
1573  if (corners[inextr].z == corners[inextl].z)
1574  {
1575  Rmax.push_back(corners[inextr].r);
1576  Rmin.push_back(corners[inextl].r);
1577  }
1578  else
1579  {
1580  Rmax.push_back(corners[inextr].r);
1581  Rmin.push_back(corners[inextl].r);
1582  }
1583 
1584  // Set original parameters Rmin,Rmax,Z
1585  //
1586  if (isConvertible)
1587  {
1589  fOriginalParameters->fZValues.resize(numPlanes);
1590  fOriginalParameters->Rmin.resize(numPlanes);
1591  fOriginalParameters->Rmax.resize(numPlanes);
1592 
1593  for (int j = 0; j < countPlanes; j++)
1594  {
1595  fOriginalParameters->fZValues[j] = Z[j];
1596  fOriginalParameters->Rmax[j] = Rmax[j];
1597  fOriginalParameters->Rmin[j] = Rmin[j];
1598  }
1601  fOriginalParameters->fNumZPlanes = countPlanes;
1602 
1603  }
1604  else // Set parameters(r,z) with Rmin==0 as convention
1605  {
1606  std::ostringstream message;
1607  message << "Polycone " << GetName() << std::endl
1608  << "cannot be converted to Polycone with (Rmin,Rmaz,Z) parameters!";
1609  UUtils::Exception("UPolycone::SetOriginalParameters()", "GeomSolids0002",
1610  Warning, 1, "can not convert");
1611 
1613 
1614  fOriginalParameters->fZValues.resize(numPlanes);
1615  fOriginalParameters->Rmin.resize(numPlanes);
1616  fOriginalParameters->Rmax.resize(numPlanes);
1617 
1618  for (int j = 0; j < numPlanes; j++)
1619  {
1621  fOriginalParameters->Rmax[j] = corners[j].r;
1622  fOriginalParameters->Rmin[j] = 0.0;
1623  }
1626  fOriginalParameters->fNumZPlanes = numPlanes;
1627  }
1628  return isConvertible;
1629 }
1630 //
1631 // Reset
1632 //
1634 {
1635  //
1636  // Clear old setup
1637  //
1638  delete enclosingCylinder;
1639 
1640  fCubicVolume = 0;
1641  fSurfaceArea = 0;
1642  double phiStart=fOriginalParameters->fStartAngle;
1643  double* Z, *R1, *R2;
1644  int num = fOriginalParameters->fNumZPlanes;
1645  Z = new double[num];
1646  R1 = new double[num];
1647  R2 = new double[num];
1648  for (int i = 0; i < num; i++)
1649  {
1650  Z[i] = fOriginalParameters->fZValues[i];
1651  R1[i] = fOriginalParameters->Rmin[i];
1652  R2[i] = fOriginalParameters->Rmax[i];
1653  }
1654 
1655  Init(phiStart, phiStart+ fOriginalParameters->fOpeningAngle, num, Z, R1, R2);
1656  delete [] R1;
1657  delete [] Z;
1658  delete [] R2;
1659 }
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:539
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:1082
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:1383
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:1189
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:1328
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: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:53
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:49
UVector3 GetPointOnSurface() const
Definition: UPolycone.cc:1209
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 GetDz() const
UPolyconeSideRZ * corners
Definition: UPolycone.hh:220
UGeometryType GetEntityType() const
Definition: UPolycone.cc:1407
void Reset()
Definition: UPolycone.cc:1633
double Bmin() const
EnumInside
Definition: VUSolid.hh:23
const G4int n
double Mag() const
Definition: UVector3.cc:48
double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const
Definition: UBox.cc:117
double SurfaceArea()
Definition: UPolycone.cc:911
T sqr(const T &x)
Definition: UUtils.hh:137
static const G4double minSafety
double Capacity()
Definition: UPolycone.cc:894
double GetZHalfLength() const
Definition: UBox.hh:121
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:245
double GetRmax1() const
static const double kPi
Definition: UUtils.hh:48
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:1350
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:1141
double Bmax() const
UVector3 GetPointOnCone(double fRmin1, double fRmax1, double fRmin2, double fRmax2, double zOne, double zTwo, double &totArea) const
Definition: UPolycone.cc:978
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:1366
double z
Definition: UVector3.hh:138
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:886
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
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