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 11.0.p4)


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