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 5.1.p1)


  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.6 2003/03/28 09:52:50 gcosmo Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-05-01-patch-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                                                    56 
 44 #if !defined(G4GEOM_USE_UPOLYHEDRA)            << 
 45                                                << 
 46 #include "G4PolyhedraSide.hh"                      57 #include "G4PolyhedraSide.hh"
 47 #include "G4PolyPhiFace.hh"                        58 #include "G4PolyPhiFace.hh"
 48                                                    59 
 49 #include "G4GeomTools.hh"                      <<  60 #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"                  61 #include "G4EnclosingCylinder.hh"
 57 #include "G4ReduciblePolygon.hh"                   62 #include "G4ReduciblePolygon.hh"
 58 #include "G4VPVParameterisation.hh"            << 
 59                                                << 
 60 namespace                                      << 
 61 {                                              << 
 62   G4Mutex surface_elementsMutex = G4MUTEX_INIT << 
 63 }                                              << 
 64                                                << 
 65 using namespace CLHEP;                         << 
 66                                                    63 
                                                   >>  64 //
 67 // Constructor (GEANT3 style parameters)           65 // Constructor (GEANT3 style parameters)
 68 //                                                 66 //
 69 // GEANT3 PGON radii are specified in the dist     67 // GEANT3 PGON radii are specified in the distance to the norm of each face.
 70 //                                             <<  68 //  
 71 G4Polyhedra::G4Polyhedra( const G4String& name <<  69 G4Polyhedra::G4Polyhedra( const G4String& name, 
 72                                 G4double phiSt     70                                 G4double phiStart,
 73                                 G4double thePh     71                                 G4double thePhiTotal,
 74                                 G4int theNumSi <<  72                                 G4int theNumSide,  
 75                                 G4int numZPlan     73                                 G4int numZPlanes,
 76                           const G4double zPlan     74                           const G4double zPlane[],
 77                           const G4double rInne     75                           const G4double rInner[],
 78                           const G4double rOute     76                           const G4double rOuter[]  )
 79   : G4VCSGfaceted( name )                          77   : G4VCSGfaceted( name )
 80 {                                                  78 {
 81   if (theNumSide <= 0)                             79   if (theNumSide <= 0)
 82   {                                                80   {
 83     std::ostringstream message;                <<  81     G4cout << "ERROR - G4Polyhedra::G4Polyhedra(): " << GetName() << G4endl
 84     message << "Solid must have at least one s <<  82            << "        No sides specified !"
 85             << "        No sides specified !"; <<  83            << G4endl;
 86     G4Exception("G4Polyhedra::G4Polyhedra()",  <<  84     G4cerr << "ERROR - G4Polyhedra::G4Polyhedra(): " << GetName() << G4endl
 87                 FatalErrorInArgument, message) <<  85            << "        No sides specified !"
                                                   >>  86            << G4endl;
                                                   >>  87     G4Exception("G4Polyhedra::G4Polyhedra() - must have at least one side");
 88   }                                                88   }
 89                                                    89 
 90   //                                               90   //
 91   // Calculate conversion factor from G3 radiu     91   // Calculate conversion factor from G3 radius to G4 radius
 92   //                                               92   //
 93   G4double phiTotal = thePhiTotal;                 93   G4double phiTotal = thePhiTotal;
 94   if ( (phiTotal <=0) || (phiTotal >= twopi*(1 <<  94   if ( (phiTotal <=0) || (phiTotal >= 2*M_PI*(1-DBL_EPSILON)) )
 95     { phiTotal = twopi; }                      <<  95     phiTotal = 2*M_PI;
 96   G4double convertRad = std::cos(0.5*phiTotal/ <<  96   G4double convertRad = cos(0.5*phiTotal/theNumSide);
 97                                                    97 
 98   //                                               98   //
 99   // Some historical stuff                         99   // Some historical stuff
100   //                                              100   //
101   original_parameters = new G4PolyhedraHistori    101   original_parameters = new G4PolyhedraHistorical;
102                                                << 102   
103   original_parameters->numSide = theNumSide;      103   original_parameters->numSide = theNumSide;
104   original_parameters->Start_angle = phiStart;    104   original_parameters->Start_angle = phiStart;
105   original_parameters->Opening_angle = phiTota    105   original_parameters->Opening_angle = phiTotal;
106   original_parameters->Num_z_planes = numZPlan    106   original_parameters->Num_z_planes = numZPlanes;
107   original_parameters->Z_values = new G4double    107   original_parameters->Z_values = new G4double[numZPlanes];
108   original_parameters->Rmin = new G4double[num    108   original_parameters->Rmin = new G4double[numZPlanes];
109   original_parameters->Rmax = new G4double[num    109   original_parameters->Rmax = new G4double[numZPlanes];
110                                                   110 
111   for (G4int i=0; i<numZPlanes; ++i)           << 111   G4int i;
                                                   >> 112   for (i=0; i<numZPlanes; i++)
112   {                                               113   {
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[    114     original_parameters->Z_values[i] = zPlane[i];
132     original_parameters->Rmin[i] = rInner[i]/c    115     original_parameters->Rmin[i] = rInner[i]/convertRad;
133     original_parameters->Rmax[i] = rOuter[i]/c    116     original_parameters->Rmax[i] = rOuter[i]/convertRad;
134   }                                               117   }
135                                                << 118   
136                                                << 119   
137   //                                              120   //
138   // Build RZ polygon using special PCON/PGON     121   // Build RZ polygon using special PCON/PGON GEANT3 constructor
139   //                                              122   //
140   auto rz = new G4ReduciblePolygon( rInner, rO << 123   G4ReduciblePolygon *rz =
                                                   >> 124     new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
141   rz->ScaleA( 1/convertRad );                     125   rz->ScaleA( 1/convertRad );
142                                                << 126   
143   //                                              127   //
144   // Do the real work                             128   // Do the real work
145   //                                              129   //
146   Create( phiStart, phiTotal, theNumSide, rz )    130   Create( phiStart, phiTotal, theNumSide, rz );
147                                                << 131   
148   delete rz;                                      132   delete rz;
149 }                                                 133 }
150                                                   134 
                                                   >> 135 
                                                   >> 136 //
151 // Constructor (generic parameters)               137 // Constructor (generic parameters)
152 //                                                138 //
153 G4Polyhedra::G4Polyhedra( const G4String& name << 139 G4Polyhedra::G4Polyhedra( const G4String& name, 
154                                 G4double phiSt    140                                 G4double phiStart,
155                                 G4double phiTo    141                                 G4double phiTotal,
156                                 G4int    theNu << 142                                 G4int    theNumSide,  
157                                 G4int    numRZ    143                                 G4int    numRZ,
158                           const G4double r[],     144                           const G4double r[],
159                           const G4double z[]      145                           const G4double z[]   )
160   : G4VCSGfaceted( name ), genericPgon(true)   << 146   : G4VCSGfaceted( name )
161 {                                                 147 {
162   if (theNumSide <= 0)                         << 148   original_parameters = 0;
163   {                                            << 149   
164     std::ostringstream message;                << 150   G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
165     message << "Solid must have at least one s << 151   
166             << "        No sides specified !"; << 
167     G4Exception("G4Polyhedra::G4Polyhedra()",  << 
168                 FatalErrorInArgument, message) << 
169   }                                            << 
170                                                << 
171   auto rz = new G4ReduciblePolygon( r, z, numR << 
172                                                << 
173   Create( phiStart, phiTotal, theNumSide, rz )    152   Create( phiStart, phiTotal, theNumSide, rz );
174                                                << 153   
175   // Set original_parameters struct for consis << 
176   //                                           << 
177   SetOriginalParameters(rz);                   << 
178                                                << 
179   delete rz;                                      154   delete rz;
180 }                                                 155 }
181                                                   156 
                                                   >> 157 
                                                   >> 158 //
182 // Create                                         159 // Create
183 //                                                160 //
184 // Generic create routine, called by each cons    161 // Generic create routine, called by each constructor
185 // after conversion of arguments                  162 // after conversion of arguments
186 //                                                163 //
187 void G4Polyhedra::Create( G4double phiStart,      164 void G4Polyhedra::Create( G4double phiStart,
188                           G4double phiTotal,      165                           G4double phiTotal,
189                           G4int    theNumSide, << 166                           G4int    theNumSide,  
190                           G4ReduciblePolygon*  << 167                           G4ReduciblePolygon *rz  )
191 {                                                 168 {
192   //                                              169   //
193   // Perform checks of rz values                  170   // Perform checks of rz values
194   //                                              171   //
195   if (rz->Amin() < 0.0)                           172   if (rz->Amin() < 0.0)
196   {                                               173   {
197     std::ostringstream message;                << 174     G4cout << "ERROR - G4Polyhedra::Create() " << GetName() << G4endl
198     message << "Illegal input parameters - " < << 175            << "        All R values must be >= 0 !"
199             << "        All R values must be > << 176            << G4endl;
200     G4Exception("G4Polyhedra::Create()", "Geom << 177     G4cerr << "ERROR - G4Polyhedra::Create() " << GetName() << G4endl
201                 FatalErrorInArgument, message) << 178            << "        All R values must be >= 0 !"
                                                   >> 179            << G4endl;
                                                   >> 180     G4Exception("G4Polyhedra::Create() - Illegal input parameters");
202   }                                               181   }
203                                                   182 
204   G4double rzArea = rz->Area();                   183   G4double rzArea = rz->Area();
205   if (rzArea < -kCarTolerance)                    184   if (rzArea < -kCarTolerance)
206   {                                            << 
207     rz->ReverseOrder();                           185     rz->ReverseOrder();
208   }                                            << 
209   else if (rzArea < kCarTolerance)             << 
210   {                                            << 
211     std::ostringstream message;                << 
212     message << "Illegal input parameters - " < << 
213             << "        R/Z cross section is z << 
214     G4Exception("G4Polyhedra::Create()", "Geom << 
215                 FatalErrorInArgument, message) << 
216   }                                            << 
217                                                   186 
218   if ( (!rz->RemoveDuplicateVertices( kCarTole << 187   else if (rzArea < -kCarTolerance)
219     || (!rz->RemoveRedundantVertices( kCarTole << 
220   {                                               188   {
221     std::ostringstream message;                << 189     G4cout << "ERROR - G4Polyhedra::Create() " << GetName() << G4endl
222     message << "Illegal input parameters - " < << 190            << "        R/Z cross section is zero or near zero: "
223             << "        Too few unique R/Z val << 191            << rzArea << G4endl;
224     G4Exception("G4Polyhedra::Create()", "Geom << 192     G4cerr << "ERROR - G4Polyhedra::Create() " << GetName() << G4endl
225                 FatalErrorInArgument, message) << 193            << "        R/Z cross section is zero or near zero: "
                                                   >> 194            << rzArea << G4endl;
                                                   >> 195     G4Exception("G4Polyhedra::Create() - Illegal input parameters");
226   }                                               196   }
227                                                << 197     
228   if (rz->CrossesItself( 1/kInfinity ))        << 198   if ( (!rz->RemoveDuplicateVertices( kCarTolerance ))
                                                   >> 199     || (!rz->RemoveRedundantVertices( kCarTolerance )) ) 
229   {                                               200   {
230     std::ostringstream message;                << 201     G4cout << "ERROR - G4Polyhedra::Create() " << GetName() << G4endl
231     message << "Illegal input parameters - " < << 202            << "        Too few unique R/Z values !"
232             << "        R/Z segments cross !"; << 203            << G4endl;
233     G4Exception("G4Polyhedra::Create()", "Geom << 204     G4cerr << "ERROR - G4Polyhedra::Create() " << GetName() << G4endl
234                 FatalErrorInArgument, message) << 205            << "        Too few unique R/Z values !"
                                                   >> 206            << G4endl;
                                                   >> 207     G4Exception("G4Polyhedra::Create() - Illegal input parameters");
                                                   >> 208   }
                                                   >> 209 
                                                   >> 210   if (rz->CrossesItself( 1/kInfinity )) 
                                                   >> 211   {
                                                   >> 212     G4cout << "ERROR - G4Polyhedra::Create() " << GetName() << G4endl
                                                   >> 213            << "        R/Z segments cross !"
                                                   >> 214            << G4endl;
                                                   >> 215     G4cerr << "ERROR - G4Polyhedra::Create() " << GetName() << G4endl
                                                   >> 216            << "        R/Z segments cross !"
                                                   >> 217            << G4endl;
                                                   >> 218     G4Exception("G4Polyhedra::Create() - Illegal input parameters");
235   }                                               219   }
236                                                   220 
237   numCorner = rz->NumVertices();                  221   numCorner = rz->NumVertices();
238                                                   222 
239                                                   223 
240   startPhi = phiStart;                            224   startPhi = phiStart;
241   while( startPhi < 0 )    // Loop checking, 1 << 225   while( startPhi < 0 ) startPhi += 2*M_PI;
242     startPhi += twopi;                         << 
243   //                                              226   //
244   // Phi opening? Account for some possible ro    227   // Phi opening? Account for some possible roundoff, and interpret
245   // nonsense value as representing no phi ope    228   // nonsense value as representing no phi opening
246   //                                              229   //
247   if ( (phiTotal <= 0) || (phiTotal > twopi*(1 << 230   if ( (phiTotal <= 0) || (phiTotal > 2.0*M_PI*(1-DBL_EPSILON)) )
248   {                                               231   {
249     phiIsOpen = false;                            232     phiIsOpen = false;
250     endPhi = startPhi + twopi;                 << 233     endPhi = phiStart+2*M_PI;
251   }                                               234   }
252   else                                            235   else
253   {                                               236   {
254     phiIsOpen = true;                             237     phiIsOpen = true;
255     endPhi = startPhi + phiTotal;              << 238     
                                                   >> 239     //
                                                   >> 240     // Convert phi into our convention
                                                   >> 241     //
                                                   >> 242     endPhi = phiStart+phiTotal;
                                                   >> 243     while( endPhi < startPhi ) endPhi += 2*M_PI;
256   }                                               244   }
257                                                << 245   
258   //                                              246   //
259   // Save number sides                            247   // Save number sides
260   //                                              248   //
261   numSide = theNumSide;                           249   numSide = theNumSide;
262                                                << 250   
263   //                                              251   //
264   // Allocate corner array.                    << 252   // Allocate corner array. 
265   //                                              253   //
266   corners = new G4PolyhedraSideRZ[numCorner];     254   corners = new G4PolyhedraSideRZ[numCorner];
267                                                   255 
268   //                                              256   //
269   // Copy corners                                 257   // Copy corners
270   //                                              258   //
271   G4ReduciblePolygonIterator iterRZ(rz);          259   G4ReduciblePolygonIterator iterRZ(rz);
272                                                << 260   
273   G4PolyhedraSideRZ *next = corners;              261   G4PolyhedraSideRZ *next = corners;
274   iterRZ.Begin();                                 262   iterRZ.Begin();
275   do    // Loop checking, 13.08.2015, G.Cosmo  << 263   do
276   {                                               264   {
277     next->r = iterRZ.GetA();                      265     next->r = iterRZ.GetA();
278     next->z = iterRZ.GetB();                      266     next->z = iterRZ.GetB();
279   } while( ++next, iterRZ.Next() );               267   } while( ++next, iterRZ.Next() );
280                                                << 268   
281   //                                              269   //
282   // Allocate face pointer array                  270   // Allocate face pointer array
283   //                                              271   //
284   numFace = phiIsOpen ? numCorner+2 : numCorne    272   numFace = phiIsOpen ? numCorner+2 : numCorner;
285   faces = new G4VCSGface*[numFace];               273   faces = new G4VCSGface*[numFace];
286                                                << 274   
287   //                                              275   //
288   // Construct side faces                         276   // Construct side faces
289   //                                              277   //
290   // To do so properly, we need to keep track     278   // To do so properly, we need to keep track of four successive RZ
291   // corners.                                     279   // corners.
292   //                                              280   //
293   // But! Don't construct a face if both point    281   // But! Don't construct a face if both points are at zero radius!
294   //                                              282   //
295   G4PolyhedraSideRZ* corner = corners,         << 283   G4PolyhedraSideRZ *corner = corners,
296                    * prev = corners + numCorne << 284                     *prev = corners + numCorner-1,
297                    * nextNext;                 << 285                     *nextNext;
298   G4VCSGface** face = faces;                   << 286   G4VCSGface   **face = faces;
299   do    // Loop checking, 13.08.2015, G.Cosmo  << 287   do
300   {                                               288   {
301     next = corner+1;                              289     next = corner+1;
302     if (next >= corners+numCorner) next = corn    290     if (next >= corners+numCorner) next = corners;
303     nextNext = next+1;                            291     nextNext = next+1;
304     if (nextNext >= corners+numCorner) nextNex    292     if (nextNext >= corners+numCorner) nextNext = corners;
305                                                << 293     
306     if (corner->r < 1/kInfinity && next->r < 1    294     if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
307 /*                                             << 295 
                                                   >> 296     //
308     // We must decide here if we can dare decl    297     // We must decide here if we can dare declare one of our faces
309     // as having a "valid" normal (i.e. allBeh    298     // as having a "valid" normal (i.e. allBehind = true). This
310     // is never possible if the face faces "in    299     // is never possible if the face faces "inward" in r *unless*
311     // we have only one side                      300     // we have only one side
312     //                                            301     //
313     G4bool allBehind;                             302     G4bool allBehind;
314     if ((corner->z > next->z) && (numSide > 1)    303     if ((corner->z > next->z) && (numSide > 1))
315     {                                             304     {
316       allBehind = false;                          305       allBehind = false;
317     }                                             306     }
318     else                                          307     else
319     {                                             308     {
320       //                                          309       //
321       // Otherwise, it is only true if the lin    310       // Otherwise, it is only true if the line passing
322       // through the two points of the segment    311       // through the two points of the segment do not
323       // split the r/z cross section              312       // split the r/z cross section
324       //                                          313       //
325       allBehind = !rz->BisectedBy( corner->r,     314       allBehind = !rz->BisectedBy( corner->r, corner->z,
326                                    next->r, ne    315                                    next->r, next->z, kCarTolerance );
327     }                                             316     }
328 */                                             << 317     
329     *face++ = new G4PolyhedraSide( prev, corne    318     *face++ = new G4PolyhedraSide( prev, corner, next, nextNext,
330                  numSide, startPhi, endPhi-sta    319                  numSide, startPhi, endPhi-startPhi, phiIsOpen );
331   } while( prev=corner, corner=next, corner >     320   } while( prev=corner, corner=next, corner > corners );
332                                                << 321   
333   if (phiIsOpen)                                  322   if (phiIsOpen)
334   {                                               323   {
335     //                                            324     //
336     // Construct phi open edges                   325     // Construct phi open edges
337     //                                            326     //
338     *face++ = new G4PolyPhiFace( rz, startPhi,    327     *face++ = new G4PolyPhiFace( rz, startPhi, phiTotal/numSide, endPhi );
339     *face++ = new G4PolyPhiFace( rz, endPhi,      328     *face++ = new G4PolyPhiFace( rz, endPhi,   phiTotal/numSide, startPhi );
340   }                                               329   }
341                                                << 330   
342   //                                              331   //
343   // We might have dropped a face or two: reca    332   // We might have dropped a face or two: recalculate numFace
344   //                                              333   //
345   numFace = (G4int)(face-faces);               << 334   numFace = face-faces;
346                                                << 335   
347   //                                              336   //
348   // Make enclosingCylinder                       337   // Make enclosingCylinder
349   //                                              338   //
350   enclosingCylinder =                             339   enclosingCylinder =
351     new G4EnclosingCylinder( rz, phiIsOpen, ph    340     new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
352 }                                                 341 }
353                                                   342 
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                                                   343 
                                                   >> 344 //
362 // Destructor                                     345 // Destructor
363 //                                                346 //
364 G4Polyhedra::~G4Polyhedra()                       347 G4Polyhedra::~G4Polyhedra()
365 {                                                 348 {
366   delete [] corners;                              349   delete [] corners;
367   delete original_parameters;                  << 350   if (original_parameters) delete original_parameters;
                                                   >> 351   
368   delete enclosingCylinder;                       352   delete enclosingCylinder;
369   delete fElements;                            << 
370   delete fpPolyhedron;                         << 
371   corners = nullptr;                           << 
372   original_parameters = nullptr;               << 
373   enclosingCylinder = nullptr;                 << 
374   fElements = nullptr;                         << 
375   fpPolyhedron = nullptr;                      << 
376 }                                                 353 }
377                                                   354 
                                                   >> 355 
                                                   >> 356 //
378 // Copy constructor                               357 // Copy constructor
379 //                                                358 //
380 G4Polyhedra::G4Polyhedra( const G4Polyhedra& s << 359 G4Polyhedra::G4Polyhedra( const G4Polyhedra &source )
381   : G4VCSGfaceted( source )                       360   : G4VCSGfaceted( source )
382 {                                                 361 {
383   CopyStuff( source );                            362   CopyStuff( source );
384 }                                                 363 }
385                                                   364 
                                                   >> 365 
                                                   >> 366 //
386 // Assignment operator                            367 // Assignment operator
387 //                                                368 //
388 G4Polyhedra &G4Polyhedra::operator=( const G4P << 369 const G4Polyhedra &G4Polyhedra::operator=( const G4Polyhedra &source )
389 {                                                 370 {
390   if (this == &source) return *this;              371   if (this == &source) return *this;
391                                                   372 
392   G4VCSGfaceted::operator=( source );             373   G4VCSGfaceted::operator=( source );
393                                                << 374   
394   delete [] corners;                              375   delete [] corners;
395   delete original_parameters;                  << 376   if (original_parameters) delete original_parameters;
                                                   >> 377   
396   delete enclosingCylinder;                       378   delete enclosingCylinder;
397                                                << 379   
398   CopyStuff( source );                            380   CopyStuff( source );
399                                                << 381   
400   return *this;                                   382   return *this;
401 }                                                 383 }
402                                                   384 
                                                   >> 385 
                                                   >> 386 //
403 // CopyStuff                                      387 // CopyStuff
404 //                                                388 //
405 void G4Polyhedra::CopyStuff( const G4Polyhedra << 389 void G4Polyhedra::CopyStuff( const G4Polyhedra &source )
406 {                                                 390 {
407   //                                              391   //
408   // Simple stuff                                 392   // Simple stuff
409   //                                              393   //
410   numSide    = source.numSide;                    394   numSide    = source.numSide;
411   startPhi   = source.startPhi;                << 395   startPhi  = source.startPhi;
412   endPhi     = source.endPhi;                  << 396   endPhi    = source.endPhi;
413   phiIsOpen  = source.phiIsOpen;                  397   phiIsOpen  = source.phiIsOpen;
414   numCorner  = source.numCorner;                  398   numCorner  = source.numCorner;
415   genericPgon= source.genericPgon;             << 
416                                                   399 
417   //                                              400   //
418   // The corner array                             401   // The corner array
419   //                                              402   //
420   corners = new G4PolyhedraSideRZ[numCorner];     403   corners = new G4PolyhedraSideRZ[numCorner];
421                                                << 404   
422   G4PolyhedraSideRZ* corn = corners,           << 405   G4PolyhedraSideRZ  *corn = corners,
423                    * sourceCorn = source.corne << 406         *sourceCorn = source.corners;
424   do    // Loop checking, 13.08.2015, G.Cosmo  << 407   do
425   {                                               408   {
426     *corn = *sourceCorn;                          409     *corn = *sourceCorn;
427   } while( ++sourceCorn, ++corn < corners+numC    410   } while( ++sourceCorn, ++corn < corners+numCorner );
428                                                << 411   
429   //                                              412   //
430   // Original parameters                          413   // Original parameters
431   //                                              414   //
432   if (source.original_parameters != nullptr)   << 415   if (source.original_parameters)
433   {                                               416   {
434     original_parameters =                         417     original_parameters =
435       new G4PolyhedraHistorical( *source.origi    418       new G4PolyhedraHistorical( *source.original_parameters );
436   }                                               419   }
437                                                << 420   
438   //                                              421   //
439   // Enclosing cylinder                           422   // Enclosing cylinder
440   //                                              423   //
441   enclosingCylinder = new G4EnclosingCylinder(    424   enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder );
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 }                                                 425 }
456                                                   426 
457 // Reset                                       << 
458 //                                             << 
459 // Recalculates and reshapes the solid, given  << 
460 // original_parameters.                        << 
461 //                                             << 
462 G4bool G4Polyhedra::Reset()                    << 
463 {                                              << 
464   if (genericPgon)                             << 
465   {                                            << 
466     std::ostringstream message;                << 
467     message << "Solid " << GetName() << " buil << 
468             << G4endl << "Not applicable to th << 
469     G4Exception("G4Polyhedra::Reset()", "GeomS << 
470                 JustWarning, message, "Paramet << 
471     return true;                               << 
472   }                                            << 
473                                                << 
474   //                                           << 
475   // Clear old setup                           << 
476   //                                           << 
477   G4VCSGfaceted::DeleteStuff();                << 
478   delete [] corners;                           << 
479   delete enclosingCylinder;                    << 
480   delete fElements;                            << 
481   corners = nullptr;                           << 
482   fElements = nullptr;                         << 
483   enclosingCylinder = nullptr;                 << 
484                                                << 
485   //                                           << 
486   // Rebuild polyhedra                         << 
487   //                                           << 
488   auto rz = new G4ReduciblePolygon( original_p << 
489                                     original_p << 
490                                     original_p << 
491                                     original_p << 
492   Create( original_parameters->Start_angle,    << 
493           original_parameters->Opening_angle,  << 
494           original_parameters->numSide, rz );  << 
495   delete rz;                                   << 
496                                                << 
497   return false;                                << 
498 }                                              << 
499                                                   427 
                                                   >> 428 //
500 // Inside                                         429 // Inside
501 //                                                430 //
502 // This is an override of G4VCSGfaceted::Insid    431 // This is an override of G4VCSGfaceted::Inside, created in order
503 // to speed things up by first checking with G    432 // to speed things up by first checking with G4EnclosingCylinder.
504 //                                                433 //
505 EInside G4Polyhedra::Inside( const G4ThreeVect << 434 EInside G4Polyhedra::Inside( const G4ThreeVector &p ) const
506 {                                                 435 {
507   //                                              436   //
508   // Quick test                                   437   // Quick test
509   //                                              438   //
510   if (enclosingCylinder->MustBeOutside(p)) ret    439   if (enclosingCylinder->MustBeOutside(p)) return kOutside;
511                                                   440 
512   //                                              441   //
513   // Long answer                                  442   // Long answer
514   //                                              443   //
515   return G4VCSGfaceted::Inside(p);                444   return G4VCSGfaceted::Inside(p);
516 }                                                 445 }
517                                                   446 
                                                   >> 447 
                                                   >> 448 //
518 // DistanceToIn                                   449 // DistanceToIn
519 //                                                450 //
520 // This is an override of G4VCSGfaceted::Insid    451 // This is an override of G4VCSGfaceted::Inside, created in order
521 // to speed things up by first checking with G    452 // to speed things up by first checking with G4EnclosingCylinder.
522 //                                                453 //
523 G4double G4Polyhedra::DistanceToIn( const G4Th << 454 G4double G4Polyhedra::DistanceToIn( const G4ThreeVector &p,
524                                     const G4Th << 455                                     const G4ThreeVector &v ) const
525 {                                                 456 {
526   //                                              457   //
527   // Quick test                                   458   // Quick test
528   //                                              459   //
529   if (enclosingCylinder->ShouldMiss(p,v))         460   if (enclosingCylinder->ShouldMiss(p,v))
530     return kInfinity;                             461     return kInfinity;
531                                                << 462   
532   //                                              463   //
533   // Long answer                                  464   // Long answer
534   //                                              465   //
535   return G4VCSGfaceted::DistanceToIn( p, v );     466   return G4VCSGfaceted::DistanceToIn( p, v );
536 }                                                 467 }
537                                                   468 
                                                   >> 469 
                                                   >> 470 //
538 // DistanceToIn                                   471 // DistanceToIn
539 //                                                472 //
540 G4double G4Polyhedra::DistanceToIn( const G4Th << 473 G4double G4Polyhedra::DistanceToIn( const G4ThreeVector &p ) const
541 {                                                 474 {
542   return G4VCSGfaceted::DistanceToIn(p);          475   return G4VCSGfaceted::DistanceToIn(p);
543 }                                                 476 }
544                                                   477 
545 // Get bounding box                            << 
546 //                                             << 
547 void G4Polyhedra::BoundingLimits(G4ThreeVector << 
548                                  G4ThreeVector << 
549 {                                              << 
550   G4double rmin = kInfinity, rmax = -kInfinity << 
551   G4double zmin = kInfinity, zmax = -kInfinity << 
552   for (G4int i=0; i<GetNumRZCorner(); ++i)     << 
553   {                                            << 
554     G4PolyhedraSideRZ corner = GetCorner(i);   << 
555     if (corner.r < rmin) rmin = corner.r;      << 
556     if (corner.r > rmax) rmax = corner.r;      << 
557     if (corner.z < zmin) zmin = corner.z;      << 
558     if (corner.z > zmax) zmax = corner.z;      << 
559   }                                            << 
560                                                   478 
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 }                                              << 
612                                                << 
613 // Calculate extent under transform and specif << 
614 //                                                479 //
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                                                << 
729 // ComputeDimensions                              480 // ComputeDimensions
730 //                                                481 //
731 void G4Polyhedra::ComputeDimensions(       G4V << 482 void G4Polyhedra::ComputeDimensions(       G4VPVParameterisation*,
732                                      const G4i << 483                                      const G4int,
733                                      const G4V << 484                                      const G4VPhysicalVolume* )
734 {                                                 485 {
735   p->ComputeDimensions(*this,n,pRep);          << 
736 }                                                 486 }
737                                                   487 
                                                   >> 488 
                                                   >> 489 //
738 // GetEntityType                                  490 // GetEntityType
739 //                                                491 //
740 G4GeometryType G4Polyhedra::GetEntityType() co    492 G4GeometryType G4Polyhedra::GetEntityType() const
741 {                                                 493 {
742   return {"G4Polyhedra"};                      << 494   return G4String("G4Polyhedra");
743 }                                                 495 }
744                                                   496 
745 // IsFaceted                                   << 
746 //                                             << 
747 G4bool G4Polyhedra::IsFaceted() const          << 
748 {                                              << 
749   return true;                                 << 
750 }                                              << 
751                                                   497 
752 // Make a clone of the object                  << 
753 //                                                498 //
754 G4VSolid* G4Polyhedra::Clone() const           << 
755 {                                              << 
756   return new G4Polyhedra(*this);               << 
757 }                                              << 
758                                                << 
759 // Stream object contents to an output stream     499 // Stream object contents to an output stream
760 //                                                500 //
761 std::ostream& G4Polyhedra::StreamInfo( std::os << 501 G4std::ostream& G4Polyhedra::StreamInfo( G4std::ostream& os ) const
762 {                                                 502 {
763   G4long oldprc = os.precision(16);            << 
764   os << "-------------------------------------    503   os << "-----------------------------------------------------------\n"
765      << "    *** Dump for solid - " << GetName    504      << "    *** Dump for solid - " << GetName() << " ***\n"
766      << "    =================================    505      << "    ===================================================\n"
767      << " Solid type: G4Polyhedra\n"              506      << " Solid type: G4Polyhedra\n"
768      << " Parameters: \n"                         507      << " Parameters: \n"
769      << "    starting phi angle : " << startPh    508      << "    starting phi angle : " << startPhi/degree << " degrees \n"
770      << "    ending phi angle   : " << endPhi/ << 509      << "    ending phi angle   : " << endPhi/degree << " degrees \n";
771      << "    number of sides    : " << numSide << 
772   G4int i=0;                                      510   G4int i=0;
773   if (!genericPgon)                            << 511   if (original_parameters)
774   {                                               512   {
775     G4int numPlanes = original_parameters->Num    513     G4int numPlanes = original_parameters->Num_z_planes;
776     os << "    number of Z planes: " << numPla    514     os << "    number of Z planes: " << numPlanes << "\n"
777        << "              Z values: \n";           515        << "              Z values: \n";
778     for (i=0; i<numPlanes; ++i)                << 516     for (i=0; i<numPlanes; i++)
779     {                                             517     {
780       os << "              Z plane " << i << "    518       os << "              Z plane " << i << ": "
781          << original_parameters->Z_values[i] <    519          << original_parameters->Z_values[i] << "\n";
782     }                                             520     }
783     os << "              Tangent distances to  << 521     os << "              Tangent distances to inner surface: \n";
784     for (i=0; i<numPlanes; ++i)                << 522     for (i=0; i<numPlanes; i++)
785     {                                             523     {
786       os << "              Z plane " << i << "    524       os << "              Z plane " << i << ": "
787          << original_parameters->Rmin[i] << "\    525          << original_parameters->Rmin[i] << "\n";
788     }                                             526     }
789     os << "              Tangent distances to  << 527     os << "              Tangent distances to outer surface: \n";
790     for (i=0; i<numPlanes; ++i)                << 528     for (i=0; i<numPlanes; i++)
791     {                                             529     {
792       os << "              Z plane " << i << "    530       os << "              Z plane " << i << ": "
793          << original_parameters->Rmax[i] << "\    531          << original_parameters->Rmax[i] << "\n";
794     }                                             532     }
795   }                                               533   }
796   os << "    number of RZ points: " << numCorn    534   os << "    number of RZ points: " << numCorner << "\n"
797      << "              RZ values (corners): \n << 535      << "              RZ values: \n";
798      for (i=0; i<numCorner; ++i)               << 536      for (i=0; i<numCorner; i++)
799      {                                            537      {
800        os << "                         "          538        os << "                         "
801           << corners[i].r << ", " << corners[i    539           << corners[i].r << ", " << corners[i].z << "\n";
802      }                                            540      }
803   os << "-------------------------------------    541   os << "-----------------------------------------------------------\n";
804   os.precision(oldprc);                        << 
805                                                   542 
806   return os;                                      543   return os;
807 }                                                 544 }
808                                                   545 
809 ////////////////////////////////////////////// << 
810 //                                             << 
811 // Return volume                               << 
812                                                   546 
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 //                                                547 //
834 // Return surface area                         << 548 // CreatePolyhedron
835                                                << 
836 G4double G4Polyhedra::GetSurfaceArea()         << 
837 {                                              << 
838   if (fSurfaceArea == 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 }                                              << 
871                                                << 
872 ////////////////////////////////////////////// << 
873 //                                                549 //
874 // Set vector of surface elements, auxiliary m << 550 G4Polyhedron* G4Polyhedra::CreatePolyhedron() const
875 // random points on surface                    << 551 { 
876                                                << 552   //
877 void G4Polyhedra::SetSurfaceElements() const   << 553   // This has to be fixed in visualization. Fake it for the moment.
878 {                                              << 554   // 
879   fElements = new std::vector<G4Polyhedra::sur << 555   if (original_parameters)
880   G4double total = 0.;                         << 556   { 
881   G4int nrz = GetNumRZCorner();                << 557     return new G4PolyhedronPgon( original_parameters->Start_angle,
882                                                << 558                                  original_parameters->Opening_angle,
883   // set lateral surface elements              << 559                                  original_parameters->numSide,
884   G4double dphi = (GetEndPhi() - GetStartPhi() << 560                                  original_parameters->Num_z_planes,
885   G4double cosa = std::cos(dphi);              << 561                                  original_parameters->Z_values,
886   G4double sina = std::sin(dphi);              << 562                                  original_parameters->Rmin,
887   G4int ia = nrz - 1;                          << 563                                  original_parameters->Rmax);
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   }                                               564   }
916                                                << 565   else
917   // set elements for phi cuts                 << 
918   if (IsOpen())                                << 
919   {                                               566   {
920     G4TwoVectorList contourRZ;                 << 567     G4cerr << "ERROR - G4Polyhedra::CreatePolyhedron() " << GetName() << G4endl
921     std::vector<G4int> triangles;              << 568            << "        Visualization of this type of G4Polycone" << G4endl
922     for (G4int i=0; i<nrz; ++i)                << 569            << "        is not supported at this time !" << G4endl;
923     {                                          << 570     return 0;
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   }                                               571   }
949 }                                              << 572 }  
950                                                   573 
951 ////////////////////////////////////////////// << 
952 //                                             << 
953 // Generate random point on surface            << 
954                                                   574 
955 G4ThreeVector G4Polyhedra::GetPointOnSurface() << 575 //
                                                   >> 576 // CreateNURBS
                                                   >> 577 //
                                                   >> 578 G4NURBS *G4Polyhedra::CreateNURBS() const
956 {                                                 579 {
957   // Set surface elements                      << 580   return 0;
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 }                                                581 }
1031                                                  582 
1032 ///////////////////////////////////////////// << 
1033 //                                            << 
1034 // CreatePolyhedron                           << 
1035                                                  583 
1036 G4Polyhedron* G4Polyhedra::CreatePolyhedron() << 584 //
                                                   >> 585 // G4Polyhedra::G4PolyhedraHistorical stuff
                                                   >> 586 //
                                                   >> 587 G4Polyhedra::
                                                   >> 588 G4PolyhedraHistorical::G4PolyhedraHistorical()
1037 {                                                589 {
1038   std::vector<G4TwoVector> rz(numCorner);     << 
1039   for (G4int i = 0; i < numCorner; ++i)       << 
1040     rz[i].set(corners[i].r, corners[i].z);    << 
1041   return new G4PolyhedronPgon(startPhi, endPh << 
1042 }                                                590 }
1043                                                  591 
1044 // SetOriginalParameters                      << 592 G4Polyhedra::
1045 //                                            << 593 G4PolyhedraHistorical::~G4PolyhedraHistorical()
1046 void G4Polyhedra::SetOriginalParameters(G4Red << 
1047 {                                                594 {
1048   G4int numPlanes = numCorner;                << 595   delete [] Z_values;
1049   G4bool isConvertible = true;                << 596   delete [] Rmin;
1050   G4double Zmax=rz->Bmax();                   << 597   delete [] Rmax;
1051   rz->StartWithZMin();                        << 598 }
1052                                               << 
1053   // Prepare vectors for storage              << 
1054   //                                          << 
1055   std::vector<G4double> Z;                    << 
1056   std::vector<G4double> Rmin;                 << 
1057   std::vector<G4double> Rmax;                 << 
1058                                               << 
1059   G4int countPlanes=1;                        << 
1060   G4int icurr=0;                              << 
1061   G4int icurl=0;                              << 
1062                                               << 
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                                                  599 
1222   }                                           << 600 G4Polyhedra::
1223   else  // Set parameters(r,z) with Rmin==0 a << 601 G4PolyhedraHistorical::
                                                   >> 602 G4PolyhedraHistorical( const G4PolyhedraHistorical &source )
                                                   >> 603 {
                                                   >> 604   Start_angle   = source.Start_angle;
                                                   >> 605   Opening_angle = source.Opening_angle;
                                                   >> 606   numSide       = source.numSide;
                                                   >> 607   Num_z_planes  = source.Num_z_planes;
                                                   >> 608   
                                                   >> 609   Z_values = new G4double[Num_z_planes];
                                                   >> 610   Rmin     = new G4double[Num_z_planes];
                                                   >> 611   Rmax     = new G4double[Num_z_planes];
                                                   >> 612   
                                                   >> 613   G4int i;
                                                   >> 614   for( i = 0; i < Num_z_planes; i++)
1224   {                                              615   {
1225 #ifdef G4SPECSDEBUG                           << 616     Z_values[i]  = source.Z_values[i];
1226     std::ostringstream message;               << 617     Rmin[i]      = source.Rmin[i];
1227     message << "Polyhedra " << GetName() << G << 618     Rmax[i]      = source.Rmax[i];
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   }                                              619   }
1248 }                                                620 }
1249                                               << 
1250 #endif                                        << 
1251                                                  621