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