G4BREPSolidPolyhedra.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 // G4BREPSolidPolyhedra.cc
00032 //
00033 // ----------------------------------------------------------------------
00034 // The polygonal solid G4BREPSolidPolyhedra is a shape defined by an inner 
00035 // and outer polygonal surface and two planes perpendicular to the Z axis. 
00036 // Each polygonal surface is created by linking a series of polygons created 
00037 // at different planes perpendicular to the Z-axis. All these polygons all 
00038 // have the same number of sides (sides) and are defined at the same Z planes 
00039 // for both inner and outer polygonal surfaces. 
00040 // ----------------------------------------------------------------------
00041 //
00042 // History
00043 // -------
00044 //
00045 // Bug-fix #266 by R.Chytracek:
00046 // - The situation when phi1 = 0 dphi1 = 2*pi and all RMINs = 0.0 is handled
00047 //   now. In this case the inner planes are not created. The fix goes even
00048 //   further this means it considers more than 2 z-planes and inner planes
00049 //   are not created whenever two consecutive RMINs are = 0.0 .
00050 // 
00051 // Corrections by S.Giani:
00052 // - Xaxis now corresponds to phi=0
00053 // - partial angle = phiTotal / Nsides
00054 // - end planes exact boundary calculation for phiTotal < 2pi 
00055 //   (also including case with RMIN=RMAX) 
00056 // - Xaxis now properly rotated to compute correct scope of vertixes
00057 // - corrected surface orientation for outer faces parallel to Z
00058 // - completed explicit setting of the orientation for all faces
00059 // - some comparison between doubles avoided by using tolerances 
00060 // - visualisation parameters made consistent with the use made by
00061 //   constructor of the input arguments (i.e. circumscribed radius).
00062 // ----------------------------------------------------------------------
00063 
00064 #include <sstream>
00065 
00066 #include "G4BREPSolidPolyhedra.hh"
00067 #include "G4PhysicalConstants.hh"
00068 #include "G4SystemOfUnits.hh"
00069 #include "G4FPlane.hh"
00070 
00071 G4BREPSolidPolyhedra::G4BREPSolidPolyhedra(const G4String& name,
00072                                                  G4double  start_angle,
00073                                                  G4double  opening_angle,
00074                                                  G4int     sides,
00075                                                  G4int     num_z_planes,      
00076                                                  G4double  z_start,
00077                                                  G4double  z_values[],
00078                                                  G4double  RMIN[],
00079                                                  G4double  RMAX[] )
00080   : G4BREPSolid(name)
00081 {    
00082   // Store the original parameters, to be used in visualisation
00083   // Note radii are not scaled because this BREP uses the radius of the
00084   // circumscribed circle and also graphics_reps/G4Polyhedron uses the
00085   // radius of the circumscribed circle.
00086   
00087   // Save contructor parameters
00088   //
00089   constructorParams.start_angle    = start_angle;
00090   constructorParams.opening_angle  = opening_angle;
00091   constructorParams.sides          = sides;
00092   constructorParams.num_z_planes   = num_z_planes;
00093   constructorParams.z_start        = z_start;
00094   constructorParams.z_values       = 0;
00095   constructorParams.RMIN           = 0;
00096   constructorParams.RMAX           = 0;
00097   
00098   if( num_z_planes > 0 )
00099   {               
00100     constructorParams.z_values       = new G4double[num_z_planes];
00101     constructorParams.RMIN           = new G4double[num_z_planes];
00102     constructorParams.RMAX           = new G4double[num_z_planes];
00103     for( G4int idx = 0; idx < num_z_planes; ++idx )
00104     {
00105       constructorParams.z_values[idx] = z_values[idx];
00106       constructorParams.RMIN[idx]     = RMIN[idx];
00107       constructorParams.RMAX[idx]     = RMAX[idx];      
00108     }
00109   }
00110 
00111   // z_values[0]  should be equal to z_start, for consistency 
00112   //   with what the constructor does.
00113   // Otherwise the z_values that are shifted by (z_values[0] - z_start) , 
00114   //   because z_values are only used in the form 
00115   //   length = z_values[d+1] - z_values[d];         // JA Apr 2, 97
00116   
00117   if( z_values[0] != z_start )
00118   {
00119     std::ostringstream message;
00120     message << "Construction Error. z_values[0] must be equal to z_start!"
00121             << G4endl
00122             << "        Wrong solid parameters: "
00123             << " z_values[0]= " << z_values[0] << " is not equal to "
00124             << " z_start= " << z_start << ".";
00125     G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00126                  "GeomSolids1002", JustWarning, message );
00127     if( num_z_planes <= 0 )  { constructorParams.z_values = new G4double[1]; }
00128     constructorParams.z_values[0]= z_start; 
00129   }
00130 
00131   active=1;
00132   InitializePolyhedra(); 
00133 }
00134 
00135 G4BREPSolidPolyhedra::G4BREPSolidPolyhedra( __void__& a )
00136   : G4BREPSolid(a)
00137 {
00138   constructorParams.start_angle    = 0.;
00139   constructorParams.opening_angle  = 0.;
00140   constructorParams.sides          = 0;
00141   constructorParams.num_z_planes   = 0;
00142   constructorParams.z_start        = 0.;
00143   constructorParams.z_values = 0;
00144   constructorParams.RMIN = 0;
00145   constructorParams.RMAX = 0;
00146 }
00147 
00148 G4BREPSolidPolyhedra::~G4BREPSolidPolyhedra()
00149 {
00150   if( constructorParams.num_z_planes > 0 )
00151   {
00152     delete [] constructorParams.z_values;
00153     delete [] constructorParams.RMIN;
00154     delete [] constructorParams.RMAX;
00155   }  
00156 }
00157 
00158 G4BREPSolidPolyhedra::G4BREPSolidPolyhedra(const G4BREPSolidPolyhedra& rhs)
00159   : G4BREPSolid(rhs)
00160 {
00161   constructorParams.start_angle    = rhs.constructorParams.start_angle;
00162   constructorParams.opening_angle  = rhs.constructorParams.opening_angle;
00163   constructorParams.sides          = rhs.constructorParams.sides;
00164   constructorParams.num_z_planes   = rhs.constructorParams.num_z_planes;
00165   constructorParams.z_start        = rhs.constructorParams.z_start;
00166   constructorParams.z_values       = 0;
00167   constructorParams.RMIN           = 0;
00168   constructorParams.RMAX           = 0;
00169   G4int num_z_planes = constructorParams.num_z_planes;
00170   if( num_z_planes > 0 )
00171   {               
00172     constructorParams.z_values       = new G4double[num_z_planes];
00173     constructorParams.RMIN           = new G4double[num_z_planes];
00174     constructorParams.RMAX           = new G4double[num_z_planes];
00175     for( G4int idx = 0; idx < num_z_planes; ++idx )
00176     {
00177       constructorParams.z_values[idx] = rhs.constructorParams.z_values[idx];
00178       constructorParams.RMIN[idx]     = rhs.constructorParams.RMIN[idx];
00179       constructorParams.RMAX[idx]     = rhs.constructorParams.RMAX[idx];      
00180     }
00181   }
00182 
00183   InitializePolyhedra();
00184 }
00185 
00186 G4BREPSolidPolyhedra&
00187 G4BREPSolidPolyhedra::operator = (const G4BREPSolidPolyhedra& rhs) 
00188 {
00189   // Check assignment to self
00190   //
00191   if (this == &rhs)  { return *this; }
00192 
00193   // Copy base class data
00194   //
00195   G4BREPSolid::operator=(rhs);
00196 
00197   // Copy data
00198   //
00199   constructorParams.start_angle    = rhs.constructorParams.start_angle;
00200   constructorParams.opening_angle  = rhs.constructorParams.opening_angle;
00201   constructorParams.sides          = rhs.constructorParams.sides;
00202   constructorParams.num_z_planes   = rhs.constructorParams.num_z_planes;
00203   constructorParams.z_start        = rhs.constructorParams.z_start;
00204   G4int num_z_planes = constructorParams.num_z_planes;
00205   if( num_z_planes > 0 )
00206   {               
00207     delete [] constructorParams.z_values;
00208     delete [] constructorParams.RMIN;
00209     delete [] constructorParams.RMAX;
00210     constructorParams.z_values       = new G4double[num_z_planes];
00211     constructorParams.RMIN           = new G4double[num_z_planes];
00212     constructorParams.RMAX           = new G4double[num_z_planes];
00213     for( G4int idx = 0; idx < num_z_planes; ++idx )
00214     {
00215       constructorParams.z_values[idx] = rhs.constructorParams.z_values[idx];
00216       constructorParams.RMIN[idx]     = rhs.constructorParams.RMIN[idx];
00217       constructorParams.RMAX[idx]     = rhs.constructorParams.RMAX[idx];      
00218     }
00219   }
00220   
00221   InitializePolyhedra();
00222 
00223   return *this;
00224 }  
00225 
00226 void G4BREPSolidPolyhedra::InitializePolyhedra()
00227 {
00228   G4double  start_angle   = constructorParams.start_angle;
00229   G4double  opening_angle = constructorParams.opening_angle;
00230   G4int     sides         = constructorParams.sides;
00231   G4int     num_z_planes  = constructorParams.num_z_planes;
00232   G4double  z_start       = constructorParams.z_start;
00233   G4double* z_values      = constructorParams.z_values;
00234   G4double* RMIN          = constructorParams.RMIN;
00235   G4double* RMAX          = constructorParams.RMAX;
00236   G4int sections          = num_z_planes - 1;
00237   
00238   if( opening_angle >= 2*pi-perMillion )
00239   {
00240     nb_of_surfaces = 2*(sections * sides) + 2;
00241   }
00242   else
00243   {
00244     nb_of_surfaces = 2*(sections * sides) + 4;
00245   }
00246 
00247   G4int       MaxNbOfSurfaces = nb_of_surfaces;
00248   G4Surface** MaxSurfaceVec   = new G4Surface*[MaxNbOfSurfaces];
00249   
00250   G4Vector3D Axis(0,0,1);
00251   G4Vector3D XAxis(1,0,0);
00252   G4Vector3D TmpAxis;
00253   G4Point3D  Origin(0,0,z_start);    
00254   G4Point3D  LocalOrigin(0,0,z_start);    
00255   G4double   Length;
00256   G4int      Count     = 0 ;
00257   G4double   PartAngle = (opening_angle)/sides;
00258 
00259 
00261   // Preconditions check
00262   
00263   // Detecting minimal required number of sides
00264   //
00265   if( sides < 3 )
00266   {
00267     G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00268                  "InvalidSetup", FatalException,
00269                  "The solid must have at least 3 sides!" );
00270   }
00271   
00272   // Detecting minimal required number of z-sections
00273   //
00274   if( num_z_planes < 2 )
00275   {
00276     G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00277                  "GeomSolids0002", FatalException,
00278                  "The solid must have at least 2 z-sections!" );
00279   }
00280 
00281   // Detect invalid configurations at the ends of polyhedra which
00282   // would not lead to a valid solid creation and likely to a crash
00283   //
00284   if( z_values[0] == z_values[1]
00285    || z_values[sections-1] == z_values[sections] )
00286   {
00287     G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00288         "GeomSolids0002", FatalException,
00289         "The solid must have the first 2 and the last 2 z-values different!" );
00290   }
00291 
00292   // Find out how the z-values sequence is ordered
00293   //
00294   G4bool increasing;
00295   if( z_values[0] < z_values[1] )
00296   {
00297     increasing = true;
00298   }
00299   else
00300   {
00301     increasing = false;
00302   }
00303   
00304   // Detecting polyhedra teeth.
00305   // It's forbidden to specify unordered, e.g. non-increasing or
00306   // non-decreasing sequence of z-values. It may be provided by a
00307   // specific solid in a future.
00308   //
00309   for( G4int idx = 0; idx < sections; idx++ )
00310   {
00311     if( ( z_values[idx] > z_values[idx+1] &&  increasing ) ||
00312         ( z_values[idx] < z_values[idx+1] && !increasing ) )
00313     {
00314       // ERROR! Invalid sequence of z-values
00315       //
00316       std::ostringstream message;
00317       message << "Unordered, non-increasing or non-decreasing sequence."
00318               << G4endl
00319               << "       Unordered z_values sequence detected !" << G4endl
00320               << "       Check z_values with indexes: "
00321               << idx << " " << (idx+1) << ".";
00322       G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00323                    "GeomSolids0002", FatalException, message );
00324     }
00325   }
00326 
00328 #ifdef G4_EXPERIMENTAL_CODE
00329   // There is one problem when sequence of z values is not increasing in a
00330   // regular way, in other words, it's not purely increasing or decreasing
00331   // Irregular sequence can be provided in order to define a polyhedra having
00332   // teeth as shown on the picture bellow
00333   // In this sequence can happen the following:
00334   //   z[a-1] > z[a] < z[a+1] && z[a+1] >= z[a-1]
00335   // One has to check the RMAX and RMIN values due to the possible
00336   // intersections.
00337   //
00338   // 1     2     3  
00339   // ___   ___   ____
00340   // 00/   00/ _ 000/
00341   // 0/    0/ |0 00|
00342   // V___  V__+0 00+--
00343   // 0000  00000 00000
00344   // ----  ----- -----
00345   // ------------------------------------ z-axis
00346   // 
00347   //
00348   // NOTE: This picture doesn't show all the possible configurations of
00349   //       a polyhedra having teeth when looking at its profile.
00350   //       The picture shows only one half of the polyhedra's profile
00352 
00353   // Experimental code! Not recommended for production, it's incomplete!
00354   // The task is to identify invalid combination of z, RMIN and RMAX values
00355   // in the case of toothydra :-)
00356   //
00357   G4int toothIdx;
00358 
00359   for( G4int idx = 1; idx < sections+1; idx++ )
00360   {
00361     if( z_values[idx-1] > z_values[idx] )
00362     {
00363       G4double toothdist = std::fabs( z_values[idx-1] - z_values[idx] );
00364       G4double aftertoothdist = std::fabs( z_values[idx+1] - z_values[idx] );
00365       if( toothdist > aftertoothdist )
00366       {
00367         // Check for possible intersection
00368         //
00369         if( RMAX[idx-1] < RMAX[idx+1] || RMIN[idx-1] > RMIN[idx+1] )
00370         {
00371           // ERROR! The surface conflict!
00372           //
00373           std::ostringstream message;
00374           message << "Unordered sequence of z_values detected." << G4endl
00375                  << "       Conflicting RMAX or RMIN values!" << G4endl
00376                  << "       Check z_values with indexes: "
00377                  << (idx-1) << " " << idx << " " << (idx+1) << ".";
00378           G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00379                        "GeomSolids0002", FatalException, message );
00380         }
00381       }
00382     }
00383   }
00384 #endif // G4_EXPERIMENTAL_CODE
00386 
00387   for(G4int a=0;a<sections;a++)
00388   {
00389     Length = z_values[a+1] - z_values[a];
00390     
00391     if( Length != 0.0 )
00392     {
00393       TmpAxis= XAxis;
00394       TmpAxis.rotateZ(start_angle);
00395       
00396       // L. Broglia: Be careful in the construction of the planes,
00397       //             see G4FPlane
00398       //
00399       for( G4int b = 0; b < sides; b++ )
00400       {
00401         // Create inner side by calculation of points for the planar surface
00402         // boundary. The order of the points gives the surface sense -> changed
00403         // to explicit sense set-up by R. Chytracek, 12/02/2002
00404         // We must check if a pair of two consecutive RMINs is not = 0.0,
00405         // this means no inner plane exists!
00406         //
00407         if( RMIN[a] != 0.0 )
00408         {
00409           if( RMIN[a+1] != 0.0 )
00410           {
00411             // Standard case
00412             //
00413             MaxSurfaceVec[Count] =
00414                CreateTrapezoidalSurface( RMIN[a], RMIN[a+1],
00415                                          LocalOrigin, Length,
00416                                          TmpAxis, PartAngle, EInverse );
00417           }
00418           else
00419           {
00420             // The special case of r1 > r2 where we end at the
00421             // point (0,0,z[a+1])
00422             //
00423             MaxSurfaceVec[Count] =
00424                CreateTriangularSurface( RMIN[a], RMIN[a+1],
00425                                         LocalOrigin, Length,
00426                                         TmpAxis, PartAngle, EInverse );
00427           }
00428         }
00429         else if( RMIN[a+1] != 0.0 )
00430         {
00431            // The special case of r1 < r2 where we start at the
00432            // point ( 0,0,z[a])
00433            //
00434            MaxSurfaceVec[Count] =
00435               CreateTriangularSurface( RMIN[a], RMIN[a+1], LocalOrigin, Length,
00436                                        TmpAxis, PartAngle, EInverse );
00437         }
00438         else
00439         {
00440           // Insert nothing into the vector of sufaces, we'll replicate
00441           // the vector anyway later
00442           //
00443           MaxSurfaceVec[Count] = 0;
00444 
00445           // We need to reduce the number of planes by 1,
00446           // one we have just skipped
00447           //
00448           nb_of_surfaces--;
00449         }      
00450 
00451         if( MaxSurfaceVec[Count] != 0 )
00452         {
00453           // Rotate axis back for the other surface point calculation
00454           // only in the case any of the Create* methods above have been
00455           // called because they modify the passed in TmpAxis
00456           //
00457           TmpAxis.rotateZ(-PartAngle);
00458         }
00459         
00460         Count++;
00461 
00462         // Create outer side
00463 
00464         if( RMAX[a] != 0.0 )
00465         {
00466           if( RMAX[a+1] != 0.0 )
00467           {
00468             // Standard case
00469             //
00470             MaxSurfaceVec[Count] =
00471                CreateTrapezoidalSurface( RMAX[a], RMAX[a+1],
00472                                          LocalOrigin, Length,
00473                                          TmpAxis, PartAngle, ENormal );
00474           }
00475           else
00476           {
00477             // The special case of r1 > r2 where we end
00478             // at the point (0,0,z[a+1])
00479             //
00480             MaxSurfaceVec[Count] =
00481                CreateTriangularSurface( RMAX[a], RMAX[a+1],
00482                                         LocalOrigin, Length,
00483                                         TmpAxis, PartAngle, ENormal );
00484           }
00485         }
00486         else if( RMAX[a+1] != 0.0 )
00487         {
00488            // The special case of r1 < r2 where we start
00489            // at the point ( 0,0,z[a])
00490            //
00491            MaxSurfaceVec[Count] =
00492               CreateTriangularSurface( RMAX[a], RMAX[a+1], LocalOrigin, Length,
00493                                        TmpAxis, PartAngle, ENormal );
00494         }
00495         else
00496         {
00497            // Two consecutive RMAX values can't be zero as
00498            // it's against the definition of BREP polyhedra
00499            //
00500            G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00501                          "GeomSolids0002", FatalException,
00502                          "Two consecutive RMAX values cannot be zero!" );
00503         }
00504         
00505         Count++;
00506       } // End of for loop over sides
00507     }
00508     else
00509     {
00510       // Create planar surfaces perpendicular to z-axis
00511 
00512       ESurfaceSense OuterSurfSense, InnerSurfSense;
00513       
00514       if( RMAX[a] != RMAX[a+1] && RMIN[a] != RMIN[a+1] )
00515       {
00516         // We're about to create a planar surface perpendicular to z-axis
00517         // We can have the 8 following configurations here:
00518         //
00519         // 1.     2.     3.     4.    
00520         // --+      +--  --+      +--   
00521         // xx|->  <-|xx  xx|      |xx   
00522         // xx+--  --+xx  --+      +--   
00523         // xxxxx  xxxxx    |      |     
00524         // xxxxx  xxxxx    +--  --+   
00525         // xx+--  --+xx    |xx  xx|   
00526         // xx|->  <-|xx    +--  --+   
00527         // --+      +--  
00528         // -------------------------- Z axis
00529         //
00532         //
00533         // 5.     6.     7.     8.
00534         // --+      +--  --+      +--
00535         // xx|->  <-|xx  xx|->  <-|xx
00536         // --+--  --+--  xx+--  --+xx
00537         // <-|xx  xx|->  xxxxx  xxxxx
00538         //   +--  --+    --+xx  xx+--
00539         //               <-|xx  xx|->
00540         //                 +--  --+  
00541         // -------------------------- Z axis
00542         //
00543         // NOTE: The pictures shows only one half of polyhedra!
00544         //       The arrows show the expected surface normal direction.
00545         //       The configuration No. 3 and 4 are not valid solids!
00546 
00547         // Eliminate the invalid cases 3 and 4.
00548         // At this point is guaranteed that each RMIN[i] < RMAX[i]
00549         // where i in in interval 0 < i < num_z_planes-1. So:
00550         //
00551         if( RMIN[a] > RMAX[a+1] || RMAX[a] < RMIN[a+1] )
00552         {
00553           std::ostringstream message;
00554           message << "The values of RMIN[" << a << "] & RMAX[" << a+1
00555                   << "] or RMAX[" << a << "] & RMIN[" << a+1 << "]" << G4endl
00556                   << "generate an invalid configuration of solid: "
00557                   << GetName() << "!";
00558           G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00559                        "GeomSolids0002", FatalException, message );
00560         }
00561 
00562         // We need to clasify all the cases in order to figure out
00563         // the planar surface sense
00564         //
00565         if( RMAX[a] > RMAX[a+1] )
00566         {
00567           // Cases 1, 5, 7
00568           //
00569           if( RMIN[a] < RMIN[a+1] )
00570           {
00571             // Case 1
00572             OuterSurfSense  = EInverse;
00573             InnerSurfSense = EInverse;
00574           }
00575           else if( RMAX[a+1] != RMIN[a])
00576           {
00577             // Case 7
00578             OuterSurfSense  = EInverse;
00579             InnerSurfSense = ENormal;
00580           }
00581           else
00582           {
00583             // Case 5
00584             OuterSurfSense  = EInverse;
00585             InnerSurfSense = ENormal;
00586           }
00587         }
00588         else
00589         {
00590           // Cases 2, 6, 8
00591           if( RMIN[a] > RMIN[a+1] )
00592           {
00593             // Case 2
00594             OuterSurfSense  = ENormal;
00595             InnerSurfSense = ENormal;
00596           }
00597           else if( RMIN[a+1] != RMAX[a] )
00598           {
00599             // Case 8
00600             OuterSurfSense  = ENormal;
00601             InnerSurfSense = EInverse;
00602           }
00603           else
00604           {
00605             // Case 6
00606             OuterSurfSense  = ENormal;
00607             InnerSurfSense = EInverse;
00608           }
00609         }
00610         
00611         TmpAxis= XAxis;
00612         TmpAxis.rotateZ(start_angle);
00613         
00614         // Compute the outer planar surface
00615         //
00616         MaxSurfaceVec[Count] =
00617            ComputePlanarSurface( RMAX[a], RMAX[a+1], LocalOrigin, TmpAxis,
00618                                  sides, PartAngle, OuterSurfSense );
00619         if( MaxSurfaceVec[Count] == 0 )
00620         {
00621           // No surface was created
00622           //
00623           nb_of_surfaces--;
00624         }
00625         Count++;
00626         
00627         TmpAxis= XAxis;
00628         TmpAxis.rotateZ(start_angle);
00629         
00630         // Compute the inner planar surface
00631         //
00632         MaxSurfaceVec[Count] =
00633            ComputePlanarSurface( RMIN[a], RMIN[a+1], LocalOrigin, TmpAxis,
00634                                  sides, PartAngle, InnerSurfSense );
00635         if( MaxSurfaceVec[Count] == 0 )
00636         {
00637           // No surface was created
00638           //
00639           nb_of_surfaces--;
00640         }        
00641         Count++;
00642         
00643         // Since we can create here at maximum 2 surfaces
00644         // we need to reflect this in the total
00645         //
00646         nb_of_surfaces -= (2*(sides-1));
00647       }
00648       else
00649       {
00650         // The case where only one of the radius values has changed
00651         //
00652         //     RMAX          RMIN
00653         //    change        change
00654         //
00655         // 1      2      3      4
00656         // --+      +--  -----  -----
00657         // 00|->  <-|00  00000  00000
00658         // 00+--  --+00  --+00  00+--
00659         // 00000  00000  <-|00  00|->
00660         //                 +--  --+
00661         // --------------------------- Z axis
00662         //
00663         // NOTE: The picture shows only one half of polyhedra!
00664         
00665         G4double      R1, R2;
00666         ESurfaceSense SurfSense;
00667         
00668         // The case by case classification
00669         //
00670         if( RMAX[a] != RMAX[a+1] )
00671         {
00672           // Cases 1, 2
00673           //
00674           R1 = RMAX[a];
00675           R2 = RMAX[a+1];
00676           if( R1 > R2 )
00677           {
00678             // Case 1
00679             //
00680             SurfSense = EInverse;
00681           }
00682           else
00683           {
00684             // Case 2
00685             //
00686             SurfSense = ENormal;
00687           }
00688         }
00689         else if(RMIN[a] != RMIN[a+1])
00690         {
00691           // Cases 3, 4
00692           //
00693           R1 = RMIN[a];
00694           R2 = RMIN[a+1];
00695           if( R1 > R2 )
00696           {
00697             // Case 3
00698             //
00699             SurfSense = ENormal;
00700           }
00701           else
00702           {
00703             // Case 4
00704             //
00705             SurfSense = EInverse;
00706           }
00707         }
00708         else
00709         {
00710           std::ostringstream message;
00711           message << "Error in construction." << G4endl
00712                   << "        Exactly the same z, rmin and rmax given for"
00713                   << G4endl
00714                   << "        consecutive indices, " << a << " and " << a+1;
00715           G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00716                        "GeomSolids1001", JustWarning, message );
00717           continue; 
00718         }
00719         TmpAxis= XAxis;
00720         TmpAxis.rotateZ(start_angle);
00721         
00722         MaxSurfaceVec[Count] =
00723            ComputePlanarSurface( R1, R2, LocalOrigin, TmpAxis,
00724                                  sides, PartAngle, SurfSense );
00725         if( MaxSurfaceVec[Count] == 0 )
00726         {
00727           // No surface was created
00728           //
00729           nb_of_surfaces--;
00730         }        
00731         Count++;
00732         
00733         // Since we can create here at maximum 1 surface
00734         // we need to reflect this in the total
00735         //
00736         nb_of_surfaces -= ((2*sides) - 1);
00737       }      
00738     } // End of if( Length != 0.0 )
00739     
00740     LocalOrigin = LocalOrigin + (Length*Axis);
00741 
00742   } // End of for loop over z sections
00743   
00744   if(opening_angle >= 2*pi-perMillion)
00745   {
00746     // Create the end planes for the configuration where delta phi >= 2*PI
00747     
00748     TmpAxis = XAxis;
00749     TmpAxis.rotateZ(start_angle);
00750     
00751     MaxSurfaceVec[Count] =
00752        ComputePlanarSurface( RMIN[0], RMAX[0], Origin, TmpAxis,
00753                              sides, PartAngle, ENormal );
00754     
00755     if( MaxSurfaceVec[Count] == 0 )
00756     {
00757       // No surface was created
00758       //
00759       nb_of_surfaces--;
00760     }
00761     Count++;
00762 
00763     // Reset plane axis
00764     //
00765     TmpAxis = XAxis;
00766     TmpAxis.rotateZ(start_angle);
00767     
00768     MaxSurfaceVec[Count] =
00769        ComputePlanarSurface( RMIN[sections], RMAX[sections],
00770                              LocalOrigin, TmpAxis,
00771                              sides, PartAngle, EInverse );
00772     
00773     if( MaxSurfaceVec[Count] == 0 )
00774     {
00775       // No surface was created
00776       //
00777       nb_of_surfaces--;
00778     }
00779     Count++;
00780   }
00781   else
00782   {
00783     // If delta phi < 2*PI then create a single boundary
00784     // (case with RMIN=0 included)
00785     
00786     // Create the lateral planars
00787     //
00788     TmpAxis             = XAxis;
00789     G4Vector3D TmpAxis2 = XAxis;
00790     TmpAxis.rotateZ(start_angle);
00791     TmpAxis2.rotateZ(start_angle);
00792     TmpAxis2.rotateZ(start_angle);
00793     
00794     LocalOrigin      = Origin;
00795     G4int points     = sections*2+2;
00796     G4int PointCount = 0;
00797     
00798     G4Point3DVector GapPointList(points);
00799     G4Point3DVector GapPointList2(points);
00800 
00801     for(G4int d=0;d<sections+1;d++)
00802     {
00803       GapPointList[PointCount] = LocalOrigin + (RMAX[d]*TmpAxis);
00804       GapPointList[points-1-PointCount] = LocalOrigin + (RMIN[d]*TmpAxis);    
00805       
00806       GapPointList2[PointCount] = LocalOrigin + (RMAX[d]*TmpAxis2);
00807       GapPointList2[points-1-PointCount] = LocalOrigin + (RMIN[d]*TmpAxis2); 
00808          
00809       PointCount++;
00810 
00811       Length = z_values[d+1] - z_values[d];
00812       LocalOrigin = LocalOrigin+(Length*Axis);
00813     }
00814     
00815     // Add the lateral planars to the surfaces list and set/reverse sense
00816     //
00817     MaxSurfaceVec[Count++] = new G4FPlane( &GapPointList,  0, ENormal );
00818     MaxSurfaceVec[Count++] = new G4FPlane( &GapPointList2, 0, EInverse );
00819     
00820     TmpAxis = XAxis;
00821     TmpAxis.rotateZ(start_angle);
00822     TmpAxis.rotateZ(opening_angle);
00823     
00824     // Create end planes
00825     //
00826     G4Point3DVector EndPointList ((sides+1)*2);
00827     G4Point3DVector EndPointList2((sides+1)*2);      
00828     
00829     for(G4int c=0;c<sides+1;c++)
00830     {
00831       // outer polylines for origin end and opposite side
00832       EndPointList[c]  = Origin + (RMAX[0] * TmpAxis);
00833       EndPointList[(sides+1)*2-1-c]  = Origin + (RMIN[0] * TmpAxis);
00834       EndPointList2[c] = LocalOrigin + (RMAX[sections] * TmpAxis);
00835       EndPointList2[(sides+1)*2-1-c] = LocalOrigin + (RMIN[sections] * TmpAxis);
00836       TmpAxis.rotateZ(-PartAngle);
00837     }
00838     
00839     // Add the end planes to the surfaces list
00840     // Note the surface sense in this case is reversed
00841     // It's because here we have created the end planes in reversed order
00842     // than it's done by ComputePlanarSurface() method
00843     //
00844     if(RMAX[0]-RMIN[0] >= perMillion)
00845     {
00846       MaxSurfaceVec[Count] = new G4FPlane( &EndPointList, 0, EInverse );
00847     }
00848     else
00849     {
00850       MaxSurfaceVec[Count] = 0;
00851       nb_of_surfaces--;
00852     }
00853     
00854     Count++;
00855     
00856     if(RMAX[sections]-RMIN[sections] >= perMillion)
00857     {
00858       MaxSurfaceVec[Count] = new G4FPlane( &EndPointList2, 0, ENormal );
00859     }
00860     else
00861     {
00862       MaxSurfaceVec[Count] = 0;
00863       nb_of_surfaces--;
00864     }    
00865   }
00866 
00867   // Now let's replicate the relevant surfaces into
00868   // G4BREPSolid's vector of surfaces
00869   //
00870   SurfaceVec = new G4Surface*[nb_of_surfaces];
00871   G4int sf = 0; G4int zeroCount = 0;
00872   for( G4int srf = 0; srf < MaxNbOfSurfaces; srf++ )
00873   {
00874     if( MaxSurfaceVec[srf] != 0 )
00875     {
00876       if( sf < nb_of_surfaces )
00877       {
00878         SurfaceVec[sf] = MaxSurfaceVec[srf];
00879       }
00880       sf++;
00881     }
00882     else
00883     {
00884       zeroCount++;
00885     }
00886   }
00887 
00888   if( sf != nb_of_surfaces )
00889   {
00890     std::ostringstream message;
00891     message << "Bad number of surfaces!" << G4endl
00892             << "          sf            : "  << sf
00893             << "          nb_of_surfaces: "  << nb_of_surfaces
00894             << "          Count         : "  << Count;
00895     G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00896                  "GeomSolids0002", FatalException, message);
00897   }
00898 
00899   // Clean up the temporary vector of surfaces
00900   //
00901   delete [] MaxSurfaceVec;
00902 
00903   Initialize();
00904 }
00905 
00906 void G4BREPSolidPolyhedra::Initialize()
00907 {
00908   // Calc bounding box for solids and surfaces
00909   // Convert concave planes to convex
00910   //
00911   ShortestDistance=1000000;
00912   CheckSurfaceNormals();
00913   if(!Box || !AxisBox)
00914     IsConvex();
00915   
00916   CalcBBoxes();
00917 }
00918 
00919 void G4BREPSolidPolyhedra::Reset() const
00920 {
00921   Active(1);
00922   ((G4BREPSolidPolyhedra*)this)->intersectionDistance=kInfinity;
00923   StartInside(0);
00924   for(register G4int a=0;a<nb_of_surfaces;a++)
00925     SurfaceVec[a]->Reset();
00926   ShortestDistance = kInfinity;
00927 }
00928 
00929 EInside G4BREPSolidPolyhedra::Inside(register const G4ThreeVector& Pt) const
00930 {
00931   // This function find if the point Pt is inside, 
00932   // outside or on the surface of the solid
00933 
00934   const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;    
00935 
00936   G4Vector3D v(1, 0, 0.01);
00937   G4Vector3D Pttmp(Pt);
00938   G4Vector3D Vtmp(v);
00939   G4Ray r(Pttmp, Vtmp);
00940   
00941   // Check if point is inside the Polyhedra bounding box
00942   //
00943   if( !GetBBox()->Inside(Pttmp) )
00944   {
00945     return kOutside;
00946   }
00947 
00948   // Set the surfaces to active again
00949   //
00950   Reset();
00951   
00952   // Test if the bounding box of each surface is intersected
00953   // by the ray. If not, the surface is deactivated.
00954   //
00955   TestSurfaceBBoxes(r);
00956   
00957   G4int hits=0, samehit=0;
00958 
00959   for(G4int a=0; a < nb_of_surfaces; a++)
00960   {
00961     G4Surface* surface = SurfaceVec[a];
00962 
00963     if(surface->IsActive())
00964     {
00965       // count the number of intersections.
00966       // if this number is odd, the start of the ray is
00967       // inside the volume bounded by the surfaces, so 
00968       // increment the number of intersection by 1 if the 
00969       // point is not on the surface and if this intersection 
00970       // was not found before
00971       //
00972       if( (surface->Intersect(r)) & 1 )
00973       {
00974         // test if the point is on the surface
00975         //
00976         if(surface->GetDistance() < sqrHalfTolerance)
00977         {
00978           return kSurface;
00979         }
00980         // test if this intersection was found before
00981         //
00982         for(G4int i=0; i<a; i++)
00983         {
00984           if(surface->GetDistance() == SurfaceVec[i]->GetDistance())
00985           {
00986             samehit++;
00987             break;
00988           }
00989         }
00990 
00991         // count the number of surfaces intersected by the ray
00992         //
00993         if(!samehit)
00994         {
00995           hits++;
00996         }
00997       }
00998     }
00999   }
01000    
01001   // if the number of surfaces intersected is odd,
01002   // the point is inside the solid
01003   //
01004   return ( (hits&1) ? kInside : kOutside );
01005 }
01006 
01007 G4ThreeVector
01008 G4BREPSolidPolyhedra::SurfaceNormal(const G4ThreeVector& Pt) const
01009 {
01010   // This function calculates the normal of the closest surface
01011   // to the given point
01012   // Note : the sense of the normal depends on the sense of the surface 
01013 
01014   G4int        iplane;
01015   G4bool       normflag = false;
01016   const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;    
01017 
01018   // Determine if the point is on the surface
01019   //
01020   G4double minDist = kInfinity;
01021   G4int normPlane = 0;
01022   for(iplane = 0; iplane < nb_of_surfaces; iplane++)
01023   {
01024     G4double dist = std::fabs(SurfaceVec[iplane]->HowNear(Pt));
01025     if( minDist > dist )
01026     {
01027       minDist = dist;
01028       normPlane = iplane;
01029     }
01030     if( dist < sqrHalfTolerance)
01031     {
01032       // the point is on this surface, so take this as the
01033       // the surface to consider for computing the normal
01034       //
01035       normflag = true;
01036       break;
01037     }
01038   }
01039 
01040   // Calculate the normal at this point, if the point is on the
01041   // surface, otherwise compute the normal to the closest surface
01042   //
01043   if ( normflag )  // point on surface
01044   {
01045     G4ThreeVector norm = SurfaceVec[iplane]->SurfaceNormal(Pt);
01046     return norm.unit();
01047   }
01048   else             // point not on surface
01049   {
01050     G4FPlane* nPlane = (G4FPlane*)(SurfaceVec[normPlane]);
01051     G4ThreeVector hitPt = nPlane->GetSrfPoint();
01052     G4ThreeVector hitNorm = nPlane->SurfaceNormal(hitPt);
01053     return hitNorm.unit();
01054   }
01055 }
01056 
01057 G4double G4BREPSolidPolyhedra::DistanceToIn(const G4ThreeVector& Pt) const
01058 {
01059   // Calculates the shortest distance ("safety") from a point
01060   // outside the solid to any boundary of this solid.
01061   // Return 0 if the point is already inside.
01062 
01063   G4double *dists = new G4double[nb_of_surfaces];
01064   G4int a;
01065 
01066   // Set the surfaces to active again
01067   //
01068   Reset();
01069    
01070   // compute the shortest distance of the point to each surfaces
01071   // Be careful : it's a signed value
01072   //
01073   for(a=0; a< nb_of_surfaces; a++)  
01074     dists[a] = SurfaceVec[a]->HowNear(Pt);
01075      
01076   G4double Dist = kInfinity;
01077   
01078   // if dists[] is positive, the point is outside
01079   // so take the shortest of the shortest positive distances
01080   // dists[] can be equal to 0 : point on a surface
01081   // ( Problem with the G4FPlane : there is no inside and no outside...
01082   //   So, to test if the point is inside to return 0, utilize the Inside
01083   //   function. But I don`t know if it is really needed because dToIn is 
01084   //   called only if the point is outside )
01085   //
01086   for(a = 0; a < nb_of_surfaces; a++)
01087     if( std::fabs(Dist) > std::fabs(dists[a]) ) 
01088       //if( dists[a] >= 0)
01089       Dist = dists[a];
01090   
01091   delete[] dists;
01092 
01093   if(Dist == kInfinity)
01094   {
01095     // the point is inside the solid or on a surface
01096     //
01097     return 0;
01098   }
01099   else 
01100   {
01101     return std::fabs(Dist);
01102   }
01103 }
01104 
01105 G4double
01106 G4BREPSolidPolyhedra::DistanceToIn(register const G4ThreeVector& Pt, 
01107                                    register const G4ThreeVector& V) const
01108 {
01109   // Calculates the distance from a point outside the solid
01110   // to the solid`s boundary along a specified direction vector.
01111   //
01112   // Note : Intersections with boundaries less than the 
01113   //        tolerance must be ignored if the direction 
01114   //        is away from the boundary
01115   
01116   G4int a;
01117   
01118   // Set the surfaces to active again
01119   //
01120   Reset();
01121   
01122   const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;    
01123   G4Vector3D Pttmp(Pt);
01124   G4Vector3D Vtmp(V);   
01125   G4Ray r(Pttmp, Vtmp);
01126 
01127   // Test if the bounding box of each surface is intersected
01128   // by the ray. If not, the surface become deactive.
01129   //
01130   TestSurfaceBBoxes(r);
01131   
01132   ShortestDistance = kInfinity;
01133   
01134   for(a=0; a< nb_of_surfaces; a++)
01135   {
01136     if( SurfaceVec[a]->IsActive() )
01137     {
01138       // test if the ray intersect the surface
01139       //
01140       if( SurfaceVec[a]->Intersect(r) )
01141       {
01142         G4double surfDistance = SurfaceVec[a]->GetDistance();
01143         
01144         // if more than 1 surface is intersected,
01145         // take the nearest one
01146         //
01147         if( surfDistance < ShortestDistance )
01148         {
01149           if( surfDistance > sqrHalfTolerance )
01150           {
01151             ShortestDistance = surfDistance;
01152           }
01153           else
01154           {
01155             // the point is within the boundary
01156             // ignore it if the direction is away from the boundary
01157             //    
01158             G4Vector3D Norm = SurfaceVec[a]->SurfaceNormal(Pttmp);
01159 
01160             if( (Norm * Vtmp) < 0 )
01161             {
01162               ShortestDistance = surfDistance;
01163 //              ShortestDistance = surfDistance==0
01164 //                                 ? sqrHalfTolerance
01165 //                                 : surfDistance;
01166             }
01167           }
01168         }
01169       }
01170     }
01171   }
01172 
01173   // Be careful !
01174   // SurfaceVec->Distance is in fact the squared distance
01175   //
01176   if(ShortestDistance != kInfinity)
01177   {
01178     return std::sqrt(ShortestDistance);
01179   }
01180   else  // no intersection
01181   {
01182     return kInfinity;
01183   }
01184 }
01185 
01186 G4double
01187 G4BREPSolidPolyhedra::DistanceToOut(register const G4ThreeVector& Pt, 
01188                                     register const G4ThreeVector& V, 
01189                                              const G4bool, 
01190                                                    G4bool *validNorm, 
01191                                                    G4ThreeVector *   ) const
01192 {
01193   // Calculates the distance from a point inside the solid
01194   // to the solid`s boundary along a specified direction vector.
01195   // Return 0 if the point is already outside (even number of
01196   // intersections greater than the tolerance).
01197   //
01198   // Note : If the shortest distance to a boundary is less 
01199   //        than the tolerance, it is ignored. This allows
01200   //        for a point within a tolerant boundary to leave
01201   //        immediately
01202 
01203   G4int parity = 0;
01204 
01205   // Set the surfaces to active again
01206   //
01207   Reset();
01208 
01209   const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;    
01210   G4Vector3D Ptv = Pt;
01211   G4int a;
01212 
01213   // I don`t understand this line
01214   //
01215   if(validNorm)
01216     *validNorm=false;
01217 
01218   G4Vector3D Pttmp(Pt);
01219   G4Vector3D Vtmp(V);   
01220   
01221   G4Ray r(Pttmp, Vtmp);
01222 
01223   // Test if the bounding box of each surface is intersected
01224   // by the ray. If not, the surface become deactive.
01225   //
01226   TestSurfaceBBoxes(r);
01227   
01228   ShortestDistance = kInfinity; // this is actually the square of the distance
01229  
01230   for(a=0; a< nb_of_surfaces; a++)
01231   {
01232     G4double surfDistance = SurfaceVec[a]->GetDistance();
01233     
01234     if(SurfaceVec[a]->IsActive())
01235     {
01236       G4int intersects = SurfaceVec[a]->Intersect(r);
01237 
01238       // test if the ray intersects the surface
01239       //
01240       if( intersects != 0 )
01241       {
01242         parity += 1;
01243 
01244         // if more than 1 surface is intersected, take the nearest one
01245         //
01246         if( surfDistance < ShortestDistance )
01247         {
01248           if( surfDistance > sqrHalfTolerance )
01249           {
01250             ShortestDistance = surfDistance;
01251           }
01252           else
01253           {
01254             // the point is within the boundary: ignore it
01255             //
01256             parity -= 1;
01257           }
01258         }
01259       }      
01260     }
01261   }
01262 
01263    G4double distance = 0.;
01264    
01265   // Be careful !
01266   // SurfaceVec->Distance is in fact the squared distance
01267   //
01268    // This condition was changed in order to give not zero answer
01269    // when particle is passing the border of two Touching Surfaces
01270    // and the distance to this surfaces is not zero.
01271    // parity is for the points on the boundary,
01272    // parity is counting only surfDistance<kCarTolerance/2. 
01273    //
01274    //  if((ShortestDistance != kInfinity) && (parity&1))
01275    //  
01276    //
01277    if((ShortestDistance != kInfinity) || (parity&1))
01278   {
01279     distance = std::sqrt(ShortestDistance);
01280   }
01281 
01282   return distance;
01283 }
01284 
01285 G4double G4BREPSolidPolyhedra::DistanceToOut(const G4ThreeVector& Pt) const
01286 {
01287   // Calculates the shortest distance ("safety") from a point
01288   // inside the solid to any boundary of this solid.
01289   // Return 0 if the point is already outside.
01290 
01291   G4double *dists = new G4double[nb_of_surfaces];
01292   G4int a;
01293 
01294   // Set the surfaces to active again
01295   //
01296   Reset();
01297   
01298   // calculate the shortest distance of the point to each surfaces
01299   // Be careful : it's a signed value
01300   //
01301   for(a=0; a< nb_of_surfaces; a++)
01302   {
01303     dists[a] = SurfaceVec[a]->HowNear(Pt);
01304   }
01305 
01306   G4double Dist = kInfinity;
01307   
01308   // if dists[] is negative, the point is inside
01309   // so take the shortest of the shortest negative distances
01310   // dists[] can be equal to 0 : point on a surface
01311   // ( Problem with the G4FPlane : there is no inside and no outside...
01312   //   So, to test if the point is outside to return 0, utilize the Inside
01313   //   function. But I don`t know if it is really needed because dToOut is 
01314   //   called only if the point is inside )
01315 
01316   for(a = 0; a < nb_of_surfaces; a++)
01317   {
01318     if( std::fabs(Dist) > std::fabs(dists[a]) )
01319     {
01320       //if( dists[a] <= 0)
01321       Dist = dists[a];
01322     }
01323   }
01324   
01325   delete[] dists;
01326 
01327   if(Dist == kInfinity)
01328   {
01329     // the point is ouside the solid or on a surface
01330     //
01331     return 0;
01332   }
01333   else
01334   {
01335     // return Dist;
01336     return std::fabs(Dist);
01337   }
01338 }
01339 
01340 G4VSolid* G4BREPSolidPolyhedra::Clone() const
01341 {
01342   return new G4BREPSolidPolyhedra(*this);
01343 }
01344 
01345 std::ostream& G4BREPSolidPolyhedra::StreamInfo(std::ostream& os) const
01346 {
01347   // Streams solid contents to output stream.
01348 
01349   G4BREPSolid::StreamInfo( os )
01350   << "\n start_angle:   " << constructorParams.start_angle
01351   << "\n opening_angle: " << constructorParams.opening_angle
01352   << "\n sides:         " << constructorParams.sides
01353   << "\n num_z_planes:  " << constructorParams.num_z_planes
01354   << "\n z_start:       " << constructorParams.z_start
01355   << "\n z_values:      ";
01356   G4int idx;
01357   for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01358   {
01359     os << constructorParams.z_values[idx] << " ";
01360   }
01361   os << "\n RMIN:          "; 
01362   for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01363   {
01364     os << constructorParams.RMIN[idx] << " ";
01365   }
01366   os << "\n RMAX:          ";
01367   for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01368   {
01369     os << constructorParams.RMAX[idx] << " ";
01370   }
01371   os << "\n-----------------------------------------------------------\n";
01372 
01373   return os;
01374 }
01375 
01376 G4Surface*
01377 G4BREPSolidPolyhedra::CreateTrapezoidalSurface( G4double r1,
01378                                                 G4double r2,
01379                                           const G4Point3D& origin,
01380                                                 G4double distance,
01381                                                 G4Vector3D& xAxis,
01382                                                 G4double partAngle,
01383                                                 ESurfaceSense sense )
01384 {
01385   // The surface to be returned
01386   //
01387   G4Surface* trapsrf = 0;
01388   G4Point3DVector PointList(4);
01389   G4Vector3D zAxis(0,0,1);
01390   
01391   PointList[0] = origin + ( r1       * xAxis);
01392   PointList[3] = origin + ( distance * zAxis)   + (r2 * xAxis);
01393   
01394   xAxis.rotateZ( partAngle );
01395   
01396   PointList[2] = origin + ( distance * zAxis)   + (r2 * xAxis);
01397   PointList[1] = origin + ( r1       * xAxis);  
01398 
01399   // Return the planar trapezoidal surface
01400   //
01401   trapsrf = new G4FPlane( &PointList, 0, sense );
01402   
01403   return trapsrf;
01404 }
01405 
01406 G4Surface*
01407 G4BREPSolidPolyhedra::CreateTriangularSurface( G4double r1,
01408                                                G4double r2,
01409                                          const G4Point3D& origin,
01410                                                G4double distance,
01411                                                G4Vector3D& xAxis,
01412                                                G4double partAngle,
01413                                                ESurfaceSense sense )
01414 {
01415   // The surface to be returned
01416   //
01417   G4Surface*      trapsrf = 0;
01418   G4Point3DVector PointList(3);
01419   G4Vector3D      zAxis(0,0,1);
01420   
01421   PointList[0] = origin + ( r1       * xAxis);
01422   PointList[2] = origin + ( distance * zAxis)   + (r2 * xAxis);
01423   
01424   xAxis.rotateZ( partAngle );
01425   
01426   if( r1 < r2 )
01427   {
01428     PointList[1] = origin + ( distance * zAxis)   + (r2 * xAxis);
01429   }
01430   else
01431   {
01432     PointList[1] = origin + ( r1       * xAxis);  
01433   }
01434 
01435   // Return the planar trapezoidal surface
01436   //
01437   trapsrf = new G4FPlane( &PointList, 0, sense );
01438   
01439   return trapsrf;
01440 }
01441 
01442 G4Surface*
01443 G4BREPSolidPolyhedra::ComputePlanarSurface( G4double r1,
01444                                             G4double r2,
01445                                       const G4Point3D& origin,
01446                                             G4Vector3D& xAxis,
01447                                             G4int sides,
01448                                             G4double partAngle,
01449                                             ESurfaceSense sense )
01450 {
01451   // This method can be called only when r1 != r2,
01452   // otherwise it returns 0 which means that no surface can be
01453   // created out of the given radius pair.
01454   // This method requires the xAxis to be pre-rotated properly.
01455 
01456   G4Point3DVector OuterPointList( sides );
01457   G4Point3DVector InnerPointList( sides );
01458     
01459   G4double   rIn, rOut;
01460   G4Surface* planarSrf = 0;
01461 
01462   if( r1 < r2 )
01463   {
01464     rIn  = r1;
01465     rOut = r2;
01466   }
01467   else if( r1 > r2 )
01468   {
01469     rIn  = r2;
01470     rOut = r1;
01471   }
01472   else
01473   {
01474     // Invalid precondition, the radius values are r1 == r2,
01475     // which means we can create only polyline but no surface
01476     //
01477     return 0;
01478   }
01479 
01480   for( G4int pidx = 0; pidx < sides; pidx++ )
01481   {
01482     // Outer polyline
01483     //
01484     OuterPointList[pidx] = origin + ( rOut * xAxis);
01485     // Inner polyline
01486     //
01487     InnerPointList[pidx] = origin + ( rIn  * xAxis);
01488     xAxis.rotateZ( partAngle );
01489   }
01490 
01491   if( rIn != 0.0 && rOut != 0.0 )
01492   {
01493     // Standard case
01494     //
01495     planarSrf = new G4FPlane( &OuterPointList, &InnerPointList, sense );
01496   }
01497   else if( rOut != 0.0 )
01498   {
01499     // Special case where inner radius is zero so no polyline
01500     // is actually created
01501     //
01502     planarSrf = new G4FPlane( &OuterPointList, 0, sense );
01503   }
01504   else
01505   {
01506     // No surface being created
01507     // This should not happen as filtered out by precondition check above
01508   }
01509   
01510   return planarSrf;
01511 }  
01512 
01513 //  In graphics_reps:
01514 
01515 #include "G4Polyhedron.hh"   
01516 
01517 G4Polyhedron* G4BREPSolidPolyhedra::CreatePolyhedron() const
01518 {
01519   return new G4PolyhedronPgon( constructorParams.start_angle, 
01520                                constructorParams.opening_angle, 
01521                                constructorParams.sides, 
01522                                constructorParams.num_z_planes, 
01523                                constructorParams.z_values,
01524                                constructorParams.RMIN,
01525                                constructorParams.RMAX);
01526 }

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