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