G4Polycone.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id: G4Polycone.cc 69790 2013-05-15 12:39:10Z gcosmo $
00028 //
00029 // 
00030 // --------------------------------------------------------------------
00031 // GEANT 4 class source file
00032 //
00033 //
00034 // G4Polycone.cc
00035 //
00036 // Implementation of a CSG polycone
00037 //
00038 // --------------------------------------------------------------------
00039 
00040 #include "G4Polycone.hh"
00041 
00042 #include "G4PolyconeSide.hh"
00043 #include "G4PolyPhiFace.hh"
00044 
00045 #include "Randomize.hh"
00046 
00047 #include "G4Polyhedron.hh"
00048 #include "G4EnclosingCylinder.hh"
00049 #include "G4ReduciblePolygon.hh"
00050 #include "G4VPVParameterisation.hh"
00051 
00052 using namespace CLHEP;
00053 
00054 //
00055 // Constructor (GEANT3 style parameters)
00056 //  
00057 G4Polycone::G4Polycone( const G4String& name, 
00058                               G4double phiStart,
00059                               G4double phiTotal,
00060                               G4int numZPlanes,
00061                         const G4double zPlane[],
00062                         const G4double rInner[],
00063                         const G4double rOuter[]  )
00064   : G4VCSGfaceted( name ), genericPcon(false)
00065 {
00066   //
00067   // Some historical ugliness
00068   //
00069   original_parameters = new G4PolyconeHistorical();
00070   
00071   original_parameters->Start_angle = phiStart;
00072   original_parameters->Opening_angle = phiTotal;
00073   original_parameters->Num_z_planes = numZPlanes;
00074   original_parameters->Z_values = new G4double[numZPlanes];
00075   original_parameters->Rmin = new G4double[numZPlanes];
00076   original_parameters->Rmax = new G4double[numZPlanes];
00077 
00078   G4int i;
00079   for (i=0; i<numZPlanes; i++)
00080   { 
00081     if(rInner[i]>rOuter[i])
00082     {
00083       DumpInfo();
00084       std::ostringstream message;
00085       message << "Cannot create a Polycone with rInner > rOuter for the same Z"
00086               << G4endl
00087               << "        rInner > rOuter for the same Z !" << G4endl
00088               << "        rMin[" << i << "] = " << rInner[i]
00089               << " -- rMax[" << i << "] = " << rOuter[i];
00090       G4Exception("G4Polycone::G4Polycone()", "GeomSolids0002",
00091                   FatalErrorInArgument, message);
00092     }
00093     if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
00094     {
00095       if( (rInner[i]   > rOuter[i+1])
00096         ||(rInner[i+1] > rOuter[i])   )
00097       {
00098         DumpInfo();
00099         std::ostringstream message;
00100         message << "Cannot create a Polycone with no contiguous segments."
00101                 << G4endl
00102                 << "        Segments are not contiguous !" << G4endl
00103                 << "        rMin[" << i << "] = " << rInner[i]
00104                 << " -- rMax[" << i+1 << "] = " << rOuter[i+1] << G4endl
00105                 << "        rMin[" << i+1 << "] = " << rInner[i+1]
00106                 << " -- rMax[" << i << "] = " << rOuter[i];
00107         G4Exception("G4Polycone::G4Polycone()", "GeomSolids0002",
00108                     FatalErrorInArgument, message);
00109       }
00110     } 
00111     original_parameters->Z_values[i] = zPlane[i];
00112     original_parameters->Rmin[i] = rInner[i];
00113     original_parameters->Rmax[i] = rOuter[i];
00114   }
00115 
00116   //
00117   // Build RZ polygon using special PCON/PGON GEANT3 constructor
00118   //
00119   G4ReduciblePolygon *rz =
00120     new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
00121   
00122   //
00123   // Do the real work
00124   //
00125   Create( phiStart, phiTotal, rz );
00126   
00127   delete rz;
00128 }
00129 
00130 
00131 //
00132 // Constructor (generic parameters)
00133 //
00134 G4Polycone::G4Polycone( const G4String& name, 
00135                               G4double phiStart,
00136                               G4double phiTotal,
00137                               G4int    numRZ,
00138                         const G4double r[],
00139                         const G4double z[]   )
00140   : G4VCSGfaceted( name ), genericPcon(true)
00141 { 
00142   
00143   G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
00144   
00145   Create( phiStart, phiTotal, rz );
00146   
00147   // Set original_parameters struct for consistency
00148   //
00149   SetOriginalParameters();
00150 
00151   delete rz;
00152 }
00153 
00154 
00155 //
00156 // Create
00157 //
00158 // Generic create routine, called by each constructor after
00159 // conversion of arguments
00160 //
00161 void G4Polycone::Create( G4double phiStart,
00162                          G4double phiTotal,
00163                          G4ReduciblePolygon *rz    )
00164 {
00165   //
00166   // Perform checks of rz values
00167   //
00168   if (rz->Amin() < 0.0)
00169   {
00170     std::ostringstream message;
00171     message << "Illegal input parameters - " << GetName() << G4endl
00172             << "        All R values must be >= 0 !";
00173     G4Exception("G4Polycone::Create()", "GeomSolids0002",
00174                 FatalErrorInArgument, message);
00175   }
00176     
00177   G4double rzArea = rz->Area();
00178   if (rzArea < -kCarTolerance)
00179     rz->ReverseOrder();
00180 
00181   else if (rzArea < -kCarTolerance)
00182   {
00183     std::ostringstream message;
00184     message << "Illegal input parameters - " << GetName() << G4endl
00185             << "        R/Z cross section is zero or near zero: " << rzArea;
00186     G4Exception("G4Polycone::Create()", "GeomSolids0002",
00187                 FatalErrorInArgument, message);
00188   }
00189     
00190   if ( (!rz->RemoveDuplicateVertices( kCarTolerance ))
00191     || (!rz->RemoveRedundantVertices( kCarTolerance ))     ) 
00192   {
00193     std::ostringstream message;
00194     message << "Illegal input parameters - " << GetName() << G4endl
00195             << "        Too few unique R/Z values !";
00196     G4Exception("G4Polycone::Create()", "GeomSolids0002",
00197                 FatalErrorInArgument, message);
00198   }
00199 
00200   if (rz->CrossesItself(1/kInfinity)) 
00201   {
00202     std::ostringstream message;
00203     message << "Illegal input parameters - " << GetName() << G4endl
00204             << "        R/Z segments cross !";
00205     G4Exception("G4Polycone::Create()", "GeomSolids0002",
00206                 FatalErrorInArgument, message);
00207   }
00208 
00209   numCorner = rz->NumVertices();
00210 
00211   //
00212   // Phi opening? Account for some possible roundoff, and interpret
00213   // nonsense value as representing no phi opening
00214   //
00215   if (phiTotal <= 0 || phiTotal > twopi-1E-10)
00216   {
00217     phiIsOpen = false;
00218     startPhi = 0;
00219     endPhi = twopi;
00220   }
00221   else
00222   {
00223     phiIsOpen = true;
00224     
00225     //
00226     // Convert phi into our convention
00227     //
00228     startPhi = phiStart;
00229     while( startPhi < 0 ) startPhi += twopi;
00230     
00231     endPhi = phiStart+phiTotal;
00232     while( endPhi < startPhi ) endPhi += twopi;
00233   }
00234   
00235   //
00236   // Allocate corner array. 
00237   //
00238   corners = new G4PolyconeSideRZ[numCorner];
00239 
00240   //
00241   // Copy corners
00242   //
00243   G4ReduciblePolygonIterator iterRZ(rz);
00244   
00245   G4PolyconeSideRZ *next = corners;
00246   iterRZ.Begin();
00247   do
00248   {
00249     next->r = iterRZ.GetA();
00250     next->z = iterRZ.GetB();
00251   } while( ++next, iterRZ.Next() );
00252   
00253   //
00254   // Allocate face pointer array
00255   //
00256   numFace = phiIsOpen ? numCorner+2 : numCorner;
00257   faces = new G4VCSGface*[numFace];
00258   
00259   //
00260   // Construct conical faces
00261   //
00262   // But! Don't construct a face if both points are at zero radius!
00263   //
00264   G4PolyconeSideRZ *corner = corners,
00265                    *prev = corners + numCorner-1,
00266                    *nextNext;
00267   G4VCSGface  **face = faces;
00268   do
00269   {
00270     next = corner+1;
00271     if (next >= corners+numCorner) next = corners;
00272     nextNext = next+1;
00273     if (nextNext >= corners+numCorner) nextNext = corners;
00274     
00275     if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
00276     
00277     //
00278     // We must decide here if we can dare declare one of our faces
00279     // as having a "valid" normal (i.e. allBehind = true). This
00280     // is never possible if the face faces "inward" in r.
00281     //
00282     G4bool allBehind;
00283     if (corner->z > next->z)
00284     {
00285       allBehind = false;
00286     }
00287     else
00288     {
00289       //
00290       // Otherwise, it is only true if the line passing
00291       // through the two points of the segment do not
00292       // split the r/z cross section
00293       //
00294       allBehind = !rz->BisectedBy( corner->r, corner->z,
00295                  next->r, next->z, kCarTolerance );
00296     }
00297     
00298     *face++ = new G4PolyconeSide( prev, corner, next, nextNext,
00299                 startPhi, endPhi-startPhi, phiIsOpen, allBehind );
00300   } while( prev=corner, corner=next, corner > corners );
00301   
00302   if (phiIsOpen)
00303   {
00304     //
00305     // Construct phi open edges
00306     //
00307     *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi  );
00308     *face++ = new G4PolyPhiFace( rz, endPhi,   0, startPhi );
00309   }
00310   
00311   //
00312   // We might have dropped a face or two: recalculate numFace
00313   //
00314   numFace = face-faces;
00315   
00316   //
00317   // Make enclosingCylinder
00318   //
00319   enclosingCylinder =
00320     new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
00321 }
00322 
00323 
00324 //
00325 // Fake default constructor - sets only member data and allocates memory
00326 //                            for usage restricted to object persistency.
00327 //
00328 G4Polycone::G4Polycone( __void__& a )
00329   : G4VCSGfaceted(a), startPhi(0.),  endPhi(0.), phiIsOpen(false),
00330     genericPcon(false), numCorner(0), corners(0),
00331     original_parameters(0), enclosingCylinder(0)
00332 {
00333 }
00334 
00335 
00336 //
00337 // Destructor
00338 //
00339 G4Polycone::~G4Polycone()
00340 {
00341   delete [] corners;
00342   delete original_parameters;
00343   delete enclosingCylinder;
00344 }
00345 
00346 
00347 //
00348 // Copy constructor
00349 //
00350 G4Polycone::G4Polycone( const G4Polycone &source )
00351   : G4VCSGfaceted( source )
00352 {
00353   CopyStuff( source );
00354 }
00355 
00356 
00357 //
00358 // Assignment operator
00359 //
00360 const G4Polycone &G4Polycone::operator=( const G4Polycone &source )
00361 {
00362   if (this == &source) return *this;
00363   
00364   G4VCSGfaceted::operator=( source );
00365   
00366   delete [] corners;
00367   if (original_parameters) delete original_parameters;
00368   
00369   delete enclosingCylinder;
00370   
00371   CopyStuff( source );
00372   
00373   return *this;
00374 }
00375 
00376 
00377 //
00378 // CopyStuff
00379 //
00380 void G4Polycone::CopyStuff( const G4Polycone &source )
00381 {
00382   //
00383   // Simple stuff
00384   //
00385   startPhi  = source.startPhi;
00386   endPhi    = source.endPhi;
00387   phiIsOpen  = source.phiIsOpen;
00388   numCorner  = source.numCorner;
00389   genericPcon= source.genericPcon;
00390 
00391   //
00392   // The corner array
00393   //
00394   corners = new G4PolyconeSideRZ[numCorner];
00395   
00396   G4PolyconeSideRZ  *corn = corners,
00397         *sourceCorn = source.corners;
00398   do
00399   {
00400     *corn = *sourceCorn;
00401   } while( ++sourceCorn, ++corn < corners+numCorner );
00402   
00403   //
00404   // Original parameters
00405   //
00406   if (source.original_parameters)
00407   {
00408     original_parameters =
00409       new G4PolyconeHistorical( *source.original_parameters );
00410   }
00411   
00412   //
00413   // Enclosing cylinder
00414   //
00415   enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder );
00416 }
00417 
00418 
00419 //
00420 // Reset
00421 //
00422 G4bool G4Polycone::Reset()
00423 {
00424   if (genericPcon)
00425   {
00426     std::ostringstream message;
00427     message << "Solid " << GetName() << " built using generic construct."
00428             << G4endl << "Not applicable to the generic construct !";
00429     G4Exception("G4Polycone::Reset()", "GeomSolids1001",
00430                 JustWarning, message, "Parameters NOT resetted.");
00431     return 1;
00432   }
00433 
00434   //
00435   // Clear old setup
00436   //
00437   G4VCSGfaceted::DeleteStuff();
00438   delete [] corners;
00439   delete enclosingCylinder;
00440 
00441   //
00442   // Rebuild polycone
00443   //
00444   G4ReduciblePolygon *rz =
00445     new G4ReduciblePolygon( original_parameters->Rmin,
00446                             original_parameters->Rmax,
00447                             original_parameters->Z_values,
00448                             original_parameters->Num_z_planes );
00449   Create( original_parameters->Start_angle,
00450           original_parameters->Opening_angle, rz );
00451   delete rz;
00452 
00453   return 0;
00454 }
00455 
00456 
00457 //
00458 // Inside
00459 //
00460 // This is an override of G4VCSGfaceted::Inside, created in order
00461 // to speed things up by first checking with G4EnclosingCylinder.
00462 //
00463 EInside G4Polycone::Inside( const G4ThreeVector &p ) const
00464 {
00465   //
00466   // Quick test
00467   //
00468   if (enclosingCylinder->MustBeOutside(p)) return kOutside;
00469 
00470   //
00471   // Long answer
00472   //
00473   return G4VCSGfaceted::Inside(p);
00474 }
00475 
00476 
00477 //
00478 // DistanceToIn
00479 //
00480 // This is an override of G4VCSGfaceted::Inside, created in order
00481 // to speed things up by first checking with G4EnclosingCylinder.
00482 //
00483 G4double G4Polycone::DistanceToIn( const G4ThreeVector &p,
00484                                    const G4ThreeVector &v ) const
00485 {
00486   //
00487   // Quick test
00488   //
00489   if (enclosingCylinder->ShouldMiss(p,v))
00490     return kInfinity;
00491   
00492   //
00493   // Long answer
00494   //
00495   return G4VCSGfaceted::DistanceToIn( p, v );
00496 }
00497 
00498 
00499 //
00500 // DistanceToIn
00501 //
00502 G4double G4Polycone::DistanceToIn( const G4ThreeVector &p ) const
00503 {
00504   return G4VCSGfaceted::DistanceToIn(p);
00505 }
00506 
00507 
00508 //
00509 // ComputeDimensions
00510 //
00511 void G4Polycone::ComputeDimensions(       G4VPVParameterisation* p,
00512                                     const G4int n,
00513                                     const G4VPhysicalVolume* pRep )
00514 {
00515   p->ComputeDimensions(*this,n,pRep);
00516 }
00517 
00518 //
00519 // GetEntityType
00520 //
00521 G4GeometryType  G4Polycone::GetEntityType() const
00522 {
00523   return G4String("G4Polycone");
00524 }
00525 
00526 //
00527 // Make a clone of the object
00528 //
00529 G4VSolid* G4Polycone::Clone() const
00530 {
00531   return new G4Polycone(*this);
00532 }
00533 
00534 //
00535 // Stream object contents to an output stream
00536 //
00537 std::ostream& G4Polycone::StreamInfo( std::ostream& os ) const
00538 {
00539   G4int oldprc = os.precision(16);
00540   os << "-----------------------------------------------------------\n"
00541      << "    *** Dump for solid - " << GetName() << " ***\n"
00542      << "    ===================================================\n"
00543      << " Solid type: G4Polycone\n"
00544      << " Parameters: \n"
00545      << "    starting phi angle : " << startPhi/degree << " degrees \n"
00546      << "    ending phi angle   : " << endPhi/degree << " degrees \n";
00547   G4int i=0;
00548   if (!genericPcon)
00549   {
00550     G4int numPlanes = original_parameters->Num_z_planes;
00551     os << "    number of Z planes: " << numPlanes << "\n"
00552        << "              Z values: \n";
00553     for (i=0; i<numPlanes; i++)
00554     {
00555       os << "              Z plane " << i << ": "
00556          << original_parameters->Z_values[i] << "\n";
00557     }
00558     os << "              Tangent distances to inner surface (Rmin): \n";
00559     for (i=0; i<numPlanes; i++)
00560     {
00561       os << "              Z plane " << i << ": "
00562          << original_parameters->Rmin[i] << "\n";
00563     }
00564     os << "              Tangent distances to outer surface (Rmax): \n";
00565     for (i=0; i<numPlanes; i++)
00566     {
00567       os << "              Z plane " << i << ": "
00568          << original_parameters->Rmax[i] << "\n";
00569     }
00570   }
00571   os << "    number of RZ points: " << numCorner << "\n"
00572      << "              RZ values (corners): \n";
00573      for (i=0; i<numCorner; i++)
00574      {
00575        os << "                         "
00576           << corners[i].r << ", " << corners[i].z << "\n";
00577      }
00578   os << "-----------------------------------------------------------\n";
00579   os.precision(oldprc);
00580 
00581   return os;
00582 }
00583 
00584 
00585 //
00586 // GetPointOnCone
00587 //
00588 // Auxiliary method for Get Point On Surface
00589 //
00590 G4ThreeVector G4Polycone::GetPointOnCone(G4double fRmin1, G4double fRmax1,
00591                                          G4double fRmin2, G4double fRmax2,
00592                                          G4double zOne,   G4double zTwo,
00593                                          G4double& totArea) const
00594 { 
00595   // declare working variables
00596   //
00597   G4double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
00598   G4double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
00599   G4double fDz=(zTwo-zOne)/2., afDz=std::fabs(fDz);
00600   G4ThreeVector point, offset=G4ThreeVector(0.,0.,0.5*(zTwo+zOne));
00601   fDPhi = endPhi - startPhi;
00602   rone = (fRmax1-fRmax2)/(2.*fDz); 
00603   rtwo = (fRmin1-fRmin2)/(2.*fDz);
00604   if(fRmax1==fRmax2){qone=0.;}
00605   else{
00606     qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
00607   }
00608   if(fRmin1==fRmin2){qtwo=0.;}
00609   else{
00610     qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
00611    }
00612   Aone   = 0.5*fDPhi*(fRmax2 + fRmax1)*(sqr(fRmin1-fRmin2)+sqr(zTwo-zOne));       
00613   Atwo   = 0.5*fDPhi*(fRmin2 + fRmin1)*(sqr(fRmax1-fRmax2)+sqr(zTwo-zOne));
00614   Afive  = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
00615   totArea = Aone+Atwo+2.*Afive;
00616   
00617   phi  = RandFlat::shoot(startPhi,endPhi);
00618   cosu = std::cos(phi);
00619   sinu = std::sin(phi);
00620 
00621 
00622   if( (startPhi == 0) && (endPhi == twopi) ) { Afive = 0; }
00623   chose = RandFlat::shoot(0.,Aone+Atwo+2.*Afive);
00624   if( (chose >= 0) && (chose < Aone) )
00625   {
00626     if(fRmax1 != fRmax2)
00627     {
00628       zRand = RandFlat::shoot(-1.*afDz,afDz); 
00629       point = G4ThreeVector (rone*cosu*(qone-zRand),
00630                              rone*sinu*(qone-zRand), zRand);
00631     }
00632     else
00633     {
00634       point = G4ThreeVector(fRmax1*cosu, fRmax1*sinu,
00635                             RandFlat::shoot(-1.*afDz,afDz));
00636      
00637     }
00638   }
00639   else if(chose >= Aone && chose < Aone + Atwo)
00640   {
00641     if(fRmin1 != fRmin2)
00642       { 
00643       zRand = RandFlat::shoot(-1.*afDz,afDz); 
00644       point = G4ThreeVector (rtwo*cosu*(qtwo-zRand),
00645                              rtwo*sinu*(qtwo-zRand), zRand);
00646       
00647     }
00648     else
00649     {
00650       point = G4ThreeVector(fRmin1*cosu, fRmin1*sinu,
00651                             RandFlat::shoot(-1.*afDz,afDz));
00652     }
00653   }
00654   else if( (chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive) )
00655   {
00656     zRand  = RandFlat::shoot(-1.*afDz,afDz);
00657     rmin   = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
00658     rmax   = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
00659     rRand1 = std::sqrt(RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
00660     point  = G4ThreeVector (rRand1*std::cos(startPhi),
00661                             rRand1*std::sin(startPhi), zRand);
00662   }
00663   else
00664   { 
00665     zRand  = RandFlat::shoot(-1.*afDz,afDz); 
00666     rmin   = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
00667     rmax   = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
00668     rRand1 = std::sqrt(RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
00669     point  = G4ThreeVector (rRand1*std::cos(endPhi),
00670                             rRand1*std::sin(endPhi), zRand);
00671    
00672   }
00673   
00674   return point+offset;
00675 }
00676 
00677 
00678 //
00679 // GetPointOnTubs
00680 //
00681 // Auxiliary method for GetPoint On Surface
00682 //
00683 G4ThreeVector G4Polycone::GetPointOnTubs(G4double fRMin, G4double fRMax,
00684                                          G4double zOne,  G4double zTwo,
00685                                          G4double& totArea) const
00686 { 
00687   G4double xRand,yRand,zRand,phi,cosphi,sinphi,chose,
00688            aOne,aTwo,aFou,rRand,fDz,fSPhi,fDPhi;
00689   fDz = std::fabs(0.5*(zTwo-zOne));
00690   fSPhi = startPhi;
00691   fDPhi = endPhi-startPhi;
00692   
00693   aOne = 2.*fDz*fDPhi*fRMax;
00694   aTwo = 2.*fDz*fDPhi*fRMin;
00695   aFou = 2.*fDz*(fRMax-fRMin);
00696   totArea = aOne+aTwo+2.*aFou;
00697   phi    = RandFlat::shoot(startPhi,endPhi);
00698   cosphi = std::cos(phi);
00699   sinphi = std::sin(phi);
00700   rRand  = fRMin + (fRMax-fRMin)*std::sqrt(RandFlat::shoot());
00701  
00702   if(startPhi == 0 && endPhi == twopi) 
00703     aFou = 0;
00704   
00705   chose  = RandFlat::shoot(0.,aOne+aTwo+2.*aFou);
00706   if( (chose >= 0) && (chose < aOne) )
00707   {
00708     xRand = fRMax*cosphi;
00709     yRand = fRMax*sinphi;
00710     zRand = RandFlat::shoot(-1.*fDz,fDz);
00711     return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
00712   }
00713   else if( (chose >= aOne) && (chose < aOne + aTwo) )
00714   {
00715     xRand = fRMin*cosphi;
00716     yRand = fRMin*sinphi;
00717     zRand = RandFlat::shoot(-1.*fDz,fDz);
00718     return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
00719   }
00720   else if( (chose >= aOne+aTwo) && (chose <aOne+aTwo+aFou) )
00721   {
00722     xRand = rRand*std::cos(fSPhi+fDPhi);
00723     yRand = rRand*std::sin(fSPhi+fDPhi);
00724     zRand = RandFlat::shoot(-1.*fDz,fDz);
00725     return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
00726   }
00727 
00728   // else
00729 
00730   xRand = rRand*std::cos(fSPhi+fDPhi);
00731   yRand = rRand*std::sin(fSPhi+fDPhi);
00732   zRand = RandFlat::shoot(-1.*fDz,fDz);
00733   return G4ThreeVector(xRand, yRand, zRand+0.5*(zTwo+zOne));
00734 }
00735 
00736 
00737 //
00738 // GetPointOnRing
00739 //
00740 // Auxiliary method for GetPoint On Surface
00741 //
00742 G4ThreeVector G4Polycone::GetPointOnRing(G4double fRMin1, G4double fRMax1,
00743                                          G4double fRMin2,G4double fRMax2,
00744                                          G4double zOne) const
00745 {
00746   G4double xRand,yRand,phi,cosphi,sinphi,rRand1,rRand2,A1,Atot,rCh;
00747   phi    = RandFlat::shoot(startPhi,endPhi);
00748   cosphi = std::cos(phi);
00749   sinphi = std::sin(phi);
00750 
00751   if(fRMin1==fRMin2)
00752   {
00753     rRand1 = fRMin1; A1=0.;
00754   }
00755   else
00756   {
00757     rRand1 = RandFlat::shoot(fRMin1,fRMin2);
00758     A1=std::fabs(fRMin2*fRMin2-fRMin1*fRMin1);
00759   }
00760   if(fRMax1==fRMax2)
00761   {
00762     rRand2=fRMax1; Atot=A1;
00763   }
00764   else
00765   {
00766     rRand2 = RandFlat::shoot(fRMax1,fRMax2);
00767     Atot   = A1+std::fabs(fRMax2*fRMax2-fRMax1*fRMax1);
00768   }
00769   rCh   = RandFlat::shoot(0.,Atot);
00770  
00771   if(rCh>A1) { rRand1=rRand2; }
00772   
00773   xRand = rRand1*cosphi;
00774   yRand = rRand1*sinphi;
00775 
00776   return G4ThreeVector(xRand, yRand, zOne);
00777 }
00778 
00779 
00780 //
00781 // GetPointOnCut
00782 //
00783 // Auxiliary method for Get Point On Surface
00784 //
00785 G4ThreeVector G4Polycone::GetPointOnCut(G4double fRMin1, G4double fRMax1,
00786                                         G4double fRMin2, G4double fRMax2,
00787                                         G4double zOne,  G4double zTwo,
00788                                         G4double& totArea) const
00789 {   if(zOne==zTwo)
00790     {
00791       return GetPointOnRing(fRMin1, fRMax1,fRMin2,fRMax2,zOne);
00792     }
00793     if( (fRMin1 == fRMin2) && (fRMax1 == fRMax2) )
00794     {
00795       return GetPointOnTubs(fRMin1, fRMax1,zOne,zTwo,totArea);
00796     }
00797     return GetPointOnCone(fRMin1,fRMax1,fRMin2,fRMax2,zOne,zTwo,totArea);
00798 }
00799 
00800 
00801 //
00802 // GetPointOnSurface
00803 //
00804 G4ThreeVector G4Polycone::GetPointOnSurface() const
00805 {
00806   if (!genericPcon)  // Polycone by faces
00807   {
00808     G4double Area=0,totArea=0,Achose1=0,Achose2=0,phi,cosphi,sinphi,rRand;
00809     G4int i=0;
00810     G4int numPlanes = original_parameters->Num_z_planes;
00811   
00812     phi = RandFlat::shoot(startPhi,endPhi);
00813     cosphi = std::cos(phi);
00814     sinphi = std::sin(phi);
00815 
00816     rRand = original_parameters->Rmin[0] +
00817       ( (original_parameters->Rmax[0]-original_parameters->Rmin[0])
00818         * std::sqrt(RandFlat::shoot()) );
00819 
00820     std::vector<G4double> areas;       // (numPlanes+1);
00821     std::vector<G4ThreeVector> points; // (numPlanes-1);
00822   
00823     areas.push_back(pi*(sqr(original_parameters->Rmax[0])
00824                        -sqr(original_parameters->Rmin[0])));
00825 
00826     for(i=0; i<numPlanes-1; i++)
00827     {
00828       Area = (original_parameters->Rmin[i]+original_parameters->Rmin[i+1])
00829            * std::sqrt(sqr(original_parameters->Rmin[i]
00830                       -original_parameters->Rmin[i+1])+
00831                        sqr(original_parameters->Z_values[i+1]
00832                       -original_parameters->Z_values[i]));
00833 
00834       Area += (original_parameters->Rmax[i]+original_parameters->Rmax[i+1])
00835             * std::sqrt(sqr(original_parameters->Rmax[i]
00836                        -original_parameters->Rmax[i+1])+
00837                         sqr(original_parameters->Z_values[i+1]
00838                        -original_parameters->Z_values[i]));
00839 
00840       Area *= 0.5*(endPhi-startPhi);
00841     
00842       if(startPhi==0.&& endPhi == twopi)
00843       {
00844         Area += std::fabs(original_parameters->Z_values[i+1]
00845                          -original_parameters->Z_values[i])*
00846                          (original_parameters->Rmax[i]
00847                          +original_parameters->Rmax[i+1]
00848                          -original_parameters->Rmin[i]
00849                          -original_parameters->Rmin[i+1]);
00850       }
00851       areas.push_back(Area);
00852       totArea += Area;
00853     }
00854   
00855     areas.push_back(pi*(sqr(original_parameters->Rmax[numPlanes-1])-
00856                         sqr(original_parameters->Rmin[numPlanes-1])));
00857   
00858     totArea += (areas[0]+areas[numPlanes]);
00859     G4double chose = RandFlat::shoot(0.,totArea);
00860 
00861     if( (chose>=0.) && (chose<areas[0]) )
00862     {
00863       return G4ThreeVector(rRand*cosphi, rRand*sinphi,
00864                            original_parameters->Z_values[0]);
00865     }
00866   
00867     for (i=0; i<numPlanes-1; i++)
00868     {
00869       Achose1 += areas[i];
00870       Achose2 = (Achose1+areas[i+1]);
00871       if(chose>=Achose1 && chose<Achose2)
00872       {
00873         return GetPointOnCut(original_parameters->Rmin[i],
00874                              original_parameters->Rmax[i],
00875                              original_parameters->Rmin[i+1],
00876                              original_parameters->Rmax[i+1],
00877                              original_parameters->Z_values[i],
00878                              original_parameters->Z_values[i+1], Area);
00879       }
00880     }
00881 
00882     rRand = original_parameters->Rmin[numPlanes-1] +
00883       ( (original_parameters->Rmax[numPlanes-1]-original_parameters->Rmin[numPlanes-1])
00884         * std::sqrt(RandFlat::shoot()) );
00885 
00886     return G4ThreeVector(rRand*cosphi,rRand*sinphi,
00887                          original_parameters->Z_values[numPlanes-1]);  
00888   }
00889   else  // Generic Polycone
00890   {
00891     return GetPointOnSurfaceGeneric();  
00892   }
00893 }
00894 
00895 //
00896 // CreatePolyhedron
00897 //
00898 G4Polyhedron* G4Polycone::CreatePolyhedron() const
00899 { 
00900   //
00901   // This has to be fixed in visualization. Fake it for the moment.
00902   // 
00903   if (!genericPcon)
00904   {
00905     return new G4PolyhedronPcon( original_parameters->Start_angle,
00906                                  original_parameters->Opening_angle,
00907                                  original_parameters->Num_z_planes,
00908                                  original_parameters->Z_values,
00909                                  original_parameters->Rmin,
00910                                  original_parameters->Rmax );
00911   }
00912   else
00913   {
00914     // The following code prepares for:
00915     // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
00916     //                                  const double xyz[][3],
00917     //                                  const int faces_vec[][4])
00918     // Here is an extract from the header file HepPolyhedron.h:
00936     const G4int numSide =
00937           G4int(G4Polyhedron::GetNumberOfRotationSteps()
00938                 * (endPhi - startPhi) / twopi) + 1;
00939     G4int nNodes;
00940     G4int nFaces;
00941     typedef G4double double3[3];
00942     double3* xyz;
00943     typedef G4int int4[4];
00944     int4* faces_vec;
00945     if (phiIsOpen)
00946     {
00947       // Triangulate open ends. Simple ear-chopping algorithm...
00948       // I'm not sure how robust this algorithm is (J.Allison).
00949       //
00950       std::vector<G4bool> chopped(numCorner, false);
00951       std::vector<G4int*> triQuads;
00952       G4int remaining = numCorner;
00953       G4int iStarter = 0;
00954       while (remaining >= 3)
00955       {
00956         // Find unchopped corners...
00957         //
00958         G4int A = -1, B = -1, C = -1;
00959         G4int iStepper = iStarter;
00960         do
00961         {
00962           if (A < 0)      { A = iStepper; }
00963           else if (B < 0) { B = iStepper; }
00964           else if (C < 0) { C = iStepper; }
00965           do
00966           {
00967             if (++iStepper >= numCorner) { iStepper = 0; }
00968           }
00969           while (chopped[iStepper]);
00970         }
00971         while (C < 0 && iStepper != iStarter);
00972 
00973         // Check triangle at B is pointing outward (an "ear").
00974         // Sign of z cross product determines...
00975         //
00976         G4double BAr = corners[A].r - corners[B].r;
00977         G4double BAz = corners[A].z - corners[B].z;
00978         G4double BCr = corners[C].r - corners[B].r;
00979         G4double BCz = corners[C].z - corners[B].z;
00980         if (BAr * BCz - BAz * BCr < kCarTolerance)
00981         {
00982           G4int* tq = new G4int[3];
00983           tq[0] = A + 1;
00984           tq[1] = B + 1;
00985           tq[2] = C + 1;
00986           triQuads.push_back(tq);
00987           chopped[B] = true;
00988           --remaining;
00989         }
00990         else
00991         {
00992           do
00993           {
00994             if (++iStarter >= numCorner) { iStarter = 0; }
00995           }
00996           while (chopped[iStarter]);
00997         }
00998       }
00999       // Transfer to faces...
01000       //
01001       nNodes = (numSide + 1) * numCorner;
01002       nFaces = numSide * numCorner + 2 * triQuads.size();
01003       faces_vec = new int4[nFaces];
01004       G4int iface = 0;
01005       G4int addition = numCorner * numSide;
01006       G4int d = numCorner - 1;
01007       for (G4int iEnd = 0; iEnd < 2; ++iEnd)
01008       {
01009         for (size_t i = 0; i < triQuads.size(); ++i)
01010         {
01011           // Negative for soft/auxiliary/normally invisible edges...
01012           //
01013           G4int a, b, c;
01014           if (iEnd == 0)
01015           {
01016             a = triQuads[i][0];
01017             b = triQuads[i][1];
01018             c = triQuads[i][2];
01019           }
01020           else
01021           {
01022             a = triQuads[i][0] + addition;
01023             b = triQuads[i][2] + addition;
01024             c = triQuads[i][1] + addition;
01025           }
01026           G4int ab = std::abs(b - a);
01027           G4int bc = std::abs(c - b);
01028           G4int ca = std::abs(a - c);
01029           faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
01030           faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
01031           faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
01032           faces_vec[iface][3] = 0;
01033           ++iface;
01034         }
01035       }
01036 
01037       // Continue with sides...
01038 
01039       xyz = new double3[nNodes];
01040       const G4double dPhi = (endPhi - startPhi) / numSide;
01041       G4double phi = startPhi;
01042       G4int ixyz = 0;
01043       for (G4int iSide = 0; iSide < numSide; ++iSide)
01044       {
01045         for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
01046         {
01047           xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
01048           xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
01049           xyz[ixyz][2] = corners[iCorner].z;
01050           if (iSide == 0)   // startPhi
01051           {
01052             if (iCorner < numCorner - 1)
01053             {
01054               faces_vec[iface][0] = ixyz + 1;
01055               faces_vec[iface][1] = -(ixyz + numCorner + 1);
01056               faces_vec[iface][2] = ixyz + numCorner + 2;
01057               faces_vec[iface][3] = ixyz + 2;
01058             }
01059             else
01060             {
01061               faces_vec[iface][0] = ixyz + 1;
01062               faces_vec[iface][1] = -(ixyz + numCorner + 1);
01063               faces_vec[iface][2] = ixyz + 2;
01064               faces_vec[iface][3] = ixyz - numCorner + 2;
01065             }
01066           }
01067           else if (iSide == numSide - 1)   // endPhi
01068           {
01069             if (iCorner < numCorner - 1)
01070               {
01071                 faces_vec[iface][0] = ixyz + 1;
01072                 faces_vec[iface][1] = ixyz + numCorner + 1;
01073                 faces_vec[iface][2] = ixyz + numCorner + 2;
01074                 faces_vec[iface][3] = -(ixyz + 2);
01075               }
01076             else
01077               {
01078                 faces_vec[iface][0] = ixyz + 1;
01079                 faces_vec[iface][1] = ixyz + numCorner + 1;
01080                 faces_vec[iface][2] = ixyz + 2;
01081                 faces_vec[iface][3] = -(ixyz - numCorner + 2);
01082               }
01083           }
01084           else
01085           {
01086             if (iCorner < numCorner - 1)
01087               {
01088                 faces_vec[iface][0] = ixyz + 1;
01089                 faces_vec[iface][1] = -(ixyz + numCorner + 1);
01090                 faces_vec[iface][2] = ixyz + numCorner + 2;
01091                 faces_vec[iface][3] = -(ixyz + 2);
01092               }
01093               else
01094               {
01095                 faces_vec[iface][0] = ixyz + 1;
01096                 faces_vec[iface][1] = -(ixyz + numCorner + 1);
01097                 faces_vec[iface][2] = ixyz + 2;
01098                 faces_vec[iface][3] = -(ixyz - numCorner + 2);
01099               }
01100             }
01101             ++iface;
01102             ++ixyz;
01103         }
01104         phi += dPhi;
01105       }
01106 
01107       // Last corners...
01108 
01109       for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
01110       {
01111         xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
01112         xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
01113         xyz[ixyz][2] = corners[iCorner].z;
01114         ++ixyz;
01115       }
01116     }
01117     else  // !phiIsOpen - i.e., a complete 360 degrees.
01118     {
01119       nNodes = numSide * numCorner;
01120       nFaces = numSide * numCorner;;
01121       xyz = new double3[nNodes];
01122       faces_vec = new int4[nFaces];
01123       const G4double dPhi = (endPhi - startPhi) / numSide;
01124       G4double phi = startPhi;
01125       G4int ixyz = 0, iface = 0;
01126       for (G4int iSide = 0; iSide < numSide; ++iSide)
01127       {
01128         for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
01129         {
01130           xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
01131           xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
01132           xyz[ixyz][2] = corners[iCorner].z;
01133 
01134           if (iSide < numSide - 1)
01135           {
01136             if (iCorner < numCorner - 1)
01137             {
01138               faces_vec[iface][0] = ixyz + 1;
01139               faces_vec[iface][1] = -(ixyz + numCorner + 1);
01140               faces_vec[iface][2] = ixyz + numCorner + 2;
01141               faces_vec[iface][3] = -(ixyz + 2);
01142             }
01143             else
01144             {
01145               faces_vec[iface][0] = ixyz + 1;
01146               faces_vec[iface][1] = -(ixyz + numCorner + 1);
01147               faces_vec[iface][2] = ixyz + 2;
01148               faces_vec[iface][3] = -(ixyz - numCorner + 2);
01149             }
01150           }
01151           else   // Last side joins ends...
01152           {
01153             if (iCorner < numCorner - 1)
01154             {
01155               faces_vec[iface][0] = ixyz + 1;
01156               faces_vec[iface][1] = -(ixyz + numCorner - nFaces + 1);
01157               faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
01158               faces_vec[iface][3] = -(ixyz + 2);
01159             }
01160             else
01161             {
01162               faces_vec[iface][0] = ixyz + 1;
01163               faces_vec[iface][1] = -(ixyz - nFaces + numCorner + 1);
01164               faces_vec[iface][2] = ixyz - nFaces + 2;
01165               faces_vec[iface][3] = -(ixyz - numCorner + 2);
01166             }
01167           }
01168           ++ixyz;
01169           ++iface;
01170         }
01171         phi += dPhi;
01172       }
01173     }
01174     G4Polyhedron* polyhedron = new G4Polyhedron;
01175     G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
01176     delete [] faces_vec;
01177     delete [] xyz;
01178     if (problem)
01179     {
01180       std::ostringstream message;
01181       message << "Problem creating G4Polyhedron for: " << GetName();
01182       G4Exception("G4Polycone::CreatePolyhedron()", "GeomSolids1002",
01183                   JustWarning, message);
01184       delete polyhedron;
01185       return 0;
01186     }
01187     else
01188     {
01189       return polyhedron;
01190     }
01191   }
01192 }
01193 
01194 //
01195 // CreateNURBS
01196 //
01197 G4NURBS *G4Polycone::CreateNURBS() const
01198 {
01199   return 0;
01200 }
01201 
01202 //
01203 // G4PolyconeHistorical stuff
01204 //
01205 
01206 G4PolyconeHistorical::G4PolyconeHistorical()
01207   : Start_angle(0.), Opening_angle(0.), Num_z_planes(0),
01208     Z_values(0), Rmin(0), Rmax(0)
01209 {
01210 }
01211 
01212 G4PolyconeHistorical::~G4PolyconeHistorical()
01213 {
01214   delete [] Z_values;
01215   delete [] Rmin;
01216   delete [] Rmax;
01217 }
01218 
01219 G4PolyconeHistorical::
01220 G4PolyconeHistorical( const G4PolyconeHistorical &source )
01221 {
01222   Start_angle   = source.Start_angle;
01223   Opening_angle = source.Opening_angle;
01224   Num_z_planes  = source.Num_z_planes;
01225   
01226   Z_values  = new G4double[Num_z_planes];
01227   Rmin      = new G4double[Num_z_planes];
01228   Rmax      = new G4double[Num_z_planes];
01229   
01230   for( G4int i = 0; i < Num_z_planes; i++)
01231   {
01232     Z_values[i] = source.Z_values[i];
01233     Rmin[i]     = source.Rmin[i];
01234     Rmax[i]     = source.Rmax[i];
01235   }
01236 }
01237 
01238 G4PolyconeHistorical&
01239 G4PolyconeHistorical::operator=( const G4PolyconeHistorical& right )
01240 {
01241   if ( &right == this ) return *this;
01242 
01243   if (&right)
01244   {
01245     Start_angle   = right.Start_angle;
01246     Opening_angle = right.Opening_angle;
01247     Num_z_planes  = right.Num_z_planes;
01248   
01249     delete [] Z_values;
01250     delete [] Rmin;
01251     delete [] Rmax;
01252     Z_values  = new G4double[Num_z_planes];
01253     Rmin      = new G4double[Num_z_planes];
01254     Rmax      = new G4double[Num_z_planes];
01255   
01256     for( G4int i = 0; i < Num_z_planes; i++)
01257     {
01258       Z_values[i] = right.Z_values[i];
01259       Rmin[i]     = right.Rmin[i];
01260       Rmax[i]     = right.Rmax[i];
01261     }
01262   }
01263   return *this;
01264 }

Generated on Mon May 27 17:49:23 2013 for Geant4 by  doxygen 1.4.7