Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/graphics_reps/src/HepPolyhedron.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 /graphics_reps/src/HepPolyhedron.cc (Version 11.3.0) and /graphics_reps/src/HepPolyhedron.cc (Version 3.0)


                                                   >>   1 // This code implementation is the intellectual property of
                                                   >>   2 // the GEANT4 collaboration.
  1 //                                                  3 //
  2 // ******************************************* <<   4 // By copying, distributing or modifying the Program (or any work
  3 // * License and Disclaimer                    <<   5 // based on the Program) you indicate your acceptance of this statement,
  4 // *                                           <<   6 // and all its terms.
  5 // * The  Geant4 software  is  copyright of th <<   7 //
  6 // * the Geant4 Collaboration.  It is provided <<   8 // $Id: HepPolyhedron.cc,v 1.4 2000/04/18 10:03:29 evc Exp $
  7 // * conditions of the Geant4 Software License <<   9 // GEANT4 tag $Name: geant4-03-00 $
  8 // * LICENSE and available at  http://cern.ch/ <<  10 //
  9 // * include a list of copyright holders.      <<  11 // 
 10 // *                                           << 
 11 // * Neither the authors of this software syst << 
 12 // * institutes,nor the agencies providing fin << 
 13 // * work  make  any representation or  warran << 
 14 // * regarding  this  software system or assum << 
 15 // * use.  Please see the license in the file  << 
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                           << 
 18 // * This  code  implementation is the result  << 
 19 // * technical work of the GEANT4 collaboratio << 
 20 // * By using,  copying,  modifying or  distri << 
 21 // * any work based  on the software)  you  ag << 
 22 // * use  in  resulting  scientific  publicati << 
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // ******************************************* << 
 25 //                                                 12 //
 26 // G4 Polyhedron library                           13 // G4 Polyhedron library
 27 //                                                 14 //
 28 // History:                                        15 // History:
 29 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@ce     16 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
 30 //                                                 17 //
 31 // 30.09.96 E.Chernyaev                            18 // 30.09.96 E.Chernyaev
 32 // - added GetNextVertexIndex, GetVertex by Ya     19 // - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
 33 // - added GetNextUnitNormal, GetNextEdgeIndic <<  20 // - added GetNextUnitNormal, GetNextEdgeIndeces, GetNextEdge
 34 //                                                 21 //
 35 // 15.12.96 E.Chernyaev                            22 // 15.12.96 E.Chernyaev
 36 // - added GetNumberOfRotationSteps, RotateEdg     23 // - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
 37 // - rewritten G4PolyhedronCons;                   24 // - rewritten G4PolyhedronCons;
 38 // - added G4PolyhedronPara, ...Trap, ...Pgon,     25 // - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
 39 //                                                 26 //
 40 // 01.06.97 E.Chernyaev                            27 // 01.06.97 E.Chernyaev
 41 // - modified RotateAroundZ, added SetSideFace     28 // - modified RotateAroundZ, added SetSideFacets
 42 //                                                 29 //
 43 // 19.03.00 E.Chernyaev                            30 // 19.03.00 E.Chernyaev
 44 // - implemented boolean operations (add, subt     31 // - implemented boolean operations (add, subtract, intersect) on polyhedra;
 45 //                                                 32 //
 46 // 25.05.01 E.Chernyaev                        <<  33   
 47 // - added GetSurfaceArea() and GetVolume()    << 
 48 //                                             << 
 49 // 05.11.02 E.Chernyaev                        << 
 50 // - added createTwistedTrap() and createPolyh << 
 51 //                                             << 
 52 // 20.06.05 G.Cosmo                            << 
 53 // - added HepPolyhedronEllipsoid              << 
 54 //                                             << 
 55 // 18.07.07 T.Nikitina                         << 
 56 // - added HepPolyhedronParaboloid             << 
 57 //                                             << 
 58 // 22.02.20 E.Chernyaev                        << 
 59 // - added HepPolyhedronTet, HepPolyhedronHybe << 
 60 //                                             << 
 61 // 12.05.21 E.Chernyaev                        << 
 62 // - added TriangulatePolygon(), RotateContour << 
 63 // - added HepPolyhedronPgon, HepPolyhedronPco << 
 64 //                                             << 
 65 // 26.03.22 E.Chernyaev                        << 
 66 // - added SetVertex(), SetFacet()             << 
 67 // - added HepPolyhedronTetMesh                << 
 68 //                                             << 
 69 // 04.04.22 E.Chernyaev                        << 
 70 // - added JoinCoplanarFacets()                << 
 71 //                                             << 
 72 // 07.04.22 E.Chernyaev                        << 
 73 // - added HepPolyhedronBoxMesh                << 
 74                                                << 
 75 #include "HepPolyhedron.h"                         34 #include "HepPolyhedron.h"
 76 #include "G4PhysicalConstants.hh"              << 
 77 #include "G4Vector3D.hh"                       << 
 78                                                << 
 79 #include <cstdlib>  // Required on some compil << 
 80 #include <cmath>                               << 
 81 #include <algorithm>                           << 
 82                                                << 
 83 using CLHEP::perMillion;                       << 
 84 using CLHEP::deg;                              << 
 85 using CLHEP::pi;                               << 
 86 using CLHEP::twopi;                            << 
 87 using CLHEP::nm;                               << 
 88 const G4double spatialTolerance = 0.01*nm;     << 
 89                                                    35 
 90 /*********************************************     36 /***********************************************************************
 91  *                                                 37  *                                                                     *
 92  * Name: HepPolyhedron operator <<                 38  * Name: HepPolyhedron operator <<                   Date:    09.05.96 *
 93  * Author: E.Chernyaev (IHEP/Protvino)             39  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
 94  *                                                 40  *                                                                     *
 95  * Function: Print contents of G4 polyhedron       41  * Function: Print contents of G4 polyhedron                           *
 96  *                                                 42  *                                                                     *
 97  *********************************************     43  ***********************************************************************/
 98 std::ostream & operator<<(std::ostream & ostr, <<  44 HepStd::ostream & operator<<(HepStd::ostream & ostr, const G4Facet & facet) {
 99   for (const auto& edge : facet.edge) {        <<  45   for (int k=0; k<4; k++) {
100     ostr << " " << edge.v << "/" << edge.f;    <<  46     ostr << " " << facet.edge[k].v << "/" << facet.edge[k].f;
101   }                                                47   }
102   return ostr;                                     48   return ostr;
103 }                                                  49 }
104                                                    50 
105 std::ostream & operator<<(std::ostream & ostr, <<  51 HepStd::ostream & operator<<(HepStd::ostream & ostr, const HepPolyhedron & ph) {
106   ostr << std::endl;                           <<  52   ostr << HepStd::endl;
107   ostr << "Nvertices=" << ph.nvert << ", Nface <<  53   ostr << "Nverteces=" << ph.nvert << ", Nfacets=" << ph.nface << HepStd::endl;
108   G4int i;                                     <<  54   int i;
109   for (i=1; i<=ph.nvert; i++) {                    55   for (i=1; i<=ph.nvert; i++) {
110      ostr << "xyz(" << i << ")="                   56      ostr << "xyz(" << i << ")="
111           << ph.pV[i].x() << ' ' << ph.pV[i].y <<  57     << ph.pV[i].x() << ' ' << ph.pV[i].y() << ' ' << ph.pV[i].z()
112           << std::endl;                        <<  58     << HepStd::endl;
113   }                                                59   }
114   for (i=1; i<=ph.nface; i++) {                    60   for (i=1; i<=ph.nface; i++) {
115     ostr << "face(" << i << ")=" << ph.pF[i] < <<  61     ostr << "face(" << i << ")=" << ph.pF[i] << HepStd::endl;
116   }                                                62   }
117   return ostr;                                     63   return ostr;
118 }                                                  64 }
119                                                    65 
120 HepPolyhedron::HepPolyhedron(G4int Nvert, G4in << 
121 /********************************************* << 
122  *                                             << 
123  * Name: HepPolyhedron constructor with        << 
124  *       allocation of memory                  << 
125  * Author: E.Tcherniaev (E.Chernyaev)          << 
126  *                                             << 
127  ********************************************* << 
128 : nvert(0), nface(0), pV(nullptr), pF(nullptr) << 
129 {                                              << 
130   AllocateMemory(Nvert, Nface);                << 
131 }                                              << 
132                                                << 
133 HepPolyhedron::HepPolyhedron(const HepPolyhedr     66 HepPolyhedron::HepPolyhedron(const HepPolyhedron &from)
134 /*********************************************     67 /***********************************************************************
135  *                                                 68  *                                                                     *
136  * Name: HepPolyhedron copy constructor            69  * Name: HepPolyhedron copy constructor             Date:    23.07.96  *
137  * Author: E.Chernyaev (IHEP/Protvino)             70  * Author: E.Chernyaev (IHEP/Protvino)              Revised:           *
138  *                                                 71  *                                                                     *
139  *********************************************     72  ***********************************************************************/
140 : nvert(0), nface(0), pV(nullptr), pF(nullptr) << 
141 {                                                  73 {
142   AllocateMemory(from.nvert, from.nface);      <<  74   if (from.nvert > 0 && from.nface > 0) {
143   for (G4int i=1; i<=nvert; i++) pV[i] = from. <<  75     nvert = from.nvert;
144   for (G4int k=1; k<=nface; k++) pF[k] = from. <<  76     nface = from.nface;
145 }                                              <<  77     pV = new HepPoint3D[nvert + 1];
146                                                <<  78     pF = new G4Facet[nface + 1];
147 HepPolyhedron::HepPolyhedron(HepPolyhedron&& f <<  79     int i;
148 /********************************************* <<  80     for (i=1; i<=nvert; i++) pV[i] = from.pV[i];
149  *                                             <<  81     for (i=1; i<=nface; i++) pF[i] = from.pF[i];
150  * Name: HepPolyhedron move constructor        <<  82   }else{
151  * Author: E.Tcherniaev (E.Chernyaev)          <<  83     nvert = 0; nface = 0; pV = 0; pF = 0;
152  *                                             <<  84   }
153  ********************************************* << 
154 : nvert(0), nface(0), pV(nullptr), pF(nullptr) << 
155 {                                              << 
156   nvert = from.nvert;                          << 
157   nface = from.nface;                          << 
158   pV = from.pV;                                << 
159   pF = from.pF;                                << 
160                                                << 
161   // Release the data from the source object   << 
162   from.nvert = 0;                              << 
163   from.nface = 0;                              << 
164   from.pV = nullptr;                           << 
165   from.pF = nullptr;                           << 
166 }                                                  85 }
167                                                    86 
168 HepPolyhedron & HepPolyhedron::operator=(const     87 HepPolyhedron & HepPolyhedron::operator=(const HepPolyhedron &from)
169 /*********************************************     88 /***********************************************************************
170  *                                                 89  *                                                                     *
171  * Name: HepPolyhedron operator =                  90  * Name: HepPolyhedron operator =                   Date:    23.07.96  *
172  * Author: E.Chernyaev (IHEP/Protvino)             91  * Author: E.Chernyaev (IHEP/Protvino)              Revised:           *
173  *                                                 92  *                                                                     *
174  * Function: Copy contents of one polyhedron t <<  93  * Function: Copy contents of one GEANT4 polyhedron to another         *
175  *                                             << 
176  ********************************************* << 
177 {                                              << 
178   if (this != &from) {                         << 
179     AllocateMemory(from.nvert, from.nface);    << 
180     for (G4int i=1; i<=nvert; i++) pV[i] = fro << 
181     for (G4int k=1; k<=nface; k++) pF[k] = fro << 
182   }                                            << 
183   return *this;                                << 
184 }                                              << 
185                                                << 
186 HepPolyhedron & HepPolyhedron::operator=(HepPo << 
187 /********************************************* << 
188  *                                             << 
189  * Name: HepPolyhedron move operator =         << 
190  * Author: E.Tcherniaev (E.Chernyaev)          << 
191  *                                             << 
192  * Function: Move contents of one polyhedron t << 
193  *                                                 94  *                                                                     *
194  *********************************************     95  ***********************************************************************/
195 {                                                  96 {
196   if (this != &from) {                         <<  97   if (this == &from) return *this;
197     delete [] pV;                              <<  98   delete [] pV;
198     delete [] pF;                              <<  99   delete [] pF;
                                                   >> 100   if (from.nvert > 0  && from.nface > 0) {
199     nvert = from.nvert;                           101     nvert = from.nvert;
200     nface = from.nface;                           102     nface = from.nface;
201     pV = from.pV;                              << 103     pV = new HepPoint3D[nvert + 1];
202     pF = from.pF;                              << 104     pF = new G4Facet[nface + 1];
203                                                << 105     int i;
204     // Release the data from the source object << 106     for (i=1; i<=nvert; i++) pV[i] = from.pV[i];
205     from.nvert = 0;                            << 107     for (i=1; i<=nface; i++) pF[i] = from.pF[i];
206     from.nface = 0;                            << 108   }else{
207     from.pV = nullptr;                         << 109     nvert = 0; nface = 0; pV = 0; pF = 0;
208     from.pF = nullptr;                         << 
209   }                                               110   }
210   return *this;                                   111   return *this;
211 }                                                 112 }
212                                                   113 
213 G4int                                          << 114 int
214 HepPolyhedron::FindNeighbour(G4int iFace, G4in << 115 HepPolyhedron::FindNeighbour(int iFace, int iNode, int iOrder) const
215 /*********************************************    116 /***********************************************************************
216  *                                                117  *                                                                     *
217  * Name: HepPolyhedron::FindNeighbour             118  * Name: HepPolyhedron::FindNeighbour                Date:    22.11.99 *
218  * Author: E.Chernyaev                            119  * Author: E.Chernyaev                               Revised:          *
219  *                                                120  *                                                                     *
220  * Function: Find neighbouring face               121  * Function: Find neighbouring face                                    *
221  *                                                122  *                                                                     *
222  *********************************************    123  ***********************************************************************/
223 {                                                 124 {
224   G4int i;                                     << 125   int i;
225   for (i=0; i<4; i++) {                           126   for (i=0; i<4; i++) {
226     if (iNode == std::abs(pF[iFace].edge[i].v) << 127     if (iNode == abs(pF[iFace].edge[i].v)) break;
227   }                                               128   }
228   if (i == 4) {                                   129   if (i == 4) {
229     std::cerr                                  << 130     HepStd::cerr
230       << "HepPolyhedron::FindNeighbour: face "    131       << "HepPolyhedron::FindNeighbour: face " << iFace
231       << " has no node " << iNode                 132       << " has no node " << iNode
232       << std::endl;                            << 133       << HepStd::endl; 
233     return 0;                                     134     return 0;
234   }                                               135   }
235   if (iOrder < 0) {                               136   if (iOrder < 0) {
236     if ( --i < 0) i = 3;                          137     if ( --i < 0) i = 3;
237     if (pF[iFace].edge[i].v == 0) i = 2;          138     if (pF[iFace].edge[i].v == 0) i = 2;
238   }                                               139   }
239   return (pF[iFace].edge[i].v > 0) ? 0 : pF[iF    140   return (pF[iFace].edge[i].v > 0) ? 0 : pF[iFace].edge[i].f;
240 }                                                 141 }
241                                                   142 
242 G4Normal3D HepPolyhedron::FindNodeNormal(G4int << 143 HepNormal3D HepPolyhedron::FindNodeNormal(int iFace, int iNode) const
243 /*********************************************    144 /***********************************************************************
244  *                                                145  *                                                                     *
245  * Name: HepPolyhedron::FindNodeNormal            146  * Name: HepPolyhedron::FindNodeNormal               Date:    22.11.99 *
246  * Author: E.Chernyaev                            147  * Author: E.Chernyaev                               Revised:          *
247  *                                                148  *                                                                     *
248  * Function: Find normal at given node            149  * Function: Find normal at given node                                 *
249  *                                                150  *                                                                     *
250  *********************************************    151  ***********************************************************************/
251 {                                                 152 {
252   G4Normal3D normal = GetUnitNormal(iFace);    << 153   HepNormal3D  normal = GetUnitNormal(iFace);
253   G4int      k = iFace, iOrder = 1;            << 154   int          k = iFace, iOrder = 1, n = 1;
254                                                   155 
255   for(;;) {                                       156   for(;;) {
256     k = FindNeighbour(k, iNode, iOrder);          157     k = FindNeighbour(k, iNode, iOrder);
257     if (k == iFace) break;                     << 158     if (k == iFace) break; 
258     if (k > 0) {                                  159     if (k > 0) {
                                                   >> 160       n++;
259       normal += GetUnitNormal(k);                 161       normal += GetUnitNormal(k);
260     }else{                                        162     }else{
261       if (iOrder < 0) break;                      163       if (iOrder < 0) break;
262       k = iFace;                                  164       k = iFace;
263       iOrder = -iOrder;                           165       iOrder = -iOrder;
264     }                                             166     }
265   }                                               167   }
266   return normal.unit();                           168   return normal.unit();
267 }                                                 169 }
268                                                   170 
269 G4int HepPolyhedron::GetNumberOfRotationSteps( << 171 void HepPolyhedron::SetNumberOfRotationSteps(int n)
270 /********************************************* << 
271  *                                             << 
272  * Name: HepPolyhedron::GetNumberOfRotationSte << 
273  * Author: J.Allison (Manchester University)   << 
274  *                                             << 
275  * Function: Get number of steps for whole cir << 
276  *                                             << 
277  ********************************************* << 
278 {                                              << 
279   return fNumberOfRotationSteps;               << 
280 }                                              << 
281                                                << 
282 void HepPolyhedron::SetVertex(G4int index, con << 
283 /********************************************* << 
284  *                                             << 
285  * Name: HepPolyhedron::SetVertex              << 
286  * Author: E.Tcherniaev (E.Chernyaev)          << 
287  *                                             << 
288  * Function: Set vertex                        << 
289  *                                             << 
290  ********************************************* << 
291 {                                              << 
292   if (index < 1 || index > nvert)              << 
293   {                                            << 
294     std::cerr                                  << 
295       << "HepPolyhedron::SetVertex: vertex ind << 
296       << " is out of range\n"                  << 
297       << "   N. of vertices = " << nvert << "\ << 
298       << "   N. of facets = " << nface << std: << 
299     return;                                    << 
300   }                                            << 
301   pV[index] = v;                               << 
302 }                                              << 
303                                                << 
304 void                                           << 
305 HepPolyhedron::SetFacet(G4int index, G4int iv1 << 
306 /********************************************* << 
307  *                                             << 
308  * Name: HepPolyhedron::SetFacet               << 
309  * Author: E.Tcherniaev (E.Chernyaev)          << 
310  *                                             << 
311  * Function: Set facet                         << 
312  *                                             << 
313  ********************************************* << 
314 {                                              << 
315   if (index < 1 || index > nface)              << 
316   {                                            << 
317     std::cerr                                  << 
318       << "HepPolyhedron::SetFacet: facet index << 
319       << " is out of range\n"                  << 
320       << "   N. of vertices = " << nvert << "\ << 
321       << "   N. of facets = " << nface << std: << 
322     return;                                    << 
323   }                                            << 
324   if (iv1 < 1 || iv1 > nvert ||                << 
325       iv2 < 1 || iv2 > nvert ||                << 
326       iv3 < 1 || iv3 > nvert ||                << 
327       iv4 < 0 || iv4 > nvert)                  << 
328   {                                            << 
329     std::cerr                                  << 
330       << "HepPolyhedron::SetFacet: incorrectly << 
331       << " (" << iv1 << ", " << iv2 << ", " << << 
332       << "   N. of vertices = " << nvert << "\ << 
333       << "   N. of facets = " << nface << std: << 
334     return;                                    << 
335   }                                            << 
336   pF[index] = G4Facet(iv1, 0, iv2, 0, iv3, 0,  << 
337 }                                              << 
338                                                << 
339 void HepPolyhedron::SetNumberOfRotationSteps(G << 
340 /*********************************************    172 /***********************************************************************
341  *                                                173  *                                                                     *
342  * Name: HepPolyhedron::SetNumberOfRotationSte    174  * Name: HepPolyhedron::SetNumberOfRotationSteps     Date:    24.06.97 *
343  * Author: J.Allison (Manchester University)      175  * Author: J.Allison (Manchester University)         Revised:          *
344  *                                                176  *                                                                     *
345  * Function: Set number of steps for whole cir    177  * Function: Set number of steps for whole circle                      *
346  *                                                178  *                                                                     *
347  *********************************************    179  ***********************************************************************/
348 {                                                 180 {
349   const G4int nMin = 3;                        << 181   const int nMin = 3;
350   if (n < nMin) {                                 182   if (n < nMin) {
351     std::cerr                                  << 183     HepStd::cerr 
352       << "HepPolyhedron::SetNumberOfRotationSt    184       << "HepPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
353       << "number of steps per circle < " << nM    185       << "number of steps per circle < " << nMin << "; forced to " << nMin
354       << std::endl;                            << 186       << HepStd::endl;
355     fNumberOfRotationSteps = nMin;                187     fNumberOfRotationSteps = nMin;
356   }else{                                          188   }else{
357     fNumberOfRotationSteps = n;                   189     fNumberOfRotationSteps = n;
358   }                                            << 190   }    
359 }                                              << 
360                                                << 
361 void HepPolyhedron::ResetNumberOfRotationSteps << 
362 /********************************************* << 
363  *                                             << 
364  * Name: HepPolyhedron::GetNumberOfRotationSte << 
365  * Author: J.Allison (Manchester University)   << 
366  *                                             << 
367  * Function: Reset number of steps for whole c << 
368  *                                             << 
369  ********************************************* << 
370 {                                              << 
371   fNumberOfRotationSteps = DEFAULT_NUMBER_OF_S << 
372 }                                                 191 }
373                                                   192 
374 void HepPolyhedron::AllocateMemory(G4int Nvert << 193 void HepPolyhedron::AllocateMemory(int Nvert, int Nface)
375 /*********************************************    194 /***********************************************************************
376  *                                                195  *                                                                     *
377  * Name: HepPolyhedron::AllocateMemory            196  * Name: HepPolyhedron::AllocateMemory               Date:    19.06.96 *
378  * Author: E.Chernyaev (IHEP/Protvino)         << 197  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
379  *                                                198  *                                                                     *
380  * Function: Allocate memory for GEANT4 polyhe    199  * Function: Allocate memory for GEANT4 polyhedron                     *
381  *                                                200  *                                                                     *
382  * Input: Nvert - number of nodes                 201  * Input: Nvert - number of nodes                                      *
383  *        Nface - number of faces                 202  *        Nface - number of faces                                      *
384  *                                                203  *                                                                     *
385  *********************************************    204  ***********************************************************************/
386 {                                                 205 {
387   if (nvert == Nvert && nface == Nface) return << 206   nvert = Nvert;
388   delete [] pV;                                << 207   nface = Nface;
389   delete [] pF;                                << 208   pV    = new HepPoint3D[nvert+1];
390   if (Nvert > 0 && Nface > 0) {                << 209   pF    = new G4Facet[nface+1];
391     nvert = Nvert;                             << 
392     nface = Nface;                             << 
393     pV    = new G4Point3D[nvert+1];            << 
394     pF    = new G4Facet[nface+1];              << 
395   }else{                                       << 
396     nvert = 0; nface = 0; pV = nullptr; pF = n << 
397   }                                            << 
398 }                                                 210 }
399                                                   211 
400 void HepPolyhedron::CreatePrism()                 212 void HepPolyhedron::CreatePrism()
401 /*********************************************    213 /***********************************************************************
402  *                                                214  *                                                                     *
403  * Name: HepPolyhedron::CreatePrism               215  * Name: HepPolyhedron::CreatePrism                  Date:    15.07.96 *
404  * Author: E.Chernyaev (IHEP/Protvino)            216  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
405  *                                                217  *                                                                     *
406  * Function: Set facets for a prism               218  * Function: Set facets for a prism                                    *
407  *                                                219  *                                                                     *
408  *********************************************    220  ***********************************************************************/
409 {                                                 221 {
410   enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRON    222   enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
411                                                   223 
412   pF[1] = G4Facet(1,LEFT,  4,BACK,  3,RIGHT,      224   pF[1] = G4Facet(1,LEFT,  4,BACK,  3,RIGHT,  2,FRONT);
413   pF[2] = G4Facet(5,TOP,   8,BACK,  4,BOTTOM,     225   pF[2] = G4Facet(5,TOP,   8,BACK,  4,BOTTOM, 1,FRONT);
414   pF[3] = G4Facet(8,TOP,   7,RIGHT, 3,BOTTOM,     226   pF[3] = G4Facet(8,TOP,   7,RIGHT, 3,BOTTOM, 4,LEFT);
415   pF[4] = G4Facet(7,TOP,   6,FRONT, 2,BOTTOM,     227   pF[4] = G4Facet(7,TOP,   6,FRONT, 2,BOTTOM, 3,BACK);
416   pF[5] = G4Facet(6,TOP,   5,LEFT,  1,BOTTOM,     228   pF[5] = G4Facet(6,TOP,   5,LEFT,  1,BOTTOM, 2,RIGHT);
417   pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK,       229   pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK,   8,LEFT);
418 }                                                 230 }
419                                                   231 
420 void HepPolyhedron::RotateEdge(G4int k1, G4int << 232 void HepPolyhedron::RotateEdge(int k1, int k2, HepDouble r1, HepDouble r2,
421                               G4int v1, G4int  << 233             int v1, int v2, int vEdge,
422                               G4bool ifWholeCi << 234                               HepBoolean ifWholeCircle, int ns, int &kface)
423 /*********************************************    235 /***********************************************************************
424  *                                                236  *                                                                     *
425  * Name: HepPolyhedron::RotateEdge                237  * Name: HepPolyhedron::RotateEdge                   Date:    05.12.96 *
426  * Author: E.Chernyaev (IHEP/Protvino)            238  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
427  *                                                239  *                                                                     *
428  * Function: Create set of facets by rotation     240  * Function: Create set of facets by rotation of an edge around Z-axis *
429  *                                                241  *                                                                     *
430  * Input: k1, k2 - end vertices of the edge       242  * Input: k1, k2 - end vertices of the edge                            *
431  *        r1, r2 - radiuses of the end vertice    243  *        r1, r2 - radiuses of the end vertices                        *
432  *        v1, v2 - visibility of edges produce    244  *        v1, v2 - visibility of edges produced by rotation of the end *
433  *                 vertices                       245  *                 vertices                                            *
434  *        vEdge  - visibility of the edge         246  *        vEdge  - visibility of the edge                              *
435  *        ifWholeCircle - is true in case of w    247  *        ifWholeCircle - is true in case of whole circle rotation     *
436  *        nds    - number of discrete steps    << 248  *        ns     - number of discrete steps                            *
437  *        r[]    - r-coordinates                  249  *        r[]    - r-coordinates                                       *
438  *        kface  - current free cell in the pF    250  *        kface  - current free cell in the pF array                   *
439  *                                                251  *                                                                     *
440  *********************************************    252  ***********************************************************************/
441 {                                                 253 {
442   if (r1 == 0. && r2 == 0.) return;            << 254   if (r1 == 0. && r2 == 0) return;
443                                                   255 
444   G4int i;                                     << 256   int i;
445   G4int i1  = k1;                              << 257   int i1  = k1;
446   G4int i2  = k2;                              << 258   int i2  = k2;
447   G4int ii1 = ifWholeCircle ? i1 : i1+nds;     << 259   int ii1 = ifWholeCircle ? i1 : i1+ns;
448   G4int ii2 = ifWholeCircle ? i2 : i2+nds;     << 260   int ii2 = ifWholeCircle ? i2 : i2+ns;
449   G4int vv  = ifWholeCircle ? vEdge : 1;       << 261   int vv  = ifWholeCircle ? vEdge : 1;
450                                                   262 
451   if (nds == 1) {                              << 263   if (ns == 1) {
452     if (r1 == 0.) {                               264     if (r1 == 0.) {
453       pF[kface++]   = G4Facet(i1,0,    v2*i2,0    265       pF[kface++]   = G4Facet(i1,0,    v2*i2,0, (i2+1),0);
454     }else if (r2 == 0.) {                         266     }else if (r2 == 0.) {
455       pF[kface++]   = G4Facet(i1,0,    i2,0,      267       pF[kface++]   = G4Facet(i1,0,    i2,0,    v1*(i1+1),0);
456     }else{                                        268     }else{
457       pF[kface++]   = G4Facet(i1,0,    v2*i2,0    269       pF[kface++]   = G4Facet(i1,0,    v2*i2,0, (i2+1),0, v1*(i1+1),0);
458     }                                             270     }
459   }else{                                          271   }else{
460     if (r1 == 0.) {                               272     if (r1 == 0.) {
461       pF[kface++]   = G4Facet(vv*i1,0,    v2*i    273       pF[kface++]   = G4Facet(vv*i1,0,    v2*i2,0, vEdge*(i2+1),0);
462       for (i2++,i=1; i<nds-1; i2++,i++) {      << 274       for (i2++,i=1; i<ns-1; i2++,i++) {
463         pF[kface++] = G4Facet(vEdge*i1,0, v2*i << 275   pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
464       }                                           276       }
465       pF[kface++]   = G4Facet(vEdge*i1,0, v2*i    277       pF[kface++]   = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
466     }else if (r2 == 0.) {                         278     }else if (r2 == 0.) {
467       pF[kface++]   = G4Facet(vv*i1,0,    vEdg    279       pF[kface++]   = G4Facet(vv*i1,0,    vEdge*i2,0, v1*(i1+1),0);
468       for (i1++,i=1; i<nds-1; i1++,i++) {      << 280       for (i1++,i=1; i<ns-1; i1++,i++) {
469         pF[kface++] = G4Facet(vEdge*i1,0, vEdg << 281   pF[kface++] = G4Facet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
470       }                                           282       }
471       pF[kface++]   = G4Facet(vEdge*i1,0, vv*i    283       pF[kface++]   = G4Facet(vEdge*i1,0, vv*i2,0,    v1*ii1,0);
472     }else{                                        284     }else{
473       pF[kface++]   = G4Facet(vv*i1,0,    v2*i    285       pF[kface++]   = G4Facet(vv*i1,0,    v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
474       for (i1++,i2++,i=1; i<nds-1; i1++,i2++,i << 286       for (i1++,i2++,i=1; i<ns-1; i1++,i2++,i++) {
475         pF[kface++] = G4Facet(vEdge*i1,0, v2*i << 287   pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
476       }                                        << 288       }  
477       pF[kface++]   = G4Facet(vEdge*i1,0, v2*i    289       pF[kface++]   = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0,      v1*ii1,0);
478     }                                             290     }
479   }                                               291   }
480 }                                                 292 }
481                                                   293 
482 void HepPolyhedron::SetSideFacets(G4int ii[4], << 294 void HepPolyhedron::SetSideFacets(int ii[4], int vv[4],
483                                  G4int *kk, G4 << 295          int *kk, HepDouble *r,
484                                  G4double dphi << 296                                  HepDouble dphi, int ns, int &kface)
485 /*********************************************    297 /***********************************************************************
486  *                                                298  *                                                                     *
487  * Name: HepPolyhedron::SetSideFacets             299  * Name: HepPolyhedron::SetSideFacets                Date:    20.05.97 *
488  * Author: E.Chernyaev (IHEP/Protvino)            300  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
489  *                                                301  *                                                                     *
490  * Function: Set side facets for the case of i    302  * Function: Set side facets for the case of incomplete rotation       *
491  *                                                303  *                                                                     *
492  * Input: ii[4] - indices of original vertices << 304  * Input: ii[4] - indeces of original verteces                         *
493  *        vv[4] - visibility of edges             305  *        vv[4] - visibility of edges                                  *
494  *        kk[]  - indices of nodes             << 306  *        kk[]  - indeces of nodes                                     *
495  *        r[]   - radiuses                        307  *        r[]   - radiuses                                             *
496  *        dphi  - delta phi                       308  *        dphi  - delta phi                                            *
497  *        nds    - number of discrete steps    << 309  *        ns     - number of discrete steps                            *
498  *        kface  - current free cell in the pF    310  *        kface  - current free cell in the pF array                   *
499  *                                                311  *                                                                     *
500  *********************************************    312  ***********************************************************************/
501 {                                                 313 {
502   G4int k1, k2, k3, k4;                        << 314   int k1, k2, k3, k4;
503                                                << 315   
504   if (std::abs(dphi-pi) < perMillion) { // hal << 316   if (abs((HepDouble)(dphi-M_PI)) < perMillion) {          // half a circle
505     for (G4int i=0; i<4; i++) {                << 317     for (int i=0; i<4; i++) {
506       k1 = ii[i];                                 318       k1 = ii[i];
507       k2 = ii[(i+1)%4];                        << 319       k2 = (i == 3) ? ii[0] : ii[i+1];
508       if (r[k1] == 0. && r[k2] == 0.) vv[i] =  << 320       if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;      
509     }                                             321     }
510   }                                               322   }
511                                                   323 
512   if (ii[1] == ii[2]) {                           324   if (ii[1] == ii[2]) {
513     k1 = kk[ii[0]];                               325     k1 = kk[ii[0]];
514     k2 = kk[ii[2]];                               326     k2 = kk[ii[2]];
515     k3 = kk[ii[3]];                               327     k3 = kk[ii[3]];
516     pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2    328     pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
517     if (r[ii[0]] != 0.) k1 += nds;             << 329     if (r[ii[0]] != 0.) k1 += ns;
518     if (r[ii[2]] != 0.) k2 += nds;             << 330     if (r[ii[2]] != 0.) k2 += ns;
519     if (r[ii[3]] != 0.) k3 += nds;             << 331     if (r[ii[3]] != 0.) k3 += ns;
520     pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2    332     pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
521   }else if (kk[ii[0]] == kk[ii[1]]) {             333   }else if (kk[ii[0]] == kk[ii[1]]) {
522     k1 = kk[ii[0]];                               334     k1 = kk[ii[0]];
523     k2 = kk[ii[2]];                               335     k2 = kk[ii[2]];
524     k3 = kk[ii[3]];                               336     k3 = kk[ii[3]];
525     pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2    337     pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
526     if (r[ii[0]] != 0.) k1 += nds;             << 338     if (r[ii[0]] != 0.) k1 += ns;
527     if (r[ii[2]] != 0.) k2 += nds;             << 339     if (r[ii[2]] != 0.) k2 += ns;
528     if (r[ii[3]] != 0.) k3 += nds;             << 340     if (r[ii[3]] != 0.) k3 += ns;
529     pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2    341     pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
530   }else if (kk[ii[2]] == kk[ii[3]]) {             342   }else if (kk[ii[2]] == kk[ii[3]]) {
531     k1 = kk[ii[0]];                               343     k1 = kk[ii[0]];
532     k2 = kk[ii[1]];                               344     k2 = kk[ii[1]];
533     k3 = kk[ii[2]];                               345     k3 = kk[ii[2]];
534     pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2    346     pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
535     if (r[ii[0]] != 0.) k1 += nds;             << 347     if (r[ii[0]] != 0.) k1 += ns;
536     if (r[ii[1]] != 0.) k2 += nds;             << 348     if (r[ii[1]] != 0.) k2 += ns;
537     if (r[ii[2]] != 0.) k3 += nds;             << 349     if (r[ii[2]] != 0.) k3 += ns;
538     pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2    350     pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
539   }else{                                          351   }else{
540     k1 = kk[ii[0]];                               352     k1 = kk[ii[0]];
541     k2 = kk[ii[1]];                               353     k2 = kk[ii[1]];
542     k3 = kk[ii[2]];                               354     k3 = kk[ii[2]];
543     k4 = kk[ii[3]];                               355     k4 = kk[ii[3]];
544     pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2    356     pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
545     if (r[ii[0]] != 0.) k1 += nds;             << 357     if (r[ii[0]] != 0.) k1 += ns;
546     if (r[ii[1]] != 0.) k2 += nds;             << 358     if (r[ii[1]] != 0.) k2 += ns;
547     if (r[ii[2]] != 0.) k3 += nds;             << 359     if (r[ii[2]] != 0.) k3 += ns;
548     if (r[ii[3]] != 0.) k4 += nds;             << 360     if (r[ii[3]] != 0.) k4 += ns;
549     pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3    361     pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
550   }                                               362   }
551 }                                                 363 }
552                                                   364 
553 void HepPolyhedron::RotateAroundZ(G4int nstep, << 365 void HepPolyhedron::RotateAroundZ(int nstep, HepDouble phi, HepDouble dphi,
554                                  G4int np1, G4 << 366                                  int np1, int np2,
555                                  const G4doubl << 367          const HepDouble *z, HepDouble *r,
556                                  G4int nodeVis << 368          int nodeVis, int edgeVis)
557 /*********************************************    369 /***********************************************************************
558  *                                                370  *                                                                     *
559  * Name: HepPolyhedron::RotateAroundZ             371  * Name: HepPolyhedron::RotateAroundZ                Date:    27.11.96 *
560  * Author: E.Chernyaev (IHEP/Protvino)            372  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
561  *                                                373  *                                                                     *
562  * Function: Create HepPolyhedron for a solid     374  * Function: Create HepPolyhedron for a solid produced by rotation of  *
563  *           two polylines around Z-axis          375  *           two polylines around Z-axis                               *
564  *                                                376  *                                                                     *
565  * Input: nstep - number of discrete steps, if    377  * Input: nstep - number of discrete steps, if 0 then default          *
566  *        phi   - starting phi angle              378  *        phi   - starting phi angle                                   *
567  *        dphi  - delta phi                       379  *        dphi  - delta phi                                            *
568  *        np1   - number of points in external    380  *        np1   - number of points in external polyline                *
569  *                (must be negative in case of    381  *                (must be negative in case of closed polyline)        *
570  *        np2   - number of points in internal    382  *        np2   - number of points in internal polyline (may be 1)     *
571  *        z[]   - z-coordinates (+z >>> -z for    383  *        z[]   - z-coordinates (+z >>> -z for both polylines)         *
572  *        r[]   - r-coordinates                   384  *        r[]   - r-coordinates                                        *
573  *        nodeVis - how to Draw edges joing co    385  *        nodeVis - how to Draw edges joing consecutive positions of   *
574  *                  node during rotation          386  *                  node during rotation                               *
575  *        edgeVis - how to Draw edges             387  *        edgeVis - how to Draw edges                                  *
576  *                                                388  *                                                                     *
577  *********************************************    389  ***********************************************************************/
578 {                                                 390 {
579   static const G4double wholeCircle   = twopi; << 391   static HepDouble wholeCircle   = 2*M_PI;
580                                                << 392     
581   //   S E T   R O T A T I O N   P A R A M E T    393   //   S E T   R O T A T I O N   P A R A M E T E R S
582                                                   394 
583   G4bool ifWholeCircle = std::abs(dphi-wholeCi << 395   HepBoolean ifWholeCircle = (abs(dphi-wholeCircle) < perMillion) ?
584   G4double delPhi = ifWholeCircle ? wholeCircl << 396     true : false;
585   G4int nSphi = nstep;                         << 397   HepDouble   delPhi  = ifWholeCircle ? wholeCircle : dphi;  
586   if (nSphi <= 0) nSphi = GetNumberOfRotationS << 398   int        nSphi    = (nstep > 0) ?
                                                   >> 399     nstep : int(delPhi*GetNumberOfRotationSteps()/wholeCircle+.5);
587   if (nSphi == 0) nSphi = 1;                      400   if (nSphi == 0) nSphi = 1;
588   G4int nVphi = ifWholeCircle ? nSphi : nSphi  << 401   int        nVphi    = ifWholeCircle ? nSphi : nSphi+1;
589   G4bool ifClosed = np1 <= 0; // true if exter << 402   HepBoolean ifClosed = np1 > 0 ? false : true;
590                                                << 403   
591   //   C O U N T   V E R T I C E S             << 404   //   C O U N T   V E R T E C E S
592                                                   405 
593   G4int absNp1 = std::abs(np1);                << 406   int absNp1 = abs(np1);
594   G4int absNp2 = std::abs(np2);                << 407   int absNp2 = abs(np2);
595   G4int i1beg = 0;                             << 408   int i1beg = 0;
596   G4int i1end = absNp1-1;                      << 409   int i1end = absNp1-1;
597   G4int i2beg = absNp1;                        << 410   int i2beg = absNp1;
598   G4int i2end = absNp1+absNp2-1;               << 411   int i2end = absNp1+absNp2-1; 
599   G4int i, j, k;                               << 412   int i, j, k;
600                                                   413 
601   for(i=i1beg; i<=i2end; i++) {                   414   for(i=i1beg; i<=i2end; i++) {
602     if (std::abs(r[i]) < spatialTolerance) r[i << 415     if (abs(r[i]) < perMillion) r[i] = 0.;
603   }                                               416   }
604                                                   417 
605   // external polyline - check position of nod << 418   j = 0;                                                // external nodes
606   //                                           << 
607   G4int Nverts = 0;                            << 
608   for (i=i1beg; i<=i1end; i++) {                  419   for (i=i1beg; i<=i1end; i++) {
609     Nverts += (r[i] == 0.) ? 1 : nVphi;        << 420     j += (r[i] == 0.) ? 1 : nVphi;
610   }                                               421   }
611                                                   422 
612   // internal polyline                         << 423   HepBoolean ifSide1 = false;                           // internal nodes
613   //                                           << 424   HepBoolean ifSide2 = false;
614   G4bool ifSide1 = false; // whether to create << 
615   G4bool ifSide2 = false; // whether to create << 
616                                                   425 
617   if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1 << 426   if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
618     Nverts += (r[i2beg] == 0.) ? 1 : nVphi;    << 427     j += (r[i2beg] == 0.) ? 1 : nVphi;
619     ifSide1 = true;                               428     ifSide1 = true;
620   }                                               429   }
621                                                   430 
622   for(i=i2beg+1; i<i2end; i++) { // intermedia << 431   for(i=i2beg+1; i<i2end; i++) {
623     Nverts += (r[i] == 0.) ? 1 : nVphi;        << 432     j += (r[i] == 0.) ? 1 : nVphi;
624   }                                               433   }
625                                                << 434   
626   if (r[i2end] != r[i1end] || z[i2end] != z[i1 << 435   if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
627     if (absNp2 > 1) Nverts += (r[i2end] == 0.) << 436     if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
628     ifSide2 = true;                               437     ifSide2 = true;
629   }                                               438   }
630                                                   439 
631   //   C O U N T   F A C E S                      440   //   C O U N T   F A C E S
632                                                   441 
633   // external lateral faces                    << 442   k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi;       // external faces
634   //                                           << 
635   G4int Nfaces = ifClosed ? absNp1*nSphi : (ab << 
636                                                   443 
637   // internal lateral faces                    << 444   if (absNp2 > 1) {                                     // internal faces
638   //                                           << 
639   if (absNp2 > 1) {                            << 
640     for(i=i2beg; i<i2end; i++) {                  445     for(i=i2beg; i<i2end; i++) {
641       if (r[i] > 0. || r[i+1] > 0.) Nfaces +=  << 446       if (r[i] > 0. || r[i+1] > 0.)       k += nSphi;
642     }                                             447     }
643                                                   448 
644     if (ifClosed) {                               449     if (ifClosed) {
645       if (r[i2end] > 0. || r[i2beg] > 0.) Nfac << 450       if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
646     }                                             451     }
647   }                                               452   }
648                                                   453 
649   // bottom and top faces                      << 454   if (!ifClosed) {                                      // side faces
650   //                                           << 455     if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
651   if (!ifClosed) {                             << 456     if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) k += nSphi;
652     if (ifSide1 && (r[i1beg] > 0. || r[i2beg]  << 
653     if (ifSide2 && (r[i1end] > 0. || r[i2end]  << 
654   }                                               457   }
655                                                   458 
656   // phi_wedge faces                           << 459   if (!ifWholeCircle) {                                 // phi_side faces
657   //                                           << 460     k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
658   if (!ifWholeCircle) {                        << 
659     Nfaces += ifClosed ? 2*absNp1 : 2*(absNp1- << 
660   }                                               461   }
661                                                   462 
662   //   A L L O C A T E   M E M O R Y              463   //   A L L O C A T E   M E M O R Y
663                                                   464 
664   AllocateMemory(Nverts, Nfaces);              << 465   AllocateMemory(j, k);
665   if (pV == nullptr || pF == nullptr) return;  << 
666                                                   466 
667   //   G E N E R A T E   V E R T I C E S       << 467   //   G E N E R A T E   V E R T E C E S
668                                                   468 
669   G4int *kk; // array of start indices along p << 469   int *kk;
670   kk = new G4int[absNp1+absNp2];               << 470   kk = new int[absNp1+absNp2];
671                                                   471 
672   // external polyline                         << 472   k = 1;
673   //                                           << 
674   k = 1; // free position in array of vertices << 
675   for(i=i1beg; i<=i1end; i++) {                   473   for(i=i1beg; i<=i1end; i++) {
676     kk[i] = k;                                    474     kk[i] = k;
677     if (r[i] == 0.)                            << 475     if (r[i] == 0.) { pV[k++] = HepPoint3D(0, 0, z[i]); } else { k += nVphi; }
678     { pV[k++] = G4Point3D(0, 0, z[i]); } else  << 
679   }                                               476   }
680                                                   477 
681   // first point of internal polyline          << 
682   //                                           << 
683   i = i2beg;                                      478   i = i2beg;
684   if (ifSide1) {                                  479   if (ifSide1) {
685     kk[i] = k;                                    480     kk[i] = k;
686     if (r[i] == 0.)                            << 481     if (r[i] == 0.) { pV[k++] = HepPoint3D(0, 0, z[i]); } else { k += nVphi; }
687     { pV[k++] = G4Point3D(0, 0, z[i]); } else  << 
688   }else{                                          482   }else{
689     kk[i] = kk[i1beg];                            483     kk[i] = kk[i1beg];
690   }                                               484   }
691                                                   485 
692   // intermediate points of internal polyline  << 
693   //                                           << 
694   for(i=i2beg+1; i<i2end; i++) {                  486   for(i=i2beg+1; i<i2end; i++) {
695     kk[i] = k;                                    487     kk[i] = k;
696     if (r[i] == 0.)                            << 488     if (r[i] == 0.) { pV[k++] = HepPoint3D(0, 0, z[i]); } else { k += nVphi; }
697     { pV[k++] = G4Point3D(0, 0, z[i]); } else  << 
698   }                                               489   }
699                                                   490 
700   // last point of internal polyline           << 
701   //                                           << 
702   if (absNp2 > 1) {                               491   if (absNp2 > 1) {
703     i = i2end;                                    492     i = i2end;
704     if (ifSide2) {                                493     if (ifSide2) {
705       kk[i] = k;                                  494       kk[i] = k;
706       if (r[i] == 0.) pV[k] = G4Point3D(0, 0,  << 495       if (r[i] == 0.) pV[k] = HepPoint3D(0, 0, z[i]);
707     }else{                                        496     }else{
708       kk[i] = kk[i1end];                          497       kk[i] = kk[i1end];
709     }                                             498     }
710   }                                               499   }
711                                                   500 
712   // set vertices                              << 501   HepDouble cosPhi, sinPhi;
713   //                                           << 
714   G4double cosPhi, sinPhi;                     << 
715                                                   502 
716   for(j=0; j<nVphi; j++) {                        503   for(j=0; j<nVphi; j++) {
717     cosPhi = std::cos(phi+j*delPhi/nSphi);     << 504     cosPhi = cos(phi+j*delPhi/nSphi);
718     sinPhi = std::sin(phi+j*delPhi/nSphi);     << 505     sinPhi = sin(phi+j*delPhi/nSphi);
719     for(i=i1beg; i<=i2end; i++) {                 506     for(i=i1beg; i<=i2end; i++) {
720       if (r[i] != 0.)                          << 507       if (r[i] != 0.) pV[kk[i]+j] = HepPoint3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
721         pV[kk[i]+j] = G4Point3D(r[i]*cosPhi,r[ << 
722     }                                             508     }
723   }                                               509   }
724                                                   510 
725   //   G E N E R A T E   F A C E S             << 511   //   G E N E R A T E   E X T E R N A L   F A C E S
726                                                   512 
727   //  external faces                           << 513   int v1,v2;
728   //                                           << 
729   G4int v1,v2;                                 << 
730                                                   514 
731   k = 1; // free position in array of faces pF << 515   k = 1;
732   v2 = ifClosed ? nodeVis : 1;                    516   v2 = ifClosed ? nodeVis : 1;
733   for(i=i1beg; i<i1end; i++) {                    517   for(i=i1beg; i<i1end; i++) {
734     v1 = v2;                                      518     v1 = v2;
735     if (!ifClosed && i == i1end-1) {              519     if (!ifClosed && i == i1end-1) {
736       v2 = 1;                                     520       v2 = 1;
737     }else{                                        521     }else{
738       v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]    522       v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
739     }                                             523     }
740     RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v    524     RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
741                edgeVis, ifWholeCircle, nSphi,  << 525          edgeVis, ifWholeCircle, nSphi, k);
742   }                                               526   }
743   if (ifClosed) {                                 527   if (ifClosed) {
744     RotateEdge(kk[i1end], kk[i1beg], r[i1end],    528     RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
745                edgeVis, ifWholeCircle, nSphi,  << 529          edgeVis, ifWholeCircle, nSphi, k);
746   }                                               530   }
747                                                   531 
748   // internal faces                            << 532   //   G E N E R A T E   I N T E R N A L   F A C E S
749   //                                           << 533 
750   if (absNp2 > 1) {                               534   if (absNp2 > 1) {
751     v2 = ifClosed ? nodeVis : 1;                  535     v2 = ifClosed ? nodeVis : 1;
752     for(i=i2beg; i<i2end; i++) {                  536     for(i=i2beg; i<i2end; i++) {
753       v1 = v2;                                    537       v1 = v2;
754       if (!ifClosed && i==i2end-1) {              538       if (!ifClosed && i==i2end-1) {
755         v2 = 1;                                << 539   v2 = 1;
756       }else{                                      540       }else{
757         v2 = (r[i] == r[i+1] && r[i+1] == r[i+ << 541   v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 :  nodeVis;
758       }                                           542       }
759       RotateEdge(kk[i+1], kk[i], r[i+1], r[i],    543       RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
760                  edgeVis, ifWholeCircle, nSphi << 544      edgeVis, ifWholeCircle, nSphi, k);
761     }                                             545     }
762     if (ifClosed) {                               546     if (ifClosed) {
763       RotateEdge(kk[i2beg], kk[i2end], r[i2beg    547       RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
764                  edgeVis, ifWholeCircle, nSphi << 548      edgeVis, ifWholeCircle, nSphi, k);
765     }                                             549     }
766   }                                               550   }
767                                                   551 
768   // bottom and top faces                      << 552   //   G E N E R A T E   S I D E   F A C E S
769   //                                           << 553 
770   if (!ifClosed) {                                554   if (!ifClosed) {
771     if (ifSide1) {                                555     if (ifSide1) {
772       RotateEdge(kk[i2beg], kk[i1beg], r[i2beg    556       RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
773                  -1, ifWholeCircle, nSphi, k); << 557      -1, ifWholeCircle, nSphi, k);
774     }                                             558     }
775     if (ifSide2) {                                559     if (ifSide2) {
776       RotateEdge(kk[i1end], kk[i2end], r[i1end    560       RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
777                  -1, ifWholeCircle, nSphi, k); << 561      -1, ifWholeCircle, nSphi, k);
778     }                                             562     }
779   }                                               563   }
780                                                   564 
781   // phi_wedge faces in case of incomplete cir << 565   //   G E N E R A T E   S I D E   F A C E S  for the case of incomplete circle
782   //                                           << 566 
783   if (!ifWholeCircle) {                           567   if (!ifWholeCircle) {
784                                                   568 
785     G4int  ii[4], vv[4];                       << 569     int  ii[4], vv[4];
786                                                   570 
787     if (ifClosed) {                               571     if (ifClosed) {
788       for (i=i1beg; i<=i1end; i++) {              572       for (i=i1beg; i<=i1end; i++) {
789         ii[0] = i;                             << 573   ii[0] = i;
790         ii[3] = (i == i1end) ? i1beg : i+1;    << 574   ii[3] = (i == i1end) ? i1beg : i+1;
791         ii[1] = (absNp2 == 1) ? i2beg : ii[0]+ << 575   ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
792         ii[2] = (absNp2 == 1) ? i2beg : ii[3]+ << 576   ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
793         vv[0] = -1;                            << 577   vv[0] = -1;
794         vv[1] = 1;                             << 578   vv[1] = 1;
795         vv[2] = -1;                            << 579   vv[2] = -1;
796         vv[3] = 1;                             << 580   vv[3] = 1;
797         SetSideFacets(ii, vv, kk, r, delPhi, n << 581   SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
798       }                                           582       }
799     }else{                                        583     }else{
800       for (i=i1beg; i<i1end; i++) {               584       for (i=i1beg; i<i1end; i++) {
801         ii[0] = i;                             << 585   ii[0] = i;
802         ii[3] = i+1;                           << 586   ii[3] = i+1;
803         ii[1] = (absNp2 == 1) ? i2beg : ii[0]+ << 587   ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
804         ii[2] = (absNp2 == 1) ? i2beg : ii[3]+ << 588   ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
805         vv[0] = (i == i1beg)   ? 1 : -1;       << 589   vv[0] = (i == i1beg)   ? 1 : -1;
806         vv[1] = 1;                             << 590   vv[1] = 1;
807         vv[2] = (i == i1end-1) ? 1 : -1;       << 591   vv[2] = (i == i1end-1) ? 1 : -1;
808         vv[3] = 1;                             << 592   vv[3] = 1;
809         SetSideFacets(ii, vv, kk, r, delPhi, n << 593   SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
810       }                                           594       }
811     }                                          << 595     }      
812   }                                               596   }
813                                                   597 
814   delete [] kk; // free memory                 << 598   delete [] kk;
815                                                   599 
816   // final check                               << 
817   //                                           << 
818   if (k-1 != nface) {                             600   if (k-1 != nface) {
819     std::cerr                                  << 601     HepStd::cerr
820       << "HepPolyhedron::RotateAroundZ: number << 602       << "Polyhedron::RotateAroundZ: number of generated faces ("
821       << k-1 << ") is not equal to the number     603       << k-1 << ") is not equal to the number of allocated faces ("
822       << nface << ")"                             604       << nface << ")"
823       << std::endl;                            << 605       << HepStd::endl;
824   }                                            << 
825 }                                              << 
826                                                << 
827 void                                           << 
828 HepPolyhedron::RotateContourAroundZ(G4int nste << 
829                                     G4double p << 
830                                     G4double d << 
831                                     const std: << 
832                                     G4int node << 
833                                     G4int edge << 
834 /********************************************* << 
835  *                                             << 
836  * Name: HepPolyhedron::RotateContourAroundZ   << 
837  * Author: E.Tcherniaev (E.Chernyaev)          << 
838  *                                             << 
839  * Function: Create HepPolyhedron for a solid  << 
840  *           a closed polyline (rz-contour) ar << 
841  *                                             << 
842  * Input: nstep - number of discrete steps, if << 
843  *        phi   - starting phi angle           << 
844  *        dphi  - delta phi                    << 
845  *        rz    - rz-contour                   << 
846  *        nodeVis - how to Draw edges joing co << 
847  *                  node during rotation       << 
848  *        edgeVis - how to Draw edges          << 
849  *                                             << 
850  ********************************************* << 
851 {                                              << 
852   //   S E T   R O T A T I O N   P A R A M E T << 
853                                                << 
854   G4bool ifWholeCircle = std::abs(dphi - twopi << 
855   G4double delPhi = (ifWholeCircle) ? twopi :  << 
856   G4int nSphi = nstep;                         << 
857   if (nSphi <= 0) nSphi = GetNumberOfRotationS << 
858   if (nSphi == 0) nSphi = 1;                   << 
859   G4int nVphi = (ifWholeCircle) ? nSphi : nSph << 
860                                                << 
861   //   C A L C U L A T E   A R E A             << 
862                                                << 
863   G4int Nrz = (G4int)rz.size();                << 
864   G4double area = 0;                           << 
865   for (G4int i = 0; i < Nrz; ++i)              << 
866   {                                            << 
867     G4int k = (i == 0) ? Nrz - 1 : i - 1;      << 
868     area += rz[k].x()*rz[i].y() - rz[i].x()*rz << 
869   }                                            << 
870                                                << 
871   //   P R E P A R E   P O L Y L I N E         << 
872                                                << 
873   auto r = new G4double[Nrz];                  << 
874   auto z = new G4double[Nrz];                  << 
875   for (G4int i = 0; i < Nrz; ++i)              << 
876   {                                            << 
877     r[i] = rz[i].x();                          << 
878     z[i] = rz[i].y();                          << 
879     if (std::abs(r[i]) < spatialTolerance) r[i << 
880   }                                            << 
881                                                << 
882   //   C O U N T   V E R T I C E S   A N D   F << 
883                                                << 
884   G4int Nverts = 0;                            << 
885   for(G4int i = 0; i < Nrz; ++i) Nverts += (r[ << 
886                                                << 
887   G4int Nedges = Nrz;                          << 
888   for (G4int i = 0; i < Nrz; ++i)              << 
889   {                                            << 
890     G4int k = (i == 0) ? Nrz - 1 : i - 1;      << 
891     Nedges -= static_cast<int>(r[k] == 0 && r[ << 
892   }                                            << 
893                                                << 
894   G4int Nfaces = Nedges*nSphi;               / << 
895   if (!ifWholeCircle) Nfaces += 2*(Nrz - 2); / << 
896                                                << 
897   //   A L L O C A T E   M E M O R Y           << 
898                                                << 
899   AllocateMemory(Nverts, Nfaces);              << 
900   if (pV == nullptr || pF == nullptr)          << 
901   {                                            << 
902     delete [] r;                               << 
903     delete [] z;                               << 
904     return;                                    << 
905   }                                            << 
906                                                << 
907   //   S E T   V E R T I C E S                 << 
908                                                << 
909   auto kk = new G4int[Nrz]; // start indices a << 
910   G4int kfree = 1; // current free position in << 
911                                                << 
912   // set start indices, set vertices for nodes << 
913   for(G4int i = 0; i < Nrz; ++i)               << 
914   {                                            << 
915     kk[i] = kfree;                             << 
916     if (r[i] == 0.) pV[kfree++] = G4Point3D(0, << 
917     if (r[i] != 0.) kfree += nVphi;            << 
918   }                                            << 
919                                                << 
920   // set vertices by rotating r                << 
921   for(G4int j = 0; j < nVphi; ++j)             << 
922   {                                            << 
923     G4double cosPhi = std::cos(phi + j*delPhi/ << 
924     G4double sinPhi = std::sin(phi + j*delPhi/ << 
925     for(G4int i = 0; i < Nrz; ++i)             << 
926     {                                          << 
927       if (r[i] != 0.)                          << 
928         pV[kk[i] + j] = G4Point3D(r[i]*cosPhi, << 
929     }                                          << 
930   }                                            << 
931                                                << 
932   //   S E T   F A C E S                       << 
933                                                << 
934   kfree = 1; // current free position in array << 
935   for(G4int i = 0; i < Nrz; ++i)               << 
936   {                                            << 
937     G4int i1 = (i < Nrz - 1) ? i + 1 : 0; // i << 
938     G4int i2 = i;                              << 
939     if (area < 0.) std::swap(i1, i2);          << 
940     RotateEdge(kk[i1], kk[i2], r[i1], r[i2], n << 
941                edgeVis, ifWholeCircle, nSphi,  << 
942   }                                            << 
943                                                << 
944   //    S E T   P H I _ W E D G E   F A C E S  << 
945                                                << 
946   if (!ifWholeCircle)                          << 
947   {                                            << 
948     std::vector<G4int> triangles;              << 
949     TriangulatePolygon(rz, triangles);         << 
950                                                << 
951     G4int ii[4], vv[4];                        << 
952     G4int ntria = G4int(triangles.size()/3);   << 
953     for (G4int i = 0; i < ntria; ++i)          << 
954     {                                          << 
955       G4int i1 = triangles[0 + i*3];           << 
956       G4int i2 = triangles[1 + i*3];           << 
957       G4int i3 = triangles[2 + i*3];           << 
958       if (area < 0.) std::swap(i1, i3);        << 
959       G4int v1 = (std::abs(i2-i1) == 1 || std: << 
960       G4int v2 = (std::abs(i3-i2) == 1 || std: << 
961       G4int v3 = (std::abs(i1-i3) == 1 || std: << 
962       ii[0] = i1; ii[1] = i2; ii[2] = i2; ii[3 << 
963       vv[0] = v1; vv[1] = -1; vv[2] = v2; vv[3 << 
964       SetSideFacets(ii, vv, kk, r, delPhi, nSp << 
965     }                                          << 
966   }                                            << 
967                                                << 
968   // free memory                               << 
969   delete [] r;                                 << 
970   delete [] z;                                 << 
971   delete [] kk;                                << 
972                                                << 
973   // final check                               << 
974   if (kfree - 1 != nface)                      << 
975   {                                            << 
976     std::cerr                                  << 
977       << "HepPolyhedron::RotateContourAroundZ: << 
978       << kfree-1 << ") is not equal to the num << 
979       << nface << ")"                          << 
980       << std::endl;                            << 
981   }                                            << 
982 }                                              << 
983                                                << 
984 G4bool                                         << 
985 HepPolyhedron::TriangulatePolygon(const std::v << 
986                                   std::vector< << 
987 /********************************************* << 
988  *                                             << 
989  * Name: HepPolyhedron::TriangulatePolygon     << 
990  * Author: E.Tcherniaev (E.Chernyaev)          << 
991  *                                             << 
992  * Function: Simple implementation of "ear cli << 
993  *           triangulation of a simple contour << 
994  *           the result in a std::vector as tr << 
995  *                                             << 
996  *           If triangulation is sucsessfull t << 
997  *           returns true, otherwise false     << 
998  *                                             << 
999  * Remark:   It's a copy of G4GeomTools::Trian << 
1000  *                                            << 
1001  ******************************************** << 
1002 {                                             << 
1003   result.resize(0);                           << 
1004   G4int n = (G4int)polygon.size();            << 
1005   if (n < 3) return false;                    << 
1006                                               << 
1007   // calculate area                           << 
1008   //                                          << 
1009   G4double area = 0.;                         << 
1010   for(G4int i = 0; i < n; ++i)                << 
1011   {                                           << 
1012     G4int k = (i == 0) ? n - 1 : i - 1;       << 
1013     area += polygon[k].x()*polygon[i].y() - p << 
1014   }                                           << 
1015                                               << 
1016   // allocate and initialize list of Vertices << 
1017   // we want a counter-clockwise polygon in V << 
1018   //                                          << 
1019   auto  V = new G4int[n];                     << 
1020   if (area > 0.)                              << 
1021     for (G4int i = 0; i < n; ++i) V[i] = i;   << 
1022   else                                        << 
1023     for (G4int i = 0; i < n; ++i) V[i] = (n - << 
1024                                               << 
1025   //  Triangulation: remove nv-2 Vertices, cr << 
1026   //                                          << 
1027   G4int nv = n;                               << 
1028   G4int count = 2*nv; // error detection coun << 
1029   for(G4int b = nv - 1; nv > 2; )             << 
1030   {                                           << 
1031     // ERROR: if we loop, it is probably a no << 
1032     if ((count--) <= 0)                       << 
1033     {                                         << 
1034       delete [] V;                            << 
1035       if (area < 0.) std::reverse(result.begi << 
1036       return false;                           << 
1037     }                                         << 
1038                                               << 
1039     // three consecutive vertices in current  << 
1040     G4int a = (b   < nv) ? b   : 0; // previo << 
1041           b = (a+1 < nv) ? a+1 : 0; // curren << 
1042     G4int c = (b+1 < nv) ? b+1 : 0; // next   << 
1043                                               << 
1044     if (CheckSnip(polygon, a,b,c, nv,V))      << 
1045     {                                         << 
1046       // output Triangle                      << 
1047       result.push_back(V[a]);                 << 
1048       result.push_back(V[b]);                 << 
1049       result.push_back(V[c]);                 << 
1050                                               << 
1051       // remove vertex b from remaining polyg << 
1052       nv--;                                   << 
1053       for(G4int i = b; i < nv; ++i) V[i] = V[ << 
1054                                               << 
1055       count = 2*nv; // resest error detection << 
1056     }                                         << 
1057   }                                           << 
1058   delete [] V;                                << 
1059   if (area < 0.) std::reverse(result.begin(), << 
1060   return true;                                << 
1061 }                                             << 
1062                                               << 
1063 G4bool HepPolyhedron::CheckSnip(const std::ve << 
1064                                 G4int a, G4in << 
1065                                 G4int n, cons << 
1066 /******************************************** << 
1067  *                                            << 
1068  * Name: HepPolyhedron::CheckSnip             << 
1069  * Author: E.Tcherniaev (E.Chernyaev)         << 
1070  *                                            << 
1071  * Function: Check for a valid snip,          << 
1072  *           it is a helper functionfor Trian << 
1073  *                                            << 
1074  ******************************************** << 
1075 {                                             << 
1076   static const G4double kCarTolerance = 1.e-9 << 
1077                                               << 
1078   // check orientation of Triangle            << 
1079   G4double Ax = contour[V[a]].x(), Ay = conto << 
1080   G4double Bx = contour[V[b]].x(), By = conto << 
1081   G4double Cx = contour[V[c]].x(), Cy = conto << 
1082   if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) < kCa << 
1083                                               << 
1084   // check that there is no point inside Tria << 
1085   G4double xmin = std::min(std::min(Ax,Bx),Cx << 
1086   G4double xmax = std::max(std::max(Ax,Bx),Cx << 
1087   G4double ymin = std::min(std::min(Ay,By),Cy << 
1088   G4double ymax = std::max(std::max(Ay,By),Cy << 
1089                                               << 
1090   for (G4int i=0; i<n; ++i)                   << 
1091   {                                           << 
1092     if((i == a) || (i == b) || (i == c)) cont << 
1093     G4double Px = contour[V[i]].x();          << 
1094     if (Px < xmin || Px > xmax) continue;     << 
1095     G4double Py = contour[V[i]].y();          << 
1096     if (Py < ymin || Py > ymax) continue;     << 
1097     // if (PointInTriangle(Ax,Ay,Bx,By,Cx,Cy, << 
1098     if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0 << 
1099     {                                         << 
1100       if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < << 
1101       if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < << 
1102       if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < << 
1103     }                                         << 
1104     else                                      << 
1105     {                                         << 
1106       if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > << 
1107       if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > << 
1108       if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > << 
1109     }                                         << 
1110     return false;                             << 
1111   }                                              606   }
1112   return true;                                << 
1113 }                                                607 }
1114                                                  608 
1115 void HepPolyhedron::SetReferences()              609 void HepPolyhedron::SetReferences()
1116 /********************************************    610 /***********************************************************************
1117  *                                               611  *                                                                     *
1118  * Name: HepPolyhedron::SetReferences            612  * Name: HepPolyhedron::SetReferences                Date:    04.12.96 *
1119  * Author: E.Chernyaev (IHEP/Protvino)           613  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
1120  *                                               614  *                                                                     *
1121  * Function: For each edge set reference to n    615  * Function: For each edge set reference to neighbouring facet         *
1122  *                                               616  *                                                                     *
1123  ********************************************    617  ***********************************************************************/
1124 {                                                618 {
1125   if (nface <= 0) return;                        619   if (nface <= 0) return;
1126                                                  620 
1127   struct edgeListMember {                        621   struct edgeListMember {
1128     edgeListMember *next;                        622     edgeListMember *next;
1129     G4int v2;                                 << 623     int v2;
1130     G4int iface;                              << 624     int iface;
1131     G4int iedge;                              << 625     int iedge;
1132   } *edgeList, *freeList, **headList;            626   } *edgeList, *freeList, **headList;
1133                                                  627 
1134                                               << 628   
1135   //   A L L O C A T E   A N D   I N I T I A     629   //   A L L O C A T E   A N D   I N I T I A T E   L I S T S
1136                                                  630 
1137   edgeList = new edgeListMember[2*nface];        631   edgeList = new edgeListMember[2*nface];
1138   headList = new edgeListMember*[nvert];         632   headList = new edgeListMember*[nvert];
1139                                               << 633   
1140   G4int i;                                    << 634   int i;
1141   for (i=0; i<nvert; i++) {                      635   for (i=0; i<nvert; i++) {
1142     headList[i] = nullptr;                    << 636     headList[i] = 0;
1143   }                                              637   }
1144   freeList = edgeList;                           638   freeList = edgeList;
1145   for (i=0; i<2*nface-1; i++) {                  639   for (i=0; i<2*nface-1; i++) {
1146     edgeList[i].next = &edgeList[i+1];           640     edgeList[i].next = &edgeList[i+1];
1147   }                                              641   }
1148   edgeList[2*nface-1].next = nullptr;         << 642   edgeList[2*nface-1].next = 0;
1149                                                  643 
1150   //   L O O P   A L O N G   E D G E S           644   //   L O O P   A L O N G   E D G E S
1151                                                  645 
1152   G4int iface, iedge, nedge, i1, i2, k1, k2;  << 646   int iface, iedge, nedge, i1, i2, k1, k2;
1153   edgeListMember *prev, *cur;                    647   edgeListMember *prev, *cur;
1154                                               << 648   
1155   for(iface=1; iface<=nface; iface++) {          649   for(iface=1; iface<=nface; iface++) {
1156     nedge = (pF[iface].edge[3].v == 0) ? 3 :     650     nedge = (pF[iface].edge[3].v == 0) ? 3 : 4;
1157     for (iedge=0; iedge<nedge; iedge++) {        651     for (iedge=0; iedge<nedge; iedge++) {
1158       i1 = iedge;                                652       i1 = iedge;
1159       i2 = (iedge < nedge-1) ? iedge+1 : 0;      653       i2 = (iedge < nedge-1) ? iedge+1 : 0;
1160       i1 = std::abs(pF[iface].edge[i1].v);    << 654       i1 = abs(pF[iface].edge[i1].v);
1161       i2 = std::abs(pF[iface].edge[i2].v);    << 655       i2 = abs(pF[iface].edge[i2].v);
1162       k1 = (i1 < i2) ? i1 : i2;          // k    656       k1 = (i1 < i2) ? i1 : i2;          // k1 = ::min(i1,i2);
1163       k2 = (i1 > i2) ? i1 : i2;          // k    657       k2 = (i1 > i2) ? i1 : i2;          // k2 = ::max(i1,i2);
1164                                               << 658       
1165       // check head of the List corresponding    659       // check head of the List corresponding to k1
1166       cur = headList[k1];                        660       cur = headList[k1];
1167       if (cur == nullptr) {                   << 661       if (cur == 0) {
1168         headList[k1] = freeList;              << 662   headList[k1] = freeList;
1169         if (freeList == nullptr) {            << 663   freeList = freeList->next;
1170           std::cerr                           << 
1171           << "Polyhedron::SetReferences: bad  << 
1172           << std::endl;                       << 
1173           break;                              << 
1174         }                                     << 
1175         freeList = freeList->next;            << 
1176         cur = headList[k1];                      664         cur = headList[k1];
1177         cur->next = nullptr;                  << 665   cur->next = 0;
1178         cur->v2 = k2;                         << 666   cur->v2 = k2;
1179         cur->iface = iface;                   << 667   cur->iface = iface;
1180         cur->iedge = iedge;                   << 668   cur->iedge = iedge;
1181         continue;                                669         continue;
1182       }                                       << 670       } 
1183                                                  671 
1184       if (cur->v2 == k2) {                       672       if (cur->v2 == k2) {
1185         headList[k1] = cur->next;                673         headList[k1] = cur->next;
1186         cur->next = freeList;                 << 674   cur->next = freeList;
1187         freeList = cur;                       << 675         freeList = cur;      
1188         pF[iface].edge[iedge].f = cur->iface;    676         pF[iface].edge[iedge].f = cur->iface;
1189         pF[cur->iface].edge[cur->iedge].f = i << 677   pF[cur->iface].edge[cur->iedge].f = iface;
1190         i1 = (pF[iface].edge[iedge].v < 0) ?     678         i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
1191         i2 = (pF[cur->iface].edge[cur->iedge] << 679   i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
1192         if (i1 != i2) {                          680         if (i1 != i2) {
1193           std::cerr                           << 681     HepStd::cerr
1194             << "Polyhedron::SetReferences: di << 682       << "Polyhedron::SetReferences: different edge visibility "
1195             << iface << "/" << iedge << "/"   << 683       << iface << "/" << iedge << "/"
1196             << pF[iface].edge[iedge].v << " a << 684       << pF[iface].edge[iedge].v << " and "
1197             << cur->iface << "/" << cur->iedg << 685       << cur->iface << "/" << cur->iedge << "/"
1198             << pF[cur->iface].edge[cur->iedge << 686       << pF[cur->iface].edge[cur->iedge].v
1199             << std::endl;                     << 687       << HepStd::endl;
1200         }                                     << 688   }
1201         continue;                             << 689   continue;
1202       }                                          690       }
1203                                                  691 
1204       // check List itself                       692       // check List itself
1205       for (;;) {                                 693       for (;;) {
1206         prev = cur;                           << 694   prev = cur;
1207         cur = prev->next;                     << 695   cur = prev->next;
1208         if (cur == nullptr) {                 << 696   if (cur == 0) {
1209           prev->next = freeList;              << 697     prev->next = freeList;
1210           if (freeList == nullptr) {          << 698     freeList = freeList->next;
1211             std::cerr                         << 699     cur = prev->next;
1212             << "Polyhedron::SetReferences: ba << 700     cur->next = 0;
1213             << std::endl;                     << 701     cur->v2 = k2;
1214             break;                            << 702     cur->iface = iface;
1215           }                                   << 703     cur->iedge = iedge;
1216           freeList = freeList->next;          << 704     break;
1217           cur = prev->next;                   << 705   }
1218           cur->next = nullptr;                << 
1219           cur->v2 = k2;                       << 
1220           cur->iface = iface;                 << 
1221           cur->iedge = iedge;                 << 
1222           break;                              << 
1223         }                                     << 
1224                                                  706 
1225         if (cur->v2 == k2) {                     707         if (cur->v2 == k2) {
1226           prev->next = cur->next;             << 708     prev->next = cur->next;
1227           cur->next = freeList;               << 709     cur->next = freeList;
1228           freeList = cur;                     << 710     freeList = cur;      
1229           pF[iface].edge[iedge].f = cur->ifac << 711     pF[iface].edge[iedge].f = cur->iface;
1230           pF[cur->iface].edge[cur->iedge].f = << 712     pF[cur->iface].edge[cur->iedge].f = iface;
1231           i1 = (pF[iface].edge[iedge].v < 0)  << 713     i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
1232           i2 = (pF[cur->iface].edge[cur->iedg << 714     i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
1233             if (i1 != i2) {                   << 715       if (i1 != i2) {
1234               std::cerr                       << 716         HepStd::cerr
1235                 << "Polyhedron::SetReferences << 717     << "Polyhedron::SetReferences: different edge visibility "
1236                 << iface << "/" << iedge << " << 718     << iface << "/" << iedge << "/"
1237                 << pF[iface].edge[iedge].v << << 719     << pF[iface].edge[iedge].v << " and "
1238                 << cur->iface << "/" << cur-> << 720     << cur->iface << "/" << cur->iedge << "/"
1239                 << pF[cur->iface].edge[cur->i << 721     << pF[cur->iface].edge[cur->iedge].v
1240                 << std::endl;                 << 722     << HepStd::endl;
1241             }                                 << 723       }
1242           break;                              << 724     break;
1243         }                                     << 725   }
1244       }                                          726       }
1245     }                                            727     }
1246   }                                              728   }
1247                                                  729 
1248   //  C H E C K   T H A T   A L L   L I S T S    730   //  C H E C K   T H A T   A L L   L I S T S   A R E   E M P T Y
1249                                                  731 
1250   for (i=0; i<nvert; i++) {                      732   for (i=0; i<nvert; i++) {
1251     if (headList[i] != nullptr) {             << 733     if (headList[i] != 0) {
1252       std::cerr                               << 734       HepStd::cerr
1253         << "Polyhedron::SetReferences: List " << 735   << "Polyhedron::SetReferences: List " << i << " is not empty"
1254         << std::endl;                         << 736   << HepStd::endl;
1255     }                                            737     }
1256   }                                              738   }
1257                                                  739 
1258   //   F R E E   M E M O R Y                     740   //   F R E E   M E M O R Y
1259                                                  741 
1260   delete [] edgeList;                            742   delete [] edgeList;
1261   delete [] headList;                            743   delete [] headList;
1262 }                                                744 }
1263                                                  745 
1264 void HepPolyhedron::JoinCoplanarFacets(G4doub << 
1265 /******************************************** << 
1266  *                                            << 
1267  * Name: HepPolyhedron::JoinCoplanarFacets    << 
1268  * Author: E.Tcherniaev (E.Chernyaev)         << 
1269  *                                            << 
1270  * Function: Join couples of triangular facet << 
1271  *           where it is possible             << 
1272  *                                            << 
1273  ******************************************** << 
1274 {                                             << 
1275   G4int njoin = 0;                            << 
1276   for (G4int icur = 1; icur <= nface; ++icur) << 
1277   {                                           << 
1278     // skip if already joined or quadrangle   << 
1279     if (pF[icur].edge[0].v == 0) continue;    << 
1280     if (pF[icur].edge[3].v != 0) continue;    << 
1281     // skip if all references point to alread << 
1282     if (pF[icur].edge[0].f < icur &&          << 
1283         pF[icur].edge[1].f < icur &&          << 
1284         pF[icur].edge[2].f < icur) continue;  << 
1285     // compute plane equation                 << 
1286     G4Normal3D norm = GetUnitNormal(icur);    << 
1287     G4double dd = norm.dot(pV[pF[icur].edge[0 << 
1288     G4int vcur0 = std::abs(pF[icur].edge[0].v << 
1289     G4int vcur1 = std::abs(pF[icur].edge[1].v << 
1290     G4int vcur2 = std::abs(pF[icur].edge[2].v << 
1291     // select neighbouring facet              << 
1292     G4int kcheck = 0, icheck = 0, vcheck = 0; << 
1293     G4double dist = DBL_MAX;                  << 
1294     for (G4int k = 0; k < 3; ++k)             << 
1295     {                                         << 
1296       G4int itmp = pF[icur].edge[k].f;        << 
1297       // skip if already checked, joined or q << 
1298       if (itmp < icur) continue;              << 
1299       if (pF[itmp].edge[0].v == 0 ||          << 
1300           pF[itmp].edge[3].v != 0) continue;  << 
1301       // get candidate vertex                 << 
1302       G4int vtmp = 0;                         << 
1303       for (G4int j = 0; j < 3; ++j)           << 
1304       {                                       << 
1305         vtmp = std::abs(pF[itmp].edge[j].v);  << 
1306   if (vtmp != vcur0 && vtmp != vcur1 && vtmp  << 
1307       }                                       << 
1308       // check distance to the plane          << 
1309       G4double dtmp = std::abs(norm.dot(pV[vt << 
1310       if (dtmp > tolerance || dtmp >= dist) c << 
1311       dist = dtmp;                            << 
1312       kcheck = k;                             << 
1313       icheck = itmp;                          << 
1314       vcheck = vtmp;                          << 
1315     }                                         << 
1316     if (icheck == 0) continue; // no facet se << 
1317     // join facets                            << 
1318     njoin++;                                  << 
1319     pF[icheck].edge[0].v = 0; // mark facet a << 
1320     if (kcheck == 0)                          << 
1321     {                                         << 
1322       pF[icur].edge[3].v = pF[icur].edge[2].v << 
1323       pF[icur].edge[2].v = pF[icur].edge[1].v << 
1324       pF[icur].edge[1].v = vcheck;            << 
1325     }                                         << 
1326     else if (kcheck == 1)                     << 
1327     {                                         << 
1328       pF[icur].edge[3].v = pF[icur].edge[2].v << 
1329       pF[icur].edge[2].v = vcheck;            << 
1330     }                                         << 
1331     else                                      << 
1332     {                                         << 
1333       pF[icur].edge[3].v = vcheck;            << 
1334     }                                         << 
1335   }                                           << 
1336   if (njoin == 0) return; // no joined facets << 
1337                                               << 
1338   // restructure facets                       << 
1339   G4int nnew = 0;                             << 
1340   for (G4int icur = 1; icur <= nface; ++icur) << 
1341   {                                           << 
1342     if (pF[icur].edge[0].v == 0) continue;    << 
1343     nnew++;                                   << 
1344     pF[nnew].edge[0].v = pF[icur].edge[0].v;  << 
1345     pF[nnew].edge[1].v = pF[icur].edge[1].v;  << 
1346     pF[nnew].edge[2].v = pF[icur].edge[2].v;  << 
1347     pF[nnew].edge[3].v = pF[icur].edge[3].v;  << 
1348   }                                           << 
1349   nface = nnew;                               << 
1350   SetReferences();                            << 
1351 }                                             << 
1352                                               << 
1353 void HepPolyhedron::InvertFacets()               746 void HepPolyhedron::InvertFacets()
1354 /********************************************    747 /***********************************************************************
1355  *                                               748  *                                                                     *
1356  * Name: HepPolyhedron::InvertFacets             749  * Name: HepPolyhedron::InvertFacets                Date:    01.12.99  *
1357  * Author: E.Chernyaev                           750  * Author: E.Chernyaev                              Revised:           *
1358  *                                               751  *                                                                     *
1359  * Function: Invert the order of the nodes in    752  * Function: Invert the order of the nodes in the facets               *
1360  *                                               753  *                                                                     *
1361  ********************************************    754  ***********************************************************************/
1362 {                                                755 {
1363   if (nface <= 0) return;                        756   if (nface <= 0) return;
1364   G4int i, k, nnode, v[4],f[4];               << 757   int i, k, nnode, v[4],f[4];
1365   for (i=1; i<=nface; i++) {                     758   for (i=1; i<=nface; i++) {
1366     nnode =  (pF[i].edge[3].v == 0) ? 3 : 4;     759     nnode =  (pF[i].edge[3].v == 0) ? 3 : 4;
1367     for (k=0; k<nnode; k++) {                    760     for (k=0; k<nnode; k++) {
1368       v[k] = (k+1 == nnode) ? pF[i].edge[0].v    761       v[k] = (k+1 == nnode) ? pF[i].edge[0].v : pF[i].edge[k+1].v;
1369       if (v[k] * pF[i].edge[k].v < 0) v[k] =     762       if (v[k] * pF[i].edge[k].v < 0) v[k] = -v[k];
1370       f[k] = pF[i].edge[k].f;                    763       f[k] = pF[i].edge[k].f;
1371     }                                            764     }
1372     for (k=0; k<nnode; k++) {                    765     for (k=0; k<nnode; k++) {
1373       pF[i].edge[nnode-1-k].v = v[k];            766       pF[i].edge[nnode-1-k].v = v[k];
1374       pF[i].edge[nnode-1-k].f = f[k];            767       pF[i].edge[nnode-1-k].f = f[k];
1375     }                                            768     }
1376   }                                              769   }
1377 }                                                770 }
1378                                                  771 
1379 HepPolyhedron & HepPolyhedron::Transform(cons << 772 HepPolyhedron & HepPolyhedron::Transform(const HepTransform3D &t)
1380 /********************************************    773 /***********************************************************************
1381  *                                               774  *                                                                     *
1382  * Name: HepPolyhedron::Transform                775  * Name: HepPolyhedron::Transform                    Date:    01.12.99  *
1383  * Author: E.Chernyaev                           776  * Author: E.Chernyaev                              Revised:           *
1384  *                                               777  *                                                                     *
1385  * Function: Make transformation of the polyh    778  * Function: Make transformation of the polyhedron                     *
1386  *                                               779  *                                                                     *
1387  ********************************************    780  ***********************************************************************/
1388 {                                                781 {
1389   if (nvert > 0) {                               782   if (nvert > 0) {
1390     for (G4int i=1; i<=nvert; i++) { pV[i] =  << 783     for (int i=1; i<=nvert; i++) { pV[i] = t * pV[i]; }
1391                                                  784 
1392     //  C H E C K   D E T E R M I N A N T   A    785     //  C H E C K   D E T E R M I N A N T   A N D
1393     //  I N V E R T   F A C E T S   I F   I T    786     //  I N V E R T   F A C E T S   I F   I T   I S   N E G A T I V E
1394                                                  787 
1395     G4Vector3D d = t * G4Vector3D(0,0,0);     << 788     HepVector3D d = t * HepVector3D(0,0,0);
1396     G4Vector3D x = t * G4Vector3D(1,0,0) - d; << 789     HepVector3D x = t * HepVector3D(1,0,0) - d;
1397     G4Vector3D y = t * G4Vector3D(0,1,0) - d; << 790     HepVector3D y = t * HepVector3D(0,1,0) - d;
1398     G4Vector3D z = t * G4Vector3D(0,0,1) - d; << 791     HepVector3D z = t * HepVector3D(0,0,1) - d;
1399     if ((x.cross(y))*z < 0) InvertFacets();      792     if ((x.cross(y))*z < 0) InvertFacets();
1400   }                                              793   }
1401   return *this;                                  794   return *this;
1402 }                                                795 }
1403                                                  796 
1404 G4bool HepPolyhedron::GetNextVertexIndex(G4in << 797 HepBoolean HepPolyhedron::GetNextVertexIndex(int &index, int &edgeFlag) const
1405 /********************************************    798 /***********************************************************************
1406  *                                               799  *                                                                     *
1407  * Name: HepPolyhedron::GetNextVertexIndex       800  * Name: HepPolyhedron::GetNextVertexIndex          Date:    03.09.96  *
1408  * Author: Yasuhide Sawada                       801  * Author: Yasuhide Sawada                          Revised:           *
1409  *                                               802  *                                                                     *
1410  * Function:                                     803  * Function:                                                           *
1411  *                                               804  *                                                                     *
1412  ********************************************    805  ***********************************************************************/
1413 {                                                806 {
1414   static G4ThreadLocal G4int iFace = 1;       << 807   static int iFace = 1;
1415   static G4ThreadLocal G4int iQVertex = 0;    << 808   static int iQVertex = 0;
1416   G4int vIndex = pF[iFace].edge[iQVertex].v;  << 809   int vIndex = pF[iFace].edge[iQVertex].v;
1417                                                  810 
1418   edgeFlag = (vIndex > 0) ? 1 : 0;               811   edgeFlag = (vIndex > 0) ? 1 : 0;
1419   index = std::abs(vIndex);                   << 812   index = abs(vIndex);
1420                                                  813 
1421   if (iQVertex >= 3 || pF[iFace].edge[iQVerte    814   if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1422     iQVertex = 0;                                815     iQVertex = 0;
1423     if (++iFace > nface) iFace = 1;              816     if (++iFace > nface) iFace = 1;
1424     return false;  // Last Edge                  817     return false;  // Last Edge
                                                   >> 818   }else{
                                                   >> 819     ++iQVertex;
                                                   >> 820     return true;  // not Last Edge
1425   }                                              821   }
1426                                               << 
1427   ++iQVertex;                                 << 
1428   return true;  // not Last Edge              << 
1429 }                                                822 }
1430                                                  823 
1431 G4Point3D HepPolyhedron::GetVertex(G4int inde << 824 HepPoint3D HepPolyhedron::GetVertex(int index) const
1432 /********************************************    825 /***********************************************************************
1433  *                                               826  *                                                                     *
1434  * Name: HepPolyhedron::GetVertex                827  * Name: HepPolyhedron::GetVertex                   Date:    03.09.96  *
1435  * Author: Yasuhide Sawada                       828  * Author: Yasuhide Sawada                          Revised: 17.11.99  *
1436  *                                               829  *                                                                     *
1437  * Function: Get vertex of the index.            830  * Function: Get vertex of the index.                                  *
1438  *                                               831  *                                                                     *
1439  ********************************************    832  ***********************************************************************/
1440 {                                                833 {
1441   if (index <= 0 || index > nvert) {             834   if (index <= 0 || index > nvert) {
1442     std::cerr                                 << 835     HepStd::cerr
1443       << "HepPolyhedron::GetVertex: irrelevan    836       << "HepPolyhedron::GetVertex: irrelevant index " << index
1444       << std::endl;                           << 837       << HepStd::endl;
1445     return G4Point3D();                       << 838     return HepPoint3D();
1446   }                                              839   }
1447   return pV[index];                              840   return pV[index];
1448 }                                                841 }
1449                                                  842 
1450 G4bool                                        << 843 HepBoolean
1451 HepPolyhedron::GetNextVertex(G4Point3D &verte << 844 HepPolyhedron::GetNextVertex(HepPoint3D &vertex, int &edgeFlag) const
1452 /********************************************    845 /***********************************************************************
1453  *                                               846  *                                                                     *
1454  * Name: HepPolyhedron::GetNextVertex            847  * Name: HepPolyhedron::GetNextVertex               Date:    22.07.96  *
1455  * Author: John Allison                          848  * Author: John Allison                             Revised:           *
1456  *                                               849  *                                                                     *
1457  * Function: Get vertices of the quadrilatera    850  * Function: Get vertices of the quadrilaterals in order for each      *
1458  *           face in face order.  Returns fal    851  *           face in face order.  Returns false when finished each     *
1459  *           face.                               852  *           face.                                                     *
1460  *                                               853  *                                                                     *
1461  ********************************************    854  ***********************************************************************/
1462 {                                                855 {
1463   G4int index;                                << 856   int index;
1464   G4bool rep = GetNextVertexIndex(index, edge << 857   HepBoolean rep = GetNextVertexIndex(index, edgeFlag);
1465   vertex = pV[index];                            858   vertex = pV[index];
1466   return rep;                                    859   return rep;
1467 }                                                860 }
1468                                                  861 
1469 G4bool HepPolyhedron::GetNextVertex(G4Point3D << 862 HepBoolean HepPolyhedron::GetNextVertex(HepPoint3D &vertex, int &edgeFlag,
1470                                   G4Normal3D  << 863                HepNormal3D &normal) const
1471 /********************************************    864 /***********************************************************************
1472  *                                               865  *                                                                     *
1473  * Name: HepPolyhedron::GetNextVertex            866  * Name: HepPolyhedron::GetNextVertex               Date:    26.11.99  *
1474  * Author: E.Chernyaev                           867  * Author: E.Chernyaev                              Revised:           *
1475  *                                               868  *                                                                     *
1476  * Function: Get vertices with normals of the    869  * Function: Get vertices with normals of the quadrilaterals in order  *
1477  *           for each face in face order.        870  *           for each face in face order.                              *
1478  *           Returns false when finished each    871  *           Returns false when finished each face.                    *
1479  *                                               872  *                                                                     *
1480  ********************************************    873  ***********************************************************************/
1481 {                                                874 {
1482   static G4ThreadLocal G4int iFace = 1;       << 875   static int iFace = 1;
1483   static G4ThreadLocal G4int iNode = 0;       << 876   static int iNode = 0;
1484                                                  877 
1485   if (nface == 0) return false;  // empty pol    878   if (nface == 0) return false;  // empty polyhedron
1486                                                  879 
1487   G4int k = pF[iFace].edge[iNode].v;          << 880   int k = pF[iFace].edge[iNode].v;
1488   if (k > 0) { edgeFlag = 1; } else { edgeFla    881   if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
1489   vertex = pV[k];                                882   vertex = pV[k];
1490   normal = FindNodeNormal(iFace,k);              883   normal = FindNodeNormal(iFace,k);
1491   if (iNode >= 3 || pF[iFace].edge[iNode+1].v    884   if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
1492     iNode = 0;                                   885     iNode = 0;
1493     if (++iFace > nface) iFace = 1;              886     if (++iFace > nface) iFace = 1;
1494     return false;                // last node    887     return false;                // last node
                                                   >> 888   }else{
                                                   >> 889     ++iNode;
                                                   >> 890     return true;                 // not last node
1495   }                                              891   }
1496   ++iNode;                                    << 
1497   return true;                 // not last no << 
1498 }                                                892 }
1499                                                  893 
1500 G4bool HepPolyhedron::GetNextEdgeIndices(G4in << 894 HepBoolean HepPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag,
1501                                        G4int  << 895               int &iface1, int &iface2) const
1502 /********************************************    896 /***********************************************************************
1503  *                                               897  *                                                                     *
1504  * Name: HepPolyhedron::GetNextEdgeIndices    << 898  * Name: HepPolyhedron::GetNextEdgeIndeces          Date:    30.09.96  *
1505  * Author: E.Chernyaev                           899  * Author: E.Chernyaev                              Revised: 17.11.99  *
1506  *                                               900  *                                                                     *
1507  * Function: Get indices of the next edge tog << 901  * Function: Get indeces of the next edge together with indeces of     *
1508  *           of the faces which share the edg    902  *           of the faces which share the edge.                        *
1509  *           Returns false when the last edge    903  *           Returns false when the last edge.                         *
1510  *                                               904  *                                                                     *
1511  ********************************************    905  ***********************************************************************/
1512 {                                                906 {
1513   static G4ThreadLocal G4int iFace    = 1;    << 907   static int iFace    = 1;
1514   static G4ThreadLocal G4int iQVertex = 0;    << 908   static int iQVertex = 0;
1515   static G4ThreadLocal G4int iOrder   = 1;    << 909   static int iOrder   = 1;
1516   G4int  k1, k2, kflag, kface1, kface2;       << 910   int  k1, k2, kflag, kface1, kface2;
1517                                                  911 
1518   if (iFace == 1 && iQVertex == 0) {             912   if (iFace == 1 && iQVertex == 0) {
1519     k2 = pF[nface].edge[0].v;                    913     k2 = pF[nface].edge[0].v;
1520     k1 = pF[nface].edge[3].v;                    914     k1 = pF[nface].edge[3].v;
1521     if (k1 == 0) k1 = pF[nface].edge[2].v;       915     if (k1 == 0) k1 = pF[nface].edge[2].v;
1522     if (std::abs(k1) > std::abs(k2)) iOrder = << 916     if (abs(k1) > abs(k2)) iOrder = -1;
1523   }                                              917   }
1524                                                  918 
1525   do {                                           919   do {
1526     k1     = pF[iFace].edge[iQVertex].v;         920     k1     = pF[iFace].edge[iQVertex].v;
1527     kflag  = k1;                                 921     kflag  = k1;
1528     k1     = std::abs(k1);                    << 922     k1     = abs(k1);
1529     kface1 = iFace;                           << 923     kface1 = iFace; 
1530     kface2 = pF[iFace].edge[iQVertex].f;         924     kface2 = pF[iFace].edge[iQVertex].f;
1531     if (iQVertex >= 3 || pF[iFace].edge[iQVer    925     if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1532       iQVertex = 0;                              926       iQVertex = 0;
1533       k2 = std::abs(pF[iFace].edge[iQVertex]. << 927       k2 = abs(pF[iFace].edge[iQVertex].v);
1534       iFace++;                                   928       iFace++;
1535     }else{                                       929     }else{
1536       iQVertex++;                                930       iQVertex++;
1537       k2 = std::abs(pF[iFace].edge[iQVertex]. << 931       k2 = abs(pF[iFace].edge[iQVertex].v);
1538     }                                            932     }
1539   } while (iOrder*k1 > iOrder*k2);               933   } while (iOrder*k1 > iOrder*k2);
1540                                                  934 
1541   i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ?     935   i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1542   iface1 = kface1; iface2 = kface2;           << 936   iface1 = kface1; iface2 = kface2; 
1543                                                  937 
1544   if (iFace > nface) {                           938   if (iFace > nface) {
1545     iFace  = 1; iOrder = 1;                      939     iFace  = 1; iOrder = 1;
1546     return false;                                940     return false;
                                                   >> 941   }else{
                                                   >> 942     return true;
1547   }                                              943   }
1548                                               << 
1549   return true;                                << 
1550 }                                                944 }
1551                                                  945 
1552 G4bool                                        << 946 HepBoolean
1553 HepPolyhedron::GetNextEdgeIndices(G4int &i1,  << 947 HepPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag) const
1554 /********************************************    948 /***********************************************************************
1555  *                                               949  *                                                                     *
1556  * Name: HepPolyhedron::GetNextEdgeIndices    << 950  * Name: HepPolyhedron::GetNextEdgeIndeces          Date:    17.11.99  *
1557  * Author: E.Chernyaev                           951  * Author: E.Chernyaev                              Revised:           *
1558  *                                               952  *                                                                     *
1559  * Function: Get indices of the next edge.    << 953  * Function: Get indeces of the next edge.                             *
1560  *           Returns false when the last edge    954  *           Returns false when the last edge.                         *
1561  *                                               955  *                                                                     *
1562  ********************************************    956  ***********************************************************************/
1563 {                                                957 {
1564   G4int kface1, kface2;                       << 958   int kface1, kface2;
1565   return GetNextEdgeIndices(i1, i2, edgeFlag, << 959   return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
1566 }                                                960 }
1567                                                  961 
1568 G4bool                                        << 962 HepBoolean
1569 HepPolyhedron::GetNextEdge(G4Point3D &p1,     << 963 HepPolyhedron::GetNextEdge(HepPoint3D &p1,
1570                            G4Point3D &p2,     << 964          HepPoint3D &p2,
1571                            G4int &edgeFlag) c << 965          int &edgeFlag) const
1572 /********************************************    966 /***********************************************************************
1573  *                                               967  *                                                                     *
1574  * Name: HepPolyhedron::GetNextEdge              968  * Name: HepPolyhedron::GetNextEdge                 Date:    30.09.96  *
1575  * Author: E.Chernyaev                           969  * Author: E.Chernyaev                              Revised:           *
1576  *                                               970  *                                                                     *
1577  * Function: Get next edge.                      971  * Function: Get next edge.                                            *
1578  *           Returns false when the last edge    972  *           Returns false when the last edge.                         *
1579  *                                               973  *                                                                     *
1580  ********************************************    974  ***********************************************************************/
1581 {                                                975 {
1582   G4int i1,i2;                                << 976   int i1,i2;
1583   G4bool rep = GetNextEdgeIndices(i1,i2,edgeF << 977   HepBoolean rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
1584   p1 = pV[i1];                                   978   p1 = pV[i1];
1585   p2 = pV[i2];                                   979   p2 = pV[i2];
1586   return rep;                                    980   return rep;
1587 }                                                981 }
1588                                                  982 
1589 G4bool                                        << 983 HepBoolean
1590 HepPolyhedron::GetNextEdge(G4Point3D &p1, G4P << 984 HepPolyhedron::GetNextEdge(HepPoint3D &p1, HepPoint3D &p2,
1591                           G4int &edgeFlag, G4 << 985         int &edgeFlag, int &iface1, int &iface2) const
1592 /********************************************    986 /***********************************************************************
1593  *                                               987  *                                                                     *
1594  * Name: HepPolyhedron::GetNextEdge              988  * Name: HepPolyhedron::GetNextEdge                 Date:    17.11.99  *
1595  * Author: E.Chernyaev                           989  * Author: E.Chernyaev                              Revised:           *
1596  *                                               990  *                                                                     *
1597  * Function: Get next edge with indices of th << 991  * Function: Get next edge with indeces of the faces which share       *
1598  *           the edge.                           992  *           the edge.                                                 *
1599  *           Returns false when the last edge    993  *           Returns false when the last edge.                         *
1600  *                                               994  *                                                                     *
1601  ********************************************    995  ***********************************************************************/
1602 {                                                996 {
1603   G4int i1,i2;                                << 997   int i1,i2;
1604   G4bool rep = GetNextEdgeIndices(i1,i2,edgeF << 998   HepBoolean rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
1605   p1 = pV[i1];                                   999   p1 = pV[i1];
1606   p2 = pV[i2];                                   1000   p2 = pV[i2];
1607   return rep;                                    1001   return rep;
1608 }                                                1002 }
1609                                                  1003 
1610 void HepPolyhedron::GetFacet(G4int iFace, G4i << 1004 void HepPolyhedron::GetFacet(int iFace, int &n, int *iNodes,
1611                             G4int *edgeFlags, << 1005           int *edgeFlags, int *iFaces) const
1612 /********************************************    1006 /***********************************************************************
1613  *                                               1007  *                                                                     *
1614  * Name: HepPolyhedron::GetFacet                 1008  * Name: HepPolyhedron::GetFacet                    Date:    15.12.99  *
1615  * Author: E.Chernyaev                           1009  * Author: E.Chernyaev                              Revised:           *
1616  *                                               1010  *                                                                     *
1617  * Function: Get face by index                   1011  * Function: Get face by index                                         *
1618  *                                               1012  *                                                                     *
1619  ********************************************    1013  ***********************************************************************/
1620 {                                                1014 {
1621   if (iFace < 1 || iFace > nface) {              1015   if (iFace < 1 || iFace > nface) {
1622     std::cerr                                 << 1016     HepStd::cerr 
1623       << "HepPolyhedron::GetFacet: irrelevant    1017       << "HepPolyhedron::GetFacet: irrelevant index " << iFace
1624       << std::endl;                           << 1018       << HepStd::endl;
1625     n = 0;                                       1019     n = 0;
1626   }else{                                         1020   }else{
1627     G4int i, k;                               << 1021     int i, k;
1628     for (i=0; i<4; i++) {                     << 1022     for (i=0; i<4; i++) { 
1629       k = pF[iFace].edge[i].v;                   1023       k = pF[iFace].edge[i].v;
1630       if (k == 0) break;                         1024       if (k == 0) break;
1631       if (iFaces != nullptr) iFaces[i] = pF[i << 1025       if (iFaces != 0) iFaces[i] = pF[iFace].edge[i].f;
1632       if (k > 0) {                            << 1026       if (k > 0) { 
1633         iNodes[i] = k;                        << 1027   iNodes[i] = k;
1634         if (edgeFlags != nullptr) edgeFlags[i << 1028   if (edgeFlags != 0) edgeFlags[i] = 1;
1635       }else{                                     1029       }else{
1636         iNodes[i] = -k;                       << 1030   iNodes[i] = -k;
1637         if (edgeFlags != nullptr) edgeFlags[i << 1031   if (edgeFlags != 0) edgeFlags[i] = -1;
1638       }                                          1032       }
1639     }                                            1033     }
1640     n = i;                                       1034     n = i;
1641   }                                              1035   }
1642 }                                                1036 }
1643                                                  1037 
1644 void HepPolyhedron::GetFacet(G4int index, G4i << 1038 void HepPolyhedron::GetFacet(int index, int &n, HepPoint3D *nodes,
1645                              G4int *edgeFlags << 1039           int *edgeFlags, HepNormal3D *normals) const
1646 /********************************************    1040 /***********************************************************************
1647  *                                               1041  *                                                                     *
1648  * Name: HepPolyhedron::GetFacet                 1042  * Name: HepPolyhedron::GetFacet                    Date:    17.11.99  *
1649  * Author: E.Chernyaev                           1043  * Author: E.Chernyaev                              Revised:           *
1650  *                                               1044  *                                                                     *
1651  * Function: Get face by index                   1045  * Function: Get face by index                                         *
1652  *                                               1046  *                                                                     *
1653  ********************************************    1047  ***********************************************************************/
1654 {                                                1048 {
1655   G4int iNodes[4];                            << 1049   int iNodes[4];
1656   GetFacet(index, n, iNodes, edgeFlags);         1050   GetFacet(index, n, iNodes, edgeFlags);
1657   if (n != 0) {                                  1051   if (n != 0) {
1658     for (G4int i=0; i<n; i++) {               << 1052     for (int i=0; i<4; i++) {
1659       nodes[i] = pV[iNodes[i]];                  1053       nodes[i] = pV[iNodes[i]];
1660       if (normals != nullptr) normals[i] = Fi << 1054       if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1661     }                                            1055     }
1662   }                                              1056   }
1663 }                                                1057 }
1664                                                  1058 
1665 G4bool                                        << 1059 HepBoolean
1666 HepPolyhedron::GetNextFacet(G4int &n, G4Point << 1060 HepPolyhedron::GetNextFacet(int &n, HepPoint3D *nodes,
1667                            G4int *edgeFlags,  << 1061          int *edgeFlags, HepNormal3D *normals) const
1668 /********************************************    1062 /***********************************************************************
1669  *                                               1063  *                                                                     *
1670  * Name: HepPolyhedron::GetNextFacet             1064  * Name: HepPolyhedron::GetNextFacet                Date:    19.11.99  *
1671  * Author: E.Chernyaev                           1065  * Author: E.Chernyaev                              Revised:           *
1672  *                                               1066  *                                                                     *
1673  * Function: Get next face with normals of un    1067  * Function: Get next face with normals of unit length at the nodes.   *
1674  *           Returns false when finished all     1068  *           Returns false when finished all faces.                    *
1675  *                                               1069  *                                                                     *
1676  ********************************************    1070  ***********************************************************************/
1677 {                                                1071 {
1678   static G4ThreadLocal G4int iFace = 1;       << 1072   static int iFace = 1;
1679                                                  1073 
1680   if (edgeFlags == nullptr) {                 << 1074   if (edgeFlags == 0) {
1681     GetFacet(iFace, n, nodes);                   1075     GetFacet(iFace, n, nodes);
1682   }else if (normals == nullptr) {             << 1076   }else if (normals == 0) {
1683     GetFacet(iFace, n, nodes, edgeFlags);        1077     GetFacet(iFace, n, nodes, edgeFlags);
1684   }else{                                         1078   }else{
1685     GetFacet(iFace, n, nodes, edgeFlags, norm    1079     GetFacet(iFace, n, nodes, edgeFlags, normals);
1686   }                                              1080   }
1687                                                  1081 
1688   if (++iFace > nface) {                         1082   if (++iFace > nface) {
1689     iFace  = 1;                                  1083     iFace  = 1;
1690     return false;                                1084     return false;
                                                   >> 1085   }else{
                                                   >> 1086     return true;
1691   }                                              1087   }
1692                                               << 
1693   return true;                                << 
1694 }                                                1088 }
1695                                                  1089 
1696 G4Normal3D HepPolyhedron::GetNormal(G4int iFa << 1090 HepNormal3D HepPolyhedron::GetNormal(int iFace) const
1697 /********************************************    1091 /***********************************************************************
1698  *                                               1092  *                                                                     *
1699  * Name: HepPolyhedron::GetNormal                1093  * Name: HepPolyhedron::GetNormal                    Date:    19.11.99 *
1700  * Author: E.Chernyaev                           1094  * Author: E.Chernyaev                               Revised:          *
1701  *                                               1095  *                                                                     *
1702  * Function: Get normal of the face given by     1096  * Function: Get normal of the face given by index                     *
1703  *                                               1097  *                                                                     *
1704  ********************************************    1098  ***********************************************************************/
1705 {                                                1099 {
1706   if (iFace < 1 || iFace > nface) {              1100   if (iFace < 1 || iFace > nface) {
1707     std::cerr                                 << 1101     HepStd::cerr 
1708       << "HepPolyhedron::GetNormal: irrelevan << 1102       << "HepPolyhedron::GetNormal: irrelevant index " << iFace 
1709       << std::endl;                           << 1103       << HepStd::endl;
1710     return G4Normal3D();                      << 1104     return HepNormal3D();
1711   }                                              1105   }
1712                                                  1106 
1713   G4int i0  = std::abs(pF[iFace].edge[0].v);  << 1107   int i0  = abs(pF[iFace].edge[0].v);
1714   G4int i1  = std::abs(pF[iFace].edge[1].v);  << 1108   int i1  = abs(pF[iFace].edge[1].v);
1715   G4int i2  = std::abs(pF[iFace].edge[2].v);  << 1109   int i2  = abs(pF[iFace].edge[2].v);
1716   G4int i3  = std::abs(pF[iFace].edge[3].v);  << 1110   int i3  = abs(pF[iFace].edge[3].v);
1717   if (i3 == 0) i3 = i0;                          1111   if (i3 == 0) i3 = i0;
1718   return (pV[i2] - pV[i0]).cross(pV[i3] - pV[    1112   return (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1]);
1719 }                                                1113 }
1720                                                  1114 
1721 G4Normal3D HepPolyhedron::GetUnitNormal(G4int << 1115 HepNormal3D HepPolyhedron::GetUnitNormal(int iFace) const
1722 /********************************************    1116 /***********************************************************************
1723  *                                               1117  *                                                                     *
1724  * Name: HepPolyhedron::GetNormal                1118  * Name: HepPolyhedron::GetNormal                    Date:    19.11.99 *
1725  * Author: E.Chernyaev                           1119  * Author: E.Chernyaev                               Revised:          *
1726  *                                               1120  *                                                                     *
1727  * Function: Get unit normal of the face give    1121  * Function: Get unit normal of the face given by index                *
1728  *                                               1122  *                                                                     *
1729  ********************************************    1123  ***********************************************************************/
1730 {                                                1124 {
1731   if (iFace < 1 || iFace > nface) {              1125   if (iFace < 1 || iFace > nface) {
1732     std::cerr                                 << 1126     HepStd::cerr 
1733       << "HepPolyhedron::GetUnitNormal: irrel    1127       << "HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
1734       << std::endl;                           << 1128       << HepStd::endl;
1735     return G4Normal3D();                      << 1129     return HepNormal3D();
1736   }                                              1130   }
1737                                                  1131 
1738   G4int i0  = std::abs(pF[iFace].edge[0].v);  << 1132   int i0  = abs(pF[iFace].edge[0].v);
1739   G4int i1  = std::abs(pF[iFace].edge[1].v);  << 1133   int i1  = abs(pF[iFace].edge[1].v);
1740   G4int i2  = std::abs(pF[iFace].edge[2].v);  << 1134   int i2  = abs(pF[iFace].edge[2].v);
1741   G4int i3  = std::abs(pF[iFace].edge[3].v);  << 1135   int i3  = abs(pF[iFace].edge[3].v);
1742   if (i3 == 0) i3 = i0;                          1136   if (i3 == 0) i3 = i0;
1743   return ((pV[i2] - pV[i0]).cross(pV[i3] - pV    1137   return ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).unit();
1744 }                                                1138 }
1745                                                  1139 
1746 G4bool HepPolyhedron::GetNextNormal(G4Normal3 << 1140 HepBoolean HepPolyhedron::GetNextNormal(HepNormal3D &normal) const
1747 /********************************************    1141 /***********************************************************************
1748  *                                               1142  *                                                                     *
1749  * Name: HepPolyhedron::GetNextNormal            1143  * Name: HepPolyhedron::GetNextNormal               Date:    22.07.96  *
1750  * Author: John Allison                          1144  * Author: John Allison                             Revised: 19.11.99  *
1751  *                                               1145  *                                                                     *
1752  * Function: Get normals of each face in face    1146  * Function: Get normals of each face in face order.  Returns false    *
1753  *           when finished all faces.            1147  *           when finished all faces.                                  *
1754  *                                               1148  *                                                                     *
1755  ********************************************    1149  ***********************************************************************/
1756 {                                                1150 {
1757   static G4ThreadLocal G4int iFace = 1;       << 1151   static int iFace = 1;
1758   normal = GetNormal(iFace);                     1152   normal = GetNormal(iFace);
1759   if (++iFace > nface) {                         1153   if (++iFace > nface) {
1760     iFace = 1;                                   1154     iFace = 1;
1761     return false;                                1155     return false;
                                                   >> 1156   }else{
                                                   >> 1157     return true;
1762   }                                              1158   }
1763   return true;                                << 
1764 }                                                1159 }
1765                                                  1160 
1766 G4bool HepPolyhedron::GetNextUnitNormal(G4Nor << 1161 HepBoolean HepPolyhedron::GetNextUnitNormal(HepNormal3D &normal) const
1767 /********************************************    1162 /***********************************************************************
1768  *                                               1163  *                                                                     *
1769  * Name: HepPolyhedron::GetNextUnitNormal        1164  * Name: HepPolyhedron::GetNextUnitNormal           Date:    16.09.96  *
1770  * Author: E.Chernyaev                           1165  * Author: E.Chernyaev                              Revised:           *
1771  *                                               1166  *                                                                     *
1772  * Function: Get normals of unit length of ea    1167  * Function: Get normals of unit length of each face in face order.    *
1773  *           Returns false when finished all     1168  *           Returns false when finished all faces.                    *
1774  *                                               1169  *                                                                     *
1775  ********************************************    1170  ***********************************************************************/
1776 {                                                1171 {
1777   G4bool rep = GetNextNormal(normal);         << 1172   HepBoolean rep = GetNextNormal(normal);
1778   normal = normal.unit();                        1173   normal = normal.unit();
1779   return rep;                                    1174   return rep;
1780 }                                                1175 }
1781                                                  1176 
1782 G4double HepPolyhedron::GetSurfaceArea() cons << 1177 HepPolyhedronTrd2::HepPolyhedronTrd2(HepDouble Dx1, HepDouble Dx2,
1783 /******************************************** << 1178              HepDouble Dy1, HepDouble Dy2,
1784  *                                            << 1179              HepDouble Dz)
1785  * Name: HepPolyhedron::GetSurfaceArea        << 
1786  * Author: E.Chernyaev                        << 
1787  *                                            << 
1788  * Function: Returns area of the surface of t << 
1789  *                                            << 
1790  ******************************************** << 
1791 {                                             << 
1792   G4double srf = 0.;                          << 
1793   for (G4int iFace=1; iFace<=nface; iFace++)  << 
1794     G4int i0 = std::abs(pF[iFace].edge[0].v); << 
1795     G4int i1 = std::abs(pF[iFace].edge[1].v); << 
1796     G4int i2 = std::abs(pF[iFace].edge[2].v); << 
1797     G4int i3 = std::abs(pF[iFace].edge[3].v); << 
1798     if (i3 == 0) i3 = i0;                     << 
1799     srf += ((pV[i2] - pV[i0]).cross(pV[i3] -  << 
1800   }                                           << 
1801   return srf/2.;                              << 
1802 }                                             << 
1803                                               << 
1804 G4double HepPolyhedron::GetVolume() const     << 
1805 /******************************************** << 
1806  *                                            << 
1807  * Name: HepPolyhedron::GetVolume             << 
1808  * Author: E.Chernyaev                        << 
1809  *                                            << 
1810  * Function: Returns volume of the polyhedron << 
1811  *                                            << 
1812  ******************************************** << 
1813 {                                             << 
1814   G4double v = 0.;                            << 
1815   for (G4int iFace=1; iFace<=nface; iFace++)  << 
1816     G4int i0 = std::abs(pF[iFace].edge[0].v); << 
1817     G4int i1 = std::abs(pF[iFace].edge[1].v); << 
1818     G4int i2 = std::abs(pF[iFace].edge[2].v); << 
1819     G4int i3 = std::abs(pF[iFace].edge[3].v); << 
1820     G4Point3D pt;                             << 
1821     if (i3 == 0) {                            << 
1822       i3 = i0;                                << 
1823       pt = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);  << 
1824     }else{                                    << 
1825       pt = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0. << 
1826     }                                         << 
1827     v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV << 
1828   }                                           << 
1829   return v/6.;                                << 
1830 }                                             << 
1831                                               << 
1832 G4int                                         << 
1833 HepPolyhedron::createTwistedTrap(G4double Dz, << 
1834                                  const G4doub << 
1835                                  const G4doub << 
1836 /******************************************** << 
1837  *                                            << 
1838  * Name: createTwistedTrap                    << 
1839  * Author: E.Chernyaev (IHEP/Protvino)        << 
1840  *                                            << 
1841  * Function: Creates polyhedron for twisted t << 
1842  *                                            << 
1843  * Input: Dz       - half-length along Z      << 
1844  *        xy1[2,4] - quadrilateral at Z=-Dz   << 
1845  *        xy2[2,4] - quadrilateral at Z=+Dz   << 
1846  *                                            << 
1847  *                                            << 
1848  ******************************************** << 
1849 {                                             << 
1850   AllocateMemory(12,18);                      << 
1851                                               << 
1852   pV[ 1] = G4Point3D(xy1[0][0],xy1[0][1],-Dz) << 
1853   pV[ 2] = G4Point3D(xy1[1][0],xy1[1][1],-Dz) << 
1854   pV[ 3] = G4Point3D(xy1[2][0],xy1[2][1],-Dz) << 
1855   pV[ 4] = G4Point3D(xy1[3][0],xy1[3][1],-Dz) << 
1856                                               << 
1857   pV[ 5] = G4Point3D(xy2[0][0],xy2[0][1], Dz) << 
1858   pV[ 6] = G4Point3D(xy2[1][0],xy2[1][1], Dz) << 
1859   pV[ 7] = G4Point3D(xy2[2][0],xy2[2][1], Dz) << 
1860   pV[ 8] = G4Point3D(xy2[3][0],xy2[3][1], Dz) << 
1861                                               << 
1862   pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;      << 
1863   pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;      << 
1864   pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;      << 
1865   pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;      << 
1866                                               << 
1867   enum {DUMMY, BOTTOM,                        << 
1868         LEFT_BOTTOM,  LEFT_FRONT,   LEFT_TOP, << 
1869         BACK_BOTTOM,  BACK_LEFT,    BACK_TOP, << 
1870         RIGHT_BOTTOM, RIGHT_BACK,   RIGHT_TOP << 
1871         FRONT_BOTTOM, FRONT_RIGHT,  FRONT_TOP << 
1872         TOP};                                 << 
1873                                               << 
1874   pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM << 
1875                                               << 
1876   pF[ 2]=G4Facet(4,BOTTOM,     -1,LEFT_FRONT, << 
1877   pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP,   << 
1878   pF[ 4]=G4Facet(5,TOP,        -8,LEFT_BACK,  << 
1879   pF[ 5]=G4Facet(8,BACK_LEFT,  -4,LEFT_BOTTOM << 
1880                                               << 
1881   pF[ 6]=G4Facet(3,BOTTOM,     -4,BACK_LEFT,  << 
1882   pF[ 7]=G4Facet(4,LEFT_BACK,  -8,BACK_TOP,   << 
1883   pF[ 8]=G4Facet(8,TOP,        -7,BACK_RIGHT, << 
1884   pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM << 
1885                                               << 
1886   pF[10]=G4Facet(2,BOTTOM,     -3,RIGHT_BACK, << 
1887   pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP,  << 
1888   pF[12]=G4Facet(7,TOP,        -6,RIGHT_FRONT << 
1889   pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTO << 
1890                                               << 
1891   pF[14]=G4Facet(1,BOTTOM,     -2,FRONT_RIGHT << 
1892   pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP,  << 
1893   pF[16]=G4Facet(6,TOP,        -5,FRONT_LEFT, << 
1894   pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTO << 
1895                                               << 
1896   pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7, << 
1897                                               << 
1898   return 0;                                   << 
1899 }                                             << 
1900                                               << 
1901 G4int                                         << 
1902 HepPolyhedron::createPolyhedron(G4int Nnodes, << 
1903                                 const G4doubl << 
1904                                 const G4int   << 
1905 /******************************************** << 
1906  *                                            << 
1907  * Name: createPolyhedron                     << 
1908  * Author: E.Chernyaev (IHEP/Protvino)        << 
1909  *                                            << 
1910  * Function: Creates user defined polyhedron  << 
1911  *                                            << 
1912  * Input: Nnodes  - number of nodes           << 
1913  *        Nfaces  - number of faces           << 
1914  *        nodes[][3] - node coordinates       << 
1915  *        faces[][4] - faces                  << 
1916  *                                            << 
1917  ******************************************** << 
1918 {                                             << 
1919   AllocateMemory(Nnodes, Nfaces);             << 
1920   if (nvert == 0) return 1;                   << 
1921                                               << 
1922   for (G4int i=0; i<Nnodes; i++) {            << 
1923     pV[i+1] = G4Point3D(xyz[i][0], xyz[i][1], << 
1924   }                                           << 
1925   for (G4int k=0; k<Nfaces; k++) {            << 
1926     pF[k+1] = G4Facet(faces[k][0],0,faces[k][ << 
1927   }                                           << 
1928   SetReferences();                            << 
1929   return 0;                                   << 
1930 }                                             << 
1931                                               << 
1932 G4Point3D HepPolyhedron::vertexUnweightedMean << 
1933   /****************************************** << 
1934    *                                          << 
1935    * Name: vertexUnweightedMean               << 
1936    * Author: S. Boogert (Manchester)          << 
1937    *                                          << 
1938    * Function: Calculate the unweighted mean  << 
1939    * in the polyhedron. Not to be confused wi << 
1940    * centre of mass                           << 
1941    ****************************************** << 
1942                                               << 
1943   auto centre = G4Point3D();                  << 
1944   for(int i=1;i<=nvert;i++) {                 << 
1945     centre += pV[i];                          << 
1946   }                                           << 
1947   centre = centre/nvert;                      << 
1948   return centre;                              << 
1949 }                                             << 
1950                                               << 
1951 HepPolyhedronTrd2::HepPolyhedronTrd2(G4double << 
1952                                      G4double << 
1953                                      G4double << 
1954 /********************************************    1180 /***********************************************************************
1955  *                                               1181  *                                                                     *
1956  * Name: HepPolyhedronTrd2                       1182  * Name: HepPolyhedronTrd2                           Date:    22.07.96 *
1957  * Author: E.Chernyaev (IHEP/Protvino)           1183  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
1958  *                                               1184  *                                                                     *
1959  * Function: Create GEANT4 TRD2-trapezoid        1185  * Function: Create GEANT4 TRD2-trapezoid                              *
1960  *                                               1186  *                                                                     *
1961  * Input: Dx1 - half-length along X at -Dz       1187  * Input: Dx1 - half-length along X at -Dz           8----7            *
1962  *        Dx2 - half-length along X ay +Dz       1188  *        Dx2 - half-length along X ay +Dz        5----6  !            *
1963  *        Dy1 - half-length along Y ay -Dz       1189  *        Dy1 - half-length along Y ay -Dz        !  4-!--3            *
1964  *        Dy2 - half-length along Y ay +Dz       1190  *        Dy2 - half-length along Y ay +Dz        1----2               *
1965  *        Dz  - half-length along Z              1191  *        Dz  - half-length along Z                                    *
1966  *                                               1192  *                                                                     *
1967  ********************************************    1193  ***********************************************************************/
1968 {                                                1194 {
1969   AllocateMemory(8,6);                           1195   AllocateMemory(8,6);
1970                                                  1196 
1971   pV[1] = G4Point3D(-Dx1,-Dy1,-Dz);           << 1197   pV[1] = HepPoint3D(-Dx1,-Dy1,-Dz);
1972   pV[2] = G4Point3D( Dx1,-Dy1,-Dz);           << 1198   pV[2] = HepPoint3D( Dx1,-Dy1,-Dz);
1973   pV[3] = G4Point3D( Dx1, Dy1,-Dz);           << 1199   pV[3] = HepPoint3D( Dx1, Dy1,-Dz);
1974   pV[4] = G4Point3D(-Dx1, Dy1,-Dz);           << 1200   pV[4] = HepPoint3D(-Dx1, Dy1,-Dz);
1975   pV[5] = G4Point3D(-Dx2,-Dy2, Dz);           << 1201   pV[5] = HepPoint3D(-Dx2,-Dy2, Dz);
1976   pV[6] = G4Point3D( Dx2,-Dy2, Dz);           << 1202   pV[6] = HepPoint3D( Dx2,-Dy2, Dz);
1977   pV[7] = G4Point3D( Dx2, Dy2, Dz);           << 1203   pV[7] = HepPoint3D( Dx2, Dy2, Dz);
1978   pV[8] = G4Point3D(-Dx2, Dy2, Dz);           << 1204   pV[8] = HepPoint3D(-Dx2, Dy2, Dz);
1979                                                  1205 
1980   CreatePrism();                                 1206   CreatePrism();
1981 }                                                1207 }
1982                                                  1208 
1983 HepPolyhedronTrd2::~HepPolyhedronTrd2() = def << 1209 HepPolyhedronTrap::HepPolyhedronTrap(HepDouble Dz,
1984                                               << 1210              HepDouble Theta,
1985 HepPolyhedronTrd1::HepPolyhedronTrd1(G4double << 1211              HepDouble Phi,
1986                                      G4double << 1212              HepDouble Dy1,
1987   : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) { << 1213              HepDouble Dx1,
1988                                               << 1214              HepDouble Dx2,
1989 HepPolyhedronTrd1::~HepPolyhedronTrd1() = def << 1215              HepDouble Alp1,
1990                                               << 1216              HepDouble Dy2,
1991 HepPolyhedronBox::HepPolyhedronBox(G4double D << 1217              HepDouble Dx3,
1992   : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}  << 1218              HepDouble Dx4,
1993                                               << 1219              HepDouble Alp2)
1994 HepPolyhedronBox::~HepPolyhedronBox() = defau << 
1995                                               << 
1996 HepPolyhedronTrap::HepPolyhedronTrap(G4double << 
1997                                      G4double << 
1998                                      G4double << 
1999                                      G4double << 
2000                                      G4double << 
2001                                      G4double << 
2002                                      G4double << 
2003                                      G4double << 
2004                                      G4double << 
2005                                      G4double << 
2006                                      G4double << 
2007 /********************************************    1220 /***********************************************************************
2008  *                                               1221  *                                                                     *
2009  * Name: HepPolyhedronTrap                       1222  * Name: HepPolyhedronTrap                           Date:    20.11.96 *
2010  * Author: E.Chernyaev                           1223  * Author: E.Chernyaev                               Revised:          *
2011  *                                               1224  *                                                                     *
2012  * Function: Create GEANT4 TRAP-trapezoid        1225  * Function: Create GEANT4 TRAP-trapezoid                              *
2013  *                                               1226  *                                                                     *
2014  * Input: DZ   - half-length in Z                1227  * Input: DZ   - half-length in Z                                      *
2015  *        Theta,Phi - polar angles of the lin    1228  *        Theta,Phi - polar angles of the line joining centres of the  *
2016  *                    faces at Z=-Dz and Z=+D    1229  *                    faces at Z=-Dz and Z=+Dz                         *
2017  *        Dy1  - half-length in Y of the face    1230  *        Dy1  - half-length in Y of the face at Z=-Dz                 *
2018  *        Dx1  - half-length in X of low edge    1231  *        Dx1  - half-length in X of low edge of the face at Z=-Dz     *
2019  *        Dx2  - half-length in X of top edge    1232  *        Dx2  - half-length in X of top edge of the face at Z=-Dz     *
2020  *        Alp1 - angle between Y-axis and the    1233  *        Alp1 - angle between Y-axis and the median joining top and   *
2021  *               low edges of the face at Z=-    1234  *               low edges of the face at Z=-Dz                        *
2022  *        Dy2  - half-length in Y of the face    1235  *        Dy2  - half-length in Y of the face at Z=+Dz                 *
2023  *        Dx3  - half-length in X of low edge    1236  *        Dx3  - half-length in X of low edge of the face at Z=+Dz     *
2024  *        Dx4  - half-length in X of top edge    1237  *        Dx4  - half-length in X of top edge of the face at Z=+Dz     *
2025  *        Alp2 - angle between Y-axis and the    1238  *        Alp2 - angle between Y-axis and the median joining top and   *
2026  *               low edges of the face at Z=+    1239  *               low edges of the face at Z=+Dz                        *
2027  *                                               1240  *                                                                     *
2028  ********************************************    1241  ***********************************************************************/
2029 {                                                1242 {
2030   G4double DzTthetaCphi = Dz*std::tan(Theta)* << 1243   HepDouble DzTthetaCphi = Dz*tan(Theta)*cos(Phi);
2031   G4double DzTthetaSphi = Dz*std::tan(Theta)* << 1244   HepDouble DzTthetaSphi = Dz*tan(Theta)*sin(Phi);
2032   G4double Dy1Talp1 = Dy1*std::tan(Alp1);     << 1245   HepDouble Dy1Talp1 = Dy1*tan(Alp1);
2033   G4double Dy2Talp2 = Dy2*std::tan(Alp2);     << 1246   HepDouble Dy2Talp2 = Dy2*tan(Alp2);
2034                                               << 1247   
2035   AllocateMemory(8,6);                           1248   AllocateMemory(8,6);
2036                                                  1249 
2037   pV[1] = G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx << 1250   pV[1] = HepPoint3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
2038   pV[2] = G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx << 1251   pV[2] = HepPoint3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
2039   pV[3] = G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx << 1252   pV[3] = HepPoint3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
2040   pV[4] = G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx << 1253   pV[4] = HepPoint3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
2041   pV[5] = G4Point3D( DzTthetaCphi-Dy2Talp2-Dx << 1254   pV[5] = HepPoint3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
2042   pV[6] = G4Point3D( DzTthetaCphi-Dy2Talp2+Dx << 1255   pV[6] = HepPoint3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
2043   pV[7] = G4Point3D( DzTthetaCphi+Dy2Talp2+Dx << 1256   pV[7] = HepPoint3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
2044   pV[8] = G4Point3D( DzTthetaCphi+Dy2Talp2-Dx << 1257   pV[8] = HepPoint3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
2045                                                  1258 
2046   CreatePrism();                                 1259   CreatePrism();
2047 }                                                1260 }
2048                                                  1261 
2049 HepPolyhedronTrap::~HepPolyhedronTrap() = def << 1262 HepPolyhedronCons::HepPolyhedronCons(HepDouble Rmn1,
2050                                               << 1263              HepDouble Rmx1,
2051 HepPolyhedronPara::HepPolyhedronPara(G4double << 1264              HepDouble Rmn2,
2052                                      G4double << 1265              HepDouble Rmx2, 
2053                                      G4double << 1266              HepDouble Dz,
2054   : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, << 1267              HepDouble Phi1,
2055                                               << 1268              HepDouble Dphi) 
2056 HepPolyhedronPara::~HepPolyhedronPara() = def << 
2057                                               << 
2058 HepPolyhedronParaboloid::HepPolyhedronParabol << 
2059                                               << 
2060                                               << 
2061                                               << 
2062                                               << 
2063 /******************************************** << 
2064  *                                            << 
2065  * Name: HepPolyhedronParaboloid              << 
2066  * Author: L.Lindroos, T.Nikitina (CERN), Jul << 
2067  *                                            << 
2068  * Function: Constructor for paraboloid       << 
2069  *                                            << 
2070  * Input: r1    - inside and outside radiuses << 
2071  *        r2    - inside and outside radiuses << 
2072  *        dz    - half length in Z            << 
2073  *        sPhi  - starting angle of the segme << 
2074  *        dPhi  - segment range               << 
2075  *                                            << 
2076  ******************************************** << 
2077 {                                             << 
2078   static const G4double wholeCircle=twopi;    << 
2079                                               << 
2080   //   C H E C K   I N P U T   P A R A M E T  << 
2081                                               << 
2082   G4int k = 0;                                << 
2083   if (r1 < 0. || r2 <= 0.)        k = 1;      << 
2084                                               << 
2085   if (dz <= 0.) k += 2;                       << 
2086                                               << 
2087   G4double phi1, phi2, dphi;                  << 
2088                                               << 
2089   if(dPhi < 0.)                               << 
2090   {                                           << 
2091     phi2 = sPhi; phi1 = phi2 + dPhi;          << 
2092   }                                           << 
2093   else if(dPhi == 0.)                         << 
2094   {                                           << 
2095     phi1 = sPhi; phi2 = phi1 + wholeCircle;   << 
2096   }                                           << 
2097   else                                        << 
2098   {                                           << 
2099     phi1 = sPhi; phi2 = phi1 + dPhi;          << 
2100   }                                           << 
2101   dphi  = phi2 - phi1;                        << 
2102                                               << 
2103   if (std::abs(dphi-wholeCircle) < perMillion << 
2104   if (dphi > wholeCircle) k += 4;             << 
2105                                               << 
2106   if (k != 0) {                               << 
2107     std::cerr << "HepPolyhedronParaboloid: er << 
2108     if ((k & 1) != 0) std::cerr << " (radiuse << 
2109     if ((k & 2) != 0) std::cerr << " (half-le << 
2110     if ((k & 4) != 0) std::cerr << " (angles) << 
2111     std::cerr << std::endl;                   << 
2112     std::cerr << " r1=" << r1;                << 
2113     std::cerr << " r2=" << r2;                << 
2114     std::cerr << " dz=" << dz << " sPhi=" <<  << 
2115               << std::endl;                   << 
2116     return;                                   << 
2117   }                                           << 
2118                                               << 
2119   //   P R E P A R E   T W O   P O L Y L I N  << 
2120                                               << 
2121   G4int n = GetNumberOfRotationSteps();       << 
2122   G4double dl = (r2 - r1) / n;                << 
2123   G4double k1 = (r2*r2 - r1*r1) / 2 / dz;     << 
2124   G4double k2 = (r2*r2 + r1*r1) / 2;          << 
2125                                               << 
2126   auto zz = new G4double[n + 2], rr = new G4d << 
2127                                               << 
2128   zz[0] = dz;                                 << 
2129   rr[0] = r2;                                 << 
2130                                               << 
2131   for(G4int i = 1; i < n - 1; i++)            << 
2132   {                                           << 
2133     rr[i] = rr[i-1] - dl;                     << 
2134     zz[i] = (rr[i]*rr[i] - k2) / k1;          << 
2135     if(rr[i] < 0)                             << 
2136     {                                         << 
2137       rr[i] = 0;                              << 
2138       zz[i] = 0;                              << 
2139     }                                         << 
2140   }                                           << 
2141                                               << 
2142   zz[n-1] = -dz;                              << 
2143   rr[n-1] = r1;                               << 
2144                                               << 
2145   zz[n] = dz;                                 << 
2146   rr[n] = 0;                                  << 
2147                                               << 
2148   zz[n+1] = -dz;                              << 
2149   rr[n+1] = 0;                                << 
2150                                               << 
2151   //   R O T A T E    P O L Y L I N E S       << 
2152                                               << 
2153   RotateAroundZ(0, phi1, dphi, n, 2, zz, rr,  << 
2154   SetReferences();                            << 
2155                                               << 
2156   delete [] zz;                               << 
2157   delete [] rr;                               << 
2158 }                                             << 
2159                                               << 
2160 HepPolyhedronParaboloid::~HepPolyhedronParabo << 
2161                                               << 
2162 HepPolyhedronHype::HepPolyhedronHype(G4double << 
2163                                      G4double << 
2164                                      G4double << 
2165                                      G4double << 
2166                                      G4double << 
2167 /******************************************** << 
2168  *                                            << 
2169  * Name: HepPolyhedronHype                    << 
2170  * Author: Tatiana Nikitina (CERN)            << 
2171  *         Evgueni Tcherniaev                 << 
2172  *                                            << 
2173  * Function: Constructor for Hype             << 
2174  *                                            << 
2175  * Input: r1       - inside radius at z=0     << 
2176  *        r2       - outside radiuses at z=0  << 
2177  *        sqrtan1  - sqr of tan of Inner Ster << 
2178  *        sqrtan2  - sqr of tan of Outer Ster << 
2179  *        halfZ    - half length in Z         << 
2180  *                                            << 
2181  ******************************************** << 
2182 {                                             << 
2183   static const G4double wholeCircle = twopi;  << 
2184                                               << 
2185   //   C H E C K   I N P U T   P A R A M E T  << 
2186                                               << 
2187   G4int k = 0;                                << 
2188   if (r1 < 0. || r2 < 0. || r1 >= r2) k = 1;  << 
2189   if (halfZ <= 0.) k += 2;                    << 
2190   if (sqrtan1 < 0.|| sqrtan2 < 0.) k += 4;    << 
2191                                               << 
2192   if (k != 0)                                 << 
2193   {                                           << 
2194     std::cerr << "HepPolyhedronHype: error in << 
2195     if ((k & 1) != 0) std::cerr << " (radiuse << 
2196     if ((k & 2) != 0) std::cerr << " (half-le << 
2197     if ((k & 4) != 0) std::cerr << " (angles) << 
2198     std::cerr << std::endl;                   << 
2199     std::cerr << " r1=" << r1 << " r2=" << r2 << 
2200     std::cerr << " halfZ=" << halfZ << " sqrT << 
2201               << " sqrTan2=" << sqrtan2       << 
2202               << std::endl;                   << 
2203     return;                                   << 
2204   }                                           << 
2205                                               << 
2206   //   P R E P A R E   T W O   P O L Y L I N  << 
2207                                               << 
2208   G4int ns = std::max(3, GetNumberOfRotationS << 
2209   G4int nz1 = (sqrtan1 == 0.) ? 2 : ns + 1;   << 
2210   G4int nz2 = (sqrtan2 == 0.) ? 2 : ns + 1;   << 
2211   auto  zz = new G4double[nz1 + nz2];         << 
2212   auto  rr = new G4double[nz1 + nz2];         << 
2213                                               << 
2214   // external polyline                        << 
2215   G4double dz2 = 2.*halfZ/(nz2 - 1);          << 
2216   for(G4int i = 0; i < nz2; ++i)              << 
2217   {                                           << 
2218     zz[i] = halfZ - dz2*i;                    << 
2219     rr[i] = std::sqrt(sqrtan2*zz[i]*zz[i] + r << 
2220   }                                           << 
2221                                               << 
2222   // internal polyline                        << 
2223   G4double dz1 = 2.*halfZ/(nz1 - 1);          << 
2224   for(G4int i = 0; i < nz1; ++i)              << 
2225   {                                           << 
2226     G4int j = nz2 + i;                        << 
2227     zz[j] = halfZ - dz1*i;                    << 
2228     rr[j] = std::sqrt(sqrtan1*zz[j]*zz[j] + r << 
2229   }                                           << 
2230                                               << 
2231   //   R O T A T E    P O L Y L I N E S       << 
2232                                               << 
2233   RotateAroundZ(0, 0., wholeCircle, nz2, nz1, << 
2234   SetReferences();                            << 
2235                                               << 
2236   delete [] zz;                               << 
2237   delete [] rr;                               << 
2238 }                                             << 
2239                                               << 
2240 HepPolyhedronHype::~HepPolyhedronHype() = def << 
2241                                               << 
2242 HepPolyhedronCons::HepPolyhedronCons(G4double << 
2243                                      G4double << 
2244                                      G4double << 
2245                                      G4double << 
2246                                      G4double << 
2247                                      G4double << 
2248                                      G4double << 
2249 /********************************************    1269 /***********************************************************************
2250  *                                               1270  *                                                                     *
2251  * Name: HepPolyhedronCons::HepPolyhedronCons    1271  * Name: HepPolyhedronCons::HepPolyhedronCons        Date:    15.12.96 *
2252  * Author: E.Chernyaev (IHEP/Protvino)           1272  * Author: E.Chernyaev (IHEP/Protvino)               Revised: 15.12.96 *
2253  *                                               1273  *                                                                     *
2254  * Function: Constructor for CONS, TUBS, CONE    1274  * Function: Constructor for CONS, TUBS, CONE, TUBE                    *
2255  *                                               1275  *                                                                     *
2256  * Input: Rmn1, Rmx1 - inside and outside rad    1276  * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz              *
2257  *        Rmn2, Rmx2 - inside and outside rad    1277  *        Rmn2, Rmx2 - inside and outside radiuses at +Dz              *
2258  *        Dz         - half length in Z          1278  *        Dz         - half length in Z                                *
2259  *        Phi1       - starting angle of the     1279  *        Phi1       - starting angle of the segment                   *
2260  *        Dphi       - segment range             1280  *        Dphi       - segment range                                   *
2261  *                                               1281  *                                                                     *
2262  ********************************************    1282  ***********************************************************************/
2263 {                                                1283 {
2264   static const G4double wholeCircle=twopi;    << 1284   static HepDouble wholeCircle=2*M_PI;
2265                                                  1285 
2266   //   C H E C K   I N P U T   P A R A M E T     1286   //   C H E C K   I N P U T   P A R A M E T E R S
2267                                                  1287 
2268   G4int k = 0;                                << 1288   int k = 0;
2269   if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. ||     1289   if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.)        k = 1;
2270   if (Rmn1 > Rmx1 || Rmn2 > Rmx2)                1290   if (Rmn1 > Rmx1 || Rmn2 > Rmx2)                              k = 1;
2271   if (Rmn1 == Rmx1 && Rmn2 == Rmx2)              1291   if (Rmn1 == Rmx1 && Rmn2 == Rmx2)                            k = 1;
2272                                                  1292 
2273   if (Dz <= 0.) k += 2;                          1293   if (Dz <= 0.) k += 2;
2274                                               << 1294  
2275   G4double phi1, phi2, dphi;                  << 1295   HepDouble phi1, phi2, dphi;
2276   if (Dphi < 0.) {                               1296   if (Dphi < 0.) {
2277     phi2 = Phi1; phi1 = phi2 - Dphi;             1297     phi2 = Phi1; phi1 = phi2 - Dphi;
2278   }else if (Dphi == 0.) {                        1298   }else if (Dphi == 0.) {
2279     phi1 = Phi1; phi2 = phi1 + wholeCircle;      1299     phi1 = Phi1; phi2 = phi1 + wholeCircle;
2280   }else{                                         1300   }else{
2281     phi1 = Phi1; phi2 = phi1 + Dphi;             1301     phi1 = Phi1; phi2 = phi1 + Dphi;
2282   }                                              1302   }
2283   dphi  = phi2 - phi1;                           1303   dphi  = phi2 - phi1;
2284   if (std::abs(dphi-wholeCircle) < perMillion << 1304   if (abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
2285   if (dphi > wholeCircle) k += 4;             << 1305   if (dphi > wholeCircle) k += 4; 
2286                                                  1306 
2287   if (k != 0) {                                  1307   if (k != 0) {
2288     std::cerr << "HepPolyhedronCone(s)/Tube(s << 1308     HepStd::cerr << "HepPolyhedronCone(s)/Tube(s): error in input parameters";
2289     if ((k & 1) != 0) std::cerr << " (radiuse << 1309     if ((k & 1) != 0) HepStd::cerr << " (radiuses)";
2290     if ((k & 2) != 0) std::cerr << " (half-le << 1310     if ((k & 2) != 0) HepStd::cerr << " (half-length)";
2291     if ((k & 4) != 0) std::cerr << " (angles) << 1311     if ((k & 4) != 0) HepStd::cerr << " (angles)";
2292     std::cerr << std::endl;                   << 1312     HepStd::cerr << HepStd::endl;
2293     std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << 1313     HepStd::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
2294     std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << 1314     HepStd::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
2295     std::cerr << " Dz=" << Dz << " Phi1=" <<  << 1315     HepStd::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
2296               << std::endl;                   << 1316         << HepStd::endl;
2297     return;                                      1317     return;
2298   }                                              1318   }
2299                                               << 1319   
2300   //   P R E P A R E   T W O   P O L Y L I N     1320   //   P R E P A R E   T W O   P O L Y L I N E S
2301                                                  1321 
2302   G4double zz[4], rr[4];                      << 1322   HepDouble zz[4], rr[4];
2303   zz[0] =  Dz;                                << 1323   zz[0] =  Dz; 
2304   zz[1] = -Dz;                                << 1324   zz[1] = -Dz; 
2305   zz[2] =  Dz;                                << 1325   zz[2] =  Dz; 
2306   zz[3] = -Dz;                                << 1326   zz[3] = -Dz; 
2307   rr[0] =  Rmx2;                                 1327   rr[0] =  Rmx2;
2308   rr[1] =  Rmx1;                                 1328   rr[1] =  Rmx1;
2309   rr[2] =  Rmn2;                                 1329   rr[2] =  Rmn2;
2310   rr[3] =  Rmn1;                                 1330   rr[3] =  Rmn1;
2311                                                  1331 
2312   //   R O T A T E    P O L Y L I N E S          1332   //   R O T A T E    P O L Y L I N E S
2313                                                  1333 
2314   RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr,  << 1334   RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1); 
2315   SetReferences();                               1335   SetReferences();
2316 }                                                1336 }
2317                                                  1337 
2318 HepPolyhedronCons::~HepPolyhedronCons() = def << 1338 HepPolyhedronPgon::HepPolyhedronPgon(HepDouble phi,
2319                                               << 1339              HepDouble dphi,
2320 HepPolyhedronCone::HepPolyhedronCone(G4double << 1340              int    npdv,
2321                                      G4double << 1341              int    nz,
2322                                      G4double << 1342              const HepDouble *z,
2323   HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, D << 1343              const HepDouble *rmin,
2324                                               << 1344              const HepDouble *rmax)
2325 HepPolyhedronCone::~HepPolyhedronCone() = def << 
2326                                               << 
2327 HepPolyhedronTubs::HepPolyhedronTubs(G4double << 
2328                                      G4double << 
2329                                      G4double << 
2330   :   HepPolyhedronCons(Rmin, Rmax, Rmin, Rma << 
2331                                               << 
2332 HepPolyhedronTubs::~HepPolyhedronTubs() = def << 
2333                                               << 
2334 HepPolyhedronTube::HepPolyhedronTube (G4doubl << 
2335                                       G4doubl << 
2336   : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, << 
2337                                               << 
2338 HepPolyhedronTube::~HepPolyhedronTube () = de << 
2339                                               << 
2340 HepPolyhedronPgon::HepPolyhedronPgon(G4double << 
2341                                      G4double << 
2342                                      G4int np << 
2343                                      G4int nz << 
2344                                      const G4 << 
2345                                      const G4 << 
2346                                      const G4 << 
2347 /********************************************    1345 /***********************************************************************
2348  *                                               1346  *                                                                     *
2349  * Name: HepPolyhedronPgon                       1347  * Name: HepPolyhedronPgon                           Date:    09.12.96 *
2350  * Author: E.Chernyaev                           1348  * Author: E.Chernyaev                               Revised:          *
2351  *                                               1349  *                                                                     *
2352  * Function: Constructor of polyhedron for PG    1350  * Function: Constructor of polyhedron for PGON, PCON                  *
2353  *                                               1351  *                                                                     *
2354  * Input: phi  - initial phi                     1352  * Input: phi  - initial phi                                           *
2355  *        dphi - delta phi                       1353  *        dphi - delta phi                                             *
2356  *        npdv - number of steps along phi       1354  *        npdv - number of steps along phi                             *
2357  *        nz   - number of z-planes (at least    1355  *        nz   - number of z-planes (at least two)                     *
2358  *        z[]  - z coordinates of the slices     1356  *        z[]  - z coordinates of the slices                           *
2359  *        rmin[] - smaller r at the slices       1357  *        rmin[] - smaller r at the slices                             *
2360  *        rmax[] - bigger  r at the slices       1358  *        rmax[] - bigger  r at the slices                             *
2361  *                                               1359  *                                                                     *
2362  ********************************************    1360  ***********************************************************************/
2363 {                                                1361 {
2364   //   C H E C K   I N P U T   P A R A M E T     1362   //   C H E C K   I N P U T   P A R A M E T E R S
2365                                                  1363 
2366   if (dphi <= 0. || dphi > twopi) {           << 1364   if (dphi <= 0. || dphi > 2*M_PI) {
2367     std::cerr                                 << 1365     HepStd::cerr
2368       << "HepPolyhedronPgon/Pcon: wrong delta    1366       << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
2369       << std::endl;                           << 1367       << HepStd::endl;
2370     return;                                      1368     return;
2371   }                                           << 1369   }    
2372                                               << 1370     
2373   if (nz < 2) {                                  1371   if (nz < 2) {
2374     std::cerr                                 << 1372     HepStd::cerr
2375       << "HepPolyhedronPgon/Pcon: number of z    1373       << "HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
2376       << std::endl;                           << 1374       << HepStd::endl;
2377     return;                                      1375     return;
2378   }                                              1376   }
2379                                                  1377 
2380   if (npdv < 0) {                                1378   if (npdv < 0) {
2381     std::cerr                                 << 1379     HepStd::cerr
2382       << "HepPolyhedronPgon/Pcon: error in nu    1380       << "HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
2383       << std::endl;                           << 1381       << HepStd::endl;
2384     return;                                      1382     return;
2385   }                                              1383   }
2386                                                  1384 
2387   G4int i;                                    << 1385   int i;
2388   for (i=0; i<nz; i++) {                         1386   for (i=0; i<nz; i++) {
2389     if (rmin[i] < 0. || rmax[i] < 0. || rmin[    1387     if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
2390       std::cerr                               << 1388       HepStd::cerr
2391         << "HepPolyhedronPgon: error in radiu << 1389   << "HepPolyhedronPgon: error in radiuses rmin[" << i << "]="
2392         << rmin[i] << " rmax[" << i << "]=" < << 1390   << rmin[i] << " rmax[" << i << "]=" << rmax[i]
2393         << std::endl;                         << 1391   << HepStd::endl;
2394       return;                                    1392       return;
2395     }                                            1393     }
2396   }                                              1394   }
2397                                                  1395 
2398   //   P R E P A R E   T W O   P O L Y L I N     1396   //   P R E P A R E   T W O   P O L Y L I N E S
2399                                                  1397 
2400   G4double *zz, *rr;                          << 1398   HepDouble *zz, *rr;
2401   zz = new G4double[2*nz];                    << 1399   zz = new HepDouble[2*nz];
2402   rr = new G4double[2*nz];                    << 1400   rr = new HepDouble[2*nz];
2403                                                  1401 
2404   if (z[0] > z[nz-1]) {                          1402   if (z[0] > z[nz-1]) {
2405     for (i=0; i<nz; i++) {                       1403     for (i=0; i<nz; i++) {
2406       zz[i]    = z[i];                           1404       zz[i]    = z[i];
2407       rr[i]    = rmax[i];                        1405       rr[i]    = rmax[i];
2408       zz[i+nz] = z[i];                           1406       zz[i+nz] = z[i];
2409       rr[i+nz] = rmin[i];                        1407       rr[i+nz] = rmin[i];
2410     }                                            1408     }
2411   }else{                                         1409   }else{
2412     for (i=0; i<nz; i++) {                       1410     for (i=0; i<nz; i++) {
2413       zz[i]    = z[nz-i-1];                      1411       zz[i]    = z[nz-i-1];
2414       rr[i]    = rmax[nz-i-1];                   1412       rr[i]    = rmax[nz-i-1];
2415       zz[i+nz] = z[nz-i-1];                      1413       zz[i+nz] = z[nz-i-1];
2416       rr[i+nz] = rmin[nz-i-1];                   1414       rr[i+nz] = rmin[nz-i-1];
2417     }                                            1415     }
2418   }                                              1416   }
2419                                                  1417 
2420   //   R O T A T E    P O L Y L I N E S          1418   //   R O T A T E    P O L Y L I N E S
2421                                                  1419 
2422   G4int nodeVis = 1;                          << 1420   RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1); 
2423   G4int edgeVis = (npdv == 0) ? -1 : 1;       << 
2424   RotateAroundZ(npdv, phi, dphi, nz, nz, zz,  << 
2425   SetReferences();                               1421   SetReferences();
2426                                               << 1422   
2427   delete [] zz;                                  1423   delete [] zz;
2428   delete [] rr;                                  1424   delete [] rr;
2429 }                                                1425 }
2430                                                  1426 
2431 HepPolyhedronPgon::HepPolyhedronPgon(G4double << 1427 HepPolyhedronSphere::HepPolyhedronSphere(HepDouble rmin, HepDouble rmax,
2432                                      G4double << 1428                HepDouble phi, HepDouble dphi,
2433                                      G4int np << 1429                HepDouble the, HepDouble dthe)
2434                                      const st << 
2435 /******************************************** << 
2436  *                                            << 
2437  * Name: HepPolyhedronPgon                    << 
2438  * Author: E.Tcherniaev (E.Chernyaev)         << 
2439  *                                            << 
2440  * Function: Constructor of polyhedron for PG << 
2441  *                                            << 
2442  * Input: phi  - initial phi                  << 
2443  *        dphi - delta phi                    << 
2444  *        npdv - number of steps along phi    << 
2445  *        rz   - rz-contour                   << 
2446  *                                            << 
2447  ******************************************** << 
2448 {                                             << 
2449   //   C H E C K   I N P U T   P A R A M E T  << 
2450                                               << 
2451   if (dphi <= 0. || dphi > twopi) {           << 
2452     std::cerr                                 << 
2453       << "HepPolyhedronPgon/Pcon: wrong delta << 
2454       << std::endl;                           << 
2455     return;                                   << 
2456   }                                           << 
2457                                               << 
2458   if (npdv < 0) {                             << 
2459     std::cerr                                 << 
2460       << "HepPolyhedronPgon/Pcon: error in nu << 
2461       << std::endl;                           << 
2462     return;                                   << 
2463   }                                           << 
2464                                               << 
2465   G4int nrz = (G4int)rz.size();               << 
2466   if (nrz < 3) {                              << 
2467     std::cerr                                 << 
2468       << "HepPolyhedronPgon/Pcon: invalid num << 
2469       << std::endl;                           << 
2470     return;                                   << 
2471   }                                           << 
2472                                               << 
2473   //   R O T A T E    P O L Y L I N E         << 
2474                                               << 
2475   G4int nodeVis = 1;                          << 
2476   G4int edgeVis = (npdv == 0) ? -1 : 1;       << 
2477   RotateContourAroundZ(npdv, phi, dphi, rz, n << 
2478   SetReferences();                            << 
2479 }                                             << 
2480                                               << 
2481 HepPolyhedronPgon::~HepPolyhedronPgon() = def << 
2482                                               << 
2483 HepPolyhedronPcon::HepPolyhedronPcon(G4double << 
2484                                      const G4 << 
2485                                      const G4 << 
2486                                      const G4 << 
2487   : HepPolyhedronPgon(phi, dphi, 0, nz, z, rm << 
2488                                               << 
2489 HepPolyhedronPcon::HepPolyhedronPcon(G4double << 
2490                                      const st << 
2491   : HepPolyhedronPgon(phi, dphi, 0, rz) {}    << 
2492                                               << 
2493 HepPolyhedronPcon::~HepPolyhedronPcon() = def << 
2494                                               << 
2495 HepPolyhedronSphere::HepPolyhedronSphere(G4do << 
2496                                          G4do << 
2497                                          G4do << 
2498 /********************************************    1430 /***********************************************************************
2499  *                                               1431  *                                                                     *
2500  * Name: HepPolyhedronSphere                     1432  * Name: HepPolyhedronSphere                         Date:    11.12.96 *
2501  * Author: E.Chernyaev (IHEP/Protvino)           1433  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
2502  *                                               1434  *                                                                     *
2503  * Function: Constructor of polyhedron for SP    1435  * Function: Constructor of polyhedron for SPHERE                      *
2504  *                                               1436  *                                                                     *
2505  * Input: rmin - internal radius                 1437  * Input: rmin - internal radius                                       *
2506  *        rmax - external radius                 1438  *        rmax - external radius                                       *
2507  *        phi  - initial phi                     1439  *        phi  - initial phi                                           *
2508  *        dphi - delta phi                       1440  *        dphi - delta phi                                             *
2509  *        the  - initial theta                   1441  *        the  - initial theta                                         *
2510  *        dthe - delta theta                     1442  *        dthe - delta theta                                           *
2511  *                                               1443  *                                                                     *
2512  ********************************************    1444  ***********************************************************************/
2513 {                                                1445 {
2514   //   C H E C K   I N P U T   P A R A M E T     1446   //   C H E C K   I N P U T   P A R A M E T E R S
2515                                                  1447 
2516   if (dphi <= 0. || dphi > twopi) {           << 1448   if (dphi <= 0. || dphi > 2*M_PI) {
2517     std::cerr                                 << 1449     HepStd::cerr
2518       << "HepPolyhedronSphere: wrong delta ph    1450       << "HepPolyhedronSphere: wrong delta phi = " << dphi
2519       << std::endl;                           << 1451       << HepStd::endl;
2520     return;                                      1452     return;
2521   }                                           << 1453   }    
2522                                                  1454 
2523   if (the < 0. || the > pi) {                 << 1455   if (the < 0. || the > M_PI) {
2524     std::cerr                                 << 1456     HepStd::cerr
2525       << "HepPolyhedronSphere: wrong theta =     1457       << "HepPolyhedronSphere: wrong theta = " << the
2526       << std::endl;                           << 1458       << HepStd::endl;
2527     return;                                      1459     return;
2528   }                                           << 1460   }    
2529                                               << 1461   
2530   if (dthe <= 0. || dthe > pi) {              << 1462   if (dthe <= 0. || dthe > M_PI) {
2531     std::cerr                                 << 1463     HepStd::cerr
2532       << "HepPolyhedronSphere: wrong delta th    1464       << "HepPolyhedronSphere: wrong delta theta = " << dthe
2533       << std::endl;                           << 1465       << HepStd::endl;
2534     return;                                      1466     return;
2535   }                                           << 1467   }    
2536                                                  1468 
2537   if (the+dthe > pi) {                        << 1469   if (the+dthe > M_PI) {
2538     std::cerr                                 << 1470     HepStd::cerr
2539       << "HepPolyhedronSphere: wrong theta +     1471       << "HepPolyhedronSphere: wrong theta + delta theta = "
2540       << the << " " << dthe                      1472       << the << " " << dthe
2541       << std::endl;                           << 1473       << HepStd::endl;
2542     return;                                      1474     return;
2543   }                                           << 1475   }    
2544                                               << 1476   
2545   if (rmin < 0. || rmin >= rmax) {               1477   if (rmin < 0. || rmin >= rmax) {
2546     std::cerr                                 << 1478     HepStd::cerr
2547       << "HepPolyhedronSphere: error in radiu    1479       << "HepPolyhedronSphere: error in radiuses"
2548       << " rmin=" << rmin << " rmax=" << rmax    1480       << " rmin=" << rmin << " rmax=" << rmax
2549       << std::endl;                           << 1481       << HepStd::endl;
2550     return;                                      1482     return;
2551   }                                              1483   }
2552                                                  1484 
2553   //   P R E P A R E   T W O   P O L Y L I N     1485   //   P R E P A R E   T W O   P O L Y L I N E S
2554                                                  1486 
2555   G4int nds = (GetNumberOfRotationSteps() + 1 << 1487   int ns = (GetNumberOfRotationSteps() + 1) / 2;
2556   G4int np1 = G4int(dthe*nds/pi+.5) + 1;      << 1488   int np1 = int(dthe*ns/M_PI+.5) + 1;
2557   if (np1 <= 1) np1 = 2;                         1489   if (np1 <= 1) np1 = 2;
2558   G4int np2 = rmin < spatialTolerance ? 1 : n << 1490   int np2 = rmin < perMillion ? 1 : np1;
2559                                                  1491 
2560   G4double *zz, *rr;                          << 1492   HepDouble *zz, *rr;
2561   zz = new G4double[np1+np2];                 << 1493   zz = new HepDouble[np1+np2];
2562   rr = new G4double[np1+np2];                 << 1494   rr = new HepDouble[np1+np2];
2563                                               << 1495 
2564   G4double a = dthe/(np1-1);                  << 1496   HepDouble a = dthe/(np1-1);
2565   G4double cosa, sina;                        << 1497   HepDouble cosa, sina;
2566   for (G4int i=0; i<np1; i++) {               << 1498   for (int i=0; i<np1; i++) {
2567     cosa  = std::cos(the+i*a);                << 1499     cosa  = cos(the+i*a);
2568     sina  = std::sin(the+i*a);                << 1500     sina  = sin(the+i*a);
2569     zz[i] = rmax*cosa;                           1501     zz[i] = rmax*cosa;
2570     rr[i] = rmax*sina;                           1502     rr[i] = rmax*sina;
2571     if (np2 > 1) {                               1503     if (np2 > 1) {
2572       zz[i+np1] = rmin*cosa;                     1504       zz[i+np1] = rmin*cosa;
2573       rr[i+np1] = rmin*sina;                     1505       rr[i+np1] = rmin*sina;
2574     }                                            1506     }
2575   }                                              1507   }
2576   if (np2 == 1) {                                1508   if (np2 == 1) {
2577     zz[np1] = 0.;                                1509     zz[np1] = 0.;
2578     rr[np1] = 0.;                                1510     rr[np1] = 0.;
2579   }                                              1511   }
2580                                                  1512 
2581   //   R O T A T E    P O L Y L I N E S          1513   //   R O T A T E    P O L Y L I N E S
2582                                                  1514 
2583   RotateAroundZ(0, phi, dphi, np1, np2, zz, r << 1515   RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1); 
2584   SetReferences();                               1516   SetReferences();
2585                                               << 1517   
2586   delete [] zz;                                  1518   delete [] zz;
2587   delete [] rr;                                  1519   delete [] rr;
2588 }                                                1520 }
2589                                                  1521 
2590 HepPolyhedronSphere::~HepPolyhedronSphere() = << 1522 HepPolyhedronTorus::HepPolyhedronTorus(HepDouble rmin,
2591                                               << 1523                HepDouble rmax,
2592 HepPolyhedronTorus::HepPolyhedronTorus(G4doub << 1524                HepDouble rtor,
2593                                        G4doub << 1525                HepDouble phi,
2594                                        G4doub << 1526                HepDouble dphi)
2595                                        G4doub << 
2596                                        G4doub << 
2597 /********************************************    1527 /***********************************************************************
2598  *                                               1528  *                                                                     *
2599  * Name: HepPolyhedronTorus                      1529  * Name: HepPolyhedronTorus                          Date:    11.12.96 *
2600  * Author: E.Chernyaev (IHEP/Protvino)           1530  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
2601  *                                               1531  *                                                                     *
2602  * Function: Constructor of polyhedron for TO    1532  * Function: Constructor of polyhedron for TORUS                       *
2603  *                                               1533  *                                                                     *
2604  * Input: rmin - internal radius                 1534  * Input: rmin - internal radius                                       *
2605  *        rmax - external radius                 1535  *        rmax - external radius                                       *
2606  *        rtor - radius of torus                 1536  *        rtor - radius of torus                                       *
2607  *        phi  - initial phi                     1537  *        phi  - initial phi                                           *
2608  *        dphi - delta phi                       1538  *        dphi - delta phi                                             *
2609  *                                               1539  *                                                                     *
2610  ********************************************    1540  ***********************************************************************/
2611 {                                                1541 {
2612   //   C H E C K   I N P U T   P A R A M E T     1542   //   C H E C K   I N P U T   P A R A M E T E R S
2613                                                  1543 
2614   if (dphi <= 0. || dphi > twopi) {           << 1544   if (dphi <= 0. || dphi > 2*M_PI) {
2615     std::cerr                                 << 1545     HepStd::cerr
2616       << "HepPolyhedronTorus: wrong delta phi    1546       << "HepPolyhedronTorus: wrong delta phi = " << dphi
2617       << std::endl;                           << 1547       << HepStd::endl;
2618     return;                                      1548     return;
2619   }                                              1549   }
2620                                                  1550 
2621   if (rmin < 0. || rmin >= rmax || rmax >= rt    1551   if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2622     std::cerr                                 << 1552     HepStd::cerr
2623       << "HepPolyhedronTorus: error in radius    1553       << "HepPolyhedronTorus: error in radiuses"
2624       << " rmin=" << rmin << " rmax=" << rmax    1554       << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
2625       << std::endl;                           << 1555       << HepStd::endl;
2626     return;                                      1556     return;
2627   }                                              1557   }
2628                                                  1558 
2629   //   P R E P A R E   T W O   P O L Y L I N     1559   //   P R E P A R E   T W O   P O L Y L I N E S
2630                                                  1560 
2631   G4int np1 = GetNumberOfRotationSteps();     << 1561   int np1 = GetNumberOfRotationSteps();
2632   G4int np2 = rmin < spatialTolerance ? 1 : n << 1562   int np2 = rmin < perMillion ? 1 : np1;
2633                                                  1563 
2634   G4double *zz, *rr;                          << 1564   HepDouble *zz, *rr;
2635   zz = new G4double[np1+np2];                 << 1565   zz = new HepDouble[np1+np2];
2636   rr = new G4double[np1+np2];                 << 1566   rr = new HepDouble[np1+np2];
2637                                               << 1567 
2638   G4double a = twopi/np1;                     << 1568   HepDouble a = 2*M_PI/np1;
2639   G4double cosa, sina;                        << 1569   HepDouble cosa, sina;
2640   for (G4int i=0; i<np1; i++) {               << 1570   for (int i=0; i<np1; i++) {
2641     cosa  = std::cos(i*a);                    << 1571     cosa  = cos(i*a);
2642     sina  = std::sin(i*a);                    << 1572     sina  = sin(i*a);
2643     zz[i] = rmax*cosa;                           1573     zz[i] = rmax*cosa;
2644     rr[i] = rtor+rmax*sina;                      1574     rr[i] = rtor+rmax*sina;
2645     if (np2 > 1) {                               1575     if (np2 > 1) {
2646       zz[i+np1] = rmin*cosa;                     1576       zz[i+np1] = rmin*cosa;
2647       rr[i+np1] = rtor+rmin*sina;                1577       rr[i+np1] = rtor+rmin*sina;
2648     }                                            1578     }
2649   }                                              1579   }
2650   if (np2 == 1) {                                1580   if (np2 == 1) {
2651     zz[np1] = 0.;                                1581     zz[np1] = 0.;
2652     rr[np1] = rtor;                              1582     rr[np1] = rtor;
2653     np2 = -1;                                    1583     np2 = -1;
2654   }                                              1584   }
2655                                                  1585 
2656   //   R O T A T E    P O L Y L I N E S          1586   //   R O T A T E    P O L Y L I N E S
2657                                                  1587 
2658   RotateAroundZ(0, phi, dphi, -np1, -np2, zz, << 1588   RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1); 
2659   SetReferences();                               1589   SetReferences();
2660                                               << 1590   
2661   delete [] zz;                               << 
2662   delete [] rr;                               << 
2663 }                                             << 
2664                                               << 
2665 HepPolyhedronTorus::~HepPolyhedronTorus() = d << 
2666                                               << 
2667 HepPolyhedronTet::HepPolyhedronTet(const G4do << 
2668                                    const G4do << 
2669                                    const G4do << 
2670                                    const G4do << 
2671 /******************************************** << 
2672  *                                            << 
2673  * Name: HepPolyhedronTet                     << 
2674  * Author: E.Tcherniaev (E.Chernyaev)         << 
2675  *                                            << 
2676  * Function: Constructor of polyhedron for TE << 
2677  *                                            << 
2678  * Input: p0,p1,p2,p3 - vertices              << 
2679  *                                            << 
2680  ******************************************** << 
2681 {                                             << 
2682   AllocateMemory(4,4);                        << 
2683                                               << 
2684   pV[1].set(p0[0], p0[1], p0[2]);             << 
2685   pV[2].set(p1[0], p1[1], p1[2]);             << 
2686   pV[3].set(p2[0], p2[1], p2[2]);             << 
2687   pV[4].set(p3[0], p3[1], p3[2]);             << 
2688                                               << 
2689   G4Vector3D v1(pV[2] - pV[1]);               << 
2690   G4Vector3D v2(pV[3] - pV[1]);               << 
2691   G4Vector3D v3(pV[4] - pV[1]);               << 
2692                                               << 
2693   if (v1.cross(v2).dot(v3) < 0.)              << 
2694   {                                           << 
2695     pV[3].set(p3[0], p3[1], p3[2]);           << 
2696     pV[4].set(p2[0], p2[1], p2[2]);           << 
2697   }                                           << 
2698                                               << 
2699   pF[1] = G4Facet(1,2,  3,4,  2,3);           << 
2700   pF[2] = G4Facet(1,3,  4,4,  3,1);           << 
2701   pF[3] = G4Facet(1,1,  2,4,  4,2);           << 
2702   pF[4] = G4Facet(2,1,  3,2,  4,3);           << 
2703 }                                             << 
2704                                               << 
2705 HepPolyhedronTet::~HepPolyhedronTet() = defau << 
2706                                               << 
2707 HepPolyhedronEllipsoid::HepPolyhedronEllipsoi << 
2708                                               << 
2709                                               << 
2710 /******************************************** << 
2711  *                                            << 
2712  * Name: HepPolyhedronEllipsoid               << 
2713  * Author: G.Guerrieri                        << 
2714  *         Evgueni Tcherniaev                 << 
2715  *                                            << 
2716  * Function: Constructor of polyhedron for EL << 
2717  *                                            << 
2718  * Input: ax - semiaxis x                     << 
2719  *        by - semiaxis y                     << 
2720  *        cz - semiaxis z                     << 
2721  *        zCut1 - lower cut plane level (soli << 
2722  *        zCut2 - upper cut plane level (soli << 
2723  *                                            << 
2724  ******************************************** << 
2725 {                                             << 
2726   //   C H E C K   I N P U T   P A R A M E T  << 
2727                                               << 
2728   if (zCut1 >= cz || zCut2 <= -cz || zCut1 >  << 
2729     std::cerr << "HepPolyhedronEllipsoid: wro << 
2730            << " zCut2 = " << zCut2            << 
2731            << " for given cz = " << cz << std << 
2732     return;                                   << 
2733   }                                           << 
2734   if (cz <= 0.0) {                            << 
2735     std::cerr << "HepPolyhedronEllipsoid: bad << 
2736       << std::endl;                           << 
2737     return;                                   << 
2738   }                                           << 
2739                                               << 
2740   //   P R E P A R E   T W O   P O L Y L I N  << 
2741   //   generate sphere of radius cz first, th << 
2742                                               << 
2743   G4double sthe = std::acos(zCut2/cz);        << 
2744   G4double dthe = std::acos(zCut1/cz) - sthe; << 
2745   G4int nds = (GetNumberOfRotationSteps() + 1 << 
2746   G4int np1 = G4int(dthe*nds/pi + 0.5) + 1;   << 
2747   if (np1 <= 1) np1 = 2;                      << 
2748   G4int np2 = 2;                              << 
2749                                               << 
2750   G4double *zz, *rr;                          << 
2751   zz = new G4double[np1 + np2];               << 
2752   rr = new G4double[np1 + np2];               << 
2753   if ((zz == nullptr) || (rr == nullptr))     << 
2754   {                                           << 
2755     G4Exception("HepPolyhedronEllipsoid::HepP << 
2756                 "greps1002", FatalException,  << 
2757   }                                           << 
2758                                               << 
2759   G4double a = dthe/(np1 - 1);                << 
2760   G4double cosa, sina;                        << 
2761   for (G4int i = 0; i < np1; ++i)             << 
2762   {                                           << 
2763     cosa  = std::cos(sthe + i*a);             << 
2764     sina  = std::sin(sthe + i*a);             << 
2765     zz[i] = cz*cosa;                          << 
2766     rr[i] = cz*sina;                          << 
2767   }                                           << 
2768   zz[np1 + 0] = zCut2;                        << 
2769   rr[np1 + 0] = 0.;                           << 
2770   zz[np1 + 1] = zCut1;                        << 
2771   rr[np1 + 1] = 0.;                           << 
2772                                               << 
2773   //   R O T A T E    P O L Y L I N E S       << 
2774                                               << 
2775   RotateAroundZ(0, 0., twopi, np1, np2, zz, r << 
2776   SetReferences();                            << 
2777                                               << 
2778   delete [] zz;                               << 
2779   delete [] rr;                               << 
2780                                               << 
2781   // rescale x and y vertex coordinates       << 
2782   G4double kx = ax/cz;                        << 
2783   G4double ky = by/cz;                        << 
2784   G4Point3D* p = pV;                          << 
2785   for (G4int i = 0; i < nvert; ++i, ++p)      << 
2786   {                                           << 
2787     p->setX(p->x()*kx);                       << 
2788     p->setY(p->y()*ky);                       << 
2789   }                                           << 
2790 }                                             << 
2791                                               << 
2792 HepPolyhedronEllipsoid::~HepPolyhedronEllipso << 
2793                                               << 
2794 HepPolyhedronEllipticalCone::HepPolyhedronEll << 
2795                                               << 
2796                                               << 
2797                                               << 
2798 /******************************************** << 
2799  *                                            << 
2800  * Name: HepPolyhedronEllipticalCone          << 
2801  * Author: D.Anninos                          << 
2802  *                                            << 
2803  * Function: Constructor for EllipticalCone   << 
2804  *                                            << 
2805  * Input: ax, ay     - X & Y semi axes at z = << 
2806  *        h          - height of full cone    << 
2807  *        zTopCut    - Top Cut in Z Axis      << 
2808  *                                            << 
2809  ******************************************** << 
2810 {                                             << 
2811   //   C H E C K   I N P U T   P A R A M E T  << 
2812                                               << 
2813   G4int k = 0;                                << 
2814   if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.)  << 
2815                                               << 
2816   if (k != 0) {                               << 
2817     std::cerr << "HepPolyhedronCone: error in << 
2818     std::cerr << std::endl;                   << 
2819     return;                                   << 
2820   }                                           << 
2821                                               << 
2822   //   P R E P A R E   T W O   P O L Y L I N  << 
2823                                               << 
2824   zTopCut = (h >= zTopCut ? zTopCut : h);     << 
2825                                               << 
2826   G4double *zz, *rr;                          << 
2827   zz = new G4double[4];                       << 
2828   rr = new G4double[4];                       << 
2829   zz[0] =   zTopCut;                          << 
2830   zz[1] =  -zTopCut;                          << 
2831   zz[2] =   zTopCut;                          << 
2832   zz[3] =  -zTopCut;                          << 
2833   rr[0] =  (h-zTopCut);                       << 
2834   rr[1] =  (h+zTopCut);                       << 
2835   rr[2] =  0.;                                << 
2836   rr[3] =  0.;                                << 
2837                                               << 
2838   //   R O T A T E    P O L Y L I N E S       << 
2839                                               << 
2840   RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, - << 
2841   SetReferences();                            << 
2842                                               << 
2843   delete [] zz;                               << 
2844   delete [] rr;                               << 
2845                                               << 
2846   // rescale x and y vertex coordinates       << 
2847  {                                            << 
2848    G4Point3D * p= pV;                         << 
2849    for (G4int i=0; i<nvert; i++, p++) {       << 
2850      p->setX( p->x() * ax );                  << 
2851      p->setY( p->y() * ay );                  << 
2852    }                                          << 
2853  }                                            << 
2854 }                                             << 
2855                                               << 
2856 HepPolyhedronEllipticalCone::~HepPolyhedronEl << 
2857                                               << 
2858 HepPolyhedronHyperbolicMirror::HepPolyhedronH << 
2859                                               << 
2860                                               << 
2861 /******************************************** << 
2862  *                                            << 
2863  * Name: HepPolyhedronHyperbolicMirror        << 
2864  * Author: E.Tcherniaev (E.Chernyaev)         << 
2865  *                                            << 
2866  * Function: Create polyhedron for Hyperbolic << 
2867  *                                            << 
2868  * Input: a - half-separation                 << 
2869  *        h - height                          << 
2870  *        r - radius                          << 
2871  *                                            << 
2872  ******************************************** << 
2873 {                                             << 
2874   G4double H = std::abs(h);                   << 
2875   G4double R = std::abs(r);                   << 
2876   G4double A = std::abs(a);                   << 
2877   G4double B = A*R/std::sqrt(2*A*H + H*H);    << 
2878                                               << 
2879   //   P R E P A R E   T W O   P O L Y L I N  << 
2880                                               << 
2881   G4int np1 = (A == 0.) ? 2 : std::max(3, Get << 
2882   G4int np2 = 2;                              << 
2883   G4double maxAng = (A == 0.) ? 0. : std::aco << 
2884   G4double delAng = maxAng/(np1 - 1);         << 
2885                                               << 
2886   auto zz = new G4double[np1 + np2];          << 
2887   auto rr = new G4double[np1 + np2];          << 
2888                                               << 
2889   // 1st polyline                             << 
2890   zz[0] = H;                                  << 
2891   rr[0] = R;                                  << 
2892   for (G4int iz = 1; iz < np1 - 1; ++iz)      << 
2893   {                                           << 
2894     G4double ang = maxAng - iz*delAng;        << 
2895     zz[iz] = A*std::cosh(ang) - A;            << 
2896     rr[iz] = B*std::sinh(ang);                << 
2897   }                                           << 
2898   zz[np1 - 1] = 0.;                           << 
2899   rr[np1 - 1] = 0.;                           << 
2900                                               << 
2901   // 2nd polyline                             << 
2902   zz[np1] = H;                                << 
2903   rr[np1] = 0.;                               << 
2904   zz[np1 + 1] = 0.;                           << 
2905   rr[np1 + 1] = 0.;                           << 
2906                                               << 
2907   //   R O T A T E    P O L Y L I N E S       << 
2908                                               << 
2909   G4double phi  = 0.;                         << 
2910   G4double dphi = CLHEP::twopi;               << 
2911   RotateAroundZ(0, phi, dphi, np1, np2, zz, r << 
2912   SetReferences();                            << 
2913                                               << 
2914   delete [] zz;                                  1591   delete [] zz;
2915   delete [] rr;                                  1592   delete [] rr;
2916 }                                                1593 }
2917                                                  1594 
2918 HepPolyhedronHyperbolicMirror::~HepPolyhedron << 1595 int HepPolyhedron::fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
2919                                               << 
2920 HepPolyhedronTetMesh::                        << 
2921 HepPolyhedronTetMesh(const std::vector<G4Thre << 
2922 /******************************************** << 
2923  *                                            << 
2924  * Name: HepPolyhedronTetMesh                 << 
2925  * Author: E.Tcherniaev (E.Chernyaev)         << 
2926  *                                            << 
2927  * Function: Create polyhedron for tetrahedro << 
2928  *                                            << 
2929  * Input: tetrahedra - array of tetrahedron v << 
2930  *                     per tetrahedron        << 
2931  *                                            << 
2932  ******************************************** << 
2933 {                                             << 
2934   // Check size of input vector               << 
2935   G4int nnodes = (G4int)tetrahedra.size();    << 
2936   if (nnodes == 0)                            << 
2937   {                                           << 
2938     std::cerr                                 << 
2939       << "HepPolyhedronTetMesh: Empty tetrahe << 
2940     return;                                   << 
2941   }                                           << 
2942   G4int ntet = nnodes/4;                      << 
2943   if (nnodes != ntet*4)                       << 
2944   {                                           << 
2945     std::cerr << "HepPolyhedronTetMesh: Numbe << 
2946               << " in tetrahedron mesh is NOT << 
2947               << std::endl;                   << 
2948     return;                                   << 
2949   }                                           << 
2950                                               << 
2951   // Find coincident vertices using hash tabl << 
2952   // This could be done using std::unordered_ << 
2953   // below runs faster.                       << 
2954   std::vector<G4int> iheads(nnodes, -1);      << 
2955   std::vector<std::pair<G4int,G4int>> ipairs( << 
2956   for (G4int i = 0; i < nnodes; ++i)          << 
2957   {                                           << 
2958     // Generate hash key                      << 
2959     G4ThreeVector point = tetrahedra[i];      << 
2960     auto key = std::hash<G4double>()(point.x( << 
2961     key ^= std::hash<G4double>()(point.y());  << 
2962     key ^= std::hash<G4double>()(point.z());  << 
2963     key %= nnodes;                            << 
2964     // Check head of the list                 << 
2965     if (iheads[key] < 0)                      << 
2966     {                                         << 
2967       iheads[key] = i;                        << 
2968       ipairs[i].first = i;                    << 
2969       continue;                               << 
2970     }                                         << 
2971     // Loop along the list                    << 
2972     for (G4int icur = iheads[key], iprev = 0; << 
2973     {                                         << 
2974       G4int icheck = ipairs[icur].first;      << 
2975       if (tetrahedra[icheck] == point)        << 
2976       {                                       << 
2977         ipairs[i].first = icheck; // coincide << 
2978         break;                                << 
2979       }                                       << 
2980       iprev = icur;                           << 
2981       icur = ipairs[icur].second;             << 
2982       // Append vertex to the list            << 
2983       if (icur < 0)                           << 
2984       {                                       << 
2985         ipairs[i].first = i;                  << 
2986         ipairs[iprev].second = i;             << 
2987         break;                                << 
2988       }                                       << 
2989     }                                         << 
2990   }                                           << 
2991                                               << 
2992   // Create vector of original facets         << 
2993   struct facet                                << 
2994   {                                           << 
2995     G4int i1, i2, i3;                         << 
2996     facet() : i1(0), i2(0), i3(0) {};         << 
2997     facet(G4int k1, G4int k2, G4int k3) : i1( << 
2998   };                                          << 
2999   G4int nfacets = nnodes;                     << 
3000   std::vector<facet> ifacets(nfacets);        << 
3001   for (G4int i = 0; i < nfacets; i += 4)      << 
3002   {                                           << 
3003     G4int i0 = ipairs[i + 0].first;           << 
3004     G4int i1 = ipairs[i + 1].first;           << 
3005     G4int i2 = ipairs[i + 2].first;           << 
3006     G4int i3 = ipairs[i + 3].first;           << 
3007     if (i0 > i1) std::swap(i0, i1);           << 
3008     if (i0 > i2) std::swap(i0, i2);           << 
3009     if (i0 > i3) std::swap(i0, i3);           << 
3010     if (i1 > i2) std::swap(i1, i2);           << 
3011     if (i1 > i3) std::swap(i1, i3);           << 
3012     G4ThreeVector e1 = tetrahedra[i1] - tetra << 
3013     G4ThreeVector e2 = tetrahedra[i2] - tetra << 
3014     G4ThreeVector e3 = tetrahedra[i3] - tetra << 
3015     G4double volume = (e1.cross(e2)).dot(e3); << 
3016     if (volume > 0.) std::swap(i2, i3);       << 
3017     ifacets[i + 0] = facet(i0, i1, i2);       << 
3018     ifacets[i + 1] = facet(i0, i2, i3);       << 
3019     ifacets[i + 2] = facet(i0, i3, i1);       << 
3020     ifacets[i + 3] = facet(i1, i3, i2);       << 
3021   }                                           << 
3022                                               << 
3023   // Find shared facets                       << 
3024   std::fill(iheads.begin(), iheads.end(), -1) << 
3025   std::fill(ipairs.begin(), ipairs.end(), std << 
3026   for (G4int i = 0; i < nfacets; ++i)         << 
3027   {                                           << 
3028     // Check head of the list                 << 
3029     G4int key = ifacets[i].i1;                << 
3030     if (iheads[key] < 0)                      << 
3031     {                                         << 
3032       iheads[key] = i;                        << 
3033       ipairs[i].first = i;                    << 
3034       continue;                               << 
3035     }                                         << 
3036     // Loop along the list                    << 
3037     G4int i2 = ifacets[i].i2, i3 = ifacets[i] << 
3038     for (G4int icur = iheads[key], iprev = -1 << 
3039     {                                         << 
3040       G4int icheck = ipairs[icur].first;      << 
3041       if (ifacets[icheck].i2 == i3 && ifacets << 
3042       {                                       << 
3043         if (iprev < 0)                        << 
3044         {                                     << 
3045           iheads[key] = ipairs[icur].second;  << 
3046         }                                     << 
3047         else                                  << 
3048         {                                     << 
3049           ipairs[iprev].second = ipairs[icur] << 
3050         }                                     << 
3051         ipairs[icur].first = -1; // shared fa << 
3052         ipairs[icur].second = -1;             << 
3053         break;                                << 
3054       }                                       << 
3055       iprev = icur;                           << 
3056       icur = ipairs[icur].second;             << 
3057       // Append facet to the list             << 
3058       if (icur < 0)                           << 
3059       {                                       << 
3060         ipairs[i].first = i;                  << 
3061         ipairs[iprev].second = i;             << 
3062         break;                                << 
3063       }                                       << 
3064     }                                         << 
3065   }                                           << 
3066                                               << 
3067   // Count vertices and facets skipping share << 
3068   std::fill(iheads.begin(), iheads.end(), -1) << 
3069   G4int nver = 0, nfac = 0;                   << 
3070   for (G4int i = 0; i < nfacets; ++i)         << 
3071   {                                           << 
3072     if (ipairs[i].first < 0) continue;        << 
3073     G4int i1 = ifacets[i].i1;                 << 
3074     G4int i2 = ifacets[i].i2;                 << 
3075     G4int i3 = ifacets[i].i3;                 << 
3076     if (iheads[i1] < 0) iheads[i1] = nver++;  << 
3077     if (iheads[i2] < 0) iheads[i2] = nver++;  << 
3078     if (iheads[i3] < 0) iheads[i3] = nver++;  << 
3079     nfac++;                                   << 
3080   }                                           << 
3081                                               << 
3082   // Construct polyhedron                     << 
3083   AllocateMemory(nver, nfac);                 << 
3084   for (G4int i = 0; i < nnodes; ++i)          << 
3085   {                                           << 
3086     G4int k = iheads[i];                      << 
3087     if (k >= 0) SetVertex(k + 1, tetrahedra[i << 
3088   }                                           << 
3089   for (G4int i = 0, k = 0; i < nfacets; ++i)  << 
3090   {                                           << 
3091     if (ipairs[i].first < 0) continue;        << 
3092     G4int i1 = iheads[ifacets[i].i1] + 1;     << 
3093     G4int i2 = iheads[ifacets[i].i2] + 1;     << 
3094     G4int i3 = iheads[ifacets[i].i3] + 1;     << 
3095     SetFacet(++k, i1, i2, i3);                << 
3096   }                                           << 
3097   SetReferences();                            << 
3098 }                                             << 
3099                                               << 
3100 HepPolyhedronTetMesh::~HepPolyhedronTetMesh() << 
3101                                               << 
3102 HepPolyhedronBoxMesh::                        << 
3103 HepPolyhedronBoxMesh(G4double sizeX, G4double << 
3104                      const std::vector<G4Thre << 
3105 /******************************************** << 
3106  *                                            << 
3107  * Name: HepPolyhedronBoxMesh                 << 
3108  * Author: E.Tcherniaev (E.Chernyaev)         << 
3109  *                                            << 
3110  * Function: Create polyhedron for box mesh   << 
3111  *                                            << 
3112  * Input: sizeX, sizeY, sizeZ - dimensions of << 
3113  *        positions - vector of cell centres  << 
3114  *                                            << 
3115  ******************************************** << 
3116 {                                             << 
3117   G4int nbox = (G4int)positions.size();       << 
3118   if (nbox == 0)                              << 
3119   {                                           << 
3120     std::cerr << "HepPolyhedronBoxMesh: Empty << 
3121     return;                                   << 
3122   }                                           << 
3123   // compute inverse dimensions               << 
3124   G4double invx = 1./sizeX, invy = 1./sizeY,  << 
3125   // find mesh bounding box                   << 
3126   G4ThreeVector pmin = positions[0], pmax = p << 
3127   for (const auto& p: positions)              << 
3128   {                                           << 
3129     if (pmin.x() > p.x()) pmin.setX(p.x());   << 
3130     if (pmin.y() > p.y()) pmin.setY(p.y());   << 
3131     if (pmin.z() > p.z()) pmin.setZ(p.z());   << 
3132     if (pmax.x() < p.x()) pmax.setX(p.x());   << 
3133     if (pmax.y() < p.y()) pmax.setY(p.y());   << 
3134     if (pmax.z() < p.z()) pmax.setZ(p.z());   << 
3135   }                                           << 
3136   // find number of voxels                    << 
3137   G4int nx = (pmax.x() - pmin.x())*invx + 1.5 << 
3138   G4int ny = (pmax.y() - pmin.y())*invy + 1.5 << 
3139   G4int nz = (pmax.z() - pmin.z())*invz + 1.5 << 
3140   // create structures for voxels and node in << 
3141   std::vector<char> voxels(nx*ny*nz, 0);      << 
3142   std::vector<G4int> indices((nx+1)*(ny+1)*(n << 
3143   // mark voxels listed in positions          << 
3144   G4int kx =  ny*nz, ky = nz;                 << 
3145   for (const auto& p: positions)              << 
3146   {                                           << 
3147     G4int ix = (p.x() - pmin.x())*invx + 0.5; << 
3148     G4int iy = (p.y() - pmin.y())*invy + 0.5; << 
3149     G4int iz = (p.z() - pmin.z())*invz + 0.5; << 
3150     G4int i = ix*kx + iy*ky + iz;             << 
3151     voxels[i] = 1;                            << 
3152   }                                           << 
3153   // count number of vertices and facets      << 
3154   // set indices                              << 
3155   G4int kvx = (ny + 1)*(nz + 1), kvy = nz + 1 << 
3156   G4int nver = 0, nfac = 0;                   << 
3157   for (const auto& p: positions)              << 
3158   {                                           << 
3159     G4int ix = (p.x() - pmin.x())*invx + 0.5; << 
3160     G4int iy = (p.y() - pmin.y())*invy + 0.5; << 
3161     G4int iz = (p.z() - pmin.z())*invz + 0.5; << 
3162     //                                        << 
3163     //    011       111                       << 
3164     //      +---–---+                       << 
3165     //      | 001   |   101                   << 
3166     //      |   +---–---+                   << 
3167     //      |   |   |   |                     << 
3168     //      +---|---+   |                     << 
3169     //    010   |   110 |                     << 
3170     //          +-------+                     << 
3171     //        000       100                   << 
3172     //                                        << 
3173     G4int vcheck = 0;                         << 
3174     // check (ix - 1) side                    << 
3175     vcheck = (ix == 0) ? 0 : voxels[(ix-1)*kx << 
3176     if (vcheck == 0)                          << 
3177     {                                         << 
3178       nfac++;                                 << 
3179       G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3180       G4int i2 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3181       G4int i3 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3182       G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3183       if (indices[i1] == 0) indices[i1] = ++n << 
3184       if (indices[i2] == 0) indices[i2] = ++n << 
3185       if (indices[i3] == 0) indices[i3] = ++n << 
3186       if (indices[i4] == 0) indices[i4] = ++n << 
3187     }                                         << 
3188     // check (ix + 1) side                    << 
3189     vcheck = (ix == nx - 1) ? 0 : voxels[(ix+ << 
3190     if (vcheck == 0)                          << 
3191     {                                         << 
3192       nfac++;                                 << 
3193       G4int i1 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3194       G4int i2 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3195       G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3196       G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3197       if (indices[i1] == 0) indices[i1] = ++n << 
3198       if (indices[i2] == 0) indices[i2] = ++n << 
3199       if (indices[i3] == 0) indices[i3] = ++n << 
3200       if (indices[i4] == 0) indices[i4] = ++n << 
3201     }                                         << 
3202     // check (iy - 1) side                    << 
3203     vcheck = (iy == 0) ? 0 : voxels[ix*kx + ( << 
3204     if (vcheck == 0)                          << 
3205     {                                         << 
3206       nfac++;                                 << 
3207       G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3208       G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3209       G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3210       G4int i4 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3211       if (indices[i1] == 0) indices[i1] = ++n << 
3212       if (indices[i2] == 0) indices[i2] = ++n << 
3213       if (indices[i3] == 0) indices[i3] = ++n << 
3214       if (indices[i4] == 0) indices[i4] = ++n << 
3215     }                                         << 
3216     // check (iy + 1) side                    << 
3217     vcheck = (iy == ny - 1) ? 0 : voxels[ix*k << 
3218     if (vcheck == 0)                          << 
3219     {                                         << 
3220       nfac++;                                 << 
3221       G4int i1 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3222       G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3223       G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3224       G4int i4 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3225       if (indices[i1] == 0) indices[i1] = ++n << 
3226       if (indices[i2] == 0) indices[i2] = ++n << 
3227       if (indices[i3] == 0) indices[i3] = ++n << 
3228       if (indices[i4] == 0) indices[i4] = ++n << 
3229     }                                         << 
3230     // check (iz - 1) side                    << 
3231     vcheck = (iz == 0) ? 0 : voxels[ix*kx + i << 
3232     if (vcheck == 0)                          << 
3233     {                                         << 
3234       nfac++;                                 << 
3235       G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3236       G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3237       G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3238       G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3239       if (indices[i1] == 0) indices[i1] = ++n << 
3240       if (indices[i2] == 0) indices[i2] = ++n << 
3241       if (indices[i3] == 0) indices[i3] = ++n << 
3242       if (indices[i4] == 0) indices[i4] = ++n << 
3243     }                                         << 
3244     // check (iz + 1) side                    << 
3245     vcheck = (iz == nz - 1) ? 0 : voxels[ix*k << 
3246     if (vcheck == 0)                          << 
3247     {                                         << 
3248       nfac++;                                 << 
3249       G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3250       G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3251       G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3252       G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3253       if (indices[i1] == 0) indices[i1] = ++n << 
3254       if (indices[i2] == 0) indices[i2] = ++n << 
3255       if (indices[i3] == 0) indices[i3] = ++n << 
3256       if (indices[i4] == 0) indices[i4] = ++n << 
3257     }                                         << 
3258   }                                           << 
3259   // Construct polyhedron                     << 
3260   AllocateMemory(nver, nfac);                 << 
3261   G4ThreeVector p0(pmin.x() - 0.5*sizeX, pmin << 
3262   for (G4int ix = 0; ix <= nx; ++ix)          << 
3263   {                                           << 
3264     for (G4int iy = 0; iy <= ny; ++iy)        << 
3265     {                                         << 
3266       for (G4int iz = 0; iz <= nz; ++iz)      << 
3267       {                                       << 
3268   G4int i = ix*kvx + iy*kvy + iz;             << 
3269   if (indices[i] == 0) continue;              << 
3270   SetVertex(indices[i], p0 + G4ThreeVector(ix << 
3271       }                                       << 
3272     }                                         << 
3273   }                                           << 
3274   nfac = 0;                                   << 
3275   for (const auto& p: positions)              << 
3276   {                                           << 
3277     G4int ix = (p.x() - pmin.x())*invx + 0.5; << 
3278     G4int iy = (p.y() - pmin.y())*invy + 0.5; << 
3279     G4int iz = (p.z() - pmin.z())*invz + 0.5; << 
3280     G4int vcheck = 0;                         << 
3281     // check (ix - 1) side                    << 
3282     vcheck = (ix == 0) ? 0 : voxels[(ix-1)*kx << 
3283     if (vcheck == 0)                          << 
3284     {                                         << 
3285       G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3286       G4int i2 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3287       G4int i3 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3288       G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3289       SetFacet(++nfac, indices[i1], indices[i << 
3290     }                                         << 
3291     // check (ix + 1) side                    << 
3292     vcheck = (ix == nx - 1) ? 0 : voxels[(ix+ << 
3293     if (vcheck == 0)                          << 
3294     {                                         << 
3295       G4int i1 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3296       G4int i2 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3297       G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3298       G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3299       SetFacet(++nfac, indices[i1], indices[i << 
3300                                               << 
3301     }                                         << 
3302     // check (iy - 1) side                    << 
3303     vcheck = (iy == 0) ? 0 : voxels[ix*kx + ( << 
3304     if (vcheck == 0)                          << 
3305     {                                         << 
3306       G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3307       G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3308       G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3309       G4int i4 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3310       SetFacet(++nfac, indices[i1], indices[i << 
3311     }                                         << 
3312     // check (iy + 1) side                    << 
3313     vcheck = (iy == ny - 1) ? 0 : voxels[ix*k << 
3314     if (vcheck == 0)                          << 
3315     {                                         << 
3316       G4int i1 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3317       G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3318       G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3319       G4int i4 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3320       SetFacet(++nfac, indices[i1], indices[i << 
3321     }                                         << 
3322     // check (iz - 1) side                    << 
3323     vcheck = (iz == 0) ? 0 : voxels[ix*kx + i << 
3324     if (vcheck == 0)                          << 
3325     {                                         << 
3326       G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3327       G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3328       G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3329       G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3330       SetFacet(++nfac, indices[i1], indices[i << 
3331     }                                         << 
3332     // check (iz + 1) side                    << 
3333     vcheck = (iz == nz - 1) ? 0 : voxels[ix*k << 
3334     if (vcheck == 0)                          << 
3335     {                                         << 
3336       G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (i << 
3337       G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (i << 
3338       G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (i << 
3339       G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (i << 
3340       SetFacet(++nfac, indices[i1], indices[i << 
3341     }                                         << 
3342   }                                           << 
3343   SetReferences();                            << 
3344 }                                             << 
3345                                               << 
3346 HepPolyhedronBoxMesh::~HepPolyhedronBoxMesh() << 
3347                                               << 
3348 G4ThreadLocal                                 << 
3349 G4int HepPolyhedron::fNumberOfRotationSteps = << 
3350 /********************************************    1596 /***********************************************************************
3351  *                                               1597  *                                                                     *
3352  * Name: HepPolyhedron::fNumberOfRotationStep    1598  * Name: HepPolyhedron::fNumberOfRotationSteps       Date:    24.06.97 *
3353  * Author: J.Allison (Manchester University)     1599  * Author: J.Allison (Manchester University)         Revised:          *
3354  *                                               1600  *                                                                     *
3355  * Function: Number of steps for whole circle    1601  * Function: Number of steps for whole circle                          *
3356  *                                               1602  *                                                                     *
3357  ********************************************    1603  ***********************************************************************/
3358                                                  1604 
3359 #include "BooleanProcessor.src"                  1605 #include "BooleanProcessor.src"
                                                   >> 1606 static BooleanProcessor processor;
3360                                                  1607 
3361 HepPolyhedron HepPolyhedron::add(const HepPol << 1608 HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const 
3362 /********************************************    1609 /***********************************************************************
3363  *                                               1610  *                                                                     *
3364  * Name: HepPolyhedron::add                      1611  * Name: HepPolyhedron::add                          Date:    19.03.00 *
3365  * Author: E.Chernyaev                           1612  * Author: E.Chernyaev                               Revised:          *
3366  *                                               1613  *                                                                     *
3367  * Function: Boolean "union" of two polyhedra    1614  * Function: Boolean "union" of two polyhedra                          *
3368  *                                               1615  *                                                                     *
3369  ********************************************    1616  ***********************************************************************/
3370 {                                                1617 {
3371   G4int ierr;                                 << 1618   return processor.execute(OP_UNION, *this, p);
3372   BooleanProcessor processor;                 << 
3373   return processor.execute(OP_UNION, *this, p << 
3374 }                                                1619 }
3375                                                  1620 
3376 HepPolyhedron HepPolyhedron::intersect(const  << 1621 HepPolyhedron HepPolyhedron::intersect(const HepPolyhedron & p) const 
3377 /********************************************    1622 /***********************************************************************
3378  *                                               1623  *                                                                     *
3379  * Name: HepPolyhedron::intersect                1624  * Name: HepPolyhedron::intersect                    Date:    19.03.00 *
3380  * Author: E.Chernyaev                           1625  * Author: E.Chernyaev                               Revised:          *
3381  *                                               1626  *                                                                     *
3382  * Function: Boolean "intersection" of two po    1627  * Function: Boolean "intersection" of two polyhedra                   *
3383  *                                               1628  *                                                                     *
3384  ********************************************    1629  ***********************************************************************/
3385 {                                                1630 {
3386   G4int ierr;                                 << 1631   return processor.execute(OP_INTERSECTION, *this, p);
3387   BooleanProcessor processor;                 << 
3388   return processor.execute(OP_INTERSECTION, * << 
3389 }                                                1632 }
3390                                                  1633 
3391 HepPolyhedron HepPolyhedron::subtract(const H << 1634 HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const 
3392 /********************************************    1635 /***********************************************************************
3393  *                                               1636  *                                                                     *
3394  * Name: HepPolyhedron::add                      1637  * Name: HepPolyhedron::add                          Date:    19.03.00 *
3395  * Author: E.Chernyaev                           1638  * Author: E.Chernyaev                               Revised:          *
3396  *                                               1639  *                                                                     *
3397  * Function: Boolean "subtraction" of "p" fro    1640  * Function: Boolean "subtraction" of "p" from "this"                  *
3398  *                                               1641  *                                                                     *
3399  ********************************************    1642  ***********************************************************************/
3400 {                                                1643 {
3401   G4int ierr;                                 << 1644   return processor.execute(OP_SUBTRACTION, *this, p);
3402   BooleanProcessor processor;                 << 
3403   return processor.execute(OP_SUBTRACTION, *t << 
3404 }                                                1645 }
3405                                               << 
3406 //NOTE : include the code of HepPolyhedronPro << 
3407 //       since there is no BooleanProcessor.h << 
3408                                               << 
3409 #undef INTERSECTION                           << 
3410                                               << 
3411 #include "HepPolyhedronProcessor.src"         << 
3412                                                  1646