Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 69790 2013-05-15 12:39:10Z 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 #include "G4PolyconeSide.hh"
43 #include "G4PolyPhiFace.hh"
44 
45 #include "Randomize.hh"
46 
47 #include "G4Polyhedron.hh"
48 #include "G4EnclosingCylinder.hh"
49 #include "G4ReduciblePolygon.hh"
50 #include "G4VPVParameterisation.hh"
51 
52 using namespace CLHEP;
53 
54 //
55 // Constructor (GEANT3 style parameters)
56 //
57 G4Polycone::G4Polycone( const G4String& name,
58  G4double phiStart,
59  G4double phiTotal,
60  G4int numZPlanes,
61  const G4double zPlane[],
62  const G4double rInner[],
63  const G4double rOuter[] )
64  : G4VCSGfaceted( name ), genericPcon(false)
65 {
66  //
67  // Some historical ugliness
68  //
69  original_parameters = new G4PolyconeHistorical();
70 
71  original_parameters->Start_angle = phiStart;
72  original_parameters->Opening_angle = phiTotal;
73  original_parameters->Num_z_planes = numZPlanes;
74  original_parameters->Z_values = new G4double[numZPlanes];
75  original_parameters->Rmin = new G4double[numZPlanes];
76  original_parameters->Rmax = new G4double[numZPlanes];
77 
78  G4int i;
79  for (i=0; i<numZPlanes; i++)
80  {
81  if(rInner[i]>rOuter[i])
82  {
83  DumpInfo();
84  std::ostringstream message;
85  message << "Cannot create a Polycone with rInner > rOuter for the same Z"
86  << G4endl
87  << " rInner > rOuter for the same Z !" << G4endl
88  << " rMin[" << i << "] = " << rInner[i]
89  << " -- rMax[" << i << "] = " << rOuter[i];
90  G4Exception("G4Polycone::G4Polycone()", "GeomSolids0002",
91  FatalErrorInArgument, message);
92  }
93  if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
94  {
95  if( (rInner[i] > rOuter[i+1])
96  ||(rInner[i+1] > rOuter[i]) )
97  {
98  DumpInfo();
99  std::ostringstream message;
100  message << "Cannot create a Polycone with no contiguous segments."
101  << G4endl
102  << " Segments are not contiguous !" << G4endl
103  << " rMin[" << i << "] = " << rInner[i]
104  << " -- rMax[" << i+1 << "] = " << rOuter[i+1] << G4endl
105  << " rMin[" << i+1 << "] = " << rInner[i+1]
106  << " -- rMax[" << i << "] = " << rOuter[i];
107  G4Exception("G4Polycone::G4Polycone()", "GeomSolids0002",
108  FatalErrorInArgument, message);
109  }
110  }
111  original_parameters->Z_values[i] = zPlane[i];
112  original_parameters->Rmin[i] = rInner[i];
113  original_parameters->Rmax[i] = rOuter[i];
114  }
115 
116  //
117  // Build RZ polygon using special PCON/PGON GEANT3 constructor
118  //
119  G4ReduciblePolygon *rz =
120  new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
121 
122  //
123  // Do the real work
124  //
125  Create( phiStart, phiTotal, rz );
126 
127  delete rz;
128 }
129 
130 
131 //
132 // Constructor (generic parameters)
133 //
134 G4Polycone::G4Polycone( const G4String& name,
135  G4double phiStart,
136  G4double phiTotal,
137  G4int numRZ,
138  const G4double r[],
139  const G4double z[] )
140  : G4VCSGfaceted( name ), genericPcon(true)
141 {
142 
143  G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
144 
145  Create( phiStart, phiTotal, rz );
146 
147  // Set original_parameters struct for consistency
148  //
149  SetOriginalParameters();
150 
151  delete rz;
152 }
153 
154 
155 //
156 // Create
157 //
158 // Generic create routine, called by each constructor after
159 // conversion of arguments
160 //
161 void G4Polycone::Create( G4double phiStart,
162  G4double phiTotal,
163  G4ReduciblePolygon *rz )
164 {
165  //
166  // Perform checks of rz values
167  //
168  if (rz->Amin() < 0.0)
169  {
170  std::ostringstream message;
171  message << "Illegal input parameters - " << GetName() << G4endl
172  << " All R values must be >= 0 !";
173  G4Exception("G4Polycone::Create()", "GeomSolids0002",
174  FatalErrorInArgument, message);
175  }
176 
177  G4double rzArea = rz->Area();
178  if (rzArea < -kCarTolerance)
179  rz->ReverseOrder();
180 
181  else if (rzArea < -kCarTolerance)
182  {
183  std::ostringstream message;
184  message << "Illegal input parameters - " << GetName() << G4endl
185  << " R/Z cross section is zero or near zero: " << rzArea;
186  G4Exception("G4Polycone::Create()", "GeomSolids0002",
187  FatalErrorInArgument, message);
188  }
189 
190  if ( (!rz->RemoveDuplicateVertices( kCarTolerance ))
191  || (!rz->RemoveRedundantVertices( kCarTolerance )) )
192  {
193  std::ostringstream message;
194  message << "Illegal input parameters - " << GetName() << G4endl
195  << " Too few unique R/Z values !";
196  G4Exception("G4Polycone::Create()", "GeomSolids0002",
197  FatalErrorInArgument, message);
198  }
199 
200  if (rz->CrossesItself(1/kInfinity))
201  {
202  std::ostringstream message;
203  message << "Illegal input parameters - " << GetName() << G4endl
204  << " R/Z segments cross !";
205  G4Exception("G4Polycone::Create()", "GeomSolids0002",
206  FatalErrorInArgument, message);
207  }
208 
209  numCorner = rz->NumVertices();
210 
211  //
212  // Phi opening? Account for some possible roundoff, and interpret
213  // nonsense value as representing no phi opening
214  //
215  if (phiTotal <= 0 || phiTotal > twopi-1E-10)
216  {
217  phiIsOpen = false;
218  startPhi = 0;
219  endPhi = twopi;
220  }
221  else
222  {
223  phiIsOpen = true;
224 
225  //
226  // Convert phi into our convention
227  //
228  startPhi = phiStart;
229  while( startPhi < 0 ) startPhi += twopi;
230 
231  endPhi = phiStart+phiTotal;
232  while( endPhi < startPhi ) endPhi += twopi;
233  }
234 
235  //
236  // Allocate corner array.
237  //
238  corners = new G4PolyconeSideRZ[numCorner];
239 
240  //
241  // Copy corners
242  //
243  G4ReduciblePolygonIterator iterRZ(rz);
244 
245  G4PolyconeSideRZ *next = corners;
246  iterRZ.Begin();
247  do
248  {
249  next->r = iterRZ.GetA();
250  next->z = iterRZ.GetB();
251  } while( ++next, iterRZ.Next() );
252 
253  //
254  // Allocate face pointer array
255  //
256  numFace = phiIsOpen ? numCorner+2 : numCorner;
257  faces = new G4VCSGface*[numFace];
258 
259  //
260  // Construct conical faces
261  //
262  // But! Don't construct a face if both points are at zero radius!
263  //
264  G4PolyconeSideRZ *corner = corners,
265  *prev = corners + numCorner-1,
266  *nextNext;
267  G4VCSGface **face = faces;
268  do
269  {
270  next = corner+1;
271  if (next >= corners+numCorner) next = corners;
272  nextNext = next+1;
273  if (nextNext >= corners+numCorner) nextNext = corners;
274 
275  if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
276 
277  //
278  // We must decide here if we can dare declare one of our faces
279  // as having a "valid" normal (i.e. allBehind = true). This
280  // is never possible if the face faces "inward" in r.
281  //
282  G4bool allBehind;
283  if (corner->z > next->z)
284  {
285  allBehind = false;
286  }
287  else
288  {
289  //
290  // Otherwise, it is only true if the line passing
291  // through the two points of the segment do not
292  // split the r/z cross section
293  //
294  allBehind = !rz->BisectedBy( corner->r, corner->z,
295  next->r, next->z, kCarTolerance );
296  }
297 
298  *face++ = new G4PolyconeSide( prev, corner, next, nextNext,
299  startPhi, endPhi-startPhi, phiIsOpen, allBehind );
300  } while( prev=corner, corner=next, corner > corners );
301 
302  if (phiIsOpen)
303  {
304  //
305  // Construct phi open edges
306  //
307  *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi );
308  *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi );
309  }
310 
311  //
312  // We might have dropped a face or two: recalculate numFace
313  //
314  numFace = face-faces;
315 
316  //
317  // Make enclosingCylinder
318  //
319  enclosingCylinder =
320  new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
321 }
322 
323 
324 //
325 // Fake default constructor - sets only member data and allocates memory
326 // for usage restricted to object persistency.
327 //
328 G4Polycone::G4Polycone( __void__& a )
329  : G4VCSGfaceted(a), startPhi(0.), endPhi(0.), phiIsOpen(false),
330  genericPcon(false), numCorner(0), corners(0),
331  original_parameters(0), enclosingCylinder(0)
332 {
333 }
334 
335 
336 //
337 // Destructor
338 //
339 G4Polycone::~G4Polycone()
340 {
341  delete [] corners;
342  delete original_parameters;
343  delete enclosingCylinder;
344 }
345 
346 
347 //
348 // Copy constructor
349 //
350 G4Polycone::G4Polycone( const G4Polycone &source )
351  : G4VCSGfaceted( source )
352 {
353  CopyStuff( source );
354 }
355 
356 
357 //
358 // Assignment operator
359 //
360 const G4Polycone &G4Polycone::operator=( const G4Polycone &source )
361 {
362  if (this == &source) return *this;
363 
364  G4VCSGfaceted::operator=( source );
365 
366  delete [] corners;
367  if (original_parameters) delete original_parameters;
368 
369  delete enclosingCylinder;
370 
371  CopyStuff( source );
372 
373  return *this;
374 }
375 
376 
377 //
378 // CopyStuff
379 //
380 void G4Polycone::CopyStuff( const G4Polycone &source )
381 {
382  //
383  // Simple stuff
384  //
385  startPhi = source.startPhi;
386  endPhi = source.endPhi;
387  phiIsOpen = source.phiIsOpen;
388  numCorner = source.numCorner;
389  genericPcon= source.genericPcon;
390 
391  //
392  // The corner array
393  //
394  corners = new G4PolyconeSideRZ[numCorner];
395 
396  G4PolyconeSideRZ *corn = corners,
397  *sourceCorn = source.corners;
398  do
399  {
400  *corn = *sourceCorn;
401  } while( ++sourceCorn, ++corn < corners+numCorner );
402 
403  //
404  // Original parameters
405  //
406  if (source.original_parameters)
407  {
408  original_parameters =
409  new G4PolyconeHistorical( *source.original_parameters );
410  }
411 
412  //
413  // Enclosing cylinder
414  //
415  enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder );
416 }
417 
418 
419 //
420 // Reset
421 //
422 G4bool G4Polycone::Reset()
423 {
424  if (genericPcon)
425  {
426  std::ostringstream message;
427  message << "Solid " << GetName() << " built using generic construct."
428  << G4endl << "Not applicable to the generic construct !";
429  G4Exception("G4Polycone::Reset()", "GeomSolids1001",
430  JustWarning, message, "Parameters NOT resetted.");
431  return 1;
432  }
433 
434  //
435  // Clear old setup
436  //
437  G4VCSGfaceted::DeleteStuff();
438  delete [] corners;
439  delete enclosingCylinder;
440 
441  //
442  // Rebuild polycone
443  //
444  G4ReduciblePolygon *rz =
445  new G4ReduciblePolygon( original_parameters->Rmin,
446  original_parameters->Rmax,
447  original_parameters->Z_values,
448  original_parameters->Num_z_planes );
449  Create( original_parameters->Start_angle,
450  original_parameters->Opening_angle, rz );
451  delete rz;
452 
453  return 0;
454 }
455 
456 
457 //
458 // Inside
459 //
460 // This is an override of G4VCSGfaceted::Inside, created in order
461 // to speed things up by first checking with G4EnclosingCylinder.
462 //
463 EInside G4Polycone::Inside( const G4ThreeVector &p ) const
464 {
465  //
466  // Quick test
467  //
468  if (enclosingCylinder->MustBeOutside(p)) return kOutside;
469 
470  //
471  // Long answer
472  //
473  return G4VCSGfaceted::Inside(p);
474 }
475 
476 
477 //
478 // DistanceToIn
479 //
480 // This is an override of G4VCSGfaceted::Inside, created in order
481 // to speed things up by first checking with G4EnclosingCylinder.
482 //
483 G4double G4Polycone::DistanceToIn( const G4ThreeVector &p,
484  const G4ThreeVector &v ) const
485 {
486  //
487  // Quick test
488  //
489  if (enclosingCylinder->ShouldMiss(p,v))
490  return kInfinity;
491 
492  //
493  // Long answer
494  //
495  return G4VCSGfaceted::DistanceToIn( p, v );
496 }
497 
498 
499 //
500 // DistanceToIn
501 //
502 G4double G4Polycone::DistanceToIn( const G4ThreeVector &p ) const
503 {
504  return G4VCSGfaceted::DistanceToIn(p);
505 }
506 
507 
508 //
509 // ComputeDimensions
510 //
511 void G4Polycone::ComputeDimensions( G4VPVParameterisation* p,
512  const G4int n,
513  const G4VPhysicalVolume* pRep )
514 {
515  p->ComputeDimensions(*this,n,pRep);
516 }
517 
518 //
519 // GetEntityType
520 //
521 G4GeometryType G4Polycone::GetEntityType() const
522 {
523  return G4String("G4Polycone");
524 }
525 
526 //
527 // Make a clone of the object
528 //
529 G4VSolid* G4Polycone::Clone() const
530 {
531  return new G4Polycone(*this);
532 }
533 
534 //
535 // Stream object contents to an output stream
536 //
537 std::ostream& G4Polycone::StreamInfo( std::ostream& os ) const
538 {
539  G4int oldprc = os.precision(16);
540  os << "-----------------------------------------------------------\n"
541  << " *** Dump for solid - " << GetName() << " ***\n"
542  << " ===================================================\n"
543  << " Solid type: G4Polycone\n"
544  << " Parameters: \n"
545  << " starting phi angle : " << startPhi/degree << " degrees \n"
546  << " ending phi angle : " << endPhi/degree << " degrees \n";
547  G4int i=0;
548  if (!genericPcon)
549  {
550  G4int numPlanes = original_parameters->Num_z_planes;
551  os << " number of Z planes: " << numPlanes << "\n"
552  << " Z values: \n";
553  for (i=0; i<numPlanes; i++)
554  {
555  os << " Z plane " << i << ": "
556  << original_parameters->Z_values[i] << "\n";
557  }
558  os << " Tangent distances to inner surface (Rmin): \n";
559  for (i=0; i<numPlanes; i++)
560  {
561  os << " Z plane " << i << ": "
562  << original_parameters->Rmin[i] << "\n";
563  }
564  os << " Tangent distances to outer surface (Rmax): \n";
565  for (i=0; i<numPlanes; i++)
566  {
567  os << " Z plane " << i << ": "
568  << original_parameters->Rmax[i] << "\n";
569  }
570  }
571  os << " number of RZ points: " << numCorner << "\n"
572  << " RZ values (corners): \n";
573  for (i=0; i<numCorner; i++)
574  {
575  os << " "
576  << corners[i].r << ", " << corners[i].z << "\n";
577  }
578  os << "-----------------------------------------------------------\n";
579  os.precision(oldprc);
580 
581  return os;
582 }
583 
584 
585 //
586 // GetPointOnCone
587 //
588 // Auxiliary method for Get Point On Surface
589 //
590 G4ThreeVector G4Polycone::GetPointOnCone(G4double fRmin1, G4double fRmax1,
591  G4double fRmin2, G4double fRmax2,
592  G4double zOne, G4double zTwo,
593  G4double& totArea) const
594 {
595  // declare working variables
596  //
597  G4double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
598  G4double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
599  G4double fDz=(zTwo-zOne)/2., afDz=std::fabs(fDz);
600  G4ThreeVector point, offset=G4ThreeVector(0.,0.,0.5*(zTwo+zOne));
601  fDPhi = endPhi - startPhi;
602  rone = (fRmax1-fRmax2)/(2.*fDz);
603  rtwo = (fRmin1-fRmin2)/(2.*fDz);
604  if(fRmax1==fRmax2){qone=0.;}
605  else{
606  qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
607  }
608  if(fRmin1==fRmin2){qtwo=0.;}
609  else{
610  qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
611  }
612  Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*(sqr(fRmin1-fRmin2)+sqr(zTwo-zOne));
613  Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*(sqr(fRmax1-fRmax2)+sqr(zTwo-zOne));
614  Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
615  totArea = Aone+Atwo+2.*Afive;
616 
617  phi = RandFlat::shoot(startPhi,endPhi);
618  cosu = std::cos(phi);
619  sinu = std::sin(phi);
620 
621 
622  if( (startPhi == 0) && (endPhi == twopi) ) { Afive = 0; }
623  chose = RandFlat::shoot(0.,Aone+Atwo+2.*Afive);
624  if( (chose >= 0) && (chose < Aone) )
625  {
626  if(fRmax1 != fRmax2)
627  {
628  zRand = RandFlat::shoot(-1.*afDz,afDz);
629  point = G4ThreeVector (rone*cosu*(qone-zRand),
630  rone*sinu*(qone-zRand), zRand);
631  }
632  else
633  {
634  point = G4ThreeVector(fRmax1*cosu, fRmax1*sinu,
635  RandFlat::shoot(-1.*afDz,afDz));
636 
637  }
638  }
639  else if(chose >= Aone && chose < Aone + Atwo)
640  {
641  if(fRmin1 != fRmin2)
642  {
643  zRand = RandFlat::shoot(-1.*afDz,afDz);
644  point = G4ThreeVector (rtwo*cosu*(qtwo-zRand),
645  rtwo*sinu*(qtwo-zRand), zRand);
646 
647  }
648  else
649  {
650  point = G4ThreeVector(fRmin1*cosu, fRmin1*sinu,
651  RandFlat::shoot(-1.*afDz,afDz));
652  }
653  }
654  else if( (chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive) )
655  {
656  zRand = RandFlat::shoot(-1.*afDz,afDz);
657  rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
658  rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
659  rRand1 = std::sqrt(RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
660  point = G4ThreeVector (rRand1*std::cos(startPhi),
661  rRand1*std::sin(startPhi), zRand);
662  }
663  else
664  {
665  zRand = RandFlat::shoot(-1.*afDz,afDz);
666  rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
667  rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
668  rRand1 = std::sqrt(RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
669  point = G4ThreeVector (rRand1*std::cos(endPhi),
670  rRand1*std::sin(endPhi), zRand);
671 
672  }
673 
674  return point+offset;
675 }
676 
677 
678 //
679 // GetPointOnTubs
680 //
681 // Auxiliary method for GetPoint On Surface
682 //
683 G4ThreeVector G4Polycone::GetPointOnTubs(G4double fRMin, G4double fRMax,
684  G4double zOne, G4double zTwo,
685  G4double& totArea) const
686 {
687  G4double xRand,yRand,zRand,phi,cosphi,sinphi,chose,
688  aOne,aTwo,aFou,rRand,fDz,fSPhi,fDPhi;
689  fDz = std::fabs(0.5*(zTwo-zOne));
690  fSPhi = startPhi;
691  fDPhi = endPhi-startPhi;
692 
693  aOne = 2.*fDz*fDPhi*fRMax;
694  aTwo = 2.*fDz*fDPhi*fRMin;
695  aFou = 2.*fDz*(fRMax-fRMin);
696  totArea = aOne+aTwo+2.*aFou;
697  phi = RandFlat::shoot(startPhi,endPhi);
698  cosphi = std::cos(phi);
699  sinphi = std::sin(phi);
700  rRand = fRMin + (fRMax-fRMin)*std::sqrt(RandFlat::shoot());
701 
702  if(startPhi == 0 && endPhi == twopi)
703  aFou = 0;
704 
705  chose = RandFlat::shoot(0.,aOne+aTwo+2.*aFou);
706  if( (chose >= 0) && (chose < aOne) )
707  {
708  xRand = fRMax*cosphi;
709  yRand = fRMax*sinphi;
710  zRand = RandFlat::shoot(-1.*fDz,fDz);
711  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
712  }
713  else if( (chose >= aOne) && (chose < aOne + aTwo) )
714  {
715  xRand = fRMin*cosphi;
716  yRand = fRMin*sinphi;
717  zRand = RandFlat::shoot(-1.*fDz,fDz);
718  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
719  }
720  else if( (chose >= aOne+aTwo) && (chose <aOne+aTwo+aFou) )
721  {
722  xRand = rRand*std::cos(fSPhi+fDPhi);
723  yRand = rRand*std::sin(fSPhi+fDPhi);
724  zRand = RandFlat::shoot(-1.*fDz,fDz);
725  return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
726  }
727 
728  // else
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 
737 //
738 // GetPointOnRing
739 //
740 // Auxiliary method for GetPoint On Surface
741 //
742 G4ThreeVector G4Polycone::GetPointOnRing(G4double fRMin1, G4double fRMax1,
743  G4double fRMin2,G4double fRMax2,
744  G4double zOne) const
745 {
746  G4double xRand,yRand,phi,cosphi,sinphi,rRand1,rRand2,A1,Atot,rCh;
747  phi = RandFlat::shoot(startPhi,endPhi);
748  cosphi = std::cos(phi);
749  sinphi = std::sin(phi);
750 
751  if(fRMin1==fRMin2)
752  {
753  rRand1 = fRMin1; A1=0.;
754  }
755  else
756  {
757  rRand1 = RandFlat::shoot(fRMin1,fRMin2);
758  A1=std::fabs(fRMin2*fRMin2-fRMin1*fRMin1);
759  }
760  if(fRMax1==fRMax2)
761  {
762  rRand2=fRMax1; Atot=A1;
763  }
764  else
765  {
766  rRand2 = RandFlat::shoot(fRMax1,fRMax2);
767  Atot = A1+std::fabs(fRMax2*fRMax2-fRMax1*fRMax1);
768  }
769  rCh = RandFlat::shoot(0.,Atot);
770 
771  if(rCh>A1) { rRand1=rRand2; }
772 
773  xRand = rRand1*cosphi;
774  yRand = rRand1*sinphi;
775 
776  return G4ThreeVector(xRand, yRand, zOne);
777 }
778 
779 
780 //
781 // GetPointOnCut
782 //
783 // Auxiliary method for Get Point On Surface
784 //
785 G4ThreeVector G4Polycone::GetPointOnCut(G4double fRMin1, G4double fRMax1,
786  G4double fRMin2, G4double fRMax2,
787  G4double zOne, G4double zTwo,
788  G4double& totArea) const
789 { if(zOne==zTwo)
790  {
791  return GetPointOnRing(fRMin1, fRMax1,fRMin2,fRMax2,zOne);
792  }
793  if( (fRMin1 == fRMin2) && (fRMax1 == fRMax2) )
794  {
795  return GetPointOnTubs(fRMin1, fRMax1,zOne,zTwo,totArea);
796  }
797  return GetPointOnCone(fRMin1,fRMax1,fRMin2,fRMax2,zOne,zTwo,totArea);
798 }
799 
800 
801 //
802 // GetPointOnSurface
803 //
804 G4ThreeVector G4Polycone::GetPointOnSurface() const
805 {
806  if (!genericPcon) // Polycone by faces
807  {
808  G4double Area=0,totArea=0,Achose1=0,Achose2=0,phi,cosphi,sinphi,rRand;
809  G4int i=0;
810  G4int numPlanes = original_parameters->Num_z_planes;
811 
812  phi = RandFlat::shoot(startPhi,endPhi);
813  cosphi = std::cos(phi);
814  sinphi = std::sin(phi);
815 
816  rRand = original_parameters->Rmin[0] +
817  ( (original_parameters->Rmax[0]-original_parameters->Rmin[0])
818  * std::sqrt(RandFlat::shoot()) );
819 
820  std::vector<G4double> areas; // (numPlanes+1);
821  std::vector<G4ThreeVector> points; // (numPlanes-1);
822 
823  areas.push_back(pi*(sqr(original_parameters->Rmax[0])
824  -sqr(original_parameters->Rmin[0])));
825 
826  for(i=0; i<numPlanes-1; i++)
827  {
828  Area = (original_parameters->Rmin[i]+original_parameters->Rmin[i+1])
829  * std::sqrt(sqr(original_parameters->Rmin[i]
830  -original_parameters->Rmin[i+1])+
831  sqr(original_parameters->Z_values[i+1]
832  -original_parameters->Z_values[i]));
833 
834  Area += (original_parameters->Rmax[i]+original_parameters->Rmax[i+1])
835  * std::sqrt(sqr(original_parameters->Rmax[i]
836  -original_parameters->Rmax[i+1])+
837  sqr(original_parameters->Z_values[i+1]
838  -original_parameters->Z_values[i]));
839 
840  Area *= 0.5*(endPhi-startPhi);
841 
842  if(startPhi==0.&& endPhi == twopi)
843  {
844  Area += std::fabs(original_parameters->Z_values[i+1]
845  -original_parameters->Z_values[i])*
846  (original_parameters->Rmax[i]
847  +original_parameters->Rmax[i+1]
848  -original_parameters->Rmin[i]
849  -original_parameters->Rmin[i+1]);
850  }
851  areas.push_back(Area);
852  totArea += Area;
853  }
854 
855  areas.push_back(pi*(sqr(original_parameters->Rmax[numPlanes-1])-
856  sqr(original_parameters->Rmin[numPlanes-1])));
857 
858  totArea += (areas[0]+areas[numPlanes]);
859  G4double chose = RandFlat::shoot(0.,totArea);
860 
861  if( (chose>=0.) && (chose<areas[0]) )
862  {
863  return G4ThreeVector(rRand*cosphi, rRand*sinphi,
864  original_parameters->Z_values[0]);
865  }
866 
867  for (i=0; i<numPlanes-1; i++)
868  {
869  Achose1 += areas[i];
870  Achose2 = (Achose1+areas[i+1]);
871  if(chose>=Achose1 && chose<Achose2)
872  {
873  return GetPointOnCut(original_parameters->Rmin[i],
874  original_parameters->Rmax[i],
875  original_parameters->Rmin[i+1],
876  original_parameters->Rmax[i+1],
877  original_parameters->Z_values[i],
878  original_parameters->Z_values[i+1], Area);
879  }
880  }
881 
882  rRand = original_parameters->Rmin[numPlanes-1] +
883  ( (original_parameters->Rmax[numPlanes-1]-original_parameters->Rmin[numPlanes-1])
884  * std::sqrt(RandFlat::shoot()) );
885 
886  return G4ThreeVector(rRand*cosphi,rRand*sinphi,
887  original_parameters->Z_values[numPlanes-1]);
888  }
889  else // Generic Polycone
890  {
891  return GetPointOnSurfaceGeneric();
892  }
893 }
894 
895 //
896 // CreatePolyhedron
897 //
898 G4Polyhedron* G4Polycone::CreatePolyhedron() const
899 {
900  //
901  // This has to be fixed in visualization. Fake it for the moment.
902  //
903  if (!genericPcon)
904  {
905  return new G4PolyhedronPcon( original_parameters->Start_angle,
906  original_parameters->Opening_angle,
907  original_parameters->Num_z_planes,
908  original_parameters->Z_values,
909  original_parameters->Rmin,
910  original_parameters->Rmax );
911  }
912  else
913  {
914  // The following code prepares for:
915  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
916  // const double xyz[][3],
917  // const int faces_vec[][4])
918  // Here is an extract from the header file HepPolyhedron.h:
936  const G4int numSide =
937  G4int(G4Polyhedron::GetNumberOfRotationSteps()
938  * (endPhi - startPhi) / twopi) + 1;
939  G4int nNodes;
940  G4int nFaces;
941  typedef G4double double3[3];
942  double3* xyz;
943  typedef G4int int4[4];
944  int4* faces_vec;
945  if (phiIsOpen)
946  {
947  // Triangulate open ends. Simple ear-chopping algorithm...
948  // I'm not sure how robust this algorithm is (J.Allison).
949  //
950  std::vector<G4bool> chopped(numCorner, false);
951  std::vector<G4int*> triQuads;
952  G4int remaining = numCorner;
953  G4int iStarter = 0;
954  while (remaining >= 3)
955  {
956  // Find unchopped corners...
957  //
958  G4int A = -1, B = -1, C = -1;
959  G4int iStepper = iStarter;
960  do
961  {
962  if (A < 0) { A = iStepper; }
963  else if (B < 0) { B = iStepper; }
964  else if (C < 0) { C = iStepper; }
965  do
966  {
967  if (++iStepper >= numCorner) { iStepper = 0; }
968  }
969  while (chopped[iStepper]);
970  }
971  while (C < 0 && iStepper != iStarter);
972 
973  // Check triangle at B is pointing outward (an "ear").
974  // Sign of z cross product determines...
975  //
976  G4double BAr = corners[A].r - corners[B].r;
977  G4double BAz = corners[A].z - corners[B].z;
978  G4double BCr = corners[C].r - corners[B].r;
979  G4double BCz = corners[C].z - corners[B].z;
980  if (BAr * BCz - BAz * BCr < kCarTolerance)
981  {
982  G4int* tq = new G4int[3];
983  tq[0] = A + 1;
984  tq[1] = B + 1;
985  tq[2] = C + 1;
986  triQuads.push_back(tq);
987  chopped[B] = true;
988  --remaining;
989  }
990  else
991  {
992  do
993  {
994  if (++iStarter >= numCorner) { iStarter = 0; }
995  }
996  while (chopped[iStarter]);
997  }
998  }
999  // Transfer to faces...
1000  //
1001  nNodes = (numSide + 1) * numCorner;
1002  nFaces = numSide * numCorner + 2 * triQuads.size();
1003  faces_vec = new int4[nFaces];
1004  G4int iface = 0;
1005  G4int addition = numCorner * numSide;
1006  G4int d = numCorner - 1;
1007  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
1008  {
1009  for (size_t i = 0; i < triQuads.size(); ++i)
1010  {
1011  // Negative for soft/auxiliary/normally invisible edges...
1012  //
1013  G4int a, b, c;
1014  if (iEnd == 0)
1015  {
1016  a = triQuads[i][0];
1017  b = triQuads[i][1];
1018  c = triQuads[i][2];
1019  }
1020  else
1021  {
1022  a = triQuads[i][0] + addition;
1023  b = triQuads[i][2] + addition;
1024  c = triQuads[i][1] + addition;
1025  }
1026  G4int ab = std::abs(b - a);
1027  G4int bc = std::abs(c - b);
1028  G4int ca = std::abs(a - c);
1029  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
1030  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
1031  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
1032  faces_vec[iface][3] = 0;
1033  ++iface;
1034  }
1035  }
1036 
1037  // Continue with sides...
1038 
1039  xyz = new double3[nNodes];
1040  const G4double dPhi = (endPhi - startPhi) / numSide;
1041  G4double phi = startPhi;
1042  G4int ixyz = 0;
1043  for (G4int iSide = 0; iSide < numSide; ++iSide)
1044  {
1045  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1046  {
1047  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1048  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1049  xyz[ixyz][2] = corners[iCorner].z;
1050  if (iSide == 0) // startPhi
1051  {
1052  if (iCorner < numCorner - 1)
1053  {
1054  faces_vec[iface][0] = ixyz + 1;
1055  faces_vec[iface][1] = -(ixyz + numCorner + 1);
1056  faces_vec[iface][2] = ixyz + numCorner + 2;
1057  faces_vec[iface][3] = ixyz + 2;
1058  }
1059  else
1060  {
1061  faces_vec[iface][0] = ixyz + 1;
1062  faces_vec[iface][1] = -(ixyz + numCorner + 1);
1063  faces_vec[iface][2] = ixyz + 2;
1064  faces_vec[iface][3] = ixyz - numCorner + 2;
1065  }
1066  }
1067  else if (iSide == numSide - 1) // endPhi
1068  {
1069  if (iCorner < numCorner - 1)
1070  {
1071  faces_vec[iface][0] = ixyz + 1;
1072  faces_vec[iface][1] = ixyz + numCorner + 1;
1073  faces_vec[iface][2] = ixyz + numCorner + 2;
1074  faces_vec[iface][3] = -(ixyz + 2);
1075  }
1076  else
1077  {
1078  faces_vec[iface][0] = ixyz + 1;
1079  faces_vec[iface][1] = ixyz + numCorner + 1;
1080  faces_vec[iface][2] = ixyz + 2;
1081  faces_vec[iface][3] = -(ixyz - numCorner + 2);
1082  }
1083  }
1084  else
1085  {
1086  if (iCorner < numCorner - 1)
1087  {
1088  faces_vec[iface][0] = ixyz + 1;
1089  faces_vec[iface][1] = -(ixyz + numCorner + 1);
1090  faces_vec[iface][2] = ixyz + numCorner + 2;
1091  faces_vec[iface][3] = -(ixyz + 2);
1092  }
1093  else
1094  {
1095  faces_vec[iface][0] = ixyz + 1;
1096  faces_vec[iface][1] = -(ixyz + numCorner + 1);
1097  faces_vec[iface][2] = ixyz + 2;
1098  faces_vec[iface][3] = -(ixyz - numCorner + 2);
1099  }
1100  }
1101  ++iface;
1102  ++ixyz;
1103  }
1104  phi += dPhi;
1105  }
1106 
1107  // Last corners...
1108 
1109  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1110  {
1111  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1112  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1113  xyz[ixyz][2] = corners[iCorner].z;
1114  ++ixyz;
1115  }
1116  }
1117  else // !phiIsOpen - i.e., a complete 360 degrees.
1118  {
1119  nNodes = numSide * numCorner;
1120  nFaces = numSide * numCorner;;
1121  xyz = new double3[nNodes];
1122  faces_vec = new int4[nFaces];
1123  const G4double dPhi = (endPhi - startPhi) / numSide;
1124  G4double phi = startPhi;
1125  G4int ixyz = 0, iface = 0;
1126  for (G4int iSide = 0; iSide < numSide; ++iSide)
1127  {
1128  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
1129  {
1130  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
1131  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
1132  xyz[ixyz][2] = corners[iCorner].z;
1133 
1134  if (iSide < numSide - 1)
1135  {
1136  if (iCorner < numCorner - 1)
1137  {
1138  faces_vec[iface][0] = ixyz + 1;
1139  faces_vec[iface][1] = -(ixyz + numCorner + 1);
1140  faces_vec[iface][2] = ixyz + numCorner + 2;
1141  faces_vec[iface][3] = -(ixyz + 2);
1142  }
1143  else
1144  {
1145  faces_vec[iface][0] = ixyz + 1;
1146  faces_vec[iface][1] = -(ixyz + numCorner + 1);
1147  faces_vec[iface][2] = ixyz + 2;
1148  faces_vec[iface][3] = -(ixyz - numCorner + 2);
1149  }
1150  }
1151  else // Last side joins ends...
1152  {
1153  if (iCorner < numCorner - 1)
1154  {
1155  faces_vec[iface][0] = ixyz + 1;
1156  faces_vec[iface][1] = -(ixyz + numCorner - nFaces + 1);
1157  faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
1158  faces_vec[iface][3] = -(ixyz + 2);
1159  }
1160  else
1161  {
1162  faces_vec[iface][0] = ixyz + 1;
1163  faces_vec[iface][1] = -(ixyz - nFaces + numCorner + 1);
1164  faces_vec[iface][2] = ixyz - nFaces + 2;
1165  faces_vec[iface][3] = -(ixyz - numCorner + 2);
1166  }
1167  }
1168  ++ixyz;
1169  ++iface;
1170  }
1171  phi += dPhi;
1172  }
1173  }
1174  G4Polyhedron* polyhedron = new G4Polyhedron;
1175  G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
1176  delete [] faces_vec;
1177  delete [] xyz;
1178  if (problem)
1179  {
1180  std::ostringstream message;
1181  message << "Problem creating G4Polyhedron for: " << GetName();
1182  G4Exception("G4Polycone::CreatePolyhedron()", "GeomSolids1002",
1183  JustWarning, message);
1184  delete polyhedron;
1185  return 0;
1186  }
1187  else
1188  {
1189  return polyhedron;
1190  }
1191  }
1192 }
1193 
1194 //
1195 // CreateNURBS
1196 //
1197 G4NURBS *G4Polycone::CreateNURBS() const
1198 {
1199  return 0;
1200 }
1201 
1202 //
1203 // G4PolyconeHistorical stuff
1204 //
1205 
1206 G4PolyconeHistorical::G4PolyconeHistorical()
1207  : Start_angle(0.), Opening_angle(0.), Num_z_planes(0),
1208  Z_values(0), Rmin(0), Rmax(0)
1209 {
1210 }
1211 
1212 G4PolyconeHistorical::~G4PolyconeHistorical()
1213 {
1214  delete [] Z_values;
1215  delete [] Rmin;
1216  delete [] Rmax;
1217 }
1218 
1219 G4PolyconeHistorical::
1220 G4PolyconeHistorical( const G4PolyconeHistorical &source )
1221 {
1222  Start_angle = source.Start_angle;
1223  Opening_angle = source.Opening_angle;
1224  Num_z_planes = source.Num_z_planes;
1225 
1226  Z_values = new G4double[Num_z_planes];
1227  Rmin = new G4double[Num_z_planes];
1228  Rmax = new G4double[Num_z_planes];
1229 
1230  for( G4int i = 0; i < Num_z_planes; i++)
1231  {
1232  Z_values[i] = source.Z_values[i];
1233  Rmin[i] = source.Rmin[i];
1234  Rmax[i] = source.Rmax[i];
1235  }
1236 }
1237 
1238 G4PolyconeHistorical&
1239 G4PolyconeHistorical::operator=( const G4PolyconeHistorical& right )
1240 {
1241  if ( &right == this ) return *this;
1242 
1243  if (&right)
1244  {
1245  Start_angle = right.Start_angle;
1246  Opening_angle = right.Opening_angle;
1247  Num_z_planes = right.Num_z_planes;
1248 
1249  delete [] Z_values;
1250  delete [] Rmin;
1251  delete [] Rmax;
1252  Z_values = new G4double[Num_z_planes];
1253  Rmin = new G4double[Num_z_planes];
1254  Rmax = new G4double[Num_z_planes];
1255 
1256  for( G4int i = 0; i < Num_z_planes; i++)
1257  {
1258  Z_values[i] = right.Z_values[i];
1259  Rmin[i] = right.Rmin[i];
1260  Rmax[i] = right.Rmax[i];
1261  }
1262  }
1263  return *this;
1264 }
Generated on Sat May 25 2013 14:33:15 for Geant4 by   doxygen 1.8.4