Geant4  10.02.p02
G4Polycone.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4Polycone.cc 93559 2015-10-26 14:14:49Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4Polycone.cc
35 //
36 // Implementation of a CSG polycone
37 //
38 // --------------------------------------------------------------------
39 
40 #include "G4Polycone.hh"
41 
42 #if !defined(G4GEOM_USE_UPOLYCONE)
43 
44 #include "G4PolyconeSide.hh"
45 #include "G4PolyPhiFace.hh"
46 
47 #include "Randomize.hh"
48 
49 #include "G4EnclosingCylinder.hh"
50 #include "G4ReduciblePolygon.hh"
51 #include "G4VPVParameterisation.hh"
52 
53 using namespace CLHEP;
54 
55 //
56 // Constructor (GEANT3 style parameters)
57 //
59  G4double phiStart,
60  G4double phiTotal,
61  G4int numZPlanes,
62  const G4double zPlane[],
63  const G4double rInner[],
64  const G4double rOuter[] )
65  : G4VCSGfaceted( name )
66 {
67  //
68  // Some historical ugliness
69  //
71 
72  original_parameters->Start_angle = phiStart;
74  original_parameters->Num_z_planes = numZPlanes;
75  original_parameters->Z_values = new G4double[numZPlanes];
76  original_parameters->Rmin = new G4double[numZPlanes];
77  original_parameters->Rmax = new G4double[numZPlanes];
78 
79  G4int i;
80  for (i=0; i<numZPlanes; i++)
81  {
82  if(rInner[i]>rOuter[i])
83  {
84  DumpInfo();
85  std::ostringstream message;
86  message << "Cannot create a Polycone with rInner > rOuter for the same Z"
87  << G4endl
88  << " rInner > rOuter for the same Z !" << G4endl
89  << " rMin[" << i << "] = " << rInner[i]
90  << " -- rMax[" << i << "] = " << rOuter[i];
91  G4Exception("G4Polycone::G4Polycone()", "GeomSolids0002",
92  FatalErrorInArgument, message);
93  }
94  if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
95  {
96  if( (rInner[i] > rOuter[i+1])
97  ||(rInner[i+1] > rOuter[i]) )
98  {
99  DumpInfo();
100  std::ostringstream message;
101  message << "Cannot create a Polycone with no contiguous segments."
102  << G4endl
103  << " Segments are not contiguous !" << G4endl
104  << " rMin[" << i << "] = " << rInner[i]
105  << " -- rMax[" << i+1 << "] = " << rOuter[i+1] << G4endl
106  << " rMin[" << i+1 << "] = " << rInner[i+1]
107  << " -- rMax[" << i << "] = " << rOuter[i];
108  G4Exception("G4Polycone::G4Polycone()", "GeomSolids0002",
109  FatalErrorInArgument, message);
110  }
111  }
112  original_parameters->Z_values[i] = zPlane[i];
113  original_parameters->Rmin[i] = rInner[i];
114  original_parameters->Rmax[i] = rOuter[i];
115  }
116 
117  //
118  // Build RZ polygon using special PCON/PGON GEANT3 constructor
119  //
120  G4ReduciblePolygon *rz =
121  new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
122 
123  //
124  // Do the real work
125  //
126  Create( phiStart, phiTotal, rz );
127 
128  delete rz;
129 }
130 
131 
132 //
133 // Constructor (generic parameters)
134 //
136  G4double phiStart,
137  G4double phiTotal,
138  G4int numRZ,
139  const G4double r[],
140  const G4double z[] )
141  : G4VCSGfaceted( name )
142 {
143 
144  G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
145 
146  Create( phiStart, phiTotal, rz );
147 
148  // Set original_parameters struct for consistency
149  //
150 
151  G4bool convertible=SetOriginalParameters(rz);
152 
153  if(!convertible)
154  {
155  std::ostringstream message;
156  message << "Polycone " << GetName() << "cannot be converted" << G4endl
157  << "to Polycone with (Rmin,Rmaz,Z) parameters!";
158  G4Exception("G4Polycone::G4Polycone()", "GeomSolids0002",
159  FatalException, message, "Use G4GenericPolycone instead!");
160  }
161  else
162  {
163  G4cout << "INFO: Converting polycone " << GetName() << G4endl
164  << "to optimized polycone with (Rmin,Rmaz,Z) parameters !"
165  << G4endl;
166  }
167  delete rz;
168 }
169 
170 
171 // Create
172 //
173 // Generic create routine, called by each constructor after
174 // conversion of arguments
175 //
177  G4double phiTotal,
178  G4ReduciblePolygon *rz )
179 {
180  //
181  // Perform checks of rz values
182  //
183  if (rz->Amin() < 0.0)
184  {
185  std::ostringstream message;
186  message << "Illegal input parameters - " << GetName() << G4endl
187  << " All R values must be >= 0 !";
188  G4Exception("G4Polycone::Create()", "GeomSolids0002",
189  FatalErrorInArgument, message);
190  }
191 
192  G4double rzArea = rz->Area();
193  if (rzArea < -kCarTolerance)
194  rz->ReverseOrder();
195 
196  else if (rzArea < -kCarTolerance)
197  {
198  std::ostringstream message;
199  message << "Illegal input parameters - " << GetName() << G4endl
200  << " R/Z cross section is zero or near zero: " << rzArea;
201  G4Exception("G4Polycone::Create()", "GeomSolids0002",
202  FatalErrorInArgument, message);
203  }
204 
206  || (!rz->RemoveRedundantVertices( kCarTolerance )) )
207  {
208  std::ostringstream message;
209  message << "Illegal input parameters - " << GetName() << G4endl
210  << " Too few unique R/Z values !";
211  G4Exception("G4Polycone::Create()", "GeomSolids0002",
212  FatalErrorInArgument, message);
213  }
214 
215  if (rz->CrossesItself(1/kInfinity))
216  {
217  std::ostringstream message;
218  message << "Illegal input parameters - " << GetName() << G4endl
219  << " R/Z segments cross !";
220  G4Exception("G4Polycone::Create()", "GeomSolids0002",
221  FatalErrorInArgument, message);
222  }
223 
224  numCorner = rz->NumVertices();
225 
226  //
227  // Phi opening? Account for some possible roundoff, and interpret
228  // nonsense value as representing no phi opening
229  //
230  if (phiTotal <= 0 || phiTotal > twopi-1E-10)
231  {
232  phiIsOpen = false;
233  startPhi = 0;
234  endPhi = twopi;
235  }
236  else
237  {
238  phiIsOpen = true;
239 
240  //
241  // Convert phi into our convention
242  //
243  startPhi = phiStart;
244  while( startPhi < 0 ) // Loop checking, 13.08.2015, G.Cosmo
245  startPhi += twopi;
246 
247  endPhi = phiStart+phiTotal;
248  while( endPhi < startPhi ) // Loop checking, 13.08.2015, G.Cosmo
249  endPhi += twopi;
250  }
251 
252  //
253  // Allocate corner array.
254  //
256 
257  //
258  // Copy corners
259  //
260  G4ReduciblePolygonIterator iterRZ(rz);
261 
262  G4PolyconeSideRZ *next = corners;
263  iterRZ.Begin();
264  do // Loop checking, 13.08.2015, G.Cosmo
265  {
266  next->r = iterRZ.GetA();
267  next->z = iterRZ.GetB();
268  } while( ++next, iterRZ.Next() );
269 
270  //
271  // Allocate face pointer array
272  //
274  faces = new G4VCSGface*[numFace];
275 
276  //
277  // Construct conical faces
278  //
279  // But! Don't construct a face if both points are at zero radius!
280  //
281  G4PolyconeSideRZ *corner = corners,
282  *prev = corners + numCorner-1,
283  *nextNext;
284  G4VCSGface **face = faces;
285  do // Loop checking, 13.08.2015, G.Cosmo
286  {
287  next = corner+1;
288  if (next >= corners+numCorner) next = corners;
289  nextNext = next+1;
290  if (nextNext >= corners+numCorner) nextNext = corners;
291 
292  if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
293 
294  //
295  // We must decide here if we can dare declare one of our faces
296  // as having a "valid" normal (i.e. allBehind = true). This
297  // is never possible if the face faces "inward" in r.
298  //
299  G4bool allBehind;
300  if (corner->z > next->z)
301  {
302  allBehind = false;
303  }
304  else
305  {
306  //
307  // Otherwise, it is only true if the line passing
308  // through the two points of the segment do not
309  // split the r/z cross section
310  //
311  allBehind = !rz->BisectedBy( corner->r, corner->z,
312  next->r, next->z, kCarTolerance );
313  }
314 
315  *face++ = new G4PolyconeSide( prev, corner, next, nextNext,
316  startPhi, endPhi-startPhi, phiIsOpen, allBehind );
317  } while( prev=corner, corner=next, corner > corners );
318 
319  if (phiIsOpen)
320  {
321  //
322  // Construct phi open edges
323  //
324  *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi );
325  *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi );
326  }
327 
328  //
329  // We might have dropped a face or two: recalculate numFace
330  //
331  numFace = face-faces;
332 
333  //
334  // Make enclosingCylinder
335  //
337  new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
338 }
339 
340 
341 //
342 // Fake default constructor - sets only member data and allocates memory
343 // for usage restricted to object persistency.
344 //
346  : G4VCSGfaceted(a), startPhi(0.), endPhi(0.), phiIsOpen(false),
347  numCorner(0), corners(0),original_parameters(0),
348  enclosingCylinder(0)
349 {
350 }
351 
352 
353 //
354 // Destructor
355 //
357 {
358  delete [] corners;
359  delete original_parameters;
360  delete enclosingCylinder;
361 }
362 
363 
364 //
365 // Copy constructor
366 //
368  : G4VCSGfaceted( source )
369 {
370  CopyStuff( source );
371 }
372 
373 
374 //
375 // Assignment operator
376 //
378 {
379  if (this == &source) return *this;
380 
381  G4VCSGfaceted::operator=( source );
382 
383  delete [] corners;
385 
386  delete enclosingCylinder;
387 
388  CopyStuff( source );
389 
390  return *this;
391 }
392 
393 
394 //
395 // CopyStuff
396 //
397 void G4Polycone::CopyStuff( const G4Polycone &source )
398 {
399  //
400  // Simple stuff
401  //
402  startPhi = source.startPhi;
403  endPhi = source.endPhi;
404  phiIsOpen = source.phiIsOpen;
405  numCorner = source.numCorner;
406 
407  //
408  // The corner array
409  //
411 
412  G4PolyconeSideRZ *corn = corners,
413  *sourceCorn = source.corners;
414  do // Loop checking, 13.08.2015, G.Cosmo
415  {
416  *corn = *sourceCorn;
417  } while( ++sourceCorn, ++corn < corners+numCorner );
418 
419  //
420  // Original parameters
421  //
422  if (source.original_parameters)
423  {
426  }
427 
428  //
429  // Enclosing cylinder
430  //
432 
433  fRebuildPolyhedron = false;
434  fpPolyhedron = 0;
435 }
436 
437 
438 //
439 // Reset
440 //
442 {
443  //
444  // Clear old setup
445  //
447  delete [] corners;
448  delete enclosingCylinder;
449 
450  //
451  // Rebuild polycone
452  //
453  G4ReduciblePolygon *rz =
460  delete rz;
461 
462  return 0;
463 }
464 
465 
466 //
467 // Inside
468 //
469 // This is an override of G4VCSGfaceted::Inside, created in order
470 // to speed things up by first checking with G4EnclosingCylinder.
471 //
473 {
474  //
475  // Quick test
476  //
477  if (enclosingCylinder->MustBeOutside(p)) return kOutside;
478 
479  //
480  // Long answer
481  //
482  return G4VCSGfaceted::Inside(p);
483 }
484 
485 
486 //
487 // DistanceToIn
488 //
489 // This is an override of G4VCSGfaceted::Inside, created in order
490 // to speed things up by first checking with G4EnclosingCylinder.
491 //
493  const G4ThreeVector &v ) const
494 {
495  //
496  // Quick test
497  //
498  if (enclosingCylinder->ShouldMiss(p,v))
499  return kInfinity;
500 
501  //
502  // Long answer
503  //
504  return G4VCSGfaceted::DistanceToIn( p, v );
505 }
506 
507 
508 //
509 // DistanceToIn
510 //
512 {
513  return G4VCSGfaceted::DistanceToIn(p);
514 }
515 
516 
517 //
518 // ComputeDimensions
519 //
521  const G4int n,
522  const G4VPhysicalVolume* pRep )
523 {
524  p->ComputeDimensions(*this,n,pRep);
525 }
526 
527 //
528 // GetEntityType
529 //
531 {
532  return G4String("G4Polycone");
533 }
534 
535 //
536 // Make a clone of the object
537 //
539 {
540  return new G4Polycone(*this);
541 }
542 
543 //
544 // Stream object contents to an output stream
545 //
546 std::ostream& G4Polycone::StreamInfo( std::ostream& os ) const
547 {
548  G4int oldprc = os.precision(16);
549  os << "-----------------------------------------------------------\n"
550  << " *** Dump for solid - " << GetName() << " ***\n"
551  << " ===================================================\n"
552  << " Solid type: G4Polycone\n"
553  << " Parameters: \n"
554  << " starting phi angle : " << startPhi/degree << " degrees \n"
555  << " ending phi angle : " << endPhi/degree << " degrees \n";
556  G4int i=0;
557 
559  os << " number of Z planes: " << numPlanes << "\n"
560  << " Z values: \n";
561  for (i=0; i<numPlanes; i++)
562  {
563  os << " Z plane " << i << ": "
564  << original_parameters->Z_values[i] << "\n";
565  }
566  os << " Tangent distances to inner surface (Rmin): \n";
567  for (i=0; i<numPlanes; i++)
568  {
569  os << " Z plane " << i << ": "
570  << original_parameters->Rmin[i] << "\n";
571  }
572  os << " Tangent distances to outer surface (Rmax): \n";
573  for (i=0; i<numPlanes; i++)
574  {
575  os << " Z plane " << i << ": "
576  << original_parameters->Rmax[i] << "\n";
577  }
578 
579  os << " number of RZ points: " << numCorner << "\n"
580  << " RZ values (corners): \n";
581  for (i=0; i<numCorner; i++)
582  {
583  os << " "
584  << corners[i].r << ", " << corners[i].z << "\n";
585  }
586  os << "-----------------------------------------------------------\n";
587  os.precision(oldprc);
588 
589  return os;
590 }
591 
592 
593 //
594 // GetPointOnCone
595 //
596 // Auxiliary method for Get Point On Surface
597 //
599  G4double fRmin2, G4double fRmax2,
600  G4double zOne, G4double zTwo,
601  G4double& totArea) const
602 {
603  // declare working variables
604  //
605  G4double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
606  G4double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
607  G4double fDz=(zTwo-zOne)/2., afDz=std::fabs(fDz);
608  G4ThreeVector point, offset=G4ThreeVector(0.,0.,0.5*(zTwo+zOne));
609  fDPhi = endPhi - startPhi;
610  rone = (fRmax1-fRmax2)/(2.*fDz);
611  rtwo = (fRmin1-fRmin2)/(2.*fDz);
612  if(fRmax1==fRmax2){qone=0.;}
613  else{
614  qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
615  }
616  if(fRmin1==fRmin2){qtwo=0.;}
617  else{
618  qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
619  }
620  Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*(sqr(fRmin1-fRmin2)+sqr(zTwo-zOne));
621  Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*(sqr(fRmax1-fRmax2)+sqr(zTwo-zOne));
622  Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
623  totArea = Aone+Atwo+2.*Afive;
624 
625  phi = RandFlat::shoot(startPhi,endPhi);
626  cosu = std::cos(phi);
627  sinu = std::sin(phi);
628 
629 
630  if( (startPhi == 0) && (endPhi == twopi) ) { Afive = 0; }
631  chose = RandFlat::shoot(0.,Aone+Atwo+2.*Afive);
632  if( (chose >= 0) && (chose < Aone) )
633  {
634  if(fRmax1 != fRmax2)
635  {
636  zRand = RandFlat::shoot(-1.*afDz,afDz);
637  point = G4ThreeVector (rone*cosu*(qone-zRand),
638  rone*sinu*(qone-zRand), zRand);
639  }
640  else
641  {
642  point = G4ThreeVector(fRmax1*cosu, fRmax1*sinu,
643  RandFlat::shoot(-1.*afDz,afDz));
644 
645  }
646  }
647  else if(chose >= Aone && chose < Aone + Atwo)
648  {
649  if(fRmin1 != fRmin2)
650  {
651  zRand = RandFlat::shoot(-1.*afDz,afDz);
652  point = G4ThreeVector (rtwo*cosu*(qtwo-zRand),
653  rtwo*sinu*(qtwo-zRand), zRand);
654 
655  }
656  else
657  {
658  point = G4ThreeVector(fRmin1*cosu, fRmin1*sinu,
659  RandFlat::shoot(-1.*afDz,afDz));
660  }
661  }
662  else if( (chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive) )
663  {
664  zRand = RandFlat::shoot(-1.*afDz,afDz);
665  rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
666  rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
667  rRand1 = std::sqrt(RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
668  point = G4ThreeVector (rRand1*std::cos(startPhi),
669  rRand1*std::sin(startPhi), zRand);
670  }
671  else
672  {
673  zRand = RandFlat::shoot(-1.*afDz,afDz);
674  rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
675  rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
676  rRand1 = std::sqrt(RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
677  point = G4ThreeVector (rRand1*std::cos(endPhi),
678  rRand1*std::sin(endPhi), zRand);
679 
680  }
681 
682  return point+offset;
683 }
684 
685 
686 //
687 // GetPointOnTubs
688 //
689 // Auxiliary method for GetPoint On Surface
690 //
692  G4double zOne, G4double zTwo,
693  G4double& totArea) const
694 {
695  G4double xRand,yRand,zRand,phi,cosphi,sinphi,chose,
696  aOne,aTwo,aFou,rRand,fDz,fSPhi,fDPhi;
697  fDz = std::fabs(0.5*(zTwo-zOne));
698  fSPhi = startPhi;
699  fDPhi = endPhi-startPhi;
700 
701  aOne = 2.*fDz*fDPhi*fRMax;
702  aTwo = 2.*fDz*fDPhi*fRMin;
703  aFou = 2.*fDz*(fRMax-fRMin);
704  totArea = aOne+aTwo+2.*aFou;
705  phi = RandFlat::shoot(startPhi,endPhi);
706  cosphi = std::cos(phi);
707  sinphi = std::sin(phi);
708  rRand = fRMin + (fRMax-fRMin)*std::sqrt(RandFlat::shoot());
709 
710  if(startPhi == 0 && endPhi == twopi)
711  aFou = 0;
712 
713  chose = RandFlat::shoot(0.,aOne+aTwo+2.*aFou);
714  if( (chose >= 0) && (chose < aOne) )
715  {
716  xRand = fRMax*cosphi;
717  yRand = fRMax*sinphi;
718  zRand = RandFlat::shoot(-1.*fDz,fDz);
719  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
720  }
721  else if( (chose >= aOne) && (chose < aOne + aTwo) )
722  {
723  xRand = fRMin*cosphi;
724  yRand = fRMin*sinphi;
725  zRand = RandFlat::shoot(-1.*fDz,fDz);
726  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
727  }
728  else if( (chose >= aOne+aTwo) && (chose <aOne+aTwo+aFou) )
729  {
730  xRand = rRand*std::cos(fSPhi+fDPhi);
731  yRand = rRand*std::sin(fSPhi+fDPhi);
732  zRand = RandFlat::shoot(-1.*fDz,fDz);
733  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
734  }
735 
736  // else
737 
738  xRand = rRand*std::cos(fSPhi+fDPhi);
739  yRand = rRand*std::sin(fSPhi+fDPhi);
740  zRand = RandFlat::shoot(-1.*fDz,fDz);
741  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
742 }
743 
744 
745 //
746 // GetPointOnRing
747 //
748 // Auxiliary method for GetPoint On Surface
749 //
751  G4double fRMin2,G4double fRMax2,
752  G4double zOne) const
753 {
754  G4double xRand,yRand,phi,cosphi,sinphi,rRand1,rRand2,A1,Atot,rCh;
756  cosphi = std::cos(phi);
757  sinphi = std::sin(phi);
758 
759  if(fRMin1==fRMin2)
760  {
761  rRand1 = fRMin1; A1=0.;
762  }
763  else
764  {
765  rRand1 = RandFlat::shoot(fRMin1,fRMin2);
766  A1=std::fabs(fRMin2*fRMin2-fRMin1*fRMin1);
767  }
768  if(fRMax1==fRMax2)
769  {
770  rRand2=fRMax1; Atot=A1;
771  }
772  else
773  {
774  rRand2 = RandFlat::shoot(fRMax1,fRMax2);
775  Atot = A1+std::fabs(fRMax2*fRMax2-fRMax1*fRMax1);
776  }
777  rCh = RandFlat::shoot(0.,Atot);
778 
779  if(rCh>A1) { rRand1=rRand2; }
780 
781  xRand = rRand1*cosphi;
782  yRand = rRand1*sinphi;
783 
784  return G4ThreeVector(xRand, yRand, zOne);
785 }
786 
787 
788 //
789 // GetPointOnCut
790 //
791 // Auxiliary method for Get Point On Surface
792 //
794  G4double fRMin2, G4double fRMax2,
795  G4double zOne, G4double zTwo,
796  G4double& totArea) const
797 { if(zOne==zTwo)
798  {
799  return GetPointOnRing(fRMin1, fRMax1,fRMin2,fRMax2,zOne);
800  }
801  if( (fRMin1 == fRMin2) && (fRMax1 == fRMax2) )
802  {
803  return GetPointOnTubs(fRMin1, fRMax1,zOne,zTwo,totArea);
804  }
805  return GetPointOnCone(fRMin1,fRMax1,fRMin2,fRMax2,zOne,zTwo,totArea);
806 }
807 
808 
809 //
810 // GetPointOnSurface
811 //
813 {
814  G4double Area=0,totArea=0,Achose1=0,Achose2=0,phi,cosphi,sinphi,rRand;
815  G4int i=0;
817 
819  cosphi = std::cos(phi);
820  sinphi = std::sin(phi);
821 
822  rRand = original_parameters->Rmin[0] +
824  * std::sqrt(RandFlat::shoot()) );
825 
826  std::vector<G4double> areas; // (numPlanes+1);
827  std::vector<G4ThreeVector> points; // (numPlanes-1);
828 
829  areas.push_back(pi*(sqr(original_parameters->Rmax[0])
830  -sqr(original_parameters->Rmin[0])));
831 
832  for(i=0; i<numPlanes-1; i++)
833  {
835  * std::sqrt(sqr(original_parameters->Rmin[i]
836  -original_parameters->Rmin[i+1])+
839 
841  * std::sqrt(sqr(original_parameters->Rmax[i]
842  -original_parameters->Rmax[i+1])+
845 
846  Area *= 0.5*(endPhi-startPhi);
847 
848  if(startPhi==0.&& endPhi == twopi)
849  {
850  Area += std::fabs(original_parameters->Z_values[i+1]
855  -original_parameters->Rmin[i+1]);
856  }
857  areas.push_back(Area);
858  totArea += Area;
859  }
860 
861  areas.push_back(pi*(sqr(original_parameters->Rmax[numPlanes-1])-
862  sqr(original_parameters->Rmin[numPlanes-1])));
863 
864  totArea += (areas[0]+areas[numPlanes]);
865  G4double chose = RandFlat::shoot(0.,totArea);
866 
867  if( (chose>=0.) && (chose<areas[0]) )
868  {
869  return G4ThreeVector(rRand*cosphi, rRand*sinphi,
871  }
872 
873  for (i=0; i<numPlanes-1; i++)
874  {
875  Achose1 += areas[i];
876  Achose2 = (Achose1+areas[i+1]);
877  if(chose>=Achose1 && chose<Achose2)
878  {
884  original_parameters->Z_values[i+1], Area);
885  }
886  }
887 
888  rRand = original_parameters->Rmin[numPlanes-1] +
889  ( (original_parameters->Rmax[numPlanes-1]-original_parameters->Rmin[numPlanes-1])
890  * std::sqrt(RandFlat::shoot()) );
891 
892  return G4ThreeVector(rRand*cosphi,rRand*sinphi,
893  original_parameters->Z_values[numPlanes-1]);
894 
895 }
896 
897 //
898 // CreatePolyhedron
899 //
901 {
902  //
903  // This has to be fixed in visualization. Fake it for the moment.
904  //
905 
912 }
913 
915 {
916  G4int numPlanes = (G4int)numCorner;
917  G4bool isConvertible=true;
918  G4double Zmax=rz->Bmax();
919  rz->StartWithZMin();
920 
921  // Prepare vectors for storage
922  //
923  std::vector<G4double> Z;
924  std::vector<G4double> Rmin;
925  std::vector<G4double> Rmax;
926 
927  G4int countPlanes=1;
928  G4int icurr=0;
929  G4int icurl=0;
930 
931  // first plane Z=Z[0]
932  //
933  Z.push_back(corners[0].z);
934  G4double Zprev=Z[0];
935  if (Zprev == corners[1].z)
936  {
937  Rmin.push_back(corners[0].r);
938  Rmax.push_back (corners[1].r);icurr=1;
939  }
940  else if (Zprev == corners[numPlanes-1].z)
941  {
942  Rmin.push_back(corners[numPlanes-1].r);
943  Rmax.push_back (corners[0].r);
944  icurl=numPlanes-1;
945  }
946  else
947  {
948  Rmin.push_back(corners[0].r);
949  Rmax.push_back (corners[0].r);
950  }
951 
952  // next planes until last
953  //
954  G4int inextr=0, inextl=0;
955  for (G4int i=0; i < numPlanes-2; i++)
956  {
957  inextr=1+icurr;
958  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
959 
960  if((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax)) { break; }
961 
962  G4double Zleft = corners[inextl].z;
963  G4double Zright = corners[inextr].z;
964  if(Zright > Zleft) // Next plane will be Zleft
965  {
966  Z.push_back(Zleft);
967  countPlanes++;
968  G4double difZr=corners[inextr].z - corners[icurr].z;
969  G4double difZl=corners[inextl].z - corners[icurl].z;
970 
971  if(std::fabs(difZl) < kCarTolerance)
972  {
973  if(std::fabs(difZr) < kCarTolerance)
974  {
975  Rmin.push_back(corners[inextl].r);
976  Rmax.push_back(corners[icurr].r);
977  }
978  else
979  {
980  Rmin.push_back(corners[inextl].r);
981  Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
982  *(corners[inextr].r - corners[icurr].r));
983  }
984  }
985  else if (difZl >= kCarTolerance)
986  {
987  if(std::fabs(difZr) < kCarTolerance)
988  {
989  Rmin.push_back(corners[icurl].r);
990  Rmax.push_back(corners[icurr].r);
991  }
992  else
993  {
994  Rmin.push_back(corners[icurl].r);
995  Rmax.push_back(corners[icurr].r + (Zleft-corners[icurr].z)/difZr
996  *(corners[inextr].r - corners[icurr].r));
997  }
998  }
999  else
1000  {
1001  isConvertible=false; break;
1002  }
1003  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1004  }
1005  else if(std::fabs(Zright-Zleft)<kCarTolerance) // Zright=Zleft
1006  {
1007  Z.push_back(Zleft);
1008  countPlanes++;
1009  icurr++;
1010 
1011  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1012 
1013  Rmin.push_back(corners[inextl].r);
1014  Rmax.push_back(corners[inextr].r);
1015  }
1016  else // Zright<Zleft
1017  {
1018  Z.push_back(Zright);
1019  countPlanes++;
1020 
1021  G4double difZr=corners[inextr].z - corners[icurr].z;
1022  G4double difZl=corners[inextl].z - corners[icurl].z;
1023  if(std::fabs(difZr) < kCarTolerance)
1024  {
1025  if(std::fabs(difZl) < kCarTolerance)
1026  {
1027  Rmax.push_back(corners[inextr].r);
1028  Rmin.push_back(corners[icurr].r);
1029  }
1030  else
1031  {
1032  Rmin.push_back(corners[icurl].r + (Zright-corners[icurl].z)/difZl
1033  *(corners[inextl].r - corners[icurl].r));
1034  Rmax.push_back(corners[inextr].r);
1035  }
1036  icurr++;
1037  } // plate
1038  else if (difZr >= kCarTolerance)
1039  {
1040  if(std::fabs(difZl) < kCarTolerance)
1041  {
1042  Rmax.push_back(corners[inextr].r);
1043  Rmin.push_back (corners[icurr].r);
1044  }
1045  else
1046  {
1047  Rmax.push_back(corners[inextr].r);
1048  Rmin.push_back (corners[icurl].r+(Zright-corners[icurl].z)/difZl
1049  * (corners[inextl].r - corners[icurl].r));
1050  }
1051  icurr++;
1052  }
1053  else
1054  {
1055  isConvertible=false; break;
1056  }
1057  }
1058  } // end for loop
1059 
1060  // last plane Z=Zmax
1061  //
1062  Z.push_back(Zmax);
1063  countPlanes++;
1064  inextr=1+icurr;
1065  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1066 
1067  Rmax.push_back(corners[inextr].r);
1068  Rmin.push_back(corners[inextl].r);
1069 
1070  // Set original parameters Rmin,Rmax,Z
1071  //
1072  if(isConvertible)
1073  {
1075  original_parameters->Z_values = new G4double[countPlanes];
1076  original_parameters->Rmin = new G4double[countPlanes];
1077  original_parameters->Rmax = new G4double[countPlanes];
1078 
1079  for(G4int j=0; j < countPlanes; j++)
1080  {
1081  original_parameters->Z_values[j] = Z[j];
1082  original_parameters->Rmax[j] = Rmax[j];
1083  original_parameters->Rmin[j] = Rmin[j];
1084  }
1087  original_parameters->Num_z_planes = countPlanes;
1088 
1089  }
1090  else // Set parameters(r,z) with Rmin==0 as convention
1091  {
1092 #ifdef G4SPECSDEBUG
1093  std::ostringstream message;
1094  message << "Polycone " << GetName() << G4endl
1095  << "cannot be converted to Polycone with (Rmin,Rmaz,Z) parameters!";
1096  G4Exception("G4Polycone::SetOriginalParameters()", "GeomSolids0002",
1097  JustWarning, message);
1098 #endif
1100  original_parameters->Z_values = new G4double[numPlanes];
1101  original_parameters->Rmin = new G4double[numPlanes];
1102  original_parameters->Rmax = new G4double[numPlanes];
1103 
1104  for(G4int j=0; j < numPlanes; j++)
1105  {
1107  original_parameters->Rmax[j] = corners[j].r;
1108  original_parameters->Rmin[j] = 0.0;
1109  }
1112  original_parameters->Num_z_planes = numPlanes;
1113  }
1114  return isConvertible;
1115 }
1116 
1117 #endif
G4String GetName() const
G4bool CrossesItself(G4double tolerance)
G4Polyhedron * CreatePolyhedron() const
Definition: G4Polycone.cc:900
ThreeVector shoot(const G4int Ap, const G4int Af)
G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Polycone.cc:520
G4bool fRebuildPolyhedron
G4Polyhedron * fpPolyhedron
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector GetPointOnSurface() const
Definition: G4Polycone.cc:812
G4double Amin() const
G4double z
Definition: TRTMaterials.hh:39
G4double endPhi
Definition: G4Polycone.hh:184
G4ThreeVector GetPointOnRing(G4double fRMin, G4double fRMax, G4double fRMin2, G4double fRMax2, G4double zOne) const
Definition: G4Polycone.cc:750
G4String name
Definition: TRTMaterials.hh:40
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Polycone.cc:492
G4double a
Definition: TRTMaterials.hh:39
G4bool MustBeOutside(const G4ThreeVector &p) const
G4PolyconeSideRZ * corners
Definition: G4Polycone.hh:187
int G4int
Definition: G4Types.hh:78
void DumpInfo() const
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Polycone.cc:546
G4ThreeVector GetPointOnCone(G4double fRmin1, G4double fRmax1, G4double fRmin2, G4double fRmax2, G4double zOne, G4double zTwo, G4double &totArea) const
Definition: G4Polycone.cc:598
G4GeometryType GetEntityType() const
Definition: G4Polycone.cc:530
G4ThreeVector GetPointOnCut(G4double fRMin1, G4double fRMax1, G4double fRMin2, G4double fRMax2, G4double zOne, G4double zTwo, G4double &totArea) const
Definition: G4Polycone.cc:793
G4bool RemoveDuplicateVertices(G4double tolerance)
G4bool RemoveRedundantVertices(G4double tolerance)
G4Polycone & operator=(const G4Polycone &source)
Definition: G4Polycone.cc:377
G4GLOB_DLL std::ostream G4cout
G4int NumVertices() const
G4double Bmax() const
G4double startPhi
Definition: G4Polycone.hh:183
bool G4bool
Definition: G4Types.hh:79
static const double twopi
Definition: G4SIunits.hh:75
virtual ~G4Polycone()
Definition: G4Polycone.cc:356
void SetOriginalParameters(G4PolyconeHistorical *pars)
EInside Inside(const G4ThreeVector &p) const
Definition: G4Polycone.cc:472
const G4int n
G4bool Reset()
Definition: G4Polycone.cc:441
G4VCSGface ** faces
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4VSolid * Clone() const
Definition: G4Polycone.cc:538
void CopyStuff(const G4Polycone &source)
Definition: G4Polycone.cc:397
static const double pi
Definition: G4SIunits.hh:74
EInside
Definition: geomdefs.hh:58
G4EnclosingCylinder * enclosingCylinder
Definition: G4Polycone.hh:192
G4ThreeVector GetPointOnTubs(G4double fRMin, G4double fRMax, G4double zOne, G4double zTwo, G4double &totArea) const
Definition: G4Polycone.cc:691
static const double degree
Definition: G4SIunits.hh:143
G4bool BisectedBy(G4double a1, G4double b1, G4double a2, G4double b2, G4double tolerance)
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:304
G4Polycone(const G4String &name, G4double phiStart, G4double phiTotal, G4int numZPlanes, const G4double zPlane[], const G4double rInner[], const G4double rOuter[])
Definition: G4Polycone.cc:58
G4bool phiIsOpen
Definition: G4Polycone.hh:185
G4int numCorner
Definition: G4Polycone.hh:186
void Create(G4double phiStart, G4double phiTotal, G4ReduciblePolygon *rz)
Definition: G4Polycone.cc:176
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4PolyconeHistorical * original_parameters
Definition: G4Polycone.hh:188
virtual EInside Inside(const G4ThreeVector &p) const
G4bool ShouldMiss(const G4ThreeVector &p, const G4ThreeVector &v) const