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


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