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


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 //
                                                   >>  27 // $Id: HepPolyhedron.cc,v 1.23.2.1 2007/08/21 14:04:32 gcosmo Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-00-patch-01 $
                                                   >>  29 //
                                                   >>  30 // 
                                                   >>  31 //
 26 // G4 Polyhedron library                           32 // G4 Polyhedron library
 27 //                                                 33 //
 28 // History:                                        34 // History:
 29 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@ce     35 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
 30 //                                                 36 //
 31 // 30.09.96 E.Chernyaev                            37 // 30.09.96 E.Chernyaev
 32 // - added GetNextVertexIndex, GetVertex by Ya     38 // - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
 33 // - added GetNextUnitNormal, GetNextEdgeIndic <<  39 // - added GetNextUnitNormal, GetNextEdgeIndeces, GetNextEdge
 34 //                                                 40 //
 35 // 15.12.96 E.Chernyaev                            41 // 15.12.96 E.Chernyaev
 36 // - added GetNumberOfRotationSteps, RotateEdg     42 // - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
 37 // - rewritten G4PolyhedronCons;                   43 // - rewritten G4PolyhedronCons;
 38 // - added G4PolyhedronPara, ...Trap, ...Pgon,     44 // - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
 39 //                                                 45 //
 40 // 01.06.97 E.Chernyaev                            46 // 01.06.97 E.Chernyaev
 41 // - modified RotateAroundZ, added SetSideFace     47 // - modified RotateAroundZ, added SetSideFacets
 42 //                                                 48 //
 43 // 19.03.00 E.Chernyaev                            49 // 19.03.00 E.Chernyaev
 44 // - implemented boolean operations (add, subt     50 // - implemented boolean operations (add, subtract, intersect) on polyhedra;
 45 //                                                 51 //
 46 // 25.05.01 E.Chernyaev                            52 // 25.05.01 E.Chernyaev
 47 // - added GetSurfaceArea() and GetVolume()    <<  53 // - added GetSurfaceArea() and GetVolume();
 48 //                                                 54 //
 49 // 05.11.02 E.Chernyaev                            55 // 05.11.02 E.Chernyaev
 50 // - added createTwistedTrap() and createPolyh <<  56 // - added createTwistedTrap() and createPolyhedron();
 51 //                                                 57 //
 52 // 20.06.05 G.Cosmo                                58 // 20.06.05 G.Cosmo
 53 // - added HepPolyhedronEllipsoid              <<  59 // - added HepPolyhedronEllipsoid;
 54 //                                             << 
 55 // 18.07.07 T.Nikitina                         << 
 56 // - added HepPolyhedronParaboloid             << 
 57 //                                             << 
 58 // 22.02.20 E.Chernyaev                        << 
 59 // - added HepPolyhedronTet, HepPolyhedronHybe << 
 60 //                                                 60 //
 61 // 12.05.21 E.Chernyaev                        <<  61   
 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"                         62 #include "HepPolyhedron.h"
 76 #include "G4PhysicalConstants.hh"              <<  63 #include <CLHEP/Units/PhysicalConstants.h>
 77 #include "G4Vector3D.hh"                       <<  64 #include <CLHEP/Geometry/Vector3D.h>
 78                                                    65 
 79 #include <cstdlib>  // Required on some compil << 
 80 #include <cmath>                                   66 #include <cmath>
 81 #include <algorithm>                           << 
 82                                                    67 
                                                   >>  68 using namespace HepGeom;
 83 using CLHEP::perMillion;                           69 using CLHEP::perMillion;
 84 using CLHEP::deg;                                  70 using CLHEP::deg;
 85 using CLHEP::pi;                                   71 using CLHEP::pi;
 86 using CLHEP::twopi;                                72 using CLHEP::twopi;
 87 using CLHEP::nm;                               << 
 88 const G4double spatialTolerance = 0.01*nm;     << 
 89                                                    73 
 90 /*********************************************     74 /***********************************************************************
 91  *                                                 75  *                                                                     *
 92  * Name: HepPolyhedron operator <<                 76  * Name: HepPolyhedron operator <<                   Date:    09.05.96 *
 93  * Author: E.Chernyaev (IHEP/Protvino)             77  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
 94  *                                                 78  *                                                                     *
 95  * Function: Print contents of G4 polyhedron       79  * Function: Print contents of G4 polyhedron                           *
 96  *                                                 80  *                                                                     *
 97  *********************************************     81  ***********************************************************************/
 98 std::ostream & operator<<(std::ostream & ostr,     82 std::ostream & operator<<(std::ostream & ostr, const G4Facet & facet) {
 99   for (const auto& edge : facet.edge) {        <<  83   for (int k=0; k<4; k++) {
100     ostr << " " << edge.v << "/" << edge.f;    <<  84     ostr << " " << facet.edge[k].v << "/" << facet.edge[k].f;
101   }                                                85   }
102   return ostr;                                     86   return ostr;
103 }                                                  87 }
104                                                    88 
105 std::ostream & operator<<(std::ostream & ostr,     89 std::ostream & operator<<(std::ostream & ostr, const HepPolyhedron & ph) {
106   ostr << std::endl;                               90   ostr << std::endl;
107   ostr << "Nvertices=" << ph.nvert << ", Nface <<  91   ostr << "Nverteces=" << ph.nvert << ", Nfacets=" << ph.nface << std::endl;
108   G4int i;                                     <<  92   int i;
109   for (i=1; i<=ph.nvert; i++) {                    93   for (i=1; i<=ph.nvert; i++) {
110      ostr << "xyz(" << i << ")="                   94      ostr << "xyz(" << i << ")="
111           << ph.pV[i].x() << ' ' << ph.pV[i].y     95           << ph.pV[i].x() << ' ' << ph.pV[i].y() << ' ' << ph.pV[i].z()
112           << std::endl;                            96           << std::endl;
113   }                                                97   }
114   for (i=1; i<=ph.nface; i++) {                    98   for (i=1; i<=ph.nface; i++) {
115     ostr << "face(" << i << ")=" << ph.pF[i] <     99     ostr << "face(" << i << ")=" << ph.pF[i] << std::endl;
116   }                                               100   }
117   return ostr;                                    101   return ostr;
118 }                                                 102 }
119                                                   103 
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    104 HepPolyhedron::HepPolyhedron(const HepPolyhedron &from)
134 /*********************************************    105 /***********************************************************************
135  *                                                106  *                                                                     *
136  * Name: HepPolyhedron copy constructor           107  * Name: HepPolyhedron copy constructor             Date:    23.07.96  *
137  * Author: E.Chernyaev (IHEP/Protvino)            108  * Author: E.Chernyaev (IHEP/Protvino)              Revised:           *
138  *                                                109  *                                                                     *
139  *********************************************    110  ***********************************************************************/
140 : nvert(0), nface(0), pV(nullptr), pF(nullptr) << 111 : nvert(0), nface(0), pV(0), pF(0)
141 {                                                 112 {
142   AllocateMemory(from.nvert, from.nface);         113   AllocateMemory(from.nvert, from.nface);
143   for (G4int i=1; i<=nvert; i++) pV[i] = from. << 114   for (int i=1; i<=nvert; i++) pV[i] = from.pV[i];
144   for (G4int k=1; k<=nface; k++) pF[k] = from. << 115   for (int k=1; k<=nface; k++) pF[k] = from.pF[k];
145 }                                              << 
146                                                << 
147 HepPolyhedron::HepPolyhedron(HepPolyhedron&& f << 
148 /********************************************* << 
149  *                                             << 
150  * Name: HepPolyhedron move constructor        << 
151  * Author: E.Tcherniaev (E.Chernyaev)          << 
152  *                                             << 
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 }                                                 116 }
167                                                   117 
168 HepPolyhedron & HepPolyhedron::operator=(const    118 HepPolyhedron & HepPolyhedron::operator=(const HepPolyhedron &from)
169 /*********************************************    119 /***********************************************************************
170  *                                                120  *                                                                     *
171  * Name: HepPolyhedron operator =                 121  * Name: HepPolyhedron operator =                   Date:    23.07.96  *
172  * Author: E.Chernyaev (IHEP/Protvino)            122  * Author: E.Chernyaev (IHEP/Protvino)              Revised:           *
173  *                                                123  *                                                                     *
174  * Function: Copy contents of one polyhedron t    124  * Function: Copy contents of one polyhedron to another                *
175  *                                                125  *                                                                     *
176  *********************************************    126  ***********************************************************************/
177 {                                                 127 {
178   if (this != &from) {                            128   if (this != &from) {
179     AllocateMemory(from.nvert, from.nface);       129     AllocateMemory(from.nvert, from.nface);
180     for (G4int i=1; i<=nvert; i++) pV[i] = fro << 130     for (int i=1; i<=nvert; i++) pV[i] = from.pV[i];
181     for (G4int k=1; k<=nface; k++) pF[k] = fro << 131     for (int k=1; k<=nface; k++) pF[k] = from.pF[k];
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  *                                             << 
194  ********************************************* << 
195 {                                              << 
196   if (this != &from) {                         << 
197     delete [] pV;                              << 
198     delete [] pF;                              << 
199     nvert = from.nvert;                        << 
200     nface = from.nface;                        << 
201     pV = from.pV;                              << 
202     pF = from.pF;                              << 
203                                                << 
204     // Release the data from the source object << 
205     from.nvert = 0;                            << 
206     from.nface = 0;                            << 
207     from.pV = nullptr;                         << 
208     from.pF = nullptr;                         << 
209   }                                               132   }
210   return *this;                                   133   return *this;
211 }                                                 134 }
212                                                   135 
213 G4int                                          << 136 int
214 HepPolyhedron::FindNeighbour(G4int iFace, G4in << 137 HepPolyhedron::FindNeighbour(int iFace, int iNode, int iOrder) const
215 /*********************************************    138 /***********************************************************************
216  *                                                139  *                                                                     *
217  * Name: HepPolyhedron::FindNeighbour             140  * Name: HepPolyhedron::FindNeighbour                Date:    22.11.99 *
218  * Author: E.Chernyaev                            141  * Author: E.Chernyaev                               Revised:          *
219  *                                                142  *                                                                     *
220  * Function: Find neighbouring face               143  * Function: Find neighbouring face                                    *
221  *                                                144  *                                                                     *
222  *********************************************    145  ***********************************************************************/
223 {                                                 146 {
224   G4int i;                                     << 147   int i;
225   for (i=0; i<4; i++) {                           148   for (i=0; i<4; i++) {
226     if (iNode == std::abs(pF[iFace].edge[i].v)    149     if (iNode == std::abs(pF[iFace].edge[i].v)) break;
227   }                                               150   }
228   if (i == 4) {                                   151   if (i == 4) {
229     std::cerr                                     152     std::cerr
230       << "HepPolyhedron::FindNeighbour: face "    153       << "HepPolyhedron::FindNeighbour: face " << iFace
231       << " has no node " << iNode                 154       << " has no node " << iNode
232       << std::endl;                            << 155       << std::endl; 
233     return 0;                                     156     return 0;
234   }                                               157   }
235   if (iOrder < 0) {                               158   if (iOrder < 0) {
236     if ( --i < 0) i = 3;                          159     if ( --i < 0) i = 3;
237     if (pF[iFace].edge[i].v == 0) i = 2;          160     if (pF[iFace].edge[i].v == 0) i = 2;
238   }                                               161   }
239   return (pF[iFace].edge[i].v > 0) ? 0 : pF[iF    162   return (pF[iFace].edge[i].v > 0) ? 0 : pF[iFace].edge[i].f;
240 }                                                 163 }
241                                                   164 
242 G4Normal3D HepPolyhedron::FindNodeNormal(G4int << 165 Normal3D<double> HepPolyhedron::FindNodeNormal(int iFace, int iNode) const
243 /*********************************************    166 /***********************************************************************
244  *                                                167  *                                                                     *
245  * Name: HepPolyhedron::FindNodeNormal            168  * Name: HepPolyhedron::FindNodeNormal               Date:    22.11.99 *
246  * Author: E.Chernyaev                            169  * Author: E.Chernyaev                               Revised:          *
247  *                                                170  *                                                                     *
248  * Function: Find normal at given node            171  * Function: Find normal at given node                                 *
249  *                                                172  *                                                                     *
250  *********************************************    173  ***********************************************************************/
251 {                                                 174 {
252   G4Normal3D normal = GetUnitNormal(iFace);    << 175   Normal3D<double>  normal = GetUnitNormal(iFace);
253   G4int      k = iFace, iOrder = 1;            << 176   int          k = iFace, iOrder = 1, n = 1;
254                                                   177 
255   for(;;) {                                       178   for(;;) {
256     k = FindNeighbour(k, iNode, iOrder);          179     k = FindNeighbour(k, iNode, iOrder);
257     if (k == iFace) break;                     << 180     if (k == iFace) break; 
258     if (k > 0) {                                  181     if (k > 0) {
                                                   >> 182       n++;
259       normal += GetUnitNormal(k);                 183       normal += GetUnitNormal(k);
260     }else{                                        184     }else{
261       if (iOrder < 0) break;                      185       if (iOrder < 0) break;
262       k = iFace;                                  186       k = iFace;
263       iOrder = -iOrder;                           187       iOrder = -iOrder;
264     }                                             188     }
265   }                                               189   }
266   return normal.unit();                           190   return normal.unit();
267 }                                                 191 }
268                                                   192 
269 G4int HepPolyhedron::GetNumberOfRotationSteps( << 193 int HepPolyhedron::GetNumberOfRotationSteps()
270 /*********************************************    194 /***********************************************************************
271  *                                                195  *                                                                     *
272  * Name: HepPolyhedron::GetNumberOfRotationSte    196  * Name: HepPolyhedron::GetNumberOfRotationSteps     Date:    24.06.97 *
273  * Author: J.Allison (Manchester University)      197  * Author: J.Allison (Manchester University)         Revised:          *
274  *                                                198  *                                                                     *
275  * Function: Get number of steps for whole cir    199  * Function: Get number of steps for whole circle                      *
276  *                                                200  *                                                                     *
277  *********************************************    201  ***********************************************************************/
278 {                                                 202 {
279   return fNumberOfRotationSteps;                  203   return fNumberOfRotationSteps;
280 }                                                 204 }
281                                                   205 
282 void HepPolyhedron::SetVertex(G4int index, con << 206 void HepPolyhedron::SetNumberOfRotationSteps(int n)
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 /*********************************************    207 /***********************************************************************
341  *                                                208  *                                                                     *
342  * Name: HepPolyhedron::SetNumberOfRotationSte    209  * Name: HepPolyhedron::SetNumberOfRotationSteps     Date:    24.06.97 *
343  * Author: J.Allison (Manchester University)      210  * Author: J.Allison (Manchester University)         Revised:          *
344  *                                                211  *                                                                     *
345  * Function: Set number of steps for whole cir    212  * Function: Set number of steps for whole circle                      *
346  *                                                213  *                                                                     *
347  *********************************************    214  ***********************************************************************/
348 {                                                 215 {
349   const G4int nMin = 3;                        << 216   const int nMin = 3;
350   if (n < nMin) {                                 217   if (n < nMin) {
351     std::cerr                                  << 218     std::cerr 
352       << "HepPolyhedron::SetNumberOfRotationSt    219       << "HepPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
353       << "number of steps per circle < " << nM    220       << "number of steps per circle < " << nMin << "; forced to " << nMin
354       << std::endl;                               221       << std::endl;
355     fNumberOfRotationSteps = nMin;                222     fNumberOfRotationSteps = nMin;
356   }else{                                          223   }else{
357     fNumberOfRotationSteps = n;                   224     fNumberOfRotationSteps = n;
358   }                                            << 225   }    
359 }                                                 226 }
360                                                   227 
361 void HepPolyhedron::ResetNumberOfRotationSteps    228 void HepPolyhedron::ResetNumberOfRotationSteps()
362 /*********************************************    229 /***********************************************************************
363  *                                                230  *                                                                     *
364  * Name: HepPolyhedron::GetNumberOfRotationSte    231  * Name: HepPolyhedron::GetNumberOfRotationSteps     Date:    24.06.97 *
365  * Author: J.Allison (Manchester University)      232  * Author: J.Allison (Manchester University)         Revised:          *
366  *                                                233  *                                                                     *
367  * Function: Reset number of steps for whole c    234  * Function: Reset number of steps for whole circle to default value   *
368  *                                                235  *                                                                     *
369  *********************************************    236  ***********************************************************************/
370 {                                                 237 {
371   fNumberOfRotationSteps = DEFAULT_NUMBER_OF_S    238   fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
372 }                                                 239 }
373                                                   240 
374 void HepPolyhedron::AllocateMemory(G4int Nvert << 241 void HepPolyhedron::AllocateMemory(int Nvert, int Nface)
375 /*********************************************    242 /***********************************************************************
376  *                                                243  *                                                                     *
377  * Name: HepPolyhedron::AllocateMemory            244  * Name: HepPolyhedron::AllocateMemory               Date:    19.06.96 *
378  * Author: E.Chernyaev (IHEP/Protvino)            245  * Author: E.Chernyaev (IHEP/Protvino)               Revised: 05.11.02 *
379  *                                                246  *                                                                     *
380  * Function: Allocate memory for GEANT4 polyhe    247  * Function: Allocate memory for GEANT4 polyhedron                     *
381  *                                                248  *                                                                     *
382  * Input: Nvert - number of nodes                 249  * Input: Nvert - number of nodes                                      *
383  *        Nface - number of faces                 250  *        Nface - number of faces                                      *
384  *                                                251  *                                                                     *
385  *********************************************    252  ***********************************************************************/
386 {                                                 253 {
387   if (nvert == Nvert && nface == Nface) return    254   if (nvert == Nvert && nface == Nface) return;
388   delete [] pV;                                << 255   if (pV != 0) delete [] pV;
389   delete [] pF;                                << 256   if (pF != 0) delete [] pF;
390   if (Nvert > 0 && Nface > 0) {                   257   if (Nvert > 0 && Nface > 0) {
391     nvert = Nvert;                                258     nvert = Nvert;
392     nface = Nface;                                259     nface = Nface;
393     pV    = new G4Point3D[nvert+1];            << 260     pV    = new Point3D<double>[nvert+1];
394     pF    = new G4Facet[nface+1];                 261     pF    = new G4Facet[nface+1];
395   }else{                                          262   }else{
396     nvert = 0; nface = 0; pV = nullptr; pF = n << 263     nvert = 0; nface = 0; pV = 0; pF = 0;
397   }                                               264   }
398 }                                                 265 }
399                                                   266 
400 void HepPolyhedron::CreatePrism()                 267 void HepPolyhedron::CreatePrism()
401 /*********************************************    268 /***********************************************************************
402  *                                                269  *                                                                     *
403  * Name: HepPolyhedron::CreatePrism               270  * Name: HepPolyhedron::CreatePrism                  Date:    15.07.96 *
404  * Author: E.Chernyaev (IHEP/Protvino)            271  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
405  *                                                272  *                                                                     *
406  * Function: Set facets for a prism               273  * Function: Set facets for a prism                                    *
407  *                                                274  *                                                                     *
408  *********************************************    275  ***********************************************************************/
409 {                                                 276 {
410   enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRON    277   enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
411                                                   278 
412   pF[1] = G4Facet(1,LEFT,  4,BACK,  3,RIGHT,      279   pF[1] = G4Facet(1,LEFT,  4,BACK,  3,RIGHT,  2,FRONT);
413   pF[2] = G4Facet(5,TOP,   8,BACK,  4,BOTTOM,     280   pF[2] = G4Facet(5,TOP,   8,BACK,  4,BOTTOM, 1,FRONT);
414   pF[3] = G4Facet(8,TOP,   7,RIGHT, 3,BOTTOM,     281   pF[3] = G4Facet(8,TOP,   7,RIGHT, 3,BOTTOM, 4,LEFT);
415   pF[4] = G4Facet(7,TOP,   6,FRONT, 2,BOTTOM,     282   pF[4] = G4Facet(7,TOP,   6,FRONT, 2,BOTTOM, 3,BACK);
416   pF[5] = G4Facet(6,TOP,   5,LEFT,  1,BOTTOM,     283   pF[5] = G4Facet(6,TOP,   5,LEFT,  1,BOTTOM, 2,RIGHT);
417   pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK,       284   pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK,   8,LEFT);
418 }                                                 285 }
419                                                   286 
420 void HepPolyhedron::RotateEdge(G4int k1, G4int << 287 void HepPolyhedron::RotateEdge(int k1, int k2, double r1, double r2,
421                               G4int v1, G4int  << 288                               int v1, int v2, int vEdge,
422                               G4bool ifWholeCi << 289                               bool ifWholeCircle, int ns, int &kface)
423 /*********************************************    290 /***********************************************************************
424  *                                                291  *                                                                     *
425  * Name: HepPolyhedron::RotateEdge                292  * Name: HepPolyhedron::RotateEdge                   Date:    05.12.96 *
426  * Author: E.Chernyaev (IHEP/Protvino)            293  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
427  *                                                294  *                                                                     *
428  * Function: Create set of facets by rotation     295  * Function: Create set of facets by rotation of an edge around Z-axis *
429  *                                                296  *                                                                     *
430  * Input: k1, k2 - end vertices of the edge       297  * Input: k1, k2 - end vertices of the edge                            *
431  *        r1, r2 - radiuses of the end vertice    298  *        r1, r2 - radiuses of the end vertices                        *
432  *        v1, v2 - visibility of edges produce    299  *        v1, v2 - visibility of edges produced by rotation of the end *
433  *                 vertices                       300  *                 vertices                                            *
434  *        vEdge  - visibility of the edge         301  *        vEdge  - visibility of the edge                              *
435  *        ifWholeCircle - is true in case of w    302  *        ifWholeCircle - is true in case of whole circle rotation     *
436  *        nds    - number of discrete steps    << 303  *        ns     - number of discrete steps                            *
437  *        r[]    - r-coordinates                  304  *        r[]    - r-coordinates                                       *
438  *        kface  - current free cell in the pF    305  *        kface  - current free cell in the pF array                   *
439  *                                                306  *                                                                     *
440  *********************************************    307  ***********************************************************************/
441 {                                                 308 {
442   if (r1 == 0. && r2 == 0.) return;            << 309   if (r1 == 0. && r2 == 0) return;
443                                                   310 
444   G4int i;                                     << 311   int i;
445   G4int i1  = k1;                              << 312   int i1  = k1;
446   G4int i2  = k2;                              << 313   int i2  = k2;
447   G4int ii1 = ifWholeCircle ? i1 : i1+nds;     << 314   int ii1 = ifWholeCircle ? i1 : i1+ns;
448   G4int ii2 = ifWholeCircle ? i2 : i2+nds;     << 315   int ii2 = ifWholeCircle ? i2 : i2+ns;
449   G4int vv  = ifWholeCircle ? vEdge : 1;       << 316   int vv  = ifWholeCircle ? vEdge : 1;
450                                                   317 
451   if (nds == 1) {                              << 318   if (ns == 1) {
452     if (r1 == 0.) {                               319     if (r1 == 0.) {
453       pF[kface++]   = G4Facet(i1,0,    v2*i2,0    320       pF[kface++]   = G4Facet(i1,0,    v2*i2,0, (i2+1),0);
454     }else if (r2 == 0.) {                         321     }else if (r2 == 0.) {
455       pF[kface++]   = G4Facet(i1,0,    i2,0,      322       pF[kface++]   = G4Facet(i1,0,    i2,0,    v1*(i1+1),0);
456     }else{                                        323     }else{
457       pF[kface++]   = G4Facet(i1,0,    v2*i2,0    324       pF[kface++]   = G4Facet(i1,0,    v2*i2,0, (i2+1),0, v1*(i1+1),0);
458     }                                             325     }
459   }else{                                          326   }else{
460     if (r1 == 0.) {                               327     if (r1 == 0.) {
461       pF[kface++]   = G4Facet(vv*i1,0,    v2*i    328       pF[kface++]   = G4Facet(vv*i1,0,    v2*i2,0, vEdge*(i2+1),0);
462       for (i2++,i=1; i<nds-1; i2++,i++) {      << 329       for (i2++,i=1; i<ns-1; i2++,i++) {
463         pF[kface++] = G4Facet(vEdge*i1,0, v2*i    330         pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
464       }                                           331       }
465       pF[kface++]   = G4Facet(vEdge*i1,0, v2*i    332       pF[kface++]   = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
466     }else if (r2 == 0.) {                         333     }else if (r2 == 0.) {
467       pF[kface++]   = G4Facet(vv*i1,0,    vEdg    334       pF[kface++]   = G4Facet(vv*i1,0,    vEdge*i2,0, v1*(i1+1),0);
468       for (i1++,i=1; i<nds-1; i1++,i++) {      << 335       for (i1++,i=1; i<ns-1; i1++,i++) {
469         pF[kface++] = G4Facet(vEdge*i1,0, vEdg    336         pF[kface++] = G4Facet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
470       }                                           337       }
471       pF[kface++]   = G4Facet(vEdge*i1,0, vv*i    338       pF[kface++]   = G4Facet(vEdge*i1,0, vv*i2,0,    v1*ii1,0);
472     }else{                                        339     }else{
473       pF[kface++]   = G4Facet(vv*i1,0,    v2*i    340       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 << 341       for (i1++,i2++,i=1; i<ns-1; i1++,i2++,i++) {
475         pF[kface++] = G4Facet(vEdge*i1,0, v2*i    342         pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
476       }                                        << 343       }  
477       pF[kface++]   = G4Facet(vEdge*i1,0, v2*i    344       pF[kface++]   = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0,      v1*ii1,0);
478     }                                             345     }
479   }                                               346   }
480 }                                                 347 }
481                                                   348 
482 void HepPolyhedron::SetSideFacets(G4int ii[4], << 349 void HepPolyhedron::SetSideFacets(int ii[4], int vv[4],
483                                  G4int *kk, G4 << 350                                  int *kk, double *r,
484                                  G4double dphi << 351                                  double dphi, int ns, int &kface)
485 /*********************************************    352 /***********************************************************************
486  *                                                353  *                                                                     *
487  * Name: HepPolyhedron::SetSideFacets             354  * Name: HepPolyhedron::SetSideFacets                Date:    20.05.97 *
488  * Author: E.Chernyaev (IHEP/Protvino)            355  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
489  *                                                356  *                                                                     *
490  * Function: Set side facets for the case of i    357  * Function: Set side facets for the case of incomplete rotation       *
491  *                                                358  *                                                                     *
492  * Input: ii[4] - indices of original vertices << 359  * Input: ii[4] - indeces of original verteces                         *
493  *        vv[4] - visibility of edges             360  *        vv[4] - visibility of edges                                  *
494  *        kk[]  - indices of nodes             << 361  *        kk[]  - indeces of nodes                                     *
495  *        r[]   - radiuses                        362  *        r[]   - radiuses                                             *
496  *        dphi  - delta phi                       363  *        dphi  - delta phi                                            *
497  *        nds    - number of discrete steps    << 364  *        ns     - number of discrete steps                            *
498  *        kface  - current free cell in the pF    365  *        kface  - current free cell in the pF array                   *
499  *                                                366  *                                                                     *
500  *********************************************    367  ***********************************************************************/
501 {                                                 368 {
502   G4int k1, k2, k3, k4;                        << 369   int k1, k2, k3, k4;
503                                                << 370   
504   if (std::abs(dphi-pi) < perMillion) { // hal << 371   if (std::abs((double)(dphi-pi)) < perMillion) {          // half a circle
505     for (G4int i=0; i<4; i++) {                << 372     for (int i=0; i<4; i++) {
506       k1 = ii[i];                                 373       k1 = ii[i];
507       k2 = ii[(i+1)%4];                        << 374       k2 = (i == 3) ? ii[0] : ii[i+1];
508       if (r[k1] == 0. && r[k2] == 0.) vv[i] =  << 375       if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;      
509     }                                             376     }
510   }                                               377   }
511                                                   378 
512   if (ii[1] == ii[2]) {                           379   if (ii[1] == ii[2]) {
513     k1 = kk[ii[0]];                               380     k1 = kk[ii[0]];
514     k2 = kk[ii[2]];                               381     k2 = kk[ii[2]];
515     k3 = kk[ii[3]];                               382     k3 = kk[ii[3]];
516     pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2    383     pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
517     if (r[ii[0]] != 0.) k1 += nds;             << 384     if (r[ii[0]] != 0.) k1 += ns;
518     if (r[ii[2]] != 0.) k2 += nds;             << 385     if (r[ii[2]] != 0.) k2 += ns;
519     if (r[ii[3]] != 0.) k3 += nds;             << 386     if (r[ii[3]] != 0.) k3 += ns;
520     pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2    387     pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
521   }else if (kk[ii[0]] == kk[ii[1]]) {             388   }else if (kk[ii[0]] == kk[ii[1]]) {
522     k1 = kk[ii[0]];                               389     k1 = kk[ii[0]];
523     k2 = kk[ii[2]];                               390     k2 = kk[ii[2]];
524     k3 = kk[ii[3]];                               391     k3 = kk[ii[3]];
525     pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2    392     pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
526     if (r[ii[0]] != 0.) k1 += nds;             << 393     if (r[ii[0]] != 0.) k1 += ns;
527     if (r[ii[2]] != 0.) k2 += nds;             << 394     if (r[ii[2]] != 0.) k2 += ns;
528     if (r[ii[3]] != 0.) k3 += nds;             << 395     if (r[ii[3]] != 0.) k3 += ns;
529     pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2    396     pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
530   }else if (kk[ii[2]] == kk[ii[3]]) {             397   }else if (kk[ii[2]] == kk[ii[3]]) {
531     k1 = kk[ii[0]];                               398     k1 = kk[ii[0]];
532     k2 = kk[ii[1]];                               399     k2 = kk[ii[1]];
533     k3 = kk[ii[2]];                               400     k3 = kk[ii[2]];
534     pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2    401     pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
535     if (r[ii[0]] != 0.) k1 += nds;             << 402     if (r[ii[0]] != 0.) k1 += ns;
536     if (r[ii[1]] != 0.) k2 += nds;             << 403     if (r[ii[1]] != 0.) k2 += ns;
537     if (r[ii[2]] != 0.) k3 += nds;             << 404     if (r[ii[2]] != 0.) k3 += ns;
538     pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2    405     pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
539   }else{                                          406   }else{
540     k1 = kk[ii[0]];                               407     k1 = kk[ii[0]];
541     k2 = kk[ii[1]];                               408     k2 = kk[ii[1]];
542     k3 = kk[ii[2]];                               409     k3 = kk[ii[2]];
543     k4 = kk[ii[3]];                               410     k4 = kk[ii[3]];
544     pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2    411     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;             << 412     if (r[ii[0]] != 0.) k1 += ns;
546     if (r[ii[1]] != 0.) k2 += nds;             << 413     if (r[ii[1]] != 0.) k2 += ns;
547     if (r[ii[2]] != 0.) k3 += nds;             << 414     if (r[ii[2]] != 0.) k3 += ns;
548     if (r[ii[3]] != 0.) k4 += nds;             << 415     if (r[ii[3]] != 0.) k4 += ns;
549     pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3    416     pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
550   }                                               417   }
551 }                                                 418 }
552                                                   419 
553 void HepPolyhedron::RotateAroundZ(G4int nstep, << 420 void HepPolyhedron::RotateAroundZ(int nstep, double phi, double dphi,
554                                  G4int np1, G4 << 421                                  int np1, int np2,
555                                  const G4doubl << 422                                  const double *z, double *r,
556                                  G4int nodeVis << 423                                  int nodeVis, int edgeVis)
557 /*********************************************    424 /***********************************************************************
558  *                                                425  *                                                                     *
559  * Name: HepPolyhedron::RotateAroundZ             426  * Name: HepPolyhedron::RotateAroundZ                Date:    27.11.96 *
560  * Author: E.Chernyaev (IHEP/Protvino)            427  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
561  *                                                428  *                                                                     *
562  * Function: Create HepPolyhedron for a solid     429  * Function: Create HepPolyhedron for a solid produced by rotation of  *
563  *           two polylines around Z-axis          430  *           two polylines around Z-axis                               *
564  *                                                431  *                                                                     *
565  * Input: nstep - number of discrete steps, if    432  * Input: nstep - number of discrete steps, if 0 then default          *
566  *        phi   - starting phi angle              433  *        phi   - starting phi angle                                   *
567  *        dphi  - delta phi                       434  *        dphi  - delta phi                                            *
568  *        np1   - number of points in external    435  *        np1   - number of points in external polyline                *
569  *                (must be negative in case of    436  *                (must be negative in case of closed polyline)        *
570  *        np2   - number of points in internal    437  *        np2   - number of points in internal polyline (may be 1)     *
571  *        z[]   - z-coordinates (+z >>> -z for    438  *        z[]   - z-coordinates (+z >>> -z for both polylines)         *
572  *        r[]   - r-coordinates                   439  *        r[]   - r-coordinates                                        *
573  *        nodeVis - how to Draw edges joing co    440  *        nodeVis - how to Draw edges joing consecutive positions of   *
574  *                  node during rotation          441  *                  node during rotation                               *
575  *        edgeVis - how to Draw edges             442  *        edgeVis - how to Draw edges                                  *
576  *                                                443  *                                                                     *
577  *********************************************    444  ***********************************************************************/
578 {                                                 445 {
579   static const G4double wholeCircle   = twopi; << 446   static double wholeCircle   = twopi;
580                                                << 447     
581   //   S E T   R O T A T I O N   P A R A M E T    448   //   S E T   R O T A T I O N   P A R A M E T E R S
582                                                   449 
583   G4bool ifWholeCircle = std::abs(dphi-wholeCi << 450   bool ifWholeCircle = (std::abs(dphi-wholeCircle) < perMillion) ? true : false;
584   G4double delPhi = ifWholeCircle ? wholeCircl << 451   double   delPhi  = ifWholeCircle ? wholeCircle : dphi;  
585   G4int nSphi = nstep;                         << 452   int        nSphi    = (nstep > 0) ?
586   if (nSphi <= 0) nSphi = GetNumberOfRotationS << 453     nstep : int(delPhi*GetNumberOfRotationSteps()/wholeCircle+.5);
587   if (nSphi == 0) nSphi = 1;                      454   if (nSphi == 0) nSphi = 1;
588   G4int nVphi = ifWholeCircle ? nSphi : nSphi  << 455   int        nVphi    = ifWholeCircle ? nSphi : nSphi+1;
589   G4bool ifClosed = np1 <= 0; // true if exter << 456   bool ifClosed = np1 > 0 ? false : true;
590                                                << 457   
591   //   C O U N T   V E R T I C E S             << 458   //   C O U N T   V E R T E C E S
592                                                   459 
593   G4int absNp1 = std::abs(np1);                << 460   int absNp1 = std::abs(np1);
594   G4int absNp2 = std::abs(np2);                << 461   int absNp2 = std::abs(np2);
595   G4int i1beg = 0;                             << 462   int i1beg = 0;
596   G4int i1end = absNp1-1;                      << 463   int i1end = absNp1-1;
597   G4int i2beg = absNp1;                        << 464   int i2beg = absNp1;
598   G4int i2end = absNp1+absNp2-1;               << 465   int i2end = absNp1+absNp2-1; 
599   G4int i, j, k;                               << 466   int i, j, k;
600                                                   467 
601   for(i=i1beg; i<=i2end; i++) {                   468   for(i=i1beg; i<=i2end; i++) {
602     if (std::abs(r[i]) < spatialTolerance) r[i << 469     if (std::abs(r[i]) < perMillion) r[i] = 0.;
603   }                                               470   }
604                                                   471 
605   // external polyline - check position of nod << 472   j = 0;                                                // external nodes
606   //                                           << 
607   G4int Nverts = 0;                            << 
608   for (i=i1beg; i<=i1end; i++) {                  473   for (i=i1beg; i<=i1end; i++) {
609     Nverts += (r[i] == 0.) ? 1 : nVphi;        << 474     j += (r[i] == 0.) ? 1 : nVphi;
610   }                                               475   }
611                                                   476 
612   // internal polyline                         << 477   bool ifSide1 = false;                           // internal nodes
613   //                                           << 478   bool ifSide2 = false;
614   G4bool ifSide1 = false; // whether to create << 
615   G4bool ifSide2 = false; // whether to create << 
616                                                   479 
617   if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1 << 480   if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
618     Nverts += (r[i2beg] == 0.) ? 1 : nVphi;    << 481     j += (r[i2beg] == 0.) ? 1 : nVphi;
619     ifSide1 = true;                               482     ifSide1 = true;
620   }                                               483   }
621                                                   484 
622   for(i=i2beg+1; i<i2end; i++) { // intermedia << 485   for(i=i2beg+1; i<i2end; i++) {
623     Nverts += (r[i] == 0.) ? 1 : nVphi;        << 486     j += (r[i] == 0.) ? 1 : nVphi;
624   }                                               487   }
625                                                << 488   
626   if (r[i2end] != r[i1end] || z[i2end] != z[i1 << 489   if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
627     if (absNp2 > 1) Nverts += (r[i2end] == 0.) << 490     if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
628     ifSide2 = true;                               491     ifSide2 = true;
629   }                                               492   }
630                                                   493 
631   //   C O U N T   F A C E S                      494   //   C O U N T   F A C E S
632                                                   495 
633   // external lateral faces                    << 496   k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi;       // external faces
634   //                                           << 
635   G4int Nfaces = ifClosed ? absNp1*nSphi : (ab << 
636                                                   497 
637   // internal lateral faces                    << 498   if (absNp2 > 1) {                                     // internal faces
638   //                                           << 
639   if (absNp2 > 1) {                            << 
640     for(i=i2beg; i<i2end; i++) {                  499     for(i=i2beg; i<i2end; i++) {
641       if (r[i] > 0. || r[i+1] > 0.) Nfaces +=  << 500       if (r[i] > 0. || r[i+1] > 0.)       k += nSphi;
642     }                                             501     }
643                                                   502 
644     if (ifClosed) {                               503     if (ifClosed) {
645       if (r[i2end] > 0. || r[i2beg] > 0.) Nfac << 504       if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
646     }                                             505     }
647   }                                               506   }
648                                                   507 
649   // bottom and top faces                      << 508   if (!ifClosed) {                                      // side faces
650   //                                           << 509     if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
651   if (!ifClosed) {                             << 510     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   }                                               511   }
655                                                   512 
656   // phi_wedge faces                           << 513   if (!ifWholeCircle) {                                 // phi_side faces
657   //                                           << 514     k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
658   if (!ifWholeCircle) {                        << 
659     Nfaces += ifClosed ? 2*absNp1 : 2*(absNp1- << 
660   }                                               515   }
661                                                   516 
662   //   A L L O C A T E   M E M O R Y              517   //   A L L O C A T E   M E M O R Y
663                                                   518 
664   AllocateMemory(Nverts, Nfaces);              << 519   AllocateMemory(j, k);
665   if (pV == nullptr || pF == nullptr) return;  << 
666                                                   520 
667   //   G E N E R A T E   V E R T I C E S       << 521   //   G E N E R A T E   V E R T E C E S
668                                                   522 
669   G4int *kk; // array of start indices along p << 523   int *kk;
670   kk = new G4int[absNp1+absNp2];               << 524   kk = new int[absNp1+absNp2];
671                                                   525 
672   // external polyline                         << 526   k = 1;
673   //                                           << 
674   k = 1; // free position in array of vertices << 
675   for(i=i1beg; i<=i1end; i++) {                   527   for(i=i1beg; i<=i1end; i++) {
676     kk[i] = k;                                    528     kk[i] = k;
677     if (r[i] == 0.)                               529     if (r[i] == 0.)
678     { pV[k++] = G4Point3D(0, 0, z[i]); } else  << 530     { pV[k++] = Point3D<double>(0, 0, z[i]); } else { k += nVphi; }
679   }                                               531   }
680                                                   532 
681   // first point of internal polyline          << 
682   //                                           << 
683   i = i2beg;                                      533   i = i2beg;
684   if (ifSide1) {                                  534   if (ifSide1) {
685     kk[i] = k;                                    535     kk[i] = k;
686     if (r[i] == 0.)                               536     if (r[i] == 0.)
687     { pV[k++] = G4Point3D(0, 0, z[i]); } else  << 537     { pV[k++] = Point3D<double>(0, 0, z[i]); } else { k += nVphi; }
688   }else{                                          538   }else{
689     kk[i] = kk[i1beg];                            539     kk[i] = kk[i1beg];
690   }                                               540   }
691                                                   541 
692   // intermediate points of internal polyline  << 
693   //                                           << 
694   for(i=i2beg+1; i<i2end; i++) {                  542   for(i=i2beg+1; i<i2end; i++) {
695     kk[i] = k;                                    543     kk[i] = k;
696     if (r[i] == 0.)                               544     if (r[i] == 0.)
697     { pV[k++] = G4Point3D(0, 0, z[i]); } else  << 545     { pV[k++] = Point3D<double>(0, 0, z[i]); } else { k += nVphi; }
698   }                                               546   }
699                                                   547 
700   // last point of internal polyline           << 
701   //                                           << 
702   if (absNp2 > 1) {                               548   if (absNp2 > 1) {
703     i = i2end;                                    549     i = i2end;
704     if (ifSide2) {                                550     if (ifSide2) {
705       kk[i] = k;                                  551       kk[i] = k;
706       if (r[i] == 0.) pV[k] = G4Point3D(0, 0,  << 552       if (r[i] == 0.) pV[k] = Point3D<double>(0, 0, z[i]);
707     }else{                                        553     }else{
708       kk[i] = kk[i1end];                          554       kk[i] = kk[i1end];
709     }                                             555     }
710   }                                               556   }
711                                                   557 
712   // set vertices                              << 558   double cosPhi, sinPhi;
713   //                                           << 
714   G4double cosPhi, sinPhi;                     << 
715                                                   559 
716   for(j=0; j<nVphi; j++) {                        560   for(j=0; j<nVphi; j++) {
717     cosPhi = std::cos(phi+j*delPhi/nSphi);        561     cosPhi = std::cos(phi+j*delPhi/nSphi);
718     sinPhi = std::sin(phi+j*delPhi/nSphi);        562     sinPhi = std::sin(phi+j*delPhi/nSphi);
719     for(i=i1beg; i<=i2end; i++) {                 563     for(i=i1beg; i<=i2end; i++) {
720       if (r[i] != 0.)                             564       if (r[i] != 0.)
721         pV[kk[i]+j] = G4Point3D(r[i]*cosPhi,r[ << 565         pV[kk[i]+j] = Point3D<double>(r[i]*cosPhi,r[i]*sinPhi,z[i]);
722     }                                             566     }
723   }                                               567   }
724                                                   568 
725   //   G E N E R A T E   F A C E S             << 569   //   G E N E R A T E   E X T E R N A L   F A C E S
726                                                   570 
727   //  external faces                           << 571   int v1,v2;
728   //                                           << 
729   G4int v1,v2;                                 << 
730                                                   572 
731   k = 1; // free position in array of faces pF << 573   k = 1;
732   v2 = ifClosed ? nodeVis : 1;                    574   v2 = ifClosed ? nodeVis : 1;
733   for(i=i1beg; i<i1end; i++) {                    575   for(i=i1beg; i<i1end; i++) {
734     v1 = v2;                                      576     v1 = v2;
735     if (!ifClosed && i == i1end-1) {              577     if (!ifClosed && i == i1end-1) {
736       v2 = 1;                                     578       v2 = 1;
737     }else{                                        579     }else{
738       v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]    580       v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
739     }                                             581     }
740     RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v    582     RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
741                edgeVis, ifWholeCircle, nSphi,     583                edgeVis, ifWholeCircle, nSphi, k);
742   }                                               584   }
743   if (ifClosed) {                                 585   if (ifClosed) {
744     RotateEdge(kk[i1end], kk[i1beg], r[i1end],    586     RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
745                edgeVis, ifWholeCircle, nSphi,     587                edgeVis, ifWholeCircle, nSphi, k);
746   }                                               588   }
747                                                   589 
748   // internal faces                            << 590   //   G E N E R A T E   I N T E R N A L   F A C E S
749   //                                           << 591 
750   if (absNp2 > 1) {                               592   if (absNp2 > 1) {
751     v2 = ifClosed ? nodeVis : 1;                  593     v2 = ifClosed ? nodeVis : 1;
752     for(i=i2beg; i<i2end; i++) {                  594     for(i=i2beg; i<i2end; i++) {
753       v1 = v2;                                    595       v1 = v2;
754       if (!ifClosed && i==i2end-1) {              596       if (!ifClosed && i==i2end-1) {
755         v2 = 1;                                   597         v2 = 1;
756       }else{                                      598       }else{
757         v2 = (r[i] == r[i+1] && r[i+1] == r[i+    599         v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 :  nodeVis;
758       }                                           600       }
759       RotateEdge(kk[i+1], kk[i], r[i+1], r[i],    601       RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
760                  edgeVis, ifWholeCircle, nSphi    602                  edgeVis, ifWholeCircle, nSphi, k);
761     }                                             603     }
762     if (ifClosed) {                               604     if (ifClosed) {
763       RotateEdge(kk[i2beg], kk[i2end], r[i2beg    605       RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
764                  edgeVis, ifWholeCircle, nSphi    606                  edgeVis, ifWholeCircle, nSphi, k);
765     }                                             607     }
766   }                                               608   }
767                                                   609 
768   // bottom and top faces                      << 610   //   G E N E R A T E   S I D E   F A C E S
769   //                                           << 611 
770   if (!ifClosed) {                                612   if (!ifClosed) {
771     if (ifSide1) {                                613     if (ifSide1) {
772       RotateEdge(kk[i2beg], kk[i1beg], r[i2beg    614       RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
773                  -1, ifWholeCircle, nSphi, k);    615                  -1, ifWholeCircle, nSphi, k);
774     }                                             616     }
775     if (ifSide2) {                                617     if (ifSide2) {
776       RotateEdge(kk[i1end], kk[i2end], r[i1end    618       RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
777                  -1, ifWholeCircle, nSphi, k);    619                  -1, ifWholeCircle, nSphi, k);
778     }                                             620     }
779   }                                               621   }
780                                                   622 
781   // phi_wedge faces in case of incomplete cir << 623   //   G E N E R A T E   S I D E   F A C E S  for the case of incomplete circle
782   //                                           << 624 
783   if (!ifWholeCircle) {                           625   if (!ifWholeCircle) {
784                                                   626 
785     G4int  ii[4], vv[4];                       << 627     int  ii[4], vv[4];
786                                                   628 
787     if (ifClosed) {                               629     if (ifClosed) {
788       for (i=i1beg; i<=i1end; i++) {              630       for (i=i1beg; i<=i1end; i++) {
789         ii[0] = i;                                631         ii[0] = i;
790         ii[3] = (i == i1end) ? i1beg : i+1;       632         ii[3] = (i == i1end) ? i1beg : i+1;
791         ii[1] = (absNp2 == 1) ? i2beg : ii[0]+    633         ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
792         ii[2] = (absNp2 == 1) ? i2beg : ii[3]+    634         ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
793         vv[0] = -1;                               635         vv[0] = -1;
794         vv[1] = 1;                                636         vv[1] = 1;
795         vv[2] = -1;                               637         vv[2] = -1;
796         vv[3] = 1;                                638         vv[3] = 1;
797         SetSideFacets(ii, vv, kk, r, delPhi, n << 639         SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
798       }                                           640       }
799     }else{                                        641     }else{
800       for (i=i1beg; i<i1end; i++) {               642       for (i=i1beg; i<i1end; i++) {
801         ii[0] = i;                                643         ii[0] = i;
802         ii[3] = i+1;                              644         ii[3] = i+1;
803         ii[1] = (absNp2 == 1) ? i2beg : ii[0]+    645         ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
804         ii[2] = (absNp2 == 1) ? i2beg : ii[3]+    646         ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
805         vv[0] = (i == i1beg)   ? 1 : -1;          647         vv[0] = (i == i1beg)   ? 1 : -1;
806         vv[1] = 1;                                648         vv[1] = 1;
807         vv[2] = (i == i1end-1) ? 1 : -1;          649         vv[2] = (i == i1end-1) ? 1 : -1;
808         vv[3] = 1;                                650         vv[3] = 1;
809         SetSideFacets(ii, vv, kk, r, delPhi, n << 651         SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
810       }                                           652       }
811     }                                          << 653     }      
812   }                                               654   }
813                                                   655 
814   delete [] kk; // free memory                 << 656   delete [] kk;
815                                                   657 
816   // final check                               << 
817   //                                           << 
818   if (k-1 != nface) {                             658   if (k-1 != nface) {
819     std::cerr                                     659     std::cerr
820       << "HepPolyhedron::RotateAroundZ: number << 660       << "Polyhedron::RotateAroundZ: number of generated faces ("
821       << k-1 << ") is not equal to the number     661       << k-1 << ") is not equal to the number of allocated faces ("
822       << nface << ")"                             662       << nface << ")"
823       << std::endl;                               663       << std::endl;
824   }                                               664   }
825 }                                                 665 }
826                                                   666 
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   }                                           << 
1112   return true;                                << 
1113 }                                             << 
1114                                               << 
1115 void HepPolyhedron::SetReferences()              667 void HepPolyhedron::SetReferences()
1116 /********************************************    668 /***********************************************************************
1117  *                                               669  *                                                                     *
1118  * Name: HepPolyhedron::SetReferences            670  * Name: HepPolyhedron::SetReferences                Date:    04.12.96 *
1119  * Author: E.Chernyaev (IHEP/Protvino)           671  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
1120  *                                               672  *                                                                     *
1121  * Function: For each edge set reference to n    673  * Function: For each edge set reference to neighbouring facet         *
1122  *                                               674  *                                                                     *
1123  ********************************************    675  ***********************************************************************/
1124 {                                                676 {
1125   if (nface <= 0) return;                        677   if (nface <= 0) return;
1126                                                  678 
1127   struct edgeListMember {                        679   struct edgeListMember {
1128     edgeListMember *next;                        680     edgeListMember *next;
1129     G4int v2;                                 << 681     int v2;
1130     G4int iface;                              << 682     int iface;
1131     G4int iedge;                              << 683     int iedge;
1132   } *edgeList, *freeList, **headList;            684   } *edgeList, *freeList, **headList;
1133                                                  685 
1134                                               << 686   
1135   //   A L L O C A T E   A N D   I N I T I A     687   //   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                                                  688 
1137   edgeList = new edgeListMember[2*nface];        689   edgeList = new edgeListMember[2*nface];
1138   headList = new edgeListMember*[nvert];         690   headList = new edgeListMember*[nvert];
1139                                               << 691   
1140   G4int i;                                    << 692   int i;
1141   for (i=0; i<nvert; i++) {                      693   for (i=0; i<nvert; i++) {
1142     headList[i] = nullptr;                    << 694     headList[i] = 0;
1143   }                                              695   }
1144   freeList = edgeList;                           696   freeList = edgeList;
1145   for (i=0; i<2*nface-1; i++) {                  697   for (i=0; i<2*nface-1; i++) {
1146     edgeList[i].next = &edgeList[i+1];           698     edgeList[i].next = &edgeList[i+1];
1147   }                                              699   }
1148   edgeList[2*nface-1].next = nullptr;         << 700   edgeList[2*nface-1].next = 0;
1149                                                  701 
1150   //   L O O P   A L O N G   E D G E S           702   //   L O O P   A L O N G   E D G E S
1151                                                  703 
1152   G4int iface, iedge, nedge, i1, i2, k1, k2;  << 704   int iface, iedge, nedge, i1, i2, k1, k2;
1153   edgeListMember *prev, *cur;                    705   edgeListMember *prev, *cur;
1154                                               << 706   
1155   for(iface=1; iface<=nface; iface++) {          707   for(iface=1; iface<=nface; iface++) {
1156     nedge = (pF[iface].edge[3].v == 0) ? 3 :     708     nedge = (pF[iface].edge[3].v == 0) ? 3 : 4;
1157     for (iedge=0; iedge<nedge; iedge++) {        709     for (iedge=0; iedge<nedge; iedge++) {
1158       i1 = iedge;                                710       i1 = iedge;
1159       i2 = (iedge < nedge-1) ? iedge+1 : 0;      711       i2 = (iedge < nedge-1) ? iedge+1 : 0;
1160       i1 = std::abs(pF[iface].edge[i1].v);       712       i1 = std::abs(pF[iface].edge[i1].v);
1161       i2 = std::abs(pF[iface].edge[i2].v);       713       i2 = std::abs(pF[iface].edge[i2].v);
1162       k1 = (i1 < i2) ? i1 : i2;          // k    714       k1 = (i1 < i2) ? i1 : i2;          // k1 = ::min(i1,i2);
1163       k2 = (i1 > i2) ? i1 : i2;          // k    715       k2 = (i1 > i2) ? i1 : i2;          // k2 = ::max(i1,i2);
1164                                               << 716       
1165       // check head of the List corresponding    717       // check head of the List corresponding to k1
1166       cur = headList[k1];                        718       cur = headList[k1];
1167       if (cur == nullptr) {                   << 719       if (cur == 0) {
1168         headList[k1] = freeList;                 720         headList[k1] = freeList;
1169         if (freeList == nullptr) {            << 
1170           std::cerr                           << 
1171           << "Polyhedron::SetReferences: bad  << 
1172           << std::endl;                       << 
1173           break;                              << 
1174         }                                     << 
1175         freeList = freeList->next;               721         freeList = freeList->next;
1176         cur = headList[k1];                      722         cur = headList[k1];
1177         cur->next = nullptr;                  << 723         cur->next = 0;
1178         cur->v2 = k2;                            724         cur->v2 = k2;
1179         cur->iface = iface;                      725         cur->iface = iface;
1180         cur->iedge = iedge;                      726         cur->iedge = iedge;
1181         continue;                                727         continue;
1182       }                                          728       }
1183                                                  729 
1184       if (cur->v2 == k2) {                       730       if (cur->v2 == k2) {
1185         headList[k1] = cur->next;                731         headList[k1] = cur->next;
1186         cur->next = freeList;                    732         cur->next = freeList;
1187         freeList = cur;                       << 733         freeList = cur;      
1188         pF[iface].edge[iedge].f = cur->iface;    734         pF[iface].edge[iedge].f = cur->iface;
1189         pF[cur->iface].edge[cur->iedge].f = i    735         pF[cur->iface].edge[cur->iedge].f = iface;
1190         i1 = (pF[iface].edge[iedge].v < 0) ?     736         i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
1191         i2 = (pF[cur->iface].edge[cur->iedge]    737         i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
1192         if (i1 != i2) {                          738         if (i1 != i2) {
1193           std::cerr                              739           std::cerr
1194             << "Polyhedron::SetReferences: di    740             << "Polyhedron::SetReferences: different edge visibility "
1195             << iface << "/" << iedge << "/"      741             << iface << "/" << iedge << "/"
1196             << pF[iface].edge[iedge].v << " a    742             << pF[iface].edge[iedge].v << " and "
1197             << cur->iface << "/" << cur->iedg    743             << cur->iface << "/" << cur->iedge << "/"
1198             << pF[cur->iface].edge[cur->iedge    744             << pF[cur->iface].edge[cur->iedge].v
1199             << std::endl;                        745             << std::endl;
1200         }                                        746         }
1201         continue;                                747         continue;
1202       }                                          748       }
1203                                                  749 
1204       // check List itself                       750       // check List itself
1205       for (;;) {                                 751       for (;;) {
1206         prev = cur;                              752         prev = cur;
1207         cur = prev->next;                        753         cur = prev->next;
1208         if (cur == nullptr) {                 << 754         if (cur == 0) {
1209           prev->next = freeList;                 755           prev->next = freeList;
1210           if (freeList == nullptr) {          << 
1211             std::cerr                         << 
1212             << "Polyhedron::SetReferences: ba << 
1213             << std::endl;                     << 
1214             break;                            << 
1215           }                                   << 
1216           freeList = freeList->next;             756           freeList = freeList->next;
1217           cur = prev->next;                      757           cur = prev->next;
1218           cur->next = nullptr;                << 758           cur->next = 0;
1219           cur->v2 = k2;                          759           cur->v2 = k2;
1220           cur->iface = iface;                    760           cur->iface = iface;
1221           cur->iedge = iedge;                    761           cur->iedge = iedge;
1222           break;                                 762           break;
1223         }                                        763         }
1224                                                  764 
1225         if (cur->v2 == k2) {                     765         if (cur->v2 == k2) {
1226           prev->next = cur->next;                766           prev->next = cur->next;
1227           cur->next = freeList;                  767           cur->next = freeList;
1228           freeList = cur;                     << 768           freeList = cur;      
1229           pF[iface].edge[iedge].f = cur->ifac    769           pF[iface].edge[iedge].f = cur->iface;
1230           pF[cur->iface].edge[cur->iedge].f =    770           pF[cur->iface].edge[cur->iedge].f = iface;
1231           i1 = (pF[iface].edge[iedge].v < 0)     771           i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
1232           i2 = (pF[cur->iface].edge[cur->iedg    772           i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
1233             if (i1 != i2) {                      773             if (i1 != i2) {
1234               std::cerr                          774               std::cerr
1235                 << "Polyhedron::SetReferences    775                 << "Polyhedron::SetReferences: different edge visibility "
1236                 << iface << "/" << iedge << "    776                 << iface << "/" << iedge << "/"
1237                 << pF[iface].edge[iedge].v <<    777                 << pF[iface].edge[iedge].v << " and "
1238                 << cur->iface << "/" << cur->    778                 << cur->iface << "/" << cur->iedge << "/"
1239                 << pF[cur->iface].edge[cur->i    779                 << pF[cur->iface].edge[cur->iedge].v
1240                 << std::endl;                    780                 << std::endl;
1241             }                                    781             }
1242           break;                                 782           break;
1243         }                                        783         }
1244       }                                          784       }
1245     }                                            785     }
1246   }                                              786   }
1247                                                  787 
1248   //  C H E C K   T H A T   A L L   L I S T S    788   //  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                                                  789 
1250   for (i=0; i<nvert; i++) {                      790   for (i=0; i<nvert; i++) {
1251     if (headList[i] != nullptr) {             << 791     if (headList[i] != 0) {
1252       std::cerr                                  792       std::cerr
1253         << "Polyhedron::SetReferences: List "    793         << "Polyhedron::SetReferences: List " << i << " is not empty"
1254         << std::endl;                            794         << std::endl;
1255     }                                            795     }
1256   }                                              796   }
1257                                                  797 
1258   //   F R E E   M E M O R Y                     798   //   F R E E   M E M O R Y
1259                                                  799 
1260   delete [] edgeList;                            800   delete [] edgeList;
1261   delete [] headList;                            801   delete [] headList;
1262 }                                                802 }
1263                                                  803 
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()               804 void HepPolyhedron::InvertFacets()
1354 /********************************************    805 /***********************************************************************
1355  *                                               806  *                                                                     *
1356  * Name: HepPolyhedron::InvertFacets             807  * Name: HepPolyhedron::InvertFacets                Date:    01.12.99  *
1357  * Author: E.Chernyaev                           808  * Author: E.Chernyaev                              Revised:           *
1358  *                                               809  *                                                                     *
1359  * Function: Invert the order of the nodes in    810  * Function: Invert the order of the nodes in the facets               *
1360  *                                               811  *                                                                     *
1361  ********************************************    812  ***********************************************************************/
1362 {                                                813 {
1363   if (nface <= 0) return;                        814   if (nface <= 0) return;
1364   G4int i, k, nnode, v[4],f[4];               << 815   int i, k, nnode, v[4],f[4];
1365   for (i=1; i<=nface; i++) {                     816   for (i=1; i<=nface; i++) {
1366     nnode =  (pF[i].edge[3].v == 0) ? 3 : 4;     817     nnode =  (pF[i].edge[3].v == 0) ? 3 : 4;
1367     for (k=0; k<nnode; k++) {                    818     for (k=0; k<nnode; k++) {
1368       v[k] = (k+1 == nnode) ? pF[i].edge[0].v    819       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] =     820       if (v[k] * pF[i].edge[k].v < 0) v[k] = -v[k];
1370       f[k] = pF[i].edge[k].f;                    821       f[k] = pF[i].edge[k].f;
1371     }                                            822     }
1372     for (k=0; k<nnode; k++) {                    823     for (k=0; k<nnode; k++) {
1373       pF[i].edge[nnode-1-k].v = v[k];            824       pF[i].edge[nnode-1-k].v = v[k];
1374       pF[i].edge[nnode-1-k].f = f[k];            825       pF[i].edge[nnode-1-k].f = f[k];
1375     }                                            826     }
1376   }                                              827   }
1377 }                                                828 }
1378                                                  829 
1379 HepPolyhedron & HepPolyhedron::Transform(cons << 830 HepPolyhedron & HepPolyhedron::Transform(const Transform3D &t)
1380 /********************************************    831 /***********************************************************************
1381  *                                               832  *                                                                     *
1382  * Name: HepPolyhedron::Transform                833  * Name: HepPolyhedron::Transform                    Date:    01.12.99  *
1383  * Author: E.Chernyaev                           834  * Author: E.Chernyaev                              Revised:           *
1384  *                                               835  *                                                                     *
1385  * Function: Make transformation of the polyh    836  * Function: Make transformation of the polyhedron                     *
1386  *                                               837  *                                                                     *
1387  ********************************************    838  ***********************************************************************/
1388 {                                                839 {
1389   if (nvert > 0) {                               840   if (nvert > 0) {
1390     for (G4int i=1; i<=nvert; i++) { pV[i] =  << 841     for (int i=1; i<=nvert; i++) { pV[i] = t * pV[i]; }
1391                                                  842 
1392     //  C H E C K   D E T E R M I N A N T   A    843     //  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    844     //  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                                                  845 
1395     G4Vector3D d = t * G4Vector3D(0,0,0);     << 846     Vector3D<double> d = t * Vector3D<double>(0,0,0);
1396     G4Vector3D x = t * G4Vector3D(1,0,0) - d; << 847     Vector3D<double> x = t * Vector3D<double>(1,0,0) - d;
1397     G4Vector3D y = t * G4Vector3D(0,1,0) - d; << 848     Vector3D<double> y = t * Vector3D<double>(0,1,0) - d;
1398     G4Vector3D z = t * G4Vector3D(0,0,1) - d; << 849     Vector3D<double> z = t * Vector3D<double>(0,0,1) - d;
1399     if ((x.cross(y))*z < 0) InvertFacets();      850     if ((x.cross(y))*z < 0) InvertFacets();
1400   }                                              851   }
1401   return *this;                                  852   return *this;
1402 }                                                853 }
1403                                                  854 
1404 G4bool HepPolyhedron::GetNextVertexIndex(G4in << 855 bool HepPolyhedron::GetNextVertexIndex(int &index, int &edgeFlag) const
1405 /********************************************    856 /***********************************************************************
1406  *                                               857  *                                                                     *
1407  * Name: HepPolyhedron::GetNextVertexIndex       858  * Name: HepPolyhedron::GetNextVertexIndex          Date:    03.09.96  *
1408  * Author: Yasuhide Sawada                       859  * Author: Yasuhide Sawada                          Revised:           *
1409  *                                               860  *                                                                     *
1410  * Function:                                     861  * Function:                                                           *
1411  *                                               862  *                                                                     *
1412  ********************************************    863  ***********************************************************************/
1413 {                                                864 {
1414   static G4ThreadLocal G4int iFace = 1;       << 865   static int iFace = 1;
1415   static G4ThreadLocal G4int iQVertex = 0;    << 866   static int iQVertex = 0;
1416   G4int vIndex = pF[iFace].edge[iQVertex].v;  << 867   int vIndex = pF[iFace].edge[iQVertex].v;
1417                                                  868 
1418   edgeFlag = (vIndex > 0) ? 1 : 0;               869   edgeFlag = (vIndex > 0) ? 1 : 0;
1419   index = std::abs(vIndex);                      870   index = std::abs(vIndex);
1420                                                  871 
1421   if (iQVertex >= 3 || pF[iFace].edge[iQVerte    872   if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1422     iQVertex = 0;                                873     iQVertex = 0;
1423     if (++iFace > nface) iFace = 1;              874     if (++iFace > nface) iFace = 1;
1424     return false;  // Last Edge                  875     return false;  // Last Edge
                                                   >> 876   }else{
                                                   >> 877     ++iQVertex;
                                                   >> 878     return true;  // not Last Edge
1425   }                                              879   }
1426                                               << 
1427   ++iQVertex;                                 << 
1428   return true;  // not Last Edge              << 
1429 }                                                880 }
1430                                                  881 
1431 G4Point3D HepPolyhedron::GetVertex(G4int inde << 882 Point3D<double> HepPolyhedron::GetVertex(int index) const
1432 /********************************************    883 /***********************************************************************
1433  *                                               884  *                                                                     *
1434  * Name: HepPolyhedron::GetVertex                885  * Name: HepPolyhedron::GetVertex                   Date:    03.09.96  *
1435  * Author: Yasuhide Sawada                       886  * Author: Yasuhide Sawada                          Revised: 17.11.99  *
1436  *                                               887  *                                                                     *
1437  * Function: Get vertex of the index.            888  * Function: Get vertex of the index.                                  *
1438  *                                               889  *                                                                     *
1439  ********************************************    890  ***********************************************************************/
1440 {                                                891 {
1441   if (index <= 0 || index > nvert) {             892   if (index <= 0 || index > nvert) {
1442     std::cerr                                    893     std::cerr
1443       << "HepPolyhedron::GetVertex: irrelevan    894       << "HepPolyhedron::GetVertex: irrelevant index " << index
1444       << std::endl;                              895       << std::endl;
1445     return G4Point3D();                       << 896     return Point3D<double>();
1446   }                                              897   }
1447   return pV[index];                              898   return pV[index];
1448 }                                                899 }
1449                                                  900 
1450 G4bool                                        << 901 bool
1451 HepPolyhedron::GetNextVertex(G4Point3D &verte << 902 HepPolyhedron::GetNextVertex(Point3D<double> &vertex, int &edgeFlag) const
1452 /********************************************    903 /***********************************************************************
1453  *                                               904  *                                                                     *
1454  * Name: HepPolyhedron::GetNextVertex            905  * Name: HepPolyhedron::GetNextVertex               Date:    22.07.96  *
1455  * Author: John Allison                          906  * Author: John Allison                             Revised:           *
1456  *                                               907  *                                                                     *
1457  * Function: Get vertices of the quadrilatera    908  * Function: Get vertices of the quadrilaterals in order for each      *
1458  *           face in face order.  Returns fal    909  *           face in face order.  Returns false when finished each     *
1459  *           face.                               910  *           face.                                                     *
1460  *                                               911  *                                                                     *
1461  ********************************************    912  ***********************************************************************/
1462 {                                                913 {
1463   G4int index;                                << 914   int index;
1464   G4bool rep = GetNextVertexIndex(index, edge << 915   bool rep = GetNextVertexIndex(index, edgeFlag);
1465   vertex = pV[index];                            916   vertex = pV[index];
1466   return rep;                                    917   return rep;
1467 }                                                918 }
1468                                                  919 
1469 G4bool HepPolyhedron::GetNextVertex(G4Point3D << 920 bool HepPolyhedron::GetNextVertex(Point3D<double> &vertex, int &edgeFlag,
1470                                   G4Normal3D  << 921                                        Normal3D<double> &normal) const
1471 /********************************************    922 /***********************************************************************
1472  *                                               923  *                                                                     *
1473  * Name: HepPolyhedron::GetNextVertex            924  * Name: HepPolyhedron::GetNextVertex               Date:    26.11.99  *
1474  * Author: E.Chernyaev                           925  * Author: E.Chernyaev                              Revised:           *
1475  *                                               926  *                                                                     *
1476  * Function: Get vertices with normals of the    927  * Function: Get vertices with normals of the quadrilaterals in order  *
1477  *           for each face in face order.        928  *           for each face in face order.                              *
1478  *           Returns false when finished each    929  *           Returns false when finished each face.                    *
1479  *                                               930  *                                                                     *
1480  ********************************************    931  ***********************************************************************/
1481 {                                                932 {
1482   static G4ThreadLocal G4int iFace = 1;       << 933   static int iFace = 1;
1483   static G4ThreadLocal G4int iNode = 0;       << 934   static int iNode = 0;
1484                                                  935 
1485   if (nface == 0) return false;  // empty pol    936   if (nface == 0) return false;  // empty polyhedron
1486                                                  937 
1487   G4int k = pF[iFace].edge[iNode].v;          << 938   int k = pF[iFace].edge[iNode].v;
1488   if (k > 0) { edgeFlag = 1; } else { edgeFla    939   if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
1489   vertex = pV[k];                                940   vertex = pV[k];
1490   normal = FindNodeNormal(iFace,k);              941   normal = FindNodeNormal(iFace,k);
1491   if (iNode >= 3 || pF[iFace].edge[iNode+1].v    942   if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
1492     iNode = 0;                                   943     iNode = 0;
1493     if (++iFace > nface) iFace = 1;              944     if (++iFace > nface) iFace = 1;
1494     return false;                // last node    945     return false;                // last node
                                                   >> 946   }else{
                                                   >> 947     ++iNode;
                                                   >> 948     return true;                 // not last node
1495   }                                              949   }
1496   ++iNode;                                    << 
1497   return true;                 // not last no << 
1498 }                                                950 }
1499                                                  951 
1500 G4bool HepPolyhedron::GetNextEdgeIndices(G4in << 952 bool HepPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag,
1501                                        G4int  << 953                                             int &iface1, int &iface2) const
1502 /********************************************    954 /***********************************************************************
1503  *                                               955  *                                                                     *
1504  * Name: HepPolyhedron::GetNextEdgeIndices    << 956  * Name: HepPolyhedron::GetNextEdgeIndeces          Date:    30.09.96  *
1505  * Author: E.Chernyaev                           957  * Author: E.Chernyaev                              Revised: 17.11.99  *
1506  *                                               958  *                                                                     *
1507  * Function: Get indices of the next edge tog << 959  * Function: Get indeces of the next edge together with indeces of     *
1508  *           of the faces which share the edg    960  *           of the faces which share the edge.                        *
1509  *           Returns false when the last edge    961  *           Returns false when the last edge.                         *
1510  *                                               962  *                                                                     *
1511  ********************************************    963  ***********************************************************************/
1512 {                                                964 {
1513   static G4ThreadLocal G4int iFace    = 1;    << 965   static int iFace    = 1;
1514   static G4ThreadLocal G4int iQVertex = 0;    << 966   static int iQVertex = 0;
1515   static G4ThreadLocal G4int iOrder   = 1;    << 967   static int iOrder   = 1;
1516   G4int  k1, k2, kflag, kface1, kface2;       << 968   int  k1, k2, kflag, kface1, kface2;
1517                                                  969 
1518   if (iFace == 1 && iQVertex == 0) {             970   if (iFace == 1 && iQVertex == 0) {
1519     k2 = pF[nface].edge[0].v;                    971     k2 = pF[nface].edge[0].v;
1520     k1 = pF[nface].edge[3].v;                    972     k1 = pF[nface].edge[3].v;
1521     if (k1 == 0) k1 = pF[nface].edge[2].v;       973     if (k1 == 0) k1 = pF[nface].edge[2].v;
1522     if (std::abs(k1) > std::abs(k2)) iOrder =    974     if (std::abs(k1) > std::abs(k2)) iOrder = -1;
1523   }                                              975   }
1524                                                  976 
1525   do {                                           977   do {
1526     k1     = pF[iFace].edge[iQVertex].v;         978     k1     = pF[iFace].edge[iQVertex].v;
1527     kflag  = k1;                                 979     kflag  = k1;
1528     k1     = std::abs(k1);                       980     k1     = std::abs(k1);
1529     kface1 = iFace;                           << 981     kface1 = iFace; 
1530     kface2 = pF[iFace].edge[iQVertex].f;         982     kface2 = pF[iFace].edge[iQVertex].f;
1531     if (iQVertex >= 3 || pF[iFace].edge[iQVer    983     if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1532       iQVertex = 0;                              984       iQVertex = 0;
1533       k2 = std::abs(pF[iFace].edge[iQVertex].    985       k2 = std::abs(pF[iFace].edge[iQVertex].v);
1534       iFace++;                                   986       iFace++;
1535     }else{                                       987     }else{
1536       iQVertex++;                                988       iQVertex++;
1537       k2 = std::abs(pF[iFace].edge[iQVertex].    989       k2 = std::abs(pF[iFace].edge[iQVertex].v);
1538     }                                            990     }
1539   } while (iOrder*k1 > iOrder*k2);               991   } while (iOrder*k1 > iOrder*k2);
1540                                                  992 
1541   i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ?     993   i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1542   iface1 = kface1; iface2 = kface2;           << 994   iface1 = kface1; iface2 = kface2; 
1543                                                  995 
1544   if (iFace > nface) {                           996   if (iFace > nface) {
1545     iFace  = 1; iOrder = 1;                      997     iFace  = 1; iOrder = 1;
1546     return false;                                998     return false;
                                                   >> 999   }else{
                                                   >> 1000     return true;
1547   }                                              1001   }
1548                                               << 
1549   return true;                                << 
1550 }                                                1002 }
1551                                                  1003 
1552 G4bool                                        << 1004 bool
1553 HepPolyhedron::GetNextEdgeIndices(G4int &i1,  << 1005 HepPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag) const
1554 /********************************************    1006 /***********************************************************************
1555  *                                               1007  *                                                                     *
1556  * Name: HepPolyhedron::GetNextEdgeIndices    << 1008  * Name: HepPolyhedron::GetNextEdgeIndeces          Date:    17.11.99  *
1557  * Author: E.Chernyaev                           1009  * Author: E.Chernyaev                              Revised:           *
1558  *                                               1010  *                                                                     *
1559  * Function: Get indices of the next edge.    << 1011  * Function: Get indeces of the next edge.                             *
1560  *           Returns false when the last edge    1012  *           Returns false when the last edge.                         *
1561  *                                               1013  *                                                                     *
1562  ********************************************    1014  ***********************************************************************/
1563 {                                                1015 {
1564   G4int kface1, kface2;                       << 1016   int kface1, kface2;
1565   return GetNextEdgeIndices(i1, i2, edgeFlag, << 1017   return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
1566 }                                                1018 }
1567                                                  1019 
1568 G4bool                                        << 1020 bool
1569 HepPolyhedron::GetNextEdge(G4Point3D &p1,     << 1021 HepPolyhedron::GetNextEdge(Point3D<double> &p1,
1570                            G4Point3D &p2,     << 1022                            Point3D<double> &p2,
1571                            G4int &edgeFlag) c << 1023                            int &edgeFlag) const
1572 /********************************************    1024 /***********************************************************************
1573  *                                               1025  *                                                                     *
1574  * Name: HepPolyhedron::GetNextEdge              1026  * Name: HepPolyhedron::GetNextEdge                 Date:    30.09.96  *
1575  * Author: E.Chernyaev                           1027  * Author: E.Chernyaev                              Revised:           *
1576  *                                               1028  *                                                                     *
1577  * Function: Get next edge.                      1029  * Function: Get next edge.                                            *
1578  *           Returns false when the last edge    1030  *           Returns false when the last edge.                         *
1579  *                                               1031  *                                                                     *
1580  ********************************************    1032  ***********************************************************************/
1581 {                                                1033 {
1582   G4int i1,i2;                                << 1034   int i1,i2;
1583   G4bool rep = GetNextEdgeIndices(i1,i2,edgeF << 1035   bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
1584   p1 = pV[i1];                                   1036   p1 = pV[i1];
1585   p2 = pV[i2];                                   1037   p2 = pV[i2];
1586   return rep;                                    1038   return rep;
1587 }                                                1039 }
1588                                                  1040 
1589 G4bool                                        << 1041 bool
1590 HepPolyhedron::GetNextEdge(G4Point3D &p1, G4P << 1042 HepPolyhedron::GetNextEdge(Point3D<double> &p1, Point3D<double> &p2,
1591                           G4int &edgeFlag, G4 << 1043                           int &edgeFlag, int &iface1, int &iface2) const
1592 /********************************************    1044 /***********************************************************************
1593  *                                               1045  *                                                                     *
1594  * Name: HepPolyhedron::GetNextEdge              1046  * Name: HepPolyhedron::GetNextEdge                 Date:    17.11.99  *
1595  * Author: E.Chernyaev                           1047  * Author: E.Chernyaev                              Revised:           *
1596  *                                               1048  *                                                                     *
1597  * Function: Get next edge with indices of th << 1049  * Function: Get next edge with indeces of the faces which share       *
1598  *           the edge.                           1050  *           the edge.                                                 *
1599  *           Returns false when the last edge    1051  *           Returns false when the last edge.                         *
1600  *                                               1052  *                                                                     *
1601  ********************************************    1053  ***********************************************************************/
1602 {                                                1054 {
1603   G4int i1,i2;                                << 1055   int i1,i2;
1604   G4bool rep = GetNextEdgeIndices(i1,i2,edgeF << 1056   bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
1605   p1 = pV[i1];                                   1057   p1 = pV[i1];
1606   p2 = pV[i2];                                   1058   p2 = pV[i2];
1607   return rep;                                    1059   return rep;
1608 }                                                1060 }
1609                                                  1061 
1610 void HepPolyhedron::GetFacet(G4int iFace, G4i << 1062 void HepPolyhedron::GetFacet(int iFace, int &n, int *iNodes,
1611                             G4int *edgeFlags, << 1063                             int *edgeFlags, int *iFaces) const
1612 /********************************************    1064 /***********************************************************************
1613  *                                               1065  *                                                                     *
1614  * Name: HepPolyhedron::GetFacet                 1066  * Name: HepPolyhedron::GetFacet                    Date:    15.12.99  *
1615  * Author: E.Chernyaev                           1067  * Author: E.Chernyaev                              Revised:           *
1616  *                                               1068  *                                                                     *
1617  * Function: Get face by index                   1069  * Function: Get face by index                                         *
1618  *                                               1070  *                                                                     *
1619  ********************************************    1071  ***********************************************************************/
1620 {                                                1072 {
1621   if (iFace < 1 || iFace > nface) {              1073   if (iFace < 1 || iFace > nface) {
1622     std::cerr                                 << 1074     std::cerr 
1623       << "HepPolyhedron::GetFacet: irrelevant    1075       << "HepPolyhedron::GetFacet: irrelevant index " << iFace
1624       << std::endl;                              1076       << std::endl;
1625     n = 0;                                       1077     n = 0;
1626   }else{                                         1078   }else{
1627     G4int i, k;                               << 1079     int i, k;
1628     for (i=0; i<4; i++) {                     << 1080     for (i=0; i<4; i++) { 
1629       k = pF[iFace].edge[i].v;                   1081       k = pF[iFace].edge[i].v;
1630       if (k == 0) break;                         1082       if (k == 0) break;
1631       if (iFaces != nullptr) iFaces[i] = pF[i << 1083       if (iFaces != 0) iFaces[i] = pF[iFace].edge[i].f;
1632       if (k > 0) {                            << 1084       if (k > 0) { 
1633         iNodes[i] = k;                           1085         iNodes[i] = k;
1634         if (edgeFlags != nullptr) edgeFlags[i << 1086         if (edgeFlags != 0) edgeFlags[i] = 1;
1635       }else{                                     1087       }else{
1636         iNodes[i] = -k;                          1088         iNodes[i] = -k;
1637         if (edgeFlags != nullptr) edgeFlags[i << 1089         if (edgeFlags != 0) edgeFlags[i] = -1;
1638       }                                          1090       }
1639     }                                            1091     }
1640     n = i;                                       1092     n = i;
1641   }                                              1093   }
1642 }                                                1094 }
1643                                                  1095 
1644 void HepPolyhedron::GetFacet(G4int index, G4i << 1096 void HepPolyhedron::GetFacet(int index, int &n, Point3D<double> *nodes,
1645                              G4int *edgeFlags << 1097                             int *edgeFlags, Normal3D<double> *normals) const
1646 /********************************************    1098 /***********************************************************************
1647  *                                               1099  *                                                                     *
1648  * Name: HepPolyhedron::GetFacet                 1100  * Name: HepPolyhedron::GetFacet                    Date:    17.11.99  *
1649  * Author: E.Chernyaev                           1101  * Author: E.Chernyaev                              Revised:           *
1650  *                                               1102  *                                                                     *
1651  * Function: Get face by index                   1103  * Function: Get face by index                                         *
1652  *                                               1104  *                                                                     *
1653  ********************************************    1105  ***********************************************************************/
1654 {                                                1106 {
1655   G4int iNodes[4];                            << 1107   int iNodes[4];
1656   GetFacet(index, n, iNodes, edgeFlags);         1108   GetFacet(index, n, iNodes, edgeFlags);
1657   if (n != 0) {                                  1109   if (n != 0) {
1658     for (G4int i=0; i<n; i++) {               << 1110     for (int i=0; i<n; i++) {
1659       nodes[i] = pV[iNodes[i]];                  1111       nodes[i] = pV[iNodes[i]];
1660       if (normals != nullptr) normals[i] = Fi << 1112       if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1661     }                                            1113     }
1662   }                                              1114   }
1663 }                                                1115 }
1664                                                  1116 
1665 G4bool                                        << 1117 bool
1666 HepPolyhedron::GetNextFacet(G4int &n, G4Point << 1118 HepPolyhedron::GetNextFacet(int &n, Point3D<double> *nodes,
1667                            G4int *edgeFlags,  << 1119                            int *edgeFlags, Normal3D<double> *normals) const
1668 /********************************************    1120 /***********************************************************************
1669  *                                               1121  *                                                                     *
1670  * Name: HepPolyhedron::GetNextFacet             1122  * Name: HepPolyhedron::GetNextFacet                Date:    19.11.99  *
1671  * Author: E.Chernyaev                           1123  * Author: E.Chernyaev                              Revised:           *
1672  *                                               1124  *                                                                     *
1673  * Function: Get next face with normals of un    1125  * Function: Get next face with normals of unit length at the nodes.   *
1674  *           Returns false when finished all     1126  *           Returns false when finished all faces.                    *
1675  *                                               1127  *                                                                     *
1676  ********************************************    1128  ***********************************************************************/
1677 {                                                1129 {
1678   static G4ThreadLocal G4int iFace = 1;       << 1130   static int iFace = 1;
1679                                                  1131 
1680   if (edgeFlags == nullptr) {                 << 1132   if (edgeFlags == 0) {
1681     GetFacet(iFace, n, nodes);                   1133     GetFacet(iFace, n, nodes);
1682   }else if (normals == nullptr) {             << 1134   }else if (normals == 0) {
1683     GetFacet(iFace, n, nodes, edgeFlags);        1135     GetFacet(iFace, n, nodes, edgeFlags);
1684   }else{                                         1136   }else{
1685     GetFacet(iFace, n, nodes, edgeFlags, norm    1137     GetFacet(iFace, n, nodes, edgeFlags, normals);
1686   }                                              1138   }
1687                                                  1139 
1688   if (++iFace > nface) {                         1140   if (++iFace > nface) {
1689     iFace  = 1;                                  1141     iFace  = 1;
1690     return false;                                1142     return false;
                                                   >> 1143   }else{
                                                   >> 1144     return true;
1691   }                                              1145   }
1692                                               << 
1693   return true;                                << 
1694 }                                                1146 }
1695                                                  1147 
1696 G4Normal3D HepPolyhedron::GetNormal(G4int iFa << 1148 Normal3D<double> HepPolyhedron::GetNormal(int iFace) const
1697 /********************************************    1149 /***********************************************************************
1698  *                                               1150  *                                                                     *
1699  * Name: HepPolyhedron::GetNormal                1151  * Name: HepPolyhedron::GetNormal                    Date:    19.11.99 *
1700  * Author: E.Chernyaev                           1152  * Author: E.Chernyaev                               Revised:          *
1701  *                                               1153  *                                                                     *
1702  * Function: Get normal of the face given by     1154  * Function: Get normal of the face given by index                     *
1703  *                                               1155  *                                                                     *
1704  ********************************************    1156  ***********************************************************************/
1705 {                                                1157 {
1706   if (iFace < 1 || iFace > nface) {              1158   if (iFace < 1 || iFace > nface) {
1707     std::cerr                                 << 1159     std::cerr 
1708       << "HepPolyhedron::GetNormal: irrelevan << 1160       << "HepPolyhedron::GetNormal: irrelevant index " << iFace 
1709       << std::endl;                              1161       << std::endl;
1710     return G4Normal3D();                      << 1162     return Normal3D<double>();
1711   }                                              1163   }
1712                                                  1164 
1713   G4int i0  = std::abs(pF[iFace].edge[0].v);  << 1165   int i0  = std::abs(pF[iFace].edge[0].v);
1714   G4int i1  = std::abs(pF[iFace].edge[1].v);  << 1166   int i1  = std::abs(pF[iFace].edge[1].v);
1715   G4int i2  = std::abs(pF[iFace].edge[2].v);  << 1167   int i2  = std::abs(pF[iFace].edge[2].v);
1716   G4int i3  = std::abs(pF[iFace].edge[3].v);  << 1168   int i3  = std::abs(pF[iFace].edge[3].v);
1717   if (i3 == 0) i3 = i0;                          1169   if (i3 == 0) i3 = i0;
1718   return (pV[i2] - pV[i0]).cross(pV[i3] - pV[    1170   return (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1]);
1719 }                                                1171 }
1720                                                  1172 
1721 G4Normal3D HepPolyhedron::GetUnitNormal(G4int << 1173 Normal3D<double> HepPolyhedron::GetUnitNormal(int iFace) const
1722 /********************************************    1174 /***********************************************************************
1723  *                                               1175  *                                                                     *
1724  * Name: HepPolyhedron::GetNormal                1176  * Name: HepPolyhedron::GetNormal                    Date:    19.11.99 *
1725  * Author: E.Chernyaev                           1177  * Author: E.Chernyaev                               Revised:          *
1726  *                                               1178  *                                                                     *
1727  * Function: Get unit normal of the face give    1179  * Function: Get unit normal of the face given by index                *
1728  *                                               1180  *                                                                     *
1729  ********************************************    1181  ***********************************************************************/
1730 {                                                1182 {
1731   if (iFace < 1 || iFace > nface) {              1183   if (iFace < 1 || iFace > nface) {
1732     std::cerr                                 << 1184     std::cerr 
1733       << "HepPolyhedron::GetUnitNormal: irrel    1185       << "HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
1734       << std::endl;                              1186       << std::endl;
1735     return G4Normal3D();                      << 1187     return Normal3D<double>();
1736   }                                              1188   }
1737                                                  1189 
1738   G4int i0  = std::abs(pF[iFace].edge[0].v);  << 1190   int i0  = std::abs(pF[iFace].edge[0].v);
1739   G4int i1  = std::abs(pF[iFace].edge[1].v);  << 1191   int i1  = std::abs(pF[iFace].edge[1].v);
1740   G4int i2  = std::abs(pF[iFace].edge[2].v);  << 1192   int i2  = std::abs(pF[iFace].edge[2].v);
1741   G4int i3  = std::abs(pF[iFace].edge[3].v);  << 1193   int i3  = std::abs(pF[iFace].edge[3].v);
1742   if (i3 == 0) i3 = i0;                          1194   if (i3 == 0) i3 = i0;
1743   return ((pV[i2] - pV[i0]).cross(pV[i3] - pV    1195   return ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).unit();
1744 }                                                1196 }
1745                                                  1197 
1746 G4bool HepPolyhedron::GetNextNormal(G4Normal3 << 1198 bool HepPolyhedron::GetNextNormal(Normal3D<double> &normal) const
1747 /********************************************    1199 /***********************************************************************
1748  *                                               1200  *                                                                     *
1749  * Name: HepPolyhedron::GetNextNormal            1201  * Name: HepPolyhedron::GetNextNormal               Date:    22.07.96  *
1750  * Author: John Allison                          1202  * Author: John Allison                             Revised: 19.11.99  *
1751  *                                               1203  *                                                                     *
1752  * Function: Get normals of each face in face    1204  * Function: Get normals of each face in face order.  Returns false    *
1753  *           when finished all faces.            1205  *           when finished all faces.                                  *
1754  *                                               1206  *                                                                     *
1755  ********************************************    1207  ***********************************************************************/
1756 {                                                1208 {
1757   static G4ThreadLocal G4int iFace = 1;       << 1209   static int iFace = 1;
1758   normal = GetNormal(iFace);                     1210   normal = GetNormal(iFace);
1759   if (++iFace > nface) {                         1211   if (++iFace > nface) {
1760     iFace = 1;                                   1212     iFace = 1;
1761     return false;                                1213     return false;
                                                   >> 1214   }else{
                                                   >> 1215     return true;
1762   }                                              1216   }
1763   return true;                                << 
1764 }                                                1217 }
1765                                                  1218 
1766 G4bool HepPolyhedron::GetNextUnitNormal(G4Nor << 1219 bool HepPolyhedron::GetNextUnitNormal(Normal3D<double> &normal) const
1767 /********************************************    1220 /***********************************************************************
1768  *                                               1221  *                                                                     *
1769  * Name: HepPolyhedron::GetNextUnitNormal        1222  * Name: HepPolyhedron::GetNextUnitNormal           Date:    16.09.96  *
1770  * Author: E.Chernyaev                           1223  * Author: E.Chernyaev                              Revised:           *
1771  *                                               1224  *                                                                     *
1772  * Function: Get normals of unit length of ea    1225  * Function: Get normals of unit length of each face in face order.    *
1773  *           Returns false when finished all     1226  *           Returns false when finished all faces.                    *
1774  *                                               1227  *                                                                     *
1775  ********************************************    1228  ***********************************************************************/
1776 {                                                1229 {
1777   G4bool rep = GetNextNormal(normal);         << 1230   bool rep = GetNextNormal(normal);
1778   normal = normal.unit();                        1231   normal = normal.unit();
1779   return rep;                                    1232   return rep;
1780 }                                                1233 }
1781                                                  1234 
1782 G4double HepPolyhedron::GetSurfaceArea() cons << 1235 double HepPolyhedron::GetSurfaceArea() const
1783 /********************************************    1236 /***********************************************************************
1784  *                                               1237  *                                                                     *
1785  * Name: HepPolyhedron::GetSurfaceArea           1238  * Name: HepPolyhedron::GetSurfaceArea              Date:    25.05.01  *
1786  * Author: E.Chernyaev                           1239  * Author: E.Chernyaev                              Revised:           *
1787  *                                               1240  *                                                                     *
1788  * Function: Returns area of the surface of t    1241  * Function: Returns area of the surface of the polyhedron.            *
1789  *                                               1242  *                                                                     *
1790  ********************************************    1243  ***********************************************************************/
1791 {                                                1244 {
1792   G4double srf = 0.;                          << 1245   double s = 0.;
1793   for (G4int iFace=1; iFace<=nface; iFace++)  << 1246   for (int iFace=1; iFace<=nface; iFace++) {
1794     G4int i0 = std::abs(pF[iFace].edge[0].v); << 1247     int i0 = std::abs(pF[iFace].edge[0].v);
1795     G4int i1 = std::abs(pF[iFace].edge[1].v); << 1248     int i1 = std::abs(pF[iFace].edge[1].v);
1796     G4int i2 = std::abs(pF[iFace].edge[2].v); << 1249     int i2 = std::abs(pF[iFace].edge[2].v);
1797     G4int i3 = std::abs(pF[iFace].edge[3].v); << 1250     int i3 = std::abs(pF[iFace].edge[3].v);
1798     if (i3 == 0) i3 = i0;                        1251     if (i3 == 0) i3 = i0;
1799     srf += ((pV[i2] - pV[i0]).cross(pV[i3] -  << 1252     s += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).mag();
1800   }                                              1253   }
1801   return srf/2.;                              << 1254   return s/2.;
1802 }                                                1255 }
1803                                                  1256 
1804 G4double HepPolyhedron::GetVolume() const     << 1257 double HepPolyhedron::GetVolume() const
1805 /********************************************    1258 /***********************************************************************
1806  *                                               1259  *                                                                     *
1807  * Name: HepPolyhedron::GetVolume                1260  * Name: HepPolyhedron::GetVolume                   Date:    25.05.01  *
1808  * Author: E.Chernyaev                           1261  * Author: E.Chernyaev                              Revised:           *
1809  *                                               1262  *                                                                     *
1810  * Function: Returns volume of the polyhedron    1263  * Function: Returns volume of the polyhedron.                         *
1811  *                                               1264  *                                                                     *
1812  ********************************************    1265  ***********************************************************************/
1813 {                                                1266 {
1814   G4double v = 0.;                            << 1267   double v = 0.;
1815   for (G4int iFace=1; iFace<=nface; iFace++)  << 1268   for (int iFace=1; iFace<=nface; iFace++) {
1816     G4int i0 = std::abs(pF[iFace].edge[0].v); << 1269     int i0 = std::abs(pF[iFace].edge[0].v);
1817     G4int i1 = std::abs(pF[iFace].edge[1].v); << 1270     int i1 = std::abs(pF[iFace].edge[1].v);
1818     G4int i2 = std::abs(pF[iFace].edge[2].v); << 1271     int i2 = std::abs(pF[iFace].edge[2].v);
1819     G4int i3 = std::abs(pF[iFace].edge[3].v); << 1272     int i3 = std::abs(pF[iFace].edge[3].v);
1820     G4Point3D pt;                             << 1273     Point3D<double> g;
1821     if (i3 == 0) {                               1274     if (i3 == 0) {
1822       i3 = i0;                                   1275       i3 = i0;
1823       pt = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);  << 1276       g  = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);
1824     }else{                                       1277     }else{
1825       pt = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0. << 1278       g  = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0.25;
1826     }                                            1279     }
1827     v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV << 1280     v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).dot(g);
1828   }                                              1281   }
1829   return v/6.;                                   1282   return v/6.;
1830 }                                                1283 }
1831                                                  1284 
1832 G4int                                         << 1285 int
1833 HepPolyhedron::createTwistedTrap(G4double Dz, << 1286 HepPolyhedron::createTwistedTrap(double Dz,
1834                                  const G4doub << 1287                                  const double xy1[][2],
1835                                  const G4doub << 1288                                  const double xy2[][2])
1836 /********************************************    1289 /***********************************************************************
1837  *                                               1290  *                                                                     *
1838  * Name: createTwistedTrap                       1291  * Name: createTwistedTrap                           Date:    05.11.02 *
1839  * Author: E.Chernyaev (IHEP/Protvino)           1292  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
1840  *                                               1293  *                                                                     *
1841  * Function: Creates polyhedron for twisted t    1294  * Function: Creates polyhedron for twisted trapezoid                  *
1842  *                                               1295  *                                                                     *
1843  * Input: Dz       - half-length along Z         1296  * Input: Dz       - half-length along Z             8----7            *
1844  *        xy1[2,4] - quadrilateral at Z=-Dz      1297  *        xy1[2,4] - quadrilateral at Z=-Dz       5----6  !            *
1845  *        xy2[2,4] - quadrilateral at Z=+Dz      1298  *        xy2[2,4] - quadrilateral at Z=+Dz       !  4-!--3            *
1846  *                                               1299  *                                                1----2               *
1847  *                                               1300  *                                                                     *
1848  ********************************************    1301  ***********************************************************************/
1849 {                                                1302 {
1850   AllocateMemory(12,18);                         1303   AllocateMemory(12,18);
1851                                                  1304 
1852   pV[ 1] = G4Point3D(xy1[0][0],xy1[0][1],-Dz) << 1305   pV[ 1] = Point3D<double>(xy1[0][0],xy1[0][1],-Dz);
1853   pV[ 2] = G4Point3D(xy1[1][0],xy1[1][1],-Dz) << 1306   pV[ 2] = Point3D<double>(xy1[1][0],xy1[1][1],-Dz);
1854   pV[ 3] = G4Point3D(xy1[2][0],xy1[2][1],-Dz) << 1307   pV[ 3] = Point3D<double>(xy1[2][0],xy1[2][1],-Dz);
1855   pV[ 4] = G4Point3D(xy1[3][0],xy1[3][1],-Dz) << 1308   pV[ 4] = Point3D<double>(xy1[3][0],xy1[3][1],-Dz);
1856                                               << 1309 
1857   pV[ 5] = G4Point3D(xy2[0][0],xy2[0][1], Dz) << 1310   pV[ 5] = Point3D<double>(xy2[0][0],xy2[0][1], Dz);
1858   pV[ 6] = G4Point3D(xy2[1][0],xy2[1][1], Dz) << 1311   pV[ 6] = Point3D<double>(xy2[1][0],xy2[1][1], Dz);
1859   pV[ 7] = G4Point3D(xy2[2][0],xy2[2][1], Dz) << 1312   pV[ 7] = Point3D<double>(xy2[2][0],xy2[2][1], Dz);
1860   pV[ 8] = G4Point3D(xy2[3][0],xy2[3][1], Dz) << 1313   pV[ 8] = Point3D<double>(xy2[3][0],xy2[3][1], Dz);
1861                                                  1314 
1862   pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;         1315   pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;
1863   pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;         1316   pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;
1864   pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;         1317   pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;
1865   pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;         1318   pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;
1866                                                  1319 
1867   enum {DUMMY, BOTTOM,                           1320   enum {DUMMY, BOTTOM,
1868         LEFT_BOTTOM,  LEFT_FRONT,   LEFT_TOP,    1321         LEFT_BOTTOM,  LEFT_FRONT,   LEFT_TOP,  LEFT_BACK,
1869         BACK_BOTTOM,  BACK_LEFT,    BACK_TOP,    1322         BACK_BOTTOM,  BACK_LEFT,    BACK_TOP,  BACK_RIGHT,
1870         RIGHT_BOTTOM, RIGHT_BACK,   RIGHT_TOP    1323         RIGHT_BOTTOM, RIGHT_BACK,   RIGHT_TOP, RIGHT_FRONT,
1871         FRONT_BOTTOM, FRONT_RIGHT,  FRONT_TOP    1324         FRONT_BOTTOM, FRONT_RIGHT,  FRONT_TOP, FRONT_LEFT,
1872         TOP};                                    1325         TOP};
1873                                                  1326 
1874   pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM    1327   pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
1875                                                  1328 
1876   pF[ 2]=G4Facet(4,BOTTOM,     -1,LEFT_FRONT,    1329   pF[ 2]=G4Facet(4,BOTTOM,     -1,LEFT_FRONT,  -12,LEFT_BACK,    0,0);
1877   pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP,      1330   pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP,    -12,LEFT_BOTTOM,  0,0);
1878   pF[ 4]=G4Facet(5,TOP,        -8,LEFT_BACK,     1331   pF[ 4]=G4Facet(5,TOP,        -8,LEFT_BACK,   -12,LEFT_FRONT,   0,0);
1879   pF[ 5]=G4Facet(8,BACK_LEFT,  -4,LEFT_BOTTOM    1332   pF[ 5]=G4Facet(8,BACK_LEFT,  -4,LEFT_BOTTOM, -12,LEFT_TOP,     0,0);
1880                                                  1333 
1881   pF[ 6]=G4Facet(3,BOTTOM,     -4,BACK_LEFT,     1334   pF[ 6]=G4Facet(3,BOTTOM,     -4,BACK_LEFT,   -11,BACK_RIGHT,   0,0);
1882   pF[ 7]=G4Facet(4,LEFT_BACK,  -8,BACK_TOP,      1335   pF[ 7]=G4Facet(4,LEFT_BACK,  -8,BACK_TOP,    -11,BACK_BOTTOM,  0,0);
1883   pF[ 8]=G4Facet(8,TOP,        -7,BACK_RIGHT,    1336   pF[ 8]=G4Facet(8,TOP,        -7,BACK_RIGHT,  -11,BACK_LEFT,    0,0);
1884   pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM    1337   pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP,     0,0);
1885                                                  1338 
1886   pF[10]=G4Facet(2,BOTTOM,     -3,RIGHT_BACK,    1339   pF[10]=G4Facet(2,BOTTOM,     -3,RIGHT_BACK,  -10,RIGHT_FRONT,  0,0);
1887   pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP,     1340   pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP,   -10,RIGHT_BOTTOM, 0,0);
1888   pF[12]=G4Facet(7,TOP,        -6,RIGHT_FRONT    1341   pF[12]=G4Facet(7,TOP,        -6,RIGHT_FRONT, -10,RIGHT_BACK,   0,0);
1889   pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTO    1342   pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP,    0,0);
1890                                                  1343 
1891   pF[14]=G4Facet(1,BOTTOM,     -2,FRONT_RIGHT    1344   pF[14]=G4Facet(1,BOTTOM,     -2,FRONT_RIGHT,  -9,FRONT_LEFT,   0,0);
1892   pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP,     1345   pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP,    -9,FRONT_BOTTOM, 0,0);
1893   pF[16]=G4Facet(6,TOP,        -5,FRONT_LEFT,    1346   pF[16]=G4Facet(6,TOP,        -5,FRONT_LEFT,   -9,FRONT_RIGHT,  0,0);
1894   pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTO    1347   pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP,    0,0);
1895                                               << 1348  
1896   pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,    1349   pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
1897                                                  1350 
1898   return 0;                                      1351   return 0;
1899 }                                                1352 }
1900                                                  1353 
1901 G4int                                         << 1354 int
1902 HepPolyhedron::createPolyhedron(G4int Nnodes, << 1355 HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
1903                                 const G4doubl << 1356                                 const double xyz[][3],
1904                                 const G4int   << 1357                                 const int  faces[][4])
1905 /********************************************    1358 /***********************************************************************
1906  *                                               1359  *                                                                     *
1907  * Name: createPolyhedron                        1360  * Name: createPolyhedron                            Date:    05.11.02 *
1908  * Author: E.Chernyaev (IHEP/Protvino)           1361  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
1909  *                                               1362  *                                                                     *
1910  * Function: Creates user defined polyhedron     1363  * Function: Creates user defined polyhedron                           *
1911  *                                               1364  *                                                                     *
1912  * Input: Nnodes  - number of nodes              1365  * Input: Nnodes  - number of nodes                                    *
1913  *        Nfaces  - number of faces              1366  *        Nfaces  - number of faces                                    *
1914  *        nodes[][3] - node coordinates          1367  *        nodes[][3] - node coordinates                                *
1915  *        faces[][4] - faces                     1368  *        faces[][4] - faces                                           *
1916  *                                               1369  *                                                                     *
1917  ********************************************    1370  ***********************************************************************/
1918 {                                                1371 {
1919   AllocateMemory(Nnodes, Nfaces);                1372   AllocateMemory(Nnodes, Nfaces);
1920   if (nvert == 0) return 1;                      1373   if (nvert == 0) return 1;
1921                                                  1374 
1922   for (G4int i=0; i<Nnodes; i++) {            << 1375   for (int i=0; i<Nnodes; i++) {
1923     pV[i+1] = G4Point3D(xyz[i][0], xyz[i][1], << 1376     pV[i+1] = Point3D<double>(xyz[i][0], xyz[i][1], xyz[i][2]);
1924   }                                              1377   }
1925   for (G4int k=0; k<Nfaces; k++) {            << 1378   for (int k=0; k<Nfaces; k++) {
1926     pF[k+1] = G4Facet(faces[k][0],0,faces[k][    1379     pF[k+1] = G4Facet(faces[k][0],0,faces[k][1],0,faces[k][2],0,faces[k][3],0);
1927   }                                              1380   }
1928   SetReferences();                               1381   SetReferences();
1929   return 0;                                      1382   return 0;
1930 }                                                1383 }
1931                                                  1384 
1932 G4Point3D HepPolyhedron::vertexUnweightedMean << 1385 HepPolyhedronTrd2::HepPolyhedronTrd2(double Dx1, double Dx2,
1933   /****************************************** << 1386                                      double Dy1, double Dy2,
1934    *                                          << 1387                                      double Dz)
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 /********************************************    1388 /***********************************************************************
1955  *                                               1389  *                                                                     *
1956  * Name: HepPolyhedronTrd2                       1390  * Name: HepPolyhedronTrd2                           Date:    22.07.96 *
1957  * Author: E.Chernyaev (IHEP/Protvino)           1391  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
1958  *                                               1392  *                                                                     *
1959  * Function: Create GEANT4 TRD2-trapezoid        1393  * Function: Create GEANT4 TRD2-trapezoid                              *
1960  *                                               1394  *                                                                     *
1961  * Input: Dx1 - half-length along X at -Dz       1395  * Input: Dx1 - half-length along X at -Dz           8----7            *
1962  *        Dx2 - half-length along X ay +Dz       1396  *        Dx2 - half-length along X ay +Dz        5----6  !            *
1963  *        Dy1 - half-length along Y ay -Dz       1397  *        Dy1 - half-length along Y ay -Dz        !  4-!--3            *
1964  *        Dy2 - half-length along Y ay +Dz       1398  *        Dy2 - half-length along Y ay +Dz        1----2               *
1965  *        Dz  - half-length along Z              1399  *        Dz  - half-length along Z                                    *
1966  *                                               1400  *                                                                     *
1967  ********************************************    1401  ***********************************************************************/
1968 {                                                1402 {
1969   AllocateMemory(8,6);                           1403   AllocateMemory(8,6);
1970                                                  1404 
1971   pV[1] = G4Point3D(-Dx1,-Dy1,-Dz);           << 1405   pV[1] = Point3D<double>(-Dx1,-Dy1,-Dz);
1972   pV[2] = G4Point3D( Dx1,-Dy1,-Dz);           << 1406   pV[2] = Point3D<double>( Dx1,-Dy1,-Dz);
1973   pV[3] = G4Point3D( Dx1, Dy1,-Dz);           << 1407   pV[3] = Point3D<double>( Dx1, Dy1,-Dz);
1974   pV[4] = G4Point3D(-Dx1, Dy1,-Dz);           << 1408   pV[4] = Point3D<double>(-Dx1, Dy1,-Dz);
1975   pV[5] = G4Point3D(-Dx2,-Dy2, Dz);           << 1409   pV[5] = Point3D<double>(-Dx2,-Dy2, Dz);
1976   pV[6] = G4Point3D( Dx2,-Dy2, Dz);           << 1410   pV[6] = Point3D<double>( Dx2,-Dy2, Dz);
1977   pV[7] = G4Point3D( Dx2, Dy2, Dz);           << 1411   pV[7] = Point3D<double>( Dx2, Dy2, Dz);
1978   pV[8] = G4Point3D(-Dx2, Dy2, Dz);           << 1412   pV[8] = Point3D<double>(-Dx2, Dy2, Dz);
1979                                                  1413 
1980   CreatePrism();                                 1414   CreatePrism();
1981 }                                                1415 }
1982                                                  1416 
1983 HepPolyhedronTrd2::~HepPolyhedronTrd2() = def << 1417 HepPolyhedronTrd2::~HepPolyhedronTrd2() {}
1984                                                  1418 
1985 HepPolyhedronTrd1::HepPolyhedronTrd1(G4double << 1419 HepPolyhedronTrd1::HepPolyhedronTrd1(double Dx1, double Dx2,
1986                                      G4double << 1420                                      double Dy, double Dz)
1987   : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {    1421   : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1988                                                  1422 
1989 HepPolyhedronTrd1::~HepPolyhedronTrd1() = def << 1423 HepPolyhedronTrd1::~HepPolyhedronTrd1() {}
1990                                                  1424 
1991 HepPolyhedronBox::HepPolyhedronBox(G4double D << 1425 HepPolyhedronBox::HepPolyhedronBox(double Dx, double Dy, double Dz)
1992   : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}     1426   : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1993                                                  1427 
1994 HepPolyhedronBox::~HepPolyhedronBox() = defau << 1428 HepPolyhedronBox::~HepPolyhedronBox() {}
1995                                                  1429 
1996 HepPolyhedronTrap::HepPolyhedronTrap(G4double << 1430 HepPolyhedronTrap::HepPolyhedronTrap(double Dz,
1997                                      G4double << 1431                                      double Theta,
1998                                      G4double << 1432                                      double Phi,
1999                                      G4double << 1433                                      double Dy1,
2000                                      G4double << 1434                                      double Dx1,
2001                                      G4double << 1435                                      double Dx2,
2002                                      G4double << 1436                                      double Alp1,
2003                                      G4double << 1437                                      double Dy2,
2004                                      G4double << 1438                                      double Dx3,
2005                                      G4double << 1439                                      double Dx4,
2006                                      G4double << 1440                                      double Alp2)
2007 /********************************************    1441 /***********************************************************************
2008  *                                               1442  *                                                                     *
2009  * Name: HepPolyhedronTrap                       1443  * Name: HepPolyhedronTrap                           Date:    20.11.96 *
2010  * Author: E.Chernyaev                           1444  * Author: E.Chernyaev                               Revised:          *
2011  *                                               1445  *                                                                     *
2012  * Function: Create GEANT4 TRAP-trapezoid        1446  * Function: Create GEANT4 TRAP-trapezoid                              *
2013  *                                               1447  *                                                                     *
2014  * Input: DZ   - half-length in Z                1448  * Input: DZ   - half-length in Z                                      *
2015  *        Theta,Phi - polar angles of the lin    1449  *        Theta,Phi - polar angles of the line joining centres of the  *
2016  *                    faces at Z=-Dz and Z=+D    1450  *                    faces at Z=-Dz and Z=+Dz                         *
2017  *        Dy1  - half-length in Y of the face    1451  *        Dy1  - half-length in Y of the face at Z=-Dz                 *
2018  *        Dx1  - half-length in X of low edge    1452  *        Dx1  - half-length in X of low edge of the face at Z=-Dz     *
2019  *        Dx2  - half-length in X of top edge    1453  *        Dx2  - half-length in X of top edge of the face at Z=-Dz     *
2020  *        Alp1 - angle between Y-axis and the    1454  *        Alp1 - angle between Y-axis and the median joining top and   *
2021  *               low edges of the face at Z=-    1455  *               low edges of the face at Z=-Dz                        *
2022  *        Dy2  - half-length in Y of the face    1456  *        Dy2  - half-length in Y of the face at Z=+Dz                 *
2023  *        Dx3  - half-length in X of low edge    1457  *        Dx3  - half-length in X of low edge of the face at Z=+Dz     *
2024  *        Dx4  - half-length in X of top edge    1458  *        Dx4  - half-length in X of top edge of the face at Z=+Dz     *
2025  *        Alp2 - angle between Y-axis and the    1459  *        Alp2 - angle between Y-axis and the median joining top and   *
2026  *               low edges of the face at Z=+    1460  *               low edges of the face at Z=+Dz                        *
2027  *                                               1461  *                                                                     *
2028  ********************************************    1462  ***********************************************************************/
2029 {                                                1463 {
2030   G4double DzTthetaCphi = Dz*std::tan(Theta)* << 1464   double DzTthetaCphi = Dz*std::tan(Theta)*std::cos(Phi);
2031   G4double DzTthetaSphi = Dz*std::tan(Theta)* << 1465   double DzTthetaSphi = Dz*std::tan(Theta)*std::sin(Phi);
2032   G4double Dy1Talp1 = Dy1*std::tan(Alp1);     << 1466   double Dy1Talp1 = Dy1*std::tan(Alp1);
2033   G4double Dy2Talp2 = Dy2*std::tan(Alp2);     << 1467   double Dy2Talp2 = Dy2*std::tan(Alp2);
2034                                               << 1468   
2035   AllocateMemory(8,6);                           1469   AllocateMemory(8,6);
2036                                                  1470 
2037   pV[1] = G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx << 1471   pV[1] = Point3D<double>(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
2038   pV[2] = G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx << 1472   pV[2] = Point3D<double>(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
2039   pV[3] = G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx << 1473   pV[3] = Point3D<double>(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
2040   pV[4] = G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx << 1474   pV[4] = Point3D<double>(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
2041   pV[5] = G4Point3D( DzTthetaCphi-Dy2Talp2-Dx << 1475   pV[5] = Point3D<double>( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
2042   pV[6] = G4Point3D( DzTthetaCphi-Dy2Talp2+Dx << 1476   pV[6] = Point3D<double>( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
2043   pV[7] = G4Point3D( DzTthetaCphi+Dy2Talp2+Dx << 1477   pV[7] = Point3D<double>( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
2044   pV[8] = G4Point3D( DzTthetaCphi+Dy2Talp2-Dx << 1478   pV[8] = Point3D<double>( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
2045                                                  1479 
2046   CreatePrism();                                 1480   CreatePrism();
2047 }                                                1481 }
2048                                                  1482 
2049 HepPolyhedronTrap::~HepPolyhedronTrap() = def << 1483 HepPolyhedronTrap::~HepPolyhedronTrap() {}
2050                                                  1484 
2051 HepPolyhedronPara::HepPolyhedronPara(G4double << 1485 HepPolyhedronPara::HepPolyhedronPara(double Dx, double Dy, double Dz,
2052                                      G4double << 1486                                      double Alpha, double Theta,
2053                                      G4double << 1487                                      double Phi)
2054   : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx,    1488   : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
2055                                                  1489 
2056 HepPolyhedronPara::~HepPolyhedronPara() = def << 1490 HepPolyhedronPara::~HepPolyhedronPara() {}
2057                                                  1491 
2058 HepPolyhedronParaboloid::HepPolyhedronParabol << 1492 HepPolyhedronCons::HepPolyhedronCons(double Rmn1,
2059                                               << 1493                                      double Rmx1,
2060                                               << 1494                                      double Rmn2,
2061                                               << 1495                                      double Rmx2, 
2062                                               << 1496                                      double Dz,
2063 /******************************************** << 1497                                      double Phi1,
2064  *                                            << 1498                                      double Dphi) 
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 /********************************************    1499 /***********************************************************************
2250  *                                               1500  *                                                                     *
2251  * Name: HepPolyhedronCons::HepPolyhedronCons    1501  * Name: HepPolyhedronCons::HepPolyhedronCons        Date:    15.12.96 *
2252  * Author: E.Chernyaev (IHEP/Protvino)           1502  * Author: E.Chernyaev (IHEP/Protvino)               Revised: 15.12.96 *
2253  *                                               1503  *                                                                     *
2254  * Function: Constructor for CONS, TUBS, CONE    1504  * Function: Constructor for CONS, TUBS, CONE, TUBE                    *
2255  *                                               1505  *                                                                     *
2256  * Input: Rmn1, Rmx1 - inside and outside rad    1506  * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz              *
2257  *        Rmn2, Rmx2 - inside and outside rad    1507  *        Rmn2, Rmx2 - inside and outside radiuses at +Dz              *
2258  *        Dz         - half length in Z          1508  *        Dz         - half length in Z                                *
2259  *        Phi1       - starting angle of the     1509  *        Phi1       - starting angle of the segment                   *
2260  *        Dphi       - segment range             1510  *        Dphi       - segment range                                   *
2261  *                                               1511  *                                                                     *
2262  ********************************************    1512  ***********************************************************************/
2263 {                                                1513 {
2264   static const G4double wholeCircle=twopi;    << 1514   static double wholeCircle=twopi;
2265                                                  1515 
2266   //   C H E C K   I N P U T   P A R A M E T     1516   //   C H E C K   I N P U T   P A R A M E T E R S
2267                                                  1517 
2268   G4int k = 0;                                << 1518   int k = 0;
2269   if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. ||     1519   if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.)        k = 1;
2270   if (Rmn1 > Rmx1 || Rmn2 > Rmx2)                1520   if (Rmn1 > Rmx1 || Rmn2 > Rmx2)                              k = 1;
2271   if (Rmn1 == Rmx1 && Rmn2 == Rmx2)              1521   if (Rmn1 == Rmx1 && Rmn2 == Rmx2)                            k = 1;
2272                                                  1522 
2273   if (Dz <= 0.) k += 2;                          1523   if (Dz <= 0.) k += 2;
2274                                               << 1524  
2275   G4double phi1, phi2, dphi;                  << 1525   double phi1, phi2, dphi;
2276   if (Dphi < 0.) {                               1526   if (Dphi < 0.) {
2277     phi2 = Phi1; phi1 = phi2 - Dphi;             1527     phi2 = Phi1; phi1 = phi2 - Dphi;
2278   }else if (Dphi == 0.) {                        1528   }else if (Dphi == 0.) {
2279     phi1 = Phi1; phi2 = phi1 + wholeCircle;      1529     phi1 = Phi1; phi2 = phi1 + wholeCircle;
2280   }else{                                         1530   }else{
2281     phi1 = Phi1; phi2 = phi1 + Dphi;             1531     phi1 = Phi1; phi2 = phi1 + Dphi;
2282   }                                              1532   }
2283   dphi  = phi2 - phi1;                           1533   dphi  = phi2 - phi1;
2284   if (std::abs(dphi-wholeCircle) < perMillion    1534   if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
2285   if (dphi > wholeCircle) k += 4;             << 1535   if (dphi > wholeCircle) k += 4; 
2286                                                  1536 
2287   if (k != 0) {                                  1537   if (k != 0) {
2288     std::cerr << "HepPolyhedronCone(s)/Tube(s    1538     std::cerr << "HepPolyhedronCone(s)/Tube(s): error in input parameters";
2289     if ((k & 1) != 0) std::cerr << " (radiuse    1539     if ((k & 1) != 0) std::cerr << " (radiuses)";
2290     if ((k & 2) != 0) std::cerr << " (half-le    1540     if ((k & 2) != 0) std::cerr << " (half-length)";
2291     if ((k & 4) != 0) std::cerr << " (angles)    1541     if ((k & 4) != 0) std::cerr << " (angles)";
2292     std::cerr << std::endl;                      1542     std::cerr << std::endl;
2293     std::cerr << " Rmn1=" << Rmn1 << " Rmx1="    1543     std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
2294     std::cerr << " Rmn2=" << Rmn2 << " Rmx2="    1544     std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
2295     std::cerr << " Dz=" << Dz << " Phi1=" <<     1545     std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
2296               << std::endl;                      1546               << std::endl;
2297     return;                                      1547     return;
2298   }                                              1548   }
2299                                               << 1549   
2300   //   P R E P A R E   T W O   P O L Y L I N     1550   //   P R E P A R E   T W O   P O L Y L I N E S
2301                                                  1551 
2302   G4double zz[4], rr[4];                      << 1552   double zz[4], rr[4];
2303   zz[0] =  Dz;                                << 1553   zz[0] =  Dz; 
2304   zz[1] = -Dz;                                << 1554   zz[1] = -Dz; 
2305   zz[2] =  Dz;                                << 1555   zz[2] =  Dz; 
2306   zz[3] = -Dz;                                << 1556   zz[3] = -Dz; 
2307   rr[0] =  Rmx2;                                 1557   rr[0] =  Rmx2;
2308   rr[1] =  Rmx1;                                 1558   rr[1] =  Rmx1;
2309   rr[2] =  Rmn2;                                 1559   rr[2] =  Rmn2;
2310   rr[3] =  Rmn1;                                 1560   rr[3] =  Rmn1;
2311                                                  1561 
2312   //   R O T A T E    P O L Y L I N E S          1562   //   R O T A T E    P O L Y L I N E S
2313                                                  1563 
2314   RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr,  << 1564   RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1); 
2315   SetReferences();                               1565   SetReferences();
2316 }                                                1566 }
2317                                                  1567 
2318 HepPolyhedronCons::~HepPolyhedronCons() = def << 1568 HepPolyhedronCons::~HepPolyhedronCons() {}
2319                                                  1569 
2320 HepPolyhedronCone::HepPolyhedronCone(G4double << 1570 HepPolyhedronCone::HepPolyhedronCone(double Rmn1, double Rmx1, 
2321                                      G4double << 1571                                      double Rmn2, double Rmx2,
2322                                      G4double << 1572                                      double Dz) :
2323   HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, D    1573   HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
2324                                                  1574 
2325 HepPolyhedronCone::~HepPolyhedronCone() = def << 1575 HepPolyhedronCone::~HepPolyhedronCone() {}
2326                                                  1576 
2327 HepPolyhedronTubs::HepPolyhedronTubs(G4double << 1577 HepPolyhedronTubs::HepPolyhedronTubs(double Rmin, double Rmax,
2328                                      G4double << 1578                                      double Dz, 
2329                                      G4double << 1579                                      double Phi1, double Dphi)
2330   :   HepPolyhedronCons(Rmin, Rmax, Rmin, Rma    1580   :   HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
2331                                                  1581 
2332 HepPolyhedronTubs::~HepPolyhedronTubs() = def << 1582 HepPolyhedronTubs::~HepPolyhedronTubs() {}
2333                                                  1583 
2334 HepPolyhedronTube::HepPolyhedronTube (G4doubl << 1584 HepPolyhedronTube::HepPolyhedronTube (double Rmin, double Rmax,
2335                                       G4doubl << 1585                                       double Dz)
2336   : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax,    1586   : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
2337                                                  1587 
2338 HepPolyhedronTube::~HepPolyhedronTube () = de << 1588 HepPolyhedronTube::~HepPolyhedronTube () {}
2339                                                  1589 
2340 HepPolyhedronPgon::HepPolyhedronPgon(G4double << 1590 HepPolyhedronPgon::HepPolyhedronPgon(double phi,
2341                                      G4double << 1591                                      double dphi,
2342                                      G4int np << 1592                                      int    npdv,
2343                                      G4int nz << 1593                                      int    nz,
2344                                      const G4 << 1594                                      const double *z,
2345                                      const G4 << 1595                                      const double *rmin,
2346                                      const G4 << 1596                                      const double *rmax)
2347 /********************************************    1597 /***********************************************************************
2348  *                                               1598  *                                                                     *
2349  * Name: HepPolyhedronPgon                       1599  * Name: HepPolyhedronPgon                           Date:    09.12.96 *
2350  * Author: E.Chernyaev                           1600  * Author: E.Chernyaev                               Revised:          *
2351  *                                               1601  *                                                                     *
2352  * Function: Constructor of polyhedron for PG    1602  * Function: Constructor of polyhedron for PGON, PCON                  *
2353  *                                               1603  *                                                                     *
2354  * Input: phi  - initial phi                     1604  * Input: phi  - initial phi                                           *
2355  *        dphi - delta phi                       1605  *        dphi - delta phi                                             *
2356  *        npdv - number of steps along phi       1606  *        npdv - number of steps along phi                             *
2357  *        nz   - number of z-planes (at least    1607  *        nz   - number of z-planes (at least two)                     *
2358  *        z[]  - z coordinates of the slices     1608  *        z[]  - z coordinates of the slices                           *
2359  *        rmin[] - smaller r at the slices       1609  *        rmin[] - smaller r at the slices                             *
2360  *        rmax[] - bigger  r at the slices       1610  *        rmax[] - bigger  r at the slices                             *
2361  *                                               1611  *                                                                     *
2362  ********************************************    1612  ***********************************************************************/
2363 {                                                1613 {
2364   //   C H E C K   I N P U T   P A R A M E T     1614   //   C H E C K   I N P U T   P A R A M E T E R S
2365                                                  1615 
2366   if (dphi <= 0. || dphi > twopi) {              1616   if (dphi <= 0. || dphi > twopi) {
2367     std::cerr                                    1617     std::cerr
2368       << "HepPolyhedronPgon/Pcon: wrong delta    1618       << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
2369       << std::endl;                              1619       << std::endl;
2370     return;                                      1620     return;
2371   }                                           << 1621   }    
2372                                               << 1622     
2373   if (nz < 2) {                                  1623   if (nz < 2) {
2374     std::cerr                                    1624     std::cerr
2375       << "HepPolyhedronPgon/Pcon: number of z    1625       << "HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
2376       << std::endl;                              1626       << std::endl;
2377     return;                                      1627     return;
2378   }                                              1628   }
2379                                                  1629 
2380   if (npdv < 0) {                                1630   if (npdv < 0) {
2381     std::cerr                                    1631     std::cerr
2382       << "HepPolyhedronPgon/Pcon: error in nu    1632       << "HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
2383       << std::endl;                              1633       << std::endl;
2384     return;                                      1634     return;
2385   }                                              1635   }
2386                                                  1636 
2387   G4int i;                                    << 1637   int i;
2388   for (i=0; i<nz; i++) {                         1638   for (i=0; i<nz; i++) {
2389     if (rmin[i] < 0. || rmax[i] < 0. || rmin[    1639     if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
2390       std::cerr                                  1640       std::cerr
2391         << "HepPolyhedronPgon: error in radiu    1641         << "HepPolyhedronPgon: error in radiuses rmin[" << i << "]="
2392         << rmin[i] << " rmax[" << i << "]=" <    1642         << rmin[i] << " rmax[" << i << "]=" << rmax[i]
2393         << std::endl;                            1643         << std::endl;
2394       return;                                    1644       return;
2395     }                                            1645     }
2396   }                                              1646   }
2397                                                  1647 
2398   //   P R E P A R E   T W O   P O L Y L I N     1648   //   P R E P A R E   T W O   P O L Y L I N E S
2399                                                  1649 
2400   G4double *zz, *rr;                          << 1650   double *zz, *rr;
2401   zz = new G4double[2*nz];                    << 1651   zz = new double[2*nz];
2402   rr = new G4double[2*nz];                    << 1652   rr = new double[2*nz];
2403                                                  1653 
2404   if (z[0] > z[nz-1]) {                          1654   if (z[0] > z[nz-1]) {
2405     for (i=0; i<nz; i++) {                       1655     for (i=0; i<nz; i++) {
2406       zz[i]    = z[i];                           1656       zz[i]    = z[i];
2407       rr[i]    = rmax[i];                        1657       rr[i]    = rmax[i];
2408       zz[i+nz] = z[i];                           1658       zz[i+nz] = z[i];
2409       rr[i+nz] = rmin[i];                        1659       rr[i+nz] = rmin[i];
2410     }                                            1660     }
2411   }else{                                         1661   }else{
2412     for (i=0; i<nz; i++) {                       1662     for (i=0; i<nz; i++) {
2413       zz[i]    = z[nz-i-1];                      1663       zz[i]    = z[nz-i-1];
2414       rr[i]    = rmax[nz-i-1];                   1664       rr[i]    = rmax[nz-i-1];
2415       zz[i+nz] = z[nz-i-1];                      1665       zz[i+nz] = z[nz-i-1];
2416       rr[i+nz] = rmin[nz-i-1];                   1666       rr[i+nz] = rmin[nz-i-1];
2417     }                                            1667     }
2418   }                                              1668   }
2419                                                  1669 
2420   //   R O T A T E    P O L Y L I N E S          1670   //   R O T A T E    P O L Y L I N E S
2421                                                  1671 
2422   G4int nodeVis = 1;                          << 1672   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();                               1673   SetReferences();
2426                                               << 1674   
2427   delete [] zz;                                  1675   delete [] zz;
2428   delete [] rr;                                  1676   delete [] rr;
2429 }                                                1677 }
2430                                                  1678 
2431 HepPolyhedronPgon::HepPolyhedronPgon(G4double << 1679 HepPolyhedronPgon::~HepPolyhedronPgon() {}
2432                                      G4double << 
2433                                      G4int np << 
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                                                  1680 
2483 HepPolyhedronPcon::HepPolyhedronPcon(G4double << 1681 HepPolyhedronPcon::HepPolyhedronPcon(double phi, double dphi, int nz,
2484                                      const G4 << 1682                                      const double *z,
2485                                      const G4 << 1683                                      const double *rmin,
2486                                      const G4 << 1684                                      const double *rmax)
2487   : HepPolyhedronPgon(phi, dphi, 0, nz, z, rm    1685   : HepPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
2488                                                  1686 
2489 HepPolyhedronPcon::HepPolyhedronPcon(G4double << 1687 HepPolyhedronPcon::~HepPolyhedronPcon() {}
2490                                      const st << 1688 
2491   : HepPolyhedronPgon(phi, dphi, 0, rz) {}    << 1689 HepPolyhedronSphere::HepPolyhedronSphere(double rmin, double rmax,
2492                                               << 1690                                        double phi, double dphi,
2493 HepPolyhedronPcon::~HepPolyhedronPcon() = def << 1691                                        double the, double dthe)
2494                                               << 
2495 HepPolyhedronSphere::HepPolyhedronSphere(G4do << 
2496                                          G4do << 
2497                                          G4do << 
2498 /********************************************    1692 /***********************************************************************
2499  *                                               1693  *                                                                     *
2500  * Name: HepPolyhedronSphere                     1694  * Name: HepPolyhedronSphere                         Date:    11.12.96 *
2501  * Author: E.Chernyaev (IHEP/Protvino)           1695  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
2502  *                                               1696  *                                                                     *
2503  * Function: Constructor of polyhedron for SP    1697  * Function: Constructor of polyhedron for SPHERE                      *
2504  *                                               1698  *                                                                     *
2505  * Input: rmin - internal radius                 1699  * Input: rmin - internal radius                                       *
2506  *        rmax - external radius                 1700  *        rmax - external radius                                       *
2507  *        phi  - initial phi                     1701  *        phi  - initial phi                                           *
2508  *        dphi - delta phi                       1702  *        dphi - delta phi                                             *
2509  *        the  - initial theta                   1703  *        the  - initial theta                                         *
2510  *        dthe - delta theta                     1704  *        dthe - delta theta                                           *
2511  *                                               1705  *                                                                     *
2512  ********************************************    1706  ***********************************************************************/
2513 {                                                1707 {
2514   //   C H E C K   I N P U T   P A R A M E T     1708   //   C H E C K   I N P U T   P A R A M E T E R S
2515                                                  1709 
2516   if (dphi <= 0. || dphi > twopi) {              1710   if (dphi <= 0. || dphi > twopi) {
2517     std::cerr                                    1711     std::cerr
2518       << "HepPolyhedronSphere: wrong delta ph    1712       << "HepPolyhedronSphere: wrong delta phi = " << dphi
2519       << std::endl;                              1713       << std::endl;
2520     return;                                      1714     return;
2521   }                                           << 1715   }    
2522                                                  1716 
2523   if (the < 0. || the > pi) {                    1717   if (the < 0. || the > pi) {
2524     std::cerr                                    1718     std::cerr
2525       << "HepPolyhedronSphere: wrong theta =     1719       << "HepPolyhedronSphere: wrong theta = " << the
2526       << std::endl;                              1720       << std::endl;
2527     return;                                      1721     return;
2528   }                                           << 1722   }    
2529                                               << 1723   
2530   if (dthe <= 0. || dthe > pi) {                 1724   if (dthe <= 0. || dthe > pi) {
2531     std::cerr                                    1725     std::cerr
2532       << "HepPolyhedronSphere: wrong delta th    1726       << "HepPolyhedronSphere: wrong delta theta = " << dthe
2533       << std::endl;                              1727       << std::endl;
2534     return;                                      1728     return;
2535   }                                           << 1729   }    
2536                                                  1730 
2537   if (the+dthe > pi) {                           1731   if (the+dthe > pi) {
2538     std::cerr                                    1732     std::cerr
2539       << "HepPolyhedronSphere: wrong theta +     1733       << "HepPolyhedronSphere: wrong theta + delta theta = "
2540       << the << " " << dthe                      1734       << the << " " << dthe
2541       << std::endl;                              1735       << std::endl;
2542     return;                                      1736     return;
2543   }                                           << 1737   }    
2544                                               << 1738   
2545   if (rmin < 0. || rmin >= rmax) {               1739   if (rmin < 0. || rmin >= rmax) {
2546     std::cerr                                    1740     std::cerr
2547       << "HepPolyhedronSphere: error in radiu    1741       << "HepPolyhedronSphere: error in radiuses"
2548       << " rmin=" << rmin << " rmax=" << rmax    1742       << " rmin=" << rmin << " rmax=" << rmax
2549       << std::endl;                              1743       << std::endl;
2550     return;                                      1744     return;
2551   }                                              1745   }
2552                                                  1746 
2553   //   P R E P A R E   T W O   P O L Y L I N     1747   //   P R E P A R E   T W O   P O L Y L I N E S
2554                                                  1748 
2555   G4int nds = (GetNumberOfRotationSteps() + 1 << 1749   int ns = (GetNumberOfRotationSteps() + 1) / 2;
2556   G4int np1 = G4int(dthe*nds/pi+.5) + 1;      << 1750   int np1 = int(dthe*ns/pi+.5) + 1;
2557   if (np1 <= 1) np1 = 2;                         1751   if (np1 <= 1) np1 = 2;
2558   G4int np2 = rmin < spatialTolerance ? 1 : n << 1752   int np2 = rmin < perMillion ? 1 : np1;
2559                                                  1753 
2560   G4double *zz, *rr;                          << 1754   double *zz, *rr;
2561   zz = new G4double[np1+np2];                 << 1755   zz = new double[np1+np2];
2562   rr = new G4double[np1+np2];                 << 1756   rr = new double[np1+np2];
2563                                               << 1757 
2564   G4double a = dthe/(np1-1);                  << 1758   double a = dthe/(np1-1);
2565   G4double cosa, sina;                        << 1759   double cosa, sina;
2566   for (G4int i=0; i<np1; i++) {               << 1760   for (int i=0; i<np1; i++) {
2567     cosa  = std::cos(the+i*a);                   1761     cosa  = std::cos(the+i*a);
2568     sina  = std::sin(the+i*a);                   1762     sina  = std::sin(the+i*a);
2569     zz[i] = rmax*cosa;                           1763     zz[i] = rmax*cosa;
2570     rr[i] = rmax*sina;                           1764     rr[i] = rmax*sina;
2571     if (np2 > 1) {                               1765     if (np2 > 1) {
2572       zz[i+np1] = rmin*cosa;                     1766       zz[i+np1] = rmin*cosa;
2573       rr[i+np1] = rmin*sina;                     1767       rr[i+np1] = rmin*sina;
2574     }                                            1768     }
2575   }                                              1769   }
2576   if (np2 == 1) {                                1770   if (np2 == 1) {
2577     zz[np1] = 0.;                                1771     zz[np1] = 0.;
2578     rr[np1] = 0.;                                1772     rr[np1] = 0.;
2579   }                                              1773   }
2580                                                  1774 
2581   //   R O T A T E    P O L Y L I N E S          1775   //   R O T A T E    P O L Y L I N E S
2582                                                  1776 
2583   RotateAroundZ(0, phi, dphi, np1, np2, zz, r << 1777   RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1); 
2584   SetReferences();                               1778   SetReferences();
2585                                               << 1779   
2586   delete [] zz;                                  1780   delete [] zz;
2587   delete [] rr;                                  1781   delete [] rr;
2588 }                                                1782 }
2589                                                  1783 
2590 HepPolyhedronSphere::~HepPolyhedronSphere() = << 1784 HepPolyhedronSphere::~HepPolyhedronSphere() {}
2591                                                  1785 
2592 HepPolyhedronTorus::HepPolyhedronTorus(G4doub << 1786 HepPolyhedronTorus::HepPolyhedronTorus(double rmin,
2593                                        G4doub << 1787                                        double rmax,
2594                                        G4doub << 1788                                        double rtor,
2595                                        G4doub << 1789                                        double phi,
2596                                        G4doub << 1790                                        double dphi)
2597 /********************************************    1791 /***********************************************************************
2598  *                                               1792  *                                                                     *
2599  * Name: HepPolyhedronTorus                      1793  * Name: HepPolyhedronTorus                          Date:    11.12.96 *
2600  * Author: E.Chernyaev (IHEP/Protvino)           1794  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
2601  *                                               1795  *                                                                     *
2602  * Function: Constructor of polyhedron for TO    1796  * Function: Constructor of polyhedron for TORUS                       *
2603  *                                               1797  *                                                                     *
2604  * Input: rmin - internal radius                 1798  * Input: rmin - internal radius                                       *
2605  *        rmax - external radius                 1799  *        rmax - external radius                                       *
2606  *        rtor - radius of torus                 1800  *        rtor - radius of torus                                       *
2607  *        phi  - initial phi                     1801  *        phi  - initial phi                                           *
2608  *        dphi - delta phi                       1802  *        dphi - delta phi                                             *
2609  *                                               1803  *                                                                     *
2610  ********************************************    1804  ***********************************************************************/
2611 {                                                1805 {
2612   //   C H E C K   I N P U T   P A R A M E T     1806   //   C H E C K   I N P U T   P A R A M E T E R S
2613                                                  1807 
2614   if (dphi <= 0. || dphi > twopi) {              1808   if (dphi <= 0. || dphi > twopi) {
2615     std::cerr                                    1809     std::cerr
2616       << "HepPolyhedronTorus: wrong delta phi    1810       << "HepPolyhedronTorus: wrong delta phi = " << dphi
2617       << std::endl;                              1811       << std::endl;
2618     return;                                      1812     return;
2619   }                                              1813   }
2620                                                  1814 
2621   if (rmin < 0. || rmin >= rmax || rmax >= rt    1815   if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2622     std::cerr                                    1816     std::cerr
2623       << "HepPolyhedronTorus: error in radius    1817       << "HepPolyhedronTorus: error in radiuses"
2624       << " rmin=" << rmin << " rmax=" << rmax    1818       << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
2625       << std::endl;                              1819       << std::endl;
2626     return;                                      1820     return;
2627   }                                              1821   }
2628                                                  1822 
2629   //   P R E P A R E   T W O   P O L Y L I N     1823   //   P R E P A R E   T W O   P O L Y L I N E S
2630                                                  1824 
2631   G4int np1 = GetNumberOfRotationSteps();     << 1825   int np1 = GetNumberOfRotationSteps();
2632   G4int np2 = rmin < spatialTolerance ? 1 : n << 1826   int np2 = rmin < perMillion ? 1 : np1;
2633                                                  1827 
2634   G4double *zz, *rr;                          << 1828   double *zz, *rr;
2635   zz = new G4double[np1+np2];                 << 1829   zz = new double[np1+np2];
2636   rr = new G4double[np1+np2];                 << 1830   rr = new double[np1+np2];
2637                                               << 1831 
2638   G4double a = twopi/np1;                     << 1832   double a = twopi/np1;
2639   G4double cosa, sina;                        << 1833   double cosa, sina;
2640   for (G4int i=0; i<np1; i++) {               << 1834   for (int i=0; i<np1; i++) {
2641     cosa  = std::cos(i*a);                       1835     cosa  = std::cos(i*a);
2642     sina  = std::sin(i*a);                       1836     sina  = std::sin(i*a);
2643     zz[i] = rmax*cosa;                           1837     zz[i] = rmax*cosa;
2644     rr[i] = rtor+rmax*sina;                      1838     rr[i] = rtor+rmax*sina;
2645     if (np2 > 1) {                               1839     if (np2 > 1) {
2646       zz[i+np1] = rmin*cosa;                     1840       zz[i+np1] = rmin*cosa;
2647       rr[i+np1] = rtor+rmin*sina;                1841       rr[i+np1] = rtor+rmin*sina;
2648     }                                            1842     }
2649   }                                              1843   }
2650   if (np2 == 1) {                                1844   if (np2 == 1) {
2651     zz[np1] = 0.;                                1845     zz[np1] = 0.;
2652     rr[np1] = rtor;                              1846     rr[np1] = rtor;
2653     np2 = -1;                                    1847     np2 = -1;
2654   }                                              1848   }
2655                                                  1849 
2656   //   R O T A T E    P O L Y L I N E S          1850   //   R O T A T E    P O L Y L I N E S
2657                                                  1851 
2658   RotateAroundZ(0, phi, dphi, -np1, -np2, zz, << 1852   RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1); 
2659   SetReferences();                               1853   SetReferences();
2660                                               << 1854   
2661   delete [] zz;                                  1855   delete [] zz;
2662   delete [] rr;                                  1856   delete [] rr;
2663 }                                                1857 }
2664                                                  1858 
2665 HepPolyhedronTorus::~HepPolyhedronTorus() = d << 1859 HepPolyhedronTorus::~HepPolyhedronTorus() {}
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                                                  1860 
2707 HepPolyhedronEllipsoid::HepPolyhedronEllipsoi << 1861 HepPolyhedronEllipsoid::HepPolyhedronEllipsoid(double ax, double by,
2708                                               << 1862                                                double cz, double zCut1,
2709                                               << 1863                                                double zCut2)
2710 /********************************************    1864 /***********************************************************************
2711  *                                               1865  *                                                                     *
2712  * Name: HepPolyhedronEllipsoid                  1866  * Name: HepPolyhedronEllipsoid                      Date:    25.02.05 *
2713  * Author: G.Guerrieri                           1867  * Author: G.Guerrieri                               Revised:          *
2714  *         Evgueni Tcherniaev                 << 
2715  *                                               1868  *                                                                     *
2716  * Function: Constructor of polyhedron for EL    1869  * Function: Constructor of polyhedron for ELLIPSOID                   *
2717  *                                               1870  *                                                                     *
2718  * Input: ax - semiaxis x                        1871  * Input: ax - semiaxis x                                              *
2719  *        by - semiaxis y                        1872  *        by - semiaxis y                                              *
2720  *        cz - semiaxis z                        1873  *        cz - semiaxis z                                              *
2721  *        zCut1 - lower cut plane level (soli    1874  *        zCut1 - lower cut plane level (solid lies above this plane)  *
2722  *        zCut2 - upper cut plane level (soli    1875  *        zCut2 - upper cut plane level (solid lies below this plane)  *
2723  *                                               1876  *                                                                     *
2724  ********************************************    1877  ***********************************************************************/
2725 {                                                1878 {
2726   //   C H E C K   I N P U T   P A R A M E T     1879   //   C H E C K   I N P U T   P A R A M E T E R S
2727                                                  1880 
2728   if (zCut1 >= cz || zCut2 <= -cz || zCut1 >     1881   if (zCut1 >= cz || zCut2 <= -cz || zCut1 > zCut2) {
2729     std::cerr << "HepPolyhedronEllipsoid: wro    1882     std::cerr << "HepPolyhedronEllipsoid: wrong zCut1 = " << zCut1
2730            << " zCut2 = " << zCut2               1883            << " zCut2 = " << zCut2
2731            << " for given cz = " << cz << std    1884            << " for given cz = " << cz << std::endl;
2732     return;                                      1885     return;
2733   }                                              1886   }
2734   if (cz <= 0.0) {                               1887   if (cz <= 0.0) {
2735     std::cerr << "HepPolyhedronEllipsoid: bad    1888     std::cerr << "HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
2736       << std::endl;                              1889       << std::endl;
2737     return;                                      1890     return;
2738   }                                              1891   }
2739                                                  1892 
                                                   >> 1893   double dthe;
                                                   >> 1894   double sthe;
                                                   >> 1895   int cutflag;
                                                   >> 1896   cutflag= 0;
                                                   >> 1897   if (zCut2 >= cz)
                                                   >> 1898     {
                                                   >> 1899       sthe= 0.0;
                                                   >> 1900     }
                                                   >> 1901   else
                                                   >> 1902     {
                                                   >> 1903       sthe= std::acos(zCut2/cz);
                                                   >> 1904       cutflag++;
                                                   >> 1905     }
                                                   >> 1906   if (zCut1 <= -cz)
                                                   >> 1907     {
                                                   >> 1908       dthe= pi - sthe;
                                                   >> 1909     }
                                                   >> 1910   else
                                                   >> 1911     {
                                                   >> 1912       dthe= std::acos(zCut1/cz)-sthe;
                                                   >> 1913       cutflag++;
                                                   >> 1914     }
                                                   >> 1915 
2740   //   P R E P A R E   T W O   P O L Y L I N     1916   //   P R E P A R E   T W O   P O L Y L I N E S
2741   //   generate sphere of radius cz first, th    1917   //   generate sphere of radius cz first, then rescale x and y later
2742                                                  1918 
2743   G4double sthe = std::acos(zCut2/cz);        << 1919   int ns = (GetNumberOfRotationSteps() + 1) / 2;
2744   G4double dthe = std::acos(zCut1/cz) - sthe; << 1920   int np1 = int(dthe*ns/pi) + 2 + cutflag;
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                                                  1921 
2750   G4double *zz, *rr;                          << 1922   double *zz, *rr;
2751   zz = new G4double[np1 + np2];               << 1923   zz = new double[np1+1];
2752   rr = new G4double[np1 + np2];               << 1924   rr = new double[np1+1];
2753   if ((zz == nullptr) || (rr == nullptr))     << 1925   if (!zz || !rr)
2754   {                                           << 1926     {
2755     G4Exception("HepPolyhedronEllipsoid::HepP << 1927       std::cerr << "Out of memory in HepPolyhedronEllipsoid!" << std::endl;
2756                 "greps1002", FatalException,  << 1928         //Exception("Out of memory in HepPolyhedronEllipsoid!");
2757   }                                           << 1929     }
2758                                                  1930 
2759   G4double a = dthe/(np1 - 1);                << 1931   double a = dthe/(np1-cutflag-1);
2760   G4double cosa, sina;                        << 1932   double cosa, sina;
2761   for (G4int i = 0; i < np1; ++i)             << 1933   int j=0;
2762   {                                           << 1934   if (sthe > 0.0)
2763     cosa  = std::cos(sthe + i*a);             << 1935     {
2764     sina  = std::sin(sthe + i*a);             << 1936       zz[j]= zCut2;
2765     zz[i] = cz*cosa;                          << 1937       rr[j]= 0.;
2766     rr[i] = cz*sina;                          << 1938       j++;
2767   }                                           << 1939     }
2768   zz[np1 + 0] = zCut2;                        << 1940   for (int i=0; i<np1-cutflag; i++) {
2769   rr[np1 + 0] = 0.;                           << 1941     cosa  = std::cos(sthe+i*a);
2770   zz[np1 + 1] = zCut1;                        << 1942     sina  = std::sin(sthe+i*a);
2771   rr[np1 + 1] = 0.;                           << 1943     zz[j] = cz*cosa;
                                                   >> 1944     rr[j] = cz*sina;
                                                   >> 1945     j++;
                                                   >> 1946   }
                                                   >> 1947   if (j < np1)
                                                   >> 1948     {
                                                   >> 1949       zz[j]= zCut1;
                                                   >> 1950       rr[j]= 0.;
                                                   >> 1951       j++;
                                                   >> 1952     }
                                                   >> 1953   if (j > np1)
                                                   >> 1954     {
                                                   >> 1955       std::cerr << "Logic error in HepPolyhedronEllipsoid, memory corrupted!"
                                                   >> 1956                 << std::endl;
                                                   >> 1957     }
                                                   >> 1958   if (j < np1)
                                                   >> 1959     {
                                                   >> 1960       std::cerr << "Warning: logic error in HepPolyhedronEllipsoid."
                                                   >> 1961                 << std::endl;
                                                   >> 1962       np1= j;
                                                   >> 1963     }
                                                   >> 1964   zz[j] = 0.;
                                                   >> 1965   rr[j] = 0.;
2772                                                  1966 
                                                   >> 1967   
2773   //   R O T A T E    P O L Y L I N E S          1968   //   R O T A T E    P O L Y L I N E S
2774                                                  1969 
2775   RotateAroundZ(0, 0., twopi, np1, np2, zz, r << 1970   RotateAroundZ(0, 0.0, twopi, np1, 1, zz, rr, -1, 1); 
2776   SetReferences();                               1971   SetReferences();
2777                                                  1972 
2778   delete [] zz;                                  1973   delete [] zz;
2779   delete [] rr;                                  1974   delete [] rr;
2780                                                  1975 
2781   // rescale x and y vertex coordinates          1976   // 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   {                                              1977   {
2787     p->setX(p->x()*kx);                       << 1978     Point3D<double> * p= pV;
2788     p->setY(p->y()*ky);                       << 1979     for (int i=0; i<nvert; i++, p++) {
                                                   >> 1980       p->setX( p->x() * ax/cz );
                                                   >> 1981       p->setY( p->y() * by/cz );
                                                   >> 1982     }
2789   }                                              1983   }
2790 }                                                1984 }
2791                                                  1985 
2792 HepPolyhedronEllipsoid::~HepPolyhedronEllipso << 1986 HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() {}
2793                                                  1987 
2794 HepPolyhedronEllipticalCone::HepPolyhedronEll << 1988 HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone(double ax,
2795                                               << 1989                                                          double ay,
2796                                               << 1990                                                          double h,
2797                                               << 1991                                                          double zTopCut) 
2798 /********************************************    1992 /***********************************************************************
2799  *                                               1993  *                                                                     *
2800  * Name: HepPolyhedronEllipticalCone             1994  * Name: HepPolyhedronEllipticalCone                 Date:    8.9.2005 *
2801  * Author: D.Anninos                             1995  * Author: D.Anninos                                 Revised: 9.9.2005 *
2802  *                                               1996  *                                                                     *
2803  * Function: Constructor for EllipticalCone      1997  * Function: Constructor for EllipticalCone                            *
2804  *                                               1998  *                                                                     *
2805  * Input: ax, ay     - X & Y semi axes at z =    1999  * Input: ax, ay     - X & Y semi axes at z = 0                        *
2806  *        h          - height of full cone       2000  *        h          - height of full cone                             *
2807  *        zTopCut    - Top Cut in Z Axis         2001  *        zTopCut    - Top Cut in Z Axis                               *
2808  *                                               2002  *                                                                     *
2809  ********************************************    2003  ***********************************************************************/
2810 {                                                2004 {
2811   //   C H E C K   I N P U T   P A R A M E T     2005   //   C H E C K   I N P U T   P A R A M E T E R S
2812                                                  2006 
2813   G4int k = 0;                                << 2007   int k = 0;
2814   if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.)     2008   if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
2815                                                  2009 
2816   if (k != 0) {                                  2010   if (k != 0) {
2817     std::cerr << "HepPolyhedronCone: error in    2011     std::cerr << "HepPolyhedronCone: error in input parameters";
2818     std::cerr << std::endl;                      2012     std::cerr << std::endl;
2819     return;                                      2013     return;
2820   }                                              2014   }
2821                                               << 2015   
2822   //   P R E P A R E   T W O   P O L Y L I N     2016   //   P R E P A R E   T W O   P O L Y L I N E S
2823                                                  2017 
2824   zTopCut = (h >= zTopCut ? zTopCut : h);        2018   zTopCut = (h >= zTopCut ? zTopCut : h);
2825                                                  2019 
2826   G4double *zz, *rr;                          << 2020   double *zz, *rr;
2827   zz = new G4double[4];                       << 2021   zz = new double[4];
2828   rr = new G4double[4];                       << 2022   rr = new double[4];
2829   zz[0] =   zTopCut;                          << 2023   zz[0] =   zTopCut; 
2830   zz[1] =  -zTopCut;                          << 2024   zz[1] =  -zTopCut; 
2831   zz[2] =   zTopCut;                          << 2025   zz[2] =   zTopCut; 
2832   zz[3] =  -zTopCut;                          << 2026   zz[3] =  -zTopCut; 
2833   rr[0] =  (h-zTopCut);                          2027   rr[0] =  (h-zTopCut);
2834   rr[1] =  (h+zTopCut);                          2028   rr[1] =  (h+zTopCut);
2835   rr[2] =  0.;                                   2029   rr[2] =  0.;
2836   rr[3] =  0.;                                   2030   rr[3] =  0.;
2837                                                  2031 
2838   //   R O T A T E    P O L Y L I N E S          2032   //   R O T A T E    P O L Y L I N E S
2839                                                  2033 
2840   RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, - << 2034   RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, -1, 1); 
2841   SetReferences();                               2035   SetReferences();
2842                                                  2036 
2843   delete [] zz;                                  2037   delete [] zz;
2844   delete [] rr;                                  2038   delete [] rr;
2845                                                  2039 
2846   // rescale x and y vertex coordinates          2040   // rescale x and y vertex coordinates
2847  {                                               2041  {
2848    G4Point3D * p= pV;                         << 2042    Point3D<double> * p= pV;
2849    for (G4int i=0; i<nvert; i++, p++) {       << 2043    for (int i=0; i<nvert; i++, p++) {
2850      p->setX( p->x() * ax );                     2044      p->setX( p->x() * ax );
2851      p->setY( p->y() * ay );                     2045      p->setY( p->y() * ay );
2852    }                                             2046    }
2853  }                                               2047  }
2854 }                                                2048 }
2855                                                  2049 
2856 HepPolyhedronEllipticalCone::~HepPolyhedronEl << 2050 HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() {}
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;                               << 
2915   delete [] rr;                               << 
2916 }                                             << 
2917                                               << 
2918 HepPolyhedronHyperbolicMirror::~HepPolyhedron << 
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                                                  2051 
3102 HepPolyhedronBoxMesh::                        << 2052 int HepPolyhedron::fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
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 /********************************************    2053 /***********************************************************************
3351  *                                               2054  *                                                                     *
3352  * Name: HepPolyhedron::fNumberOfRotationStep    2055  * Name: HepPolyhedron::fNumberOfRotationSteps       Date:    24.06.97 *
3353  * Author: J.Allison (Manchester University)     2056  * Author: J.Allison (Manchester University)         Revised:          *
3354  *                                               2057  *                                                                     *
3355  * Function: Number of steps for whole circle    2058  * Function: Number of steps for whole circle                          *
3356  *                                               2059  *                                                                     *
3357  ********************************************    2060  ***********************************************************************/
3358                                                  2061 
3359 #include "BooleanProcessor.src"                  2062 #include "BooleanProcessor.src"
                                                   >> 2063 static BooleanProcessor processor;
3360                                                  2064 
3361 HepPolyhedron HepPolyhedron::add(const HepPol << 2065 HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const 
3362 /********************************************    2066 /***********************************************************************
3363  *                                               2067  *                                                                     *
3364  * Name: HepPolyhedron::add                      2068  * Name: HepPolyhedron::add                          Date:    19.03.00 *
3365  * Author: E.Chernyaev                           2069  * Author: E.Chernyaev                               Revised:          *
3366  *                                               2070  *                                                                     *
3367  * Function: Boolean "union" of two polyhedra    2071  * Function: Boolean "union" of two polyhedra                          *
3368  *                                               2072  *                                                                     *
3369  ********************************************    2073  ***********************************************************************/
3370 {                                                2074 {
3371   G4int ierr;                                 << 2075   return processor.execute(OP_UNION, *this, p);
3372   BooleanProcessor processor;                 << 
3373   return processor.execute(OP_UNION, *this, p << 
3374 }                                                2076 }
3375                                                  2077 
3376 HepPolyhedron HepPolyhedron::intersect(const  << 2078 HepPolyhedron HepPolyhedron::intersect(const HepPolyhedron & p) const 
3377 /********************************************    2079 /***********************************************************************
3378  *                                               2080  *                                                                     *
3379  * Name: HepPolyhedron::intersect                2081  * Name: HepPolyhedron::intersect                    Date:    19.03.00 *
3380  * Author: E.Chernyaev                           2082  * Author: E.Chernyaev                               Revised:          *
3381  *                                               2083  *                                                                     *
3382  * Function: Boolean "intersection" of two po    2084  * Function: Boolean "intersection" of two polyhedra                   *
3383  *                                               2085  *                                                                     *
3384  ********************************************    2086  ***********************************************************************/
3385 {                                                2087 {
3386   G4int ierr;                                 << 2088   return processor.execute(OP_INTERSECTION, *this, p);
3387   BooleanProcessor processor;                 << 
3388   return processor.execute(OP_INTERSECTION, * << 
3389 }                                                2089 }
3390                                                  2090 
3391 HepPolyhedron HepPolyhedron::subtract(const H << 2091 HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const 
3392 /********************************************    2092 /***********************************************************************
3393  *                                               2093  *                                                                     *
3394  * Name: HepPolyhedron::add                      2094  * Name: HepPolyhedron::add                          Date:    19.03.00 *
3395  * Author: E.Chernyaev                           2095  * Author: E.Chernyaev                               Revised:          *
3396  *                                               2096  *                                                                     *
3397  * Function: Boolean "subtraction" of "p" fro    2097  * Function: Boolean "subtraction" of "p" from "this"                  *
3398  *                                               2098  *                                                                     *
3399  ********************************************    2099  ***********************************************************************/
3400 {                                                2100 {
3401   G4int ierr;                                 << 2101   return processor.execute(OP_SUBTRACTION, *this, p);
3402   BooleanProcessor processor;                 << 
3403   return processor.execute(OP_SUBTRACTION, *t << 
3404 }                                                2102 }
3405                                                  2103 
3406 //NOTE : include the code of HepPolyhedronPro << 2104 bool HepPolyhedron::IsErrorBooleanProcess() const {
3407 //       since there is no BooleanProcessor.h << 2105   return processor.get_processor_error();
3408                                               << 2106 }
3409 #undef INTERSECTION                           << 
3410                                               << 
3411 #include "HepPolyhedronProcessor.src"         << 
3412                                                  2107