Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4UPolyhedra.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/G4UPolyhedra.cc (Version 11.3.0) and /geometry/solids/specific/src/G4UPolyhedra.cc (Version 10.1.p2)


  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 G4UPolyhedra wrapper clas << 
 27 //                                                 26 //
 28 // 31.10.13 G.Cosmo, CERN                      <<  27 // $Id:$
                                                   >>  28 //
                                                   >>  29 // Implementation of G4UPolycone wrapper class
 29 // -------------------------------------------     30 // --------------------------------------------------------------------
 30                                                    31 
 31 #include "G4Polyhedra.hh"                          32 #include "G4Polyhedra.hh"
 32 #include "G4UPolyhedra.hh"                         33 #include "G4UPolyhedra.hh"
 33                                                << 
 34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G << 
 35                                                << 
 36 #include "G4GeomTools.hh"                      << 
 37 #include "G4GeometryTolerance.hh"              << 
 38 #include "G4AffineTransform.hh"                << 
 39 #include "G4VPVParameterisation.hh"                34 #include "G4VPVParameterisation.hh"
 40 #include "G4BoundingEnvelope.hh"               << 
 41                                                    35 
 42 using namespace CLHEP;                         <<  36 using CLHEP::twopi;
 43                                                    37 
 44 //////////////////////////////////////////////     38 ////////////////////////////////////////////////////////////////////////
 45 //                                                 39 //
 46 // Constructor (GEANT3 style parameters)           40 // Constructor (GEANT3 style parameters)
 47 //                                                 41 //
 48 // GEANT3 PGON radii are specified in the dist     42 // GEANT3 PGON radii are specified in the distance to the norm of each face.
 49 //                                             <<  43 //  
 50 G4UPolyhedra::G4UPolyhedra(const G4String& nam <<  44 G4UPolyhedra::G4UPolyhedra(const G4String& name, 
 51                                  G4double phiS     45                                  G4double phiStart,
 52                                  G4double phiT     46                                  G4double phiTotal,
 53                                  G4int numSide <<  47                                  G4int numSide,  
 54                                  G4int numZPla     48                                  G4int numZPlanes,
 55                            const G4double zPla     49                            const G4double zPlane[],
 56                            const G4double rInn     50                            const G4double rInner[],
 57                            const G4double rOut     51                            const G4double rOuter[]  )
 58   : Base_t(name, phiStart, phiTotal, numSide,  <<  52   : G4USolid(name, new UPolyhedra(name,phiStart, phiTotal, numSide,
 59            numZPlanes, zPlane, rInner, rOuter) <<  53                                   numZPlanes, zPlane, rInner, rOuter))
 60 {                                                  54 {
 61   fGenericPgon = false;                        << 
 62   SetOriginalParameters();                     << 
 63   wrStart = phiStart;                          << 
 64   while (wrStart < 0)                          << 
 65   {                                            << 
 66     wrStart += twopi;                          << 
 67   }                                            << 
 68   wrDelta = phiTotal;                          << 
 69   if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL << 
 70   {                                            << 
 71     wrDelta = twopi;                           << 
 72   }                                            << 
 73   wrNumSide = numSide;                         << 
 74   G4double convertRad = 1./std::cos(0.5*wrDelt << 
 75   rzcorners.resize(0);                         << 
 76   for (G4int i=0; i<numZPlanes; ++i)           << 
 77   {                                            << 
 78     G4double z = zPlane[i];                    << 
 79     G4double r = rOuter[i]*convertRad;         << 
 80     rzcorners.emplace_back(r,z);               << 
 81   }                                            << 
 82   for (G4int i=numZPlanes-1; i>=0; --i)        << 
 83   {                                            << 
 84     G4double z = zPlane[i];                    << 
 85     G4double r = rInner[i]*convertRad;         << 
 86     rzcorners.emplace_back(r,z);               << 
 87   }                                            << 
 88   std::vector<G4int> iout;                     << 
 89   G4GeomTools::RemoveRedundantVertices(rzcorne << 
 90 }                                                  55 }
 91                                                    56 
 92                                                    57 
 93 //////////////////////////////////////////////     58 ////////////////////////////////////////////////////////////////////////
 94 //                                                 59 //
 95 // Constructor (generic parameters)                60 // Constructor (generic parameters)
 96 //                                                 61 //
 97 G4UPolyhedra::G4UPolyhedra(const G4String& nam <<  62 G4UPolyhedra::G4UPolyhedra(const G4String& name, 
 98                                  G4double phiS     63                                  G4double phiStart,
 99                                  G4double phiT     64                                  G4double phiTotal,
100                                  G4int    numS <<  65                                  G4int    numSide,  
101                                  G4int    numR     66                                  G4int    numRZ,
102                            const G4double r[],     67                            const G4double r[],
103                            const G4double z[]      68                            const G4double z[]   )
104   : Base_t(name, phiStart, phiTotal, numSide,  <<  69   : G4USolid(name, new UPolyhedra(name, phiStart, phiTotal, numSide,
105 {                                              <<  70                                   numRZ, r, z))
106   fGenericPgon = true;                         <<  71 { 
107   SetOriginalParameters();                     << 
108   wrStart = phiStart;                          << 
109   while (wrStart < 0.)                         << 
110   {                                            << 
111     wrStart += twopi;                          << 
112   }                                            << 
113   wrDelta = phiTotal;                          << 
114   if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL << 
115   {                                            << 
116     wrDelta = twopi;                           << 
117   }                                            << 
118   wrNumSide = numSide;                         << 
119   rzcorners.resize(0);                         << 
120   for (G4int i=0; i<numRZ; ++i)                << 
121   {                                            << 
122     rzcorners.emplace_back(r[i],z[i]);         << 
123   }                                            << 
124   std::vector<G4int> iout;                     << 
125   G4GeomTools::RemoveRedundantVertices(rzcorne << 
126 }                                                  72 }
127                                                    73 
128                                                    74 
129 //////////////////////////////////////////////     75 ////////////////////////////////////////////////////////////////////////
130 //                                                 76 //
131 // Fake default constructor - sets only member     77 // Fake default constructor - sets only member data and allocates memory
132 //                            for usage restri     78 //                            for usage restricted to object persistency.
133 //                                                 79 //
134 G4UPolyhedra::G4UPolyhedra( __void__& a )          80 G4UPolyhedra::G4UPolyhedra( __void__& a )
135   : Base_t(a)                                  <<  81   : G4USolid(a)
136 {                                                  82 {
137 }                                                  83 }
138                                                    84 
139                                                    85 
140 //////////////////////////////////////////////     86 ////////////////////////////////////////////////////////////////////////
141 //                                                 87 //
142 // Destructor                                      88 // Destructor
143 //                                                 89 //
144 G4UPolyhedra::~G4UPolyhedra() = default;       <<  90 G4UPolyhedra::~G4UPolyhedra()
                                                   >>  91 {
                                                   >>  92 }
145                                                    93 
146                                                    94 
147 //////////////////////////////////////////////     95 ////////////////////////////////////////////////////////////////////////
148 //                                                 96 //
149 // Copy constructor                                97 // Copy constructor
150 //                                                 98 //
151 G4UPolyhedra::G4UPolyhedra( const G4UPolyhedra <<  99 G4UPolyhedra::G4UPolyhedra( const G4UPolyhedra &source )
152   : Base_t( source )                           << 100   : G4USolid( source )
153 {                                                 101 {
154   fGenericPgon = source.fGenericPgon;          << 
155   fOriginalParameters = source.fOriginalParame << 
156   wrStart   = source.wrStart;                  << 
157   wrDelta   = source.wrDelta;                  << 
158   wrNumSide = source.wrNumSide;                << 
159   rzcorners = source.rzcorners;                << 
160 }                                                 102 }
161                                                   103 
162                                                   104 
163 //////////////////////////////////////////////    105 ////////////////////////////////////////////////////////////////////////
164 //                                                106 //
165 // Assignment operator                            107 // Assignment operator
166 //                                                108 //
167 G4UPolyhedra& G4UPolyhedra::operator=( const G << 109 G4UPolyhedra& G4UPolyhedra::operator=( const G4UPolyhedra &source )
168 {                                                 110 {
169   if (this == &source) return *this;              111   if (this == &source) return *this;
170                                                   112 
171   Base_t::operator=( source );                 << 113   G4USolid::operator=( source );
172   fGenericPgon = source.fGenericPgon;          << 114     
173   fOriginalParameters = source.fOriginalParame << 
174   wrStart   = source.wrStart;                  << 
175   wrDelta   = source.wrDelta;                  << 
176   wrNumSide = source.wrNumSide;                << 
177   rzcorners = source.rzcorners;                << 
178                                                << 
179   return *this;                                   115   return *this;
180 }                                                 116 }
181                                                   117 
182                                                   118 
183 //////////////////////////////////////////////    119 ////////////////////////////////////////////////////////////////////////
184 //                                                120 //
185 // Accessors & modifiers                       << 
186 //                                             << 
187 G4int G4UPolyhedra::GetNumSide() const         << 
188 {                                              << 
189   return wrNumSide;                            << 
190 }                                              << 
191 G4double G4UPolyhedra::GetStartPhi() const     << 
192 {                                              << 
193   return wrStart;                              << 
194 }                                              << 
195 G4double G4UPolyhedra::GetEndPhi() const       << 
196 {                                              << 
197   return (wrStart + wrDelta);                  << 
198 }                                              << 
199 G4double G4UPolyhedra::GetSinStartPhi() const  << 
200 {                                              << 
201   G4double phi = GetStartPhi();                << 
202   return std::sin(phi);                        << 
203 }                                              << 
204 G4double G4UPolyhedra::GetCosStartPhi() const  << 
205 {                                              << 
206   G4double phi = GetStartPhi();                << 
207   return std::cos(phi);                        << 
208 }                                              << 
209 G4double G4UPolyhedra::GetSinEndPhi() const    << 
210 {                                              << 
211   G4double phi = GetEndPhi();                  << 
212   return std::sin(phi);                        << 
213 }                                              << 
214 G4double G4UPolyhedra::GetCosEndPhi() const    << 
215 {                                              << 
216   G4double phi = GetEndPhi();                  << 
217   return std::cos(phi);                        << 
218 }                                              << 
219 G4bool G4UPolyhedra::IsOpen() const            << 
220 {                                              << 
221   return (wrDelta < twopi);                    << 
222 }                                              << 
223 G4bool G4UPolyhedra::IsGeneric() const         << 
224 {                                              << 
225   return fGenericPgon;                         << 
226 }                                              << 
227 G4int G4UPolyhedra::GetNumRZCorner() const     << 
228 {                                              << 
229   return rzcorners.size();                     << 
230 }                                              << 
231 G4PolyhedraSideRZ G4UPolyhedra::GetCorner(G4in << 
232 {                                              << 
233   G4TwoVector rz = rzcorners.at(index);        << 
234   G4PolyhedraSideRZ psiderz = { rz.x(), rz.y() << 
235                                                << 
236   return psiderz;                              << 
237 }                                              << 
238 G4PolyhedraHistorical* G4UPolyhedra::GetOrigin << 
239 {                                              << 
240   return new G4PolyhedraHistorical(fOriginalPa << 
241 }                                              << 
242 void G4UPolyhedra::SetOriginalParameters()     << 
243 {                                              << 
244   G4double startPhi = GetPhiStart();           << 
245   G4double deltaPhi = GetPhiDelta();           << 
246   G4int numPlanes   = GetZSegmentCount() + 1;  << 
247   G4int numSides    = GetSideCount();          << 
248                                                << 
249   fOriginalParameters.Start_angle   = startPhi << 
250   fOriginalParameters.Opening_angle = deltaPhi << 
251   fOriginalParameters.Num_z_planes  = numPlane << 
252   fOriginalParameters.numSide       = numSides << 
253                                                << 
254   delete [] fOriginalParameters.Z_values;      << 
255   delete [] fOriginalParameters.Rmin;          << 
256   delete [] fOriginalParameters.Rmax;          << 
257   fOriginalParameters.Z_values = new G4double[ << 
258   fOriginalParameters.Rmin = new G4double[numP << 
259   fOriginalParameters.Rmax = new G4double[numP << 
260                                                << 
261   G4double convertRad = fGenericPgon           << 
262                       ? 1.0 : std::cos(0.5*del << 
263   for (G4int i=0; i<numPlanes; ++i)            << 
264   {                                            << 
265     fOriginalParameters.Z_values[i] = GetZPlan << 
266     fOriginalParameters.Rmax[i]     = GetRMax( << 
267     fOriginalParameters.Rmin[i]     = GetRMin( << 
268   }                                            << 
269 }                                              << 
270 void G4UPolyhedra::SetOriginalParameters(G4Pol << 
271 {                                              << 
272   fOriginalParameters = *pars;                 << 
273   fRebuildPolyhedron = true;                   << 
274   Reset();                                     << 
275 }                                              << 
276                                                << 
277 G4bool G4UPolyhedra::Reset()                   << 
278 {                                              << 
279   if (fGenericPgon)                            << 
280   {                                            << 
281     std::ostringstream message;                << 
282     message << "Solid " << GetName() << " buil << 
283             << G4endl << "Not applicable to th << 
284     G4Exception("G4UPolyhedra::Reset()", "Geom << 
285                 JustWarning, message, "Paramet << 
286     return true;  // error code set            << 
287   }                                            << 
288                                                << 
289   //                                           << 
290   // Rebuild polyhedra based on original param << 
291   //                                           << 
292   wrStart = fOriginalParameters.Start_angle;   << 
293   while (wrStart < 0.)                         << 
294   {                                            << 
295     wrStart += twopi;                          << 
296   }                                            << 
297   wrDelta = fOriginalParameters.Opening_angle; << 
298   if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL << 
299   {                                            << 
300     wrDelta = twopi;                           << 
301   }                                            << 
302   wrNumSide = fOriginalParameters.numSide;     << 
303   rzcorners.resize(0);                         << 
304   for (G4int i=0; i<fOriginalParameters.Num_z_ << 
305   {                                            << 
306     G4double z = fOriginalParameters.Z_values[ << 
307     G4double r = fOriginalParameters.Rmax[i];  << 
308     rzcorners.emplace_back(r,z);               << 
309   }                                            << 
310   for (G4int i=fOriginalParameters.Num_z_plane << 
311   {                                            << 
312     G4double z = fOriginalParameters.Z_values[ << 
313     G4double r = fOriginalParameters.Rmin[i];  << 
314     rzcorners.emplace_back(r,z);               << 
315   }                                            << 
316   std::vector<G4int> iout;                     << 
317   G4GeomTools::RemoveRedundantVertices(rzcorne << 
318                                                << 
319   return false;  // error code unset           << 
320 }                                              << 
321                                                << 
322                                                << 
323 ////////////////////////////////////////////// << 
324 //                                             << 
325 // Dispatch to parameterisation for replicatio    121 // Dispatch to parameterisation for replication mechanism dimension
326 // computation & modification.                    122 // computation & modification.
327 //                                                123 //
328 void G4UPolyhedra::ComputeDimensions(G4VPVPara    124 void G4UPolyhedra::ComputeDimensions(G4VPVParameterisation* p,
329                                      const G4i    125                                      const G4int n,
330                                      const G4V    126                                      const G4VPhysicalVolume* pRep)
331 {                                                 127 {
332   p->ComputeDimensions(*(G4Polyhedra*)this,n,p    128   p->ComputeDimensions(*(G4Polyhedra*)this,n,pRep);
333 }                                                 129 }
334                                                   130 
335                                                   131 
336 //////////////////////////////////////////////    132 //////////////////////////////////////////////////////////////////////////
337 //                                                133 //
338 // Make a clone of the object                     134 // Make a clone of the object
339                                                   135 
340 G4VSolid* G4UPolyhedra::Clone() const             136 G4VSolid* G4UPolyhedra::Clone() const
341 {                                                 137 {
342   return new G4UPolyhedra(*this);                 138   return new G4UPolyhedra(*this);
343 }                                                 139 }
344                                                   140 
345                                                   141 
346 ////////////////////////////////////////////// << 142 ////////////////////////////////////////////////////////////////////////
347 //                                                143 //
348 // Get bounding box                            << 144 // CreatePolyhedron
349                                                << 145 //
350 void G4UPolyhedra::BoundingLimits(G4ThreeVecto << 146 G4Polyhedron* G4UPolyhedra::CreatePolyhedron() const
351                                   G4ThreeVecto << 
352 {                                                 147 {
353   static G4bool checkBBox = true;              << 148   if (!IsGeneric())
354   static G4bool checkPhi  = true;              << 
355                                                << 
356   G4double rmin = kInfinity, rmax = -kInfinity << 
357   G4double zmin = kInfinity, zmax = -kInfinity << 
358   for (G4int i=0; i<GetNumRZCorner(); ++i)     << 
359   {                                               149   {
360     G4PolyhedraSideRZ corner = GetCorner(i);   << 150     G4PolyhedraHistorical* original_parameters = GetOriginalParameters();
361     if (corner.r < rmin) rmin = corner.r;      << 151     G4PolyhedronPgon*
362     if (corner.r > rmax) rmax = corner.r;      << 152       polyhedron = new G4PolyhedronPgon( GetOriginalParameters()->Start_angle,
363     if (corner.z < zmin) zmin = corner.z;      << 153                                  GetOriginalParameters()->Opening_angle,
364     if (corner.z > zmax) zmax = corner.z;      << 154                                  GetOriginalParameters()->numSide,
365   }                                            << 155                                  GetOriginalParameters()->Num_z_planes,
366                                                << 156                                  GetOriginalParameters()->Z_values,
367   G4double sphi    = GetStartPhi();            << 157                                  GetOriginalParameters()->Rmin,
368   G4double ephi    = GetEndPhi();              << 158                                  GetOriginalParameters()->Rmax);
369   G4double dphi    = IsOpen() ? ephi-sphi : tw << 159     delete original_parameters;  // delete local copy 
370   G4int    ksteps  = GetNumSide();             << 160     return polyhedron;
371   G4double astep   = dphi/ksteps;              << 161   }
372   G4double sinStep = std::sin(astep);          << 162   else
373   G4double cosStep = std::cos(astep);          << 163   {
374                                                << 164     // The following code prepares for:
375   G4double sinCur = GetSinStartPhi();          << 165     // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
376   G4double cosCur = GetCosStartPhi();          << 166     //                                 const double xyz[][3],
377   if (!IsOpen()) rmin = 0.;                    << 167     //                                 const int faces_vec[][4])
378   G4double xmin = rmin*cosCur, xmax = xmin;    << 168     // Here is an extract from the header file HepPolyhedron.h:
379   G4double ymin = rmin*sinCur, ymax = ymin;    << 169     /**
380   for (G4int k=0; k<ksteps+1; ++k)             << 170      * Creates user defined polyhedron.
381   {                                            << 171      * This function allows to the user to define arbitrary polyhedron.
382     G4double x = rmax*cosCur;                  << 172      * The faces of the polyhedron should be either triangles or planar
383     if (x < xmin) xmin = x;                    << 173      * quadrilateral. Nodes of a face are defined by indexes pointing to
384     if (x > xmax) xmax = x;                    << 174      * the elements in the xyz array. Numeration of the elements in the
385     G4double y = rmax*sinCur;                  << 175      * array starts from 1 (like in fortran). The indexes can be positive
386     if (y < ymin) ymin = y;                    << 176      * or negative. Negative sign means that the corresponding edge is
387     if (y > ymax) ymax = y;                    << 177      * invisible. The normal of the face should be directed to exterior
388     if (rmin > 0.)                             << 178      * of the polyhedron. 
                                                   >> 179      * 
                                                   >> 180      * @param  Nnodes number of nodes
                                                   >> 181      * @param  Nfaces number of faces
                                                   >> 182      * @param  xyz    nodes
                                                   >> 183      * @param  faces_vec  faces (quadrilaterals or triangles)
                                                   >> 184      * @return status of the operation - is non-zero in case of problem
                                                   >> 185      */
                                                   >> 186     G4int nNodes;
                                                   >> 187     G4int nFaces;
                                                   >> 188     typedef G4double double3[3];
                                                   >> 189     double3* xyz;
                                                   >> 190     typedef G4int int4[4];
                                                   >> 191     int4* faces_vec;
                                                   >> 192     if (IsOpen())
389     {                                             193     {
390       G4double xx = rmin*cosCur;               << 194       // Triangulate open ends.  Simple ear-chopping algorithm...
391       if (xx < xmin) xmin = xx;                << 195       // I'm not sure how robust this algorithm is (J.Allison).
392       if (xx > xmax) xmax = xx;                << 196       //
393       G4double yy = rmin*sinCur;               << 197       std::vector<G4bool> chopped(GetNumRZCorner(), false);
394       if (yy < ymin) ymin = yy;                << 198       std::vector<G4int*> triQuads;
395       if (yy > ymax) ymax = yy;                << 199       G4int remaining = GetNumRZCorner();
                                                   >> 200       G4int iStarter = 0;
                                                   >> 201       while (remaining >= 3)
                                                   >> 202       {
                                                   >> 203         // Find unchopped corners...
                                                   >> 204         //
                                                   >> 205         G4int A = -1, B = -1, C = -1;
                                                   >> 206         G4int iStepper = iStarter;
                                                   >> 207         do
                                                   >> 208         {
                                                   >> 209           if (A < 0)      { A = iStepper; }
                                                   >> 210           else if (B < 0) { B = iStepper; }
                                                   >> 211           else if (C < 0) { C = iStepper; }
                                                   >> 212           do
                                                   >> 213           {
                                                   >> 214             if (++iStepper >= GetNumRZCorner()) iStepper = 0;
                                                   >> 215           }
                                                   >> 216           while (chopped[iStepper]);
                                                   >> 217         }
                                                   >> 218         while (C < 0 && iStepper != iStarter);
                                                   >> 219 
                                                   >> 220         // Check triangle at B is pointing outward (an "ear").
                                                   >> 221         // Sign of z cross product determines...
                                                   >> 222 
                                                   >> 223         G4double BAr = GetCorner(A).r - GetCorner(B).r;
                                                   >> 224         G4double BAz = GetCorner(A).z - GetCorner(B).z;
                                                   >> 225         G4double BCr = GetCorner(C).r - GetCorner(B).r;
                                                   >> 226         G4double BCz = GetCorner(C).z - GetCorner(B).z;
                                                   >> 227         if (BAr * BCz - BAz * BCr < kCarTolerance)
                                                   >> 228         {
                                                   >> 229           G4int* tq = new G4int[3];
                                                   >> 230           tq[0] = A + 1;
                                                   >> 231           tq[1] = B + 1;
                                                   >> 232           tq[2] = C + 1;
                                                   >> 233           triQuads.push_back(tq);
                                                   >> 234           chopped[B] = true;
                                                   >> 235           --remaining;
                                                   >> 236         }
                                                   >> 237         else
                                                   >> 238         {
                                                   >> 239           do
                                                   >> 240           {
                                                   >> 241             if (++iStarter >= GetNumRZCorner()) { iStarter = 0; }
                                                   >> 242           }
                                                   >> 243           while (chopped[iStarter]);
                                                   >> 244         }
                                                   >> 245       }
                                                   >> 246 
                                                   >> 247       // Transfer to faces...
                                                   >> 248       G4int numSide=GetNumSide();
                                                   >> 249       nNodes = (numSide + 1) * GetNumRZCorner();
                                                   >> 250       nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size();
                                                   >> 251       faces_vec = new int4[nFaces];
                                                   >> 252       G4int iface = 0;
                                                   >> 253       G4int addition = GetNumRZCorner() * numSide;
                                                   >> 254       G4int d = GetNumRZCorner() - 1;
                                                   >> 255       for (G4int iEnd = 0; iEnd < 2; ++iEnd)
                                                   >> 256       {
                                                   >> 257         for (size_t i = 0; i < triQuads.size(); ++i)
                                                   >> 258         {
                                                   >> 259           // Negative for soft/auxiliary/normally invisible edges...
                                                   >> 260           //
                                                   >> 261           G4int a, b, c;
                                                   >> 262           if (iEnd == 0)
                                                   >> 263           {
                                                   >> 264             a = triQuads[i][0];
                                                   >> 265             b = triQuads[i][1];
                                                   >> 266             c = triQuads[i][2];
                                                   >> 267           }
                                                   >> 268           else
                                                   >> 269           {
                                                   >> 270             a = triQuads[i][0] + addition;
                                                   >> 271             b = triQuads[i][2] + addition;
                                                   >> 272             c = triQuads[i][1] + addition;
                                                   >> 273           }
                                                   >> 274           G4int ab = std::abs(b - a);
                                                   >> 275           G4int bc = std::abs(c - b);
                                                   >> 276           G4int ca = std::abs(a - c);
                                                   >> 277           faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
                                                   >> 278           faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
                                                   >> 279           faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
                                                   >> 280           faces_vec[iface][3] = 0;
                                                   >> 281           ++iface;
                                                   >> 282         }
                                                   >> 283       }
                                                   >> 284 
                                                   >> 285       // Continue with sides...
                                                   >> 286 
                                                   >> 287       xyz = new double3[nNodes];
                                                   >> 288       const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide;
                                                   >> 289       G4double phi = GetStartPhi();
                                                   >> 290       G4int ixyz = 0;
                                                   >> 291       for (G4int iSide = 0; iSide < numSide; ++iSide)
                                                   >> 292       {
                                                   >> 293         for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
                                                   >> 294         {
                                                   >> 295           xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
                                                   >> 296           xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
                                                   >> 297           xyz[ixyz][2] = GetCorner(iCorner).z;
                                                   >> 298           if (iCorner < GetNumRZCorner() - 1)
                                                   >> 299           {
                                                   >> 300             faces_vec[iface][0] = ixyz + 1;
                                                   >> 301             faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
                                                   >> 302             faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
                                                   >> 303             faces_vec[iface][3] = ixyz + 2;
                                                   >> 304           }
                                                   >> 305           else
                                                   >> 306           {
                                                   >> 307             faces_vec[iface][0] = ixyz + 1;
                                                   >> 308             faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
                                                   >> 309             faces_vec[iface][2] = ixyz + 2;
                                                   >> 310             faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
                                                   >> 311           }
                                                   >> 312           ++iface;
                                                   >> 313           ++ixyz;
                                                   >> 314         }
                                                   >> 315         phi += dPhi;
                                                   >> 316       }
                                                   >> 317 
                                                   >> 318       // Last GetCorner...
                                                   >> 319 
                                                   >> 320       for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
                                                   >> 321       {
                                                   >> 322         xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
                                                   >> 323         xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
                                                   >> 324         xyz[ixyz][2] = GetCorner(iCorner).z;
                                                   >> 325         ++ixyz;
                                                   >> 326       }
396     }                                             327     }
397     G4double sinTmp = sinCur;                  << 328     else  // !phiIsOpen - i.e., a complete 360 degrees.
398     sinCur = sinCur*cosStep + cosCur*sinStep;  << 
399     cosCur = cosCur*cosStep - sinTmp*sinStep;  << 
400   }                                            << 
401   pMin.set(xmin,ymin,zmin);                    << 
402   pMax.set(xmax,ymax,zmax);                    << 
403                                                << 
404   // Check correctness of the bounding box     << 
405   //                                           << 
406   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
407   {                                            << 
408     std::ostringstream message;                << 
409     message << "Bad bounding box (min >= max)  << 
410             << GetName() << " !"               << 
411             << "\npMin = " << pMin             << 
412             << "\npMax = " << pMax;            << 
413     G4Exception("G4UPolyhedra::BoundingLimits( << 
414                 JustWarning, message);         << 
415     StreamInfo(G4cout);                        << 
416   }                                            << 
417                                                << 
418   // Check consistency of bounding boxes       << 
419   //                                           << 
420   if (checkBBox)                               << 
421   {                                            << 
422     U3Vector vmin, vmax;                       << 
423     Extent(vmin,vmax);                         << 
424     if (std::abs(pMin.x()-vmin.x()) > kCarTole << 
425         std::abs(pMin.y()-vmin.y()) > kCarTole << 
426         std::abs(pMin.z()-vmin.z()) > kCarTole << 
427         std::abs(pMax.x()-vmax.x()) > kCarTole << 
428         std::abs(pMax.y()-vmax.y()) > kCarTole << 
429         std::abs(pMax.z()-vmax.z()) > kCarTole << 
430     {                                             329     {
431       std::ostringstream message;              << 330       nNodes = GetNumSide() * GetNumRZCorner();
432       message << "Inconsistency in bounding bo << 331       nFaces = GetNumSide() * GetNumRZCorner();;
433               << GetName() << " !"             << 332       xyz = new double3[nNodes];
434               << "\nBBox min: wrapper = " << p << 333       faces_vec = new int4[nFaces];
435               << "\nBBox max: wrapper = " << p << 334       // const G4double dPhi = (endPhi - startPhi) / numSide;
436       G4Exception("G4UPolyhedra::BoundingLimit << 335       const G4double dPhi = twopi / GetNumSide(); // !phiIsOpen endPhi-startPhi = 360 degrees.
437                   JustWarning, message);       << 336       G4double phi = GetStartPhi();
438       checkBBox = false;                       << 337       G4int ixyz = 0, iface = 0;
                                                   >> 338       for (G4int iSide = 0; iSide < GetNumSide(); ++iSide)
                                                   >> 339       {
                                                   >> 340         for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner)
                                                   >> 341         {
                                                   >> 342           xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi);
                                                   >> 343           xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi);
                                                   >> 344           xyz[ixyz][2] = GetCorner(iCorner).z;
                                                   >> 345           if (iSide < GetNumSide() - 1)
                                                   >> 346           {
                                                   >> 347             if (iCorner < GetNumRZCorner() - 1)
                                                   >> 348             {
                                                   >> 349               faces_vec[iface][0] = ixyz + 1;
                                                   >> 350               faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
                                                   >> 351               faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2;
                                                   >> 352               faces_vec[iface][3] = ixyz + 2;
                                                   >> 353             }
                                                   >> 354             else
                                                   >> 355             {
                                                   >> 356               faces_vec[iface][0] = ixyz + 1;
                                                   >> 357               faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1;
                                                   >> 358               faces_vec[iface][2] = ixyz + 2;
                                                   >> 359               faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
                                                   >> 360             }
                                                   >> 361           }
                                                   >> 362           else   // Last side joins ends...
                                                   >> 363           {
                                                   >> 364             if (iCorner < GetNumRZCorner() - 1)
                                                   >> 365             {
                                                   >> 366               faces_vec[iface][0] = ixyz + 1;
                                                   >> 367               faces_vec[iface][1] = ixyz + GetNumRZCorner() - nFaces + 1;
                                                   >> 368               faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2;
                                                   >> 369               faces_vec[iface][3] = ixyz + 2;
                                                   >> 370             }
                                                   >> 371             else
                                                   >> 372             {
                                                   >> 373               faces_vec[iface][0] = ixyz + 1;
                                                   >> 374               faces_vec[iface][1] = ixyz - nFaces + GetNumRZCorner() + 1;
                                                   >> 375               faces_vec[iface][2] = ixyz - nFaces + 2;
                                                   >> 376               faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2;
                                                   >> 377             }
                                                   >> 378           }
                                                   >> 379           ++ixyz;
                                                   >> 380           ++iface;
                                                   >> 381         }
                                                   >> 382         phi += dPhi;
                                                   >> 383       }
439     }                                             384     }
440   }                                            << 385     G4Polyhedron* polyhedron = new G4Polyhedron;
441                                                << 386     G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
442   // Check consistency of angles               << 387     delete [] faces_vec;
443   //                                           << 388     delete [] xyz;
444   if (checkPhi)                                << 389     if (problem)
445   {                                            << 
446     if (GetStartPhi() != GetPhiStart()  ||     << 
447         GetEndPhi()   != GetPhiEnd()    ||     << 
448         GetNumSide()  != GetSideCount() ||     << 
449         IsOpen()      != (Base_t::GetPhiDelta( << 
450     {                                             390     {
451       std::ostringstream message;                 391       std::ostringstream message;
452       message << "Inconsistency in Phi angles  << 392       message << "Problem creating G4Polyhedron for: " << GetName();
453               << GetName() << " !"             << 393       G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
454               << "\nPhi start  : wrapper = " < << 
455               << " solid = " <<     GetPhiStar << 
456               << "\nPhi end    : wrapper = " < << 
457               << " solid = " <<     GetPhiEnd( << 
458               << "\nPhi # sides: wrapper = " < << 
459               << " solid = " <<     GetSideCou << 
460               << "\nPhi is open: wrapper = " < << 
461               << " solid = "                   << 
462               << ((Base_t::GetPhiDelta() < two << 
463       G4Exception("G4UPolyhedra::BoundingLimit << 
464                   JustWarning, message);          394                   JustWarning, message);
465       checkPhi = false;                        << 395       delete polyhedron;
                                                   >> 396       return 0;
466     }                                             397     }
467   }                                            << 398     else
468 }                                              << 
469                                                << 
470 ////////////////////////////////////////////// << 
471 //                                             << 
472 // Calculate extent under transform and specif << 
473                                                << 
474 G4bool                                         << 
475 G4UPolyhedra::CalculateExtent(const EAxis pAxi << 
476                               const G4VoxelLim << 
477                               const G4AffineTr << 
478                                     G4double&  << 
479 {                                              << 
480   G4ThreeVector bmin, bmax;                    << 
481   G4bool exist;                                << 
482                                                << 
483   // Check bounding box (bbox)                 << 
484   //                                           << 
485   BoundingLimits(bmin,bmax);                   << 
486   G4BoundingEnvelope bbox(bmin,bmax);          << 
487 #ifdef G4BBOX_EXTENT                           << 
488   if (true) return bbox.CalculateExtent(pAxis, << 
489 #endif                                         << 
490   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 
491   {                                            << 
492     return exist = pMin < pMax;                << 
493   }                                            << 
494                                                << 
495   // To find the extent, RZ contour of the pol << 
496   // in triangles. The extent is calculated as << 
497   // all sub-polycones formed by rotation of t << 
498   //                                           << 
499   G4TwoVectorList contourRZ;                   << 
500   G4TwoVectorList triangles;                   << 
501   std::vector<G4int> iout;                     << 
502   G4double eminlim = pVoxelLimit.GetMinExtent( << 
503   G4double emaxlim = pVoxelLimit.GetMaxExtent( << 
504                                                << 
505   // get RZ contour, ensure anticlockwise orde << 
506   for (G4int i=0; i<GetNumRZCorner(); ++i)     << 
507   {                                            << 
508     G4PolyhedraSideRZ corner = GetCorner(i);   << 
509     contourRZ.emplace_back(corner.r,corner.z); << 
510   }                                            << 
511   G4GeomTools::RemoveRedundantVertices(contour << 
512   G4double area = G4GeomTools::PolygonArea(con << 
513   if (area < 0.) std::reverse(contourRZ.begin( << 
514                                                << 
515   // triangulate RZ countour                   << 
516   if (!G4GeomTools::TriangulatePolygon(contour << 
517   {                                            << 
518     std::ostringstream message;                << 
519     message << "Triangulation of RZ contour ha << 
520             << GetName() << " !"               << 
521             << "\nExtent has been calculated u << 
522     G4Exception("G4UPolyhedra::CalculateExtent << 
523                 "GeomMgt1002",JustWarning,mess << 
524     return bbox.CalculateExtent(pAxis,pVoxelLi << 
525   }                                            << 
526                                                << 
527   // set trigonometric values                  << 
528   G4double sphi     = GetStartPhi();           << 
529   G4double ephi     = GetEndPhi();             << 
530   G4double dphi     = IsOpen() ? ephi-sphi : t << 
531   G4int    ksteps   = GetNumSide();            << 
532   G4double astep    = dphi/ksteps;             << 
533   G4double sinStep  = std::sin(astep);         << 
534   G4double cosStep  = std::cos(astep);         << 
535   G4double sinStart = GetSinStartPhi();        << 
536   G4double cosStart = GetCosStartPhi();        << 
537                                                << 
538   // allocate vector lists                     << 
539   std::vector<const G4ThreeVectorList *> polyg << 
540   polygons.resize(ksteps+1);                   << 
541   for (G4int k=0; k<ksteps+1; ++k)             << 
542   {                                            << 
543     polygons[k] = new G4ThreeVectorList(3);    << 
544   }                                            << 
545                                                << 
546   // main loop along triangles                 << 
547   pMin =  kInfinity;                           << 
548   pMax = -kInfinity;                           << 
549   G4int ntria = triangles.size()/3;            << 
550   for (G4int i=0; i<ntria; ++i)                << 
551   {                                            << 
552     G4double sinCur = sinStart;                << 
553     G4double cosCur = cosStart;                << 
554     G4int i3 = i*3;                            << 
555     for (G4int k=0; k<ksteps+1; ++k) // rotate << 
556     {                                             399     {
557       auto ptr = const_cast<G4ThreeVectorList* << 400       return polyhedron;
558       auto iter = ptr->begin();                << 
559       iter->set(triangles[i3+0].x()*cosCur,    << 
560                 triangles[i3+0].x()*sinCur,    << 
561                 triangles[i3+0].y());          << 
562       iter++;                                  << 
563       iter->set(triangles[i3+1].x()*cosCur,    << 
564                 triangles[i3+1].x()*sinCur,    << 
565                 triangles[i3+1].y());          << 
566       iter++;                                  << 
567       iter->set(triangles[i3+2].x()*cosCur,    << 
568                 triangles[i3+2].x()*sinCur,    << 
569                 triangles[i3+2].y());          << 
570                                                << 
571       G4double sinTmp = sinCur;                << 
572       sinCur = sinCur*cosStep + cosCur*sinStep << 
573       cosCur = cosCur*cosStep - sinTmp*sinStep << 
574     }                                             401     }
575                                                << 
576     // set sub-envelope and adjust extent      << 
577     G4double emin,emax;                        << 
578     G4BoundingEnvelope benv(polygons);         << 
579     if (!benv.CalculateExtent(pAxis,pVoxelLimi << 
580     if (emin < pMin) pMin = emin;              << 
581     if (emax > pMax) pMax = emax;              << 
582     if (eminlim > pMin && emaxlim < pMax) brea << 
583   }                                               402   }
584   // free memory                               << 
585   for (G4int k=0; k<ksteps+1; ++k) { delete po << 
586   return (pMin < pMax);                        << 
587 }                                                 403 }
588                                                << 
589                                                << 
590 ////////////////////////////////////////////// << 
591 //                                             << 
592 // CreatePolyhedron                            << 
593 //                                             << 
594 G4Polyhedron* G4UPolyhedra::CreatePolyhedron() << 
595 {                                              << 
596   return new G4PolyhedronPgon(wrStart, wrDelta << 
597 }                                              << 
598                                                << 
599 #endif  // G4GEOM_USE_USOLIDS                  << 
600                                                   404