Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4Polyhedra.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/specific/src/G4Polyhedra.cc (Version 11.3.0) and /geometry/solids/specific/src/G4Polyhedra.cc (Version 4.1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 // Implementation of G4Polyhedra, a CSG polyhe <<  23 //
 27 // as an inherited class of G4VCSGfaceted.     <<  24 // $Id: G4Polyhedra.cc,v 1.3 2001/07/11 10:00:16 gunter Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-04-01 $
                                                   >>  26 //
                                                   >>  27 // 
                                                   >>  28 // --------------------------------------------------------------------
                                                   >>  29 // GEANT 4 class source file
                                                   >>  30 //
                                                   >>  31 //
                                                   >>  32 // G4Polyhedra.cc
                                                   >>  33 //
                                                   >>  34 // Implementation of a CSG polyhedra, as an inherited class of G4VCSGfaceted.
                                                   >>  35 //
                                                   >>  36 // To be done:
                                                   >>  37 //    * Cracks: there are probably small cracks in the seams between the
                                                   >>  38 //      phi face (G4PolyPhiFace) and sides (G4PolyhedraSide) that are not
                                                   >>  39 //      entirely leakproof. Also, I am not sure all vertices are leak proof.
                                                   >>  40 //    * Many optimizations are possible, but not implemented.
                                                   >>  41 //    * Visualization needs to be updated outside of this routine.
 28 //                                                 42 //
 29 // Utility classes:                                43 // Utility classes:
 30 //    * G4EnclosingCylinder: decided a quick c <<  44 //    * G4EnclosingCylinder: I decided a quick check of geometry would be a
 31 //      good idea (for CPU speed). If the quic     45 //      good idea (for CPU speed). If the quick check fails, the regular
 32 //      full-blown G4VCSGfaceted version is in     46 //      full-blown G4VCSGfaceted version is invoked.
 33 //    * G4ReduciblePolygon: Really meant as a      47 //    * G4ReduciblePolygon: Really meant as a check of input parameters,
 34 //      this utility class also "converts" the     48 //      this utility class also "converts" the GEANT3-like PGON/PCON
 35 //      arguments into the newer ones.             49 //      arguments into the newer ones.
 36 // Both these classes are implemented outside      50 // Both these classes are implemented outside this file because they are
 37 // shared with G4Polycone.                         51 // shared with G4Polycone.
 38 //                                                 52 //
 39 // Author: David C. Williams (davidw@scipp.ucs << 
 40 // -------------------------------------------     53 // --------------------------------------------------------------------
 41                                                    54 
 42 #include "G4Polyhedra.hh"                          55 #include "G4Polyhedra.hh"
 43                                                << 
 44 #if !defined(G4GEOM_USE_UPOLYHEDRA)            << 
 45                                                << 
 46 #include "G4PolyhedraSide.hh"                      56 #include "G4PolyhedraSide.hh"
 47 #include "G4PolyPhiFace.hh"                        57 #include "G4PolyPhiFace.hh"
 48                                                    58 
 49 #include "G4GeomTools.hh"                      <<  59 #include "G4Polyhedron.hh"
 50 #include "G4VoxelLimits.hh"                    << 
 51 #include "G4AffineTransform.hh"                << 
 52 #include "G4BoundingEnvelope.hh"               << 
 53                                                << 
 54 #include "G4QuickRand.hh"                      << 
 55                                                << 
 56 #include "G4EnclosingCylinder.hh"                  60 #include "G4EnclosingCylinder.hh"
 57 #include "G4ReduciblePolygon.hh"                   61 #include "G4ReduciblePolygon.hh"
 58 #include "G4VPVParameterisation.hh"            << 
 59                                                << 
 60 namespace                                      << 
 61 {                                              << 
 62   G4Mutex surface_elementsMutex = G4MUTEX_INIT << 
 63 }                                              << 
 64                                                << 
 65 using namespace CLHEP;                         << 
 66                                                    62 
                                                   >>  63 //
 67 // Constructor (GEANT3 style parameters)           64 // Constructor (GEANT3 style parameters)
 68 //                                                 65 //
 69 // GEANT3 PGON radii are specified in the dist     66 // GEANT3 PGON radii are specified in the distance to the norm of each face.
 70 //                                             <<  67 //  
 71 G4Polyhedra::G4Polyhedra( const G4String& name <<  68 G4Polyhedra::G4Polyhedra( const G4String& name, 
 72                                 G4double phiSt     69                                 G4double phiStart,
 73                                 G4double thePh     70                                 G4double thePhiTotal,
 74                                 G4int theNumSi <<  71                     G4int theNumSide, 
 75                                 G4int numZPlan     72                                 G4int numZPlanes,
 76                           const G4double zPlan     73                           const G4double zPlane[],
 77                           const G4double rInne     74                           const G4double rInner[],
 78                           const G4double rOute <<  75                           const G4double rOuter[]  ) : G4VCSGfaceted( name )
 79   : G4VCSGfaceted( name )                      << 
 80 {                                                  76 {
 81   if (theNumSide <= 0)                         <<  77   if (theNumSide <= 0) G4Exception( "G4Polyhedra:: must have at least one side" );
 82   {                                            << 
 83     std::ostringstream message;                << 
 84     message << "Solid must have at least one s << 
 85             << "        No sides specified !"; << 
 86     G4Exception("G4Polyhedra::G4Polyhedra()",  << 
 87                 FatalErrorInArgument, message) << 
 88   }                                            << 
 89                                                << 
 90   //                                           << 
 91   // Calculate conversion factor from G3 radiu << 
 92   //                                           << 
 93   G4double phiTotal = thePhiTotal;             << 
 94   if ( (phiTotal <=0) || (phiTotal >= twopi*(1 << 
 95     { phiTotal = twopi; }                      << 
 96   G4double convertRad = std::cos(0.5*phiTotal/ << 
 97                                                << 
 98   //                                           << 
 99   // Some historical stuff                     << 
100   //                                           << 
101   original_parameters = new G4PolyhedraHistori << 
102                                                << 
103   original_parameters->numSide = theNumSide;   << 
104   original_parameters->Start_angle = phiStart; << 
105   original_parameters->Opening_angle = phiTota << 
106   original_parameters->Num_z_planes = numZPlan << 
107   original_parameters->Z_values = new G4double << 
108   original_parameters->Rmin = new G4double[num << 
109   original_parameters->Rmax = new G4double[num << 
110                                                << 
111   for (G4int i=0; i<numZPlanes; ++i)           << 
112   {                                            << 
113     if (( i < numZPlanes-1) && ( zPlane[i] ==  << 
114     {                                          << 
115       if( (rInner[i]   > rOuter[i+1])          << 
116         ||(rInner[i+1] > rOuter[i])   )        << 
117       {                                        << 
118         DumpInfo();                            << 
119         std::ostringstream message;            << 
120         message << "Cannot create a Polyhedra  << 
121                 << G4endl                      << 
122                 << "        Segments are not c << 
123                 << "        rMin[" << i << "]  << 
124                 << " -- rMax[" << i+1 << "] =  << 
125                 << "        rMin[" << i+1 << " << 
126                 << " -- rMax[" << i << "] = "  << 
127         G4Exception("G4Polyhedra::G4Polyhedra( << 
128                     FatalErrorInArgument, mess << 
129       }                                        << 
130     }                                          << 
131     original_parameters->Z_values[i] = zPlane[ << 
132     original_parameters->Rmin[i] = rInner[i]/c << 
133     original_parameters->Rmax[i] = rOuter[i]/c << 
134   }                                            << 
135                                                << 
136                                                << 
137   //                                           << 
138   // Build RZ polygon using special PCON/PGON  << 
139   //                                           << 
140   auto rz = new G4ReduciblePolygon( rInner, rO << 
141   rz->ScaleA( 1/convertRad );                  << 
142                                                << 
143   //                                           << 
144   // Do the real work                          << 
145   //                                           << 
146   Create( phiStart, phiTotal, theNumSide, rz ) << 
147                                                    78 
148   delete rz;                                   <<  79   //
                                                   >>  80   // Calculate conversion factor from G3 radius to G4 radius
                                                   >>  81   //
                                                   >>  82   G4double phiTotal = thePhiTotal;
                                                   >>  83   if (phiTotal <=0 || phiTotal >= 2*M_PI*(1-DBL_EPSILON)) phiTotal = 2*M_PI;
                                                   >>  84   G4double convertRad = cos(0.5*phiTotal/theNumSide);
                                                   >>  85 
                                                   >>  86   //
                                                   >>  87   // Some historical stuff
                                                   >>  88   //
                                                   >>  89   original_parameters = new G4PolyhedraHistorical;
                                                   >>  90   
                                                   >>  91   original_parameters->numSide = theNumSide;
                                                   >>  92   original_parameters->Start_angle = phiStart;
                                                   >>  93   original_parameters->Opening_angle = phiTotal;
                                                   >>  94   original_parameters->Num_z_planes = numZPlanes;
                                                   >>  95   original_parameters->Z_values = new G4double[numZPlanes];
                                                   >>  96   original_parameters->Rmin = new G4double[numZPlanes];
                                                   >>  97   original_parameters->Rmax = new G4double[numZPlanes];
                                                   >>  98 
                                                   >>  99         G4int i;
                                                   >> 100   for (i=0; i<numZPlanes; i++) {
                                                   >> 101     original_parameters->Z_values[i] = zPlane[i];
                                                   >> 102     original_parameters->Rmin[i] = rInner[i]/convertRad;
                                                   >> 103     original_parameters->Rmax[i] = rOuter[i]/convertRad;
                                                   >> 104   }
                                                   >> 105   
                                                   >> 106   
                                                   >> 107   //
                                                   >> 108   // Build RZ polygon using special PCON/PGON GEANT3 constructor
                                                   >> 109   //
                                                   >> 110   G4ReduciblePolygon *rz = new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
                                                   >> 111   rz->ScaleA( 1/convertRad );
                                                   >> 112   
                                                   >> 113   //
                                                   >> 114   // Do the real work
                                                   >> 115   //
                                                   >> 116   Create( phiStart, phiTotal, theNumSide, rz );
                                                   >> 117   
                                                   >> 118   delete rz;
149 }                                                 119 }
150                                                   120 
                                                   >> 121 
                                                   >> 122 //
151 // Constructor (generic parameters)               123 // Constructor (generic parameters)
152 //                                                124 //
153 G4Polyhedra::G4Polyhedra( const G4String& name << 125 G4Polyhedra::G4Polyhedra( const G4String& name, 
154                                 G4double phiSt << 126                   G4double phiStart,
155                                 G4double phiTo << 127                               G4double phiTotal,
156                                 G4int    theNu << 128                     G4int    theNumSide,  
157                                 G4int    numRZ << 129                   G4int    numRZ,
158                           const G4double r[],  << 130             const G4double r[],
159                           const G4double z[]   << 131             const G4double z[]   ) : G4VCSGfaceted( name )
160   : G4VCSGfaceted( name ), genericPgon(true)   << 132 {
161 {                                              << 133   original_parameters = 0;
162   if (theNumSide <= 0)                         << 134   
163   {                                            << 135   G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
164     std::ostringstream message;                << 136   
165     message << "Solid must have at least one s << 137   Create( phiStart, phiTotal, theNumSide, rz );
166             << "        No sides specified !"; << 138   
167     G4Exception("G4Polyhedra::G4Polyhedra()",  << 139   delete rz;
168                 FatalErrorInArgument, message) << 
169   }                                            << 
170                                                << 
171   auto rz = new G4ReduciblePolygon( r, z, numR << 
172                                                << 
173   Create( phiStart, phiTotal, theNumSide, rz ) << 
174                                                << 
175   // Set original_parameters struct for consis << 
176   //                                           << 
177   SetOriginalParameters(rz);                   << 
178                                                << 
179   delete rz;                                   << 
180 }                                                 140 }
181                                                   141 
                                                   >> 142 
                                                   >> 143 //
182 // Create                                         144 // Create
183 //                                                145 //
184 // Generic create routine, called by each cons << 146 // Generic create routine, called by each constructor after conversion of arguments
185 // after conversion of arguments               << 
186 //                                                147 //
187 void G4Polyhedra::Create( G4double phiStart,      148 void G4Polyhedra::Create( G4double phiStart,
188                           G4double phiTotal,   << 149                       G4double phiTotal,
189                           G4int    theNumSide, << 150               G4int    theNumSide,  
190                           G4ReduciblePolygon*  << 151             G4ReduciblePolygon *rz  )
191 {                                              << 152 {
192   //                                           << 153   //
193   // Perform checks of rz values               << 154   // Perform checks of rz values
194   //                                           << 155   //
195   if (rz->Amin() < 0.0)                        << 156   if (rz->Amin() < 0.0) 
196   {                                            << 157     G4Exception( "G4Polyhedra: Illegal input parameters: All R values must be >= 0" );
197     std::ostringstream message;                << 158     
198     message << "Illegal input parameters - " < << 159   G4double rzArea = rz->Area();
199             << "        All R values must be > << 160   if (rzArea < -kCarTolerance) rz->ReverseOrder();
200     G4Exception("G4Polyhedra::Create()", "Geom << 161 
201                 FatalErrorInArgument, message) << 162   else if (rzArea < -kCarTolerance)
202   }                                            << 163     G4Exception( "G4Polyhedra: Illegal input parameters: R/Z cross section is zero or near zero" );
203                                                << 164     
204   G4double rzArea = rz->Area();                << 165   if ((!rz->RemoveDuplicateVertices( kCarTolerance )) || 
205   if (rzArea < -kCarTolerance)                 << 166       (!rz->RemoveRedundantVertices( kCarTolerance ))     ) 
206   {                                            << 167         G4Exception( "G4Polyhedra: Illegal input parameters: Too few unique R/Z values" );
207     rz->ReverseOrder();                        << 168 
208   }                                            << 169   if (rz->CrossesItself( 1/kInfinity )) 
209   else if (rzArea < kCarTolerance)             << 170     G4Exception( "G4Polyhedra: Illegal input parameters: R/Z segments cross" );
210   {                                            << 171 
211     std::ostringstream message;                << 172   numCorner = rz->NumVertices();
212     message << "Illegal input parameters - " < << 173 
213             << "        R/Z cross section is z << 174 
214     G4Exception("G4Polyhedra::Create()", "Geom << 175   startPhi = phiStart;
215                 FatalErrorInArgument, message) << 176   while( startPhi < 0 ) startPhi += 2*M_PI;
216   }                                            << 177         //
217                                                << 178         // Phi opening? Account for some possible roundoff, and interpret
218   if ( (!rz->RemoveDuplicateVertices( kCarTole << 179         // nonsense value as representing no phi opening
219     || (!rz->RemoveRedundantVertices( kCarTole << 180         //
220   {                                            << 181         if (phiTotal <= 0 || phiTotal > 2.0*M_PI*(1-DBL_EPSILON)) {
221     std::ostringstream message;                << 182                 phiIsOpen = false;
222     message << "Illegal input parameters - " < << 183                 endPhi = phiStart+2*M_PI;
223             << "        Too few unique R/Z val << 184         }
224     G4Exception("G4Polyhedra::Create()", "Geom << 185         else {
225                 FatalErrorInArgument, message) << 186                 phiIsOpen = true;
226   }                                            << 187     
227                                                << 188     //
228   if (rz->CrossesItself( 1/kInfinity ))        << 189     // Convert phi into our convention
229   {                                            << 190     //
230     std::ostringstream message;                << 191     endPhi = phiStart+phiTotal;
231     message << "Illegal input parameters - " < << 192     while( endPhi < startPhi ) endPhi += 2*M_PI;
232             << "        R/Z segments cross !"; << 193         }
233     G4Exception("G4Polyhedra::Create()", "Geom << 194   
234                 FatalErrorInArgument, message) << 195   //
235   }                                            << 196   // Save number sides
236                                                << 197   //
237   numCorner = rz->NumVertices();               << 198   numSide = theNumSide;
238                                                << 199   
239                                                << 200   //
240   startPhi = phiStart;                         << 201   // Allocate corner array. 
241   while( startPhi < 0 )    // Loop checking, 1 << 202   //
242     startPhi += twopi;                         << 203   corners = new G4PolyhedraSideRZ[numCorner];
243   //                                           << 204 
244   // Phi opening? Account for some possible ro << 205   //
245   // nonsense value as representing no phi ope << 206   // Copy corners
246   //                                           << 207   //
247   if ( (phiTotal <= 0) || (phiTotal > twopi*(1 << 208   G4ReduciblePolygonIterator iterRZ(rz);
248   {                                            << 209   
249     phiIsOpen = false;                         << 210   G4PolyhedraSideRZ *next = corners;
250     endPhi = startPhi + twopi;                 << 211   iterRZ.Begin();
251   }                                            << 212   do {
252   else                                         << 213     next->r = iterRZ.GetA();
253   {                                            << 214     next->z = iterRZ.GetB();
254     phiIsOpen = true;                          << 215   } while( ++next, iterRZ.Next() );
255     endPhi = startPhi + phiTotal;              << 216   
256   }                                            << 217   //
257                                                << 218   // Allocate face pointer array
258   //                                           << 219   //
259   // Save number sides                         << 220   numFace = phiIsOpen ? numCorner+2 : numCorner;
260   //                                           << 221   faces = new G4VCSGface*[numFace];
261   numSide = theNumSide;                        << 222   
262                                                << 223   //
263   //                                           << 224   // Construct side faces
264   // Allocate corner array.                    << 225   //
265   //                                           << 226   // To do so properly, we need to keep track of four successive RZ
266   corners = new G4PolyhedraSideRZ[numCorner];  << 227   // corners.
267                                                << 228   //
268   //                                           << 229   // But! Don't construct a face if both points are at zero radius!
269   // Copy corners                              << 230   //
270   //                                           << 231   G4PolyhedraSideRZ *corner = corners,
271   G4ReduciblePolygonIterator iterRZ(rz);       << 232         *prev = corners + numCorner-1,
272                                                << 233         *nextNext;
273   G4PolyhedraSideRZ *next = corners;           << 234   G4VCSGface   **face = faces;
274   iterRZ.Begin();                              << 235   do {
275   do    // Loop checking, 13.08.2015, G.Cosmo  << 236     next = corner+1;
276   {                                            << 237     if (next >= corners+numCorner) next = corners;
277     next->r = iterRZ.GetA();                   << 238     nextNext = next+1;
278     next->z = iterRZ.GetB();                   << 239     if (nextNext >= corners+numCorner) nextNext = corners;
279   } while( ++next, iterRZ.Next() );            << 240     
280                                                << 241     if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
281   //                                           << 242 
282   // Allocate face pointer array               << 243     //
283   //                                           << 244     // We must decide here if we can dare declare one of our faces
284   numFace = phiIsOpen ? numCorner+2 : numCorne << 245     // as having a "valid" normal (i.e. allBehind = true). This
285   faces = new G4VCSGface*[numFace];            << 246     // is never possible if the face faces "inward" in r *unless*
286                                                << 247     // we have only one side
287   //                                           << 248     //
288   // Construct side faces                      << 249     G4bool allBehind;
289   //                                           << 250     if ((corner->z > next->z) && (numSide > 1)) {
290   // To do so properly, we need to keep track  << 251       allBehind = false;
291   // corners.                                  << 252     }
292   //                                           << 253     else {
293   // But! Don't construct a face if both point << 254       //
294   //                                           << 255       // Otherwise, it is only true if the line passing
295   G4PolyhedraSideRZ* corner = corners,         << 256       // through the two points of the segment do not
296                    * prev = corners + numCorne << 257       // split the r/z cross section
297                    * nextNext;                 << 258       //
298   G4VCSGface** face = faces;                   << 259       allBehind = !rz->BisectedBy( corner->r, corner->z,
299   do    // Loop checking, 13.08.2015, G.Cosmo  << 260                  next->r, next->z, kCarTolerance );
300   {                                            << 261     }
301     next = corner+1;                           << 262     
302     if (next >= corners+numCorner) next = corn << 263     *face++ = new G4PolyhedraSide( prev, corner, next, nextNext,
303     nextNext = next+1;                         << 264                  numSide, startPhi, endPhi-startPhi, phiIsOpen );
304     if (nextNext >= corners+numCorner) nextNex << 265   } while( prev=corner, corner=next, corner > corners );
305                                                << 266   
306     if (corner->r < 1/kInfinity && next->r < 1 << 267   if (phiIsOpen) {
307 /*                                             << 268     //
308     // We must decide here if we can dare decl << 269     // Construct phi open edges
309     // as having a "valid" normal (i.e. allBeh << 270     //
310     // is never possible if the face faces "in << 271     *face++ = new G4PolyPhiFace( rz, startPhi, phiTotal/numSide, endPhi );
311     // we have only one side                   << 272     *face++ = new G4PolyPhiFace( rz, endPhi,   phiTotal/numSide, startPhi );
312     //                                         << 273   }
313     G4bool allBehind;                          << 274   
314     if ((corner->z > next->z) && (numSide > 1) << 275   //
315     {                                          << 276   // We might have dropped a face or two: recalculate numFace
316       allBehind = false;                       << 277   //
317     }                                          << 278   numFace = face-faces;
318     else                                       << 279   
319     {                                          << 280   //
320       //                                       << 281   // Make enclosingCylinder
321       // Otherwise, it is only true if the lin << 282   //
322       // through the two points of the segment << 283   enclosingCylinder = new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
323       // split the r/z cross section           << 
324       //                                       << 
325       allBehind = !rz->BisectedBy( corner->r,  << 
326                                    next->r, ne << 
327     }                                          << 
328 */                                             << 
329     *face++ = new G4PolyhedraSide( prev, corne << 
330                  numSide, startPhi, endPhi-sta << 
331   } while( prev=corner, corner=next, corner >  << 
332                                                << 
333   if (phiIsOpen)                               << 
334   {                                            << 
335     //                                         << 
336     // Construct phi open edges                << 
337     //                                         << 
338     *face++ = new G4PolyPhiFace( rz, startPhi, << 
339     *face++ = new G4PolyPhiFace( rz, endPhi,   << 
340   }                                            << 
341                                                << 
342   //                                           << 
343   // We might have dropped a face or two: reca << 
344   //                                           << 
345   numFace = (G4int)(face-faces);               << 
346                                                << 
347   //                                           << 
348   // Make enclosingCylinder                    << 
349   //                                           << 
350   enclosingCylinder =                          << 
351     new G4EnclosingCylinder( rz, phiIsOpen, ph << 
352 }                                                 284 }
353                                                   285 
354 // Fake default constructor - sets only member << 
355 //                            for usage restri << 
356 //                                             << 
357 G4Polyhedra::G4Polyhedra( __void__& a )        << 
358   : G4VCSGfaceted(a), startPhi(0.), endPhi(0.) << 
359 {                                              << 
360 }                                              << 
361                                                   286 
                                                   >> 287 //
362 // Destructor                                     288 // Destructor
363 //                                                289 //
364 G4Polyhedra::~G4Polyhedra()                       290 G4Polyhedra::~G4Polyhedra()
365 {                                                 291 {
366   delete [] corners;                           << 292   delete [] corners;
367   delete original_parameters;                  << 293   if (original_parameters) delete original_parameters;
368   delete enclosingCylinder;                    << 294   
369   delete fElements;                            << 295   delete enclosingCylinder;
370   delete fpPolyhedron;                         << 
371   corners = nullptr;                           << 
372   original_parameters = nullptr;               << 
373   enclosingCylinder = nullptr;                 << 
374   fElements = nullptr;                         << 
375   fpPolyhedron = nullptr;                      << 
376 }                                                 296 }
377                                                   297 
                                                   >> 298 
                                                   >> 299 //
378 // Copy constructor                               300 // Copy constructor
379 //                                                301 //
380 G4Polyhedra::G4Polyhedra( const G4Polyhedra& s << 302 G4Polyhedra::G4Polyhedra( const G4Polyhedra &source ) : G4VCSGfaceted( source )
381   : G4VCSGfaceted( source )                    << 
382 {                                                 303 {
383   CopyStuff( source );                         << 304   CopyStuff( source );
384 }                                                 305 }
385                                                   306 
                                                   >> 307 
                                                   >> 308 //
386 // Assignment operator                            309 // Assignment operator
387 //                                                310 //
388 G4Polyhedra &G4Polyhedra::operator=( const G4P << 311 const G4Polyhedra &G4Polyhedra::operator=( const G4Polyhedra &source )
389 {                                                 312 {
390   if (this == &source) return *this;           << 313   if (this == &source) return *this;
391                                                << 
392   G4VCSGfaceted::operator=( source );          << 
393                                                   314 
394   delete [] corners;                           << 315   G4VCSGfaceted::operator=( source );
395   delete original_parameters;                  << 316   
396   delete enclosingCylinder;                    << 317   delete [] corners;
397                                                << 318   if (original_parameters) delete original_parameters;
398   CopyStuff( source );                         << 319   
399                                                << 320   delete enclosingCylinder;
400   return *this;                                << 321   
                                                   >> 322   CopyStuff( source );
                                                   >> 323   
                                                   >> 324   return *this;
401 }                                                 325 }
402                                                   326 
403 // CopyStuff                                   << 
404 //                                             << 
405 void G4Polyhedra::CopyStuff( const G4Polyhedra << 
406 {                                              << 
407   //                                           << 
408   // Simple stuff                              << 
409   //                                           << 
410   numSide    = source.numSide;                 << 
411   startPhi   = source.startPhi;                << 
412   endPhi     = source.endPhi;                  << 
413   phiIsOpen  = source.phiIsOpen;               << 
414   numCorner  = source.numCorner;               << 
415   genericPgon= source.genericPgon;             << 
416                                                << 
417   //                                           << 
418   // The corner array                          << 
419   //                                           << 
420   corners = new G4PolyhedraSideRZ[numCorner];  << 
421                                                << 
422   G4PolyhedraSideRZ* corn = corners,           << 
423                    * sourceCorn = source.corne << 
424   do    // Loop checking, 13.08.2015, G.Cosmo  << 
425   {                                            << 
426     *corn = *sourceCorn;                       << 
427   } while( ++sourceCorn, ++corn < corners+numC << 
428                                                << 
429   //                                           << 
430   // Original parameters                       << 
431   //                                           << 
432   if (source.original_parameters != nullptr)   << 
433   {                                            << 
434     original_parameters =                      << 
435       new G4PolyhedraHistorical( *source.origi << 
436   }                                            << 
437                                                << 
438   //                                           << 
439   // Enclosing cylinder                        << 
440   //                                           << 
441   enclosingCylinder = new G4EnclosingCylinder( << 
442                                                << 
443   //                                           << 
444   // Surface elements                          << 
445   //                                           << 
446   delete fElements;                            << 
447   fElements = nullptr;                         << 
448                                                << 
449   //                                           << 
450   // Polyhedron                                << 
451   //                                           << 
452   fRebuildPolyhedron = false;                  << 
453   delete fpPolyhedron;                         << 
454   fpPolyhedron = nullptr;                      << 
455 }                                              << 
456                                                   327 
457 // Reset                                       << 
458 //                                                328 //
459 // Recalculates and reshapes the solid, given  << 329 // CopyStuff
460 // original_parameters.                        << 
461 //                                                330 //
462 G4bool G4Polyhedra::Reset()                    << 331 void G4Polyhedra::CopyStuff( const G4Polyhedra &source )
463 {                                                 332 {
464   if (genericPgon)                             << 333   //
465   {                                            << 334   // Simple stuff
466     std::ostringstream message;                << 335   //
467     message << "Solid " << GetName() << " buil << 336   numSide   = source.numSide;
468             << G4endl << "Not applicable to th << 337   startPhi  = source.startPhi;
469     G4Exception("G4Polyhedra::Reset()", "GeomS << 338   endPhi    = source.endPhi;
470                 JustWarning, message, "Paramet << 339   phiIsOpen = source.phiIsOpen;
471     return true;                               << 340   numCorner = source.numCorner;
472   }                                            << 341 
473                                                << 342   //
474   //                                           << 343   // The corner array
475   // Clear old setup                           << 344   //
476   //                                           << 345   corners = new G4PolyhedraSideRZ[numCorner];
477   G4VCSGfaceted::DeleteStuff();                << 346   
478   delete [] corners;                           << 347   G4PolyhedraSideRZ *corn = corners,
479   delete enclosingCylinder;                    << 348         *sourceCorn = source.corners;
480   delete fElements;                            << 349   do {
481   corners = nullptr;                           << 350     *corn = *sourceCorn;
482   fElements = nullptr;                         << 351   } while( ++sourceCorn, ++corn < corners+numCorner );
483   enclosingCylinder = nullptr;                 << 352   
484                                                << 353   //
485   //                                           << 354   // Original parameters
486   // Rebuild polyhedra                         << 355   //
487   //                                           << 356   if (source.original_parameters) {
488   auto rz = new G4ReduciblePolygon( original_p << 357     original_parameters = new G4PolyhedraHistorical( *source.original_parameters );
489                                     original_p << 358   }
490                                     original_p << 359   
491                                     original_p << 360   //
492   Create( original_parameters->Start_angle,    << 361   // Enclosing cylinder
493           original_parameters->Opening_angle,  << 362   //
494           original_parameters->numSide, rz );  << 363   enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder );
495   delete rz;                                   << 
496                                                << 
497   return false;                                << 
498 }                                                 364 }
499                                                   365 
                                                   >> 366 
                                                   >> 367 //
500 // Inside                                         368 // Inside
501 //                                                369 //
502 // This is an override of G4VCSGfaceted::Insid << 370 // This is an override of G4VCSGfaceted::Inside, created in order to speed things
503 // to speed things up by first checking with G << 371 // up by first checking with G4EnclosingCylinder.
504 //                                                372 //
505 EInside G4Polyhedra::Inside( const G4ThreeVect << 373 EInside G4Polyhedra::Inside( const G4ThreeVector &p ) const
506 {                                                 374 {
507   //                                           << 375   //
508   // Quick test                                << 376   // Quick test
509   //                                           << 377   //
510   if (enclosingCylinder->MustBeOutside(p)) ret << 378   if (enclosingCylinder->MustBeOutside(p)) return kOutside;
511                                                << 379 
512   //                                           << 380   //
513   // Long answer                               << 381   // Long answer
514   //                                           << 382   //
515   return G4VCSGfaceted::Inside(p);             << 383   return G4VCSGfaceted::Inside(p);
516 }                                                 384 }
517                                                   385 
518 // DistanceToIn                                << 
519 //                                             << 
520 // This is an override of G4VCSGfaceted::Insid << 
521 // to speed things up by first checking with G << 
522 //                                             << 
523 G4double G4Polyhedra::DistanceToIn( const G4Th << 
524                                     const G4Th << 
525 {                                              << 
526   //                                           << 
527   // Quick test                                << 
528   //                                           << 
529   if (enclosingCylinder->ShouldMiss(p,v))      << 
530     return kInfinity;                          << 
531                                                << 
532   //                                           << 
533   // Long answer                               << 
534   //                                           << 
535   return G4VCSGfaceted::DistanceToIn( p, v );  << 
536 }                                              << 
537                                                   386 
                                                   >> 387 //
538 // DistanceToIn                                   388 // DistanceToIn
539 //                                                389 //
540 G4double G4Polyhedra::DistanceToIn( const G4Th << 390 // This is an override of G4VCSGfaceted::Inside, created in order to speed things
541 {                                              << 391 // up by first checking with G4EnclosingCylinder.
542   return G4VCSGfaceted::DistanceToIn(p);       << 
543 }                                              << 
544                                                << 
545 // Get bounding box                            << 
546 //                                                392 //
547 void G4Polyhedra::BoundingLimits(G4ThreeVector << 393 G4double G4Polyhedra::DistanceToIn( const G4ThreeVector &p, const G4ThreeVector &v ) const
548                                  G4ThreeVector << 
549 {                                                 394 {
550   G4double rmin = kInfinity, rmax = -kInfinity << 395   //
551   G4double zmin = kInfinity, zmax = -kInfinity << 396   // Quick test
552   for (G4int i=0; i<GetNumRZCorner(); ++i)     << 397   //
553   {                                            << 398   if (enclosingCylinder->ShouldMiss(p,v)) return kInfinity;
554     G4PolyhedraSideRZ corner = GetCorner(i);   << 399   
555     if (corner.r < rmin) rmin = corner.r;      << 400   //
556     if (corner.r > rmax) rmax = corner.r;      << 401   // Long answer
557     if (corner.z < zmin) zmin = corner.z;      << 402   //
558     if (corner.z > zmax) zmax = corner.z;      << 403   return G4VCSGfaceted::DistanceToIn( p, v );
559   }                                            << 
560                                                << 
561   G4double sphi    = GetStartPhi();            << 
562   G4double ephi    = GetEndPhi();              << 
563   G4double dphi    = IsOpen() ? ephi-sphi : tw << 
564   G4int    ksteps  = GetNumSide();             << 
565   G4double astep   = dphi/ksteps;              << 
566   G4double sinStep = std::sin(astep);          << 
567   G4double cosStep = std::cos(astep);          << 
568                                                << 
569   G4double sinCur = GetSinStartPhi();          << 
570   G4double cosCur = GetCosStartPhi();          << 
571   if (!IsOpen()) rmin = 0.;                    << 
572   G4double xmin = rmin*cosCur, xmax = xmin;    << 
573   G4double ymin = rmin*sinCur, ymax = ymin;    << 
574   for (G4int k=0; k<ksteps+1; ++k)             << 
575   {                                            << 
576     G4double x = rmax*cosCur;                  << 
577     if (x < xmin) xmin = x;                    << 
578     if (x > xmax) xmax = x;                    << 
579     G4double y = rmax*sinCur;                  << 
580     if (y < ymin) ymin = y;                    << 
581     if (y > ymax) ymax = y;                    << 
582     if (rmin > 0)                              << 
583     {                                          << 
584       G4double xx = rmin*cosCur;               << 
585       if (xx < xmin) xmin = xx;                << 
586       if (xx > xmax) xmax = xx;                << 
587       G4double yy = rmin*sinCur;               << 
588       if (yy < ymin) ymin = yy;                << 
589       if (yy > ymax) ymax = yy;                << 
590     }                                          << 
591     G4double sinTmp = sinCur;                  << 
592     sinCur = sinCur*cosStep + cosCur*sinStep;  << 
593     cosCur = cosCur*cosStep - sinTmp*sinStep;  << 
594   }                                            << 
595   pMin.set(xmin,ymin,zmin);                    << 
596   pMax.set(xmax,ymax,zmax);                    << 
597                                                << 
598   // Check correctness of the bounding box     << 
599   //                                           << 
600   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
601   {                                            << 
602     std::ostringstream message;                << 
603     message << "Bad bounding box (min >= max)  << 
604             << GetName() << " !"               << 
605             << "\npMin = " << pMin             << 
606             << "\npMax = " << pMax;            << 
607     G4Exception("G4Polyhedra::BoundingLimits() << 
608                 JustWarning, message);         << 
609     DumpInfo();                                << 
610   }                                            << 
611 }                                                 404 }
612                                                   405 
613 // Calculate extent under transform and specif << 
614 //                                             << 
615 G4bool G4Polyhedra::CalculateExtent(const EAxi << 
616                                     const G4Vo << 
617                                     const G4Af << 
618                                     G4double&  << 
619 {                                              << 
620   G4ThreeVector bmin, bmax;                    << 
621   G4bool exist;                                << 
622                                                << 
623   // Check bounding box (bbox)                 << 
624   //                                           << 
625   BoundingLimits(bmin,bmax);                   << 
626   G4BoundingEnvelope bbox(bmin,bmax);          << 
627 #ifdef G4BBOX_EXTENT                           << 
628   return bbox.CalculateExtent(pAxis,pVoxelLimi << 
629 #endif                                         << 
630   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 
631   {                                            << 
632     return exist = pMin < pMax;                << 
633   }                                            << 
634                                                << 
635   // To find the extent, RZ contour of the pol << 
636   // in triangles. The extent is calculated as << 
637   // all sub-polycones formed by rotation of t << 
638   //                                           << 
639   G4TwoVectorList contourRZ;                   << 
640   G4TwoVectorList triangles;                   << 
641   std::vector<G4int> iout;                     << 
642   G4double eminlim = pVoxelLimit.GetMinExtent( << 
643   G4double emaxlim = pVoxelLimit.GetMaxExtent( << 
644                                                << 
645   // get RZ contour, ensure anticlockwise orde << 
646   for (G4int i=0; i<GetNumRZCorner(); ++i)     << 
647   {                                            << 
648     G4PolyhedraSideRZ corner = GetCorner(i);   << 
649     contourRZ.emplace_back(corner.r,corner.z); << 
650   }                                            << 
651   G4GeomTools::RemoveRedundantVertices(contour << 
652   G4double area = G4GeomTools::PolygonArea(con << 
653   if (area < 0.) std::reverse(contourRZ.begin( << 
654                                                << 
655   // triangulate RZ countour                   << 
656   if (!G4GeomTools::TriangulatePolygon(contour << 
657   {                                            << 
658     std::ostringstream message;                << 
659     message << "Triangulation of RZ contour ha << 
660             << GetName() << " !"               << 
661             << "\nExtent has been calculated u << 
662     G4Exception("G4Polyhedra::CalculateExtent( << 
663                 "GeomMgt1002",JustWarning,mess << 
664     return bbox.CalculateExtent(pAxis,pVoxelLi << 
665   }                                            << 
666                                                << 
667   // set trigonometric values                  << 
668   G4double sphi     = GetStartPhi();           << 
669   G4double ephi     = GetEndPhi();             << 
670   G4double dphi     = IsOpen() ? ephi-sphi : t << 
671   G4int    ksteps   = GetNumSide();            << 
672   G4double astep    = dphi/ksteps;             << 
673   G4double sinStep  = std::sin(astep);         << 
674   G4double cosStep  = std::cos(astep);         << 
675   G4double sinStart = GetSinStartPhi();        << 
676   G4double cosStart = GetCosStartPhi();        << 
677                                                << 
678   // allocate vector lists                     << 
679   std::vector<const G4ThreeVectorList *> polyg << 
680   polygons.resize(ksteps+1);                   << 
681   for (G4int k=0; k<ksteps+1; ++k)             << 
682   {                                            << 
683     polygons[k] = new G4ThreeVectorList(3);    << 
684   }                                            << 
685                                                << 
686   // main loop along triangles                 << 
687   pMin =  kInfinity;                           << 
688   pMax = -kInfinity;                           << 
689   G4int ntria = (G4int)triangles.size()/3;     << 
690   for (G4int i=0; i<ntria; ++i)                << 
691   {                                            << 
692     G4double sinCur = sinStart;                << 
693     G4double cosCur = cosStart;                << 
694     G4int i3 = i*3;                            << 
695     for (G4int k=0; k<ksteps+1; ++k) // rotate << 
696     {                                          << 
697       auto ptr = const_cast<G4ThreeVectorList* << 
698       auto iter = ptr->begin();                << 
699       iter->set(triangles[i3+0].x()*cosCur,    << 
700                 triangles[i3+0].x()*sinCur,    << 
701                 triangles[i3+0].y());          << 
702       iter++;                                  << 
703       iter->set(triangles[i3+1].x()*cosCur,    << 
704                 triangles[i3+1].x()*sinCur,    << 
705                 triangles[i3+1].y());          << 
706       iter++;                                  << 
707       iter->set(triangles[i3+2].x()*cosCur,    << 
708                 triangles[i3+2].x()*sinCur,    << 
709                 triangles[i3+2].y());          << 
710                                                << 
711       G4double sinTmp = sinCur;                << 
712       sinCur = sinCur*cosStep + cosCur*sinStep << 
713       cosCur = cosCur*cosStep - sinTmp*sinStep << 
714     }                                          << 
715                                                << 
716     // set sub-envelope and adjust extent      << 
717     G4double emin,emax;                        << 
718     G4BoundingEnvelope benv(polygons);         << 
719     if (!benv.CalculateExtent(pAxis,pVoxelLimi << 
720     if (emin < pMin) pMin = emin;              << 
721     if (emax > pMax) pMax = emax;              << 
722     if (eminlim > pMin && emaxlim < pMax) brea << 
723   }                                            << 
724   // free memory                               << 
725   for (G4int k=0; k<ksteps+1; ++k) { delete po << 
726   return (pMin < pMax);                        << 
727 }                                              << 
728                                                   406 
                                                   >> 407 //
729 // ComputeDimensions                              408 // ComputeDimensions
730 //                                                409 //
731 void G4Polyhedra::ComputeDimensions(       G4V << 410 void G4Polyhedra::ComputeDimensions( G4VPVParameterisation* p,
732                                      const G4i << 411             const G4int n,
733                                      const G4V << 412             const G4VPhysicalVolume* pRep)
734 {                                                 413 {
735   p->ComputeDimensions(*this,n,pRep);          << 
736 }                                                 414 }
737                                                   415 
738 // GetEntityType                               << 
739 //                                             << 
740 G4GeometryType G4Polyhedra::GetEntityType() co << 
741 {                                              << 
742   return {"G4Polyhedra"};                      << 
743 }                                              << 
744                                                   416 
745 // IsFaceted                                   << 
746 //                                                417 //
747 G4bool G4Polyhedra::IsFaceted() const          << 418 // CreatePolyhedron
748 {                                              << 
749   return true;                                 << 
750 }                                              << 
751                                                << 
752 // Make a clone of the object                  << 
753 //                                                419 //
754 G4VSolid* G4Polyhedra::Clone() const           << 420 G4Polyhedron *G4Polyhedra::CreatePolyhedron() const
755 {                                              << 421 { 
756   return new G4Polyhedra(*this);               << 422   //
757 }                                              << 423   // This has to be fixed in visualization. Fake it for the moment.
                                                   >> 424   // 
                                                   >> 425   if (original_parameters) {
                                                   >> 426   
                                                   >> 427     return new G4PolyhedronPgon( original_parameters->Start_angle,
                                                   >> 428                original_parameters->Opening_angle,
                                                   >> 429                original_parameters->numSide,
                                                   >> 430                original_parameters->Num_z_planes,
                                                   >> 431                original_parameters->Z_values,
                                                   >> 432                original_parameters->Rmin,
                                                   >> 433                original_parameters->Rmax);
                                                   >> 434   }
                                                   >> 435   else {
                                                   >> 436     G4cerr << "G4Polyhedra: visualization of this type of G4Polyhedra is not supported at this time" << G4endl;
                                                   >> 437     return 0;
                                                   >> 438   }
758                                                   439 
759 // Stream object contents to an output stream  << 440 } 
760 //                                             << 
761 std::ostream& G4Polyhedra::StreamInfo( std::os << 
762 {                                              << 
763   G4long oldprc = os.precision(16);            << 
764   os << "------------------------------------- << 
765      << "    *** Dump for solid - " << GetName << 
766      << "    ================================= << 
767      << " Solid type: G4Polyhedra\n"           << 
768      << " Parameters: \n"                      << 
769      << "    starting phi angle : " << startPh << 
770      << "    ending phi angle   : " << endPhi/ << 
771      << "    number of sides    : " << numSide << 
772   G4int i=0;                                   << 
773   if (!genericPgon)                            << 
774   {                                            << 
775     G4int numPlanes = original_parameters->Num << 
776     os << "    number of Z planes: " << numPla << 
777        << "              Z values: \n";        << 
778     for (i=0; i<numPlanes; ++i)                << 
779     {                                          << 
780       os << "              Z plane " << i << " << 
781          << original_parameters->Z_values[i] < << 
782     }                                          << 
783     os << "              Tangent distances to  << 
784     for (i=0; i<numPlanes; ++i)                << 
785     {                                          << 
786       os << "              Z plane " << i << " << 
787          << original_parameters->Rmin[i] << "\ << 
788     }                                          << 
789     os << "              Tangent distances to  << 
790     for (i=0; i<numPlanes; ++i)                << 
791     {                                          << 
792       os << "              Z plane " << i << " << 
793          << original_parameters->Rmax[i] << "\ << 
794     }                                          << 
795   }                                            << 
796   os << "    number of RZ points: " << numCorn << 
797      << "              RZ values (corners): \n << 
798      for (i=0; i<numCorner; ++i)               << 
799      {                                         << 
800        os << "                         "       << 
801           << corners[i].r << ", " << corners[i << 
802      }                                         << 
803   os << "------------------------------------- << 
804   os.precision(oldprc);                        << 
805                                                   441 
806   return os;                                   << 
807 }                                              << 
808                                                   442 
809 ////////////////////////////////////////////// << 
810 //                                                443 //
811 // Return volume                               << 444 // CreateNURBS
812                                                << 
813 G4double G4Polyhedra::GetCubicVolume()         << 
814 {                                              << 
815   if (fCubicVolume == 0.)                      << 
816   {                                            << 
817     G4double total = 0.;                       << 
818     G4int nrz = GetNumRZCorner();              << 
819     G4PolyhedraSideRZ a = GetCorner(nrz - 1);  << 
820     for (G4int i=0; i<nrz; ++i)                << 
821     {                                          << 
822       G4PolyhedraSideRZ b = GetCorner(i);      << 
823       total += (b.r*b.r + b.r*a.r + a.r*a.r)*( << 
824       a = b;                                   << 
825     }                                          << 
826     fCubicVolume = std::abs(total)*            << 
827       std::sin((GetEndPhi() - GetStartPhi())/G << 
828   }                                            << 
829   return fCubicVolume;                         << 
830 }                                              << 
831                                                << 
832 ////////////////////////////////////////////// << 
833 //                                                445 //
834 // Return surface area                         << 446 G4NURBS *G4Polyhedra::CreateNURBS() const
835                                                << 
836 G4double G4Polyhedra::GetSurfaceArea()         << 
837 {                                                 447 {
838   if (fSurfaceArea == 0.)                      << 448   return 0;
839   {                                            << 
840     G4double total = 0.;                       << 
841     G4int nrz = GetNumRZCorner();              << 
842     if (IsOpen())                              << 
843     {                                          << 
844       G4PolyhedraSideRZ a = GetCorner(nrz - 1) << 
845       for (G4int i=0; i<nrz; ++i)              << 
846       {                                        << 
847         G4PolyhedraSideRZ b = GetCorner(i);    << 
848         total += a.r*b.z - a.z*b.r;            << 
849         a = b;                                 << 
850       }                                        << 
851       total = std::abs(total);                 << 
852     }                                          << 
853     G4double alp = (GetEndPhi() - GetStartPhi( << 
854     G4double cosa = std::cos(alp);             << 
855     G4double sina = std::sin(alp);             << 
856     G4PolyhedraSideRZ a = GetCorner(nrz - 1);  << 
857     for (G4int i=0; i<nrz; ++i)                << 
858     {                                          << 
859       G4PolyhedraSideRZ b = GetCorner(i);      << 
860       G4ThreeVector p1(a.r, 0, a.z);           << 
861       G4ThreeVector p2(a.r*cosa, a.r*sina, a.z << 
862       G4ThreeVector p3(b.r*cosa, b.r*sina, b.z << 
863       G4ThreeVector p4(b.r, 0, b.z);           << 
864       total += GetNumSide()*(G4GeomTools::Quad << 
865       a = b;                                   << 
866     }                                          << 
867     fSurfaceArea = total;                      << 
868   }                                            << 
869   return fSurfaceArea;                         << 
870 }                                                 449 }
871                                                   450 
872 ////////////////////////////////////////////// << 
873 //                                             << 
874 // Set vector of surface elements, auxiliary m << 
875 // random points on surface                    << 
876                                                << 
877 void G4Polyhedra::SetSurfaceElements() const   << 
878 {                                              << 
879   fElements = new std::vector<G4Polyhedra::sur << 
880   G4double total = 0.;                         << 
881   G4int nrz = GetNumRZCorner();                << 
882                                                << 
883   // set lateral surface elements              << 
884   G4double dphi = (GetEndPhi() - GetStartPhi() << 
885   G4double cosa = std::cos(dphi);              << 
886   G4double sina = std::sin(dphi);              << 
887   G4int ia = nrz - 1;                          << 
888   for (G4int ib=0; ib<nrz; ++ib)               << 
889   {                                            << 
890     G4PolyhedraSideRZ a = GetCorner(ia);       << 
891     G4PolyhedraSideRZ b = GetCorner(ib);       << 
892     G4Polyhedra::surface_element selem;        << 
893     selem.i0 = ia;                             << 
894     selem.i1 = ib;                             << 
895     ia = ib;                                   << 
896     if (a.r == 0. && b.r == 0.) continue;      << 
897     G4ThreeVector p1(a.r, 0, a.z);             << 
898     G4ThreeVector p2(a.r*cosa, a.r*sina, a.z); << 
899     G4ThreeVector p3(b.r*cosa, b.r*sina, b.z); << 
900     G4ThreeVector p4(b.r, 0, b.z);             << 
901     if (a.r > 0.)                              << 
902     {                                          << 
903       selem.i2 = -1;                           << 
904       total += GetNumSide()*(G4GeomTools::Tria << 
905       selem.area = total;                      << 
906       fElements->push_back(selem);             << 
907     }                                          << 
908     if (b.r > 0.)                              << 
909     {                                          << 
910       selem.i2 = -2;                           << 
911       total += GetNumSide()*(G4GeomTools::Tria << 
912       selem.area = total;                      << 
913       fElements->push_back(selem);             << 
914     }                                          << 
915   }                                            << 
916                                                << 
917   // set elements for phi cuts                 << 
918   if (IsOpen())                                << 
919   {                                            << 
920     G4TwoVectorList contourRZ;                 << 
921     std::vector<G4int> triangles;              << 
922     for (G4int i=0; i<nrz; ++i)                << 
923     {                                          << 
924       G4PolyhedraSideRZ corner = GetCorner(i); << 
925       contourRZ.emplace_back(corner.r, corner. << 
926     }                                          << 
927     G4GeomTools::TriangulatePolygon(contourRZ, << 
928     auto ntria = (G4int)triangles.size();      << 
929     for (G4int i=0; i<ntria; i+=3)             << 
930     {                                          << 
931       G4Polyhedra::surface_element selem;      << 
932       selem.i0 = triangles[i];                 << 
933       selem.i1 = triangles[i+1];               << 
934       selem.i2 = triangles[i+2];               << 
935       G4PolyhedraSideRZ a = GetCorner(selem.i0 << 
936       G4PolyhedraSideRZ b = GetCorner(selem.i1 << 
937       G4PolyhedraSideRZ c = GetCorner(selem.i2 << 
938       G4double stria =                         << 
939         std::abs(G4GeomTools::TriangleArea(a.r << 
940       total += stria;                          << 
941       selem.area = total;                      << 
942       fElements->push_back(selem); // start ph << 
943       total += stria;                          << 
944       selem.area = total;                      << 
945       selem.i0 += nrz;                         << 
946       fElements->push_back(selem); // end phi  << 
947     }                                          << 
948   }                                            << 
949 }                                              << 
950                                                   451 
951 ////////////////////////////////////////////// << 
952 //                                             << 
953 // Generate random point on surface            << 
954                                                   452 
955 G4ThreeVector G4Polyhedra::GetPointOnSurface() << 
956 {                                              << 
957   // Set surface elements                      << 
958   if (fElements == nullptr)                    << 
959   {                                            << 
960     G4AutoLock l(&surface_elementsMutex);      << 
961     SetSurfaceElements();                      << 
962     l.unlock();                                << 
963   }                                            << 
964                                                << 
965   // Select surface element                    << 
966   G4Polyhedra::surface_element selem;          << 
967   selem = fElements->back();                   << 
968   G4double select = selem.area*G4QuickRand();  << 
969   auto it = std::lower_bound(fElements->begin( << 
970                              [](const G4Polyhe << 
971                              -> G4bool { retur << 
972                                                << 
973   // Generate random point                     << 
974   G4double x = 0, y = 0, z = 0;                << 
975   G4double u = G4QuickRand();                  << 
976   G4double v = G4QuickRand();                  << 
977   if (u + v > 1.) { u = 1. - u; v = 1. - v; }  << 
978   G4int i0 = (*it).i0;                         << 
979   G4int i1 = (*it).i1;                         << 
980   G4int i2 = (*it).i2;                         << 
981   if (i2 < 0) // lateral surface               << 
982   {                                            << 
983     // sample point                            << 
984     G4int nside = GetNumSide();                << 
985     G4double dphi = (GetEndPhi() - GetStartPhi << 
986     G4double cosa = std::cos(dphi);            << 
987     G4double sina = std::sin(dphi);            << 
988     G4PolyhedraSideRZ a = GetCorner(i0);       << 
989     G4PolyhedraSideRZ b = GetCorner(i1);       << 
990     G4ThreeVector p0(a.r, 0, a.z);             << 
991     G4ThreeVector p1(b.r, 0, b.z);             << 
992     G4ThreeVector p2(b.r*cosa, b.r*sina, b.z); << 
993     if (i2 == -1) p1.set(a.r*cosa, a.r*sina, a << 
994     p0 += (p1 - p0)*u + (p2 - p0)*v;           << 
995     // find selected side and rotate point     << 
996     G4double scurr = (*it).area;               << 
997     G4double sprev = (it == fElements->begin() << 
998     G4int iside = nside*(select - sprev)/(scur << 
999     if (iside == 0 && GetStartPhi() == 0.)     << 
1000     {                                         << 
1001       x = p0.x();                             << 
1002       y = p0.y();                             << 
1003       z = p0.z();                             << 
1004     }                                         << 
1005     else                                      << 
1006     {                                         << 
1007       if (iside == nside) --iside; // iside m << 
1008       G4double phi = iside*dphi + GetStartPhi << 
1009       G4double cosphi = std::cos(phi);        << 
1010       G4double sinphi = std::sin(phi);        << 
1011       x = p0.x()*cosphi - p0.y()*sinphi;      << 
1012       y = p0.x()*sinphi + p0.y()*cosphi;      << 
1013       z = p0.z();                             << 
1014     }                                         << 
1015   }                                           << 
1016   else // phi cut                             << 
1017   {                                           << 
1018     G4int nrz = GetNumRZCorner();             << 
1019     G4double phi = (i0 < nrz) ? GetStartPhi() << 
1020     if (i0 >= nrz) { i0 -= nrz; }             << 
1021     G4PolyhedraSideRZ p0 = GetCorner(i0);     << 
1022     G4PolyhedraSideRZ p1 = GetCorner(i1);     << 
1023     G4PolyhedraSideRZ p2 = GetCorner(i2);     << 
1024     G4double r = (p1.r - p0.r)*u + (p2.r - p0 << 
1025     x = r*std::cos(phi);                      << 
1026     y = r*std::sin(phi);                      << 
1027     z = (p1.z - p0.z)*u + (p2.z - p0.z)*v + p << 
1028   }                                           << 
1029   return {x, y, z};                           << 
1030 }                                             << 
1031                                                  453 
1032 ///////////////////////////////////////////// << 
1033 //                                               454 //
1034 // CreatePolyhedron                           << 455 // G4Polyhedra::G4PolyhedraHistorical stuff
1035                                               << 456 //
1036 G4Polyhedron* G4Polyhedra::CreatePolyhedron() << 457 G4Polyhedra::G4PolyhedraHistorical::~G4PolyhedraHistorical()
1037 {                                                458 {
1038   std::vector<G4TwoVector> rz(numCorner);     << 459   delete [] Z_values;
1039   for (G4int i = 0; i < numCorner; ++i)       << 460   delete [] Rmin;
1040     rz[i].set(corners[i].r, corners[i].z);    << 461   delete [] Rmax;
1041   return new G4PolyhedronPgon(startPhi, endPh << 
1042 }                                                462 }
1043                                                  463 
1044 // SetOriginalParameters                      << 464 G4Polyhedra::G4PolyhedraHistorical::G4PolyhedraHistorical( const G4PolyhedraHistorical &source )
1045 //                                            << 
1046 void G4Polyhedra::SetOriginalParameters(G4Red << 
1047 {                                                465 {
1048   G4int numPlanes = numCorner;                << 466   Start_angle   = source.Start_angle;
1049   G4bool isConvertible = true;                << 467   Opening_angle = source.Opening_angle;
1050   G4double Zmax=rz->Bmax();                   << 468   numSide   = source.numSide;
1051   rz->StartWithZMin();                        << 469   Num_z_planes  = source.Num_z_planes;
1052                                               << 470   
1053   // Prepare vectors for storage              << 471   Z_values  = new G4double[Num_z_planes];
1054   //                                          << 472   Rmin    = new G4double[Num_z_planes];
1055   std::vector<G4double> Z;                    << 473   Rmax    = new G4double[Num_z_planes];
1056   std::vector<G4double> Rmin;                 << 474   
1057   std::vector<G4double> Rmax;                 << 475   G4int i;
1058                                               << 476   for( i = 0; i < Num_z_planes; i++) {
1059   G4int countPlanes=1;                        << 477     Z_values[i] = source.Z_values[i];
1060   G4int icurr=0;                              << 478     Rmin[i]     = source.Rmin[i];
1061   G4int icurl=0;                              << 479     Rmax[i]     = source.Rmax[i];
1062                                               << 480   }
1063   // first plane Z=Z[0]                       << 
1064   //                                          << 
1065   Z.push_back(corners[0].z);                  << 
1066   G4double Zprev=Z[0];                        << 
1067   if (Zprev == corners[1].z)                  << 
1068   {                                           << 
1069     Rmin.push_back(corners[0].r);             << 
1070     Rmax.push_back (corners[1].r);icurr=1;    << 
1071   }                                           << 
1072   else if (Zprev == corners[numPlanes-1].z)   << 
1073   {                                           << 
1074     Rmin.push_back(corners[numPlanes-1].r);   << 
1075     Rmax.push_back (corners[0].r);            << 
1076     icurl=numPlanes-1;                        << 
1077   }                                           << 
1078   else                                        << 
1079   {                                           << 
1080     Rmin.push_back(corners[0].r);             << 
1081     Rmax.push_back (corners[0].r);            << 
1082   }                                           << 
1083                                               << 
1084   // next planes until last                   << 
1085   //                                          << 
1086   G4int inextr=0, inextl=0;                   << 
1087   for (G4int i=0; i < numPlanes-2; ++i)       << 
1088   {                                           << 
1089     inextr=1+icurr;                           << 
1090     inextl=(icurl <= 0)? numPlanes-1 : icurl- << 
1091                                               << 
1092     if((corners[inextr].z >= Zmax) & (corners << 
1093                                               << 
1094     G4double Zleft = corners[inextl].z;       << 
1095     G4double Zright = corners[inextr].z;      << 
1096     if(Zright>Zleft)                          << 
1097     {                                         << 
1098       Z.push_back(Zleft);                     << 
1099       countPlanes++;                          << 
1100       G4double difZr=corners[inextr].z - corn << 
1101       G4double difZl=corners[inextl].z - corn << 
1102                                               << 
1103       if(std::fabs(difZl) < kCarTolerance)    << 
1104       {                                       << 
1105         if(std::fabs(difZr) < kCarTolerance)  << 
1106         {                                     << 
1107           Rmin.push_back(corners[inextl].r);  << 
1108           Rmax.push_back(corners[icurr].r);   << 
1109         }                                     << 
1110         else                                  << 
1111         {                                     << 
1112           Rmin.push_back(corners[inextl].r);  << 
1113           Rmax.push_back(corners[icurr].r + ( << 
1114                                 *(corners[ine << 
1115         }                                     << 
1116       }                                       << 
1117       else if (difZl >= kCarTolerance)        << 
1118       {                                       << 
1119         if(std::fabs(difZr) < kCarTolerance)  << 
1120         {                                     << 
1121           Rmin.push_back(corners[icurl].r);   << 
1122           Rmax.push_back(corners[icurr].r);   << 
1123         }                                     << 
1124         else                                  << 
1125         {                                     << 
1126           Rmin.push_back(corners[icurl].r);   << 
1127           Rmax.push_back(corners[icurr].r + ( << 
1128                                 *(corners[ine << 
1129         }                                     << 
1130       }                                       << 
1131       else                                    << 
1132       {                                       << 
1133         isConvertible=false; break;           << 
1134       }                                       << 
1135       icurl=(icurl == 0)? numPlanes-1 : icurl << 
1136     }                                         << 
1137     else if(std::fabs(Zright-Zleft)<kCarToler << 
1138     {                                         << 
1139       Z.push_back(Zleft);                     << 
1140       ++countPlanes;                          << 
1141       ++icurr;                                << 
1142                                               << 
1143       icurl=(icurl == 0)? numPlanes-1 : icurl << 
1144                                               << 
1145       Rmin.push_back(corners[inextl].r);      << 
1146       Rmax.push_back (corners[inextr].r);     << 
1147     }                                         << 
1148     else  // Zright<Zleft                     << 
1149     {                                         << 
1150       Z.push_back(Zright);                    << 
1151       ++countPlanes;                          << 
1152                                               << 
1153       G4double difZr=corners[inextr].z - corn << 
1154       G4double difZl=corners[inextl].z - corn << 
1155       if(std::fabs(difZr) < kCarTolerance)    << 
1156       {                                       << 
1157         if(std::fabs(difZl) < kCarTolerance)  << 
1158         {                                     << 
1159           Rmax.push_back(corners[inextr].r);  << 
1160           Rmin.push_back(corners[icurr].r);   << 
1161         }                                     << 
1162         else                                  << 
1163         {                                     << 
1164           Rmin.push_back(corners[icurl].r + ( << 
1165                                 * (corners[in << 
1166           Rmax.push_back(corners[inextr].r);  << 
1167         }                                     << 
1168         ++icurr;                              << 
1169       }           // plate                    << 
1170       else if (difZr >= kCarTolerance)        << 
1171       {                                       << 
1172         if(std::fabs(difZl) < kCarTolerance)  << 
1173         {                                     << 
1174           Rmax.push_back(corners[inextr].r);  << 
1175           Rmin.push_back (corners[icurr].r);  << 
1176         }                                     << 
1177         else                                  << 
1178         {                                     << 
1179           Rmax.push_back(corners[inextr].r);  << 
1180           Rmin.push_back (corners[icurl].r+(Z << 
1181                                   * (corners[ << 
1182         }                                     << 
1183         ++icurr;                              << 
1184       }                                       << 
1185       else                                    << 
1186       {                                       << 
1187         isConvertible=false; break;           << 
1188       }                                       << 
1189     }                                         << 
1190   }   // end for loop                         << 
1191                                               << 
1192   // last plane Z=Zmax                        << 
1193   //                                          << 
1194   Z.push_back(Zmax);                          << 
1195   ++countPlanes;                              << 
1196   inextr=1+icurr;                             << 
1197   inextl=(icurl <= 0)? numPlanes-1 : icurl-1; << 
1198                                               << 
1199   Rmax.push_back(corners[inextr].r);          << 
1200   Rmin.push_back(corners[inextl].r);          << 
1201                                               << 
1202   // Set original parameters Rmin,Rmax,Z      << 
1203   //                                          << 
1204   if(isConvertible)                           << 
1205   {                                           << 
1206    original_parameters = new G4PolyhedraHisto << 
1207    original_parameters->numSide = numSide;    << 
1208    original_parameters->Z_values = new G4doub << 
1209    original_parameters->Rmin = new G4double[c << 
1210    original_parameters->Rmax = new G4double[c << 
1211                                               << 
1212    for(G4int j=0; j < countPlanes; ++j)       << 
1213    {                                          << 
1214      original_parameters->Z_values[j] = Z[j]; << 
1215      original_parameters->Rmax[j] = Rmax[j];  << 
1216      original_parameters->Rmin[j] = Rmin[j];  << 
1217    }                                          << 
1218    original_parameters->Start_angle = startPh << 
1219    original_parameters->Opening_angle = endPh << 
1220    original_parameters->Num_z_planes = countP << 
1221                                               << 
1222   }                                           << 
1223   else  // Set parameters(r,z) with Rmin==0 a << 
1224   {                                           << 
1225 #ifdef G4SPECSDEBUG                           << 
1226     std::ostringstream message;               << 
1227     message << "Polyhedra " << GetName() << G << 
1228       << "cannot be converted to Polyhedra wi << 
1229     G4Exception("G4Polyhedra::SetOriginalPara << 
1230                 "GeomSolids0002", JustWarning << 
1231 #endif                                        << 
1232     original_parameters = new G4PolyhedraHist << 
1233     original_parameters->numSide = numSide;   << 
1234     original_parameters->Z_values = new G4dou << 
1235     original_parameters->Rmin = new G4double[ << 
1236     original_parameters->Rmax = new G4double[ << 
1237                                               << 
1238     for(G4int j=0; j < numPlanes; ++j)        << 
1239     {                                         << 
1240       original_parameters->Z_values[j] = corn << 
1241       original_parameters->Rmax[j] = corners[ << 
1242       original_parameters->Rmin[j] = 0.0;     << 
1243     }                                         << 
1244     original_parameters->Start_angle = startP << 
1245     original_parameters->Opening_angle = endP << 
1246     original_parameters->Num_z_planes = numPl << 
1247   }                                           << 
1248 }                                                481 }
1249                                                  482 
1250 #endif                                        << 
1251                                                  483