g4tools  5.4.0
polyhedron.icc
Go to the documentation of this file.
1 // Copyright (C) 2010, Guy Barrand. All rights reserved.
2 // See the file tools.license for terms.
3 
4 
5 #include <cfloat> //G.Barrand : to have DBL_EPSILON on Windows.
6 #include <list>
7 
8 #include "../mnmx"
9 
10 namespace tools {
11 namespace hep {
12 
13 // G.Barrand : introduce iabs to avoid a mess with cmath and some compiler.
14 inline int iabs(int a) {
15  return a < 0 ? -a : a;
16 }
17 
18 //--------------------------------------------------------------------//
19 // JFB: //
20 // polyhedron was polyhedron, retrofitted to Open Inventor //
21 // infrastructure: //
22 //--------------------------------------------------------------------//
23 
24 
25 //
26 // ********************************************************************
27 // * DISCLAIMER *
28 // * *
29 // * The following disclaimer summarizes all the specific disclaimers *
30 // * of contributors to this software. The specific disclaimers,which *
31 // * govern, are listed with their locations in: *
32 // * http://cern.ch/geant4/license *
33 // * *
34 // * Neither the authors of this software system, nor their employing *
35 // * institutes,nor the agencies providing financial support for this *
36 // * work make any representation or warranty, express or implied, *
37 // * regarding this software system or assume any liability for its *
38 // * use. *
39 // * *
40 // * This code implementation is the intellectual property of the *
41 // * GEANT4 collaboration. *
42 // * By copying, distributing or modifying the Program (or any work *
43 // * based on the Program) you indicate your acceptance of this *
44 // * statement, and all its terms. *
45 // ********************************************************************
46 //
47 //
48 //
49 //
50 //
51 // G4 Polyhedron library
52 //
53 // History:
54 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
55 //
56 // 30.09.96 E.Chernyaev
57 // - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
58 // - added GetNextUnitNormal, GetNextEdgeIndeces, GetNextEdge
59 //
60 // 15.12.96 E.Chernyaev
61 // - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
62 // - rewritten G4PolyhedronCons;
63 // - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
64 //
65 // 01.06.97 E.Chernyaev
66 // - modified RotateAroundZ, added SetSideFacets
67 //
68 // 19.03.00 E.Chernyaev
69 // - implemented boolean operations (add, subtract, intersect) on polyhedra;
70 //
71 // 25.05.01 E.Chernyaev
72 // - added GetSurfaceArea() and GetVolume();
73 //
74 
75 
76 /***********************************************************************
77  * *
78  * Name: polyhedron operator << Date: 09.05.96 *
79  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
80  * *
81  * Function: Print contents of G4 polyhedron *
82  * *
83  ***********************************************************************/
84 inline int operator==(const SbFacet& v1, const SbFacet& v2) { //G.Barrand
85  for(int i=0;i<4;i++) {
86  if(v1.edge[i].v != v2.edge[i].v) return false;
87  if(v1.edge[i].f != v2.edge[i].f) return false;
88  }
89  return true;
90 }
91 inline int operator!=(const SbFacet& v1, const SbFacet& v2) { //G.Barrand
92  return !(v1 == v2);
93 }
94 inline std::ostream& operator<<(std::ostream & ostr, const SbFacet & facet) {
95  for (int k=0; k<4; k++) {
96  ostr << " " << facet.edge[k].v << "/" << facet.edge[k].f;
97  }
98  return ostr;
99 }
100 inline std::ostream& operator<<(std::ostream & ostr, const polyhedron & ph) {
101  ostr << std::endl;
102  ostr << "Nverteces=" << ph.nvert << ", Nfacets=" << ph.nface << std::endl;
103  int i;
104  for (i=1; i<=ph.nvert; i++) {
105  ostr << "xyz(" << i << ")="
106  << ph.pV[i].v0() << ' ' << ph.pV[i].v1() << ' ' << ph.pV[i].v2()
107  << std::endl;
108  }
109  for (i=1; i<=ph.nface; i++) {
110  ostr << "face(" << i << ")=" << ph.pF[i] << std::endl;
111  }
112  return ostr;
113 }
114 
115 inline polyhedron::polyhedron(const polyhedron &from)
116 /***********************************************************************
117  * *
118  * Name: polyhedron copy constructor Date: 23.07.96 *
119  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
120  * *
121  ***********************************************************************/
122 {
123 #ifdef TOOLS_MEM
124  mem::increment(s_class().c_str());
125 #endif
126  //m_name = 0; //G.Barrand
127  //if(from.m_name) m_name = new std::string(*from.m_name); //G.Barrand
128 
129  if (from.nvert > 0 && from.nface > 0) {
130  nvert = from.nvert;
131  nface = from.nface;
132  pV = new HVPoint3D[nvert + 1];
133  pF = new SbFacet[nface + 1];
134  int i;
135  for (i=1; i<=nvert; i++) pV[i] = from.pV[i];
136  for (i=1; i<=nface; i++) pF[i] = from.pF[i];
137  }else{
138  nvert = 0; nface = 0; pV = 0; pF = 0;
139  }
140  fNumberOfRotationSteps = from.fNumberOfRotationSteps;
141 }
142 
143 inline int operator==(const polyhedron& v1,const polyhedron& v2) { //G.Barrand
144  return v1.isEqual(v2);
145 }
146 inline int operator!=(const polyhedron& v1,const polyhedron& v2) { //G.Barrand
147  return !(v1 == v2);
148 }
149 
150 inline polyhedron & polyhedron::operator=(const polyhedron &from)
151 /***********************************************************************
152  * *
153  * Name: polyhedron operator = Date: 23.07.96 *
154  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
155  * *
156  * Function: Copy contents of one GEANT4 polyhedron to another *
157  * *
158  ***********************************************************************/
159 {
160  if (this == &from) {
161 // ::printf("debug : polyhedron::operaot=() : on same object\n");
162 // std::cerr << "polyhedron::operaot=() :"
163 // << " on same object !"
164 // << std::endl;
165  return *this;
166  }
167 
168  //delete m_name;m_name = 0; //G.Barrand
169  //if(from.m_name) m_name = new std::string(*from.m_name); //G.Barrand
170 
171  delete [] pV;
172  delete [] pF;
173  if (from.nvert > 0 && from.nface > 0) {
174  nvert = from.nvert;
175  nface = from.nface;
176  pV = new HVPoint3D[nvert + 1];
177  pF = new SbFacet[nface + 1];
178  int i;
179  for (i=1; i<=nvert; i++) pV[i] = from.pV[i];
180  for (i=1; i<=nface; i++) pF[i] = from.pF[i];
181  }else{
182  nvert = 0; nface = 0; pV = 0; pF = 0;
183  }
184  fNumberOfRotationSteps = from.fNumberOfRotationSteps;
185  return *this;
186 }
187 
188 //G.Barrand
189 inline bool polyhedron::isEqual(const polyhedron &from) const {
190  if (this == &from) return true;
191  if(nvert!=from.nvert) return false;
192  if(nface!=from.nface) return false;
193  int i;
194  for (i=1; i<=nvert; i++) {
195  if(pV[i]!=from.pV[i]) return false;
196  }
197  for (i=1; i<=nface; i++) {
198  if(!pF[i].isEqual(from.pF[i])) return false;
199  }
200  return true;
201 }
202 inline bool polyhedron::isConsistent(const char*) const {
203  for(int i=1;i<=nface;i++) {
204  const SbFacet& facet = pF[i];
205  for(int j=0;j<4;j++) {
206  int v,f;
207  facet.GetEdge(j,v,f);
208  if(iabs(v)>nvert) return false;
209  if(iabs(f)>nvert) return false;
210  }
211  }
212  return true;
213 }
214 inline void polyhedron::dump(std::ostream& a_out) const {
215  a_out << " nface = " << nface << std::endl;
216  for(int i=1;i<=nface;i++) {
217  const SbFacet& facet = pF[i];
218  for(int j=0;j<4;j++) {
219  int v,f;
220  facet.GetEdge(j,v,f);
221  a_out << " " << v << " " << f << std::endl;
222  }
223  }
224 }
225 
226 inline int polyhedron::FindNeighbour(int iFace, int iNode, int iOrder) const
227 /***********************************************************************
228  * *
229  * Name: polyhedron::FindNeighbour Date: 22.11.99 *
230  * Author: E.Chernyaev Revised: *
231  * *
232  * Function: Find neighbouring face *
233  * *
234  ***********************************************************************/
235 {
236  int i;
237  for (i=0; i<4; i++) {
238  if (iNode == iabs(pF[iFace].edge[i].v)) break;
239  }
240  if (i == 4) {
241 #ifdef TOOLS_HEP_PH_OUT_ERR
242  std::cerr
243  << "polyhedron::FindNeighbour: face " << iFace
244  << " has no node " << iNode
245  << std::endl;
246 #endif
247  return 0;
248  }
249  if (iOrder < 0) {
250  if ( --i < 0) i = 3;
251  if (pF[iFace].edge[i].v == 0) i = 2;
252  }
253  return (pF[iFace].edge[i].v > 0) ? 0 : pF[iFace].edge[i].f;
254 }
255 
256 inline HVNormal3D polyhedron::FindNodeNormal(int iFace, int iNode) const
257 /***********************************************************************
258  * *
259  * Name: polyhedron::FindNodeNormal Date: 22.11.99 *
260  * Author: E.Chernyaev Revised: *
261  * *
262  * Function: Find normal at given node *
263  * *
264  ***********************************************************************/
265 {
266  HVNormal3D normal = GetUnitNormal(iFace);
267  int k = iFace, iOrder = 1, n = 1;
268 
269  for(;;) {
270  k = FindNeighbour(k, iNode, iOrder);
271  if (k == iFace) break;
272  if (k > 0) {
273  n++;
274  normal += GetUnitNormal(k);
275  }else{
276  if (iOrder < 0) break;
277  k = iFace;
278  iOrder = -iOrder;
279  }
280  }
281  normal.normalize();
282  return normal;
283 }
284 
285 inline void polyhedron::SetNumberOfRotationSteps(int n)
286 /***********************************************************************
287  * *
288  * Name: polyhedron::SetNumberOfRotationSteps Date: 24.06.97 *
289  * Author: J.Allison (Manchester University) Revised: *
290  * *
291  * Function: Set number of steps for whole circle *
292  * *
293  ***********************************************************************/
294 {
295  const int nMin = 3;
296  if (n < nMin) {
297 #ifdef TOOLS_HEP_PH_OUT_ERR
298  std::cerr
299  << "polyhedron::SetNumberOfRotationSteps: attempt to set the\n"
300  << "number of steps per circle < " << nMin << "; forced to " << nMin
301  << std::endl;
302 #endif
303  fNumberOfRotationSteps = nMin;
304  }else{
305  fNumberOfRotationSteps = n;
306  }
307 }
308 
309 inline void polyhedron::AllocateMemory(int Nvert, int Nface)
310 /***********************************************************************
311  * *
312  * Name: polyhedron::AllocateMemory Date: 19.06.96 *
313  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
314  * *
315  * Function: Allocate memory for GEANT4 polyhedron *
316  * *
317  * Input: Nvert - number of nodes *
318  * Nface - number of faces *
319  * *
320  ***********************************************************************/
321 {
322  nvert = Nvert;
323  nface = Nface;
324  pV = new HVPoint3D[nvert+1];
325  pF = new SbFacet[nface+1];
326 }
327 
328 inline void polyhedron::CreatePrism()
329 /***********************************************************************
330  * *
331  * Name: polyhedron::CreatePrism Date: 15.07.96 *
332  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
333  * *
334  * Function: Set facets for a prism *
335  * *
336  ***********************************************************************/
337 {
338  enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
339 
340  pF[1] = SbFacet(1,LEFT, 4,BACK, 3,RIGHT, 2,FRONT);
341  pF[2] = SbFacet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
342  pF[3] = SbFacet(8,TOP, 7,RIGHT, 3,BOTTOM, 4,LEFT);
343  pF[4] = SbFacet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
344  pF[5] = SbFacet(6,TOP, 5,LEFT, 1,BOTTOM, 2,RIGHT);
345  pF[6] = SbFacet(5,FRONT, 6,RIGHT, 7,BACK, 8,LEFT);
346 }
347 
348 inline void polyhedron::RotateEdge(int k1, int k2, double r1, double r2,
349  int v1, int v2, int vEdge,
350  bool ifWholeCircle, int ns, int &kface)
351 /***********************************************************************
352  * *
353  * Name: polyhedron::RotateEdge Date: 05.12.96 *
354  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
355  * *
356  * Function: Create set of facets by rotation of an edge around Z-axis *
357  * *
358  * Input: k1, k2 - end vertices of the edge *
359  * r1, r2 - radiuses of the end vertices *
360  * v1, v2 - visibility of edges produced by rotation of the end *
361  * vertices *
362  * vEdge - visibility of the edge *
363  * ifWholeCircle - is true in case of whole circle rotation *
364  * ns - number of discrete steps *
365  * r[] - r-coordinates *
366  * kface - current free cell in the pF array *
367  * *
368  ***********************************************************************/
369 {
370  if (r1 == 0. && r2 == 0) return;
371 
372  int i;
373  int i1 = k1;
374  int i2 = k2;
375  int ii1 = ifWholeCircle ? i1 : i1+ns;
376  int ii2 = ifWholeCircle ? i2 : i2+ns;
377  int vv = ifWholeCircle ? vEdge : 1;
378 
379  if (ns == 1) {
380  if (r1 == 0.) {
381  pF[kface++] = SbFacet(i1,0, v2*i2,0, (i2+1),0);
382  }else if (r2 == 0.) {
383  pF[kface++] = SbFacet(i1,0, i2,0, v1*(i1+1),0);
384  }else{
385  pF[kface++] = SbFacet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),0);
386  }
387  }else{
388  if (r1 == 0.) {
389  pF[kface++] = SbFacet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0);
390  for (i2++,i=1; i<ns-1; i2++,i++) {
391  pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
392  }
393  pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
394  }else if (r2 == 0.) {
395  pF[kface++] = SbFacet(vv*i1,0, vEdge*i2,0, v1*(i1+1),0);
396  for (i1++,i=1; i<ns-1; i1++,i++) {
397  pF[kface++] = SbFacet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
398  }
399  pF[kface++] = SbFacet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
400  }else{
401  pF[kface++] = SbFacet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
402  for (i1++,i2++,i=1; i<ns-1; i1++,i2++,i++) {
403  pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
404  }
405  pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
406  }
407  }
408 }
409 
410 inline void polyhedron::SetSideFacets(int ii[4], int vv[4],
411  int *kk, double *r,
412  double dphi, int ns, int &kface)
413 /***********************************************************************
414  * *
415  * Name: polyhedron::SetSideFacets Date: 20.05.97 *
416  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
417  * *
418  * Function: Set side facets for the case of incomplete rotation *
419  * *
420  * Input: ii[4] - indeces of original verteces *
421  * vv[4] - visibility of edges *
422  * kk[] - indeces of nodes *
423  * r[] - radiuses *
424  * dphi - delta phi *
425  * ns - number of discrete steps *
426  * kface - current free cell in the pF array *
427  * *
428  ***********************************************************************/
429 {
430  const double perMillion = 0.000001;
431 
432  int k1, k2, k3, k4;
433  if (fabs(dphi-_M_PI()) < perMillion) { // half a circle
434  for (int i=0; i<4; i++) {
435  k1 = ii[i];
436  k2 = (i == 3) ? ii[0] : ii[i+1];
437  if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
438  }
439  }
440 
441  if (ii[1] == ii[2]) {
442  k1 = kk[ii[0]];
443  k2 = kk[ii[2]];
444  k3 = kk[ii[3]];
445  pF[kface++] = SbFacet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
446  if (r[ii[0]] != 0.) k1 += ns;
447  if (r[ii[2]] != 0.) k2 += ns;
448  if (r[ii[3]] != 0.) k3 += ns;
449  pF[kface++] = SbFacet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
450  }else if (kk[ii[0]] == kk[ii[1]]) {
451  k1 = kk[ii[0]];
452  k2 = kk[ii[2]];
453  k3 = kk[ii[3]];
454  pF[kface++] = SbFacet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
455  if (r[ii[0]] != 0.) k1 += ns;
456  if (r[ii[2]] != 0.) k2 += ns;
457  if (r[ii[3]] != 0.) k3 += ns;
458  pF[kface++] = SbFacet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
459  }else if (kk[ii[2]] == kk[ii[3]]) {
460  k1 = kk[ii[0]];
461  k2 = kk[ii[1]];
462  k3 = kk[ii[2]];
463  pF[kface++] = SbFacet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
464  if (r[ii[0]] != 0.) k1 += ns;
465  if (r[ii[1]] != 0.) k2 += ns;
466  if (r[ii[2]] != 0.) k3 += ns;
467  pF[kface++] = SbFacet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
468  }else{
469  k1 = kk[ii[0]];
470  k2 = kk[ii[1]];
471  k3 = kk[ii[2]];
472  k4 = kk[ii[3]];
473  pF[kface++] = SbFacet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
474  if (r[ii[0]] != 0.) k1 += ns;
475  if (r[ii[1]] != 0.) k2 += ns;
476  if (r[ii[2]] != 0.) k3 += ns;
477  if (r[ii[3]] != 0.) k4 += ns;
478  pF[kface++] = SbFacet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
479  }
480 }
481 
482 inline void polyhedron::RotateAroundZ(int nstep, double phi, double dphi,
483  int np1, int np2,
484  const double *z, double *r,
485  int nodeVis, int edgeVis)
486 /***********************************************************************
487  * *
488  * Name: polyhedron::RotateAroundZ Date: 27.11.96 *
489  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
490  * *
491  * Function: Create polyhedron for a solid produced by rotation of *
492  * two polylines around Z-axis *
493  * *
494  * Input: nstep - number of discrete steps, if 0 then default *
495  * phi - starting phi angle *
496  * dphi - delta phi *
497  * np1 - number of points in external polyline *
498  * (must be negative in case of closed polyline) *
499  * np2 - number of points in internal polyline (may be 1) *
500  * z[] - z-coordinates (+z >>> -z for both polylines) *
501  * r[] - r-coordinates *
502  * nodeVis - how to Draw edges joing consecutive positions of *
503  * node during rotation *
504  * edgeVis - how to Draw edges *
505  * *
506  ***********************************************************************/
507 {
508  static const double wholeCircle = 2*_M_PI(); //G.Barrand : const
509 
510  // S E T R O T A T I O N P A R A M E T E R S
511 
512  const double perMillion = 0.000001;
513 
514  bool ifWholeCircle = (fabs(dphi-wholeCircle) < perMillion) ?
515  true : false;
516  double delPhi = ifWholeCircle ? wholeCircle : dphi;
517 
518  int n_step = (nstep > 0) ? nstep : GetNumberOfRotationSteps(); //G.Barrand
519  int nSphi = int(delPhi*n_step/wholeCircle+.5);
520 
521  if (nSphi == 0) nSphi = 1;
522  int nVphi = ifWholeCircle ? nSphi : nSphi+1;
523  bool ifClosed = np1 > 0 ? false : true;
524 
525  // C O U N T V E R T E C E S
526 
527  int absNp1 = iabs(np1);
528  int absNp2 = iabs(np2);
529  int i1beg = 0;
530  int i1end = absNp1-1;
531  int i2beg = absNp1;
532  int i2end = absNp1+absNp2-1;
533  int i, j, k;
534 
535  for(i=i1beg; i<=i2end; i++) {
536  if (fabs(r[i]) < perMillion) r[i] = 0.;
537  }
538 
539  j = 0; // external nodes
540  for (i=i1beg; i<=i1end; i++) {
541  j += (r[i] == 0.) ? 1 : nVphi;
542  }
543 
544  bool ifSide1 = false; // internal nodes
545  bool ifSide2 = false;
546 
547  if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
548  j += (r[i2beg] == 0.) ? 1 : nVphi;
549  ifSide1 = true;
550  }
551 
552  for(i=i2beg+1; i<i2end; i++) {
553  j += (r[i] == 0.) ? 1 : nVphi;
554  }
555 
556  if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
557  if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
558  ifSide2 = true;
559  }
560 
561  // C O U N T F A C E S
562 
563  k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi; // external faces
564 
565  if (absNp2 > 1) { // internal faces
566  for(i=i2beg; i<i2end; i++) {
567  if (r[i] > 0. || r[i+1] > 0.) k += nSphi;
568  }
569 
570  if (ifClosed) {
571  if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
572  }
573  }
574 
575  if (!ifClosed) { // side faces
576  if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
577  if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) k += nSphi;
578  }
579 
580  if (!ifWholeCircle) { // phi_side faces
581  k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
582  }
583 
584  // A L L O C A T E M E M O R Y
585 
586  AllocateMemory(j, k);
587 
588  // G E N E R A T E V E R T E C E S
589 
590  int *kk;
591  kk = new int[absNp1+absNp2];
592 
593  k = 1;
594  for(i=i1beg; i<=i1end; i++) {
595  kk[i] = k;
596  if (r[i] == 0.) { pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
597  }
598 
599  i = i2beg;
600  if (ifSide1) {
601  kk[i] = k;
602  if (r[i] == 0.) { pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
603  }else{
604  kk[i] = kk[i1beg];
605  }
606 
607  for(i=i2beg+1; i<i2end; i++) {
608  kk[i] = k;
609  if (r[i] == 0.) { pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
610  }
611 
612  if (absNp2 > 1) {
613  i = i2end;
614  if (ifSide2) {
615  kk[i] = k;
616  if (r[i] == 0.) pV[k] = HVPoint3D(0, 0, z[i]);
617  }else{
618  kk[i] = kk[i1end];
619  }
620  }
621 
622  double cosPhi, sinPhi;
623 
624  for(j=0; j<nVphi; j++) {
625  cosPhi = cos(phi+j*delPhi/nSphi);
626  sinPhi = sin(phi+j*delPhi/nSphi);
627  for(i=i1beg; i<=i2end; i++) {
628  if (r[i] != 0.) pV[kk[i]+j] = HVPoint3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
629  }
630  }
631 
632  // G E N E R A T E E X T E R N A L F A C E S
633 
634  int v1,v2;
635 
636  k = 1;
637  v2 = ifClosed ? nodeVis : 1;
638  for(i=i1beg; i<i1end; i++) {
639  v1 = v2;
640  if (!ifClosed && i == i1end-1) {
641  v2 = 1;
642  }else{
643  v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
644  }
645  RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
646  edgeVis, ifWholeCircle, nSphi, k);
647  }
648  if (ifClosed) {
649  RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
650  edgeVis, ifWholeCircle, nSphi, k);
651  }
652 
653  // G E N E R A T E I N T E R N A L F A C E S
654 
655  if (absNp2 > 1) {
656  v2 = ifClosed ? nodeVis : 1;
657  for(i=i2beg; i<i2end; i++) {
658  v1 = v2;
659  if (!ifClosed && i==i2end-1) {
660  v2 = 1;
661  }else{
662  v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
663  }
664  RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
665  edgeVis, ifWholeCircle, nSphi, k);
666  }
667  if (ifClosed) {
668  RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
669  edgeVis, ifWholeCircle, nSphi, k);
670  }
671  }
672 
673  // G E N E R A T E S I D E F A C E S
674 
675  if (!ifClosed) {
676  if (ifSide1) {
677  RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
678  -1, ifWholeCircle, nSphi, k);
679  }
680  if (ifSide2) {
681  RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
682  -1, ifWholeCircle, nSphi, k);
683  }
684  }
685 
686  // G E N E R A T E S I D E F A C E S for the case of incomplete circle
687 
688  if (!ifWholeCircle) {
689 
690  int ii[4], vv[4];
691 
692  if (ifClosed) {
693  for (i=i1beg; i<=i1end; i++) {
694  ii[0] = i;
695  ii[3] = (i == i1end) ? i1beg : i+1;
696  ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
697  ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
698  vv[0] = -1;
699  vv[1] = 1;
700  vv[2] = -1;
701  vv[3] = 1;
702  SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
703  }
704  }else{
705  for (i=i1beg; i<i1end; i++) {
706  ii[0] = i;
707  ii[3] = i+1;
708  ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
709  ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
710  vv[0] = (i == i1beg) ? 1 : -1;
711  vv[1] = 1;
712  vv[2] = (i == i1end-1) ? 1 : -1;
713  vv[3] = 1;
714  SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
715  }
716  }
717  }
718 
719  delete [] kk;
720 
721  if (k-1 != nface) {
722 #ifdef TOOLS_HEP_PH_OUT_ERR
723  std::cerr
724  << "Polyhedron::RotateAroundZ: number of generated faces ("
725  << k-1 << ") is not equal to the number of allocated faces ("
726  << nface << ")"
727  << std::endl;
728 #endif
729  }
730 }
731 
732 inline void polyhedron::SetReferences()
733 /***********************************************************************
734  * *
735  * Name: polyhedron::SetReferences Date: 04.12.96 *
736  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
737  * *
738  * Function: For each edge set reference to neighbouring facet *
739  * *
740  ***********************************************************************/
741 {
742  if (nface <= 0) return;
743 
744  struct edgeListMember {
745  edgeListMember *next;
746  int v2;
747  int iface;
748  int iedge;
749  } *edgeList, *freeList, **headList;
750 
751 
752  // A L L O C A T E A N D I N I T I A T E L I S T S
753 
754  edgeList = new edgeListMember[2*nface];
755  headList = new edgeListMember*[nvert];
756 
757  int i;
758 #ifdef TOOLS_HEP_PH_NOT_OPT
759  for (i=0; i<nvert; i++) {
760  headList[i] = 0;
761  }
762  freeList = edgeList;
763  for (i=0; i<2*nface-1; i++) {
764  edgeList[i].next = &edgeList[i+1];
765  }
766  edgeList[2*nface-1].next = 0;
767 #else
768  {edgeListMember** hpos = headList;
769  for (i=0; i<nvert; i++,hpos++) *hpos = 0;
770  freeList = edgeList;
771  int num = 2*nface-1;
772  edgeListMember* epos = edgeList;
773  for (i=0; i<num; i++,epos++) {
774  (*epos).next = epos+1;
775  }
776  (*epos).next = 0;}
777 #endif
778 
779  // L O O P A L O N G E D G E S
780 
781  int iface, iedge, nedge, i1, i2, k1, k2;
782  edgeListMember *prev, *cur;
783 
784 #ifdef TOOLS_HEP_PH_NOT_OPT
785 #else
786  edgeListMember** hpos;
787  SbFacet* pF_iface;
788  SbFacet* pF_cur_iface;
789 #endif
790 
791  for(iface=1; iface<=nface; iface++) {
792 #ifdef TOOLS_HEP_PH_NOT_OPT
793  nedge = (pF[iface].edge[3].v == 0) ? 3 : 4;
794 #else
795  pF_iface = pF+iface;
796  nedge = (pF_iface->edge[3].v == 0) ? 3 : 4;
797 #endif
798  for (iedge=0; iedge<nedge; iedge++) {
799  i1 = iedge;
800  i2 = (iedge < nedge-1) ? iedge+1 : 0;
801 #ifdef TOOLS_HEP_PH_NOT_OPT
802  i1 = iabs(pF[iface].edge[i1].v);
803  i2 = iabs(pF[iface].edge[i2].v);
804  k1 = (i1 < i2) ? i1 : i2; // k1 = ::min(i1,i2);
805  k2 = (i1 > i2) ? i1 : i2; // k2 = ::max(i1,i2);
806 #else
807  i1 = iabs(pF_iface->edge[i1].v);
808  i2 = iabs(pF_iface->edge[i2].v);
809  if(i1<i2) {
810  k1 = i1;
811  k2 = i2;
812  } else {
813  k1 = i2;
814  k2 = i1;
815  }
816 #endif
817  // check head of the List corresponding to k1
818 #ifdef TOOLS_HEP_PH_NOT_OPT
819  cur = headList[k1];
820 #else
821  hpos = headList+k1;
822  cur = *hpos;
823 #endif
824  if (cur == 0) {
825 #ifdef TOOLS_HEP_PH_NOT_OPT
826  headList[k1] = freeList;
827  freeList = freeList->next;
828  cur = headList[k1];
829 #else
830  *hpos = freeList;
831  freeList = freeList->next;
832  cur = *hpos;
833 #endif
834  cur->next = 0;
835  cur->v2 = k2;
836  cur->iface = iface;
837  cur->iedge = iedge;
838  continue;
839  }
840 
841  if (cur->v2 == k2) {
842 #ifdef TOOLS_HEP_PH_NOT_OPT
843  headList[k1] = cur->next;
844 #else
845  *hpos = cur->next;
846 #endif
847  cur->next = freeList;
848  freeList = cur;
849 #ifdef TOOLS_HEP_PH_NOT_OPT
850  pF[iface].edge[iedge].f = cur->iface;
851  pF[cur->iface].edge[cur->iedge].f = iface;
852  i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
853  i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
854 #else
855  pF_iface->edge[iedge].f = cur->iface;
856  pF_cur_iface = pF+cur->iface;
857  pF_cur_iface->edge[cur->iedge].f = iface;
858  i1 = (pF_iface->edge[iedge].v < 0) ? -1 : 1;
859  i2 = (pF_cur_iface->edge[cur->iedge].v < 0) ? -1 : 1;
860 #endif
861  if (i1 != i2) {
862 #ifdef TOOLS_HEP_PH_OUT_ERR
863  std::cerr
864  << "Polyhedron::SetReferences: different edge visibility "
865  << iface << "/" << iedge << "/"
866  << pF[iface].edge[iedge].v << " and "
867  << cur->iface << "/" << cur->iedge << "/"
868  << pF[cur->iface].edge[cur->iedge].v
869  << std::endl;
870 #endif
871  }
872  continue;
873  }
874 
875  // check List itself
876  for (;;) {
877  prev = cur;
878  cur = prev->next;
879  if (cur == 0) {
880  prev->next = freeList;
881  freeList = freeList->next;
882  cur = prev->next;
883  cur->next = 0;
884  cur->v2 = k2;
885  cur->iface = iface;
886  cur->iedge = iedge;
887  break;
888  }
889 
890  if (cur->v2 == k2) {
891  prev->next = cur->next;
892  cur->next = freeList;
893  freeList = cur;
894 #ifdef TOOLS_HEP_PH_NOT_OPT
895  pF[iface].edge[iedge].f = cur->iface;
896  pF[cur->iface].edge[cur->iedge].f = iface;
897  i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
898  i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
899 #else
900  pF_iface->edge[iedge].f = cur->iface;
901  pF_cur_iface = pF+cur->iface;
902  pF_cur_iface->edge[cur->iedge].f = iface;
903  i1 = (pF_iface->edge[iedge].v < 0) ? -1 : 1;
904  i2 = (pF_cur_iface->edge[cur->iedge].v < 0) ? -1 : 1;
905 #endif
906  if (i1 != i2) {
907 #ifdef TOOLS_HEP_PH_OUT_ERR
908  std::cerr
909  << "Polyhedron::SetReferences: different edge visibility "
910  << iface << "/" << iedge << "/"
911  << pF[iface].edge[iedge].v << " and "
912  << cur->iface << "/" << cur->iedge << "/"
913  << pF[cur->iface].edge[cur->iedge].v
914  << std::endl;
915 #endif
916  }
917  break;
918  }
919  }
920  }
921  }
922 
923  // 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
924 
925 #ifdef TOOLS_HEP_PH_OUT_ERR
926 #ifdef TOOLS_HEP_PH_NOT_OPT
927  for (i=0; i<nvert; i++) {
928  if (headList[i] != 0) {
929  std::cerr
930  << "Polyhedron::SetReferences: List " << i << " is not empty"
931  << std::endl;
932  }
933  }
934 #else
935  {edgeListMember** hpos = headList;
936  for (i=0; i<nvert; i++,hpos++) {
937  if (*hpos != 0) {
938  std::cerr
939  << "Polyhedron::SetReferences: List " << i << " is not empty"
940  << std::endl;
941  }
942  }}
943 #endif
944 #endif
945 
946  // F R E E M E M O R Y
947 
948  delete [] edgeList;
949  delete [] headList;
950 }
951 
952 inline void polyhedron::InvertFacets()
953 /***********************************************************************
954  * *
955  * Name: polyhedron::InvertFacets Date: 01.12.99 *
956  * Author: E.Chernyaev Revised: *
957  * *
958  * Function: Invert the order of the nodes in the facets *
959  * *
960  ***********************************************************************/
961 {
962  if (nface <= 0) return;
963  int i, k, nnode, v[4],f[4];
964  for (i=1; i<=nface; i++) {
965  nnode = (pF[i].edge[3].v == 0) ? 3 : 4;
966  for (k=0; k<nnode; k++) {
967  v[k] = (k+1 == nnode) ? pF[i].edge[0].v : pF[i].edge[k+1].v;
968  if (v[k] * pF[i].edge[k].v < 0) v[k] = -v[k];
969  f[k] = pF[i].edge[k].f;
970  }
971  for (k=0; k<nnode; k++) {
972  pF[i].edge[nnode-1-k].v = v[k];
973  pF[i].edge[nnode-1-k].f = f[k];
974  }
975  }
976 }
977 
978 /*
979 inline polyhedron & polyhedron::Transform(
980  const SbRotation& rotation
981 ,const SbVec3f& translation
982 )
983 {
984  if (nvert > 0) {
985  for (int i=1; i<=nvert; i++) {
986  const HVPoint3D& pv = pV[i];
987  SbVec3f tmp;
988  rotation.multVec(SbVec3f(pv[0],pv[1],pv[2]),tmp);
989  pV[i] = HVPoint3D(tmp[0],tmp[1],tmp[2])
990  +HVPoint3D(translation[0],translation[1],translation[2]);
991  }
992 
993  // C H E C K D E T E R M I N A N T A N D
994  // 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
995 
996  //FIXME : have the below done in double.
997  SbVec3f x; rotation.multVec(SbVec3f(1,0,0),x);
998  SbVec3f y; rotation.multVec(SbVec3f(0,1,0),y);
999  SbVec3f z; rotation.multVec(SbVec3f(0,0,1),z);
1000  if ((x.cross(y)).dot(z) < 0) InvertFacets();
1001  }
1002  return *this;
1003 }
1004 */
1005 
1006 /*G.Barrand : optimized version.*/
1007 inline polyhedron & polyhedron::Translate(
1008  double a_x
1009 ,double a_y
1010 ,double a_z
1011 )
1012 {
1013  if (nvert > 0) {
1014  for (int i=1; i<=nvert; i++) {
1015  const HVPoint3D& p = pV[i];
1016  pV[i].set_value(p.x()+a_x,p.y()+a_y,p.z()+a_z);
1017  }
1018  }
1019  return *this;
1020 }
1021 
1022 inline polyhedron & polyhedron::Transform(const rotd& rotation,double a_x,double a_y,double a_z) {
1023  if (nvert > 0) {
1024  {double x,y,z;
1025  for (int i=1; i<=nvert; i++) {
1026  HVPoint3D& p = pV[i];
1027  x = p.x();
1028  y = p.y();
1029  z = p.z();
1030  rotation.mul_3(x,y,z); //tmp = R*pV[i]
1031  p.set_value(x+a_x,y+a_y,z+a_z);
1032  }}
1033 
1034  double x_x = 1;
1035  double x_y = 0;
1036  double x_z = 0;
1037  rotation.mul_3(x_x,x_y,x_z);
1038 
1039  double y_x = 0;
1040  double y_y = 1;
1041  double y_z = 0;
1042  rotation.mul_3(y_x,y_y,y_z);
1043 
1044  double z_x = 0;
1045  double z_y = 0;
1046  double z_z = 1;
1047  rotation.mul_3(z_x,z_y,z_z);
1048 
1049  // x_x y_x
1050  // x_y y_y
1051  // x_z y_z
1052  double x_cross_y_x = x_y*y_z-x_z*y_y;
1053  double x_cross_y_y = x_z*y_x-x_x*y_z;
1054  double x_cross_y_z = x_x*y_y-x_y*y_x;
1055 
1056  double x_cross_y_dot_z =
1057  x_cross_y_x*z_x + x_cross_y_y*z_y + x_cross_y_z*z_z;
1058  if (x_cross_y_dot_z < 0) InvertFacets();
1059  }
1060  return *this;
1061 }
1062 
1063 inline polyhedron & polyhedron::Transform(const rotd& rotation,const vec3d& translation) {
1064 /*
1065  if (nvert > 0) {
1066  vec3d tmp;
1067  for (int i=1; i<=nvert; i++) {
1068  rotation.mul_vec(pV[i],tmp); //tmp = R*pV[i]
1069  pV[i] = tmp+translation;
1070  }
1071  vec3d x; rotation.mul_vec(vec3d::s_x(),x);
1072  vec3d y; rotation.mul_vec(vec3d::s_y(),y);
1073  vec3d z; rotation.mul_vec(vec3d::s_z(),z);
1074  if ((x.cross(y)).dot(z) < 0) InvertFacets();
1075  }
1076 */
1077  return Transform(rotation,translation.x(),translation.y(),translation.z());
1078 }
1079 
1080 /*G.Barrand : inline
1081 bool polyhedron::GetNextVertexIndex(int &index, int &edgeFlag) const
1082 // ***********************************************************************
1083 // * *
1084 // * Name: polyhedron::GetNextVertexIndex Date: 03.09.96 *
1085 // * Author: Yasuhide Sawada Revised: *
1086 // * *
1087 // * Function: *
1088 // * *
1089 // ***********************************************************************
1090 {
1091  static int iFace = 1;
1092  static int iQVertex = 0;
1093  int vIndex = pF[iFace].edge[iQVertex].v;
1094 
1095  edgeFlag = (vIndex > 0) ? 1 : 0;
1096  index = iabs(vIndex);
1097 
1098  if(index>nvert) {
1099 #ifdef TOOLS_HEP_PH_OUT_ERR
1100  std::cerr << "polyhedron::GetNextVertexIndex: pV index problem "
1101  << index << " exceed " << nvert << std::endl;
1102 #endif
1103  index = 0;
1104  }
1105 
1106  if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1107  iQVertex = 0;
1108  if (++iFace > nface) iFace = 1;
1109  return false; // Last Edge
1110  }else{
1111  ++iQVertex;
1112  return true; // not Last Edge
1113  }
1114 }
1115 */
1116 
1117 inline HVPoint3D polyhedron::GetVertex(int index) const
1118 /***********************************************************************
1119  * *
1120  * Name: polyhedron::GetVertex Date: 03.09.96 *
1121  * Author: Yasuhide Sawada Revised: 17.11.99 *
1122  * *
1123  * Function: Get vertex of the index. *
1124  * *
1125  ***********************************************************************/
1126 {
1127  if (index <= 0 || index > nvert) {
1128 #ifdef TOOLS_HEP_PH_OUT_ERR
1129  std::cerr
1130  << "polyhedron::GetVertex: irrelevant index " << index
1131  << std::endl;
1132 #endif
1133  return HVPoint3D();
1134  }
1135  return pV[index];
1136 }
1137 
1138 inline const HVPoint3D& polyhedron::GetVertexFast(int index) const {//G.Barrand
1139  return pV[index];
1140 }
1141 
1142 inline bool polyhedron::GetNextVertex(HVPoint3D &vertex, int &edgeFlag) const
1143 /***********************************************************************
1144  * *
1145  * Name: polyhedron::GetNextVertex Date: 22.07.96 *
1146  * Author: John Allison Revised: *
1147  * *
1148  * Function: Get vertices of the quadrilaterals in order for each *
1149  * face in face order. Returns false when finished each *
1150  * face. *
1151  * *
1152  ***********************************************************************/
1153 {
1154  int index;
1155  bool rep = GetNextVertexIndex(index, edgeFlag);
1156  vertex = pV[index];
1157  return rep;
1158 }
1159 
1160 inline bool polyhedron::GetNextVertex(HVPoint3D &vertex, int &edgeFlag,HVNormal3D &normal) const
1161 /***********************************************************************
1162  * *
1163  * Name: polyhedron::GetNextVertex Date: 26.11.99 *
1164  * Author: E.Chernyaev Revised: *
1165  * *
1166  * Function: Get vertices with normals of the quadrilaterals in order *
1167  * for each face in face order. *
1168  * Returns false when finished each face. *
1169  * *
1170  ***********************************************************************/
1171 {
1172  static int iFace = 1;
1173  static int iNode = 0;
1174 
1175  if (nface == 0) return false; // empty polyhedron
1176 
1177  int k = pF[iFace].edge[iNode].v;
1178  if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
1179  vertex = pV[k];
1180  normal = FindNodeNormal(iFace,k);
1181  if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
1182  iNode = 0;
1183  if (++iFace > nface) iFace = 1;
1184  return false; // last node
1185  }else{
1186  ++iNode;
1187  return true; // not last node
1188  }
1189 }
1190 
1191 inline bool polyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag,
1192  int &iface1, int &iface2) const
1193 /***********************************************************************
1194  * *
1195  * Name: polyhedron::GetNextEdgeIndeces Date: 30.09.96 *
1196  * Author: E.Chernyaev Revised: 17.11.99 *
1197  * *
1198  * Function: Get indeces of the next edge together with indeces of *
1199  * of the faces which share the edge. *
1200  * Returns false when the last edge. *
1201  * *
1202  ***********************************************************************/
1203 {
1204  static int iFace = 1;
1205  static int iQVertex = 0;
1206  static int iOrder = 1;
1207  int k1, k2, kflag, kface1, kface2;
1208 
1209  if (iFace == 1 && iQVertex == 0) {
1210  k2 = pF[nface].edge[0].v;
1211  k1 = pF[nface].edge[3].v;
1212  if (k1 == 0) k1 = pF[nface].edge[2].v;
1213  if (iabs(k1) > iabs(k2)) iOrder = -1;
1214  }
1215 
1216  do {
1217  k1 = pF[iFace].edge[iQVertex].v;
1218  kflag = k1;
1219  k1 = iabs(k1);
1220  kface1 = iFace;
1221  kface2 = pF[iFace].edge[iQVertex].f;
1222  if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1223  iQVertex = 0;
1224  k2 = iabs(pF[iFace].edge[iQVertex].v);
1225  iFace++;
1226  }else{
1227  iQVertex++;
1228  k2 = iabs(pF[iFace].edge[iQVertex].v);
1229  }
1230  } while (iOrder*k1 > iOrder*k2);
1231 
1232  i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1233  iface1 = kface1; iface2 = kface2;
1234 
1235  if (iFace > nface) {
1236  iFace = 1; iOrder = 1;
1237  return false;
1238  }else{
1239  return true;
1240  }
1241 }
1242 
1243 inline bool polyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag) const
1244 /***********************************************************************
1245  * *
1246  * Name: polyhedron::GetNextEdgeIndeces Date: 17.11.99 *
1247  * Author: E.Chernyaev Revised: *
1248  * *
1249  * Function: Get indeces of the next edge. *
1250  * Returns false when the last edge. *
1251  * *
1252  ***********************************************************************/
1253 {
1254  int kface1, kface2;
1255  return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
1256 }
1257 
1258 inline bool polyhedron::GetNextEdge(HVPoint3D &p1,HVPoint3D &p2,int &edgeFlag) const
1259 /***********************************************************************
1260  * *
1261  * Name: polyhedron::GetNextEdge Date: 30.09.96 *
1262  * Author: E.Chernyaev Revised: *
1263  * *
1264  * Function: Get next edge. *
1265  * Returns false when the last edge. *
1266  * *
1267  ***********************************************************************/
1268 {
1269  int i1,i2;
1270  bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
1271  p1 = pV[i1];
1272  p2 = pV[i2];
1273  return rep;
1274 }
1275 
1276 inline bool polyhedron::GetNextEdge(HVPoint3D &p1, HVPoint3D &p2,int &edgeFlag, int &iface1, int &iface2) const
1277 /***********************************************************************
1278  * *
1279  * Name: polyhedron::GetNextEdge Date: 17.11.99 *
1280  * Author: E.Chernyaev Revised: *
1281  * *
1282  * Function: Get next edge with indeces of the faces which share *
1283  * the edge. *
1284  * Returns false when the last edge. *
1285  * *
1286  ***********************************************************************/
1287 {
1288  int i1,i2;
1289  bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
1290  p1 = pV[i1];
1291  p2 = pV[i2];
1292  return rep;
1293 }
1294 
1295 inline void polyhedron::GetFacet(int iFace, int &n, int *iNodes,int *edgeFlags, int *iFaces) const
1296 /***********************************************************************
1297  * *
1298  * Name: polyhedron::GetFacet Date: 15.12.99 *
1299  * Author: E.Chernyaev Revised: *
1300  * *
1301  * Function: Get face by index *
1302  * *
1303  ***********************************************************************/
1304 {
1305  if (iFace < 1 || iFace > nface) {
1306 #ifdef TOOLS_HEP_PH_OUT_ERR
1307  std::cerr
1308  << "polyhedron::GetFacet: irrelevant index " << iFace
1309  << std::endl;
1310 #endif
1311  n = 0;
1312  }else{
1313  int i, k;
1314  for (i=0; i<4; i++) {
1315  k = pF[iFace].edge[i].v;
1316  if (k == 0) break;
1317  if (iFaces != 0) iFaces[i] = pF[iFace].edge[i].f;
1318  if (k > 0) {
1319  iNodes[i] = k;
1320  if (edgeFlags != 0) edgeFlags[i] = 1;
1321  }else{
1322  iNodes[i] = -k;
1323  if (edgeFlags != 0) edgeFlags[i] = -1;
1324  }
1325  }
1326  n = i;
1327  }
1328 }
1329 
1330 inline void polyhedron::GetFacet(int index, int &n, HVPoint3D *nodes,int *edgeFlags, HVNormal3D *normals) const
1331 /***********************************************************************
1332  * *
1333  * Name: polyhedron::GetFacet Date: 17.11.99 *
1334  * Author: E.Chernyaev Revised: *
1335  * *
1336  * Function: Get face by index *
1337  * *
1338  ***********************************************************************/
1339 {
1340  int iNodes[4];
1341  GetFacet(index, n, iNodes, edgeFlags);
1342  if (n != 0) {
1343  for (int i=0; i<4; i++) {
1344  nodes[i] = pV[iNodes[i]];
1345  if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1346  }
1347  }
1348 }
1349 
1350 inline bool polyhedron::GetNextFacet(int &n, HVPoint3D *nodes,int *edgeFlags, HVNormal3D *normals) const
1351 /***********************************************************************
1352  * *
1353  * Name: polyhedron::GetNextFacet Date: 19.11.99 *
1354  * Author: E.Chernyaev Revised: *
1355  * *
1356  * Function: Get next face with normals of unit length at the nodes. *
1357  * Returns false when finished all faces. *
1358  * *
1359  ***********************************************************************/
1360 {
1361  static int iFace = 1;
1362 
1363  if (edgeFlags == 0) {
1364  GetFacet(iFace, n, nodes);
1365  }else if (normals == 0) {
1366  GetFacet(iFace, n, nodes, edgeFlags);
1367  }else{
1368  GetFacet(iFace, n, nodes, edgeFlags, normals);
1369  }
1370 
1371  if (++iFace > nface) {
1372  iFace = 1;
1373  return false;
1374  }else{
1375  return true;
1376  }
1377 }
1378 
1379 //G.Barrand
1380 #ifdef TOOLS_HEP_PH_OUT_ERR
1381 inline bool polyhedron::CHECK_INDEX(const char* a_method,int a_index) const {
1382  if(a_index>nvert) {
1383  std::cerr << "polyhedron::" << a_method << " :"
1384  << " index problem. "
1385  << a_index << " exceed " << nvert << std::endl;
1386  return false;
1387  }
1388  return true;
1389 }
1390 #else
1391 inline bool polyhedron::CHECK_INDEX(const char*,int a_index) const {
1392  if(a_index>nvert) {
1393  return false;
1394  }
1395  return true;
1396 }
1397 #endif
1398 
1399 inline HVNormal3D polyhedron::GetNormal(int iFace) const
1400 /***********************************************************************
1401  * *
1402  * Name: polyhedron::GetNormal Date: 19.11.99 *
1403  * Author: E.Chernyaev Revised: *
1404  * *
1405  * Function: Get normal of the face given by index *
1406  * *
1407  ***********************************************************************/
1408 {
1409  if (iFace < 1 || iFace > nface) {
1410 #ifdef TOOLS_HEP_PH_OUT_ERR
1411  std::cerr
1412  << "polyhedron::GetNormal: irrelevant index " << iFace
1413  << std::endl;
1414 #endif
1415  return HVNormal3D();
1416  }
1417 
1418  int i0 = iabs(pF[iFace].edge[0].v);
1419  int i1 = iabs(pF[iFace].edge[1].v);
1420  int i2 = iabs(pF[iFace].edge[2].v);
1421  int i3 = iabs(pF[iFace].edge[3].v);
1422  if (i3 == 0) i3 = i0;
1423 
1424  if(!CHECK_INDEX("GetNormal",i0)) return HVNormal3D();
1425  if(!CHECK_INDEX("GetNormal",i1)) return HVNormal3D();
1426  if(!CHECK_INDEX("GetNormal",i2)) return HVNormal3D();
1427  if(!CHECK_INDEX("GetNormal",i3)) return HVNormal3D();
1428 
1429  HVNormal3D nm;
1430  (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1],nm);
1431  return nm;
1432 }
1433 
1434 inline HVNormal3D polyhedron::GetUnitNormal(int iFace) const
1435 /***********************************************************************
1436  * *
1437  * Name: polyhedron::GetNormal Date: 19.11.99 *
1438  * Author: E.Chernyaev Revised: *
1439  * *
1440  * Function: Get unit normal of the face given by index *
1441  * *
1442  ***********************************************************************/
1443 {
1444  if (iFace < 1 || iFace > nface) {
1445 #ifdef TOOLS_HEP_PH_OUT_ERR
1446  std::cerr
1447  << "polyhedron::GetUnitNormal: irrelevant index " << iFace
1448  << std::endl;
1449 #endif
1450  return HVNormal3D();
1451  }
1452 
1453  int i0 = iabs(pF[iFace].edge[0].v);
1454  int i1 = iabs(pF[iFace].edge[1].v);
1455  int i2 = iabs(pF[iFace].edge[2].v);
1456  int i3 = iabs(pF[iFace].edge[3].v);
1457  if (i3 == 0) i3 = i0;
1458 
1459  if(!CHECK_INDEX("GetUnitNormal",i0)) return HVNormal3D();
1460  if(!CHECK_INDEX("GetUnitNormal",i1)) return HVNormal3D();
1461  if(!CHECK_INDEX("GetUnitNormal",i2)) return HVNormal3D();
1462  if(!CHECK_INDEX("GetUnitNormal",i3)) return HVNormal3D();
1463 
1464  HVNormal3D nm;
1465  (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1],nm);
1466  nm.normalize();
1467  return nm;
1468 }
1469 
1470 inline bool polyhedron::GetNextNormal(HVNormal3D &normal) const
1471 /***********************************************************************
1472  * *
1473  * Name: polyhedron::GetNextNormal Date: 22.07.96 *
1474  * Author: John Allison Revised: 19.11.99 *
1475  * *
1476  * Function: Get normals of each face in face order. Returns false *
1477  * when finished all faces. *
1478  * *
1479  ***********************************************************************/
1480 {
1481  static int iFace = 1;
1482  normal = GetNormal(iFace);
1483  if (++iFace > nface) {
1484  iFace = 1;
1485  return false;
1486  }else{
1487  return true;
1488  }
1489 }
1490 
1491 inline bool polyhedron::GetNextUnitNormal(HVNormal3D &normal) const
1492 /***********************************************************************
1493  * *
1494  * Name: polyhedron::GetNextUnitNormal Date: 16.09.96 *
1495  * Author: E.Chernyaev Revised: *
1496  * *
1497  * Function: Get normals of unit length of each face in face order. *
1498  * Returns false when finished all faces. *
1499  * *
1500  ***********************************************************************/
1501 {
1502  bool rep = GetNextNormal(normal);
1503  normal.normalize();
1504  return rep;
1505 }
1506 
1507 inline double polyhedron::GetSurfaceArea() const
1508 /***********************************************************************
1509  * *
1510  * Name: polyhedron::GetSurfaceArea Date: 25.05.01 *
1511  * Author: E.Chernyaev Revised: *
1512  * *
1513  * Function: Returns area of the surface of the polyhedron. *
1514  * *
1515  ***********************************************************************/
1516 {
1517  double s = 0.;
1518  HVPoint3D p;
1519  for (int iFace=1; iFace<=nface; iFace++) {
1520  int i0 = iabs(pF[iFace].edge[0].v);
1521  int i1 = iabs(pF[iFace].edge[1].v);
1522  int i2 = iabs(pF[iFace].edge[2].v);
1523  int i3 = iabs(pF[iFace].edge[3].v);
1524  if (i3 == 0) i3 = i0;
1525  (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1],p);
1526  s += p.length();
1527  }
1528  return s/2.;
1529 }
1530 
1531 inline double polyhedron::GetVolume() const
1532 /***********************************************************************
1533  * *
1534  * Name: polyhedron::GetVolume Date: 25.05.01 *
1535  * Author: E.Chernyaev Revised: *
1536  * *
1537  * Function: Returns volume of the polyhedron. *
1538  * *
1539  ***********************************************************************/
1540 {
1541  double v = 0.;
1542  HVPoint3D p;
1543  for (int iFace=1; iFace<=nface; iFace++) {
1544  int i0 = iabs(pF[iFace].edge[0].v);
1545  int i1 = iabs(pF[iFace].edge[1].v);
1546  int i2 = iabs(pF[iFace].edge[2].v);
1547  int i3 = iabs(pF[iFace].edge[3].v);
1548  HVPoint3D g;
1549  if (i3 == 0) {
1550  i3 = i0;
1551  g = (pV[i0]+pV[i1]+pV[i2]) * (1.0/3.0);
1552  }else{
1553  g = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0.25;
1554  }
1555  (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1],p);
1556  v += p.dot(g);
1557  }
1558  return v/6.;
1559 }
1560 
1561 //G.Barrand
1562 inline
1563 bool polyhedron::set_polyhedron_trd2(double Dx1, double Dx2,
1564  double Dy1, double Dy2,
1565  double Dz)
1566 /***********************************************************************
1567  * *
1568  * Name: polyhedron_trd2 Date: 22.07.96 *
1569  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1570  * *
1571  * Function: Create GEANT4 TRD2-trapezoid *
1572  * *
1573  * Input: Dx1 - half-length along X at -Dz 8----7 *
1574  * Dx2 - half-length along X ay +Dz 5----6 ! *
1575  * Dy1 - half-length along Y ay -Dz ! 4-!--3 *
1576  * Dy2 - half-length along Y ay +Dz 1----2 *
1577  * Dz - half-length along Z *
1578  * *
1579  ***********************************************************************/
1580 {
1581  _clear(); //G.Barrand
1582 
1583  //From TGeoTrd2::Capacity() :
1584  double capacity = 2*(Dx1+Dx2)*(Dy1+Dy2)*Dz +
1585  (2./3.)*(Dx1-Dx2)*(Dy1-Dy2)*Dz;
1586 
1587  //if(//(Dx1<=0.0)||(Dx2<=0.0)||
1588  // //(Dy1<=0.0)||(Dy2<=0.0)||
1589  // (Dz<=0.0)){ //G.Barrand
1590  if(capacity<=0) {
1591 #if defined(TOOLS_HEP_PH_OUT_ERR) || defined(TOOLS_HEP_PH_OUT_ERR_TRD2)
1592  std::cerr << "set_polyhedron_trd2: error in input parameters";
1593  std::cerr << " Dx1=" << Dx1
1594  << " Dx2=" << Dx2
1595  << " Dy1=" << Dy1
1596  << " Dy2=" << Dy2
1597  << " Dz=" << Dz
1598  << std::endl;
1599 #endif
1600  return false;
1601  }
1602 
1603  AllocateMemory(8,6);
1604 
1605  pV[1] = HVPoint3D(-Dx1,-Dy1,-Dz);
1606  pV[2] = HVPoint3D( Dx1,-Dy1,-Dz);
1607  pV[3] = HVPoint3D( Dx1, Dy1,-Dz);
1608  pV[4] = HVPoint3D(-Dx1, Dy1,-Dz);
1609  pV[5] = HVPoint3D(-Dx2,-Dy2, Dz);
1610  pV[6] = HVPoint3D( Dx2,-Dy2, Dz);
1611  pV[7] = HVPoint3D( Dx2, Dy2, Dz);
1612  pV[8] = HVPoint3D(-Dx2, Dy2, Dz);
1613 
1614  CreatePrism();
1615 
1616  return true;
1617 }
1618 
1619 inline
1620 polyhedron_trd2::polyhedron_trd2(double Dx1, double Dx2,
1621  double Dy1, double Dy2,
1622  double Dz)
1623 {
1624  set_polyhedron_trd2(Dx1,Dx2,Dy1,Dy2,Dz); //G.Barrand
1625 }
1626 
1627 //G.Barrand
1628 inline
1629 bool polyhedron::set_polyhedron_arb8(double Dz,const double* xy) {
1630  _clear(); //G.Barrand
1631 
1632  // xy as if xy[8][2]
1633  // then xy[i][j] = xy[i*2+j]
1634 
1635  // from TGeoArb8::Capacity() :
1636  double capacity = 0;
1637  {int j;
1638  for(int i=0; i<4; i++) {
1639  j = (i+1)%4;
1640  capacity +=
1641  0.25*Dz*((vxy(xy,i,0)+vxy(xy,i+4,0))*(vxy(xy,j,1)+vxy(xy,j+4,1)) -
1642  (vxy(xy,j,0)+vxy(xy,j+4,0))*(vxy(xy,i,1)+vxy(xy,i+4,1)) +
1643  (1./3)*((vxy(xy,i+4,0)-vxy(xy,i,0))*(vxy(xy,j+4,1)-vxy(xy,j,1)) -
1644  (vxy(xy,j,0)-vxy(xy,j+4,0))*(vxy(xy,i,1)-vxy(xy,i+4,1))));
1645  }
1646  capacity = ::fabs(capacity);}
1647  if(capacity<=0) {
1648  return false;
1649  }
1650 
1651  AllocateMemory(8,6);
1652 
1653  pV[1] = HVPoint3D( vxy(xy,0,0), vxy(xy,0,1),-Dz);
1654  pV[2] = HVPoint3D( vxy(xy,1,0), vxy(xy,1,1),-Dz);
1655  pV[3] = HVPoint3D( vxy(xy,2,0), vxy(xy,2,1),-Dz);
1656  pV[4] = HVPoint3D( vxy(xy,3,0), vxy(xy,3,1),-Dz);
1657  pV[5] = HVPoint3D( vxy(xy,4,0), vxy(xy,4,1), Dz);
1658  pV[6] = HVPoint3D( vxy(xy,5,0), vxy(xy,5,1), Dz);
1659  pV[7] = HVPoint3D( vxy(xy,6,0), vxy(xy,6,1), Dz);
1660  pV[8] = HVPoint3D( vxy(xy,7,0), vxy(xy,7,1), Dz);
1661 
1662  CreatePrism();
1663 
1664  return true;
1665 }
1666 
1667 inline
1668 polyhedron_arb8::polyhedron_arb8(double Dz,const double* xy)
1669 {
1670  set_polyhedron_arb8(Dz,xy); //G.Barrand
1671 }
1672 
1673 //G.Barrand
1674 inline
1675 int polyhedron::_ixy(
1676  int a_ixy
1677 ,int a_npts
1678 ,int a_iz
1679 ,int a_nz
1680 ,bool a_acw
1681 ,bool a_zfb
1682 ){
1683  if(a_acw) {
1684  if(a_zfb) {
1685  return a_iz*a_npts+a_ixy;
1686  } else {
1687  return (a_nz-1-a_iz)*a_npts+a_ixy;
1688  }
1689  } else {
1690  if(a_zfb) {
1691  return a_iz*a_npts+(a_npts-1-a_ixy);
1692  } else {
1693  return (a_nz-1-a_iz)*a_npts+(a_npts-1-a_ixy);
1694  }
1695  }
1696 }
1697 
1698 inline
1699 bool polyhedron::set_polyhedron_xtru(
1700  int a_npts // number of vertices of the 2D polygon (at least 3)
1701 ,int a_nz // number of z planes (at least two)
1702 ,double* a_xs //[a_nz][a_npts] X positions for polygon vertices
1703 ,double* a_ys //[a_nz][a_npts] Y positions for polygon vertices
1704 ,double* a_zs //[a_nz] Z positions
1705 //default orientations :
1706 ,bool a_acw //= true
1707 ,bool a_zfb //= true
1708 ){
1709  _clear(); //G.Barrand
1710 
1711  if(a_npts<=2) return false;
1712  if(a_nz<=1) return false;
1713 
1714  //check if convex :
1715  bool convex = true;
1716  {double xv = a_xs[1]-a_xs[0];
1717  double yv = a_ys[1]-a_ys[0];
1718  double xw,yw,cross_z;
1719  for(int j=2;j<a_npts;j++) {
1720  xw = a_xs[j]-a_xs[j-1];
1721  yw = a_ys[j]-a_ys[j-1];
1722  //z of cross :
1723  // xv xw
1724  // yv yw
1725  // 0 0
1726  cross_z = xv*yw-yv*xw;
1727  if(a_acw) {
1728  if(cross_z<0) {convex = false;break;}
1729  } else {
1730  if(cross_z>0) {convex = false;break;}
1731  }
1732  xv = xw;
1733  yv = yw;
1734  }
1735  if(convex) {
1736  // have to check seg(n-2,n-1) to seg(n-1,0)
1737  xw = a_xs[0]-a_xs[a_npts-1];
1738  yw = a_ys[0]-a_ys[a_npts-1];
1739  cross_z = xv*yw-yv*xw;
1740  if(a_acw) {
1741  if(cross_z<0) convex = false;
1742  } else {
1743  if(cross_z>0) convex = false;
1744  }
1745  xv = xw;
1746  yv = yw;
1747  }
1748  if(convex) {
1749  // have to check seg(n-1,0) to seg(0,1)
1750  xw = a_xs[1]-a_xs[0];
1751  yw = a_ys[1]-a_ys[0];
1752  cross_z = xv*yw-yv*xw;
1753  if(a_acw) {
1754  if(cross_z<0) convex = false;
1755  } else {
1756  if(cross_z>0) convex = false;
1757  }
1758  }
1759  if(!convex) {
1760 #ifdef TOOLS_HEP_PH_OUT_ERR
1761  std::cerr << "tools::hep::set_polyhedron_xtru :"
1762  << " not convex polygon."
1763  << std::endl;
1764 #endif
1765  return false;
1766  }}
1767 
1768  // C O U N T V E R T E C E S
1769 
1770  int i;
1771 
1772  int j = a_nz*a_npts+2;
1773 
1774  // C O U N T F A C E S
1775 
1776  int k = ((a_nz-1)+1+1)*a_npts;
1777 
1778  // A L L O C A T E M E M O R Y
1779 
1780  AllocateMemory(j, k);
1781 
1782  // G E N E R A T E V E R T E C E S
1783 
1784  int* kk = new int[a_nz+2];
1785 
1786  k = 1;
1787  for(i=0; i<a_nz; i++) {
1788  kk[i] = k;
1789  k += a_npts;
1790  }
1791 
1792  kk[a_nz] = k;
1793  kk[a_nz+1] = k+1;
1794 
1795  int ixy,iz;
1796  if(a_acw) {
1797  if(a_zfb) {
1798  for(j=0; j<a_npts; j++) {
1799  for(i=0; i<a_nz; i++) {
1800  iz = i;
1801  ixy = iz*a_npts+j;
1802  pV[kk[i]+j] = HVPoint3D(a_xs[ixy],a_ys[ixy],a_zs[iz]);
1803  }
1804  }
1805  } else {
1806  for(j=0; j<a_npts; j++) {
1807  for(i=0; i<a_nz; i++) {
1808  iz = a_nz-1-i;
1809  ixy = iz*a_npts+j;
1810  pV[kk[i]+j] = HVPoint3D(a_xs[ixy],a_ys[ixy],a_zs[iz]);
1811  }
1812  }
1813  }
1814  } else {
1815  if(a_zfb) {
1816  for(j=0; j<a_npts; j++) {
1817  for(i=0; i<a_nz; i++) {
1818  iz = i;
1819  ixy = iz*a_npts+a_npts-1-j;
1820  pV[kk[i]+j] = HVPoint3D(a_xs[ixy],a_ys[ixy],a_zs[iz]);
1821  }
1822  }
1823  } else {
1824  for(j=0; j<a_npts; j++) {
1825  for(i=0; i<a_nz; i++) {
1826  iz = a_nz-1-i;
1827  ixy = iz*a_npts+a_npts-1-j;
1828  pV[kk[i]+j] = HVPoint3D(a_xs[ixy],a_ys[ixy],a_zs[iz]);
1829  }
1830  }
1831  }
1832  }
1833 
1834  {double xcbeg = 0;
1835  double ycbeg = 0;
1836  for(j=0; j<a_npts; j++) {
1837  ixy = _ixy(j,a_npts,0,a_nz,a_acw,a_zfb);
1838  xcbeg += a_xs[ixy];
1839  ycbeg += a_ys[ixy];
1840  }
1841  xcbeg /= a_npts;
1842  ycbeg /= a_npts;
1843  if(a_zfb) {
1844  pV[kk[a_nz]] = HVPoint3D(xcbeg,ycbeg,a_zs[0]);
1845  } else {
1846  pV[kk[a_nz]] = HVPoint3D(xcbeg,ycbeg,a_zs[a_nz-1]);
1847  }}
1848 
1849  {double xcend = 0;
1850  double ycend = 0;
1851  for(j=0; j<a_npts; j++) {
1852  ixy = _ixy(j,a_npts,a_nz-1,a_nz,a_acw,a_zfb);
1853  xcend += a_xs[ixy];
1854  ycend += a_ys[ixy];
1855  }
1856  xcend /= a_npts;
1857  ycend /= a_npts;
1858  if(a_zfb) {
1859  pV[kk[a_nz+1]] = HVPoint3D(xcend,ycend,a_zs[a_nz-1]);
1860  } else {
1861  pV[kk[a_nz+1]] = HVPoint3D(xcend,ycend,a_zs[0]);
1862  }}
1863 
1864  // G E N E R A T E E X T E R N A L F A C E S
1865  int edgeVis = 1;
1866 
1867  int v1 = 1;
1868  int v2 = 1;
1869  k = 1;
1870  for(i=0; i<(a_nz-1); i++) {
1871  RotateEdge(kk[i], kk[i+1], 1, 1, v1, v2,
1872  edgeVis, true, a_npts, k);
1873  }
1874 
1875  // G E N E R A T E S I D E F A C E S
1876  //begin side
1877  RotateEdge(kk[a_nz], kk[0], 0, 1, 1, 1,
1878  -1, true, a_npts, k);
1879  //end side
1880  RotateEdge(kk[a_nz-1], kk[a_nz+1], 1, 0, 1, 1,
1881  -1, true, a_npts, k);
1882 
1883  delete [] kk;
1884 
1885  if ((k-1) != nface) {
1886 #ifdef TOOLS_HEP_PH_OUT_ERR
1887  std::cerr
1888  << "Polyhedron::RotateAroundZ: number of generated faces ("
1889  << k-1 << ") is not equal to the number of allocated faces ("
1890  << nface << ")"
1891  << std::endl;
1892 #endif
1893  }
1894 
1895  SetReferences();
1896  return true;
1897 }
1898 
1899 inline
1900 polyhedron_xtru::polyhedron_xtru(int a_npts,int a_nz,
1901  double* a_xs,double* a_ys,double* a_zs,
1902  bool a_acw, //= true
1903  bool a_zfb){ //= true
1904  set_polyhedron_xtru(a_npts,a_nz,a_xs,a_ys,a_zs,a_acw,a_zfb);
1905 }
1906 
1907 //G.Barrand :
1908 inline
1909 bool polyhedron::set_polyhedron_hype(double a_st_in,double a_st_out,
1910  double a_rmin,double a_rmax,double a_dz,
1911  int a_nz,int a_nphi) //G.Barrand
1912 {
1913  static const double wholeCircle = 2*_M_PI(); //G.Barrand : const
1914 
1915  _clear(); //G.Barrand
1916 
1917  // C H E C K I N P U T P A R A M E T E R S
1918 
1919  int k = 0;
1920  if (a_rmin < 0. || a_rmax < 0.) k = 1;
1921  if (a_rmin > a_rmax) k = 1;
1922  if (a_rmin == a_rmax) k = 1;
1923 
1924  if (a_dz <= 0.) k += 2;
1925 
1926  if (a_nz <= 0) k += 4;
1927  if (a_nphi <= 0) k += 4;
1928 
1929  double tout = ::tan(a_st_in);
1930  double tin = ::tan(a_st_out);
1931 
1932  double d_in = a_rmin*a_rmin + tin*tin*a_dz*a_dz;
1933  double d_out = a_rmax*a_rmax + tout*tout*a_dz*a_dz;
1934  if(d_in>d_out) k += 8;
1935 
1936  double phi1 = 0;
1937  double dphi = wholeCircle;
1938 
1939  if (k != 0) {
1940 #ifdef TOOLS_HEP_PH_OUT_ERR
1941  std::cerr << "polyhedron_cone(s)/Tube(s): error in input parameters";
1942  if ((k & 1) != 0) std::cerr << " (radiuses)";
1943  if ((k & 2) != 0) std::cerr << " (half-length)";
1944  if ((k & 4) != 0) std::cerr << " (steps)";
1945  if ((k & 8) != 0) std::cerr << " (angles)";
1946  std::cerr << std::endl;
1947  //std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
1948  //std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
1949  //std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
1950  // << std::endl;
1951 #endif
1952  return false;
1953  }
1954 
1955  // P R E P A R E T W O P O L Y L I N E S
1956 
1957  double* zz = new double[2*(a_nz+1)];
1958  double* rr = new double[2*(a_nz+1)];
1959 
1960  double dz = 2.0f*a_dz/double(a_nz);
1961 
1962  double z;
1963 
1964  // r^2 - (tout*z)^2 = rout^2
1965  double rout_sq = a_rmax*a_rmax;
1966  double tout_sq = tout*tout;
1967  for(int iz=0;iz<=a_nz;iz++) {
1968  z = a_dz-iz*dz;
1969  zz[iz] = z;
1970  rr[iz] = ::sqrt(rout_sq+tout_sq*z*z);
1971  }
1972 
1973  // r^2 - (tin*z)^2 = rin^2
1974  double rin_sq = a_rmin*a_rmin;
1975  double tin_sq = tin*tin;
1976  for(int iz=0;iz<a_nz;iz++) {
1977  z = a_dz-iz*dz;
1978  zz[a_nz+iz] = z;
1979  rr[a_nz+iz] = ::sqrt(rin_sq+tin_sq*z*z);
1980  }
1981 
1982  // R O T A T E P O L Y L I N E S
1983 
1984  RotateAroundZ(a_nphi, phi1, dphi, a_nz, a_nz, zz, rr, -1, -1);
1985  SetReferences();
1986 
1987  delete [] zz;
1988  delete [] rr;
1989  return true;
1990 }
1991 
1992 inline
1993 polyhedron_hype::polyhedron_hype(double a_st_in,double a_st_out,
1994  double a_rmin,double a_rmax,double a_dz,
1995  int a_nz,int a_nphi) {
1996  set_polyhedron_hype(a_st_in,a_st_out,a_rmin,a_rmax,a_dz,a_nz,a_nphi);
1997 }
1998 
1999 //G.Barrand :
2000 inline
2001 bool polyhedron::set_polyhedron_eltu(double a_dx,double a_dy,double a_dz,
2002  int a_nz,int a_nphi) //G.Barrand
2003 {
2004  // elliptical tube class. An elliptical tube has 3 parameters :
2005  // a_dx - semi-axis of the ellipse along x
2006  // a_dy - semi-axis of the ellipse along y
2007  // a_dz - half length in z
2008  static const double wholeCircle = 2*_M_PI(); //G.Barrand : const
2009 
2010  _clear(); //G.Barrand
2011 
2012  // C H E C K I N P U T P A R A M E T E R S
2013 
2014  int k = 0;
2015 
2016  if (a_dx <= 0.) k += 1;
2017  if (a_dy <= 0.) k += 1;
2018  if (a_dz <= 0.) k += 1;
2019 
2020  if (a_nz <= 0) k += 2;
2021  if (a_nphi <= 0) k += 2;
2022 
2023  if (k != 0) {
2024 #ifdef TOOLS_HEP_PH_OUT_ERR
2025  std::cerr << "polyhedron_cone(s)/Tube(s): error in input parameters";
2026  if ((k & 1) != 0) std::cerr << " (half-length)";
2027  if ((k & 2) != 0) std::cerr << " (steps)";
2028  std::cerr << std::endl;
2029  std::cerr << " Dx=" << a_dx
2030  << " Dy=" << a_dy
2031  << " Dz=" << a_dz
2032  << " nz=" << a_nz
2033  << " nphi=" << a_nphi
2034  << std::endl;
2035 #endif
2036  return false;
2037  }
2038 
2039  // C O U N T V E R T E C E S
2040 
2041  int i;
2042 
2043  int j = a_nz*a_nphi+2;
2044 
2045  // C O U N T F A C E S
2046 
2047  k = ((a_nz-1)+1+1)*a_nphi;
2048 
2049  // A L L O C A T E M E M O R Y
2050 
2051  AllocateMemory(j, k);
2052 
2053  // G E N E R A T E V E R T E C E S
2054 
2055  int* kk = new int[a_nz+2];
2056 
2057  k = 1;
2058  for(i=0; i<a_nz; i++) {
2059  kk[i] = k;
2060  k += a_nphi;
2061  }
2062 
2063  kk[a_nz] = k;
2064  kk[a_nz+1] = k+1;
2065 
2066 
2067  {double a2 = a_dx*a_dx;
2068  double b2 = a_dy*a_dy;
2069  double dphi = wholeCircle/a_nphi;
2070  double phi = 0;
2071  double sph,cph,r2,r,x,y;
2072  double sz = (2.0*a_dz)/a_nz;
2073  for(j=0; j<a_nphi; j++) {
2074  phi = j*dphi;
2075  sph=::sin(phi);
2076  cph=::cos(phi);
2077  r2=(a2*b2)/(b2+(a2-b2)*sph*sph);
2078  r=::sqrt(r2);
2079  x = r*cph;
2080  y = r*sph;
2081  for(i=0; i<a_nz; i++) {
2082  pV[kk[i]+j] = HVPoint3D(x,y,a_dz-sz*i);
2083  }
2084  }}
2085 
2086  pV[kk[a_nz]] = HVPoint3D(0,0,a_dz);
2087  pV[kk[a_nz+1]] = HVPoint3D(0,0,-a_dz);
2088 
2089  // G E N E R A T E E X T E R N A L F A C E S
2090  int edgeVis = 1;
2091 
2092  int v1 = 1;
2093  int v2 = 1;
2094  k = 1;
2095  for(i=0; i<(a_nz-1); i++) {
2096  RotateEdge(kk[i], kk[i+1], 1, 1, v1, v2,
2097  edgeVis, true, a_nphi, k);
2098  }
2099 
2100  // G E N E R A T E S I D E F A C E S
2101  //begin side
2102  RotateEdge(kk[a_nz], kk[0], 0, 1, 1, 1,
2103  -1, true, a_nphi, k);
2104  //end side
2105  RotateEdge(kk[a_nz-1], kk[a_nz+1], 1, 0, 1, 1,
2106  -1, true, a_nphi, k);
2107 
2108  delete [] kk;
2109 
2110  if ((k-1) != nface) {
2111 #ifdef TOOLS_HEP_PH_OUT_ERR
2112  std::cerr
2113  << "Polyhedron::RotateAroundZ: number of generated faces ("
2114  << k-1 << ") is not equal to the number of allocated faces ("
2115  << nface << ")"
2116  << std::endl;
2117 #endif
2118  }
2119 
2120  SetReferences();
2121  return true;
2122 }
2123 
2124 inline
2125 polyhedron_trd1::polyhedron_trd1(double Dx1, double Dx2,
2126  double Dy, double Dz)
2127  : polyhedron_trd2(Dx1, Dx2, Dy, Dy, Dz) {}
2128 
2129 inline
2130 polyhedron_box::polyhedron_box(double Dx, double Dy, double Dz)
2131  : polyhedron_trd2(Dx, Dx, Dy, Dy, Dz) {}
2132 
2133 //G.Barrand
2134 inline
2135 bool polyhedron::set_polyhedron_trap(double Dz,
2136  double Theta,
2137  double Phi,
2138  double Dy1,
2139  double Dx1,
2140  double Dx2,
2141  double Alp1,
2142  double Dy2,
2143  double Dx3,
2144  double Dx4,
2145  double Alp2)
2146 /***********************************************************************
2147  * *
2148  * Name: polyhedron_trap Date: 20.11.96 *
2149  * Author: E.Chernyaev Revised: *
2150  * *
2151  * Function: Create GEANT4 TRAP-trapezoid *
2152  * *
2153  * Input: DZ - half-length in Z *
2154  * Theta,Phi - polar angles of the line joining centres of the *
2155  * faces at Z=-Dz and Z=+Dz *
2156  * Dy1 - half-length in Y of the face at Z=-Dz *
2157  * Dx1 - half-length in X of low edge of the face at Z=-Dz *
2158  * Dx2 - half-length in X of top edge of the face at Z=-Dz *
2159  * Alp1 - angle between Y-axis and the median joining top and *
2160  * low edges of the face at Z=-Dz *
2161  * Dy2 - half-length in Y of the face at Z=+Dz *
2162  * Dx3 - half-length in X of low edge of the face at Z=+Dz *
2163  * Dx4 - half-length in X of top edge of the face at Z=+Dz *
2164  * Alp2 - angle between Y-axis and the median joining top and *
2165  * low edges of the face at Z=+Dz *
2166  * *
2167  ***********************************************************************/
2168 {
2169  _clear(); //G.Barrand
2170 
2171  //FIXME : check capacity = 0; //see TGeoTrap to set a arb8.
2172 
2173  if(//(Dx1<=0.0)||(Dx2<=0.0)||(Dx3<=0.0)||(Dx4<=0.0)||
2174  //(Dy1<=0.0)||(Dy2<=0.0)||
2175  (Dz<=0.0)){ //G.Barrand
2176 #ifdef TOOLS_HEP_PH_OUT_ERR
2177  std::cerr << "set_polyhedron_trap: error in input parameters";
2178  std::cerr << " Dx1=" << Dx1
2179  << " Dx2=" << Dx2
2180  << " Dy1=" << Dy1
2181  << " Dy2=" << Dy2
2182  << " Dz=" << Dz
2183  << std::endl;
2184 #endif
2185  return false;
2186  }
2187 
2188  double DzTthetaCphi = Dz*tan(Theta)*cos(Phi);
2189  double DzTthetaSphi = Dz*tan(Theta)*sin(Phi);
2190  double Dy1Talp1 = Dy1*tan(Alp1);
2191  double Dy2Talp2 = Dy2*tan(Alp2);
2192 
2193  AllocateMemory(8,6);
2194 
2195  pV[1] = HVPoint3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
2196  pV[2] = HVPoint3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
2197  pV[3] = HVPoint3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
2198  pV[4] = HVPoint3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
2199  pV[5] = HVPoint3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
2200  pV[6] = HVPoint3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
2201  pV[7] = HVPoint3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
2202  pV[8] = HVPoint3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
2203 
2204  CreatePrism();
2205 
2206  return true;
2207 }
2208 
2209 inline
2210 polyhedron_trap::polyhedron_trap(double Dz,
2211  double Theta,
2212  double Phi,
2213  double Dy1,
2214  double Dx1,
2215  double Dx2,
2216  double Alp1,
2217  double Dy2,
2218  double Dx3,
2219  double Dx4,
2220  double Alp2)
2221 {
2222  //G.Barrand
2223  set_polyhedron_trap(Dz,Theta,Phi,Dy1,Dx1,Dx2,Alp1,Dy2,Dx3,Dx4,Alp2);
2224 }
2225 
2226 inline
2227 polyhedron_para::polyhedron_para(double Dx, double Dy, double Dz,
2228  double Alpha, double Theta,
2229  double Phi)
2230  : polyhedron_trap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
2231 
2232 //G.Barrand : have the below set_ to optimize exlib/sg/polyhedron setup.
2233 inline
2234 bool polyhedron::set_polyhedron_cons(double Rmn1,
2235  double Rmx1,
2236  double Rmn2,
2237  double Rmx2,
2238  double Dz,
2239  double Phi1,
2240  double Dphi,
2241  int nstep) //G.Barrand
2242 /***********************************************************************
2243  * *
2244  * Name: polyhedron_cons::polyhedron_cons Date: 15.12.96 *
2245  * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
2246  * *
2247  * Function: Constructor for CONS, TUBS, CONE, TUBE *
2248  * *
2249  * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
2250  * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
2251  * Dz - half length in Z *
2252  * Phi1 - starting angle of the segment *
2253  * Dphi - segment range *
2254  * *
2255  ***********************************************************************/
2256 {
2257  static const double wholeCircle = 2*_M_PI(); //G.Barrand : const
2258  _clear(); //G.Barrand
2259 
2260  // C H E C K I N P U T P A R A M E T E R S
2261 
2262  int k = 0;
2263  if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
2264  if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
2265  if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
2266 
2267  //G.Barrand : for atlas.root.
2268  //if (Dz <= 0.) k += 2;
2269  const double perMillion = 0.000001;
2270  if (Dz <= 0.) Dz = perMillion*mx<double>(Rmx1,Rmx2);
2271 
2272  double phi1, phi2, dphi;
2273  if (Dphi < 0.) {
2274  phi2 = Phi1; phi1 = phi2 - Dphi;
2275  }else if (Dphi == 0.) {
2276  phi1 = Phi1; phi2 = phi1 + wholeCircle;
2277  }else{
2278  phi1 = Phi1; phi2 = phi1 + Dphi;
2279  }
2280  dphi = phi2 - phi1;
2281 
2282  if (fabs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
2283  if (dphi > wholeCircle) k += 4;
2284 
2285  if (k != 0) {
2286 #ifdef TOOLS_HEP_PH_OUT_ERR
2287  std::cerr << "polyhedron_cone(s)/Tube(s): error in input parameters";
2288  if ((k & 1) != 0) std::cerr << " (radiuses)";
2289  if ((k & 2) != 0) std::cerr << " (half-length)";
2290  if ((k & 4) != 0) std::cerr << " (angles)";
2291  std::cerr << std::endl;
2292  std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
2293  std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
2294  std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
2295  << std::endl;
2296 #endif
2297  return false;
2298  }
2299 
2300  // P R E P A R E T W O P O L Y L I N E S
2301 
2302  double zz[4], rr[4];
2303  zz[0] = Dz;
2304  zz[1] = -Dz;
2305  zz[2] = Dz;
2306  zz[3] = -Dz;
2307  rr[0] = Rmx2;
2308  rr[1] = Rmx1;
2309  rr[2] = Rmn2;
2310  rr[3] = Rmn1;
2311 
2312  // R O T A T E P O L Y L I N E S
2313 
2314  //G.Barrand : nstep
2315  RotateAroundZ(nstep, phi1, dphi, 2, 2, zz, rr, -1, -1);
2316  SetReferences();
2317 
2318  return true;
2319 }
2320 
2321 inline
2322 polyhedron_cons::polyhedron_cons(double Rmn1,
2323  double Rmx1,
2324  double Rmn2,
2325  double Rmx2,
2326  double Dz,
2327  double Phi1,
2328  double Dphi,
2329  int nstep) //G.Barrand
2330 {
2331  set_polyhedron_cons(Rmn1,Rmx1,Rmn2,Rmx2,Dz,Phi1,Dphi,nstep); //G.Barrand
2332 }
2333 
2334 inline
2335 polyhedron_cone::polyhedron_cone(double Rmn1, double Rmx1,
2336  double Rmn2, double Rmx2,
2337  double Dz,
2338  int nstep) //G.Barrand
2339 : polyhedron_cons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0, 2*_M_PI(), nstep) {}
2340 
2341 inline
2342 polyhedron_tubs::polyhedron_tubs(double Rmin, double Rmax,
2343  double Dz,
2344  double Phi1, double Dphi,
2345  int nstep) //G.Barrand
2346 : polyhedron_cons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi, nstep) {}
2347 
2348 inline
2349 polyhedron_tube::polyhedron_tube (double Rmin, double Rmax,
2350  double Dz,
2351  int nstep) //G.Barrand
2352 : polyhedron_cons(Rmin, Rmax, Rmin, Rmax, Dz, 0, 2*_M_PI(), nstep) {}
2353 
2354 //G.Barrand
2355 inline
2356 bool polyhedron::set_polyhedron_pgon(double phi,
2357  double dphi,
2358  int npdv,
2359  int nz,
2360  const double *z,
2361  const double *rmin,
2362  const double *rmax)
2363 /***********************************************************************
2364  * *
2365  * Name: polyhedron_pgon Date: 09.12.96 *
2366  * Author: E.Chernyaev Revised: *
2367  * *
2368  * Function: Constructor of polyhedron for PGON, PCON *
2369  * *
2370  * Input: phi - initial phi *
2371  * dphi - delta phi *
2372  * npdv - number of steps along phi *
2373  * nz - number of z-planes (at least two) *
2374  * z[] - z coordinates of the slices *
2375  * rmin[] - smaller r at the slices *
2376  * rmax[] - bigger r at the slices *
2377  * *
2378  ***********************************************************************/
2379 {
2380  _clear();
2381 
2382  // C H E C K I N P U T P A R A M E T E R S
2383 
2384  if (dphi <= 0. || dphi > 2*_M_PI()) {
2385 #ifdef TOOLS_HEP_PH_OUT_ERR
2386  std::cerr
2387  << "polyhedron_pgon/Pcon: wrong delta phi = " << dphi
2388  << std::endl;
2389 #endif
2390  return false;
2391  }
2392 
2393  if (nz < 2) {
2394 #ifdef TOOLS_HEP_PH_OUT_ERR
2395  std::cerr
2396  << "polyhedron_pgon/Pcon: number of z-planes less than two = " << nz
2397  << std::endl;
2398 #endif
2399  return false;
2400  }
2401 
2402  if (npdv < 0) {
2403 #ifdef TOOLS_HEP_PH_OUT_ERR
2404  std::cerr
2405  << "polyhedron_pgon/Pcon: error in number of phi-steps =" << npdv
2406  << std::endl;
2407 #endif
2408  return false;
2409  }
2410 
2411  int i;
2412  for (i=0; i<nz; i++) {
2413  if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
2414 #ifdef TOOLS_HEP_PH_OUT_ERR
2415  std::cerr
2416  << "polyhedron_pgon: error in radiuses rmin[" << i << "]="
2417  << rmin[i] << " rmax[" << i << "]=" << rmax[i]
2418  << std::endl;
2419 #endif
2420  return false;
2421  }
2422  }
2423 
2424  // P R E P A R E T W O P O L Y L I N E S
2425 
2426  double *zz, *rr;
2427  zz = new double[2*nz];
2428  rr = new double[2*nz];
2429 
2430  if (z[0] > z[nz-1]) {
2431  for (i=0; i<nz; i++) {
2432  zz[i] = z[i];
2433  rr[i] = rmax[i];
2434  zz[i+nz] = z[i];
2435  rr[i+nz] = rmin[i];
2436  }
2437  }else{
2438  for (i=0; i<nz; i++) {
2439  zz[i] = z[nz-i-1];
2440  rr[i] = rmax[nz-i-1];
2441  zz[i+nz] = z[nz-i-1];
2442  rr[i+nz] = rmin[nz-i-1];
2443  }
2444  }
2445 
2446  // R O T A T E P O L Y L I N E S
2447 
2448  RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1);
2449  SetReferences();
2450 
2451  delete [] zz;
2452  delete [] rr;
2453 
2454  return true;
2455 }
2456 
2457 inline
2458 polyhedron_pgon::polyhedron_pgon(double phi,
2459  double dphi,
2460  int npdv,
2461  int nz,
2462  const double *z,
2463  const double *rmin,
2464  const double *rmax)
2465 {
2466  set_polyhedron_pgon(phi,dphi,npdv,nz,z,rmin,rmax);
2467 }
2468 
2469 inline
2470 polyhedron_pcon::polyhedron_pcon(double phi, double dphi, int nz,
2471  const double *z,
2472  const double *rmin,
2473  const double *rmax)
2474  : polyhedron_pgon(phi, dphi, 0, nz, z, rmin, rmax) {}
2475 
2476 inline
2477 bool polyhedron::set_polyhedron_sphere(double rmin, double rmax,
2478  double phi, double dphi,
2479  double the, double dthe,
2480  int nphi, //G.Barrand
2481  int nthe) //G.Barrand
2482 /***********************************************************************
2483  * *
2484  * Name: polyhedron_sphere Date: 11.12.96 *
2485  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
2486  * *
2487  * Function: Constructor of polyhedron for SPHERE *
2488  * *
2489  * Input: rmin - internal radius *
2490  * rmax - external radius *
2491  * phi - initial phi *
2492  * dphi - delta phi *
2493  * the - initial theta *
2494  * dthe - delta theta *
2495  * *
2496  ***********************************************************************/
2497 {
2498  _clear();
2499 
2500  // C H E C K I N P U T P A R A M E T E R S
2501 
2502  if (dphi <= 0. || dphi > 2*_M_PI()) {
2503 #ifdef TOOLS_HEP_PH_OUT_ERR
2504  std::cerr
2505  << "polyhedron_sphere: wrong delta phi = " << dphi
2506  << std::endl;
2507 #endif
2508  return false;
2509  }
2510 
2511  if (the < 0. || the > _M_PI()) {
2512 #ifdef TOOLS_HEP_PH_OUT_ERR
2513  std::cerr
2514  << "polyhedron_sphere: wrong theta = " << the
2515  << std::endl;
2516 #endif
2517  return false;
2518  }
2519 
2520  if (dthe <= 0. || dthe > _M_PI()) {
2521 #ifdef TOOLS_HEP_PH_OUT_ERR
2522  std::cerr
2523  << "polyhedron_sphere: wrong delta theta = " << dthe
2524  << std::endl;
2525 #endif
2526  return false;
2527  }
2528 
2529  if ( (the+dthe >= _M_PI()) && (the+dthe < _M_PI() + 2*DBL_EPSILON) )
2530  dthe = _M_PI() - the; //G.Barrand : coming from LHCb/S.Ponce.
2531 
2532  if (the+dthe > _M_PI()) {
2533 #ifdef TOOLS_HEP_PH_OUT_ERR
2534  std::cerr
2535  << "polyhedron_sphere: wrong theta + delta theta = "
2536  << the << " " << dthe
2537  << std::endl;
2538 #endif
2539  return false;
2540  }
2541 
2542  if (rmin < 0. || rmin >= rmax) {
2543 #ifdef TOOLS_HEP_PH_OUT_ERR
2544  std::cerr
2545  << "polyhedron_sphere: error in radiuses"
2546  << " rmin=" << rmin << " rmax=" << rmax
2547  << std::endl;
2548 #endif
2549  return false;
2550  }
2551 
2552  // P R E P A R E T W O P O L Y L I N E S
2553 
2554  int n_the = (nthe>0) ? nthe : GetNumberOfRotationSteps(); //G.Barrand.
2555  int ns = (n_the + 1) / 2;
2556 
2557  int np1 = int(dthe*ns/_M_PI()+.5) + 1;
2558  if (np1 <= 1) np1 = 2;
2559  const double perMillion = 0.000001;
2560  int np2 = rmin < perMillion ? 1 : np1;
2561 
2562  double *zz, *rr;
2563  zz = new double[np1+np2];
2564  rr = new double[np1+np2];
2565 
2566  double a = dthe/(np1-1);
2567  double cosa, sina;
2568  for (int i=0; i<np1; i++) {
2569  cosa = cos(the+i*a);
2570  sina = sin(the+i*a);
2571  zz[i] = rmax*cosa;
2572  rr[i] = rmax*sina;
2573  if (np2 > 1) {
2574  zz[i+np1] = rmin*cosa;
2575  rr[i+np1] = rmin*sina;
2576  }
2577  }
2578  if (np2 == 1) {
2579  zz[np1] = 0.;
2580  rr[np1] = 0.;
2581  }
2582 
2583  // R O T A T E P O L Y L I N E S
2584 
2585  //G.Barrand : nphi.
2586  RotateAroundZ(nphi, phi, dphi, np1, np2, zz, rr, -1, -1);
2587  SetReferences();
2588 
2589  delete [] zz;
2590  delete [] rr;
2591 
2592  return true;
2593 }
2594 
2595 inline
2596 polyhedron_sphere::polyhedron_sphere(double rmin, double rmax,
2597  double phi, double dphi,
2598  double the, double dthe,
2599  int nphi, //G.Barrand
2600  int nthe) //G.Barrand
2601 {
2602  set_polyhedron_sphere(rmin,rmax,
2603  phi,dphi,
2604  the,dthe,
2605  nphi,nthe);
2606 }
2607 
2608 
2609 inline
2610 bool polyhedron::set_polyhedron_torus(double rmin,
2611  double rmax,
2612  double rtor,
2613  double phi,
2614  double dphi,
2615  int nphi, //G.Barrand
2616  int nthe) //G.Barrand
2617 /***********************************************************************
2618  * *
2619  * Name: polyhedron_torus Date: 11.12.96 *
2620  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
2621  * *
2622  * Function: Constructor of polyhedron for TORUS *
2623  * *
2624  * Input: rmin - internal radius *
2625  * rmax - external radius *
2626  * rtor - radius of torus *
2627  * phi - initial phi *
2628  * dphi - delta phi *
2629  * *
2630  ***********************************************************************/
2631 {
2632  _clear();
2633 
2634  // C H E C K I N P U T P A R A M E T E R S
2635 
2636  if (dphi <= 0. || dphi > 2*_M_PI()) {
2637 #ifdef TOOLS_HEP_PH_OUT_ERR
2638  std::cerr
2639  << "polyhedron_torus: wrong delta phi = " << dphi
2640  << std::endl;
2641 #endif
2642  return false;
2643  }
2644 
2645  if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2646 #ifdef TOOLS_HEP_PH_OUT_ERR
2647  std::cerr
2648  << "polyhedron_torus: error in radiuses"
2649  << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
2650  << std::endl;
2651 #endif
2652  return false;
2653  }
2654 
2655  // P R E P A R E T W O P O L Y L I N E S
2656 
2657  int n_the = (nthe>0) ? nthe : GetNumberOfRotationSteps(); //G.Barrand.
2658  int np1 = n_the;
2659 
2660  const double perMillion = 0.000001;
2661  int np2 = rmin < perMillion ? 1 : np1;
2662 
2663  double *zz, *rr;
2664  zz = new double[np1+np2];
2665  rr = new double[np1+np2];
2666 
2667  double a = 2*_M_PI()/np1;
2668  double cosa, sina;
2669  for (int i=0; i<np1; i++) {
2670  cosa = cos(i*a);
2671  sina = sin(i*a);
2672  zz[i] = rmax*cosa;
2673  rr[i] = rtor+rmax*sina;
2674  if (np2 > 1) {
2675  zz[i+np1] = rmin*cosa;
2676  rr[i+np1] = rtor+rmin*sina;
2677  }
2678  }
2679  if (np2 == 1) {
2680  zz[np1] = 0.;
2681  rr[np1] = rtor;
2682  np2 = -1;
2683  }
2684 
2685  // R O T A T E P O L Y L I N E S
2686 
2687  //G.Barrand : nphi.
2688  RotateAroundZ(nphi, phi, dphi, -np1, -np2, zz, rr, -1,-1);
2689  SetReferences();
2690 
2691  delete [] zz;
2692  delete [] rr;
2693 
2694  return true;
2695 }
2696 
2697 inline
2698 polyhedron_torus::polyhedron_torus(double rmin,
2699  double rmax,
2700  double rtor,
2701  double phi,
2702  double dphi,
2703  int nphi, //G.Barrand
2704  int nthe) //G.Barrand
2705 {
2706  set_polyhedron_torus(rmin,rmax,rtor,phi,dphi,nphi,nthe); //G.Barrand
2707 }
2708 
2709 //int polyhedron::fNumberOfRotationSteps = NUMBER_OF_STEPS;
2710 
2711 // G.Barrand : begin.
2712 inline int polyhedron::GetNumberOfRotationSteps() {
2713  return fNumberOfRotationSteps;
2714 }
2715 inline void polyhedron::ResetNumberOfRotationSteps() {
2716  fNumberOfRotationSteps = NUMBER_OF_STEPS();
2717 }
2718 // G.Barrand : end.
2719 /***********************************************************************
2720  * *
2721  * Name: polyhedron::fNumberOfRotationSteps Date: 24.06.97 *
2722  * Author: J.Allison (Manchester University) Revised: *
2723  * *
2724  * Function: Number of steps for whole circle *
2725  * *
2726  ***********************************************************************/
2727 
2728 }}
2729 
2730 #include "pbp.icc" //BooleanProcessor
2731 
2732 namespace tools {
2733 namespace hep {
2734 
2735 //G.Barrand : static BooleanProcessor processor;
2736 
2737 inline polyhedron polyhedron::add(const polyhedron & p) const
2738 /***********************************************************************
2739  * *
2740  * Name: polyhedron::add Date: 19.03.00 *
2741  * Author: E.Chernyaev Revised: *
2742  * *
2743  * Function: Boolean "union" of two polyhedra *
2744  * *
2745  ***********************************************************************/
2746 {
2747  BooleanProcessor processor; //G.Barrand
2748  int err;
2749  return processor.execute(OP_UNION, *this, p, err);
2750 }
2751 
2752 inline polyhedron polyhedron::intersect(const polyhedron & p) const
2753 /***********************************************************************
2754  * *
2755  * Name: polyhedron::intersect Date: 19.03.00 *
2756  * Author: E.Chernyaev Revised: *
2757  * *
2758  * Function: Boolean "intersection" of two polyhedra *
2759  * *
2760  ***********************************************************************/
2761 {
2762  BooleanProcessor processor; //G.Barrand
2763  int err;
2764  return processor.execute(OP_INTERSECTION, *this, p, err);
2765 }
2766 
2767 inline polyhedron polyhedron::subtract(const polyhedron & p) const
2768 /***********************************************************************
2769  * *
2770  * Name: polyhedron::add Date: 19.03.00 *
2771  * Author: E.Chernyaev Revised: *
2772  * *
2773  * Function: Boolean "subtraction" of "p" from "this" *
2774  * *
2775  ***********************************************************************/
2776 {
2777  BooleanProcessor processor; //G.Barrand
2778  int err;
2779  return processor.execute(OP_SUBTRACTION, *this, p, err);
2780 }
2781 
2782 
2783 //G.Barrand : begin
2784 
2785 inline bool is_in(unsigned int a_index,
2786  const std::list<unsigned int>& a_is) {
2787  std::list<unsigned int>::const_iterator it;
2788  for(it=a_is.begin();it!=a_is.end();++it) {
2789  if(*it==a_index) return true;
2790  }
2791  return false;
2792 }
2793 
2794 class bijection_visitor {
2795 #ifdef TOOLS_MEM
2796  TOOLS_SCLASS(tools::hep::bijection_visitor)
2797 #endif
2798 public:
2799  typedef std::vector<unsigned int> is_t;
2800  virtual bool visit(const is_t&) = 0;
2801 public:
2802  bijection_visitor(unsigned int a_number):m_number(a_number){
2803 #ifdef TOOLS_MEM
2804  mem::increment(s_class().c_str());
2805 #endif
2806  }
2807  virtual ~bijection_visitor(){
2808 #ifdef TOOLS_MEM
2809  mem::decrement(s_class().c_str());
2810 #endif
2811  }
2812 protected:
2813  bijection_visitor(const bijection_visitor&){
2814 #ifdef TOOLS_MEM
2815  mem::increment(s_class().c_str());
2816 #endif
2817  }
2818  bijection_visitor& operator=(const bijection_visitor&){return *this;}
2819 public:
2820  bool visitx() {
2821  m_is.clear();
2822  m_is.resize(m_number,0);
2823  std::list<unsigned int> is;
2824  return visit(0,is);
2825  }
2826 private:
2827  bool visit(unsigned int a_level,std::list<unsigned int>& a_is) {
2828  for(unsigned int index=0;index<m_number;index++) {
2829  if(is_in(index,a_is)) {
2830  } else {
2831  a_is.push_back(index);
2832  m_is[a_level] = index;
2833  if(a_level==m_number-1) {
2834  if(!visit(m_is)) return false;
2835  } else {
2836  if(!visit(a_level+1,a_is)) return false;
2837  }
2838  a_is.pop_back();
2839  }
2840  }
2841  return true;
2842  }
2843 private:
2844  unsigned int m_number;
2845  is_t m_is;
2846 };
2847 
2848 //inline void dump(const std::vector<unsigned int>& a_is) {
2849 // unsigned int number = a_is.size();
2850 // for(unsigned int index=0;index<number;index++) {
2851 // printf("%d ",a_is[index]);
2852 // }
2853 // printf("\n");
2854 //}
2855 
2856 //class bijection_dump : public bijection_visitor {
2857 //public:
2858 // bijection_dump(unsigned int a_number)
2859 // : bijection_visitor(a_number)
2860 // {}
2861 // virtual bool visit(const is_t& a_is) {
2862 // dump(a_is);
2863 // return true;//continue
2864 // }
2865 //};
2866 
2867 class polyhedron_exec : public bijection_visitor {
2868 #ifdef TOOLS_MEM
2869  TOOLS_SCLASS(tools::hep::polyhedron_exec)
2870 #endif
2871 public:
2872  polyhedron_exec(unsigned int a_number,
2873  polyhedronProcessor& a_proc,
2874  polyhedron& a_poly)
2875  : bijection_visitor(a_number)
2876  ,m_proc(a_proc)
2877  ,m_poly(a_poly)
2878  {
2879 #ifdef TOOLS_MEM
2880  mem::increment(s_class().c_str());
2881 #endif
2882  }
2883  virtual ~polyhedron_exec(){
2884 #ifdef TOOLS_MEM
2885  mem::decrement(s_class().c_str());
2886 #endif
2887  }
2888 protected:
2889  polyhedron_exec(const polyhedron_exec& a_from)
2890  :bijection_visitor(a_from)
2891  ,m_proc(a_from.m_proc)
2892  ,m_poly(a_from.m_poly)
2893  {
2894 #ifdef TOOLS_MEM
2895  mem::increment(s_class().c_str());
2896 #endif
2897  }
2898  polyhedron_exec& operator=(const polyhedron_exec&){return *this;}
2899 public:
2900  virtual bool visit(const is_t& a_is) {
2901  if(m_proc.execute1(m_poly,a_is)==true) return false; //stop
2902  return true;//continue
2903  }
2904 private:
2905  polyhedronProcessor& m_proc;
2906  polyhedron& m_poly;
2907 };
2908 
2909 inline bool polyhedronProcessor::execute(polyhedron& a_poly) {
2910  //{for(unsigned int index=0;index<5;index++) {
2911  // printf("debug : bijection : %d\n",index);
2912  // bijection_dump bd(index);
2913  // bd.visitx();
2914  //}}
2915 
2916  polyhedron_exec e((unsigned int)m_ops.size(),*this,a_poly);
2917  if(!e.visitx()) return true;
2918 #ifdef TOOLS_HEP_PH_OUT_ERR
2919  //std::cerr << "polyhedronProcessor::execute :"
2920  // << " all shifts and combinatory tried."
2921  // << " Boolean operations failed."
2922  // << std::endl;
2923 #endif
2924  return false;
2925 }
2926 inline bool polyhedronProcessor::execute1(
2927  polyhedron& a_poly
2928 ,const std::vector<unsigned int>& a_is
2929 ) {
2930  polyhedron result(a_poly);
2931  size_t number = m_ops.size();
2932  int num_shift = BooleanProcessor::get_num_shift();
2933  for(int ishift=0;ishift<num_shift;ishift++) {
2934  BooleanProcessor::set_shift(ishift);
2935 
2936  result = a_poly;
2937  bool done = true;
2938  for(size_t index=0;index<number;index++) {
2939  BooleanProcessor processor; //take a fresh one.
2940  const op_t& elem = m_ops[a_is[index]];
2941  int err;
2942  result = processor.execute(elem.first,result,elem.second,err);
2943  if(err) {
2944  done = false;
2945  break;
2946  }
2947  }
2948  if(done) {
2949  a_poly = result;
2950  return true;
2951  }
2952  }
2953 
2954 #ifdef TOOLS_HEP_PH_OUT_ERR
2955  //std::cerr << "polyhedronProcessor::execute :"
2956  // << " all shifts tried. Boolean operations failed."
2957  // << std::endl;
2958 #endif
2959 
2960  //a_poly = result;
2961  return false;
2962 }
2963 
2964 }}
2965 //G.Barrand : end
Generated by doxygen 1.8.20