ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HepPolyhedron.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HepPolyhedron.cc
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 //
27 //
28 //
29 //
30 // G4 Polyhedron library
31 //
32 // History:
33 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
34 //
35 // 30.09.96 E.Chernyaev
36 // - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
37 // - added GetNextUnitNormal, GetNextEdgeIndices, GetNextEdge
38 //
39 // 15.12.96 E.Chernyaev
40 // - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
41 // - rewritten G4PolyhedronCons;
42 // - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
43 //
44 // 01.06.97 E.Chernyaev
45 // - modified RotateAroundZ, added SetSideFacets
46 //
47 // 19.03.00 E.Chernyaev
48 // - implemented boolean operations (add, subtract, intersect) on polyhedra;
49 //
50 // 25.05.01 E.Chernyaev
51 // - added GetSurfaceArea() and GetVolume();
52 //
53 // 05.11.02 E.Chernyaev
54 // - added createTwistedTrap() and createPolyhedron();
55 //
56 // 20.06.05 G.Cosmo
57 // - added HepPolyhedronEllipsoid;
58 //
59 // 18.07.07 T.Nikitin
60 // - added HepParaboloid;
61 
62 #include "HepPolyhedron.h"
63 #include "G4PhysicalConstants.hh"
64 #include "G4Vector3D.hh"
65 
66 #include <cstdlib> // Required on some compilers for std::abs(int) ...
67 #include <cmath>
68 #include <cassert>
69 
70 using CLHEP::perMillion;
71 using CLHEP::deg;
72 using CLHEP::pi;
73 using CLHEP::twopi;
74 using CLHEP::nm;
75 const G4double spatialTolerance = 0.01*nm;
76 
77 /***********************************************************************
78  * *
79  * Name: HepPolyhedron operator << Date: 09.05.96 *
80  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
81  * *
82  * Function: Print contents of G4 polyhedron *
83  * *
84  ***********************************************************************/
85 std::ostream & operator<<(std::ostream & ostr, const G4Facet & facet) {
86  for (G4int k=0; k<4; k++) {
87  ostr << " " << facet.edge[k].v << "/" << facet.edge[k].f;
88  }
89  return ostr;
90 }
91 
92 std::ostream & operator<<(std::ostream & ostr, const HepPolyhedron & ph) {
93  ostr << std::endl;
94  ostr << "Nvertices=" << ph.nvert << ", Nfacets=" << ph.nface << std::endl;
95  G4int i;
96  for (i=1; i<=ph.nvert; i++) {
97  ostr << "xyz(" << i << ")="
98  << ph.pV[i].x() << ' ' << ph.pV[i].y() << ' ' << ph.pV[i].z()
99  << std::endl;
100  }
101  for (i=1; i<=ph.nface; i++) {
102  ostr << "face(" << i << ")=" << ph.pF[i] << std::endl;
103  }
104  return ostr;
105 }
106 
107 HepPolyhedron::HepPolyhedron(const HepPolyhedron &from)
108 /***********************************************************************
109  * *
110  * Name: HepPolyhedron copy constructor Date: 23.07.96 *
111  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
112  * *
113  ***********************************************************************/
114 : nvert(0), nface(0), pV(0), pF(0)
115 {
116  AllocateMemory(from.nvert, from.nface);
117  for (G4int i=1; i<=nvert; i++) pV[i] = from.pV[i];
118  for (G4int k=1; k<=nface; k++) pF[k] = from.pF[k];
119 }
120 
121 HepPolyhedron::HepPolyhedron(HepPolyhedron&& from)
122 /***********************************************************************
123  * *
124  * Name: HepPolyhedron move constructor Date: 04.11.2019 *
125  * Author: E.Tcherniaev (E.Cheryaev) Revised: *
126  * *
127  ***********************************************************************/
128 : nvert(0), nface(0), pV(nullptr), pF(nullptr)
129 {
130  nvert = from.nvert;
131  nface = from.nface;
132  pV = from.pV;
133  pF = from.pF;
134 
135  // Release the data from the source object
136  from.nvert = 0;
137  from.nface = 0;
138  from.pV = nullptr;
139  from.pF = nullptr;
140 }
141 
142 HepPolyhedron & HepPolyhedron::operator=(const HepPolyhedron &from)
143 /***********************************************************************
144  * *
145  * Name: HepPolyhedron operator = Date: 23.07.96 *
146  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
147  * *
148  * Function: Copy contents of one polyhedron to another *
149  * *
150  ***********************************************************************/
151 {
152  if (this != &from) {
153  AllocateMemory(from.nvert, from.nface);
154  for (G4int i=1; i<=nvert; i++) pV[i] = from.pV[i];
155  for (G4int k=1; k<=nface; k++) pF[k] = from.pF[k];
156  }
157  return *this;
158 }
159 
160 HepPolyhedron & HepPolyhedron::operator=(HepPolyhedron&& from)
161 /***********************************************************************
162  * *
163  * Name: HepPolyhedron move operator = Date: 04.11.2019 *
164  * Author: E.Tcherniaev (E.Chernyaev) Revised: *
165  * *
166  * Function: Move contents of one polyhedron to another *
167  * *
168  ***********************************************************************/
169 {
170  if (this != &from) {
171  delete [] pV;
172  delete [] pF;
173  nvert = from.nvert;
174  nface = from.nface;
175  pV = from.pV;
176  pF = from.pF;
177 
178  // Release the data from the source object
179  from.nvert = 0;
180  from.nface = 0;
181  from.pV = nullptr;
182  from.pF = nullptr;
183  }
184  return *this;
185 }
186 
187 G4int
188 HepPolyhedron::FindNeighbour(G4int iFace, G4int iNode, G4int iOrder) const
189 /***********************************************************************
190  * *
191  * Name: HepPolyhedron::FindNeighbour Date: 22.11.99 *
192  * Author: E.Chernyaev Revised: *
193  * *
194  * Function: Find neighbouring face *
195  * *
196  ***********************************************************************/
197 {
198  G4int i;
199  for (i=0; i<4; i++) {
200  if (iNode == std::abs(pF[iFace].edge[i].v)) break;
201  }
202  if (i == 4) {
203  std::cerr
204  << "HepPolyhedron::FindNeighbour: face " << iFace
205  << " has no node " << iNode
206  << std::endl;
207  return 0;
208  }
209  if (iOrder < 0) {
210  if ( --i < 0) i = 3;
211  if (pF[iFace].edge[i].v == 0) i = 2;
212  }
213  return (pF[iFace].edge[i].v > 0) ? 0 : pF[iFace].edge[i].f;
214 }
215 
216 G4Normal3D HepPolyhedron::FindNodeNormal(G4int iFace, G4int iNode) const
217 /***********************************************************************
218  * *
219  * Name: HepPolyhedron::FindNodeNormal Date: 22.11.99 *
220  * Author: E.Chernyaev Revised: *
221  * *
222  * Function: Find normal at given node *
223  * *
224  ***********************************************************************/
225 {
226  G4Normal3D normal = GetUnitNormal(iFace);
227  G4int k = iFace, iOrder = 1, n = 1;
228 
229  for(;;) {
230  k = FindNeighbour(k, iNode, iOrder);
231  if (k == iFace) break;
232  if (k > 0) {
233  n++;
234  normal += GetUnitNormal(k);
235  }else{
236  if (iOrder < 0) break;
237  k = iFace;
238  iOrder = -iOrder;
239  }
240  }
241  return normal.unit();
242 }
243 
244 G4int HepPolyhedron::GetNumberOfRotationSteps()
245 /***********************************************************************
246  * *
247  * Name: HepPolyhedron::GetNumberOfRotationSteps Date: 24.06.97 *
248  * Author: J.Allison (Manchester University) Revised: *
249  * *
250  * Function: Get number of steps for whole circle *
251  * *
252  ***********************************************************************/
253 {
254  return fNumberOfRotationSteps;
255 }
256 
257 void HepPolyhedron::SetNumberOfRotationSteps(G4int n)
258 /***********************************************************************
259  * *
260  * Name: HepPolyhedron::SetNumberOfRotationSteps Date: 24.06.97 *
261  * Author: J.Allison (Manchester University) Revised: *
262  * *
263  * Function: Set number of steps for whole circle *
264  * *
265  ***********************************************************************/
266 {
267  const G4int nMin = 3;
268  if (n < nMin) {
269  std::cerr
270  << "HepPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
271  << "number of steps per circle < " << nMin << "; forced to " << nMin
272  << std::endl;
273  fNumberOfRotationSteps = nMin;
274  }else{
275  fNumberOfRotationSteps = n;
276  }
277 }
278 
279 void HepPolyhedron::ResetNumberOfRotationSteps()
280 /***********************************************************************
281  * *
282  * Name: HepPolyhedron::GetNumberOfRotationSteps Date: 24.06.97 *
283  * Author: J.Allison (Manchester University) Revised: *
284  * *
285  * Function: Reset number of steps for whole circle to default value *
286  * *
287  ***********************************************************************/
288 {
289  fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
290 }
291 
292 void HepPolyhedron::AllocateMemory(G4int Nvert, G4int Nface)
293 /***********************************************************************
294  * *
295  * Name: HepPolyhedron::AllocateMemory Date: 19.06.96 *
296  * Author: E.Chernyaev (IHEP/Protvino) Revised: 05.11.02 *
297  * *
298  * Function: Allocate memory for GEANT4 polyhedron *
299  * *
300  * Input: Nvert - number of nodes *
301  * Nface - number of faces *
302  * *
303  ***********************************************************************/
304 {
305  if (nvert == Nvert && nface == Nface) return;
306  if (pV != 0) delete [] pV;
307  if (pF != 0) delete [] pF;
308  if (Nvert > 0 && Nface > 0) {
309  nvert = Nvert;
310  nface = Nface;
311  pV = new G4Point3D[nvert+1];
312  pF = new G4Facet[nface+1];
313  }else{
314  nvert = 0; nface = 0; pV = 0; pF = 0;
315  }
316 }
317 
318 void HepPolyhedron::CreatePrism()
319 /***********************************************************************
320  * *
321  * Name: HepPolyhedron::CreatePrism Date: 15.07.96 *
322  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
323  * *
324  * Function: Set facets for a prism *
325  * *
326  ***********************************************************************/
327 {
328  enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
329 
330  pF[1] = G4Facet(1,LEFT, 4,BACK, 3,RIGHT, 2,FRONT);
331  pF[2] = G4Facet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
332  pF[3] = G4Facet(8,TOP, 7,RIGHT, 3,BOTTOM, 4,LEFT);
333  pF[4] = G4Facet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
334  pF[5] = G4Facet(6,TOP, 5,LEFT, 1,BOTTOM, 2,RIGHT);
335  pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK, 8,LEFT);
336 }
337 
338 void HepPolyhedron::RotateEdge(G4int k1, G4int k2, G4double r1, G4double r2,
339  G4int v1, G4int v2, G4int vEdge,
340  G4bool ifWholeCircle, G4int nds, G4int &kface)
341 /***********************************************************************
342  * *
343  * Name: HepPolyhedron::RotateEdge Date: 05.12.96 *
344  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
345  * *
346  * Function: Create set of facets by rotation of an edge around Z-axis *
347  * *
348  * Input: k1, k2 - end vertices of the edge *
349  * r1, r2 - radiuses of the end vertices *
350  * v1, v2 - visibility of edges produced by rotation of the end *
351  * vertices *
352  * vEdge - visibility of the edge *
353  * ifWholeCircle - is true in case of whole circle rotation *
354  * nds - number of discrete steps *
355  * r[] - r-coordinates *
356  * kface - current free cell in the pF array *
357  * *
358  ***********************************************************************/
359 {
360  if (r1 == 0. && r2 == 0) return;
361 
362  G4int i;
363  G4int i1 = k1;
364  G4int i2 = k2;
365  G4int ii1 = ifWholeCircle ? i1 : i1+nds;
366  G4int ii2 = ifWholeCircle ? i2 : i2+nds;
367  G4int vv = ifWholeCircle ? vEdge : 1;
368 
369  if (nds == 1) {
370  if (r1 == 0.) {
371  pF[kface++] = G4Facet(i1,0, v2*i2,0, (i2+1),0);
372  }else if (r2 == 0.) {
373  pF[kface++] = G4Facet(i1,0, i2,0, v1*(i1+1),0);
374  }else{
375  pF[kface++] = G4Facet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),0);
376  }
377  }else{
378  if (r1 == 0.) {
379  pF[kface++] = G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0);
380  for (i2++,i=1; i<nds-1; i2++,i++) {
381  pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
382  }
383  pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
384  }else if (r2 == 0.) {
385  pF[kface++] = G4Facet(vv*i1,0, vEdge*i2,0, v1*(i1+1),0);
386  for (i1++,i=1; i<nds-1; i1++,i++) {
387  pF[kface++] = G4Facet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
388  }
389  pF[kface++] = G4Facet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
390  }else{
391  pF[kface++] = G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
392  for (i1++,i2++,i=1; i<nds-1; i1++,i2++,i++) {
393  pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
394  }
395  pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
396  }
397  }
398 }
399 
400 void HepPolyhedron::SetSideFacets(G4int ii[4], G4int vv[4],
401  G4int *kk, G4double *r,
402  G4double dphi, G4int nds, G4int &kface)
403 /***********************************************************************
404  * *
405  * Name: HepPolyhedron::SetSideFacets Date: 20.05.97 *
406  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
407  * *
408  * Function: Set side facets for the case of incomplete rotation *
409  * *
410  * Input: ii[4] - indices of original vertices *
411  * vv[4] - visibility of edges *
412  * kk[] - indices of nodes *
413  * r[] - radiuses *
414  * dphi - delta phi *
415  * nds - number of discrete steps *
416  * kface - current free cell in the pF array *
417  * *
418  ***********************************************************************/
419 {
420  G4int k1, k2, k3, k4;
421 
422  if (std::abs((G4double)(dphi-pi)) < perMillion) { // half a circle
423  for (G4int i=0; i<4; i++) {
424  k1 = ii[i];
425  k2 = ii[(i+1)%4];
426  if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
427  }
428  }
429 
430  if (ii[1] == ii[2]) {
431  k1 = kk[ii[0]];
432  k2 = kk[ii[2]];
433  k3 = kk[ii[3]];
434  pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
435  if (r[ii[0]] != 0.) k1 += nds;
436  if (r[ii[2]] != 0.) k2 += nds;
437  if (r[ii[3]] != 0.) k3 += nds;
438  pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
439  }else if (kk[ii[0]] == kk[ii[1]]) {
440  k1 = kk[ii[0]];
441  k2 = kk[ii[2]];
442  k3 = kk[ii[3]];
443  pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
444  if (r[ii[0]] != 0.) k1 += nds;
445  if (r[ii[2]] != 0.) k2 += nds;
446  if (r[ii[3]] != 0.) k3 += nds;
447  pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
448  }else if (kk[ii[2]] == kk[ii[3]]) {
449  k1 = kk[ii[0]];
450  k2 = kk[ii[1]];
451  k3 = kk[ii[2]];
452  pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
453  if (r[ii[0]] != 0.) k1 += nds;
454  if (r[ii[1]] != 0.) k2 += nds;
455  if (r[ii[2]] != 0.) k3 += nds;
456  pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
457  }else{
458  k1 = kk[ii[0]];
459  k2 = kk[ii[1]];
460  k3 = kk[ii[2]];
461  k4 = kk[ii[3]];
462  pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
463  if (r[ii[0]] != 0.) k1 += nds;
464  if (r[ii[1]] != 0.) k2 += nds;
465  if (r[ii[2]] != 0.) k3 += nds;
466  if (r[ii[3]] != 0.) k4 += nds;
467  pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
468  }
469 }
470 
471 void HepPolyhedron::RotateAroundZ(G4int nstep, G4double phi, G4double dphi,
472  G4int np1, G4int np2,
473  const G4double *z, G4double *r,
474  G4int nodeVis, G4int edgeVis)
475 /***********************************************************************
476  * *
477  * Name: HepPolyhedron::RotateAroundZ Date: 27.11.96 *
478  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
479  * *
480  * Function: Create HepPolyhedron for a solid produced by rotation of *
481  * two polylines around Z-axis *
482  * *
483  * Input: nstep - number of discrete steps, if 0 then default *
484  * phi - starting phi angle *
485  * dphi - delta phi *
486  * np1 - number of points in external polyline *
487  * (must be negative in case of closed polyline) *
488  * np2 - number of points in internal polyline (may be 1) *
489  * z[] - z-coordinates (+z >>> -z for both polylines) *
490  * r[] - r-coordinates *
491  * nodeVis - how to Draw edges joing consecutive positions of *
492  * node during rotation *
493  * edgeVis - how to Draw edges *
494  * *
495  ***********************************************************************/
496 {
497  static const G4double wholeCircle = twopi;
498 
499  // S E T R O T A T I O N P A R A M E T E R S
500 
501  G4bool ifWholeCircle = (std::abs(dphi-wholeCircle) < perMillion) ? true : false;
502  G4double delPhi = ifWholeCircle ? wholeCircle : dphi;
503  G4int nSphi = (nstep > 0) ?
504  nstep : G4int(delPhi*GetNumberOfRotationSteps()/wholeCircle+.5);
505  if (nSphi == 0) nSphi = 1;
506  G4int nVphi = ifWholeCircle ? nSphi : nSphi+1;
507  G4bool ifClosed = np1 > 0 ? false : true;
508 
509  // C O U N T V E R T E C E S
510 
511  G4int absNp1 = std::abs(np1);
512  G4int absNp2 = std::abs(np2);
513  G4int i1beg = 0;
514  G4int i1end = absNp1-1;
515  G4int i2beg = absNp1;
516  G4int i2end = absNp1+absNp2-1;
517  G4int i, j, k;
518 
519  for(i=i1beg; i<=i2end; i++) {
520  if (std::abs(r[i]) < spatialTolerance) r[i] = 0.;
521  }
522 
523  j = 0; // external nodes
524  for (i=i1beg; i<=i1end; i++) {
525  j += (r[i] == 0.) ? 1 : nVphi;
526  }
527 
528  G4bool ifSide1 = false; // internal nodes
529  G4bool ifSide2 = false;
530 
531  if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
532  j += (r[i2beg] == 0.) ? 1 : nVphi;
533  ifSide1 = true;
534  }
535 
536  for(i=i2beg+1; i<i2end; i++) {
537  j += (r[i] == 0.) ? 1 : nVphi;
538  }
539 
540  if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
541  if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
542  ifSide2 = true;
543  }
544 
545  // C O U N T F A C E S
546 
547  k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi; // external faces
548 
549  if (absNp2 > 1) { // internal faces
550  for(i=i2beg; i<i2end; i++) {
551  if (r[i] > 0. || r[i+1] > 0.) k += nSphi;
552  }
553 
554  if (ifClosed) {
555  if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
556  }
557  }
558 
559  if (!ifClosed) { // side faces
560  if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
561  if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) k += nSphi;
562  }
563 
564  if (!ifWholeCircle) { // phi_side faces
565  k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
566  }
567 
568  // A L L O C A T E M E M O R Y
569 
570  AllocateMemory(j, k);
571 
572  // G E N E R A T E V E R T E C E S
573 
574  G4int *kk;
575  kk = new G4int[absNp1+absNp2];
576 
577  k = 1;
578  for(i=i1beg; i<=i1end; i++) {
579  kk[i] = k;
580  if (r[i] == 0.)
581  { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
582  }
583 
584  i = i2beg;
585  if (ifSide1) {
586  kk[i] = k;
587  if (r[i] == 0.)
588  { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
589  }else{
590  kk[i] = kk[i1beg];
591  }
592 
593  for(i=i2beg+1; i<i2end; i++) {
594  kk[i] = k;
595  if (r[i] == 0.)
596  { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
597  }
598 
599  if (absNp2 > 1) {
600  i = i2end;
601  if (ifSide2) {
602  kk[i] = k;
603  if (r[i] == 0.) pV[k] = G4Point3D(0, 0, z[i]);
604  }else{
605  kk[i] = kk[i1end];
606  }
607  }
608 
609  G4double cosPhi, sinPhi;
610 
611  for(j=0; j<nVphi; j++) {
612  cosPhi = std::cos(phi+j*delPhi/nSphi);
613  sinPhi = std::sin(phi+j*delPhi/nSphi);
614  for(i=i1beg; i<=i2end; i++) {
615  if (r[i] != 0.)
616  pV[kk[i]+j] = G4Point3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
617  }
618  }
619 
620  // G E N E R A T E E X T E R N A L F A C E S
621 
622  G4int v1,v2;
623 
624  k = 1;
625  v2 = ifClosed ? nodeVis : 1;
626  for(i=i1beg; i<i1end; i++) {
627  v1 = v2;
628  if (!ifClosed && i == i1end-1) {
629  v2 = 1;
630  }else{
631  v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
632  }
633  RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
634  edgeVis, ifWholeCircle, nSphi, k);
635  }
636  if (ifClosed) {
637  RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
638  edgeVis, ifWholeCircle, nSphi, k);
639  }
640 
641  // G E N E R A T E I N T E R N A L F A C E S
642 
643  if (absNp2 > 1) {
644  v2 = ifClosed ? nodeVis : 1;
645  for(i=i2beg; i<i2end; i++) {
646  v1 = v2;
647  if (!ifClosed && i==i2end-1) {
648  v2 = 1;
649  }else{
650  v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
651  }
652  RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
653  edgeVis, ifWholeCircle, nSphi, k);
654  }
655  if (ifClosed) {
656  RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
657  edgeVis, ifWholeCircle, nSphi, k);
658  }
659  }
660 
661  // G E N E R A T E S I D E F A C E S
662 
663  if (!ifClosed) {
664  if (ifSide1) {
665  RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
666  -1, ifWholeCircle, nSphi, k);
667  }
668  if (ifSide2) {
669  RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
670  -1, ifWholeCircle, nSphi, k);
671  }
672  }
673 
674  // G E N E R A T E S I D E F A C E S for the case of incomplete circle
675 
676  if (!ifWholeCircle) {
677 
678  G4int ii[4], vv[4];
679 
680  if (ifClosed) {
681  for (i=i1beg; i<=i1end; i++) {
682  ii[0] = i;
683  ii[3] = (i == i1end) ? i1beg : i+1;
684  ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
685  ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
686  vv[0] = -1;
687  vv[1] = 1;
688  vv[2] = -1;
689  vv[3] = 1;
690  SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
691  }
692  }else{
693  for (i=i1beg; i<i1end; i++) {
694  ii[0] = i;
695  ii[3] = i+1;
696  ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
697  ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
698  vv[0] = (i == i1beg) ? 1 : -1;
699  vv[1] = 1;
700  vv[2] = (i == i1end-1) ? 1 : -1;
701  vv[3] = 1;
702  SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
703  }
704  }
705  }
706 
707  delete [] kk;
708 
709  if (k-1 != nface) {
710  std::cerr
711  << "Polyhedron::RotateAroundZ: number of generated faces ("
712  << k-1 << ") is not equal to the number of allocated faces ("
713  << nface << ")"
714  << std::endl;
715  }
716 }
717 
718 void HepPolyhedron::SetReferences()
719 /***********************************************************************
720  * *
721  * Name: HepPolyhedron::SetReferences Date: 04.12.96 *
722  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
723  * *
724  * Function: For each edge set reference to neighbouring facet *
725  * *
726  ***********************************************************************/
727 {
728  if (nface <= 0) return;
729 
730  struct edgeListMember {
731  edgeListMember *next;
732  G4int v2;
733  G4int iface;
734  G4int iedge;
735  } *edgeList, *freeList, **headList;
736 
737 
738  // A L L O C A T E A N D I N I T I A T E L I S T S
739 
740  edgeList = new edgeListMember[2*nface];
741  headList = new edgeListMember*[nvert];
742 
743  G4int i;
744  for (i=0; i<nvert; i++) {
745  headList[i] = 0;
746  }
747  freeList = edgeList;
748  for (i=0; i<2*nface-1; i++) {
749  edgeList[i].next = &edgeList[i+1];
750  }
751  edgeList[2*nface-1].next = 0;
752 
753  // L O O P A L O N G E D G E S
754 
755  G4int iface, iedge, nedge, i1, i2, k1, k2;
756  edgeListMember *prev, *cur;
757 
758  for(iface=1; iface<=nface; iface++) {
759  nedge = (pF[iface].edge[3].v == 0) ? 3 : 4;
760  for (iedge=0; iedge<nedge; iedge++) {
761  i1 = iedge;
762  i2 = (iedge < nedge-1) ? iedge+1 : 0;
763  i1 = std::abs(pF[iface].edge[i1].v);
764  i2 = std::abs(pF[iface].edge[i2].v);
765  k1 = (i1 < i2) ? i1 : i2; // k1 = ::min(i1,i2);
766  k2 = (i1 > i2) ? i1 : i2; // k2 = ::max(i1,i2);
767 
768  // check head of the List corresponding to k1
769  cur = headList[k1];
770  if (cur == 0) {
771  headList[k1] = freeList;
772  if (!freeList) {
773  std::cerr
774  << "Polyhedron::SetReferences: bad link "
775  << std::endl;
776  break;
777  }
778  freeList = freeList->next;
779  cur = headList[k1];
780  cur->next = 0;
781  cur->v2 = k2;
782  cur->iface = iface;
783  cur->iedge = iedge;
784  continue;
785  }
786 
787  if (cur->v2 == k2) {
788  headList[k1] = cur->next;
789  cur->next = freeList;
790  freeList = cur;
791  pF[iface].edge[iedge].f = cur->iface;
792  pF[cur->iface].edge[cur->iedge].f = iface;
793  i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
794  i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
795  if (i1 != i2) {
796  std::cerr
797  << "Polyhedron::SetReferences: different edge visibility "
798  << iface << "/" << iedge << "/"
799  << pF[iface].edge[iedge].v << " and "
800  << cur->iface << "/" << cur->iedge << "/"
801  << pF[cur->iface].edge[cur->iedge].v
802  << std::endl;
803  }
804  continue;
805  }
806 
807  // check List itself
808  for (;;) {
809  prev = cur;
810  cur = prev->next;
811  if (cur == 0) {
812  prev->next = freeList;
813  if (!freeList) {
814  std::cerr
815  << "Polyhedron::SetReferences: bad link "
816  << std::endl;
817  break;
818  }
819  freeList = freeList->next;
820  cur = prev->next;
821  cur->next = 0;
822  cur->v2 = k2;
823  cur->iface = iface;
824  cur->iedge = iedge;
825  break;
826  }
827 
828  if (cur->v2 == k2) {
829  prev->next = cur->next;
830  cur->next = freeList;
831  freeList = cur;
832  pF[iface].edge[iedge].f = cur->iface;
833  pF[cur->iface].edge[cur->iedge].f = iface;
834  i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
835  i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
836  if (i1 != i2) {
837  std::cerr
838  << "Polyhedron::SetReferences: different edge visibility "
839  << iface << "/" << iedge << "/"
840  << pF[iface].edge[iedge].v << " and "
841  << cur->iface << "/" << cur->iedge << "/"
842  << pF[cur->iface].edge[cur->iedge].v
843  << std::endl;
844  }
845  break;
846  }
847  }
848  }
849  }
850 
851  // 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
852 
853  for (i=0; i<nvert; i++) {
854  if (headList[i] != 0) {
855  std::cerr
856  << "Polyhedron::SetReferences: List " << i << " is not empty"
857  << std::endl;
858  }
859  }
860 
861  // F R E E M E M O R Y
862 
863  delete [] edgeList;
864  delete [] headList;
865 }
866 
867 void HepPolyhedron::InvertFacets()
868 /***********************************************************************
869  * *
870  * Name: HepPolyhedron::InvertFacets Date: 01.12.99 *
871  * Author: E.Chernyaev Revised: *
872  * *
873  * Function: Invert the order of the nodes in the facets *
874  * *
875  ***********************************************************************/
876 {
877  if (nface <= 0) return;
878  G4int i, k, nnode, v[4],f[4];
879  for (i=1; i<=nface; i++) {
880  nnode = (pF[i].edge[3].v == 0) ? 3 : 4;
881  for (k=0; k<nnode; k++) {
882  v[k] = (k+1 == nnode) ? pF[i].edge[0].v : pF[i].edge[k+1].v;
883  if (v[k] * pF[i].edge[k].v < 0) v[k] = -v[k];
884  f[k] = pF[i].edge[k].f;
885  }
886  for (k=0; k<nnode; k++) {
887  pF[i].edge[nnode-1-k].v = v[k];
888  pF[i].edge[nnode-1-k].f = f[k];
889  }
890  }
891 }
892 
893 HepPolyhedron & HepPolyhedron::Transform(const G4Transform3D &t)
894 /***********************************************************************
895  * *
896  * Name: HepPolyhedron::Transform Date: 01.12.99 *
897  * Author: E.Chernyaev Revised: *
898  * *
899  * Function: Make transformation of the polyhedron *
900  * *
901  ***********************************************************************/
902 {
903  if (nvert > 0) {
904  for (G4int i=1; i<=nvert; i++) { pV[i] = t * pV[i]; }
905 
906  // C H E C K D E T E R M I N A N T A N D
907  // 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
908 
909  G4Vector3D d = t * G4Vector3D(0,0,0);
910  G4Vector3D x = t * G4Vector3D(1,0,0) - d;
911  G4Vector3D y = t * G4Vector3D(0,1,0) - d;
912  G4Vector3D z = t * G4Vector3D(0,0,1) - d;
913  if ((x.cross(y))*z < 0) InvertFacets();
914  }
915  return *this;
916 }
917 
918 G4bool HepPolyhedron::GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
919 /***********************************************************************
920  * *
921  * Name: HepPolyhedron::GetNextVertexIndex Date: 03.09.96 *
922  * Author: Yasuhide Sawada Revised: *
923  * *
924  * Function: *
925  * *
926  ***********************************************************************/
927 {
928  static G4ThreadLocal G4int iFace = 1;
929  static G4ThreadLocal G4int iQVertex = 0;
930  G4int vIndex = pF[iFace].edge[iQVertex].v;
931 
932  edgeFlag = (vIndex > 0) ? 1 : 0;
933  index = std::abs(vIndex);
934 
935  if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
936  iQVertex = 0;
937  if (++iFace > nface) iFace = 1;
938  return false; // Last Edge
939  }else{
940  ++iQVertex;
941  return true; // not Last Edge
942  }
943 }
944 
945 G4Point3D HepPolyhedron::GetVertex(G4int index) const
946 /***********************************************************************
947  * *
948  * Name: HepPolyhedron::GetVertex Date: 03.09.96 *
949  * Author: Yasuhide Sawada Revised: 17.11.99 *
950  * *
951  * Function: Get vertex of the index. *
952  * *
953  ***********************************************************************/
954 {
955  if (index <= 0 || index > nvert) {
956  std::cerr
957  << "HepPolyhedron::GetVertex: irrelevant index " << index
958  << std::endl;
959  return G4Point3D();
960  }
961  return pV[index];
962 }
963 
964 G4bool
965 HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag) const
966 /***********************************************************************
967  * *
968  * Name: HepPolyhedron::GetNextVertex Date: 22.07.96 *
969  * Author: John Allison Revised: *
970  * *
971  * Function: Get vertices of the quadrilaterals in order for each *
972  * face in face order. Returns false when finished each *
973  * face. *
974  * *
975  ***********************************************************************/
976 {
977  G4int index;
978  G4bool rep = GetNextVertexIndex(index, edgeFlag);
979  vertex = pV[index];
980  return rep;
981 }
982 
983 G4bool HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag,
984  G4Normal3D &normal) const
985 /***********************************************************************
986  * *
987  * Name: HepPolyhedron::GetNextVertex Date: 26.11.99 *
988  * Author: E.Chernyaev Revised: *
989  * *
990  * Function: Get vertices with normals of the quadrilaterals in order *
991  * for each face in face order. *
992  * Returns false when finished each face. *
993  * *
994  ***********************************************************************/
995 {
996  static G4ThreadLocal G4int iFace = 1;
997  static G4ThreadLocal G4int iNode = 0;
998 
999  if (nface == 0) return false; // empty polyhedron
1000 
1001  G4int k = pF[iFace].edge[iNode].v;
1002  if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
1003  vertex = pV[k];
1004  normal = FindNodeNormal(iFace,k);
1005  if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
1006  iNode = 0;
1007  if (++iFace > nface) iFace = 1;
1008  return false; // last node
1009  }else{
1010  ++iNode;
1011  return true; // not last node
1012  }
1013 }
1014 
1015 G4bool HepPolyhedron::GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag,
1016  G4int &iface1, G4int &iface2) const
1017 /***********************************************************************
1018  * *
1019  * Name: HepPolyhedron::GetNextEdgeIndices Date: 30.09.96 *
1020  * Author: E.Chernyaev Revised: 17.11.99 *
1021  * *
1022  * Function: Get indices of the next edge together with indices of *
1023  * of the faces which share the edge. *
1024  * Returns false when the last edge. *
1025  * *
1026  ***********************************************************************/
1027 {
1028  static G4ThreadLocal G4int iFace = 1;
1029  static G4ThreadLocal G4int iQVertex = 0;
1030  static G4ThreadLocal G4int iOrder = 1;
1031  G4int k1, k2, kflag, kface1, kface2;
1032 
1033  if (iFace == 1 && iQVertex == 0) {
1034  k2 = pF[nface].edge[0].v;
1035  k1 = pF[nface].edge[3].v;
1036  if (k1 == 0) k1 = pF[nface].edge[2].v;
1037  if (std::abs(k1) > std::abs(k2)) iOrder = -1;
1038  }
1039 
1040  do {
1041  k1 = pF[iFace].edge[iQVertex].v;
1042  kflag = k1;
1043  k1 = std::abs(k1);
1044  kface1 = iFace;
1045  kface2 = pF[iFace].edge[iQVertex].f;
1046  if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1047  iQVertex = 0;
1048  k2 = std::abs(pF[iFace].edge[iQVertex].v);
1049  iFace++;
1050  }else{
1051  iQVertex++;
1052  k2 = std::abs(pF[iFace].edge[iQVertex].v);
1053  }
1054  } while (iOrder*k1 > iOrder*k2);
1055 
1056  i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1057  iface1 = kface1; iface2 = kface2;
1058 
1059  if (iFace > nface) {
1060  iFace = 1; iOrder = 1;
1061  return false;
1062  }else{
1063  return true;
1064  }
1065 }
1066 
1067 G4bool
1068 HepPolyhedron::GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag) const
1069 /***********************************************************************
1070  * *
1071  * Name: HepPolyhedron::GetNextEdgeIndices Date: 17.11.99 *
1072  * Author: E.Chernyaev Revised: *
1073  * *
1074  * Function: Get indices of the next edge. *
1075  * Returns false when the last edge. *
1076  * *
1077  ***********************************************************************/
1078 {
1079  G4int kface1, kface2;
1080  return GetNextEdgeIndices(i1, i2, edgeFlag, kface1, kface2);
1081 }
1082 
1083 G4bool
1084 HepPolyhedron::GetNextEdge(G4Point3D &p1,
1085  G4Point3D &p2,
1086  G4int &edgeFlag) const
1087 /***********************************************************************
1088  * *
1089  * Name: HepPolyhedron::GetNextEdge Date: 30.09.96 *
1090  * Author: E.Chernyaev Revised: *
1091  * *
1092  * Function: Get next edge. *
1093  * Returns false when the last edge. *
1094  * *
1095  ***********************************************************************/
1096 {
1097  G4int i1,i2;
1098  G4bool rep = GetNextEdgeIndices(i1,i2,edgeFlag);
1099  p1 = pV[i1];
1100  p2 = pV[i2];
1101  return rep;
1102 }
1103 
1104 G4bool
1105 HepPolyhedron::GetNextEdge(G4Point3D &p1, G4Point3D &p2,
1106  G4int &edgeFlag, G4int &iface1, G4int &iface2) const
1107 /***********************************************************************
1108  * *
1109  * Name: HepPolyhedron::GetNextEdge Date: 17.11.99 *
1110  * Author: E.Chernyaev Revised: *
1111  * *
1112  * Function: Get next edge with indices of the faces which share *
1113  * the edge. *
1114  * Returns false when the last edge. *
1115  * *
1116  ***********************************************************************/
1117 {
1118  G4int i1,i2;
1119  G4bool rep = GetNextEdgeIndices(i1,i2,edgeFlag,iface1,iface2);
1120  p1 = pV[i1];
1121  p2 = pV[i2];
1122  return rep;
1123 }
1124 
1125 void HepPolyhedron::GetFacet(G4int iFace, G4int &n, G4int *iNodes,
1126  G4int *edgeFlags, G4int *iFaces) const
1127 /***********************************************************************
1128  * *
1129  * Name: HepPolyhedron::GetFacet Date: 15.12.99 *
1130  * Author: E.Chernyaev Revised: *
1131  * *
1132  * Function: Get face by index *
1133  * *
1134  ***********************************************************************/
1135 {
1136  if (iFace < 1 || iFace > nface) {
1137  std::cerr
1138  << "HepPolyhedron::GetFacet: irrelevant index " << iFace
1139  << std::endl;
1140  n = 0;
1141  }else{
1142  G4int i, k;
1143  for (i=0; i<4; i++) {
1144  k = pF[iFace].edge[i].v;
1145  if (k == 0) break;
1146  if (iFaces != 0) iFaces[i] = pF[iFace].edge[i].f;
1147  if (k > 0) {
1148  iNodes[i] = k;
1149  if (edgeFlags != 0) edgeFlags[i] = 1;
1150  }else{
1151  iNodes[i] = -k;
1152  if (edgeFlags != 0) edgeFlags[i] = -1;
1153  }
1154  }
1155  n = i;
1156  }
1157 }
1158 
1159 void HepPolyhedron::GetFacet(G4int index, G4int &n, G4Point3D *nodes,
1160  G4int *edgeFlags, G4Normal3D *normals) const
1161 /***********************************************************************
1162  * *
1163  * Name: HepPolyhedron::GetFacet Date: 17.11.99 *
1164  * Author: E.Chernyaev Revised: *
1165  * *
1166  * Function: Get face by index *
1167  * *
1168  ***********************************************************************/
1169 {
1170  G4int iNodes[4];
1171  GetFacet(index, n, iNodes, edgeFlags);
1172  if (n != 0) {
1173  for (G4int i=0; i<n; i++) {
1174  nodes[i] = pV[iNodes[i]];
1175  if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1176  }
1177  }
1178 }
1179 
1180 G4bool
1181 HepPolyhedron::GetNextFacet(G4int &n, G4Point3D *nodes,
1182  G4int *edgeFlags, G4Normal3D *normals) const
1183 /***********************************************************************
1184  * *
1185  * Name: HepPolyhedron::GetNextFacet Date: 19.11.99 *
1186  * Author: E.Chernyaev Revised: *
1187  * *
1188  * Function: Get next face with normals of unit length at the nodes. *
1189  * Returns false when finished all faces. *
1190  * *
1191  ***********************************************************************/
1192 {
1193  static G4ThreadLocal G4int iFace = 1;
1194 
1195  if (edgeFlags == 0) {
1196  GetFacet(iFace, n, nodes);
1197  }else if (normals == 0) {
1198  GetFacet(iFace, n, nodes, edgeFlags);
1199  }else{
1200  GetFacet(iFace, n, nodes, edgeFlags, normals);
1201  }
1202 
1203  if (++iFace > nface) {
1204  iFace = 1;
1205  return false;
1206  }else{
1207  return true;
1208  }
1209 }
1210 
1211 G4Normal3D HepPolyhedron::GetNormal(G4int iFace) const
1212 /***********************************************************************
1213  * *
1214  * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1215  * Author: E.Chernyaev Revised: *
1216  * *
1217  * Function: Get normal of the face given by index *
1218  * *
1219  ***********************************************************************/
1220 {
1221  if (iFace < 1 || iFace > nface) {
1222  std::cerr
1223  << "HepPolyhedron::GetNormal: irrelevant index " << iFace
1224  << std::endl;
1225  return G4Normal3D();
1226  }
1227 
1228  G4int i0 = std::abs(pF[iFace].edge[0].v);
1229  G4int i1 = std::abs(pF[iFace].edge[1].v);
1230  G4int i2 = std::abs(pF[iFace].edge[2].v);
1231  G4int i3 = std::abs(pF[iFace].edge[3].v);
1232  if (i3 == 0) i3 = i0;
1233  return (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1]);
1234 }
1235 
1236 G4Normal3D HepPolyhedron::GetUnitNormal(G4int iFace) const
1237 /***********************************************************************
1238  * *
1239  * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1240  * Author: E.Chernyaev Revised: *
1241  * *
1242  * Function: Get unit normal of the face given by index *
1243  * *
1244  ***********************************************************************/
1245 {
1246  if (iFace < 1 || iFace > nface) {
1247  std::cerr
1248  << "HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
1249  << std::endl;
1250  return G4Normal3D();
1251  }
1252 
1253  G4int i0 = std::abs(pF[iFace].edge[0].v);
1254  G4int i1 = std::abs(pF[iFace].edge[1].v);
1255  G4int i2 = std::abs(pF[iFace].edge[2].v);
1256  G4int i3 = std::abs(pF[iFace].edge[3].v);
1257  if (i3 == 0) i3 = i0;
1258  return ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).unit();
1259 }
1260 
1261 G4bool HepPolyhedron::GetNextNormal(G4Normal3D &normal) const
1262 /***********************************************************************
1263  * *
1264  * Name: HepPolyhedron::GetNextNormal Date: 22.07.96 *
1265  * Author: John Allison Revised: 19.11.99 *
1266  * *
1267  * Function: Get normals of each face in face order. Returns false *
1268  * when finished all faces. *
1269  * *
1270  ***********************************************************************/
1271 {
1272  static G4ThreadLocal G4int iFace = 1;
1273  normal = GetNormal(iFace);
1274  if (++iFace > nface) {
1275  iFace = 1;
1276  return false;
1277  }else{
1278  return true;
1279  }
1280 }
1281 
1282 G4bool HepPolyhedron::GetNextUnitNormal(G4Normal3D &normal) const
1283 /***********************************************************************
1284  * *
1285  * Name: HepPolyhedron::GetNextUnitNormal Date: 16.09.96 *
1286  * Author: E.Chernyaev Revised: *
1287  * *
1288  * Function: Get normals of unit length of each face in face order. *
1289  * Returns false when finished all faces. *
1290  * *
1291  ***********************************************************************/
1292 {
1293  G4bool rep = GetNextNormal(normal);
1294  normal = normal.unit();
1295  return rep;
1296 }
1297 
1298 G4double HepPolyhedron::GetSurfaceArea() const
1299 /***********************************************************************
1300  * *
1301  * Name: HepPolyhedron::GetSurfaceArea Date: 25.05.01 *
1302  * Author: E.Chernyaev Revised: *
1303  * *
1304  * Function: Returns area of the surface of the polyhedron. *
1305  * *
1306  ***********************************************************************/
1307 {
1308  G4double srf = 0.;
1309  for (G4int iFace=1; iFace<=nface; iFace++) {
1310  G4int i0 = std::abs(pF[iFace].edge[0].v);
1311  G4int i1 = std::abs(pF[iFace].edge[1].v);
1312  G4int i2 = std::abs(pF[iFace].edge[2].v);
1313  G4int i3 = std::abs(pF[iFace].edge[3].v);
1314  if (i3 == 0) i3 = i0;
1315  srf += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).mag();
1316  }
1317  return srf/2.;
1318 }
1319 
1320 G4double HepPolyhedron::GetVolume() const
1321 /***********************************************************************
1322  * *
1323  * Name: HepPolyhedron::GetVolume Date: 25.05.01 *
1324  * Author: E.Chernyaev Revised: *
1325  * *
1326  * Function: Returns volume of the polyhedron. *
1327  * *
1328  ***********************************************************************/
1329 {
1330  G4double v = 0.;
1331  for (G4int iFace=1; iFace<=nface; iFace++) {
1332  G4int i0 = std::abs(pF[iFace].edge[0].v);
1333  G4int i1 = std::abs(pF[iFace].edge[1].v);
1334  G4int i2 = std::abs(pF[iFace].edge[2].v);
1335  G4int i3 = std::abs(pF[iFace].edge[3].v);
1336  G4Point3D pt;
1337  if (i3 == 0) {
1338  i3 = i0;
1339  pt = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);
1340  }else{
1341  pt = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0.25;
1342  }
1343  v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).dot(pt);
1344  }
1345  return v/6.;
1346 }
1347 
1348 G4int
1349 HepPolyhedron::createTwistedTrap(G4double Dz,
1350  const G4double xy1[][2],
1351  const G4double xy2[][2])
1352 /***********************************************************************
1353  * *
1354  * Name: createTwistedTrap Date: 05.11.02 *
1355  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1356  * *
1357  * Function: Creates polyhedron for twisted trapezoid *
1358  * *
1359  * Input: Dz - half-length along Z 8----7 *
1360  * xy1[2,4] - quadrilateral at Z=-Dz 5----6 ! *
1361  * xy2[2,4] - quadrilateral at Z=+Dz ! 4-!--3 *
1362  * 1----2 *
1363  * *
1364  ***********************************************************************/
1365 {
1366  AllocateMemory(12,18);
1367 
1368  pV[ 1] = G4Point3D(xy1[0][0],xy1[0][1],-Dz);
1369  pV[ 2] = G4Point3D(xy1[1][0],xy1[1][1],-Dz);
1370  pV[ 3] = G4Point3D(xy1[2][0],xy1[2][1],-Dz);
1371  pV[ 4] = G4Point3D(xy1[3][0],xy1[3][1],-Dz);
1372 
1373  pV[ 5] = G4Point3D(xy2[0][0],xy2[0][1], Dz);
1374  pV[ 6] = G4Point3D(xy2[1][0],xy2[1][1], Dz);
1375  pV[ 7] = G4Point3D(xy2[2][0],xy2[2][1], Dz);
1376  pV[ 8] = G4Point3D(xy2[3][0],xy2[3][1], Dz);
1377 
1378  pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;
1379  pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;
1380  pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;
1381  pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;
1382 
1383  enum {DUMMY, BOTTOM,
1384  LEFT_BOTTOM, LEFT_FRONT, LEFT_TOP, LEFT_BACK,
1385  BACK_BOTTOM, BACK_LEFT, BACK_TOP, BACK_RIGHT,
1386  RIGHT_BOTTOM, RIGHT_BACK, RIGHT_TOP, RIGHT_FRONT,
1387  FRONT_BOTTOM, FRONT_RIGHT, FRONT_TOP, FRONT_LEFT,
1388  TOP};
1389 
1390  pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
1391 
1392  pF[ 2]=G4Facet(4,BOTTOM, -1,LEFT_FRONT, -12,LEFT_BACK, 0,0);
1393  pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP, -12,LEFT_BOTTOM, 0,0);
1394  pF[ 4]=G4Facet(5,TOP, -8,LEFT_BACK, -12,LEFT_FRONT, 0,0);
1395  pF[ 5]=G4Facet(8,BACK_LEFT, -4,LEFT_BOTTOM, -12,LEFT_TOP, 0,0);
1396 
1397  pF[ 6]=G4Facet(3,BOTTOM, -4,BACK_LEFT, -11,BACK_RIGHT, 0,0);
1398  pF[ 7]=G4Facet(4,LEFT_BACK, -8,BACK_TOP, -11,BACK_BOTTOM, 0,0);
1399  pF[ 8]=G4Facet(8,TOP, -7,BACK_RIGHT, -11,BACK_LEFT, 0,0);
1400  pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP, 0,0);
1401 
1402  pF[10]=G4Facet(2,BOTTOM, -3,RIGHT_BACK, -10,RIGHT_FRONT, 0,0);
1403  pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP, -10,RIGHT_BOTTOM, 0,0);
1404  pF[12]=G4Facet(7,TOP, -6,RIGHT_FRONT, -10,RIGHT_BACK, 0,0);
1405  pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP, 0,0);
1406 
1407  pF[14]=G4Facet(1,BOTTOM, -2,FRONT_RIGHT, -9,FRONT_LEFT, 0,0);
1408  pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP, -9,FRONT_BOTTOM, 0,0);
1409  pF[16]=G4Facet(6,TOP, -5,FRONT_LEFT, -9,FRONT_RIGHT, 0,0);
1410  pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP, 0,0);
1411 
1412  pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
1413 
1414  return 0;
1415 }
1416 
1417 G4int
1418 HepPolyhedron::createPolyhedron(G4int Nnodes, G4int Nfaces,
1419  const G4double xyz[][3],
1420  const G4int faces[][4])
1421 /***********************************************************************
1422  * *
1423  * Name: createPolyhedron Date: 05.11.02 *
1424  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1425  * *
1426  * Function: Creates user defined polyhedron *
1427  * *
1428  * Input: Nnodes - number of nodes *
1429  * Nfaces - number of faces *
1430  * nodes[][3] - node coordinates *
1431  * faces[][4] - faces *
1432  * *
1433  ***********************************************************************/
1434 {
1435  AllocateMemory(Nnodes, Nfaces);
1436  if (nvert == 0) return 1;
1437 
1438  for (G4int i=0; i<Nnodes; i++) {
1439  pV[i+1] = G4Point3D(xyz[i][0], xyz[i][1], xyz[i][2]);
1440  }
1441  for (G4int k=0; k<Nfaces; k++) {
1442  pF[k+1] = G4Facet(faces[k][0],0,faces[k][1],0,faces[k][2],0,faces[k][3],0);
1443  }
1444  SetReferences();
1445  return 0;
1446 }
1447 
1448 HepPolyhedronTrd2::HepPolyhedronTrd2(G4double Dx1, G4double Dx2,
1449  G4double Dy1, G4double Dy2,
1450  G4double Dz)
1451 /***********************************************************************
1452  * *
1453  * Name: HepPolyhedronTrd2 Date: 22.07.96 *
1454  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1455  * *
1456  * Function: Create GEANT4 TRD2-trapezoid *
1457  * *
1458  * Input: Dx1 - half-length along X at -Dz 8----7 *
1459  * Dx2 - half-length along X ay +Dz 5----6 ! *
1460  * Dy1 - half-length along Y ay -Dz ! 4-!--3 *
1461  * Dy2 - half-length along Y ay +Dz 1----2 *
1462  * Dz - half-length along Z *
1463  * *
1464  ***********************************************************************/
1465 {
1466  AllocateMemory(8,6);
1467 
1468  pV[1] = G4Point3D(-Dx1,-Dy1,-Dz);
1469  pV[2] = G4Point3D( Dx1,-Dy1,-Dz);
1470  pV[3] = G4Point3D( Dx1, Dy1,-Dz);
1471  pV[4] = G4Point3D(-Dx1, Dy1,-Dz);
1472  pV[5] = G4Point3D(-Dx2,-Dy2, Dz);
1473  pV[6] = G4Point3D( Dx2,-Dy2, Dz);
1474  pV[7] = G4Point3D( Dx2, Dy2, Dz);
1475  pV[8] = G4Point3D(-Dx2, Dy2, Dz);
1476 
1477  CreatePrism();
1478 }
1479 
1480 HepPolyhedronTrd2::~HepPolyhedronTrd2() {}
1481 
1482 HepPolyhedronTrd1::HepPolyhedronTrd1(G4double Dx1, G4double Dx2,
1483  G4double Dy, G4double Dz)
1484  : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1485 
1486 HepPolyhedronTrd1::~HepPolyhedronTrd1() {}
1487 
1488 HepPolyhedronBox::HepPolyhedronBox(G4double Dx, G4double Dy, G4double Dz)
1489  : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1490 
1491 HepPolyhedronBox::~HepPolyhedronBox() {}
1492 
1493 HepPolyhedronTrap::HepPolyhedronTrap(G4double Dz,
1494  G4double Theta,
1495  G4double Phi,
1496  G4double Dy1,
1497  G4double Dx1,
1498  G4double Dx2,
1499  G4double Alp1,
1500  G4double Dy2,
1501  G4double Dx3,
1502  G4double Dx4,
1503  G4double Alp2)
1504 /***********************************************************************
1505  * *
1506  * Name: HepPolyhedronTrap Date: 20.11.96 *
1507  * Author: E.Chernyaev Revised: *
1508  * *
1509  * Function: Create GEANT4 TRAP-trapezoid *
1510  * *
1511  * Input: DZ - half-length in Z *
1512  * Theta,Phi - polar angles of the line joining centres of the *
1513  * faces at Z=-Dz and Z=+Dz *
1514  * Dy1 - half-length in Y of the face at Z=-Dz *
1515  * Dx1 - half-length in X of low edge of the face at Z=-Dz *
1516  * Dx2 - half-length in X of top edge of the face at Z=-Dz *
1517  * Alp1 - angle between Y-axis and the median joining top and *
1518  * low edges of the face at Z=-Dz *
1519  * Dy2 - half-length in Y of the face at Z=+Dz *
1520  * Dx3 - half-length in X of low edge of the face at Z=+Dz *
1521  * Dx4 - half-length in X of top edge of the face at Z=+Dz *
1522  * Alp2 - angle between Y-axis and the median joining top and *
1523  * low edges of the face at Z=+Dz *
1524  * *
1525  ***********************************************************************/
1526 {
1527  G4double DzTthetaCphi = Dz*std::tan(Theta)*std::cos(Phi);
1528  G4double DzTthetaSphi = Dz*std::tan(Theta)*std::sin(Phi);
1529  G4double Dy1Talp1 = Dy1*std::tan(Alp1);
1530  G4double Dy2Talp2 = Dy2*std::tan(Alp2);
1531 
1532  AllocateMemory(8,6);
1533 
1534  pV[1] = G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
1535  pV[2] = G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
1536  pV[3] = G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
1537  pV[4] = G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
1538  pV[5] = G4Point3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
1539  pV[6] = G4Point3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
1540  pV[7] = G4Point3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
1541  pV[8] = G4Point3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
1542 
1543  CreatePrism();
1544 }
1545 
1546 HepPolyhedronTrap::~HepPolyhedronTrap() {}
1547 
1548 HepPolyhedronPara::HepPolyhedronPara(G4double Dx, G4double Dy, G4double Dz,
1549  G4double Alpha, G4double Theta,
1550  G4double Phi)
1551  : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
1552 
1553 HepPolyhedronPara::~HepPolyhedronPara() {}
1554 
1555 HepPolyhedronParaboloid::HepPolyhedronParaboloid(G4double r1,
1556  G4double r2,
1557  G4double dz,
1558  G4double sPhi,
1559  G4double dPhi)
1560 /***********************************************************************
1561  * *
1562  * Name: HepPolyhedronParaboloid Date: 28.06.07 *
1563  * Author: L.Lindroos, T.Nikitina (CERN), July 2007 Revised: 28.06.07 *
1564  * *
1565  * Function: Constructor for paraboloid *
1566  * *
1567  * Input: r1 - inside and outside radiuses at -Dz *
1568  * r2 - inside and outside radiuses at +Dz *
1569  * dz - half length in Z *
1570  * sPhi - starting angle of the segment *
1571  * dPhi - segment range *
1572  * *
1573  ***********************************************************************/
1574 {
1575  static const G4double wholeCircle=twopi;
1576 
1577  // C H E C K I N P U T P A R A M E T E R S
1578 
1579  G4int k = 0;
1580  if (r1 < 0. || r2 <= 0.) k = 1;
1581 
1582  if (dz <= 0.) k += 2;
1583 
1584  G4double phi1, phi2, dphi;
1585 
1586  if(dPhi < 0.)
1587  {
1588  phi2 = sPhi; phi1 = phi2 + dPhi;
1589  }
1590  else if(dPhi == 0.)
1591  {
1592  phi1 = sPhi; phi2 = phi1 + wholeCircle;
1593  }
1594  else
1595  {
1596  phi1 = sPhi; phi2 = phi1 + dPhi;
1597  }
1598  dphi = phi2 - phi1;
1599 
1600  if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1601  if (dphi > wholeCircle) k += 4;
1602 
1603  if (k != 0) {
1604  std::cerr << "HepPolyhedronParaboloid: error in input parameters";
1605  if ((k & 1) != 0) std::cerr << " (radiuses)";
1606  if ((k & 2) != 0) std::cerr << " (half-length)";
1607  if ((k & 4) != 0) std::cerr << " (angles)";
1608  std::cerr << std::endl;
1609  std::cerr << " r1=" << r1;
1610  std::cerr << " r2=" << r2;
1611  std::cerr << " dz=" << dz << " sPhi=" << sPhi << " dPhi=" << dPhi
1612  << std::endl;
1613  return;
1614  }
1615 
1616  // P R E P A R E T W O P O L Y L I N E S
1617 
1618  G4int n = GetNumberOfRotationSteps();
1619  G4double dl = (r2 - r1) / n;
1620  G4double k1 = (r2*r2 - r1*r1) / 2 / dz;
1621  G4double k2 = (r2*r2 + r1*r1) / 2;
1622 
1623  G4double *zz = new G4double[n + 2], *rr = new G4double[n + 2];
1624 
1625  zz[0] = dz;
1626  rr[0] = r2;
1627 
1628  for(G4int i = 1; i < n - 1; i++)
1629  {
1630  rr[i] = rr[i-1] - dl;
1631  zz[i] = (rr[i]*rr[i] - k2) / k1;
1632  if(rr[i] < 0)
1633  {
1634  rr[i] = 0;
1635  zz[i] = 0;
1636  }
1637  }
1638 
1639  zz[n-1] = -dz;
1640  rr[n-1] = r1;
1641 
1642  zz[n] = dz;
1643  rr[n] = 0;
1644 
1645  zz[n+1] = -dz;
1646  rr[n+1] = 0;
1647 
1648  // R O T A T E P O L Y L I N E S
1649 
1650  RotateAroundZ(0, phi1, dphi, n, 2, zz, rr, -1, -1);
1651  SetReferences();
1652 
1653  delete [] zz;
1654  delete [] rr;
1655 }
1656 
1657 HepPolyhedronParaboloid::~HepPolyhedronParaboloid() {}
1658 
1659 HepPolyhedronHype::HepPolyhedronHype(G4double r1,
1660  G4double r2,
1661  G4double sqrtan1,
1662  G4double sqrtan2,
1663  G4double halfZ)
1664 /***********************************************************************
1665  * *
1666  * Name: HepPolyhedronHype Date: 14.04.08 *
1667  * Author: Tatiana Nikitina (CERN) Revised: 14.04.08 *
1668  * *
1669  * Function: Constructor for Hype *
1670  * *
1671  * Input: r1 - inside radius at z=0 *
1672  * r2 - outside radiuses at z=0 *
1673  * sqrtan1 - sqr of tan of Inner Stereo Angle *
1674  * sqrtan2 - sqr of tan of Outer Stereo Angle *
1675  * halfZ - half length in Z *
1676  * *
1677  ***********************************************************************/
1678 {
1679  static const G4double wholeCircle=twopi;
1680 
1681  // C H E C K I N P U T P A R A M E T E R S
1682 
1683  G4int k = 0;
1684  if (r2 < 0. || r1 < 0. ) k = 1;
1685  if (r1 > r2 ) k = 1;
1686  if (r1 == r2) k = 1;
1687 
1688  if (halfZ <= 0.) k += 2;
1689 
1690  if (sqrtan1<0.||sqrtan2<0.) k += 4;
1691 
1692  if (k != 0)
1693  {
1694  std::cerr << "HepPolyhedronHype: error in input parameters";
1695  if ((k & 1) != 0) std::cerr << " (radiuses)";
1696  if ((k & 2) != 0) std::cerr << " (half-length)";
1697  if ((k & 4) != 0) std::cerr << " (angles)";
1698  std::cerr << std::endl;
1699  std::cerr << " r1=" << r1 << " r2=" << r2;
1700  std::cerr << " halfZ=" << halfZ << " sqrTan1=" << sqrtan1
1701  << " sqrTan2=" << sqrtan2
1702  << std::endl;
1703  return;
1704  }
1705 
1706  // P R E P A R E T W O P O L Y L I N E S
1707 
1708  G4int n = GetNumberOfRotationSteps();
1709  G4double dz = 2.*halfZ / n;
1710  G4double k1 = r1*r1;
1711  G4double k2 = r2*r2;
1712 
1713  G4double *zz = new G4double[n+n+1], *rr = new G4double[n+n+1];
1714 
1715  zz[0] = halfZ;
1716  rr[0] = std::sqrt(sqrtan2*halfZ*halfZ+k2);
1717 
1718  for(G4int i = 1; i < n-1; i++)
1719  {
1720  zz[i] = zz[i-1] - dz;
1721  rr[i] =std::sqrt(sqrtan2*zz[i]*zz[i]+k2);
1722  }
1723 
1724  zz[n-1] = -halfZ;
1725  rr[n-1] = rr[0];
1726 
1727  zz[n] = halfZ;
1728  rr[n] = std::sqrt(sqrtan1*halfZ*halfZ+k1);
1729 
1730  for(G4int i = n+1; i < n+n; i++)
1731  {
1732  zz[i] = zz[i-1] - dz;
1733  rr[i] =std::sqrt(sqrtan1*zz[i]*zz[i]+k1);
1734  }
1735  zz[n+n] = -halfZ;
1736  rr[n+n] = rr[n];
1737 
1738  // R O T A T E P O L Y L I N E S
1739 
1740  RotateAroundZ(0, 0., wholeCircle, n, n, zz, rr, -1, -1);
1741  SetReferences();
1742 
1743  delete [] zz;
1744  delete [] rr;
1745 }
1746 
1747 HepPolyhedronHype::~HepPolyhedronHype() {}
1748 
1749 HepPolyhedronCons::HepPolyhedronCons(G4double Rmn1,
1750  G4double Rmx1,
1751  G4double Rmn2,
1752  G4double Rmx2,
1753  G4double Dz,
1754  G4double Phi1,
1755  G4double Dphi)
1756 /***********************************************************************
1757  * *
1758  * Name: HepPolyhedronCons::HepPolyhedronCons Date: 15.12.96 *
1759  * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
1760  * *
1761  * Function: Constructor for CONS, TUBS, CONE, TUBE *
1762  * *
1763  * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
1764  * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
1765  * Dz - half length in Z *
1766  * Phi1 - starting angle of the segment *
1767  * Dphi - segment range *
1768  * *
1769  ***********************************************************************/
1770 {
1771  static const G4double wholeCircle=twopi;
1772 
1773  // C H E C K I N P U T P A R A M E T E R S
1774 
1775  G4int k = 0;
1776  if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
1777  if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
1778  if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
1779 
1780  if (Dz <= 0.) k += 2;
1781 
1782  G4double phi1, phi2, dphi;
1783  if (Dphi < 0.) {
1784  phi2 = Phi1; phi1 = phi2 - Dphi;
1785  }else if (Dphi == 0.) {
1786  phi1 = Phi1; phi2 = phi1 + wholeCircle;
1787  }else{
1788  phi1 = Phi1; phi2 = phi1 + Dphi;
1789  }
1790  dphi = phi2 - phi1;
1791  if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1792  if (dphi > wholeCircle) k += 4;
1793 
1794  if (k != 0) {
1795  std::cerr << "HepPolyhedronCone(s)/Tube(s): error in input parameters";
1796  if ((k & 1) != 0) std::cerr << " (radiuses)";
1797  if ((k & 2) != 0) std::cerr << " (half-length)";
1798  if ((k & 4) != 0) std::cerr << " (angles)";
1799  std::cerr << std::endl;
1800  std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
1801  std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
1802  std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
1803  << std::endl;
1804  return;
1805  }
1806 
1807  // P R E P A R E T W O P O L Y L I N E S
1808 
1809  G4double zz[4], rr[4];
1810  zz[0] = Dz;
1811  zz[1] = -Dz;
1812  zz[2] = Dz;
1813  zz[3] = -Dz;
1814  rr[0] = Rmx2;
1815  rr[1] = Rmx1;
1816  rr[2] = Rmn2;
1817  rr[3] = Rmn1;
1818 
1819  // R O T A T E P O L Y L I N E S
1820 
1821  RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1);
1822  SetReferences();
1823 }
1824 
1825 HepPolyhedronCons::~HepPolyhedronCons() {}
1826 
1827 HepPolyhedronCone::HepPolyhedronCone(G4double Rmn1, G4double Rmx1,
1828  G4double Rmn2, G4double Rmx2,
1829  G4double Dz) :
1830  HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
1831 
1832 HepPolyhedronCone::~HepPolyhedronCone() {}
1833 
1834 HepPolyhedronTubs::HepPolyhedronTubs(G4double Rmin, G4double Rmax,
1835  G4double Dz,
1836  G4double Phi1, G4double Dphi)
1837  : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
1838 
1839 HepPolyhedronTubs::~HepPolyhedronTubs() {}
1840 
1841 HepPolyhedronTube::HepPolyhedronTube (G4double Rmin, G4double Rmax,
1842  G4double Dz)
1843  : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
1844 
1845 HepPolyhedronTube::~HepPolyhedronTube () {}
1846 
1847 HepPolyhedronPgon::HepPolyhedronPgon(G4double phi,
1848  G4double dphi,
1849  G4int npdv,
1850  G4int nz,
1851  const G4double *z,
1852  const G4double *rmin,
1853  const G4double *rmax)
1854 /***********************************************************************
1855  * *
1856  * Name: HepPolyhedronPgon Date: 09.12.96 *
1857  * Author: E.Chernyaev Revised: *
1858  * *
1859  * Function: Constructor of polyhedron for PGON, PCON *
1860  * *
1861  * Input: phi - initial phi *
1862  * dphi - delta phi *
1863  * npdv - number of steps along phi *
1864  * nz - number of z-planes (at least two) *
1865  * z[] - z coordinates of the slices *
1866  * rmin[] - smaller r at the slices *
1867  * rmax[] - bigger r at the slices *
1868  * *
1869  ***********************************************************************/
1870 {
1871  // C H E C K I N P U T P A R A M E T E R S
1872 
1873  if (dphi <= 0. || dphi > twopi) {
1874  std::cerr
1875  << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
1876  << std::endl;
1877  return;
1878  }
1879 
1880  if (nz < 2) {
1881  std::cerr
1882  << "HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
1883  << std::endl;
1884  return;
1885  }
1886 
1887  if (npdv < 0) {
1888  std::cerr
1889  << "HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
1890  << std::endl;
1891  return;
1892  }
1893 
1894  G4int i;
1895  for (i=0; i<nz; i++) {
1896  if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
1897  std::cerr
1898  << "HepPolyhedronPgon: error in radiuses rmin[" << i << "]="
1899  << rmin[i] << " rmax[" << i << "]=" << rmax[i]
1900  << std::endl;
1901  return;
1902  }
1903  }
1904 
1905  // P R E P A R E T W O P O L Y L I N E S
1906 
1907  G4double *zz, *rr;
1908  zz = new G4double[2*nz];
1909  rr = new G4double[2*nz];
1910 
1911  if (z[0] > z[nz-1]) {
1912  for (i=0; i<nz; i++) {
1913  zz[i] = z[i];
1914  rr[i] = rmax[i];
1915  zz[i+nz] = z[i];
1916  rr[i+nz] = rmin[i];
1917  }
1918  }else{
1919  for (i=0; i<nz; i++) {
1920  zz[i] = z[nz-i-1];
1921  rr[i] = rmax[nz-i-1];
1922  zz[i+nz] = z[nz-i-1];
1923  rr[i+nz] = rmin[nz-i-1];
1924  }
1925  }
1926 
1927  // R O T A T E P O L Y L I N E S
1928 
1929  RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1);
1930  SetReferences();
1931 
1932  delete [] zz;
1933  delete [] rr;
1934 }
1935 
1936 HepPolyhedronPgon::~HepPolyhedronPgon() {}
1937 
1938 HepPolyhedronPcon::HepPolyhedronPcon(G4double phi, G4double dphi, G4int nz,
1939  const G4double *z,
1940  const G4double *rmin,
1941  const G4double *rmax)
1942  : HepPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
1943 
1944 HepPolyhedronPcon::~HepPolyhedronPcon() {}
1945 
1946 HepPolyhedronSphere::HepPolyhedronSphere(G4double rmin, G4double rmax,
1947  G4double phi, G4double dphi,
1948  G4double the, G4double dthe)
1949 /***********************************************************************
1950  * *
1951  * Name: HepPolyhedronSphere Date: 11.12.96 *
1952  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1953  * *
1954  * Function: Constructor of polyhedron for SPHERE *
1955  * *
1956  * Input: rmin - internal radius *
1957  * rmax - external radius *
1958  * phi - initial phi *
1959  * dphi - delta phi *
1960  * the - initial theta *
1961  * dthe - delta theta *
1962  * *
1963  ***********************************************************************/
1964 {
1965  // C H E C K I N P U T P A R A M E T E R S
1966 
1967  if (dphi <= 0. || dphi > twopi) {
1968  std::cerr
1969  << "HepPolyhedronSphere: wrong delta phi = " << dphi
1970  << std::endl;
1971  return;
1972  }
1973 
1974  if (the < 0. || the > pi) {
1975  std::cerr
1976  << "HepPolyhedronSphere: wrong theta = " << the
1977  << std::endl;
1978  return;
1979  }
1980 
1981  if (dthe <= 0. || dthe > pi) {
1982  std::cerr
1983  << "HepPolyhedronSphere: wrong delta theta = " << dthe
1984  << std::endl;
1985  return;
1986  }
1987 
1988  if (the+dthe > pi) {
1989  std::cerr
1990  << "HepPolyhedronSphere: wrong theta + delta theta = "
1991  << the << " " << dthe
1992  << std::endl;
1993  return;
1994  }
1995 
1996  if (rmin < 0. || rmin >= rmax) {
1997  std::cerr
1998  << "HepPolyhedronSphere: error in radiuses"
1999  << " rmin=" << rmin << " rmax=" << rmax
2000  << std::endl;
2001  return;
2002  }
2003 
2004  // P R E P A R E T W O P O L Y L I N E S
2005 
2006  G4int nds = (GetNumberOfRotationSteps() + 1) / 2;
2007  G4int np1 = G4int(dthe*nds/pi+.5) + 1;
2008  if (np1 <= 1) np1 = 2;
2009  G4int np2 = rmin < spatialTolerance ? 1 : np1;
2010 
2011  G4double *zz, *rr;
2012  zz = new G4double[np1+np2];
2013  rr = new G4double[np1+np2];
2014 
2015  G4double a = dthe/(np1-1);
2016  G4double cosa, sina;
2017  for (G4int i=0; i<np1; i++) {
2018  cosa = std::cos(the+i*a);
2019  sina = std::sin(the+i*a);
2020  zz[i] = rmax*cosa;
2021  rr[i] = rmax*sina;
2022  if (np2 > 1) {
2023  zz[i+np1] = rmin*cosa;
2024  rr[i+np1] = rmin*sina;
2025  }
2026  }
2027  if (np2 == 1) {
2028  zz[np1] = 0.;
2029  rr[np1] = 0.;
2030  }
2031 
2032  // R O T A T E P O L Y L I N E S
2033 
2034  RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
2035  SetReferences();
2036 
2037  delete [] zz;
2038  delete [] rr;
2039 }
2040 
2041 HepPolyhedronSphere::~HepPolyhedronSphere() {}
2042 
2043 HepPolyhedronTorus::HepPolyhedronTorus(G4double rmin,
2044  G4double rmax,
2045  G4double rtor,
2046  G4double phi,
2047  G4double dphi)
2048 /***********************************************************************
2049  * *
2050  * Name: HepPolyhedronTorus Date: 11.12.96 *
2051  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
2052  * *
2053  * Function: Constructor of polyhedron for TORUS *
2054  * *
2055  * Input: rmin - internal radius *
2056  * rmax - external radius *
2057  * rtor - radius of torus *
2058  * phi - initial phi *
2059  * dphi - delta phi *
2060  * *
2061  ***********************************************************************/
2062 {
2063  // C H E C K I N P U T P A R A M E T E R S
2064 
2065  if (dphi <= 0. || dphi > twopi) {
2066  std::cerr
2067  << "HepPolyhedronTorus: wrong delta phi = " << dphi
2068  << std::endl;
2069  return;
2070  }
2071 
2072  if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2073  std::cerr
2074  << "HepPolyhedronTorus: error in radiuses"
2075  << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
2076  << std::endl;
2077  return;
2078  }
2079 
2080  // P R E P A R E T W O P O L Y L I N E S
2081 
2082  G4int np1 = GetNumberOfRotationSteps();
2083  assert(np1>0);
2084  G4int np2 = rmin < spatialTolerance ? 1 : np1;
2085 
2086  G4double *zz, *rr;
2087  zz = new G4double[np1+np2];
2088  rr = new G4double[np1+np2];
2089 
2090  G4double a = twopi/np1;
2091  G4double cosa, sina;
2092  for (G4int i=0; i<np1; i++) {
2093  cosa = std::cos(i*a);
2094  sina = std::sin(i*a);
2095  zz[i] = rmax*cosa;
2096  rr[i] = rtor+rmax*sina;
2097  if (np2 > 1) {
2098  zz[i+np1] = rmin*cosa;
2099  rr[i+np1] = rtor+rmin*sina;
2100  }
2101  }
2102  if (np2 == 1) {
2103  zz[np1] = 0.;
2104  rr[np1] = rtor;
2105  np2 = -1;
2106  }
2107 
2108  // R O T A T E P O L Y L I N E S
2109 
2110  RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1);
2111  SetReferences();
2112 
2113  delete [] zz;
2114  delete [] rr;
2115 }
2116 
2117 HepPolyhedronTorus::~HepPolyhedronTorus() {}
2118 
2119 HepPolyhedronEllipsoid::HepPolyhedronEllipsoid(G4double ax, G4double by,
2120  G4double cz, G4double zCut1,
2121  G4double zCut2)
2122 /***********************************************************************
2123  * *
2124  * Name: HepPolyhedronEllipsoid Date: 25.02.05 *
2125  * Author: G.Guerrieri Revised: *
2126  * *
2127  * Function: Constructor of polyhedron for ELLIPSOID *
2128  * *
2129  * Input: ax - semiaxis x *
2130  * by - semiaxis y *
2131  * cz - semiaxis z *
2132  * zCut1 - lower cut plane level (solid lies above this plane) *
2133  * zCut2 - upper cut plane level (solid lies below this plane) *
2134  * *
2135  ***********************************************************************/
2136 {
2137  // C H E C K I N P U T P A R A M E T E R S
2138 
2139  if (zCut1 >= cz || zCut2 <= -cz || zCut1 > zCut2) {
2140  std::cerr << "HepPolyhedronEllipsoid: wrong zCut1 = " << zCut1
2141  << " zCut2 = " << zCut2
2142  << " for given cz = " << cz << std::endl;
2143  return;
2144  }
2145  if (cz <= 0.0) {
2146  std::cerr << "HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
2147  << std::endl;
2148  return;
2149  }
2150 
2151  G4double dthe;
2152  G4double sthe;
2153  G4int cutflag;
2154  cutflag= 0;
2155  if (zCut2 >= cz)
2156  {
2157  sthe= 0.0;
2158  }
2159  else
2160  {
2161  sthe= std::acos(zCut2/cz);
2162  cutflag++;
2163  }
2164  if (zCut1 <= -cz)
2165  {
2166  dthe= pi - sthe;
2167  }
2168  else
2169  {
2170  dthe= std::acos(zCut1/cz)-sthe;
2171  cutflag++;
2172  }
2173 
2174  // P R E P A R E T W O P O L Y L I N E S
2175  // generate sphere of radius cz first, then rescale x and y later
2176 
2177  G4int nds = (GetNumberOfRotationSteps() + 1) / 2;
2178  G4int np1 = G4int(dthe*nds/pi) + 2 + cutflag;
2179 
2180  G4double *zz, *rr;
2181  zz = new G4double[np1+1];
2182  rr = new G4double[np1+1];
2183  if (!zz || !rr)
2184  {
2185  G4Exception("HepPolyhedronEllipsoid::HepPolyhedronEllipsoid",
2186  "greps1002", FatalException, "Out of memory");
2187  }
2188 
2189  G4double a = dthe/(np1-cutflag-1);
2190  G4double cosa, sina;
2191  G4int j=0;
2192  if (sthe > 0.0)
2193  {
2194  zz[j]= zCut2;
2195  rr[j]= 0.;
2196  j++;
2197  }
2198  for (G4int i=0; i<np1-cutflag; i++) {
2199  cosa = std::cos(sthe+i*a);
2200  sina = std::sin(sthe+i*a);
2201  zz[j] = cz*cosa;
2202  rr[j] = cz*sina;
2203  j++;
2204  }
2205  if (j < np1)
2206  {
2207  zz[j]= zCut1;
2208  rr[j]= 0.;
2209  j++;
2210  }
2211  if (j > np1)
2212  {
2213  std::cerr << "Logic error in HepPolyhedronEllipsoid, memory corrupted!"
2214  << std::endl;
2215  }
2216  if (j < np1)
2217  {
2218  std::cerr << "Warning: logic error in HepPolyhedronEllipsoid."
2219  << std::endl;
2220  np1= j;
2221  }
2222  zz[j] = 0.;
2223  rr[j] = 0.;
2224 
2225 
2226  // R O T A T E P O L Y L I N E S
2227 
2228  RotateAroundZ(0, 0.0, twopi, np1, 1, zz, rr, -1, 1);
2229  SetReferences();
2230 
2231  delete [] zz;
2232  delete [] rr;
2233 
2234  // rescale x and y vertex coordinates
2235  {
2236  G4Point3D * p= pV;
2237  for (G4int i=0; i<nvert; i++, p++) {
2238  p->setX( p->x() * ax/cz );
2239  p->setY( p->y() * by/cz );
2240  }
2241  }
2242 }
2243 
2244 HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() {}
2245 
2246 HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone(G4double ax,
2247  G4double ay,
2248  G4double h,
2249  G4double zTopCut)
2250 /***********************************************************************
2251  * *
2252  * Name: HepPolyhedronEllipticalCone Date: 8.9.2005 *
2253  * Author: D.Anninos Revised: 9.9.2005 *
2254  * *
2255  * Function: Constructor for EllipticalCone *
2256  * *
2257  * Input: ax, ay - X & Y semi axes at z = 0 *
2258  * h - height of full cone *
2259  * zTopCut - Top Cut in Z Axis *
2260  * *
2261  ***********************************************************************/
2262 {
2263  // C H E C K I N P U T P A R A M E T E R S
2264 
2265  G4int k = 0;
2266  if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
2267 
2268  if (k != 0) {
2269  std::cerr << "HepPolyhedronCone: error in input parameters";
2270  std::cerr << std::endl;
2271  return;
2272  }
2273 
2274  // P R E P A R E T W O P O L Y L I N E S
2275 
2276  zTopCut = (h >= zTopCut ? zTopCut : h);
2277 
2278  G4double *zz, *rr;
2279  zz = new G4double[4];
2280  rr = new G4double[4];
2281  zz[0] = zTopCut;
2282  zz[1] = -zTopCut;
2283  zz[2] = zTopCut;
2284  zz[3] = -zTopCut;
2285  rr[0] = (h-zTopCut);
2286  rr[1] = (h+zTopCut);
2287  rr[2] = 0.;
2288  rr[3] = 0.;
2289 
2290  // R O T A T E P O L Y L I N E S
2291 
2292  RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, -1, -1);
2293  SetReferences();
2294 
2295  delete [] zz;
2296  delete [] rr;
2297 
2298  // rescale x and y vertex coordinates
2299  {
2300  G4Point3D * p= pV;
2301  for (G4int i=0; i<nvert; i++, p++) {
2302  p->setX( p->x() * ax );
2303  p->setY( p->y() * ay );
2304  }
2305  }
2306 }
2307 
2308 HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() {}
2309 
2310 G4ThreadLocal G4int HepPolyhedron::fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
2311 /***********************************************************************
2312  * *
2313  * Name: HepPolyhedron::fNumberOfRotationSteps Date: 24.06.97 *
2314  * Author: J.Allison (Manchester University) Revised: *
2315  * *
2316  * Function: Number of steps for whole circle *
2317  * *
2318  ***********************************************************************/
2319 
2320 #include "BooleanProcessor.src"
2321 
2322 HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const
2323 /***********************************************************************
2324  * *
2325  * Name: HepPolyhedron::add Date: 19.03.00 *
2326  * Author: E.Chernyaev Revised: *
2327  * *
2328  * Function: Boolean "union" of two polyhedra *
2329  * *
2330  ***********************************************************************/
2331 {
2332  G4int ierr;
2333  BooleanProcessor processor;
2334  return processor.execute(OP_UNION, *this, p,ierr);
2335 }
2336 
2337 HepPolyhedron HepPolyhedron::intersect(const HepPolyhedron & p) const
2338 /***********************************************************************
2339  * *
2340  * Name: HepPolyhedron::intersect Date: 19.03.00 *
2341  * Author: E.Chernyaev Revised: *
2342  * *
2343  * Function: Boolean "intersection" of two polyhedra *
2344  * *
2345  ***********************************************************************/
2346 {
2347  G4int ierr;
2348  BooleanProcessor processor;
2349  return processor.execute(OP_INTERSECTION, *this, p,ierr);
2350 }
2351 
2352 HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const
2353 /***********************************************************************
2354  * *
2355  * Name: HepPolyhedron::add Date: 19.03.00 *
2356  * Author: E.Chernyaev Revised: *
2357  * *
2358  * Function: Boolean "subtraction" of "p" from "this" *
2359  * *
2360  ***********************************************************************/
2361 {
2362  G4int ierr;
2363  BooleanProcessor processor;
2364  return processor.execute(OP_SUBTRACTION, *this, p,ierr);
2365 }
2366 
2367 //NOTE : include the code of HepPolyhedronProcessor here
2368 // since there is no BooleanProcessor.h
2369 
2370 #undef INTERSECTION
2371 
2372 #include "HepPolyhedronProcessor.src"
2373