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 ]

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