Geant4  10.00.p03
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 84624 2014-10-17 09:56:00Z 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 ) startPhi += twopi;
245 
246  endPhi = phiStart+phiTotal;
247  while( endPhi < startPhi ) endPhi += twopi;
248  }
249 
250  //
251  // Allocate corner array.
252  //
254 
255  //
256  // Copy corners
257  //
258  G4ReduciblePolygonIterator iterRZ(rz);
259 
260  G4PolyconeSideRZ *next = corners;
261  iterRZ.Begin();
262  do
263  {
264  next->r = iterRZ.GetA();
265  next->z = iterRZ.GetB();
266  } while( ++next, iterRZ.Next() );
267 
268  //
269  // Allocate face pointer array
270  //
272  faces = new G4VCSGface*[numFace];
273 
274  //
275  // Construct conical faces
276  //
277  // But! Don't construct a face if both points are at zero radius!
278  //
279  G4PolyconeSideRZ *corner = corners,
280  *prev = corners + numCorner-1,
281  *nextNext;
282  G4VCSGface **face = faces;
283  do
284  {
285  next = corner+1;
286  if (next >= corners+numCorner) next = corners;
287  nextNext = next+1;
288  if (nextNext >= corners+numCorner) nextNext = corners;
289 
290  if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
291 
292  //
293  // We must decide here if we can dare declare one of our faces
294  // as having a "valid" normal (i.e. allBehind = true). This
295  // is never possible if the face faces "inward" in r.
296  //
297  G4bool allBehind;
298  if (corner->z > next->z)
299  {
300  allBehind = false;
301  }
302  else
303  {
304  //
305  // Otherwise, it is only true if the line passing
306  // through the two points of the segment do not
307  // split the r/z cross section
308  //
309  allBehind = !rz->BisectedBy( corner->r, corner->z,
310  next->r, next->z, kCarTolerance );
311  }
312 
313  *face++ = new G4PolyconeSide( prev, corner, next, nextNext,
314  startPhi, endPhi-startPhi, phiIsOpen, allBehind );
315  } while( prev=corner, corner=next, corner > corners );
316 
317  if (phiIsOpen)
318  {
319  //
320  // Construct phi open edges
321  //
322  *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi );
323  *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi );
324  }
325 
326  //
327  // We might have dropped a face or two: recalculate numFace
328  //
329  numFace = face-faces;
330 
331  //
332  // Make enclosingCylinder
333  //
335  new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
336 }
337 
338 
339 //
340 // Fake default constructor - sets only member data and allocates memory
341 // for usage restricted to object persistency.
342 //
344  : G4VCSGfaceted(a), startPhi(0.), endPhi(0.), phiIsOpen(false),
345  numCorner(0), corners(0),original_parameters(0),
346  enclosingCylinder(0)
347 {
348 }
349 
350 
351 //
352 // Destructor
353 //
355 {
356  delete [] corners;
357  delete original_parameters;
358  delete enclosingCylinder;
359 }
360 
361 
362 //
363 // Copy constructor
364 //
366  : G4VCSGfaceted( source )
367 {
368  CopyStuff( source );
369 }
370 
371 
372 //
373 // Assignment operator
374 //
376 {
377  if (this == &source) return *this;
378 
379  G4VCSGfaceted::operator=( source );
380 
381  delete [] corners;
383 
384  delete enclosingCylinder;
385 
386  CopyStuff( source );
387 
388  return *this;
389 }
390 
391 
392 //
393 // CopyStuff
394 //
395 void G4Polycone::CopyStuff( const G4Polycone &source )
396 {
397  //
398  // Simple stuff
399  //
400  startPhi = source.startPhi;
401  endPhi = source.endPhi;
402  phiIsOpen = source.phiIsOpen;
403  numCorner = source.numCorner;
404 
405  //
406  // The corner array
407  //
409 
410  G4PolyconeSideRZ *corn = corners,
411  *sourceCorn = source.corners;
412  do
413  {
414  *corn = *sourceCorn;
415  } while( ++sourceCorn, ++corn < corners+numCorner );
416 
417  //
418  // Original parameters
419  //
420  if (source.original_parameters)
421  {
424  }
425 
426  //
427  // Enclosing cylinder
428  //
430 
431  fRebuildPolyhedron = false;
432  fpPolyhedron = 0;
433 }
434 
435 
436 //
437 // Reset
438 //
440 {
441  //
442  // Clear old setup
443  //
445  delete [] corners;
446  delete enclosingCylinder;
447 
448  //
449  // Rebuild polycone
450  //
451  G4ReduciblePolygon *rz =
458  delete rz;
459 
460  return 0;
461 }
462 
463 
464 //
465 // Inside
466 //
467 // This is an override of G4VCSGfaceted::Inside, created in order
468 // to speed things up by first checking with G4EnclosingCylinder.
469 //
471 {
472  //
473  // Quick test
474  //
475  if (enclosingCylinder->MustBeOutside(p)) return kOutside;
476 
477  //
478  // Long answer
479  //
480  return G4VCSGfaceted::Inside(p);
481 }
482 
483 
484 //
485 // DistanceToIn
486 //
487 // This is an override of G4VCSGfaceted::Inside, created in order
488 // to speed things up by first checking with G4EnclosingCylinder.
489 //
491  const G4ThreeVector &v ) const
492 {
493  //
494  // Quick test
495  //
496  if (enclosingCylinder->ShouldMiss(p,v))
497  return kInfinity;
498 
499  //
500  // Long answer
501  //
502  return G4VCSGfaceted::DistanceToIn( p, v );
503 }
504 
505 
506 //
507 // DistanceToIn
508 //
510 {
511  return G4VCSGfaceted::DistanceToIn(p);
512 }
513 
514 
515 //
516 // ComputeDimensions
517 //
519  const G4int n,
520  const G4VPhysicalVolume* pRep )
521 {
522  p->ComputeDimensions(*this,n,pRep);
523 }
524 
525 //
526 // GetEntityType
527 //
529 {
530  return G4String("G4Polycone");
531 }
532 
533 //
534 // Make a clone of the object
535 //
537 {
538  return new G4Polycone(*this);
539 }
540 
541 //
542 // Stream object contents to an output stream
543 //
544 std::ostream& G4Polycone::StreamInfo( std::ostream& os ) const
545 {
546  G4int oldprc = os.precision(16);
547  os << "-----------------------------------------------------------\n"
548  << " *** Dump for solid - " << GetName() << " ***\n"
549  << " ===================================================\n"
550  << " Solid type: G4Polycone\n"
551  << " Parameters: \n"
552  << " starting phi angle : " << startPhi/degree << " degrees \n"
553  << " ending phi angle : " << endPhi/degree << " degrees \n";
554  G4int i=0;
555 
557  os << " number of Z planes: " << numPlanes << "\n"
558  << " Z values: \n";
559  for (i=0; i<numPlanes; i++)
560  {
561  os << " Z plane " << i << ": "
562  << original_parameters->Z_values[i] << "\n";
563  }
564  os << " Tangent distances to inner surface (Rmin): \n";
565  for (i=0; i<numPlanes; i++)
566  {
567  os << " Z plane " << i << ": "
568  << original_parameters->Rmin[i] << "\n";
569  }
570  os << " Tangent distances to outer surface (Rmax): \n";
571  for (i=0; i<numPlanes; i++)
572  {
573  os << " Z plane " << i << ": "
574  << original_parameters->Rmax[i] << "\n";
575  }
576 
577  os << " number of RZ points: " << numCorner << "\n"
578  << " RZ values (corners): \n";
579  for (i=0; i<numCorner; i++)
580  {
581  os << " "
582  << corners[i].r << ", " << corners[i].z << "\n";
583  }
584  os << "-----------------------------------------------------------\n";
585  os.precision(oldprc);
586 
587  return os;
588 }
589 
590 
591 //
592 // GetPointOnCone
593 //
594 // Auxiliary method for Get Point On Surface
595 //
597  G4double fRmin2, G4double fRmax2,
598  G4double zOne, G4double zTwo,
599  G4double& totArea) const
600 {
601  // declare working variables
602  //
603  G4double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
604  G4double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
605  G4double fDz=(zTwo-zOne)/2., afDz=std::fabs(fDz);
606  G4ThreeVector point, offset=G4ThreeVector(0.,0.,0.5*(zTwo+zOne));
607  fDPhi = endPhi - startPhi;
608  rone = (fRmax1-fRmax2)/(2.*fDz);
609  rtwo = (fRmin1-fRmin2)/(2.*fDz);
610  if(fRmax1==fRmax2){qone=0.;}
611  else{
612  qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
613  }
614  if(fRmin1==fRmin2){qtwo=0.;}
615  else{
616  qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
617  }
618  Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*(sqr(fRmin1-fRmin2)+sqr(zTwo-zOne));
619  Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*(sqr(fRmax1-fRmax2)+sqr(zTwo-zOne));
620  Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
621  totArea = Aone+Atwo+2.*Afive;
622 
623  phi = RandFlat::shoot(startPhi,endPhi);
624  cosu = std::cos(phi);
625  sinu = std::sin(phi);
626 
627 
628  if( (startPhi == 0) && (endPhi == twopi) ) { Afive = 0; }
629  chose = RandFlat::shoot(0.,Aone+Atwo+2.*Afive);
630  if( (chose >= 0) && (chose < Aone) )
631  {
632  if(fRmax1 != fRmax2)
633  {
634  zRand = RandFlat::shoot(-1.*afDz,afDz);
635  point = G4ThreeVector (rone*cosu*(qone-zRand),
636  rone*sinu*(qone-zRand), zRand);
637  }
638  else
639  {
640  point = G4ThreeVector(fRmax1*cosu, fRmax1*sinu,
641  RandFlat::shoot(-1.*afDz,afDz));
642 
643  }
644  }
645  else if(chose >= Aone && chose < Aone + Atwo)
646  {
647  if(fRmin1 != fRmin2)
648  {
649  zRand = RandFlat::shoot(-1.*afDz,afDz);
650  point = G4ThreeVector (rtwo*cosu*(qtwo-zRand),
651  rtwo*sinu*(qtwo-zRand), zRand);
652 
653  }
654  else
655  {
656  point = G4ThreeVector(fRmin1*cosu, fRmin1*sinu,
657  RandFlat::shoot(-1.*afDz,afDz));
658  }
659  }
660  else if( (chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive) )
661  {
662  zRand = RandFlat::shoot(-1.*afDz,afDz);
663  rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
664  rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
665  rRand1 = std::sqrt(RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
666  point = G4ThreeVector (rRand1*std::cos(startPhi),
667  rRand1*std::sin(startPhi), zRand);
668  }
669  else
670  {
671  zRand = RandFlat::shoot(-1.*afDz,afDz);
672  rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
673  rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
674  rRand1 = std::sqrt(RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
675  point = G4ThreeVector (rRand1*std::cos(endPhi),
676  rRand1*std::sin(endPhi), zRand);
677 
678  }
679 
680  return point+offset;
681 }
682 
683 
684 //
685 // GetPointOnTubs
686 //
687 // Auxiliary method for GetPoint On Surface
688 //
690  G4double zOne, G4double zTwo,
691  G4double& totArea) const
692 {
693  G4double xRand,yRand,zRand,phi,cosphi,sinphi,chose,
694  aOne,aTwo,aFou,rRand,fDz,fSPhi,fDPhi;
695  fDz = std::fabs(0.5*(zTwo-zOne));
696  fSPhi = startPhi;
697  fDPhi = endPhi-startPhi;
698 
699  aOne = 2.*fDz*fDPhi*fRMax;
700  aTwo = 2.*fDz*fDPhi*fRMin;
701  aFou = 2.*fDz*(fRMax-fRMin);
702  totArea = aOne+aTwo+2.*aFou;
703  phi = RandFlat::shoot(startPhi,endPhi);
704  cosphi = std::cos(phi);
705  sinphi = std::sin(phi);
706  rRand = fRMin + (fRMax-fRMin)*std::sqrt(RandFlat::shoot());
707 
708  if(startPhi == 0 && endPhi == twopi)
709  aFou = 0;
710 
711  chose = RandFlat::shoot(0.,aOne+aTwo+2.*aFou);
712  if( (chose >= 0) && (chose < aOne) )
713  {
714  xRand = fRMax*cosphi;
715  yRand = fRMax*sinphi;
716  zRand = RandFlat::shoot(-1.*fDz,fDz);
717  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
718  }
719  else if( (chose >= aOne) && (chose < aOne + aTwo) )
720  {
721  xRand = fRMin*cosphi;
722  yRand = fRMin*sinphi;
723  zRand = RandFlat::shoot(-1.*fDz,fDz);
724  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
725  }
726  else if( (chose >= aOne+aTwo) && (chose <aOne+aTwo+aFou) )
727  {
728  xRand = rRand*std::cos(fSPhi+fDPhi);
729  yRand = rRand*std::sin(fSPhi+fDPhi);
730  zRand = RandFlat::shoot(-1.*fDz,fDz);
731  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
732  }
733 
734  // else
735 
736  xRand = rRand*std::cos(fSPhi+fDPhi);
737  yRand = rRand*std::sin(fSPhi+fDPhi);
738  zRand = RandFlat::shoot(-1.*fDz,fDz);
739  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
740 }
741 
742 
743 //
744 // GetPointOnRing
745 //
746 // Auxiliary method for GetPoint On Surface
747 //
749  G4double fRMin2,G4double fRMax2,
750  G4double zOne) const
751 {
752  G4double xRand,yRand,phi,cosphi,sinphi,rRand1,rRand2,A1,Atot,rCh;
754  cosphi = std::cos(phi);
755  sinphi = std::sin(phi);
756 
757  if(fRMin1==fRMin2)
758  {
759  rRand1 = fRMin1; A1=0.;
760  }
761  else
762  {
763  rRand1 = RandFlat::shoot(fRMin1,fRMin2);
764  A1=std::fabs(fRMin2*fRMin2-fRMin1*fRMin1);
765  }
766  if(fRMax1==fRMax2)
767  {
768  rRand2=fRMax1; Atot=A1;
769  }
770  else
771  {
772  rRand2 = RandFlat::shoot(fRMax1,fRMax2);
773  Atot = A1+std::fabs(fRMax2*fRMax2-fRMax1*fRMax1);
774  }
775  rCh = RandFlat::shoot(0.,Atot);
776 
777  if(rCh>A1) { rRand1=rRand2; }
778 
779  xRand = rRand1*cosphi;
780  yRand = rRand1*sinphi;
781 
782  return G4ThreeVector(xRand, yRand, zOne);
783 }
784 
785 
786 //
787 // GetPointOnCut
788 //
789 // Auxiliary method for Get Point On Surface
790 //
792  G4double fRMin2, G4double fRMax2,
793  G4double zOne, G4double zTwo,
794  G4double& totArea) const
795 { if(zOne==zTwo)
796  {
797  return GetPointOnRing(fRMin1, fRMax1,fRMin2,fRMax2,zOne);
798  }
799  if( (fRMin1 == fRMin2) && (fRMax1 == fRMax2) )
800  {
801  return GetPointOnTubs(fRMin1, fRMax1,zOne,zTwo,totArea);
802  }
803  return GetPointOnCone(fRMin1,fRMax1,fRMin2,fRMax2,zOne,zTwo,totArea);
804 }
805 
806 
807 //
808 // GetPointOnSurface
809 //
811 {
812  G4double Area=0,totArea=0,Achose1=0,Achose2=0,phi,cosphi,sinphi,rRand;
813  G4int i=0;
815 
817  cosphi = std::cos(phi);
818  sinphi = std::sin(phi);
819 
820  rRand = original_parameters->Rmin[0] +
822  * std::sqrt(RandFlat::shoot()) );
823 
824  std::vector<G4double> areas; // (numPlanes+1);
825  std::vector<G4ThreeVector> points; // (numPlanes-1);
826 
827  areas.push_back(pi*(sqr(original_parameters->Rmax[0])
828  -sqr(original_parameters->Rmin[0])));
829 
830  for(i=0; i<numPlanes-1; i++)
831  {
833  * std::sqrt(sqr(original_parameters->Rmin[i]
834  -original_parameters->Rmin[i+1])+
837 
839  * std::sqrt(sqr(original_parameters->Rmax[i]
840  -original_parameters->Rmax[i+1])+
843 
844  Area *= 0.5*(endPhi-startPhi);
845 
846  if(startPhi==0.&& endPhi == twopi)
847  {
848  Area += std::fabs(original_parameters->Z_values[i+1]
853  -original_parameters->Rmin[i+1]);
854  }
855  areas.push_back(Area);
856  totArea += Area;
857  }
858 
859  areas.push_back(pi*(sqr(original_parameters->Rmax[numPlanes-1])-
860  sqr(original_parameters->Rmin[numPlanes-1])));
861 
862  totArea += (areas[0]+areas[numPlanes]);
863  G4double chose = RandFlat::shoot(0.,totArea);
864 
865  if( (chose>=0.) && (chose<areas[0]) )
866  {
867  return G4ThreeVector(rRand*cosphi, rRand*sinphi,
869  }
870 
871  for (i=0; i<numPlanes-1; i++)
872  {
873  Achose1 += areas[i];
874  Achose2 = (Achose1+areas[i+1]);
875  if(chose>=Achose1 && chose<Achose2)
876  {
882  original_parameters->Z_values[i+1], Area);
883  }
884  }
885 
886  rRand = original_parameters->Rmin[numPlanes-1] +
887  ( (original_parameters->Rmax[numPlanes-1]-original_parameters->Rmin[numPlanes-1])
888  * std::sqrt(RandFlat::shoot()) );
889 
890  return G4ThreeVector(rRand*cosphi,rRand*sinphi,
891  original_parameters->Z_values[numPlanes-1]);
892 
893 }
894 
895 //
896 // CreatePolyhedron
897 //
899 {
900  //
901  // This has to be fixed in visualization. Fake it for the moment.
902  //
903 
910 }
911 
913 {
914  G4int numPlanes = (G4int)numCorner;
915  G4bool isConvertible=true;
916  G4double Zmax=rz->Bmax();
917  rz->StartWithZMin();
918 
919  // Prepare vectors for storage
920  //
921  std::vector<G4double> Z;
922  std::vector<G4double> Rmin;
923  std::vector<G4double> Rmax;
924 
925  G4int countPlanes=1;
926  G4int icurr=0;
927  G4int icurl=0;
928 
929  // first plane Z=Z[0]
930  //
931  Z.push_back(corners[0].z);
932  G4double Zprev=Z[0];
933  if (Zprev == corners[1].z)
934  {
935  Rmin.push_back(corners[0].r);
936  Rmax.push_back (corners[1].r);icurr=1;
937  }
938  else if (Zprev == corners[numPlanes-1].z)
939  {
940  Rmin.push_back(corners[numPlanes-1].r);
941  Rmax.push_back (corners[0].r);
942  icurl=numPlanes-1;
943  }
944  else
945  {
946  Rmin.push_back(corners[0].r);
947  Rmax.push_back (corners[0].r);
948  }
949 
950  // next planes until last
951  //
952  G4int inextr=0, inextl=0;
953  for (G4int i=0; i < numPlanes-2; i++)
954  {
955  inextr=1+icurr;
956  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
957 
958  if((corners[inextr].z >= Zmax) & (corners[inextl].z >= Zmax)) { break; }
959 
960  G4double Zleft = corners[inextl].z;
961  G4double Zright = corners[inextr].z;
962  if(Zright>Zleft)
963  {
964  Z.push_back(Zleft);
965  countPlanes++;
966  G4double difZr=corners[inextr].z - corners[icurr].z;
967  G4double difZl=corners[inextl].z - corners[icurl].z;
968 
969  if(std::fabs(difZl) < kCarTolerance)
970  {
971  if(corners[inextl].r >= corners[icurl].r)
972  {
973  Rmin.push_back(corners[icurl].r);
974  Rmax.push_back(Rmax[countPlanes-2]);
975  Rmax[countPlanes-2]=corners[icurl].r;
976  }
977  else
978  {
979  Rmin.push_back(corners[inextl].r);
980  Rmax.push_back(corners[icurl].r);
981  }
982  }
983  else if (difZl >= kCarTolerance)
984  {
985  Rmin.push_back(corners[inextl].r);
986  Rmax.push_back (corners[icurr].r + (Zleft-corners[icurr].z)/difZr
987  *(corners[inextr].r - corners[icurr].r));
988  }
989  else
990  {
991  isConvertible=false; break;
992  }
993  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
994  }
995  else if(std::fabs(Zright-Zleft)<kCarTolerance) // Zright=Zleft
996  {
997  Z.push_back(Zleft);
998  countPlanes++;
999  icurr++;
1000 
1001  icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1002 
1003  Rmin.push_back(corners[inextl].r);
1004  Rmax.push_back (corners[inextr].r);
1005  }
1006  else // Zright<Zleft
1007  {
1008  Z.push_back(Zright);
1009  countPlanes++;
1010 
1011  G4double difZr=corners[inextr].z - corners[icurr].z;
1012  G4double difZl=corners[inextl].z - corners[icurl].z;
1013  if(std::fabs(difZr) < kCarTolerance)
1014  {
1015  if(corners[inextr].r >= corners[icurr].r)
1016  {
1017  Rmin.push_back(corners[icurr].r);
1018  Rmax.push_back(corners[inextr].r);
1019  }
1020  else
1021  {
1022  Rmin.push_back(corners[inextr].r);
1023  Rmax.push_back(corners[icurr].r);
1024  Rmax[countPlanes-2]=corners[inextr].r;
1025  }
1026  icurr++;
1027  } // plate
1028  else if (difZr >= kCarTolerance)
1029  {
1030  if(std::fabs(difZl)<kCarTolerance)
1031  {
1032  Rmax.push_back(corners[inextr].r);
1033  Rmin.push_back (corners[icurr].r);
1034  }
1035  else
1036  {
1037  Rmax.push_back(corners[inextr].r);
1038  Rmin.push_back (corners[icurl].r+(Zright-corners[icurl].z)/difZl
1039  * (corners[inextl].r - corners[icurl].r));
1040  }
1041  icurr++;
1042  }
1043  else
1044  {
1045  isConvertible=false; break;
1046  }
1047  }
1048  } // end for loop
1049 
1050  // last plane Z=Zmax
1051  //
1052  Z.push_back(Zmax);
1053  countPlanes++;
1054  inextr=1+icurr;
1055  inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1056 
1057  if(corners[inextr].z==corners[inextl].z)
1058  {
1059  Rmax.push_back(corners[inextr].r);
1060  Rmin.push_back(corners[inextl].r);
1061  }
1062  else
1063  {
1064  Rmax.push_back(corners[inextr].r);
1065  Rmin.push_back(corners[inextl].r);
1066  }
1067 
1068  // Set original parameters Rmin,Rmax,Z
1069  //
1070  if(isConvertible)
1071  {
1073  original_parameters->Z_values = new G4double[countPlanes];
1074  original_parameters->Rmin = new G4double[countPlanes];
1075  original_parameters->Rmax = new G4double[countPlanes];
1076 
1077  for(G4int j=0; j < countPlanes; j++)
1078  {
1079  original_parameters->Z_values[j] = Z[j];
1080  original_parameters->Rmax[j] = Rmax[j];
1081  original_parameters->Rmin[j] = Rmin[j];
1082  }
1085  original_parameters->Num_z_planes = countPlanes;
1086 
1087  }
1088  else // Set parameters(r,z) with Rmin==0 as convention
1089  {
1090 #ifdef G4SPECSDEBUG
1091  std::ostringstream message;
1092  message << "Polycone " << GetName() << G4endl
1093  << "cannot be converted to Polycone with (Rmin,Rmaz,Z) parameters!";
1094  G4Exception("G4Polycone::SetOriginalParameters()", "GeomSolids0002",
1095  JustWarning, message);
1096 #endif
1098  original_parameters->Z_values = new G4double[numPlanes];
1099  original_parameters->Rmin = new G4double[numPlanes];
1100  original_parameters->Rmax = new G4double[numPlanes];
1101 
1102  for(G4int j=0; j < numPlanes; j++)
1103  {
1105  original_parameters->Rmax[j] = corners[j].r;
1106  original_parameters->Rmin[j] = 0.0;
1107  }
1110  original_parameters->Num_z_planes = numPlanes;
1111  }
1112  return isConvertible;
1113 }
1114 
1115 #endif
G4String GetName() const
G4bool CrossesItself(G4double tolerance)
G4Polyhedron * CreatePolyhedron() const
Definition: G4Polycone.cc:898
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:518
G4bool fRebuildPolyhedron
G4Polyhedron * fpPolyhedron
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector GetPointOnSurface() const
Definition: G4Polycone.cc:810
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:748
G4String name
Definition: TRTMaterials.hh:40
const G4double pi
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Polycone.cc:490
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:544
G4ThreeVector GetPointOnCone(G4double fRmin1, G4double fRmax1, G4double fRmin2, G4double fRmax2, G4double zOne, G4double zTwo, G4double &totArea) const
Definition: G4Polycone.cc:596
G4GeometryType GetEntityType() const
Definition: G4Polycone.cc:528
G4ThreeVector GetPointOnCut(G4double fRMin1, G4double fRMax1, G4double fRMin2, G4double fRMax2, G4double zOne, G4double zTwo, G4double &totArea) const
Definition: G4Polycone.cc:791
G4bool RemoveDuplicateVertices(G4double tolerance)
G4bool RemoveRedundantVertices(G4double tolerance)
G4Polycone & operator=(const G4Polycone &source)
Definition: G4Polycone.cc:375
G4GLOB_DLL std::ostream G4cout
G4int NumVertices() const
G4double Bmax() const
G4double startPhi
Definition: G4Polycone.hh:183
bool G4bool
Definition: G4Types.hh:79
virtual ~G4Polycone()
Definition: G4Polycone.cc:354
void SetOriginalParameters(G4PolyconeHistorical *pars)
EInside Inside(const G4ThreeVector &p) const
Definition: G4Polycone.cc:470
const G4int n
G4bool Reset()
Definition: G4Polycone.cc:439
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:536
void CopyStuff(const G4Polycone &source)
Definition: G4Polycone.cc:395
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:689
static const double degree
Definition: G4SIunits.hh:125
G4bool BisectedBy(G4double a1, G4double b1, G4double a2, G4double b2, G4double tolerance)
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:305
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