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