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