G4BREPSolidPCone.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 // $Id$
00027 //
00028 // ----------------------------------------------------------------------
00029 // GEANT 4 class source file
00030 //
00031 // G4BREPSolidPCone.cc
00032 //
00033 // ----------------------------------------------------------------------
00034 // The polyconical solid G4BREPSolidPCone is a shape defined by a set of 
00035 // inner and outer conical or cylindrical surface sections and two planes 
00036 // perpendicular to the Z axis. Each conical surface is defined by its 
00037 // radius at two different planes perpendicular to the Z-axis. Inner and 
00038 // outer conical surfaces are defined using common Z planes. 
00039 // ----------------------------------------------------------------------
00040 
00041 #include "G4Types.hh"
00042 #include <sstream>
00043 
00044 #include "G4BREPSolidPCone.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4FCylindricalSurface.hh"
00048 #include "G4FConicalSurface.hh"
00049 #include "G4CircularCurve.hh"
00050 #include "G4FPlane.hh"
00051 
00052 typedef enum
00053 {
00054   EInverse = 0,
00055   ENormal = 1
00056 } ESurfaceSense;
00057 
00058 G4BREPSolidPCone::G4BREPSolidPCone(const G4String& name,
00059                                          G4double  start_angle,
00060                                          G4double  opening_angle,
00061                                          G4int     num_z_planes, // sections,
00062                                          G4double  z_start,   
00063                                          G4double  z_values[],
00064                                          G4double  RMIN[],
00065                                          G4double  RMAX[] )
00066   : G4BREPSolid(name)
00067 {
00068   // Save contructor parameters
00069   //
00070   constructorParams.start_angle    = start_angle;
00071   constructorParams.opening_angle  = opening_angle;
00072   constructorParams.num_z_planes   = num_z_planes;
00073   constructorParams.z_start        = z_start;
00074   constructorParams.z_values       = 0;
00075   constructorParams.RMIN           = 0;
00076   constructorParams.RMAX           = 0;
00077   
00078   if( num_z_planes > 0 )
00079   {               
00080     constructorParams.z_values       = new G4double[num_z_planes];
00081     constructorParams.RMIN           = new G4double[num_z_planes];
00082     constructorParams.RMAX           = new G4double[num_z_planes];
00083     for( G4int idx = 0; idx < num_z_planes; ++idx )
00084     {
00085       constructorParams.z_values[idx] = z_values[idx];
00086       constructorParams.RMIN[idx]     = RMIN[idx];
00087       constructorParams.RMAX[idx]     = RMAX[idx];      
00088     }
00089   }
00090 
00091   active=1;
00092   InitializePCone();
00093 }
00094 
00095 G4BREPSolidPCone::G4BREPSolidPCone( __void__& a )
00096   : G4BREPSolid(a)
00097 {
00098   constructorParams.start_angle    = 0.;
00099   constructorParams.opening_angle  = 0.;
00100   constructorParams.num_z_planes   = 0;
00101   constructorParams.z_start        = 0.;
00102   constructorParams.z_values = 0;
00103   constructorParams.RMIN = 0;
00104   constructorParams.RMAX = 0;
00105 }
00106 
00107 G4BREPSolidPCone::~G4BREPSolidPCone()
00108 {
00109   if( constructorParams.num_z_planes > 0 )
00110   {
00111     delete [] constructorParams.z_values;
00112     delete [] constructorParams.RMIN;
00113     delete [] constructorParams.RMAX;
00114   }
00115 }
00116 
00117 G4BREPSolidPCone::G4BREPSolidPCone(const G4BREPSolidPCone& rhs)
00118   : G4BREPSolid(rhs)
00119 {
00120   constructorParams.start_angle   = rhs.constructorParams.start_angle;
00121   constructorParams.opening_angle = rhs.constructorParams.opening_angle;
00122   constructorParams.num_z_planes  = rhs.constructorParams.num_z_planes;
00123   constructorParams.z_start       = rhs.constructorParams.z_start;
00124   constructorParams.z_values = 0;
00125   constructorParams.RMIN     = 0;
00126   constructorParams.RMAX     = 0;
00127   G4int nplanes = constructorParams.num_z_planes;
00128   if( nplanes > 0 )
00129   {
00130     constructorParams.z_values = new G4double[nplanes];
00131     constructorParams.RMIN     = new G4double[nplanes];
00132     constructorParams.RMAX     = new G4double[nplanes];
00133     for( G4int idx = 0; idx < nplanes; ++idx )
00134     {
00135       constructorParams.z_values[idx] = rhs.constructorParams.z_values[idx];
00136       constructorParams.RMIN[idx]     = rhs.constructorParams.RMIN[idx];
00137       constructorParams.RMAX[idx]     = rhs.constructorParams.RMAX[idx];      
00138     }
00139   }
00140   
00141   InitializePCone();
00142 }
00143 
00144 G4BREPSolidPCone&
00145 G4BREPSolidPCone::operator = (const G4BREPSolidPCone& rhs) 
00146 {
00147   // Check assignment to self
00148   //
00149   if (this == &rhs)  { return *this; }
00150 
00151   // Copy base class data
00152   //
00153   G4BREPSolid::operator=(rhs);
00154 
00155   // Copy data
00156   //
00157   constructorParams.start_angle   = rhs.constructorParams.start_angle;
00158   constructorParams.opening_angle = rhs.constructorParams.opening_angle;
00159   constructorParams.num_z_planes  = rhs.constructorParams.num_z_planes;
00160   constructorParams.z_start       = rhs.constructorParams.z_start;
00161   G4int nplanes = constructorParams.num_z_planes;
00162   if( nplanes > 0 )
00163   {
00164     delete [] constructorParams.z_values;
00165     delete [] constructorParams.RMIN;
00166     delete [] constructorParams.RMAX;
00167     constructorParams.z_values = new G4double[nplanes];
00168     constructorParams.RMIN     = new G4double[nplanes];
00169     constructorParams.RMAX     = new G4double[nplanes];
00170     for( G4int idx = 0; idx < nplanes; ++idx )
00171     {
00172       constructorParams.z_values[idx] = rhs.constructorParams.z_values[idx];
00173       constructorParams.RMIN[idx]     = rhs.constructorParams.RMIN[idx];
00174       constructorParams.RMAX[idx]     = rhs.constructorParams.RMAX[idx];      
00175     }
00176   }
00177   
00178   InitializePCone();
00179 
00180   return *this;
00181 }  
00182 
00183 void G4BREPSolidPCone::InitializePCone()
00184 {
00185   G4double  opening_angle = constructorParams.opening_angle;
00186   G4int     num_z_planes  = constructorParams.num_z_planes;
00187   G4double  z_start       = constructorParams.z_start;
00188   G4double* z_values      = constructorParams.z_values;
00189   G4double* RMIN          = constructorParams.RMIN;
00190   G4double* RMAX          = constructorParams.RMAX;
00191 
00192   G4int sections= constructorParams.num_z_planes-1;
00193   nb_of_surfaces = 2*sections+2;
00194 
00195   SurfaceVec = new G4Surface*[nb_of_surfaces];
00196   G4ThreeVector Axis(0,0,1);
00197   G4ThreeVector Origin(0,0,z_start);    
00198   G4double Length;
00199   G4ThreeVector LocalOrigin(0,0,z_start);    
00200   G4int a, b = 0;
00201 
00202   G4ThreeVector PlaneAxis(0, 0, 1);
00203   G4ThreeVector PlaneDir (0, 1, 0);   
00204 
00206   // Test delta phi
00207   
00208   // At the moment (11/03/2002) the phi section is not implemented
00209   // so we take a G4 application down if there is a request for such
00210   // a configuration
00211   //
00212   if( opening_angle < 2*pi-perMillion )
00213   {
00214     G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00215        "GeomSolids0001", FatalException,
00216        "Sorry, phi-section not supported yet, try to use G4Polycone instead!");
00217   }
00218 
00220   // Test the validity of the R values
00221   
00222   // RMIN[0] and RMIN[num_z_planes-1] cannot be = 0
00223   // when RMIN[0] or RMIN[num_z_planes-1] are = 0
00224   //
00225   if(  ((RMIN[0] == 0)              && (RMAX[0] == 0))             || 
00226        ((RMIN[num_z_planes-1] == 0) && (RMAX[num_z_planes-1] == 0))   )
00227       G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00228                   "GeomSolids0002", FatalException,
00229                   "RMIN at the extremities cannot be 0 when RMAX=0 !");  
00230   
00231   // only RMAX[0] and RMAX[num_z_planes-1] can be = 0
00232   //
00233   for(a = 1; a < num_z_planes-1; a++)
00234     if (RMAX[a] == 0)
00235       G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00236                   "GeomSolids0002", FatalException,
00237                   "RMAX inside the solid cannot be 0 !");
00238   
00239   // RMAX[a] must be greater than RMIN[a]
00240   //
00241   for(a = 1; a < num_z_planes-1; a++)
00242     if (RMIN[a] >= RMAX[a])
00243       G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00244                   "GeomSolids0002", FatalException,
00245                   "RMAX must be greater than RMIN in the middle Z planes !");
00246   
00247   if( (RMIN[num_z_planes-1] > RMAX[num_z_planes-1] )
00248       || (RMIN[0] > RMAX[0]) ) 
00249       G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00250                   "GeomSolids0002", FatalException,
00251                   "RMAX must be greater or equal than RMIN at the ends !");
00252 
00253   // Create surfaces
00254 
00255   for( a = 0; a < sections; a++)
00256   {
00257     // Surface length
00258     //
00259     Length = z_values[a+1] - z_values[a];
00260 
00261     if( Length == 0 )
00262     {
00263       // We need to create planar surface(s)
00264       //
00265       if( RMAX[a] != RMAX[a+1] && RMIN[a] != RMIN[a+1] )
00266       {
00267         // We can have the 8 following configurations here:
00268         //
00269         // 1.     2.     3.     4.    
00270         // --+      +--  --+      +--   
00271         // xx|->  <-|xx  xx|      |xx   
00272         // xx+--  --+xx  --+      +--   
00273         // xxxxx  xxxxx    |      |     
00274         // xxxxx  xxxxx    +--  --+   
00275         // xx+--  --+xx    |xx  xx|   
00276         // xx|->  <-|xx    +--  --+   
00277         // --+      +--  
00278         // -------------------------- Z axis
00279         //
00282         //
00283         // 5.     6.     7.     8.
00284         // --+      +--  --+      +--
00285         // xx|->  <-|xx  xx|->  <-|xx
00286         // --+--  --+--  xx+--  --+xx
00287         // <-|xx  xx|->  xxxxx  xxxxx
00288         //   +--  --+    --+xx  xx+--
00289         //               <-|xx  xx|->
00290         //                 +--  --+  
00291         // -------------------------- Z axis
00292         //
00293         // NOTE: The pictures show only one half of polycone!
00294         //       The arrows show the expected surface normal direction.
00295         //       The configuration No. 3 and 4 are not valid solids!
00296         
00297         // Eliminate the invalid cases 3 and 4.
00298         // At this point is guaranteed that each RMIN[i] < RMAX[i]
00299         // where i in in interval 0 < i < num_z_planes-1. So:
00300         //
00301         if( RMIN[a] > RMAX[a+1] || RMAX[a] < RMIN[a+1] )
00302         {
00303           std::ostringstream message;
00304           message << "The values of RMIN[" << a
00305                   << "] & RMAX[" << a+1 << "] or RMAX[" << a
00306                   << "] & RMIN[" << a+1 << "]" << G4endl
00307                   << "generate an invalid configuration for solid: "
00308                   << GetName() << "!";
00309           G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00310                       "GeomSolids0002", FatalException, message );
00311         }
00312         
00313         ESurfaceSense UpSurfSense, LowSurfSense;
00314         
00315         // We need to classify all the cases in order to figure out
00316         // the planar surface sense
00317         //
00318         if( RMAX[a] > RMAX[a+1] )
00319         {
00320           // Cases 1, 5, 7
00321           //
00322           if( RMIN[a] < RMIN[a+1] )
00323           {
00324             // Case 1
00325             //
00326             UpSurfSense  = ENormal;
00327             LowSurfSense = ENormal;
00328           }
00329           else if( RMAX[a+1] != RMIN[a])
00330           {
00331             // Case 7
00332             //
00333             UpSurfSense  = ENormal;
00334             LowSurfSense = EInverse;
00335           }
00336           else
00337           {
00338             // Case 5
00339             //
00340             UpSurfSense  = ENormal;
00341             LowSurfSense = EInverse;
00342           }
00343         }
00344         else
00345         {
00346           // Cases 2, 6, 8
00347           //
00348           if( RMIN[a] > RMIN[a+1] )
00349           {
00350             // Case 2
00351             //
00352             UpSurfSense  = EInverse;
00353             LowSurfSense = EInverse;
00354           }
00355           else if( RMIN[a+1] != RMAX[a] )
00356           {
00357             // Case 8
00358             //
00359             UpSurfSense  = EInverse;
00360             LowSurfSense = ENormal;
00361           }
00362           else
00363           {
00364             // Case 6
00365             UpSurfSense  = EInverse;
00366             LowSurfSense = ENormal;
00367           }
00368         }
00369             
00370         SurfaceVec[b] =
00371            ComputePlanarSurface( RMAX[a], RMAX[a+1], LocalOrigin, PlaneAxis,
00372                                  PlaneDir, UpSurfSense );
00373         //SurfaceVec[b]->SetSameSense( UpSurfSense );
00374         b++;
00375         
00376         SurfaceVec[b] =
00377            ComputePlanarSurface( RMIN[a], RMIN[a+1], LocalOrigin, PlaneAxis,
00378                                  PlaneDir, LowSurfSense );
00379         //SurfaceVec[b]->SetSameSense( LowSurfSense );
00380         b++;
00381       }
00382       else
00383       {
00384         // The original code creating single planar surface
00385         // in case where only either RMAX or RMIN have changed
00386         // at the same Z value; e.g.:
00387         //
00388         //     RMAX          RMIN
00389         //    change        change
00390         //
00391         // 1      2      3      4
00392         // --+      +--  -----  -----
00393         // 00|->  <-|00  00000  00000
00394         // 00+--  --+00  --+00  00+--
00395         // 00000  00000  <-|00  00|->
00396         //                 +--  --+
00397         // --------------------------- Z axis
00398         //
00399         // NOTE: The picture shows only one half of polycone!
00400         
00401         G4double R1, R2;
00402         ESurfaceSense SurfSense;
00403         
00404         // test where is the plane surface
00405         // if( RMAX[a] != RMAX[a+1] )
00406         // {
00407         //   R1 = RMAX[a];
00408         //   R2 = RMAX[a+1];
00409         // }
00410         // else if(RMIN[a] != RMIN[a+1])
00411         // {
00412         //   R1 = RMIN[a];
00413         //   R2 = RMIN[a+1];
00414         // }
00415         // else
00416         // {
00417         //   std::ostringstream message;
00418         //   message << "Error in construction." << G4endl
00419         //           << "Exactly the same z, rmin and rmax given for "
00420         //           << "consecutive indices, " << a << " and " << a+1;
00421         //   G4Exception("G4BREPSolidPCone::InitializePCone()",
00422         //               "GeomSolids1002", JustWarning, message);
00423         //   continue; 
00424         // }
00425         if( RMAX[a] != RMAX[a+1] )
00426         {
00427           // Cases 1, 2
00428           //
00429           R1 = RMAX[a];
00430           R2 = RMAX[a+1];
00431           if( R1 > R2 )
00432           {
00433             // Case 1
00434             //
00435             SurfSense = ENormal;
00436           }
00437           else
00438           {
00439             // Case 2
00440             //
00441             SurfSense = EInverse;
00442           }
00443         }
00444         else if(RMIN[a] != RMIN[a+1])
00445         {
00446           // Cases 1, 2
00447           //
00448           R1 = RMIN[a];
00449           R2 = RMIN[a+1];
00450           if( R1 > R2 )
00451           {
00452             // Case 3
00453             //
00454             SurfSense = EInverse;
00455           }
00456           else
00457           {
00458             // Case 4
00459             //
00460             SurfSense = ENormal;
00461           }
00462         }
00463         else
00464         {
00465           std::ostringstream message;
00466           message << "Error in construction." << G4endl
00467                   << "Exactly the same z, rmin and rmax given for"
00468                   << "consecutive indices, " << a << " and " << a+1;
00469           G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00470                       "GeomSolids1001", JustWarning, message);
00471           continue; 
00472         }
00473         
00474         SurfaceVec[b] =
00475            ComputePlanarSurface( R1, R2, LocalOrigin, PlaneAxis,
00476                                  PlaneDir, SurfSense );
00477         //SurfaceVec[b]->SetSameSense( SurfSense );
00478         b++;
00479         
00480         //      SurfaceVec[b]->SetSameSense(1);
00481         nb_of_surfaces--;
00482       }
00483     }
00484     else
00485     {
00486       // The surface to create is conical or cylindrical
00487 
00488       // Inner PCone
00489       //
00490       if(RMIN[a] != RMIN[a+1])
00491       {
00492         // Create cone
00493         //
00494         if(RMIN[a] > RMIN[a+1])
00495         {
00496           G4Vector3D ConeOrigin = G4Vector3D(LocalOrigin) ; 
00497           SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis, Length,
00498                                                 RMIN[a+1], RMIN[a]);
00499           SurfaceVec[b]->SetSameSense(0);
00500         }
00501         else
00502         {      
00503           G4Vector3D Axis2 = G4Vector3D( -1*Axis );
00504           G4Vector3D LocalOrigin2 = G4Vector3D( LocalOrigin + (Length*Axis) );
00505           G4Vector3D ConeOrigin = LocalOrigin2; 
00506           SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis2, Length,
00507                                                 RMIN[a], RMIN[a+1]); 
00508           SurfaceVec[b]->SetSameSense(0);
00509         }
00510         b++;  
00511       }
00512       else
00513       {
00514         if( RMIN[a] == 0 )
00515         {
00516           // Do not create any surface and decrease nb_of_surfaces
00517           //
00518           nb_of_surfaces--;
00519         }
00520         else
00521         {
00522           // Create cylinder
00523           //
00524           G4Vector3D CylOrigin = G4Vector3D( LocalOrigin ); 
00525           SurfaceVec[b] = new G4FCylindricalSurface(CylOrigin, Axis,
00526                                                     RMIN[a], Length );
00527           SurfaceVec[b]->SetSameSense(0);
00528           b++;
00529         }    
00530       }
00531 
00532       // Outer PCone
00533       //
00534       if(RMAX[a] != RMAX[a+1])
00535       {
00536         // Create cone
00537         //
00538         if(RMAX[a] > RMAX[a+1])
00539         {
00540           G4Vector3D ConeOrigin = G4Vector3D( LocalOrigin );
00541           SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis, Length, RMAX[a+1], RMAX[a]);
00542           SurfaceVec[b]->SetSameSense(1);
00543         }
00544         else
00545         {
00546           G4Vector3D Axis2 = G4Vector3D( -1*Axis );
00547           G4Vector3D LocalOrigin2 = G4Vector3D( LocalOrigin + (Length*Axis) );
00548           G4Vector3D ConeOrigin = LocalOrigin2;
00549           SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis2, Length,
00550                                                 RMAX[a], RMAX[a+1]);
00551           SurfaceVec[b]->SetSameSense(1);
00552         }
00553         b++;
00554       }
00555       else
00556       {
00557         // Create cylinder
00558         //
00559         G4Vector3D CylOrigin = G4Vector3D( LocalOrigin );
00560 
00561         if (RMAX[a] == 0)
00562         {
00563           // Do not create any surface and decrease nb_of_surfaces
00564           //
00565           nb_of_surfaces--;
00566         }
00567         else
00568         {
00569           SurfaceVec[b] = new G4FCylindricalSurface(CylOrigin, Axis,
00570                                                     RMAX[a], Length );
00571           SurfaceVec[b]->SetSameSense(1);
00572           b++;
00573         }
00574       }
00575     }
00576 
00577     // Move surface origin to next section
00578     //
00579     LocalOrigin = LocalOrigin + (Length*Axis);
00580   }
00581   
00583   // Create two end planes  
00584 
00585   G4CurveVector cv;
00586   G4CircularCurve* tmp;
00587 
00588   if(RMIN[0] < RMAX[0])
00589   {
00590     // Create start G4Plane & boundaries    
00591     //
00592     G4Point3D ArcStart1a = G4Point3D( Origin + (RMIN[0]*PlaneDir) );
00593     G4Point3D ArcStart1b = G4Point3D( Origin + (RMAX[0]*PlaneDir) );
00594  
00595     if( RMIN[0] > 0.0 )
00596     {
00597       tmp = new G4CircularCurve;
00598       tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, Origin), RMIN[0]);
00599       tmp->SetBounds(ArcStart1a, ArcStart1a);
00600       tmp->SetSameSense(0);
00601       cv.push_back(tmp);
00602     }
00603   
00604     tmp = new G4CircularCurve;
00605     tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, Origin), RMAX[0]);
00606     tmp->SetBounds(ArcStart1b, ArcStart1b);
00607     tmp->SetSameSense(1);
00608     cv.push_back(tmp);
00609 
00610     SurfaceVec[nb_of_surfaces-2] = new G4FPlane(PlaneDir, -PlaneAxis, Origin);
00611     SurfaceVec[nb_of_surfaces-2]->SetBoundaries(&cv);
00612     SurfaceVec[nb_of_surfaces-2]->SetSameSense(0);
00613   }
00614   else
00615   {
00616     // RMIN[0] == RMAX[0] so no surface is needed, it is a line!
00617     //
00618     nb_of_surfaces--;    
00619   }
00620 
00621   if(RMIN[sections] < RMAX[sections])
00622   {
00623     // Create end G4Plane & boundaries
00624     //
00625     G4Point3D ArcStart2a = G4Point3D( LocalOrigin+(RMIN[sections]*PlaneDir) );  
00626     G4Point3D ArcStart2b = G4Point3D( LocalOrigin+(RMAX[sections]*PlaneDir) );    
00627   
00628     cv.clear();
00629 
00630     if( RMIN[sections] > 0.0 )
00631     {
00632       tmp = new G4CircularCurve;
00633       tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, LocalOrigin),
00634                 RMIN[sections]);
00635       tmp->SetBounds(ArcStart2a, ArcStart2a);
00636       tmp->SetSameSense(0);
00637       cv.push_back(tmp);
00638     }
00639     
00640     tmp = new G4CircularCurve;
00641     tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, LocalOrigin), 
00642               RMAX[sections]);
00643     tmp->SetBounds(ArcStart2b, ArcStart2b);
00644     tmp->SetSameSense(1);
00645     cv.push_back(tmp);
00646   
00647     SurfaceVec[nb_of_surfaces-1] = new G4FPlane(PlaneDir, PlaneAxis,
00648                                                 LocalOrigin);
00649     SurfaceVec[nb_of_surfaces-1]->SetBoundaries(&cv);
00650   
00651     // set sense of the surface
00652     //
00653     SurfaceVec[nb_of_surfaces-1]->SetSameSense(0);
00654   }
00655   else
00656   {
00657     // RMIN[0] == RMAX[0] so no surface is needed, it is a line!
00658     //
00659     nb_of_surfaces--;    
00660   }
00661 
00662   Initialize();
00663 }
00664 
00665 void G4BREPSolidPCone::Initialize()
00666 { 
00667   // Computes the bounding box for solids and surfaces
00668   // Converts concave planes to convex
00669   //    
00670   ShortestDistance=1000000;
00671   CheckSurfaceNormals();
00672   
00673   if(!Box || !AxisBox)
00674     IsConvex();
00675   
00676   CalcBBoxes();
00677 }
00678 
00679 EInside G4BREPSolidPCone::Inside(register const G4ThreeVector& Pt) const
00680 {
00681   // This function find if the point Pt is inside, 
00682   // outside or on the surface of the solid
00683 
00684   G4Vector3D v(1, 0, 0.01);
00685   //G4Vector3D v(1, 0, 0); // This will miss the planar surface perp. to Z axis
00686   //G4Vector3D v(0, 0, 1); // This works, however considered as hack not a fix
00687   G4Vector3D Pttmp(Pt);
00688   G4Vector3D Vtmp(v);
00689   G4Ray r(Pttmp, Vtmp);
00690   
00691   // Check if point is inside the PCone bounding box
00692   //
00693   if( !GetBBox()->Inside(Pttmp) )
00694   {
00695     return kOutside;
00696   }
00697 
00698   // Set the surfaces to active again
00699   //
00700   Reset();
00701   
00702   // Test if the bounding box of each surface is intersected
00703   // by the ray. If not, the surface is deactivated.
00704   //
00705   TestSurfaceBBoxes(r);
00706   
00707   G4double dist = kInfinity;
00708   G4bool isIntersected = false;
00709   G4int WhichSurface = 0;
00710 
00711   const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
00712 
00713   // Chech if the point is on the surface, otherwise
00714   // find the nearest intersected suface. If there are not intersections the
00715   // point is outside
00716 
00717   for(G4int a=0; a < nb_of_surfaces; a++)
00718   {
00719     G4Surface* surface = SurfaceVec[a];
00720     
00721     if( surface->IsActive() )
00722     {
00723       G4double hownear = surface->HowNear(Pt);
00724       
00725       if( std::fabs( hownear ) < sqrHalfTolerance )
00726       {
00727         return kSurface;
00728       }
00729 
00730       if( surface->Intersect(r) )
00731       {
00732         isIntersected = true;
00733         hownear = surface->GetDistance();
00734         
00735         if ( std::fabs( hownear ) < dist )
00736         {
00737           dist         = hownear;
00738           WhichSurface = a;
00739         }
00740       }
00741     }
00742   }
00743   
00744   if ( !isIntersected )
00745   {
00746     return kOutside;
00747   }
00748 
00749   // Find the point of intersection on the surface and the normal
00750   // !!!! be carefull the distance is std::sqrt(dist) !!!!
00751 
00752   dist = std::sqrt( dist );
00753   G4Vector3D IntersectionPoint = Pttmp + dist*Vtmp;
00754   G4Vector3D Normal =
00755              SurfaceVec[WhichSurface]->SurfaceNormal( IntersectionPoint );
00756   
00757   G4double dot = Normal*Vtmp;
00758   
00759   return ( (dot > 0) ? kInside : kOutside );
00760 }
00761 
00762 G4ThreeVector
00763 G4BREPSolidPCone::SurfaceNormal(const G4ThreeVector& Pt) const
00764 {  
00765   // This function calculates the normal of the closest surface
00766   // at a point on the surface
00767   // Note : the sense of the normal depends on the sense of the surface 
00768 
00769   G4int        isurf;
00770   G4bool       normflag = false;
00771 
00772   const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
00773  
00774   // Determine if the point is on the surface
00775   //
00776   G4double minDist = kInfinity;
00777   G4int normSurface = 0;
00778   for(isurf = 0; isurf < nb_of_surfaces; isurf++)
00779   {
00780     G4double dist = std::fabs(SurfaceVec[isurf]->HowNear(Pt));
00781     if( minDist > dist )
00782     {
00783       minDist = dist;
00784       normSurface = isurf;
00785     }
00786     if( dist < sqrHalfTolerance)
00787     {
00788       // the point is on this surface
00789       //
00790       normflag = true;
00791       break;
00792     }
00793   }
00794   
00795   // Calculate the normal at this point, if the point is on the
00796   // surface, otherwise compute the normal to the closest surface
00797   //
00798   if ( normflag )  // point on surface
00799   {
00800     G4ThreeVector norm = SurfaceVec[isurf]->SurfaceNormal(Pt);
00801     return norm.unit();
00802   }
00803   else             // point not on surface
00804   {
00805     G4Surface* nSurface = SurfaceVec[normSurface];
00806     G4ThreeVector hitPt = nSurface->GetClosestHit();
00807     G4ThreeVector hitNorm = nSurface->SurfaceNormal(hitPt);
00808     return hitNorm.unit();
00809   }
00810 }
00811 
00812 G4double G4BREPSolidPCone::DistanceToIn(const G4ThreeVector& Pt) const
00813 {
00814   // Calculates the shortest distance ("safety") from a point
00815   // outside the solid to any boundary of this solid.
00816   // Return 0 if the point is already inside.
00817 
00818   G4double *dists = new G4double[nb_of_surfaces];
00819   G4int a;
00820 
00821   // Set the surfaces to active again
00822   //
00823   Reset();
00824   
00825   // Compute the shortest distance of the point to each surfaces
00826   // Be careful : it's a signed value
00827   //
00828   for(a=0; a< nb_of_surfaces; a++)  
00829     dists[a] = SurfaceVec[a]->HowNear(Pt);
00830      
00831   G4double Dist = kInfinity;
00832   
00833   // if dists[] is positive, the point is outside
00834   // so take the shortest of the shortest positive distances
00835   // dists[] can be equal to 0 : point on a surface
00836   // ( Problem with the G4FPlane : there is no inside and no outside...
00837   //   So, to test if the point is inside to return 0, utilize the Inside
00838   //   function. But I don`t know if it is really needed because dToIn is 
00839   //   called only if the point is outside )
00840   //
00841   for(a = 0; a < nb_of_surfaces; a++)
00842     if( std::fabs(Dist) > std::fabs(dists[a]) ) 
00843       //if( dists[a] >= 0)
00844       Dist = dists[a];
00845   
00846   delete[] dists;
00847 
00848   if(Dist == kInfinity)
00849     // the point is inside the solid or on a surface
00850     //
00851     return 0;
00852   else 
00853     return std::fabs(Dist);
00854 }
00855 
00856 G4double G4BREPSolidPCone::DistanceToIn(register const G4ThreeVector& Pt, 
00857                                         register const G4ThreeVector& V) const
00858 {
00859   // Calculates the distance from a point outside the solid
00860   // to the solid`s boundary along a specified direction vector.
00861   // 
00862   // Note : Intersections with boundaries less than the 
00863   //        tolerance must be ignored if the direction 
00864   //        is away from the boundary
00865   
00866   G4int a;
00867   
00868   // Set the surfaces to active again
00869   //
00870   Reset();
00871   
00872   G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;    
00873   G4Vector3D Pttmp(Pt);
00874   G4Vector3D Vtmp(V);   
00875   G4Ray r(Pttmp, Vtmp);
00876 
00877   // Test if the bounding box of each surface is intersected
00878   // by the ray. If not, the surface become deactive.
00879   //
00880   TestSurfaceBBoxes(r);
00881   
00882   ShortestDistance = kInfinity;
00883   
00884   for(a=0; a< nb_of_surfaces; a++)
00885   {
00886     if(SurfaceVec[a]->IsActive())
00887     {
00888       // test if the ray intersect the surface
00889       //
00890       G4Vector3D Norm = SurfaceVec[a]->SurfaceNormal(Pttmp);
00891       G4double hownear = SurfaceVec[a]->HowNear(Pt);
00892       
00893       if( (Norm * Vtmp) < 0 && std::fabs(hownear) < sqrHalfTolerance )
00894         return 0;
00895       
00896       if( (SurfaceVec[a]->Intersect(r)) )
00897       {
00898         // if more than 1 surface is intersected,
00899         // take the nearest one
00900         G4double distance = SurfaceVec[a]->GetDistance();
00901         
00902         if( distance < ShortestDistance )
00903         {
00904           if( distance > sqrHalfTolerance )
00905           {
00906             ShortestDistance = distance;
00907           }
00908         }
00909       }
00910     }
00911   }
00912   
00913   // Be careful !
00914   // SurfaceVec->Distance is in fact the squared distance
00915   //
00916   if(ShortestDistance != kInfinity)
00917     return( std::sqrt(ShortestDistance) );
00918   else
00919     // no intersection
00920     //
00921     return kInfinity; 
00922 }
00923 
00924 G4double G4BREPSolidPCone::DistanceToOut(register const G4ThreeVector& Pt, 
00925                                          register const G4ThreeVector& V, 
00926                                                   const G4bool, 
00927                                                         G4bool *validNorm, 
00928                                                         G4ThreeVector *) const
00929 {
00930   // Calculates the distance from a point inside the solid
00931   // to the solid`s boundary along a specified direction vector.
00932   // Return 0 if the point is already outside.
00933   //
00934   // Note : If the shortest distance to a boundary is less 
00935   //        than the tolerance, it is ignored. This allows
00936   //        for a point within a tolerant boundary to leave
00937   //        immediately
00938 
00939   // Set the surfaces to active again
00940   //
00941   Reset();
00942 
00943   const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;    
00944   G4Vector3D Ptv = G4Vector3D( Pt );
00945   G4int a;
00946 
00947   // I don`t understand this line
00948   //
00949   if(validNorm)
00950     *validNorm=false;
00951 
00952   G4Vector3D Pttmp(Pt);
00953   G4Vector3D Vtmp(V);   
00954   
00955   G4Ray r(Pttmp, Vtmp);
00956 
00957   // Test if the bounding box of each surface is intersected
00958   // by the ray. If not, the surface become deactive.
00959   //
00960   TestSurfaceBBoxes(r);
00961   
00962   ShortestDistance = kInfinity;
00963  
00964   for(a=0; a< nb_of_surfaces; a++)
00965   {
00966     if( SurfaceVec[a]->IsActive() )
00967     {
00968       G4Vector3D Norm = SurfaceVec[a]->SurfaceNormal(Pttmp);
00969       G4double hownear = SurfaceVec[a]->HowNear(Pt);
00970       
00971       if( (Norm * Vtmp) > 0 && std::fabs( hownear ) < sqrHalfTolerance )
00972         return 0;
00973       
00974       // test if the ray intersect the surface
00975       //
00976       if( SurfaceVec[a]->Intersect(r) )
00977       {
00978         // if more than 1 surface is intersected,
00979         // take the nearest one
00980         //
00981         G4double distance = SurfaceVec[a]->GetDistance();
00982 
00983         if( distance < ShortestDistance )
00984         {
00985           if( distance > sqrHalfTolerance )
00986           {
00987             ShortestDistance = distance;
00988           }
00989         }
00990       }
00991     }
00992   }
00993 
00994   // Be careful !
00995   // SurfaceVec->Distance is in fact the squared distance
00996   //
00997   if(ShortestDistance != kInfinity)
00998     return std::sqrt(ShortestDistance);
00999   else
01000     // if no intersection is founded, the point is outside
01001     //
01002     return 0; 
01003 }
01004 
01005 G4double G4BREPSolidPCone::DistanceToOut(const G4ThreeVector& Pt) const
01006 {
01007   // Calculates the shortest distance ("safety") from a point
01008   // inside the solid to any boundary of this solid.
01009   // Return 0 if the point is already outside.
01010 
01011   G4double *dists = new G4double[nb_of_surfaces];
01012   G4int a;
01013 
01014   // Set the surfaces to active again
01015   Reset();
01016   
01017   // calcul of the shortest distance of the point to each surfaces
01018   // Be careful : it's a signed value
01019   //
01020   for(a=0; a< nb_of_surfaces; a++)
01021     dists[a] = SurfaceVec[a]->HowNear(Pt);  
01022 
01023   G4double Dist = kInfinity;
01024   
01025   // if dists[] is negative, the point is inside
01026   // so take the shortest of the shortest negative distances
01027   // dists[] can be equal to 0 : point on a surface
01028   // ( Problem with the G4FPlane : there is no inside and no outside...
01029   //   So, to test if the point is outside to return 0, utilize the Inside
01030   //   function. But I don`t know if it is really needed because dToOut is 
01031   //   called only if the point is inside )
01032 
01033   for(a = 0; a < nb_of_surfaces; a++)
01034     if( std::fabs(Dist) > std::fabs(dists[a]) ) 
01035       //if( dists[a] <= 0)
01036       Dist = dists[a];
01037   
01038   delete[] dists;
01039 
01040   if(Dist == kInfinity)
01041     // the point is ouside the solid or on a surface
01042     //
01043     return 0;
01044   else
01045     return std::fabs(Dist);
01046 }
01047 
01048 G4VSolid* G4BREPSolidPCone::Clone() const
01049 {
01050   return new G4BREPSolidPCone(*this);
01051 }
01052 
01053 std::ostream& G4BREPSolidPCone::StreamInfo(std::ostream& os) const
01054 {  
01055   // Streams solid contents to output stream.
01056 
01057   G4BREPSolid::StreamInfo( os )
01058   << "\n start_angle:   " << constructorParams.start_angle
01059   << "\n opening_angle: " << constructorParams.opening_angle
01060   << "\n num_z_planes:  " << constructorParams.num_z_planes
01061   << "\n z_start:       " << constructorParams.z_start
01062   << "\n z_values:      ";
01063   G4int idx;
01064   for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01065   {
01066     os << constructorParams.z_values[idx] << " ";
01067   }
01068   os << "\n RMIN:          "; 
01069   for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01070   {
01071     os << constructorParams.RMIN[idx] << " ";
01072   }
01073   os << "\n RMAX:          ";
01074   for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01075   {
01076     os << constructorParams.RMAX[idx] << " ";
01077   }
01078   os << "\n-----------------------------------------------------------\n";
01079 
01080   return os;
01081 }
01082 
01083 void G4BREPSolidPCone::Reset() const
01084 {
01085   Active(1);
01086   ((G4BREPSolidPCone*)this)->intersectionDistance=kInfinity;
01087   StartInside(0);
01088   for(register int a=0;a<nb_of_surfaces;a++)
01089     SurfaceVec[a]->Reset();
01090   ShortestDistance = kInfinity;
01091 }
01092 
01093 G4Surface*
01094 G4BREPSolidPCone::ComputePlanarSurface( G4double r1,
01095                                         G4double r2,
01096                                         G4ThreeVector& origin,
01097                                         G4ThreeVector& planeAxis,
01098                                         G4ThreeVector& planeDirection,
01099                                         G4int surfSense )
01100 {
01101   // The planar surface to return
01102   G4Surface* planarFace = 0;
01103   
01104   G4CurveVector    cv1;
01105   G4CircularCurve  *tmp1, *tmp2;
01106   
01107   // Create plane surface
01108   G4Point3D ArcStart1 = G4Point3D( origin + (r1 * planeDirection) );
01109   G4Point3D ArcStart2 = G4Point3D( origin + (r2 * planeDirection) );
01110 
01111   if(r1 != 0)  
01112   {
01113     tmp1 = new G4CircularCurve;
01114     tmp1->Init(G4Axis2Placement3D( planeDirection, planeAxis, origin), r1);
01115     tmp1->SetBounds(ArcStart1, ArcStart1);
01116     
01117     if( surfSense )
01118       tmp1->SetSameSense(1);
01119     else
01120       tmp1->SetSameSense(0);
01121 
01122     cv1.push_back(tmp1);
01123   }
01124 
01125   if(r2 != 0)  
01126   {
01127     tmp2 = new G4CircularCurve;
01128     tmp2->Init(G4Axis2Placement3D( planeDirection, planeAxis, origin), r2);
01129     tmp2->SetBounds(ArcStart2, ArcStart2);
01130     
01131     if( surfSense )
01132       tmp2->SetSameSense(0);
01133     else
01134       tmp2->SetSameSense(1);
01135     
01136     cv1.push_back(tmp2);
01137   }
01138 
01139   planarFace = new G4FPlane( planeDirection, planeAxis, origin, surfSense );
01140   
01141   planarFace->SetBoundaries(&cv1);
01142 
01143   return planarFace;
01144 }              
01145 
01146 //  In graphics_reps:   
01147 
01148 #include "G4Polyhedron.hh"   
01149 
01150 G4Polyhedron* G4BREPSolidPCone::CreatePolyhedron() const
01151 {
01152   return new G4PolyhedronPcon( constructorParams.start_angle, 
01153                                constructorParams.opening_angle, 
01154                                constructorParams.num_z_planes, 
01155                                constructorParams.z_values,
01156                                constructorParams.RMIN,
01157                                constructorParams.RMAX );
01158 }

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