Geant4  10.01.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 83392 2014-08-21 14:36:35Z 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  if (!freeList) {
725  std::cerr
726  << "Polyhedron::SetReferences: bad link "
727  << std::endl;
728  break;
729  }
730  freeList = freeList->next;
731  cur = headList[k1];
732  cur->next = 0;
733  cur->v2 = k2;
734  cur->iface = iface;
735  cur->iedge = iedge;
736  continue;
737  }
738 
739  if (cur->v2 == k2) {
740  headList[k1] = cur->next;
741  cur->next = freeList;
742  freeList = cur;
743  pF[iface].edge[iedge].f = cur->iface;
744  pF[cur->iface].edge[cur->iedge].f = iface;
745  i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
746  i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
747  if (i1 != i2) {
748  std::cerr
749  << "Polyhedron::SetReferences: different edge visibility "
750  << iface << "/" << iedge << "/"
751  << pF[iface].edge[iedge].v << " and "
752  << cur->iface << "/" << cur->iedge << "/"
753  << pF[cur->iface].edge[cur->iedge].v
754  << std::endl;
755  }
756  continue;
757  }
758 
759  // check List itself
760  for (;;) {
761  prev = cur;
762  cur = prev->next;
763  if (cur == 0) {
764  prev->next = freeList;
765  if (!freeList) {
766  std::cerr
767  << "Polyhedron::SetReferences: bad link "
768  << std::endl;
769  break;
770  }
771  freeList = freeList->next;
772  cur = prev->next;
773  cur->next = 0;
774  cur->v2 = k2;
775  cur->iface = iface;
776  cur->iedge = iedge;
777  break;
778  }
779 
780  if (cur->v2 == k2) {
781  prev->next = cur->next;
782  cur->next = freeList;
783  freeList = cur;
784  pF[iface].edge[iedge].f = cur->iface;
785  pF[cur->iface].edge[cur->iedge].f = iface;
786  i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
787  i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
788  if (i1 != i2) {
789  std::cerr
790  << "Polyhedron::SetReferences: different edge visibility "
791  << iface << "/" << iedge << "/"
792  << pF[iface].edge[iedge].v << " and "
793  << cur->iface << "/" << cur->iedge << "/"
794  << pF[cur->iface].edge[cur->iedge].v
795  << std::endl;
796  }
797  break;
798  }
799  }
800  }
801  }
802 
803  // 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
804 
805  for (i=0; i<nvert; i++) {
806  if (headList[i] != 0) {
807  std::cerr
808  << "Polyhedron::SetReferences: List " << i << " is not empty"
809  << std::endl;
810  }
811  }
812 
813  // F R E E M E M O R Y
814 
815  delete [] edgeList;
816  delete [] headList;
817 }
818 
819 void HepPolyhedron::InvertFacets()
820 /***********************************************************************
821  * *
822  * Name: HepPolyhedron::InvertFacets Date: 01.12.99 *
823  * Author: E.Chernyaev Revised: *
824  * *
825  * Function: Invert the order of the nodes in the facets *
826  * *
827  ***********************************************************************/
828 {
829  if (nface <= 0) return;
830  G4int i, k, nnode, v[4],f[4];
831  for (i=1; i<=nface; i++) {
832  nnode = (pF[i].edge[3].v == 0) ? 3 : 4;
833  for (k=0; k<nnode; k++) {
834  v[k] = (k+1 == nnode) ? pF[i].edge[0].v : pF[i].edge[k+1].v;
835  if (v[k] * pF[i].edge[k].v < 0) v[k] = -v[k];
836  f[k] = pF[i].edge[k].f;
837  }
838  for (k=0; k<nnode; k++) {
839  pF[i].edge[nnode-1-k].v = v[k];
840  pF[i].edge[nnode-1-k].f = f[k];
841  }
842  }
843 }
844 
845 HepPolyhedron & HepPolyhedron::Transform(const G4Transform3D &t)
846 /***********************************************************************
847  * *
848  * Name: HepPolyhedron::Transform Date: 01.12.99 *
849  * Author: E.Chernyaev Revised: *
850  * *
851  * Function: Make transformation of the polyhedron *
852  * *
853  ***********************************************************************/
854 {
855  if (nvert > 0) {
856  for (G4int i=1; i<=nvert; i++) { pV[i] = t * pV[i]; }
857 
858  // C H E C K D E T E R M I N A N T A N D
859  // 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
860 
861  G4Vector3D d = t * G4Vector3D(0,0,0);
862  G4Vector3D x = t * G4Vector3D(1,0,0) - d;
863  G4Vector3D y = t * G4Vector3D(0,1,0) - d;
864  G4Vector3D z = t * G4Vector3D(0,0,1) - d;
865  if ((x.cross(y))*z < 0) InvertFacets();
866  }
867  return *this;
868 }
869 
870 G4bool HepPolyhedron::GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
871 /***********************************************************************
872  * *
873  * Name: HepPolyhedron::GetNextVertexIndex Date: 03.09.96 *
874  * Author: Yasuhide Sawada Revised: *
875  * *
876  * Function: *
877  * *
878  ***********************************************************************/
879 {
880  static G4ThreadLocal G4int iFace = 1;
881  static G4ThreadLocal G4int iQVertex = 0;
882  G4int vIndex = pF[iFace].edge[iQVertex].v;
883 
884  edgeFlag = (vIndex > 0) ? 1 : 0;
885  index = std::abs(vIndex);
886 
887  if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
888  iQVertex = 0;
889  if (++iFace > nface) iFace = 1;
890  return false; // Last Edge
891  }else{
892  ++iQVertex;
893  return true; // not Last Edge
894  }
895 }
896 
897 G4Point3D HepPolyhedron::GetVertex(G4int index) const
898 /***********************************************************************
899  * *
900  * Name: HepPolyhedron::GetVertex Date: 03.09.96 *
901  * Author: Yasuhide Sawada Revised: 17.11.99 *
902  * *
903  * Function: Get vertex of the index. *
904  * *
905  ***********************************************************************/
906 {
907  if (index <= 0 || index > nvert) {
908  std::cerr
909  << "HepPolyhedron::GetVertex: irrelevant index " << index
910  << std::endl;
911  return G4Point3D();
912  }
913  return pV[index];
914 }
915 
916 G4bool
917 HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag) const
918 /***********************************************************************
919  * *
920  * Name: HepPolyhedron::GetNextVertex Date: 22.07.96 *
921  * Author: John Allison Revised: *
922  * *
923  * Function: Get vertices of the quadrilaterals in order for each *
924  * face in face order. Returns false when finished each *
925  * face. *
926  * *
927  ***********************************************************************/
928 {
929  G4int index;
930  G4bool rep = GetNextVertexIndex(index, edgeFlag);
931  vertex = pV[index];
932  return rep;
933 }
934 
935 G4bool HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag,
936  G4Normal3D &normal) const
937 /***********************************************************************
938  * *
939  * Name: HepPolyhedron::GetNextVertex Date: 26.11.99 *
940  * Author: E.Chernyaev Revised: *
941  * *
942  * Function: Get vertices with normals of the quadrilaterals in order *
943  * for each face in face order. *
944  * Returns false when finished each face. *
945  * *
946  ***********************************************************************/
947 {
948  static G4ThreadLocal G4int iFace = 1;
949  static G4ThreadLocal G4int iNode = 0;
950 
951  if (nface == 0) return false; // empty polyhedron
952 
953  G4int k = pF[iFace].edge[iNode].v;
954  if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
955  vertex = pV[k];
956  normal = FindNodeNormal(iFace,k);
957  if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
958  iNode = 0;
959  if (++iFace > nface) iFace = 1;
960  return false; // last node
961  }else{
962  ++iNode;
963  return true; // not last node
964  }
965 }
966 
967 G4bool HepPolyhedron::GetNextEdgeIndeces(G4int &i1, G4int &i2, G4int &edgeFlag,
968  G4int &iface1, G4int &iface2) const
969 /***********************************************************************
970  * *
971  * Name: HepPolyhedron::GetNextEdgeIndeces Date: 30.09.96 *
972  * Author: E.Chernyaev Revised: 17.11.99 *
973  * *
974  * Function: Get indeces of the next edge together with indeces of *
975  * of the faces which share the edge. *
976  * Returns false when the last edge. *
977  * *
978  ***********************************************************************/
979 {
980  static G4ThreadLocal G4int iFace = 1;
981  static G4ThreadLocal G4int iQVertex = 0;
982  static G4ThreadLocal G4int iOrder = 1;
983  G4int k1, k2, kflag, kface1, kface2;
984 
985  if (iFace == 1 && iQVertex == 0) {
986  k2 = pF[nface].edge[0].v;
987  k1 = pF[nface].edge[3].v;
988  if (k1 == 0) k1 = pF[nface].edge[2].v;
989  if (std::abs(k1) > std::abs(k2)) iOrder = -1;
990  }
991 
992  do {
993  k1 = pF[iFace].edge[iQVertex].v;
994  kflag = k1;
995  k1 = std::abs(k1);
996  kface1 = iFace;
997  kface2 = pF[iFace].edge[iQVertex].f;
998  if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
999  iQVertex = 0;
1000  k2 = std::abs(pF[iFace].edge[iQVertex].v);
1001  iFace++;
1002  }else{
1003  iQVertex++;
1004  k2 = std::abs(pF[iFace].edge[iQVertex].v);
1005  }
1006  } while (iOrder*k1 > iOrder*k2);
1007 
1008  i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1009  iface1 = kface1; iface2 = kface2;
1010 
1011  if (iFace > nface) {
1012  iFace = 1; iOrder = 1;
1013  return false;
1014  }else{
1015  return true;
1016  }
1017 }
1018 
1019 G4bool
1020 HepPolyhedron::GetNextEdgeIndeces(G4int &i1, G4int &i2, G4int &edgeFlag) const
1021 /***********************************************************************
1022  * *
1023  * Name: HepPolyhedron::GetNextEdgeIndeces Date: 17.11.99 *
1024  * Author: E.Chernyaev Revised: *
1025  * *
1026  * Function: Get indeces of the next edge. *
1027  * Returns false when the last edge. *
1028  * *
1029  ***********************************************************************/
1030 {
1031  G4int kface1, kface2;
1032  return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
1033 }
1034 
1035 G4bool
1036 HepPolyhedron::GetNextEdge(G4Point3D &p1,
1037  G4Point3D &p2,
1038  G4int &edgeFlag) const
1039 /***********************************************************************
1040  * *
1041  * Name: HepPolyhedron::GetNextEdge Date: 30.09.96 *
1042  * Author: E.Chernyaev Revised: *
1043  * *
1044  * Function: Get next edge. *
1045  * Returns false when the last edge. *
1046  * *
1047  ***********************************************************************/
1048 {
1049  G4int i1,i2;
1050  G4bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
1051  p1 = pV[i1];
1052  p2 = pV[i2];
1053  return rep;
1054 }
1055 
1056 G4bool
1057 HepPolyhedron::GetNextEdge(G4Point3D &p1, G4Point3D &p2,
1058  G4int &edgeFlag, G4int &iface1, G4int &iface2) const
1059 /***********************************************************************
1060  * *
1061  * Name: HepPolyhedron::GetNextEdge Date: 17.11.99 *
1062  * Author: E.Chernyaev Revised: *
1063  * *
1064  * Function: Get next edge with indeces of the faces which share *
1065  * the edge. *
1066  * Returns false when the last edge. *
1067  * *
1068  ***********************************************************************/
1069 {
1070  G4int i1,i2;
1071  G4bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
1072  p1 = pV[i1];
1073  p2 = pV[i2];
1074  return rep;
1075 }
1076 
1077 void HepPolyhedron::GetFacet(G4int iFace, G4int &n, G4int *iNodes,
1078  G4int *edgeFlags, G4int *iFaces) const
1079 /***********************************************************************
1080  * *
1081  * Name: HepPolyhedron::GetFacet Date: 15.12.99 *
1082  * Author: E.Chernyaev Revised: *
1083  * *
1084  * Function: Get face by index *
1085  * *
1086  ***********************************************************************/
1087 {
1088  if (iFace < 1 || iFace > nface) {
1089  std::cerr
1090  << "HepPolyhedron::GetFacet: irrelevant index " << iFace
1091  << std::endl;
1092  n = 0;
1093  }else{
1094  G4int i, k;
1095  for (i=0; i<4; i++) {
1096  k = pF[iFace].edge[i].v;
1097  if (k == 0) break;
1098  if (iFaces != 0) iFaces[i] = pF[iFace].edge[i].f;
1099  if (k > 0) {
1100  iNodes[i] = k;
1101  if (edgeFlags != 0) edgeFlags[i] = 1;
1102  }else{
1103  iNodes[i] = -k;
1104  if (edgeFlags != 0) edgeFlags[i] = -1;
1105  }
1106  }
1107  n = i;
1108  }
1109 }
1110 
1111 void HepPolyhedron::GetFacet(G4int index, G4int &n, G4Point3D *nodes,
1112  G4int *edgeFlags, G4Normal3D *normals) const
1113 /***********************************************************************
1114  * *
1115  * Name: HepPolyhedron::GetFacet Date: 17.11.99 *
1116  * Author: E.Chernyaev Revised: *
1117  * *
1118  * Function: Get face by index *
1119  * *
1120  ***********************************************************************/
1121 {
1122  G4int iNodes[4];
1123  GetFacet(index, n, iNodes, edgeFlags);
1124  if (n != 0) {
1125  for (G4int i=0; i<n; i++) {
1126  nodes[i] = pV[iNodes[i]];
1127  if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1128  }
1129  }
1130 }
1131 
1132 G4bool
1133 HepPolyhedron::GetNextFacet(G4int &n, G4Point3D *nodes,
1134  G4int *edgeFlags, G4Normal3D *normals) const
1135 /***********************************************************************
1136  * *
1137  * Name: HepPolyhedron::GetNextFacet Date: 19.11.99 *
1138  * Author: E.Chernyaev Revised: *
1139  * *
1140  * Function: Get next face with normals of unit length at the nodes. *
1141  * Returns false when finished all faces. *
1142  * *
1143  ***********************************************************************/
1144 {
1145  static G4ThreadLocal G4int iFace = 1;
1146 
1147  if (edgeFlags == 0) {
1148  GetFacet(iFace, n, nodes);
1149  }else if (normals == 0) {
1150  GetFacet(iFace, n, nodes, edgeFlags);
1151  }else{
1152  GetFacet(iFace, n, nodes, edgeFlags, normals);
1153  }
1154 
1155  if (++iFace > nface) {
1156  iFace = 1;
1157  return false;
1158  }else{
1159  return true;
1160  }
1161 }
1162 
1163 G4Normal3D HepPolyhedron::GetNormal(G4int iFace) const
1164 /***********************************************************************
1165  * *
1166  * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1167  * Author: E.Chernyaev Revised: *
1168  * *
1169  * Function: Get normal of the face given by index *
1170  * *
1171  ***********************************************************************/
1172 {
1173  if (iFace < 1 || iFace > nface) {
1174  std::cerr
1175  << "HepPolyhedron::GetNormal: irrelevant index " << iFace
1176  << std::endl;
1177  return G4Normal3D();
1178  }
1179 
1180  G4int i0 = std::abs(pF[iFace].edge[0].v);
1181  G4int i1 = std::abs(pF[iFace].edge[1].v);
1182  G4int i2 = std::abs(pF[iFace].edge[2].v);
1183  G4int i3 = std::abs(pF[iFace].edge[3].v);
1184  if (i3 == 0) i3 = i0;
1185  return (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1]);
1186 }
1187 
1188 G4Normal3D HepPolyhedron::GetUnitNormal(G4int iFace) const
1189 /***********************************************************************
1190  * *
1191  * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1192  * Author: E.Chernyaev Revised: *
1193  * *
1194  * Function: Get unit normal of the face given by index *
1195  * *
1196  ***********************************************************************/
1197 {
1198  if (iFace < 1 || iFace > nface) {
1199  std::cerr
1200  << "HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
1201  << std::endl;
1202  return G4Normal3D();
1203  }
1204 
1205  G4int i0 = std::abs(pF[iFace].edge[0].v);
1206  G4int i1 = std::abs(pF[iFace].edge[1].v);
1207  G4int i2 = std::abs(pF[iFace].edge[2].v);
1208  G4int i3 = std::abs(pF[iFace].edge[3].v);
1209  if (i3 == 0) i3 = i0;
1210  return ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).unit();
1211 }
1212 
1213 G4bool HepPolyhedron::GetNextNormal(G4Normal3D &normal) const
1214 /***********************************************************************
1215  * *
1216  * Name: HepPolyhedron::GetNextNormal Date: 22.07.96 *
1217  * Author: John Allison Revised: 19.11.99 *
1218  * *
1219  * Function: Get normals of each face in face order. Returns false *
1220  * when finished all faces. *
1221  * *
1222  ***********************************************************************/
1223 {
1224  static G4ThreadLocal G4int iFace = 1;
1225  normal = GetNormal(iFace);
1226  if (++iFace > nface) {
1227  iFace = 1;
1228  return false;
1229  }else{
1230  return true;
1231  }
1232 }
1233 
1234 G4bool HepPolyhedron::GetNextUnitNormal(G4Normal3D &normal) const
1235 /***********************************************************************
1236  * *
1237  * Name: HepPolyhedron::GetNextUnitNormal Date: 16.09.96 *
1238  * Author: E.Chernyaev Revised: *
1239  * *
1240  * Function: Get normals of unit length of each face in face order. *
1241  * Returns false when finished all faces. *
1242  * *
1243  ***********************************************************************/
1244 {
1245  G4bool rep = GetNextNormal(normal);
1246  normal = normal.unit();
1247  return rep;
1248 }
1249 
1250 G4double HepPolyhedron::GetSurfaceArea() const
1251 /***********************************************************************
1252  * *
1253  * Name: HepPolyhedron::GetSurfaceArea Date: 25.05.01 *
1254  * Author: E.Chernyaev Revised: *
1255  * *
1256  * Function: Returns area of the surface of the polyhedron. *
1257  * *
1258  ***********************************************************************/
1259 {
1260  G4double srf = 0.;
1261  for (G4int iFace=1; iFace<=nface; iFace++) {
1262  G4int i0 = std::abs(pF[iFace].edge[0].v);
1263  G4int i1 = std::abs(pF[iFace].edge[1].v);
1264  G4int i2 = std::abs(pF[iFace].edge[2].v);
1265  G4int i3 = std::abs(pF[iFace].edge[3].v);
1266  if (i3 == 0) i3 = i0;
1267  srf += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).mag();
1268  }
1269  return srf/2.;
1270 }
1271 
1272 G4double HepPolyhedron::GetVolume() const
1273 /***********************************************************************
1274  * *
1275  * Name: HepPolyhedron::GetVolume Date: 25.05.01 *
1276  * Author: E.Chernyaev Revised: *
1277  * *
1278  * Function: Returns volume of the polyhedron. *
1279  * *
1280  ***********************************************************************/
1281 {
1282  G4double v = 0.;
1283  for (G4int iFace=1; iFace<=nface; iFace++) {
1284  G4int i0 = std::abs(pF[iFace].edge[0].v);
1285  G4int i1 = std::abs(pF[iFace].edge[1].v);
1286  G4int i2 = std::abs(pF[iFace].edge[2].v);
1287  G4int i3 = std::abs(pF[iFace].edge[3].v);
1288  G4Point3D pt;
1289  if (i3 == 0) {
1290  i3 = i0;
1291  pt = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);
1292  }else{
1293  pt = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0.25;
1294  }
1295  v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).dot(pt);
1296  }
1297  return v/6.;
1298 }
1299 
1300 G4int
1301 HepPolyhedron::createTwistedTrap(G4double Dz,
1302  const G4double xy1[][2],
1303  const G4double xy2[][2])
1304 /***********************************************************************
1305  * *
1306  * Name: createTwistedTrap Date: 05.11.02 *
1307  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1308  * *
1309  * Function: Creates polyhedron for twisted trapezoid *
1310  * *
1311  * Input: Dz - half-length along Z 8----7 *
1312  * xy1[2,4] - quadrilateral at Z=-Dz 5----6 ! *
1313  * xy2[2,4] - quadrilateral at Z=+Dz ! 4-!--3 *
1314  * 1----2 *
1315  * *
1316  ***********************************************************************/
1317 {
1318  AllocateMemory(12,18);
1319 
1320  pV[ 1] = G4Point3D(xy1[0][0],xy1[0][1],-Dz);
1321  pV[ 2] = G4Point3D(xy1[1][0],xy1[1][1],-Dz);
1322  pV[ 3] = G4Point3D(xy1[2][0],xy1[2][1],-Dz);
1323  pV[ 4] = G4Point3D(xy1[3][0],xy1[3][1],-Dz);
1324 
1325  pV[ 5] = G4Point3D(xy2[0][0],xy2[0][1], Dz);
1326  pV[ 6] = G4Point3D(xy2[1][0],xy2[1][1], Dz);
1327  pV[ 7] = G4Point3D(xy2[2][0],xy2[2][1], Dz);
1328  pV[ 8] = G4Point3D(xy2[3][0],xy2[3][1], Dz);
1329 
1330  pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;
1331  pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;
1332  pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;
1333  pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;
1334 
1335  enum {DUMMY, BOTTOM,
1336  LEFT_BOTTOM, LEFT_FRONT, LEFT_TOP, LEFT_BACK,
1337  BACK_BOTTOM, BACK_LEFT, BACK_TOP, BACK_RIGHT,
1338  RIGHT_BOTTOM, RIGHT_BACK, RIGHT_TOP, RIGHT_FRONT,
1339  FRONT_BOTTOM, FRONT_RIGHT, FRONT_TOP, FRONT_LEFT,
1340  TOP};
1341 
1342  pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
1343 
1344  pF[ 2]=G4Facet(4,BOTTOM, -1,LEFT_FRONT, -12,LEFT_BACK, 0,0);
1345  pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP, -12,LEFT_BOTTOM, 0,0);
1346  pF[ 4]=G4Facet(5,TOP, -8,LEFT_BACK, -12,LEFT_FRONT, 0,0);
1347  pF[ 5]=G4Facet(8,BACK_LEFT, -4,LEFT_BOTTOM, -12,LEFT_TOP, 0,0);
1348 
1349  pF[ 6]=G4Facet(3,BOTTOM, -4,BACK_LEFT, -11,BACK_RIGHT, 0,0);
1350  pF[ 7]=G4Facet(4,LEFT_BACK, -8,BACK_TOP, -11,BACK_BOTTOM, 0,0);
1351  pF[ 8]=G4Facet(8,TOP, -7,BACK_RIGHT, -11,BACK_LEFT, 0,0);
1352  pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP, 0,0);
1353 
1354  pF[10]=G4Facet(2,BOTTOM, -3,RIGHT_BACK, -10,RIGHT_FRONT, 0,0);
1355  pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP, -10,RIGHT_BOTTOM, 0,0);
1356  pF[12]=G4Facet(7,TOP, -6,RIGHT_FRONT, -10,RIGHT_BACK, 0,0);
1357  pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP, 0,0);
1358 
1359  pF[14]=G4Facet(1,BOTTOM, -2,FRONT_RIGHT, -9,FRONT_LEFT, 0,0);
1360  pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP, -9,FRONT_BOTTOM, 0,0);
1361  pF[16]=G4Facet(6,TOP, -5,FRONT_LEFT, -9,FRONT_RIGHT, 0,0);
1362  pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP, 0,0);
1363 
1364  pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
1365 
1366  return 0;
1367 }
1368 
1369 G4int
1370 HepPolyhedron::createPolyhedron(G4int Nnodes, G4int Nfaces,
1371  const G4double xyz[][3],
1372  const G4int faces[][4])
1373 /***********************************************************************
1374  * *
1375  * Name: createPolyhedron Date: 05.11.02 *
1376  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1377  * *
1378  * Function: Creates user defined polyhedron *
1379  * *
1380  * Input: Nnodes - number of nodes *
1381  * Nfaces - number of faces *
1382  * nodes[][3] - node coordinates *
1383  * faces[][4] - faces *
1384  * *
1385  ***********************************************************************/
1386 {
1387  AllocateMemory(Nnodes, Nfaces);
1388  if (nvert == 0) return 1;
1389 
1390  for (G4int i=0; i<Nnodes; i++) {
1391  pV[i+1] = G4Point3D(xyz[i][0], xyz[i][1], xyz[i][2]);
1392  }
1393  for (G4int k=0; k<Nfaces; k++) {
1394  pF[k+1] = G4Facet(faces[k][0],0,faces[k][1],0,faces[k][2],0,faces[k][3],0);
1395  }
1396  SetReferences();
1397  return 0;
1398 }
1399 
1400 HepPolyhedronTrd2::HepPolyhedronTrd2(G4double Dx1, G4double Dx2,
1401  G4double Dy1, G4double Dy2,
1402  G4double Dz)
1403 /***********************************************************************
1404  * *
1405  * Name: HepPolyhedronTrd2 Date: 22.07.96 *
1406  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1407  * *
1408  * Function: Create GEANT4 TRD2-trapezoid *
1409  * *
1410  * Input: Dx1 - half-length along X at -Dz 8----7 *
1411  * Dx2 - half-length along X ay +Dz 5----6 ! *
1412  * Dy1 - half-length along Y ay -Dz ! 4-!--3 *
1413  * Dy2 - half-length along Y ay +Dz 1----2 *
1414  * Dz - half-length along Z *
1415  * *
1416  ***********************************************************************/
1417 {
1418  AllocateMemory(8,6);
1419 
1420  pV[1] = G4Point3D(-Dx1,-Dy1,-Dz);
1421  pV[2] = G4Point3D( Dx1,-Dy1,-Dz);
1422  pV[3] = G4Point3D( Dx1, Dy1,-Dz);
1423  pV[4] = G4Point3D(-Dx1, Dy1,-Dz);
1424  pV[5] = G4Point3D(-Dx2,-Dy2, Dz);
1425  pV[6] = G4Point3D( Dx2,-Dy2, Dz);
1426  pV[7] = G4Point3D( Dx2, Dy2, Dz);
1427  pV[8] = G4Point3D(-Dx2, Dy2, Dz);
1428 
1429  CreatePrism();
1430 }
1431 
1432 HepPolyhedronTrd2::~HepPolyhedronTrd2() {}
1433 
1434 HepPolyhedronTrd1::HepPolyhedronTrd1(G4double Dx1, G4double Dx2,
1435  G4double Dy, G4double Dz)
1436  : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1437 
1438 HepPolyhedronTrd1::~HepPolyhedronTrd1() {}
1439 
1440 HepPolyhedronBox::HepPolyhedronBox(G4double Dx, G4double Dy, G4double Dz)
1441  : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1442 
1443 HepPolyhedronBox::~HepPolyhedronBox() {}
1444 
1445 HepPolyhedronTrap::HepPolyhedronTrap(G4double Dz,
1446  G4double Theta,
1447  G4double Phi,
1448  G4double Dy1,
1449  G4double Dx1,
1450  G4double Dx2,
1451  G4double Alp1,
1452  G4double Dy2,
1453  G4double Dx3,
1454  G4double Dx4,
1455  G4double Alp2)
1456 /***********************************************************************
1457  * *
1458  * Name: HepPolyhedronTrap Date: 20.11.96 *
1459  * Author: E.Chernyaev Revised: *
1460  * *
1461  * Function: Create GEANT4 TRAP-trapezoid *
1462  * *
1463  * Input: DZ - half-length in Z *
1464  * Theta,Phi - polar angles of the line joining centres of the *
1465  * faces at Z=-Dz and Z=+Dz *
1466  * Dy1 - half-length in Y of the face at Z=-Dz *
1467  * Dx1 - half-length in X of low edge of the face at Z=-Dz *
1468  * Dx2 - half-length in X of top edge of the face at Z=-Dz *
1469  * Alp1 - angle between Y-axis and the median joining top and *
1470  * low edges of the face at Z=-Dz *
1471  * Dy2 - half-length in Y of the face at Z=+Dz *
1472  * Dx3 - half-length in X of low edge of the face at Z=+Dz *
1473  * Dx4 - half-length in X of top edge of the face at Z=+Dz *
1474  * Alp2 - angle between Y-axis and the median joining top and *
1475  * low edges of the face at Z=+Dz *
1476  * *
1477  ***********************************************************************/
1478 {
1479  G4double DzTthetaCphi = Dz*std::tan(Theta)*std::cos(Phi);
1480  G4double DzTthetaSphi = Dz*std::tan(Theta)*std::sin(Phi);
1481  G4double Dy1Talp1 = Dy1*std::tan(Alp1);
1482  G4double Dy2Talp2 = Dy2*std::tan(Alp2);
1483 
1484  AllocateMemory(8,6);
1485 
1486  pV[1] = G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
1487  pV[2] = G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
1488  pV[3] = G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
1489  pV[4] = G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
1490  pV[5] = G4Point3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
1491  pV[6] = G4Point3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
1492  pV[7] = G4Point3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
1493  pV[8] = G4Point3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
1494 
1495  CreatePrism();
1496 }
1497 
1498 HepPolyhedronTrap::~HepPolyhedronTrap() {}
1499 
1500 HepPolyhedronPara::HepPolyhedronPara(G4double Dx, G4double Dy, G4double Dz,
1501  G4double Alpha, G4double Theta,
1502  G4double Phi)
1503  : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
1504 
1505 HepPolyhedronPara::~HepPolyhedronPara() {}
1506 
1507 HepPolyhedronParaboloid::HepPolyhedronParaboloid(G4double r1,
1508  G4double r2,
1509  G4double dz,
1510  G4double sPhi,
1511  G4double dPhi)
1512 /***********************************************************************
1513  * *
1514  * Name: HepPolyhedronParaboloid Date: 28.06.07 *
1515  * Author: L.Lindroos, T.Nikitina (CERN), July 2007 Revised: 28.06.07 *
1516  * *
1517  * Function: Constructor for paraboloid *
1518  * *
1519  * Input: r1 - inside and outside radiuses at -Dz *
1520  * r2 - inside and outside radiuses at +Dz *
1521  * dz - half length in Z *
1522  * sPhi - starting angle of the segment *
1523  * dPhi - segment range *
1524  * *
1525  ***********************************************************************/
1526 {
1527  static const G4double wholeCircle=twopi;
1528 
1529  // C H E C K I N P U T P A R A M E T E R S
1530 
1531  G4int k = 0;
1532  if (r1 < 0. || r2 <= 0.) k = 1;
1533 
1534  if (dz <= 0.) k += 2;
1535 
1536  G4double phi1, phi2, dphi;
1537 
1538  if(dPhi < 0.)
1539  {
1540  phi2 = sPhi; phi1 = phi2 + dPhi;
1541  }
1542  else if(dPhi == 0.)
1543  {
1544  phi1 = sPhi; phi2 = phi1 + wholeCircle;
1545  }
1546  else
1547  {
1548  phi1 = sPhi; phi2 = phi1 + dPhi;
1549  }
1550  dphi = phi2 - phi1;
1551 
1552  if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1553  if (dphi > wholeCircle) k += 4;
1554 
1555  if (k != 0) {
1556  std::cerr << "HepPolyhedronParaboloid: error in input parameters";
1557  if ((k & 1) != 0) std::cerr << " (radiuses)";
1558  if ((k & 2) != 0) std::cerr << " (half-length)";
1559  if ((k & 4) != 0) std::cerr << " (angles)";
1560  std::cerr << std::endl;
1561  std::cerr << " r1=" << r1;
1562  std::cerr << " r2=" << r2;
1563  std::cerr << " dz=" << dz << " sPhi=" << sPhi << " dPhi=" << dPhi
1564  << std::endl;
1565  return;
1566  }
1567 
1568  // P R E P A R E T W O P O L Y L I N E S
1569 
1570  G4int n = GetNumberOfRotationSteps();
1571  G4double dl = (r2 - r1) / n;
1572  G4double k1 = (r2*r2 - r1*r1) / 2 / dz;
1573  G4double k2 = (r2*r2 + r1*r1) / 2;
1574 
1575  G4double *zz = new G4double[n + 2], *rr = new G4double[n + 2];
1576 
1577  zz[0] = dz;
1578  rr[0] = r2;
1579 
1580  for(G4int i = 1; i < n - 1; i++)
1581  {
1582  rr[i] = rr[i-1] - dl;
1583  zz[i] = (rr[i]*rr[i] - k2) / k1;
1584  if(rr[i] < 0)
1585  {
1586  rr[i] = 0;
1587  zz[i] = 0;
1588  }
1589  }
1590 
1591  zz[n-1] = -dz;
1592  rr[n-1] = r1;
1593 
1594  zz[n] = dz;
1595  rr[n] = 0;
1596 
1597  zz[n+1] = -dz;
1598  rr[n+1] = 0;
1599 
1600  // R O T A T E P O L Y L I N E S
1601 
1602  RotateAroundZ(0, phi1, dphi, n, 2, zz, rr, -1, -1);
1603  SetReferences();
1604 
1605  delete [] zz;
1606  delete [] rr;
1607 }
1608 
1609 HepPolyhedronParaboloid::~HepPolyhedronParaboloid() {}
1610 
1611 HepPolyhedronHype::HepPolyhedronHype(G4double r1,
1612  G4double r2,
1613  G4double sqrtan1,
1614  G4double sqrtan2,
1615  G4double halfZ)
1616 /***********************************************************************
1617  * *
1618  * Name: HepPolyhedronHype Date: 14.04.08 *
1619  * Author: Tatiana Nikitina (CERN) Revised: 14.04.08 *
1620  * *
1621  * Function: Constructor for Hype *
1622  * *
1623  * Input: r1 - inside radius at z=0 *
1624  * r2 - outside radiuses at z=0 *
1625  * sqrtan1 - sqr of tan of Inner Stereo Angle *
1626  * sqrtan2 - sqr of tan of Outer Stereo Angle *
1627  * halfZ - half length in Z *
1628  * *
1629  ***********************************************************************/
1630 {
1631  static const G4double wholeCircle=twopi;
1632 
1633  // C H E C K I N P U T P A R A M E T E R S
1634 
1635  G4int k = 0;
1636  if (r2 < 0. || r1 < 0. ) k = 1;
1637  if (r1 > r2 ) k = 1;
1638  if (r1 == r2) k = 1;
1639 
1640  if (halfZ <= 0.) k += 2;
1641 
1642  if (sqrtan1<0.||sqrtan2<0.) k += 4;
1643 
1644  if (k != 0)
1645  {
1646  std::cerr << "HepPolyhedronHype: error in input parameters";
1647  if ((k & 1) != 0) std::cerr << " (radiuses)";
1648  if ((k & 2) != 0) std::cerr << " (half-length)";
1649  if ((k & 4) != 0) std::cerr << " (angles)";
1650  std::cerr << std::endl;
1651  std::cerr << " r1=" << r1 << " r2=" << r2;
1652  std::cerr << " halfZ=" << halfZ << " sqrTan1=" << sqrtan1
1653  << " sqrTan2=" << sqrtan2
1654  << std::endl;
1655  return;
1656  }
1657 
1658  // P R E P A R E T W O P O L Y L I N E S
1659 
1660  G4int n = GetNumberOfRotationSteps();
1661  G4double dz = 2.*halfZ / n;
1662  G4double k1 = r1*r1;
1663  G4double k2 = r2*r2;
1664 
1665  G4double *zz = new G4double[n+n+1], *rr = new G4double[n+n+1];
1666 
1667  zz[0] = halfZ;
1668  rr[0] = std::sqrt(sqrtan2*halfZ*halfZ+k2);
1669 
1670  for(G4int i = 1; i < n-1; i++)
1671  {
1672  zz[i] = zz[i-1] - dz;
1673  rr[i] =std::sqrt(sqrtan2*zz[i]*zz[i]+k2);
1674  }
1675 
1676  zz[n-1] = -halfZ;
1677  rr[n-1] = rr[0];
1678 
1679  zz[n] = halfZ;
1680  rr[n] = std::sqrt(sqrtan1*halfZ*halfZ+k1);
1681 
1682  for(G4int i = n+1; i < n+n; i++)
1683  {
1684  zz[i] = zz[i-1] - dz;
1685  rr[i] =std::sqrt(sqrtan1*zz[i]*zz[i]+k1);
1686  }
1687  zz[n+n] = -halfZ;
1688  rr[n+n] = rr[n];
1689 
1690  // R O T A T E P O L Y L I N E S
1691 
1692  RotateAroundZ(0, 0., wholeCircle, n, n, zz, rr, -1, -1);
1693  SetReferences();
1694 
1695  delete [] zz;
1696  delete [] rr;
1697 }
1698 
1699 HepPolyhedronHype::~HepPolyhedronHype() {}
1700 
1701 HepPolyhedronCons::HepPolyhedronCons(G4double Rmn1,
1702  G4double Rmx1,
1703  G4double Rmn2,
1704  G4double Rmx2,
1705  G4double Dz,
1706  G4double Phi1,
1707  G4double Dphi)
1708 /***********************************************************************
1709  * *
1710  * Name: HepPolyhedronCons::HepPolyhedronCons Date: 15.12.96 *
1711  * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
1712  * *
1713  * Function: Constructor for CONS, TUBS, CONE, TUBE *
1714  * *
1715  * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
1716  * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
1717  * Dz - half length in Z *
1718  * Phi1 - starting angle of the segment *
1719  * Dphi - segment range *
1720  * *
1721  ***********************************************************************/
1722 {
1723  static const G4double wholeCircle=twopi;
1724 
1725  // C H E C K I N P U T P A R A M E T E R S
1726 
1727  G4int k = 0;
1728  if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
1729  if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
1730  if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
1731 
1732  if (Dz <= 0.) k += 2;
1733 
1734  G4double phi1, phi2, dphi;
1735  if (Dphi < 0.) {
1736  phi2 = Phi1; phi1 = phi2 - Dphi;
1737  }else if (Dphi == 0.) {
1738  phi1 = Phi1; phi2 = phi1 + wholeCircle;
1739  }else{
1740  phi1 = Phi1; phi2 = phi1 + Dphi;
1741  }
1742  dphi = phi2 - phi1;
1743  if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1744  if (dphi > wholeCircle) k += 4;
1745 
1746  if (k != 0) {
1747  std::cerr << "HepPolyhedronCone(s)/Tube(s): error in input parameters";
1748  if ((k & 1) != 0) std::cerr << " (radiuses)";
1749  if ((k & 2) != 0) std::cerr << " (half-length)";
1750  if ((k & 4) != 0) std::cerr << " (angles)";
1751  std::cerr << std::endl;
1752  std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
1753  std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
1754  std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
1755  << std::endl;
1756  return;
1757  }
1758 
1759  // P R E P A R E T W O P O L Y L I N E S
1760 
1761  G4double zz[4], rr[4];
1762  zz[0] = Dz;
1763  zz[1] = -Dz;
1764  zz[2] = Dz;
1765  zz[3] = -Dz;
1766  rr[0] = Rmx2;
1767  rr[1] = Rmx1;
1768  rr[2] = Rmn2;
1769  rr[3] = Rmn1;
1770 
1771  // R O T A T E P O L Y L I N E S
1772 
1773  RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1);
1774  SetReferences();
1775 }
1776 
1777 HepPolyhedronCons::~HepPolyhedronCons() {}
1778 
1779 HepPolyhedronCone::HepPolyhedronCone(G4double Rmn1, G4double Rmx1,
1780  G4double Rmn2, G4double Rmx2,
1781  G4double Dz) :
1782  HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
1783 
1784 HepPolyhedronCone::~HepPolyhedronCone() {}
1785 
1786 HepPolyhedronTubs::HepPolyhedronTubs(G4double Rmin, G4double Rmax,
1787  G4double Dz,
1788  G4double Phi1, G4double Dphi)
1789  : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
1790 
1791 HepPolyhedronTubs::~HepPolyhedronTubs() {}
1792 
1793 HepPolyhedronTube::HepPolyhedronTube (G4double Rmin, G4double Rmax,
1794  G4double Dz)
1795  : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
1796 
1797 HepPolyhedronTube::~HepPolyhedronTube () {}
1798 
1799 HepPolyhedronPgon::HepPolyhedronPgon(G4double phi,
1800  G4double dphi,
1801  G4int npdv,
1802  G4int nz,
1803  const G4double *z,
1804  const G4double *rmin,
1805  const G4double *rmax)
1806 /***********************************************************************
1807  * *
1808  * Name: HepPolyhedronPgon Date: 09.12.96 *
1809  * Author: E.Chernyaev Revised: *
1810  * *
1811  * Function: Constructor of polyhedron for PGON, PCON *
1812  * *
1813  * Input: phi - initial phi *
1814  * dphi - delta phi *
1815  * npdv - number of steps along phi *
1816  * nz - number of z-planes (at least two) *
1817  * z[] - z coordinates of the slices *
1818  * rmin[] - smaller r at the slices *
1819  * rmax[] - bigger r at the slices *
1820  * *
1821  ***********************************************************************/
1822 {
1823  // C H E C K I N P U T P A R A M E T E R S
1824 
1825  if (dphi <= 0. || dphi > twopi) {
1826  std::cerr
1827  << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
1828  << std::endl;
1829  return;
1830  }
1831 
1832  if (nz < 2) {
1833  std::cerr
1834  << "HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
1835  << std::endl;
1836  return;
1837  }
1838 
1839  if (npdv < 0) {
1840  std::cerr
1841  << "HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
1842  << std::endl;
1843  return;
1844  }
1845 
1846  G4int i;
1847  for (i=0; i<nz; i++) {
1848  if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
1849  std::cerr
1850  << "HepPolyhedronPgon: error in radiuses rmin[" << i << "]="
1851  << rmin[i] << " rmax[" << i << "]=" << rmax[i]
1852  << std::endl;
1853  return;
1854  }
1855  }
1856 
1857  // P R E P A R E T W O P O L Y L I N E S
1858 
1859  G4double *zz, *rr;
1860  zz = new G4double[2*nz];
1861  rr = new G4double[2*nz];
1862 
1863  if (z[0] > z[nz-1]) {
1864  for (i=0; i<nz; i++) {
1865  zz[i] = z[i];
1866  rr[i] = rmax[i];
1867  zz[i+nz] = z[i];
1868  rr[i+nz] = rmin[i];
1869  }
1870  }else{
1871  for (i=0; i<nz; i++) {
1872  zz[i] = z[nz-i-1];
1873  rr[i] = rmax[nz-i-1];
1874  zz[i+nz] = z[nz-i-1];
1875  rr[i+nz] = rmin[nz-i-1];
1876  }
1877  }
1878 
1879  // R O T A T E P O L Y L I N E S
1880 
1881  RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1);
1882  SetReferences();
1883 
1884  delete [] zz;
1885  delete [] rr;
1886 }
1887 
1888 HepPolyhedronPgon::~HepPolyhedronPgon() {}
1889 
1890 HepPolyhedronPcon::HepPolyhedronPcon(G4double phi, G4double dphi, G4int nz,
1891  const G4double *z,
1892  const G4double *rmin,
1893  const G4double *rmax)
1894  : HepPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
1895 
1896 HepPolyhedronPcon::~HepPolyhedronPcon() {}
1897 
1898 HepPolyhedronSphere::HepPolyhedronSphere(G4double rmin, G4double rmax,
1899  G4double phi, G4double dphi,
1900  G4double the, G4double dthe)
1901 /***********************************************************************
1902  * *
1903  * Name: HepPolyhedronSphere Date: 11.12.96 *
1904  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1905  * *
1906  * Function: Constructor of polyhedron for SPHERE *
1907  * *
1908  * Input: rmin - internal radius *
1909  * rmax - external radius *
1910  * phi - initial phi *
1911  * dphi - delta phi *
1912  * the - initial theta *
1913  * dthe - delta theta *
1914  * *
1915  ***********************************************************************/
1916 {
1917  // C H E C K I N P U T P A R A M E T E R S
1918 
1919  if (dphi <= 0. || dphi > twopi) {
1920  std::cerr
1921  << "HepPolyhedronSphere: wrong delta phi = " << dphi
1922  << std::endl;
1923  return;
1924  }
1925 
1926  if (the < 0. || the > pi) {
1927  std::cerr
1928  << "HepPolyhedronSphere: wrong theta = " << the
1929  << std::endl;
1930  return;
1931  }
1932 
1933  if (dthe <= 0. || dthe > pi) {
1934  std::cerr
1935  << "HepPolyhedronSphere: wrong delta theta = " << dthe
1936  << std::endl;
1937  return;
1938  }
1939 
1940  if (the+dthe > pi) {
1941  std::cerr
1942  << "HepPolyhedronSphere: wrong theta + delta theta = "
1943  << the << " " << dthe
1944  << std::endl;
1945  return;
1946  }
1947 
1948  if (rmin < 0. || rmin >= rmax) {
1949  std::cerr
1950  << "HepPolyhedronSphere: error in radiuses"
1951  << " rmin=" << rmin << " rmax=" << rmax
1952  << std::endl;
1953  return;
1954  }
1955 
1956  // P R E P A R E T W O P O L Y L I N E S
1957 
1958  G4int nds = (GetNumberOfRotationSteps() + 1) / 2;
1959  G4int np1 = G4int(dthe*nds/pi+.5) + 1;
1960  if (np1 <= 1) np1 = 2;
1961  G4int np2 = rmin < spatialTolerance ? 1 : np1;
1962 
1963  G4double *zz, *rr;
1964  zz = new G4double[np1+np2];
1965  rr = new G4double[np1+np2];
1966 
1967  G4double a = dthe/(np1-1);
1968  G4double cosa, sina;
1969  for (G4int i=0; i<np1; i++) {
1970  cosa = std::cos(the+i*a);
1971  sina = std::sin(the+i*a);
1972  zz[i] = rmax*cosa;
1973  rr[i] = rmax*sina;
1974  if (np2 > 1) {
1975  zz[i+np1] = rmin*cosa;
1976  rr[i+np1] = rmin*sina;
1977  }
1978  }
1979  if (np2 == 1) {
1980  zz[np1] = 0.;
1981  rr[np1] = 0.;
1982  }
1983 
1984  // R O T A T E P O L Y L I N E S
1985 
1986  RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
1987  SetReferences();
1988 
1989  delete [] zz;
1990  delete [] rr;
1991 }
1992 
1993 HepPolyhedronSphere::~HepPolyhedronSphere() {}
1994 
1995 HepPolyhedronTorus::HepPolyhedronTorus(G4double rmin,
1996  G4double rmax,
1997  G4double rtor,
1998  G4double phi,
1999  G4double dphi)
2000 /***********************************************************************
2001  * *
2002  * Name: HepPolyhedronTorus Date: 11.12.96 *
2003  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
2004  * *
2005  * Function: Constructor of polyhedron for TORUS *
2006  * *
2007  * Input: rmin - internal radius *
2008  * rmax - external radius *
2009  * rtor - radius of torus *
2010  * phi - initial phi *
2011  * dphi - delta phi *
2012  * *
2013  ***********************************************************************/
2014 {
2015  // C H E C K I N P U T P A R A M E T E R S
2016 
2017  if (dphi <= 0. || dphi > twopi) {
2018  std::cerr
2019  << "HepPolyhedronTorus: wrong delta phi = " << dphi
2020  << std::endl;
2021  return;
2022  }
2023 
2024  if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2025  std::cerr
2026  << "HepPolyhedronTorus: error in radiuses"
2027  << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
2028  << std::endl;
2029  return;
2030  }
2031 
2032  // P R E P A R E T W O P O L Y L I N E S
2033 
2034  G4int np1 = GetNumberOfRotationSteps();
2035  G4int np2 = rmin < spatialTolerance ? 1 : np1;
2036 
2037  G4double *zz, *rr;
2038  zz = new G4double[np1+np2];
2039  rr = new G4double[np1+np2];
2040 
2041  G4double a = twopi/np1;
2042  G4double cosa, sina;
2043  for (G4int i=0; i<np1; i++) {
2044  cosa = std::cos(i*a);
2045  sina = std::sin(i*a);
2046  zz[i] = rmax*cosa;
2047  rr[i] = rtor+rmax*sina;
2048  if (np2 > 1) {
2049  zz[i+np1] = rmin*cosa;
2050  rr[i+np1] = rtor+rmin*sina;
2051  }
2052  }
2053  if (np2 == 1) {
2054  zz[np1] = 0.;
2055  rr[np1] = rtor;
2056  np2 = -1;
2057  }
2058 
2059  // R O T A T E P O L Y L I N E S
2060 
2061  RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1);
2062  SetReferences();
2063 
2064  delete [] zz;
2065  delete [] rr;
2066 }
2067 
2068 HepPolyhedronTorus::~HepPolyhedronTorus() {}
2069 
2070 HepPolyhedronEllipsoid::HepPolyhedronEllipsoid(G4double ax, G4double by,
2071  G4double cz, G4double zCut1,
2072  G4double zCut2)
2073 /***********************************************************************
2074  * *
2075  * Name: HepPolyhedronEllipsoid Date: 25.02.05 *
2076  * Author: G.Guerrieri Revised: *
2077  * *
2078  * Function: Constructor of polyhedron for ELLIPSOID *
2079  * *
2080  * Input: ax - semiaxis x *
2081  * by - semiaxis y *
2082  * cz - semiaxis z *
2083  * zCut1 - lower cut plane level (solid lies above this plane) *
2084  * zCut2 - upper cut plane level (solid lies below this plane) *
2085  * *
2086  ***********************************************************************/
2087 {
2088  // C H E C K I N P U T P A R A M E T E R S
2089 
2090  if (zCut1 >= cz || zCut2 <= -cz || zCut1 > zCut2) {
2091  std::cerr << "HepPolyhedronEllipsoid: wrong zCut1 = " << zCut1
2092  << " zCut2 = " << zCut2
2093  << " for given cz = " << cz << std::endl;
2094  return;
2095  }
2096  if (cz <= 0.0) {
2097  std::cerr << "HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
2098  << std::endl;
2099  return;
2100  }
2101 
2102  G4double dthe;
2103  G4double sthe;
2104  G4int cutflag;
2105  cutflag= 0;
2106  if (zCut2 >= cz)
2107  {
2108  sthe= 0.0;
2109  }
2110  else
2111  {
2112  sthe= std::acos(zCut2/cz);
2113  cutflag++;
2114  }
2115  if (zCut1 <= -cz)
2116  {
2117  dthe= pi - sthe;
2118  }
2119  else
2120  {
2121  dthe= std::acos(zCut1/cz)-sthe;
2122  cutflag++;
2123  }
2124 
2125  // P R E P A R E T W O P O L Y L I N E S
2126  // generate sphere of radius cz first, then rescale x and y later
2127 
2128  G4int nds = (GetNumberOfRotationSteps() + 1) / 2;
2129  G4int np1 = G4int(dthe*nds/pi) + 2 + cutflag;
2130 
2131  G4double *zz, *rr;
2132  zz = new G4double[np1+1];
2133  rr = new G4double[np1+1];
2134  if (!zz || !rr)
2135  {
2136  G4Exception("HepPolyhedronEllipsoid::HepPolyhedronEllipsoid",
2137  "greps1002", FatalException, "Out of memory");
2138  }
2139 
2140  G4double a = dthe/(np1-cutflag-1);
2141  G4double cosa, sina;
2142  G4int j=0;
2143  if (sthe > 0.0)
2144  {
2145  zz[j]= zCut2;
2146  rr[j]= 0.;
2147  j++;
2148  }
2149  for (G4int i=0; i<np1-cutflag; i++) {
2150  cosa = std::cos(sthe+i*a);
2151  sina = std::sin(sthe+i*a);
2152  zz[j] = cz*cosa;
2153  rr[j] = cz*sina;
2154  j++;
2155  }
2156  if (j < np1)
2157  {
2158  zz[j]= zCut1;
2159  rr[j]= 0.;
2160  j++;
2161  }
2162  if (j > np1)
2163  {
2164  std::cerr << "Logic error in HepPolyhedronEllipsoid, memory corrupted!"
2165  << std::endl;
2166  }
2167  if (j < np1)
2168  {
2169  std::cerr << "Warning: logic error in HepPolyhedronEllipsoid."
2170  << std::endl;
2171  np1= j;
2172  }
2173  zz[j] = 0.;
2174  rr[j] = 0.;
2175 
2176 
2177  // R O T A T E P O L Y L I N E S
2178 
2179  RotateAroundZ(0, 0.0, twopi, np1, 1, zz, rr, -1, 1);
2180  SetReferences();
2181 
2182  delete [] zz;
2183  delete [] rr;
2184 
2185  // rescale x and y vertex coordinates
2186  {
2187  G4Point3D * p= pV;
2188  for (G4int i=0; i<nvert; i++, p++) {
2189  p->setX( p->x() * ax/cz );
2190  p->setY( p->y() * by/cz );
2191  }
2192  }
2193 }
2194 
2195 HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() {}
2196 
2197 HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone(G4double ax,
2198  G4double ay,
2199  G4double h,
2200  G4double zTopCut)
2201 /***********************************************************************
2202  * *
2203  * Name: HepPolyhedronEllipticalCone Date: 8.9.2005 *
2204  * Author: D.Anninos Revised: 9.9.2005 *
2205  * *
2206  * Function: Constructor for EllipticalCone *
2207  * *
2208  * Input: ax, ay - X & Y semi axes at z = 0 *
2209  * h - height of full cone *
2210  * zTopCut - Top Cut in Z Axis *
2211  * *
2212  ***********************************************************************/
2213 {
2214  // C H E C K I N P U T P A R A M E T E R S
2215 
2216  G4int k = 0;
2217  if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
2218 
2219  if (k != 0) {
2220  std::cerr << "HepPolyhedronCone: error in input parameters";
2221  std::cerr << std::endl;
2222  return;
2223  }
2224 
2225  // P R E P A R E T W O P O L Y L I N E S
2226 
2227  zTopCut = (h >= zTopCut ? zTopCut : h);
2228 
2229  G4double *zz, *rr;
2230  zz = new G4double[4];
2231  rr = new G4double[4];
2232  zz[0] = zTopCut;
2233  zz[1] = -zTopCut;
2234  zz[2] = zTopCut;
2235  zz[3] = -zTopCut;
2236  rr[0] = (h-zTopCut);
2237  rr[1] = (h+zTopCut);
2238  rr[2] = 0.;
2239  rr[3] = 0.;
2240 
2241  // R O T A T E P O L Y L I N E S
2242 
2243  RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, -1, -1);
2244  SetReferences();
2245 
2246  delete [] zz;
2247  delete [] rr;
2248 
2249  // rescale x and y vertex coordinates
2250  {
2251  G4Point3D * p= pV;
2252  for (G4int i=0; i<nvert; i++, p++) {
2253  p->setX( p->x() * ax );
2254  p->setY( p->y() * ay );
2255  }
2256  }
2257 }
2258 
2259 HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() {}
2260 
2261 G4ThreadLocal G4int HepPolyhedron::fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
2262 /***********************************************************************
2263  * *
2264  * Name: HepPolyhedron::fNumberOfRotationSteps Date: 24.06.97 *
2265  * Author: J.Allison (Manchester University) Revised: *
2266  * *
2267  * Function: Number of steps for whole circle *
2268  * *
2269  ***********************************************************************/
2270 
2271 #include "BooleanProcessor.src"
2272 
2273 HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const
2274 /***********************************************************************
2275  * *
2276  * Name: HepPolyhedron::add Date: 19.03.00 *
2277  * Author: E.Chernyaev Revised: *
2278  * *
2279  * Function: Boolean "union" of two polyhedra *
2280  * *
2281  ***********************************************************************/
2282 {
2283  G4int ierr;
2284  BooleanProcessor processor;
2285  return processor.execute(OP_UNION, *this, p,ierr);
2286 }
2287 
2288 HepPolyhedron HepPolyhedron::intersect(const HepPolyhedron & p) const
2289 /***********************************************************************
2290  * *
2291  * Name: HepPolyhedron::intersect Date: 19.03.00 *
2292  * Author: E.Chernyaev Revised: *
2293  * *
2294  * Function: Boolean "intersection" of two polyhedra *
2295  * *
2296  ***********************************************************************/
2297 {
2298  G4int ierr;
2299  BooleanProcessor processor;
2300  return processor.execute(OP_INTERSECTION, *this, p,ierr);
2301 }
2302 
2303 HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const
2304 /***********************************************************************
2305  * *
2306  * Name: HepPolyhedron::add Date: 19.03.00 *
2307  * Author: E.Chernyaev Revised: *
2308  * *
2309  * Function: Boolean "subtraction" of "p" from "this" *
2310  * *
2311  ***********************************************************************/
2312 {
2313  G4int ierr;
2314  BooleanProcessor processor;
2315  return processor.execute(OP_SUBTRACTION, *this, p,ierr);
2316 }
2317 
2318 //NOTE : include the code of HepPolyhedronProcessor here
2319 // since there is no BooleanProcessor.h
2320 
2321 #undef INTERSECTION
2322 
2323 #include "HepPolyhedronProcessor.src"
2324 
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:68
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:89
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