source: trunk/source/geometry/solids/specific/src/G4Polycone.cc @ 1350

Last change on this file since 1350 was 1337, checked in by garnier, 14 years ago

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