ATLAS Offline Software
SbPolyhedron.cxx
Go to the documentation of this file.
1 // //
3 // Source file for class SbPolyhedron //
4 // //
5 // Update:
6 // Riccardo-Maria BIANCHI (rbianchi@cern.ch) //
7 // 13.12.2012 //
8 // //
10 
11 
12 // this :
13 #include <VP1HEPVis/SbPolyhedron.h>
14 #include <cassert>
15 
16 #define perMillion 0.000001
17 #define deg (M_PI/180.0)
18 
19 #include <float.h> //G.Barrand : to have DBL_EPSILON on Windows.
20 
21 // G.Barrand : introduce iabs to avoid a mess with cmath and some compiler.
22 inline int iabs(int a) {
23  return a < 0 ? -a : a;
24 }
25 inline float ffabs(float a) {
26  return a < 0.0f ? -a : a;
27 }
28 
29 
30 // rbianchi
31 HVPoint3D::HVPoint3D(): SbVec3d(0,0,0){}
32 HVPoint3D::HVPoint3D(double x,double y,double z): SbVec3d(x,y,z){}
33 HVPoint3D::HVPoint3D(const HVPoint3D& v): SbVec3d(v){}
34 HVPoint3D::HVPoint3D(const SbVec3d& v): SbVec3d(v){}
36  setValue(v[0],v[1],v[2]);
37  return *this;
38 }
39 HVPoint3D& HVPoint3D::operator=(const SbVec3d& v) {
40  setValue(v[0],v[1],v[2]);
41  return *this;
42 }
44  return HVPoint3D(v1[0] + v2[0],v1[1] + v2[1],v1[2] + v2[2]);
45 }
46 //---
47 
48 
49 
50 //--------------------------------------------------------------------//
51 // JFB: //
52 // SbPolyhedron was SbPolyhedron, retrofitted to Open Inventor //
53 // infrastructure: //
54 //--------------------------------------------------------------------//
55 
56 
57 //
58 // ********************************************************************
59 // * DISCLAIMER *
60 // * *
61 // * The following disclaimer summarizes all the specific disclaimers *
62 // * of contributors to this software. The specific disclaimers,which *
63 // * govern, are listed with their locations in: *
64 // * http://cern.ch/geant4/license *
65 // * *
66 // * Neither the authors of this software system, nor their employing *
67 // * institutes,nor the agencies providing financial support for this *
68 // * work make any representation or warranty, express or implied, *
69 // * regarding this software system or assume any liability for its *
70 // * use. *
71 // * *
72 // * This code implementation is the intellectual property of the *
73 // * GEANT4 collaboration. *
74 // * By copying, distributing or modifying the Program (or any work *
75 // * based on the Program) you indicate your acceptance of this *
76 // * statement, and all its terms. *
77 // ********************************************************************
78 //
79 //
80 // $Id: SbPolyhedron.cxx,v 1.3 2008-09-14 19:04:40 tkittel Exp $
81 // GEANT4 tag $Name: not supported by cvs2svn $
82 //
83 //
84 //
85 // G4 Polyhedron library
86 //
87 // History:
88 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
89 //
90 // 30.09.96 E.Chernyaev
91 // - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
92 // - added GetNextUnitNormal, GetNextEdgeIndeces, GetNextEdge
93 //
94 // 15.12.96 E.Chernyaev
95 // - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
96 // - rewritten G4PolyhedronCons;
97 // - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
98 //
99 // 01.06.97 E.Chernyaev
100 // - modified RotateAroundZ, added SetSideFacets
101 //
102 // 19.03.00 E.Chernyaev
103 // - implemented boolean operations (add, subtract, intersect) on polyhedra;
104 //
105 // 25.05.01 E.Chernyaev
106 // - added GetSurfaceArea() and GetVolume();
107 //
108 
109 
110 /***********************************************************************
111  * *
112  * Name: SbPolyhedron operator << Date: 09.05.96 *
113  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
114  * *
115  * Function: Print contents of G4 polyhedron *
116  * *
117  ***********************************************************************/
118 std::ostream & operator<<(std::ostream & ostr, const SbFacet & facet) {
119  for (int k=0; k<4; k++) {
120  ostr << " " << facet.m_edge[k].v << "/" << facet.m_edge[k].f;
121  }
122  return ostr;
123 }
124 
125 std::ostream & operator<<(std::ostream & ostr, const SbPolyhedron & ph) {
126  ostr << std::endl;
127  ostr << "Nverteces=" << ph.m_nvert << ", Nfacets=" << ph.m_nface << std::endl;
128  int i;
129  for (i=1; i<=ph.m_nvert; i++) {
130  ostr << "xyz(" << i << ")="
131  << ph.m_pV[i][0] << ' ' << ph.m_pV[i][1] << ' ' << ph.m_pV[i][2]
132  << std::endl;
133  }
134  for (i=1; i<=ph.m_nface; i++) {
135  ostr << "face(" << i << ")=" << ph.m_pF[i] << std::endl;
136  }
137  return ostr;
138 }
139 
141 /***********************************************************************
142  * *
143  * Name: SbPolyhedron copy constructor Date: 23.07.96 *
144  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
145  * *
146  ***********************************************************************/
147 {
148  if (from.m_nvert > 0 && from.m_nface > 0) {
149  m_nvert = from.m_nvert;
150  m_nface = from.m_nface;
151  m_pV = new HVPoint3D[m_nvert + 1];
152  m_pF = new SbFacet[m_nface + 1];
153  int i;
154  for (i=1; i<=m_nvert; i++) m_pV[i] = from.m_pV[i];
155  for (i=1; i<=m_nface; i++) m_pF[i] = from.m_pF[i];
156  }else{
157  m_nvert = 0; m_nface = 0; m_pV = 0; m_pF = 0;
158  }
159 }
160 
162 /***********************************************************************
163  * *
164  * Name: SbPolyhedron operator = Date: 23.07.96 *
165  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
166  * *
167  * Function: Copy contents of one GEANT4 polyhedron to another *
168  * *
169  ***********************************************************************/
170 {
171  if (this == &from) return *this;
172  delete [] m_pV;
173  delete [] m_pF;
174  if (from.m_nvert > 0 && from.m_nface > 0) {
175  m_nvert = from.m_nvert;
176  m_nface = from.m_nface;
177  m_pV = new HVPoint3D[m_nvert + 1];
178  m_pF = new SbFacet[m_nface + 1];
179  int i;
180  for (i=1; i<=m_nvert; i++) m_pV[i] = from.m_pV[i];
181  for (i=1; i<=m_nface; i++) m_pF[i] = from.m_pF[i];
182  }else{
183  m_nvert = 0; m_nface = 0; m_pV = 0; m_pF = 0;
184  }
185  return *this;
186 }
187 
188 int
189 SbPolyhedron::FindNeighbour(int iFace, int iNode, int iOrder) const
190 /***********************************************************************
191  * *
192  * Name: SbPolyhedron::FindNeighbour Date: 22.11.99 *
193  * Author: E.Chernyaev Revised: *
194  * *
195  * Function: Find neighbouring face *
196  * *
197  ***********************************************************************/
198 {
199  int i;
200  for (i=0; i<4; i++) {
201  if (iNode == iabs(m_pF[iFace].m_edge[i].v)) break;
202  }
203  if (i == 4) {
204  std::cerr
205  << "SbPolyhedron::FindNeighbour: face " << iFace
206  << " has no node " << iNode
207  << std::endl;
208  return 0;
209  }
210  if (iOrder < 0) {
211  if ( --i < 0) i = 3;
212  if (m_pF[iFace].m_edge[i].v == 0) i = 2;
213  }
214  return (m_pF[iFace].m_edge[i].v > 0) ? 0 : m_pF[iFace].m_edge[i].f;
215 }
216 
217 HVNormal3D SbPolyhedron::FindNodeNormal(int iFace, int iNode) const
218 /***********************************************************************
219  * *
220  * Name: SbPolyhedron::FindNodeNormal Date: 22.11.99 *
221  * Author: E.Chernyaev Revised: *
222  * *
223  * Function: Find normal at given node *
224  * *
225  ***********************************************************************/
226 {
227  HVNormal3D normal = GetUnitNormal(iFace);
228  int k = iFace, iOrder = 1;
229 
230  for(;;) {
231  k = FindNeighbour(k, iNode, iOrder);
232  if (k == iFace) break;
233  if (k > 0) {
234  normal += GetUnitNormal(k);
235  }else{
236  if (iOrder < 0) break;
237  k = iFace;
238  iOrder = -iOrder;
239  }
240  }
241  normal.normalize();
242  return normal;
243 }
244 
246 /***********************************************************************
247  * *
248  * Name: SbPolyhedron::SetNumberOfRotationSteps Date: 24.06.97 *
249  * Author: J.Allison (Manchester University) Revised: *
250  * *
251  * Function: Set number of steps for whole circle *
252  * *
253  ***********************************************************************/
254 {
255  const int nMin = 3;
256  if (n < nMin) {
257  std::cerr
258  << "SbPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
259  << "number of steps per circle < " << nMin << "; forced to " << nMin
260  << std::endl;
262  }else{
264  }
265 }
266 
267 void SbPolyhedron::AllocateMemory(int Nvert, int Nface)
268 /***********************************************************************
269  * *
270  * Name: SbPolyhedron::AllocateMemory Date: 19.06.96 *
271  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
272  * *
273  * Function: Allocate memory for GEANT4 polyhedron *
274  * *
275  * Input: Nvert - number of nodes *
276  * Nface - number of faces *
277  * *
278  ***********************************************************************/
279 {
280  m_nvert = Nvert;
281  m_nface = Nface;
282  m_pV = new HVPoint3D[m_nvert+1];
283  m_pF = new SbFacet[m_nface+1];
284 }
285 
287 /***********************************************************************
288  * *
289  * Name: SbPolyhedron::CreatePrism Date: 15.07.96 *
290  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
291  * *
292  * Function: Set facets for a prism *
293  * *
294  ***********************************************************************/
295 {
296  enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
297 
298  m_pF[1] = SbFacet(1,LEFT, 4,BACK, 3,RIGHT, 2,FRONT);
299  m_pF[2] = SbFacet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
300  m_pF[3] = SbFacet(8,TOP, 7,RIGHT, 3,BOTTOM, 4,LEFT);
301  m_pF[4] = SbFacet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
302  m_pF[5] = SbFacet(6,TOP, 5,LEFT, 1,BOTTOM, 2,RIGHT);
303  m_pF[6] = SbFacet(5,FRONT, 6,RIGHT, 7,BACK, 8,LEFT);
304 }
305 
306 void SbPolyhedron::RotateEdge(int k1, int k2, double r1, double r2,
307  int v1, int v2, int vEdge,
308  bool ifWholeCircle, int ns, int &kface)
309 /***********************************************************************
310  * *
311  * Name: SbPolyhedron::RotateEdge Date: 05.12.96 *
312  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
313  * *
314  * Function: Create set of facets by rotation of an edge around Z-axis *
315  * *
316  * Input: k1, k2 - end vertices of the edge *
317  * r1, r2 - radiuses of the end vertices *
318  * v1, v2 - visibility of edges produced by rotation of the end *
319  * vertices *
320  * vEdge - visibility of the edge *
321  * ifWholeCircle - is true in case of whole circle rotation *
322  * ns - number of discrete steps *
323  * r[] - r-coordinates *
324  * kface - current free cell in the m_pF array *
325  * *
326  ***********************************************************************/
327 {
328  if (r1 == 0. && r2 == 0) return;
329 
330  int i;
331  int i1 = k1;
332  int i2 = k2;
333  int ii1 = ifWholeCircle ? i1 : i1+ns;
334  int ii2 = ifWholeCircle ? i2 : i2+ns;
335  int vv = ifWholeCircle ? vEdge : 1;
336 
337  if (ns == 1) {
338  if (r1 == 0.) {
339  m_pF[kface++] = SbFacet(i1,0, v2*i2,0, (i2+1),0);
340  }else if (r2 == 0.) {
341  m_pF[kface++] = SbFacet(i1,0, i2,0, v1*(i1+1),0);
342  }else{
343  m_pF[kface++] = SbFacet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),0);
344  }
345  }else{
346  if (r1 == 0.) {
347  m_pF[kface++] = SbFacet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0);
348  for (i2++,i=1; i<ns-1; i2++,i++) {
349  m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
350  }
351  m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
352  }else if (r2 == 0.) {
353  m_pF[kface++] = SbFacet(vv*i1,0, vEdge*i2,0, v1*(i1+1),0);
354  for (i1++,i=1; i<ns-1; i1++,i++) {
355  m_pF[kface++] = SbFacet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
356  }
357  m_pF[kface++] = SbFacet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
358  }else{
359  m_pF[kface++] = SbFacet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
360  for (i1++,i2++,i=1; i<ns-1; i1++,i2++,i++) {
361  m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
362  }
363  m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
364  }
365  }
366 }
367 
368 void SbPolyhedron::SetSideFacets(int ii[4], int vv[4],
369  int *kk, double *r,
370  double dphi, int ns, int &kface)
371 /***********************************************************************
372  * *
373  * Name: SbPolyhedron::SetSideFacets Date: 20.05.97 *
374  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
375  * *
376  * Function: Set side facets for the case of incomplete rotation *
377  * *
378  * Input: ii[4] - indeces of original verteces *
379  * vv[4] - visibility of edges *
380  * kk[] - indeces of nodes *
381  * r[] - radiuses *
382  * dphi - delta phi *
383  * ns - number of discrete steps *
384  * kface - current free cell in the m_pF array *
385  * *
386  ***********************************************************************/
387 {
388  int k1, k2, k3, k4;
389  if (fabs(dphi-M_PI) < perMillion) { // half a circle
390  for (int i=0; i<4; i++) {
391  k1 = ii[i];
392  k2 = (i == 3) ? ii[0] : ii[i+1];
393  if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
394  }
395  }
396 
397  if (ii[1] == ii[2]) {
398  k1 = kk[ii[0]];
399  k2 = kk[ii[2]];
400  k3 = kk[ii[3]];
401  m_pF[kface++] = SbFacet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
402  if (r[ii[0]] != 0.) k1 += ns;
403  if (r[ii[2]] != 0.) k2 += ns;
404  if (r[ii[3]] != 0.) k3 += ns;
405  m_pF[kface++] = SbFacet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
406  }else if (kk[ii[0]] == kk[ii[1]]) {
407  k1 = kk[ii[0]];
408  k2 = kk[ii[2]];
409  k3 = kk[ii[3]];
410  m_pF[kface++] = SbFacet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
411  if (r[ii[0]] != 0.) k1 += ns;
412  if (r[ii[2]] != 0.) k2 += ns;
413  if (r[ii[3]] != 0.) k3 += ns;
414  m_pF[kface++] = SbFacet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
415  }else if (kk[ii[2]] == kk[ii[3]]) {
416  k1 = kk[ii[0]];
417  k2 = kk[ii[1]];
418  k3 = kk[ii[2]];
419  m_pF[kface++] = SbFacet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
420  if (r[ii[0]] != 0.) k1 += ns;
421  if (r[ii[1]] != 0.) k2 += ns;
422  if (r[ii[2]] != 0.) k3 += ns;
423  m_pF[kface++] = SbFacet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
424  }else{
425  k1 = kk[ii[0]];
426  k2 = kk[ii[1]];
427  k3 = kk[ii[2]];
428  k4 = kk[ii[3]];
429  m_pF[kface++] = SbFacet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
430  if (r[ii[0]] != 0.) k1 += ns;
431  if (r[ii[1]] != 0.) k2 += ns;
432  if (r[ii[2]] != 0.) k3 += ns;
433  if (r[ii[3]] != 0.) k4 += ns;
434  m_pF[kface++] = SbFacet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
435  }
436 }
437 
438 void SbPolyhedron::RotateAroundZ(int nstep, double phi, double dphi,
439  int np1, int np2,
440  const double *z, double *r,
441  int nodeVis, int edgeVis)
442 /***********************************************************************
443  * *
444  * Name: SbPolyhedron::RotateAroundZ Date: 27.11.96 *
445  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
446  * *
447  * Function: Create SbPolyhedron for a solid produced by rotation of *
448  * two polylines around Z-axis *
449  * *
450  * Input: nstep - number of discrete steps, if 0 then default *
451  * phi - starting phi angle *
452  * dphi - delta phi *
453  * np1 - number of points in external polyline *
454  * (must be negative in case of closed polyline) *
455  * np2 - number of points in internal polyline (may be 1) *
456  * z[] - z-coordinates (+z >>> -z for both polylines) *
457  * r[] - r-coordinates *
458  * nodeVis - how to Draw edges joing consecutive positions of *
459  * node during rotation *
460  * edgeVis - how to Draw edges *
461  * *
462  ***********************************************************************/
463 {
464  static double wholeCircle = 2*M_PI;
465 
466  // S E T R O T A T I O N P A R A M E T E R S
467 
468  bool ifWholeCircle = (fabs(dphi-wholeCircle) < perMillion) ?
469  true : false;
470  double delPhi = ifWholeCircle ? wholeCircle : dphi;
471  int nSphi = (nstep > 0) ?
472  nstep : int(delPhi*GetNumberOfRotationSteps()/wholeCircle+.5);
473  if (nSphi == 0) nSphi = 1;
474  int nVphi = ifWholeCircle ? nSphi : nSphi+1;
475  bool ifClosed = np1 > 0 ? false : true;
476 
477  // C O U N T V E R T E C E S
478 
479  int absNp1 = iabs(np1);
480  int absNp2 = iabs(np2);
481  int i1beg = 0;
482  int i1end = absNp1-1;
483  int i2beg = absNp1;
484  int i2end = absNp1+absNp2-1;
485  int i, j, k;
486 
487  for(i=i1beg; i<=i2end; i++) {
488  if (fabs(r[i]) < perMillion) r[i] = 0.;
489  }
490 
491  j = 0; // external nodes
492  for (i=i1beg; i<=i1end; i++) {
493  j += (r[i] == 0.) ? 1 : nVphi;
494  }
495 
496  bool ifSide1 = false; // internal nodes
497  bool ifSide2 = false;
498 
499  if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
500  j += (r[i2beg] == 0.) ? 1 : nVphi;
501  ifSide1 = true;
502  }
503 
504  for(i=i2beg+1; i<i2end; i++) {
505  j += (r[i] == 0.) ? 1 : nVphi;
506  }
507 
508  if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
509  if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
510  ifSide2 = true;
511  }
512 
513  // C O U N T F A C E S
514 
515  k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi; // external faces
516 
517  if (absNp2 > 1) { // internal faces
518  for(i=i2beg; i<i2end; i++) {
519  if (r[i] > 0. || r[i+1] > 0.) k += nSphi;
520  }
521 
522  if (ifClosed) {
523  if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
524  }
525  }
526 
527  if (!ifClosed) { // side faces
528  if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
529  if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) k += nSphi;
530  }
531 
532  if (!ifWholeCircle) { // phi_side faces
533  k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
534  }
535 
536  // A L L O C A T E M E M O R Y
537 
538  AllocateMemory(j, k);
539 
540  // G E N E R A T E V E R T E C E S
541 
542  int *kk;
543  kk = new int[absNp1+absNp2];
544 
545  k = 1;
546  for(i=i1beg; i<=i1end; i++) {
547  kk[i] = k;
548  if (r[i] == 0.) { m_pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
549  }
550 
551  i = i2beg;
552  if (ifSide1) {
553  kk[i] = k;
554  if (r[i] == 0.) { m_pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
555  }else{
556  kk[i] = kk[i1beg];
557  }
558 
559  for(i=i2beg+1; i<i2end; i++) {
560  kk[i] = k;
561  if (r[i] == 0.) { m_pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
562  }
563 
564  if (absNp2 > 1) {
565  i = i2end;
566  if (ifSide2) {
567  kk[i] = k;
568  if (r[i] == 0.) m_pV[k] = HVPoint3D(0, 0, z[i]);
569  }else{
570  kk[i] = kk[i1end];
571  }
572  }
573 
574  double cosPhi, sinPhi;
575 
576  for(j=0; j<nVphi; j++) {
577  cosPhi = cos(phi+j*delPhi/nSphi);
578  sinPhi = sin(phi+j*delPhi/nSphi);
579  for(i=i1beg; i<=i2end; i++) {
580  if (r[i] != 0.) m_pV[kk[i]+j] = HVPoint3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
581  }
582  }
583 
584  // G E N E R A T E E X T E R N A L F A C E S
585 
586  int v1,v2;
587 
588  k = 1;
589  v2 = ifClosed ? nodeVis : 1;
590  for(i=i1beg; i<i1end; i++) {
591  v1 = v2;
592  if (!ifClosed && i == i1end-1) {
593  v2 = 1;
594  }else{
595  v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
596  }
597  RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
598  edgeVis, ifWholeCircle, nSphi, k);
599  }
600  if (ifClosed) {
601  RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
602  edgeVis, ifWholeCircle, nSphi, k);
603  }
604 
605  // G E N E R A T E I N T E R N A L F A C E S
606 
607  if (absNp2 > 1) {
608  v2 = ifClosed ? nodeVis : 1;
609  for(i=i2beg; i<i2end; i++) {
610  v1 = v2;
611  if (!ifClosed && i==i2end-1) {
612  v2 = 1;
613  }else{
614  v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
615  }
616  RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
617  edgeVis, ifWholeCircle, nSphi, k);
618  }
619  if (ifClosed) {
620  RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
621  edgeVis, ifWholeCircle, nSphi, k);
622  }
623  }
624 
625  // G E N E R A T E S I D E F A C E S
626 
627  if (!ifClosed) {
628  if (ifSide1) {
629  RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
630  -1, ifWholeCircle, nSphi, k);
631  }
632  if (ifSide2) {
633  RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
634  -1, ifWholeCircle, nSphi, k);
635  }
636  }
637 
638  // G E N E R A T E S I D E F A C E S for the case of incomplete circle
639 
640  if (!ifWholeCircle) {
641 
642  int ii[4], vv[4];
643 
644  if (ifClosed) {
645  for (i=i1beg; i<=i1end; i++) {
646  ii[0] = i;
647  ii[3] = (i == i1end) ? i1beg : i+1;
648  ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
649  ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
650  vv[0] = -1;
651  vv[1] = 1;
652  vv[2] = -1;
653  vv[3] = 1;
654  SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
655  }
656  }else{
657  for (i=i1beg; i<i1end; i++) {
658  ii[0] = i;
659  ii[3] = i+1;
660  ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
661  ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
662  vv[0] = (i == i1beg) ? 1 : -1;
663  vv[1] = 1;
664  vv[2] = (i == i1end-1) ? 1 : -1;
665  vv[3] = 1;
666  SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
667  }
668  }
669  }
670 
671  delete [] kk;
672 
673  if (k-1 != m_nface) {
674  std::cerr
675  << "Polyhedron::RotateAroundZ: number of generated faces ("
676  << k-1 << ") is not equal to the number of allocated faces ("
677  << m_nface << ")"
678  << std::endl;
679  }
680 }
681 
683 /***********************************************************************
684  * *
685  * Name: SbPolyhedron::SetReferences Date: 04.12.96 *
686  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
687  * *
688  * Function: For each edge set reference to neighbouring facet *
689  * *
690  ***********************************************************************/
691 {
692  if (m_nface <= 0) return;
693 
694  struct edgeListMember {
695  edgeListMember *next;
696  int v2;
697  int iface;
698  int iedge;
699  } *edgeList, *freeList, **headList;
700 
701 
702  // A L L O C A T E A N D I N I T I A T E L I S T S
703 
704  edgeList = new edgeListMember[2*m_nface];
705  headList = new edgeListMember*[m_nvert];
706 
707  int i;
708  for (i=0; i<m_nvert; i++) {
709  headList[i] = 0;
710  }
711  freeList = edgeList;
712  for (i=0; i<2*m_nface-1; i++) {
713  edgeList[i].next = &edgeList[i+1];
714  }
715  edgeList[2*m_nface-1].next = 0;
716 
717  // L O O P A L O N G E D G E S
718 
719  int iface, iedge, nedge, i1, i2, k1, k2;
720  edgeListMember *prev, *cur;
721 
722  for(iface=1; iface<=m_nface; iface++) {
723  nedge = (m_pF[iface].m_edge[3].v == 0) ? 3 : 4;
724  for (iedge=0; iedge<nedge; iedge++) {
725  i1 = iedge;
726  i2 = (iedge < nedge-1) ? iedge+1 : 0;
727  i1 = iabs(m_pF[iface].m_edge[i1].v);
728  i2 = iabs(m_pF[iface].m_edge[i2].v);
729  k1 = (i1 < i2) ? i1 : i2; // k1 = ::min(i1,i2);
730  k2 = (i1 > i2) ? i1 : i2; // k2 = ::max(i1,i2);
731 
732  // check head of the List corresponding to k1
733  cur = headList[k1];
734  if (cur == 0) {
735  headList[k1] = freeList;
736  freeList = freeList->next;
737  cur = headList[k1];
738  cur->next = 0;
739  cur->v2 = k2;
740  cur->iface = iface;
741  cur->iedge = iedge;
742  continue;
743  }
744 
745  if (cur->v2 == k2) {
746  headList[k1] = cur->next;
747  cur->next = freeList;
748  freeList = cur;
749  m_pF[iface].m_edge[iedge].f = cur->iface;
750  m_pF[cur->iface].m_edge[cur->iedge].f = iface;
751  i1 = (m_pF[iface].m_edge[iedge].v < 0) ? -1 : 1;
752  i2 = (m_pF[cur->iface].m_edge[cur->iedge].v < 0) ? -1 : 1;
753  if (i1 != i2) {
754  std::cerr
755  << "Polyhedron::SetReferences: different edge visibility "
756  << iface << "/" << iedge << "/"
757  << m_pF[iface].m_edge[iedge].v << " and "
758  << cur->iface << "/" << cur->iedge << "/"
759  << m_pF[cur->iface].m_edge[cur->iedge].v
760  << std::endl;
761  }
762  continue;
763  }
764 
765  // check List itself
766  for (;;) {
767  prev = cur;
768  cur = prev->next;
769  if (cur == 0) {
770  prev->next = freeList;
771  freeList = freeList->next;
772  cur = prev->next;
773  cur->next = 0;
774  cur->v2 = k2;
775  cur->iface = iface;
776  cur->iedge = iedge;
777  break;
778  }
779 
780  if (cur->v2 == k2) {
781  prev->next = cur->next;
782  cur->next = freeList;
783  freeList = cur;
784  m_pF[iface].m_edge[iedge].f = cur->iface;
785  m_pF[cur->iface].m_edge[cur->iedge].f = iface;
786  i1 = (m_pF[iface].m_edge[iedge].v < 0) ? -1 : 1;
787  i2 = (m_pF[cur->iface].m_edge[cur->iedge].v < 0) ? -1 : 1;
788  if (i1 != i2) {
789  std::cerr
790  << "Polyhedron::SetReferences: different edge visibility "
791  << iface << "/" << iedge << "/"
792  << m_pF[iface].m_edge[iedge].v << " and "
793  << cur->iface << "/" << cur->iedge << "/"
794  << m_pF[cur->iface].m_edge[cur->iedge].v
795  << std::endl;
796  }
797  break;
798  }
799  }
800  }
801  }
802 
803  // C H E C K T H A T A L L L I S T S A R E E M P T Y
804 
805  for (i=0; i<m_nvert; i++) {
806  if (headList[i] != 0) {
807  std::cerr
808  << "Polyhedron::SetReferences: List " << i << " is not empty"
809  << std::endl;
810  }
811  }
812 
813  // F R E E M E M O R Y
814 
815  delete [] edgeList;
816  delete [] headList;
817 }
818 
820 /***********************************************************************
821  * *
822  * Name: SbPolyhedron::InvertFacets Date: 01.12.99 *
823  * Author: E.Chernyaev Revised: *
824  * *
825  * Function: Invert the order of the nodes in the facets *
826  * *
827  ***********************************************************************/
828 {
829  if (m_nface <= 0) return;
830  int i, k, nnode, v[4],f[4];
831  for (i=1; i<=m_nface; i++) {
832  nnode = (m_pF[i].m_edge[3].v == 0) ? 3 : 4;
833  for (k=0; k<nnode; k++) {
834  v[k] = (k+1 == nnode) ? m_pF[i].m_edge[0].v : m_pF[i].m_edge[k+1].v;
835  if (v[k] * m_pF[i].m_edge[k].v < 0) v[k] = -v[k];
836  f[k] = m_pF[i].m_edge[k].f;
837  }
838  for (k=0; k<nnode; k++) {
839  m_pF[i].m_edge[nnode-1-k].v = v[k];
840  m_pF[i].m_edge[nnode-1-k].f = f[k];
841  }
842  }
843 }
844 
845 
846 // rbianchi change
847 //- original transform method
848 /*
849 SbPolyhedron & SbPolyhedron::Transform(const HVRotation & rotation, const HVVector3D & translation)
850 // ***********************************************************************
851 // * *
852 // * Name: SbPolyhedron::Transform Date: 01.12.99 *
853 // * Author: E.Chernyaev Revised: *
854 // * *
855 // * Function: Make transformation of the polyhedron *
856 // * *
857 // ***********************************************************************
858 {
859  if (m_nvert > 0) {
860  for (int i=1; i<=m_nvert; i++) {
861  HVVector3D tmp;
862  rotation.multVec(m_pV[i],tmp);
863  m_pV[i] = tmp+translation;
864  }
865 
866  // C H E C K D E T E R M I N A N T A N D
867  // 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
868 
869  HVVector3D x; rotation.multVec(HVVector3D(1,0,0),x);
870  HVVector3D y; rotation.multVec(HVVector3D(0,1,0),y);
871  HVVector3D z; rotation.multVec(HVVector3D(0,0,1),z);
872  if ((x.cross(y)).dot(z) < 0) InvertFacets();
873  }
874  return *this;
875 }
876 */
877 //- new transform method from the newest OpenScientist/HEPVis
880 ,const SbVec3d& translation
881 )
882 {
883  if (m_nvert > 0) {
884  SbVec3d tmp;
885  for (int i=1; i<=m_nvert; i++) {
886  rotation.multVec(m_pV[i],tmp);
887  m_pV[i] = tmp+translation;
888  }
889  SbVec3d x; rotation.multVec(SbVec3d(1,0,0),x);
890  SbVec3d y; rotation.multVec(SbVec3d(0,1,0),y);
891  SbVec3d z; rotation.multVec(SbVec3d(0,0,1),z);
892  if ((x.cross(y)).dot(z) < 0) InvertFacets();
893  }
894  return *this;
895 }
896 //---
897 
898 
899 
900 
901 
902 bool SbPolyhedron::GetNextVertexIndex(int &index, int &edgeFlag) const
903 /***********************************************************************
904  * *
905  * Name: SbPolyhedron::GetNextVertexIndex Date: 03.09.96 *
906  * Author: Yasuhide Sawada Revised: *
907  * *
908  * Function: *
909  * *
910  ***********************************************************************/
911 {
912  static int iFace = 1;
913  static int iQVertex = 0;
914  int vIndex = m_pF[iFace].m_edge[iQVertex].v;
915 
916  edgeFlag = (vIndex > 0) ? 1 : 0;
917  index = iabs(vIndex);
918 
919  if (iQVertex >= 3 || m_pF[iFace].m_edge[iQVertex+1].v == 0) {
920  iQVertex = 0;
921  if (++iFace > m_nface) iFace = 1;
922  return false; // Last Edge
923  }else{
924  ++iQVertex;
925  return true; // not Last Edge
926  }
927 }
928 
930 /***********************************************************************
931  * *
932  * Name: SbPolyhedron::GetVertex Date: 03.09.96 *
933  * Author: Yasuhide Sawada Revised: 17.11.99 *
934  * *
935  * Function: Get vertex of the index. *
936  * *
937  ***********************************************************************/
938 {
939  if (index <= 0 || index > m_nvert) {
940  std::cerr
941  << "SbPolyhedron::GetVertex: irrelevant index " << index
942  << std::endl;
943  return HVPoint3D();
944  }
945  return m_pV[index];
946 }
947 
948 
949 // rbianchi - 14.12.2012
950 const HVPoint3D& SbPolyhedron::GetVertexFast(int index) const { //G.Barrand
951  return m_pV[index];
952 }
953 //---
954 
955 
956 bool
958 /***********************************************************************
959  * *
960  * Name: SbPolyhedron::GetNextVertex Date: 22.07.96 *
961  * Author: John Allison Revised: *
962  * *
963  * Function: Get vertices of the quadrilaterals in order for each *
964  * face in face order. Returns false when finished each *
965  * face. *
966  * *
967  ***********************************************************************/
968 {
969  int index;
970  bool rep = GetNextVertexIndex(index, edgeFlag);
971  vertex = m_pV[index];
972  return rep;
973 }
974 
976  HVNormal3D &normal) const
977 /***********************************************************************
978  * *
979  * Name: SbPolyhedron::GetNextVertex Date: 26.11.99 *
980  * Author: E.Chernyaev Revised: *
981  * *
982  * Function: Get vertices with normals of the quadrilaterals in order *
983  * for each face in face order. *
984  * Returns false when finished each face. *
985  * *
986  ***********************************************************************/
987 {
988  static int iFace = 1;
989  static int iNode = 0;
990 
991  if (m_nface == 0) return false; // empty polyhedron
992 
993  int k = m_pF[iFace].m_edge[iNode].v;
994  if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
995  vertex = m_pV[k];
996  normal = FindNodeNormal(iFace,k);
997  if (iNode >= 3 || m_pF[iFace].m_edge[iNode+1].v == 0) {
998  iNode = 0;
999  if (++iFace > m_nface) iFace = 1;
1000  return false; // last node
1001  }else{
1002  ++iNode;
1003  return true; // not last node
1004  }
1005 }
1006 
1007 bool SbPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag,
1008  int &iface1, int &iface2) const
1009 /***********************************************************************
1010  * *
1011  * Name: SbPolyhedron::GetNextEdgeIndeces Date: 30.09.96 *
1012  * Author: E.Chernyaev Revised: 17.11.99 *
1013  * *
1014  * Function: Get indeces of the next edge together with indeces of *
1015  * of the faces which share the edge. *
1016  * Returns false when the last edge. *
1017  * *
1018  ***********************************************************************/
1019 {
1020  static int iFace = 1;
1021  static int iQVertex = 0;
1022  static int iOrder = 1;
1023  int k1, k2, kflag, kface1, kface2;
1024 
1025  if (iFace == 1 && iQVertex == 0) {
1026  k2 = m_pF[m_nface].m_edge[0].v;
1027  k1 = m_pF[m_nface].m_edge[3].v;
1028  if (k1 == 0) k1 = m_pF[m_nface].m_edge[2].v;
1029  if (iabs(k1) > iabs(k2)) iOrder = -1;
1030  }
1031 
1032  do {
1033  k1 = m_pF[iFace].m_edge[iQVertex].v;
1034  kflag = k1;
1035  k1 = iabs(k1);
1036  kface1 = iFace;
1037  kface2 = m_pF[iFace].m_edge[iQVertex].f;
1038  if (iQVertex >= 3 || m_pF[iFace].m_edge[iQVertex+1].v == 0) {
1039  iQVertex = 0;
1040  k2 = iabs(m_pF[iFace].m_edge[iQVertex].v);
1041  iFace++;
1042  }else{
1043  iQVertex++;
1044  k2 = iabs(m_pF[iFace].m_edge[iQVertex].v);
1045  }
1046  } while (iOrder*k1 > iOrder*k2);
1047 
1048  i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1049  iface1 = kface1; iface2 = kface2;
1050 
1051  if (iFace > m_nface) {
1052  iFace = 1; iOrder = 1;
1053  return false;
1054  }else{
1055  return true;
1056  }
1057 }
1058 
1059 bool
1060 SbPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag) const
1061 /***********************************************************************
1062  * *
1063  * Name: SbPolyhedron::GetNextEdgeIndeces Date: 17.11.99 *
1064  * Author: E.Chernyaev Revised: *
1065  * *
1066  * Function: Get indeces of the next edge. *
1067  * Returns false when the last edge. *
1068  * *
1069  ***********************************************************************/
1070 {
1071  int kface1, kface2;
1072  return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
1073 }
1074 
1075 bool
1077  HVPoint3D &p2,
1078  int &edgeFlag) const
1079 /***********************************************************************
1080  * *
1081  * Name: SbPolyhedron::GetNextEdge Date: 30.09.96 *
1082  * Author: E.Chernyaev Revised: *
1083  * *
1084  * Function: Get next edge. *
1085  * Returns false when the last edge. *
1086  * *
1087  ***********************************************************************/
1088 {
1089  int i1,i2;
1090  bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
1091  p1 = m_pV[i1];
1092  p2 = m_pV[i2];
1093  return rep;
1094 }
1095 
1096 bool
1098  int &edgeFlag, int &iface1, int &iface2) const
1099 /***********************************************************************
1100  * *
1101  * Name: SbPolyhedron::GetNextEdge Date: 17.11.99 *
1102  * Author: E.Chernyaev Revised: *
1103  * *
1104  * Function: Get next edge with indeces of the faces which share *
1105  * the edge. *
1106  * Returns false when the last edge. *
1107  * *
1108  ***********************************************************************/
1109 {
1110  int i1,i2;
1111  bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
1112  p1 = m_pV[i1];
1113  p2 = m_pV[i2];
1114  return rep;
1115 }
1116 
1117 void SbPolyhedron::GetFacet(int iFace, int &n, int *iNodes,
1118  int *edgeFlags, int *iFaces) const
1119 /***********************************************************************
1120  * *
1121  * Name: SbPolyhedron::GetFacet Date: 15.12.99 *
1122  * Author: E.Chernyaev Revised: *
1123  * *
1124  * Function: Get face by index *
1125  * *
1126  ***********************************************************************/
1127 {
1128  if (iFace < 1 || iFace > m_nface) {
1129  std::cerr
1130  << "SbPolyhedron::GetFacet: irrelevant index " << iFace
1131  << std::endl;
1132  n = 0;
1133  }else{
1134  int i, k;
1135  for (i=0; i<4; i++) {
1136  k = m_pF[iFace].m_edge[i].v;
1137  if (k == 0) break;
1138  if (iFaces != 0) iFaces[i] = m_pF[iFace].m_edge[i].f;
1139  if (k > 0) {
1140  iNodes[i] = k;
1141  if (edgeFlags != 0) edgeFlags[i] = 1;
1142  }else{
1143  iNodes[i] = -k;
1144  if (edgeFlags != 0) edgeFlags[i] = -1;
1145  }
1146  }
1147  n = i;
1148  }
1149 }
1150 
1151 void SbPolyhedron::GetFacet(int index, int &n, HVPoint3D *nodes,
1152  int *edgeFlags, HVNormal3D *normals) const
1153 /***********************************************************************
1154  * *
1155  * Name: SbPolyhedron::GetFacet Date: 17.11.99 *
1156  * Author: E.Chernyaev Revised: *
1157  * *
1158  * Function: Get face by index *
1159  * *
1160  ***********************************************************************/
1161 {
1162  int iNodes[4];
1163  GetFacet(index, n, iNodes, edgeFlags);
1164  if (n != 0) {
1165  for (int i=0; i<4; i++) {
1166  nodes[i] = m_pV[iNodes[i]];
1167  if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1168  }
1169  }
1170 }
1171 
1172 bool
1174  int *edgeFlags, HVNormal3D *normals) const
1175 /***********************************************************************
1176  * *
1177  * Name: SbPolyhedron::GetNextFacet Date: 19.11.99 *
1178  * Author: E.Chernyaev Revised: *
1179  * *
1180  * Function: Get next face with normals of unit length at the nodes. *
1181  * Returns false when finished all faces. *
1182  * *
1183  ***********************************************************************/
1184 {
1185  static int iFace = 1;
1186 
1187  if (edgeFlags == 0) {
1188  GetFacet(iFace, n, nodes);
1189  }else if (normals == 0) {
1190  GetFacet(iFace, n, nodes, edgeFlags);
1191  }else{
1192  GetFacet(iFace, n, nodes, edgeFlags, normals);
1193  }
1194 
1195  if (++iFace > m_nface) {
1196  iFace = 1;
1197  return false;
1198  }else{
1199  return true;
1200  }
1201 }
1202 
1204 /***********************************************************************
1205  * *
1206  * Name: SbPolyhedron::GetNormal Date: 19.11.99 *
1207  * Author: E.Chernyaev Revised: *
1208  * *
1209  * Function: Get normal of the face given by index *
1210  * *
1211  ***********************************************************************/
1212 {
1213  if (iFace < 1 || iFace > m_nface) {
1214  std::cerr
1215  << "SbPolyhedron::GetNormal: irrelevant index " << iFace
1216  << std::endl;
1217  return HVNormal3D();
1218  }
1219 
1220  int i0 = iabs(m_pF[iFace].m_edge[0].v);
1221  int i1 = iabs(m_pF[iFace].m_edge[1].v);
1222  int i2 = iabs(m_pF[iFace].m_edge[2].v);
1223  int i3 = iabs(m_pF[iFace].m_edge[3].v);
1224  if (i3 == 0) i3 = i0;
1225  return (m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1]);
1226 }
1227 
1229 /***********************************************************************
1230  * *
1231  * Name: SbPolyhedron::GetNormal Date: 19.11.99 *
1232  * Author: E.Chernyaev Revised: *
1233  * *
1234  * Function: Get unit normal of the face given by index *
1235  * *
1236  ***********************************************************************/
1237 {
1238  if (iFace < 1 || iFace > m_nface) {
1239  std::cerr
1240  << "SbPolyhedron::GetUnitNormal: irrelevant index " << iFace
1241  << std::endl;
1242  return HVNormal3D();
1243  }
1244 
1245  int i0 = iabs(m_pF[iFace].m_edge[0].v);
1246  int i1 = iabs(m_pF[iFace].m_edge[1].v);
1247  int i2 = iabs(m_pF[iFace].m_edge[2].v);
1248  int i3 = iabs(m_pF[iFace].m_edge[3].v);
1249  if (i3 == 0) i3 = i0;
1250  HVNormal3D nm = (m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1]);
1251  nm.normalize();
1252  return nm;
1253 }
1254 
1256 /***********************************************************************
1257  * *
1258  * Name: SbPolyhedron::GetNextNormal Date: 22.07.96 *
1259  * Author: John Allison Revised: 19.11.99 *
1260  * *
1261  * Function: Get normals of each face in face order. Returns false *
1262  * when finished all faces. *
1263  * *
1264  ***********************************************************************/
1265 {
1266  static int iFace = 1;
1267  normal = GetNormal(iFace);
1268  if (++iFace > m_nface) {
1269  iFace = 1;
1270  return false;
1271  }else{
1272  return true;
1273  }
1274 }
1275 
1277 /***********************************************************************
1278  * *
1279  * Name: SbPolyhedron::GetNextUnitNormal Date: 16.09.96 *
1280  * Author: E.Chernyaev Revised: *
1281  * *
1282  * Function: Get normals of unit length of each face in face order. *
1283  * Returns false when finished all faces. *
1284  * *
1285  ***********************************************************************/
1286 {
1287  bool rep = GetNextNormal(normal);
1288  normal.normalize();
1289  return rep;
1290 }
1291 
1293 /***********************************************************************
1294  * *
1295  * Name: SbPolyhedron::GetSurfaceArea Date: 25.05.01 *
1296  * Author: E.Chernyaev Revised: *
1297  * *
1298  * Function: Returns area of the surface of the polyhedron. *
1299  * *
1300  ***********************************************************************/
1301 {
1302  double s = 0.;
1303  for (int iFace=1; iFace<=m_nface; iFace++) {
1304  int i0 = iabs(m_pF[iFace].m_edge[0].v);
1305  int i1 = iabs(m_pF[iFace].m_edge[1].v);
1306  int i2 = iabs(m_pF[iFace].m_edge[2].v);
1307  int i3 = iabs(m_pF[iFace].m_edge[3].v);
1308  if (i3 == 0) i3 = i0;
1309  s += ((m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1])).length();
1310  }
1311  return s/2.;
1312 }
1313 
1315 /***********************************************************************
1316  * *
1317  * Name: SbPolyhedron::GetVolume Date: 25.05.01 *
1318  * Author: E.Chernyaev Revised: *
1319  * *
1320  * Function: Returns volume of the polyhedron. *
1321  * *
1322  ***********************************************************************/
1323 {
1324  double v = 0.;
1325  for (int iFace=1; iFace<=m_nface; iFace++) {
1326  int i0 = iabs(m_pF[iFace].m_edge[0].v);
1327  int i1 = iabs(m_pF[iFace].m_edge[1].v);
1328  int i2 = iabs(m_pF[iFace].m_edge[2].v);
1329  int i3 = iabs(m_pF[iFace].m_edge[3].v);
1330  HVPoint3D g;
1331  if (i3 == 0) {
1332  i3 = i0;
1333  g = (m_pV[i0]+m_pV[i1]+m_pV[i2]) * (1.0f/3.0f);
1334  }else{
1335  g = (m_pV[i0]+m_pV[i1]+m_pV[i2]+m_pV[i3]) * 0.25f;
1336  }
1337  v += ((m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1])).dot(g);
1338  }
1339  return v/6.;
1340 }
1341 
1342 SbPolyhedronTrd2::SbPolyhedronTrd2(double Dx1, double Dx2,
1343  double Dy1, double Dy2,
1344  double Dz)
1345 /***********************************************************************
1346  * *
1347  * Name: SbPolyhedronTrd2 Date: 22.07.96 *
1348  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1349  * *
1350  * Function: Create GEANT4 TRD2-trapezoid *
1351  * *
1352  * Input: Dx1 - half-length along X at -Dz 8----7 *
1353  * Dx2 - half-length along X ay +Dz 5----6 ! *
1354  * Dy1 - half-length along Y ay -Dz ! 4-!--3 *
1355  * Dy2 - half-length along Y ay +Dz 1----2 *
1356  * Dz - half-length along Z *
1357  * *
1358  ***********************************************************************/
1359 {
1360  AllocateMemory(8,6);
1361 
1362  m_pV[1] = HVPoint3D(-Dx1,-Dy1,-Dz);
1363  m_pV[2] = HVPoint3D( Dx1,-Dy1,-Dz);
1364  m_pV[3] = HVPoint3D( Dx1, Dy1,-Dz);
1365  m_pV[4] = HVPoint3D(-Dx1, Dy1,-Dz);
1366  m_pV[5] = HVPoint3D(-Dx2,-Dy2, Dz);
1367  m_pV[6] = HVPoint3D( Dx2,-Dy2, Dz);
1368  m_pV[7] = HVPoint3D( Dx2, Dy2, Dz);
1369  m_pV[8] = HVPoint3D(-Dx2, Dy2, Dz);
1370 
1371  CreatePrism();
1372 }
1373 
1375 
1376 SbPolyhedronTrd1::SbPolyhedronTrd1(double Dx1, double Dx2,
1377  double Dy, double Dz)
1378  : SbPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1379 
1381 
1382 SbPolyhedronBox::SbPolyhedronBox(double Dx, double Dy, double Dz)
1383  : SbPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1384 
1386 
1388  double Theta,
1389  double Phi,
1390  double Dy1,
1391  double Dx1,
1392  double Dx2,
1393  double Alp1,
1394  double Dy2,
1395  double Dx3,
1396  double Dx4,
1397  double Alp2)
1398 /***********************************************************************
1399  * *
1400  * Name: SbPolyhedronTrap Date: 20.11.96 *
1401  * Author: E.Chernyaev Revised: *
1402  * *
1403  * Function: Create GEANT4 TRAP-trapezoid *
1404  * *
1405  * Input: DZ - half-length in Z *
1406  * Theta,Phi - polar angles of the line joining centres of the *
1407  * faces at Z=-Dz and Z=+Dz *
1408  * Dy1 - half-length in Y of the face at Z=-Dz *
1409  * Dx1 - half-length in X of low edge of the face at Z=-Dz *
1410  * Dx2 - half-length in X of top edge of the face at Z=-Dz *
1411  * Alp1 - angle between Y-axis and the median joining top and *
1412  * low edges of the face at Z=-Dz *
1413  * Dy2 - half-length in Y of the face at Z=+Dz *
1414  * Dx3 - half-length in X of low edge of the face at Z=+Dz *
1415  * Dx4 - half-length in X of top edge of the face at Z=+Dz *
1416  * Alp2 - angle between Y-axis and the median joining top and *
1417  * low edges of the face at Z=+Dz *
1418  * *
1419  ***********************************************************************/
1420 {
1421  double DzTthetaCphi = Dz*tan(Theta)*cos(Phi);
1422  double DzTthetaSphi = Dz*tan(Theta)*sin(Phi);
1423  double Dy1Talp1 = Dy1*tan(Alp1);
1424  double Dy2Talp2 = Dy2*tan(Alp2);
1425 
1426  AllocateMemory(8,6);
1427 
1428  m_pV[1] = HVPoint3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
1429  m_pV[2] = HVPoint3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
1430  m_pV[3] = HVPoint3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
1431  m_pV[4] = HVPoint3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
1432  m_pV[5] = HVPoint3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
1433  m_pV[6] = HVPoint3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
1434  m_pV[7] = HVPoint3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
1435  m_pV[8] = HVPoint3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
1436 
1437  CreatePrism();
1438 }
1439 
1441 
1442 SbPolyhedronPara::SbPolyhedronPara(double Dx, double Dy, double Dz,
1443  double Alpha, double Theta,
1444  double Phi)
1445  : SbPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
1446 
1448 
1450  double Rmx1,
1451  double Rmn2,
1452  double Rmx2,
1453  double Dz,
1454  double Phi1,
1455  double Dphi)
1456 /***********************************************************************
1457  * *
1458  * Name: SbPolyhedronCons::SbPolyhedronCons Date: 15.12.96 *
1459  * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
1460  * *
1461  * Function: Constructor for CONS, TUBS, CONE, TUBE *
1462  * *
1463  * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
1464  * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
1465  * Dz - half length in Z *
1466  * Phi1 - starting angle of the segment *
1467  * Dphi - segment range *
1468  * *
1469  ***********************************************************************/
1470 {
1471  static double wholeCircle=2*M_PI;
1472 
1473  // C H E C K I N P U T P A R A M E T E R S
1474 
1475  int k = 0;
1476  if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
1477  if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
1478  if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
1479 
1480  // We can get this from the tracking geometry. Don't complain.
1481  if (Dz == 0) return;
1482 
1483  if (Dz < 0.) k += 2;
1484 
1485  double phi1, phi2, dphi;
1486  if (Dphi < 0.) {
1487  phi2 = Phi1; phi1 = phi2 - Dphi;
1488  }else if (Dphi == 0.) {
1489  phi1 = Phi1; phi2 = phi1 + wholeCircle;
1490  }else{
1491  phi1 = Phi1; phi2 = phi1 + Dphi;
1492  }
1493  dphi = phi2 - phi1;
1494  if (fabs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1495  if (dphi > wholeCircle) k += 4;
1496 
1497  if (k != 0) {
1498  std::cerr << "SbPolyhedronCone(s)/Tube(s): error in input parameters";
1499  if ((k & 1) != 0) std::cerr << " (radiuses)";
1500  if ((k & 2) != 0) std::cerr << " (half-length)";
1501  if ((k & 4) != 0) std::cerr << " (angles)";
1502  std::cerr << std::endl;
1503  std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
1504  std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
1505  std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
1506  << std::endl;
1507  return;
1508  }
1509 
1510  // P R E P A R E T W O P O L Y L I N E S
1511 
1512  double zz[4], rr[4];
1513  zz[0] = Dz;
1514  zz[1] = -Dz;
1515  zz[2] = Dz;
1516  zz[3] = -Dz;
1517  rr[0] = Rmx2;
1518  rr[1] = Rmx1;
1519  rr[2] = Rmn2;
1520  rr[3] = Rmn1;
1521 
1522  // R O T A T E P O L Y L I N E S
1523 
1524  RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1);
1525  SetReferences();
1526 }
1527 
1529 
1530 SbPolyhedronCone::SbPolyhedronCone(double Rmn1, double Rmx1,
1531  double Rmn2, double Rmx2,
1532  double Dz) :
1533  SbPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
1534 
1536 
1538  double Dz,
1539  double Phi1, double Dphi)
1540  : SbPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
1541 
1543 
1545  double Dz)
1546  : SbPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
1547 
1549 
1551  double dphi,
1552  int npdv,
1553  int nz,
1554  const double *z,
1555  const double *rmin,
1556  const double *rmax)
1557 /***********************************************************************
1558  * *
1559  * Name: SbPolyhedronPgon Date: 09.12.96 *
1560  * Author: E.Chernyaev Revised: *
1561  * *
1562  * Function: Constructor of polyhedron for PGON, PCON *
1563  * *
1564  * Input: phi - initial phi *
1565  * dphi - delta phi *
1566  * npdv - number of steps along phi *
1567  * nz - number of z-planes (at least two) *
1568  * z[] - z coordinates of the slices *
1569  * rmin[] - smaller r at the slices *
1570  * rmax[] - bigger r at the slices *
1571  * *
1572  ***********************************************************************/
1573 {
1574  // C H E C K I N P U T P A R A M E T E R S
1575 
1576  if (dphi <= 0. || dphi > 2*M_PI) {
1577  std::cerr
1578  << "SbPolyhedronPgon/Pcon: wrong delta phi = " << dphi
1579  << std::endl;
1580  return;
1581  }
1582 
1583  if (nz < 2) {
1584  std::cerr
1585  << "SbPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
1586  << std::endl;
1587  return;
1588  }
1589 
1590  if (npdv < 0) {
1591  std::cerr
1592  << "SbPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
1593  << std::endl;
1594  return;
1595  }
1596 
1597  int i;
1598  for (i=0; i<nz; i++) {
1599  if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
1600  std::cerr
1601  << "SbPolyhedronPgon: error in radiuses rmin[" << i << "]="
1602  << rmin[i] << " rmax[" << i << "]=" << rmax[i]
1603  << std::endl;
1604  return;
1605  }
1606  }
1607 
1608  // P R E P A R E T W O P O L Y L I N E S
1609 
1610  double *zz, *rr;
1611  zz = new double[2*nz];
1612  rr = new double[2*nz];
1613 
1614  if (z[0] > z[nz-1]) {
1615  for (i=0; i<nz; i++) {
1616  zz[i] = z[i];
1617  rr[i] = rmax[i];
1618  zz[i+nz] = z[i];
1619  rr[i+nz] = rmin[i];
1620  }
1621  }else{
1622  for (i=0; i<nz; i++) {
1623  zz[i] = z[nz-i-1];
1624  rr[i] = rmax[nz-i-1];
1625  zz[i+nz] = z[nz-i-1];
1626  rr[i+nz] = rmin[nz-i-1];
1627  }
1628  }
1629 
1630  // R O T A T E P O L Y L I N E S
1631 
1632  RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1);
1633  SetReferences();
1634 
1635  delete [] zz;
1636  delete [] rr;
1637 }
1638 
1640 
1641 SbPolyhedronPcon::SbPolyhedronPcon(double phi, double dphi, int nz,
1642  const double *z,
1643  const double *rmin,
1644  const double *rmax)
1645  : SbPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
1646 
1648 
1650  double phi, double dphi,
1651  double the, double dthe)
1652 /***********************************************************************
1653  * *
1654  * Name: SbPolyhedronSphere Date: 11.12.96 *
1655  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1656  * *
1657  * Function: Constructor of polyhedron for SPHERE *
1658  * *
1659  * Input: rmin - internal radius *
1660  * rmax - external radius *
1661  * phi - initial phi *
1662  * dphi - delta phi *
1663  * the - initial theta *
1664  * dthe - delta theta *
1665  * *
1666  ***********************************************************************/
1667 {
1668  // C H E C K I N P U T P A R A M E T E R S
1669 
1670  if (dphi <= 0. || dphi > 2*M_PI) {
1671  std::cerr
1672  << "SbPolyhedronSphere: wrong delta phi = " << dphi
1673  << std::endl;
1674  return;
1675  }
1676 
1677  if (the < 0. || the > M_PI) {
1678  std::cerr
1679  << "SbPolyhedronSphere: wrong theta = " << the
1680  << std::endl;
1681  return;
1682  }
1683 
1684  if (dthe <= 0. || dthe > M_PI) {
1685  std::cerr
1686  << "SbPolyhedronSphere: wrong delta theta = " << dthe
1687  << std::endl;
1688  return;
1689  }
1690 
1691  if ( (the+dthe >= M_PI) && (the+dthe < M_PI + 2*DBL_EPSILON) )
1692  dthe = M_PI - the; //G.Barrand : coming from LHCb/S.Ponce.
1693 
1694  if (the+dthe > M_PI) {
1695  std::cerr
1696  << "SbPolyhedronSphere: wrong theta + delta theta = "
1697  << the << " " << dthe
1698  << std::endl;
1699  return;
1700  }
1701 
1702  if (rmin < 0. || rmin >= rmax) {
1703  std::cerr
1704  << "SbPolyhedronSphere: error in radiuses"
1705  << " rmin=" << rmin << " rmax=" << rmax
1706  << std::endl;
1707  return;
1708  }
1709 
1710  // P R E P A R E T W O P O L Y L I N E S
1711 
1712  int ns = (GetNumberOfRotationSteps() + 1) / 2;
1713  int np1 = int(dthe*ns/M_PI+.5) + 1;
1714  if (np1 <= 1) np1 = 2;
1715  int np2 = rmin < perMillion ? 1 : np1;
1716 
1717  double *zz, *rr;
1718  zz = new double[np1+np2];
1719  rr = new double[np1+np2];
1720 
1721  double a = dthe/(np1-1);
1722  double cosa, sina;
1723  for (int i=0; i<np1; i++) {
1724  cosa = cos(the+i*a);
1725  sina = sin(the+i*a);
1726  zz[i] = rmax*cosa;
1727  rr[i] = rmax*sina;
1728  if (np2 > 1) {
1729  zz[i+np1] = rmin*cosa;
1730  rr[i+np1] = rmin*sina;
1731  }
1732  }
1733  if (np2 == 1) {
1734  zz[np1] = 0.;
1735  rr[np1] = 0.;
1736  }
1737 
1738  // R O T A T E P O L Y L I N E S
1739 
1740  RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
1741  SetReferences();
1742 
1743  delete [] zz;
1744  delete [] rr;
1745 }
1746 
1748 
1750  double rmax,
1751  double rtor,
1752  double phi,
1753  double dphi)
1754 /***********************************************************************
1755  * *
1756  * Name: SbPolyhedronTorus Date: 11.12.96 *
1757  * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1758  * *
1759  * Function: Constructor of polyhedron for TORUS *
1760  * *
1761  * Input: rmin - internal radius *
1762  * rmax - external radius *
1763  * rtor - radius of torus *
1764  * phi - initial phi *
1765  * dphi - delta phi *
1766  * *
1767  ***********************************************************************/
1768 {
1769  // C H E C K I N P U T P A R A M E T E R S
1770 
1771  if (dphi <= 0. || dphi > 2*M_PI) {
1772  std::cerr
1773  << "SbPolyhedronTorus: wrong delta phi = " << dphi
1774  << std::endl;
1775  return;
1776  }
1777 
1778  if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
1779  std::cerr
1780  << "SbPolyhedronTorus: error in radiuses"
1781  << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
1782  << std::endl;
1783  return;
1784  }
1785 
1786  // P R E P A R E T W O P O L Y L I N E S
1787 
1788  int np1 = GetNumberOfRotationSteps();
1789  int np2 = rmin < perMillion ? 1 : np1;
1790 
1791  //double *zz(nullptr), *rr;
1792  const auto totNumPts = np1+np2;
1793  auto zz = new double[totNumPts]{};
1794  auto rr = new double[totNumPts]{};
1795  for (unsigned int i(0);i!=(unsigned int) totNumPts;++i){
1796  rr[i]=0.0;
1797  zz[i]=0.0;
1798  }
1799 
1800  double a = 2*M_PI/np1;
1801  double cosa, sina;
1802  for (int i=0; i<np1; i++) {
1803  cosa = cos(i*a);
1804  sina = sin(i*a);
1805  zz[i] = rmax*cosa;
1806  rr[i] = rtor+rmax*sina;
1807  if (np2 > 1) {
1808  zz[i+np1] = rmin*cosa;
1809  rr[i+np1] = rtor+rmin*sina;
1810  }
1811  }
1812  if (np2 == 1) {
1813  zz[np1] = 0.;
1814  rr[np1] = rtor;
1815  np2 = -1;
1816  }
1817 
1818  // R O T A T E P O L Y L I N E S
1819 
1820  RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1);
1821  SetReferences();
1822 
1823  delete [] zz;
1824  delete [] rr;
1825 }
1826 
1828 
1830 /***********************************************************************
1831  * *
1832  * Name: SbPolyhedron::s_numberOfRotationSteps Date: 24.06.97 *
1833  * Author: J.Allison (Manchester University) Revised: *
1834  * *
1835  * Function: Number of steps for whole circle *
1836  * *
1837  ***********************************************************************/
1838 
1839 #include "BooleanProcessor.h"
1840 static HEPVis_BooleanProcessor processor;
1841 
1843 /***********************************************************************
1844  * *
1845  * Name: SbPolyhedron::add Date: 19.03.00 *
1846  * Author: E.Chernyaev Revised: *
1847  * *
1848  * Function: Boolean "union" of two polyhedra *
1849  * *
1850  ***********************************************************************/
1851 {
1852  // rbianchi change
1853  //return processor.execute(OP_UNION, *this, p);
1854  HEPVis_BooleanProcessor processor; //G.Barrand
1855  int err;
1856  return processor.execute(OP_UNION, *this, p, err);
1857  //---
1858 }
1859 
1861 /***********************************************************************
1862  * *
1863  * Name: SbPolyhedron::intersect Date: 19.03.00 *
1864  * Author: E.Chernyaev Revised: *
1865  * *
1866  * Function: Boolean "intersection" of two polyhedra *
1867  * *
1868  ***********************************************************************/
1869 {
1870  // rbianchi change
1871  //return processor.execute(OP_INTERSECTION, *this, p);
1872  HEPVis_BooleanProcessor processor; //G.Barrand
1873  int err;
1874  return processor.execute(OP_INTERSECTION, *this, p, err);
1875  //---
1876 }
1877 
1879 /***********************************************************************
1880  * *
1881  * Name: SbPolyhedron::add Date: 19.03.00 *
1882  * Author: E.Chernyaev Revised: *
1883  * *
1884  * Function: Boolean "subtraction" of "p" from "this" *
1885  * *
1886  ***********************************************************************/
1887 {
1888  // rbianchi change
1889  //return processor.execute(OP_SUBTRACTION, *this, p);
1890  HEPVis_BooleanProcessor processor; //G.Barrand
1891  int err;
1892  return processor.execute(OP_SUBTRACTION, *this, p, err);
1893  //---
1894 }
1895 
1896 
1901 
1902 #include "PolygonTriangulator.h"
1903 #include <map>
1904 #include <set>
1905 
1907 public:
1908 
1909  typedef std::pair<int,int> Edge;
1912 
1913  // Methods:
1917  ~Internals() { delete poly; }
1918 
1919  void setData(const std::vector<double> * xx,const std::vector<double>* yy, const double& dz);
1920  void triangulate();
1922  void addExtraVertices();
1924  void defineFacesTopology();
1925  //Data:
1926  double dz;
1927  const std::vector<double> * x;
1928  const std::vector<double> * y;
1929  unsigned n;
1930  unsigned ntriangles;
1932 
1933  //The triangulation:
1935  //Various helper maps, sets, ...:
1936  std::set<Edge> edges_internal;
1937  std::set<Edge> edges_external;
1938  std::map<Edge,const Triangle*> neighbourmap;//In this map the Edges are oriented!
1939  std::set<Edge> edges_with_extra_vertex;
1942  std::map<Edge,unsigned> edgewithextravertex_2_id;
1944 
1945  //Helper methods:
1946 
1947  bool isinternaledge(const Edge& edge, const std::map<Edge,std::vector<Triangles::const_iterator> >& edge2triangles_map) const;
1948  Edge GetEdge(const Triangle * tr, const unsigned& i, const bool& oriented=false) const;//If oriented, 2-4 and 4-2 are different. If not, 4-2 is returned as 2-4.
1949  void setupface(const unsigned& face_id, const std::vector<unsigned>&v) const;
1950  void setupExternalFace(const unsigned& id_extedge,const unsigned& prev_edgeface,const unsigned& next_edgeface,
1951  const unsigned& idfirstvertex,const unsigned& idsecondvertex,
1952  const unsigned& id_triangleface,const unsigned& n, std::vector<unsigned>& faceinfo) const;
1953 
1954  void getSurroundingValues(const std::vector<unsigned>& numbercycle, const unsigned& centralValue,
1955  unsigned& prevValue, unsigned& nextValue) const;
1956 
1957  unsigned top2bottomvertexid(const unsigned& topid) const;
1958  unsigned top2bottomfaceid(const unsigned& topid) const;
1959 };
1960 
1961 
1962 //_______________________________________________________________________________________________________________________
1963 SbPolyhedronPolygonXSect::SbPolyhedronPolygonXSect(const std::vector<double>& x,const std::vector<double>& y, const double& dz)
1964 {
1965 
1966  Internals * internals = new Internals(this);
1967  internals->setData(&x,&y,dz);
1968  internals->triangulate();
1969 
1971  internals->addExtraVertices();
1973  internals->defineFacesTopology();
1974 
1975  delete internals;
1976 }
1977 
1978 //_______________________________________________________________________________________________________________________
1980 
1981 //_______________________________________________________________________________________________________________________
1982 void SbPolyhedronPolygonXSect::Internals::setData(const std::vector<double> * xx,const std::vector<double>* yy, const double& the_dz)
1983 {
1984  n = xx->size();
1985  ntriangles = n-2;
1986  assert (n==yy->size()&&n>2);//fixme n>2, and special code for n==3.
1987  dz = the_dz;
1988  x = xx;
1989  y = yy;
1991 }
1992 
1993 //_______________________________________________________________________________________________________________________
1995 {
1996  poly = new PolygonTriangulator(*x,*y);
1997 }
1998 
1999 //_______________________________________________________________________________________________________________________
2001  //Goal: classify edges as external/internal, get ability to find neighbours.
2002 
2003  //As a help, we create a map of edges->triangles. We define an edge
2004  //as a pair of int's, where the first should always be lower than
2005  //the second (otherwise the edges 2-4 and 4-2 would seem different):
2006  std::map<Edge,std::vector<PolygonTriangulator::Triangles::const_iterator> > edge2triangles_map;
2007  std::set<Edge> alledges;
2008  Triangles::const_iterator itt=poly->triangles()->begin();
2009  Triangles::const_iterator ittE=poly->triangles()->end();
2010  for(; itt!=ittE; ++itt)
2011  for (unsigned iedge=0;iedge<3;++iedge) {
2012  Edge eun = GetEdge(&(*itt),iedge,false/*unoriented!*/);
2013  if (edge2triangles_map.find(eun)==edge2triangles_map.end()) {
2014  alledges.insert(eun);
2015  }
2016  edge2triangles_map[eun].push_back(itt);
2017  }
2018 
2019 
2020  //set of external/internal edges:
2021  std::set<Edge>::const_iterator itedges = alledges.begin();
2022  std::set<Edge>::const_iterator itedgesE = alledges.end();
2023  for (;itedges!=itedgesE;++itedges) {
2024  if (isinternaledge((*itedges),edge2triangles_map))
2025  edges_internal.insert(*itedges);
2026  else
2027  edges_external.insert(*itedges);
2028  }
2029 
2030  //Neighbourmap:
2031  itt=poly->triangles()->begin();
2032  for(; itt!=ittE; ++itt)
2033  for (unsigned iedge=0;iedge<3;++iedge) {
2034  Edge e = GetEdge(&(*itt),iedge,true/*oriented!*/);
2035  Edge eun = GetEdge(&(*itt),iedge,false/*unoriented!*/);
2036  assert(edge2triangles_map.find(eun)!=edge2triangles_map.end());//fixme
2037  if (edges_external.find(eun)!=edges_external.end()) {
2038  neighbourmap[e]=0;
2039  neighbourmap[eun]=0;
2040  } else {
2041  assert(edge2triangles_map[eun].size()==2);//fixme
2042  const Triangle* triangle = static_cast<const Triangle*>( &(*itt) );
2043  const Triangle* triangle1 = static_cast<const Triangle*>( &(*(edge2triangles_map[eun].front())) );
2044  const Triangle* triangle2 = static_cast<const Triangle*>( &(*(edge2triangles_map[eun].back())) );
2045  neighbourmap[e]=(triangle==triangle1?triangle2:triangle1);
2046  }
2047  }
2048 }
2049 
2050 //_______________________________________________________________________________________________________________________
2052 
2062 
2063  std::set<const Triangle*> triangles_with_extra_vertex;//Eventually, all triangles should end up here.
2064 
2066  Triangles::const_iterator itt=poly->triangles()->begin();
2067  Triangles::const_iterator ittE=poly->triangles()->end();
2068  for(; itt!=ittE; ++itt) {
2069  unsigned ninternaledges(0);
2070  for (unsigned iedge=0;iedge<3;++iedge) {
2071  Edge eun = GetEdge(&(*itt),iedge,false/*unoriented!*/);
2072  if (edges_internal.find(eun)!=edges_internal.end())
2073  ++ninternaledges;
2074  }
2075  if (ninternaledges<3)
2076  continue;
2077 
2078 
2079  //Check the three edges to find a neighbour that has not already become a quadrilateral.
2080  unsigned iedge(0);
2081  for (;iedge<3;++iedge) {
2082  Edge e = GetEdge(&(*itt),iedge,true/*oriented!*/);
2083  const Triangle* trneighbour = neighbourmap[e];
2084  assert(trneighbour);//fixme
2085  if (triangles_with_extra_vertex.insert(trneighbour).second) {
2086  triangles_with_extra_vertex.insert(&(*itt));
2087  Edge eun = GetEdge(&(*itt),iedge,false/*unoriented*/);
2088  edges_with_extra_vertex.insert(eun);
2089  ++nextrainternalvertices;
2090  break;
2091  }
2092  }
2093  assert(iedge<3);//If we cannot add a vertex to a triangle with 3 internal edges, then we have a problem...
2094  }
2095 
2097  itt=poly->triangles()->begin();
2098  for(; itt!=ittE; ++itt) {
2099  Triangle* triangle = const_cast<Triangle*>( &(*itt) );
2100  if (triangles_with_extra_vertex.find(triangle)!=triangles_with_extra_vertex.end())
2101  continue;
2102  for (unsigned iedge=0;iedge<3;++iedge) {
2103  Edge eun = GetEdge(&(*itt),iedge,false/*unoriented!*/);
2104  //only look at external edges
2105  if (edges_external.find(eun)==edges_external.end())
2106  continue;
2107  assert(edges_with_extra_vertex.find(eun)==edges_with_extra_vertex.end());//fixme
2108  triangles_with_extra_vertex.insert(triangle);
2109  edges_with_extra_vertex.insert(eun);
2110  ++nextraexternalvertices;
2111  break;
2112  }
2113  }
2114  assert(ntriangles==triangles_with_extra_vertex.size());
2115 }
2116 
2117 //_______________________________________________________________________________________________________________________
2119 
2120  const unsigned nfaces = 2*ntriangles+n+nextraexternalvertices;
2121  const unsigned nvertices = 2*(n+nextraexternalvertices+nextrainternalvertices);
2122  sbpolyhedron->AllocateMemory(nvertices,nfaces);
2123 
2125 
2126  //Make edgewithextravertex_2_id map.
2127  std::set<Edge>::const_iterator itl = edges_with_extra_vertex.begin();
2128  std::set<Edge>::const_iterator itlE = edges_with_extra_vertex.end();
2129  unsigned il = 2*n;
2130  unsigned ielev = 3*n-4;
2131  for (;itl!=itlE;++itl) {
2132  edgewithextravertex_2_id[*itl]=++il;
2133  if (edges_external.find(*itl)!=edges_external.end()) {
2134  externaledgewithextravertex_2_2ndedgeid[*itl]=++ielev;
2135  }
2136  }
2137  //Now place the vertices:
2138 
2139  //First the original n vertices at +dz
2140  for (unsigned i = 0; i<n; ++i)
2141  sbpolyhedron->m_pV[i+1] = HVPoint3D(x->at(i),y->at(i),dz);
2142  //Then the original n vertices at -dz
2143  for (unsigned i = 0; i<n; ++i)
2144  sbpolyhedron->m_pV[i+1+n] = HVPoint3D(x->at(i),y->at(i),-dz);
2145  //Extra vertices at +dz.
2146  std::map<Edge,unsigned>::const_iterator itlid = edgewithextravertex_2_id.begin();
2147  std::map<Edge,unsigned>::const_iterator itlidE = edgewithextravertex_2_id.end();
2148  for (; itlid!=itlidE; ++itlid) {
2149  double xx = 0.5*(x->at(itlid->first.first-1)+x->at(itlid->first.second-1));
2150  double yy = 0.5*(y->at(itlid->first.first-1)+y->at(itlid->first.second-1));
2151  sbpolyhedron->m_pV[itlid->second] = HVPoint3D(xx,yy,dz);
2152  }
2153  //Extra vertices at -dz.
2154  itlid = edgewithextravertex_2_id.begin();
2155  for (; itlid!=itlidE; ++itlid) {
2156  double xx = 0.5*(x->at(itlid->first.first-1)+x->at(itlid->first.second-1));
2157  double yy = 0.5*(y->at(itlid->first.first-1)+y->at(itlid->first.second-1));
2158  sbpolyhedron->m_pV[itlid->second+edgewithextravertex_2_id.size()] = HVPoint3D(xx,yy,-dz);
2159  }
2160 
2161 }
2162 
2163 //_______________________________________________________________________________________________________________________
2165  const std::map<Edge,std::vector<Triangles::const_iterator> >& edge2triangles_map) const {
2166 
2167  Edge reverse_edge(edge.second,edge.first);
2168  unsigned k = 0;
2169  std::map<Edge,std::vector<Triangles::const_iterator> >::const_iterator it = edge2triangles_map.find(edge);
2170  if (it!=edge2triangles_map.end()) k += it->second.size();
2171  it = edge2triangles_map.find(reverse_edge);
2172  if (it!=edge2triangles_map.end()) k += it->second.size();
2173  return k==2;
2174 
2175 }
2176 
2177 //_______________________________________________________________________________________________________________________
2179  assert(i<=2);
2180  if (i==2)
2181  return oriented ? Edge((*tr)[2],(*tr)[0]) : Edge((*tr)[0],(*tr)[2]);
2182  else
2183  return Edge((*tr)[i],(*tr)[i+1]);
2184 }
2185 
2186 //____________________________________________________________________________________________________
2187 void SbPolyhedronPolygonXSect::Internals::setupface(const unsigned& face_id, const std::vector<unsigned>&v) const {
2188  assert(v.size()==8);
2189  sbpolyhedron->m_pF[face_id] = SbFacet(v.at(0),v.at(1),v.at(2),v.at(3),v.at(4),v.at(5),v.at(6),v.at(7));
2190 }
2191 
2192 
2193 //____________________________________________________________________________________________________
2194 void SbPolyhedronPolygonXSect::Internals::setupExternalFace(const unsigned& id_extedge,const unsigned& prev_edgeface,const unsigned& next_edgeface,
2195  const unsigned& idfirstvertex,const unsigned& idsecondvertex,
2196  const unsigned& id_triangleface,const unsigned& n, std::vector<unsigned>& faceinfo) const {
2197  faceinfo.clear();
2198  faceinfo.push_back(idfirstvertex);//point
2199  faceinfo.push_back(prev_edgeface);//neighbour face
2200  faceinfo.push_back(top2bottomvertexid(idfirstvertex));//point
2201  faceinfo.push_back(id_triangleface+n-2);//neighbour face
2202  faceinfo.push_back(top2bottomvertexid(idsecondvertex));//point
2203  faceinfo.push_back(next_edgeface);//neighbour face
2204  faceinfo.push_back(idsecondvertex);//point
2205  faceinfo.push_back(id_triangleface);//neighbour face
2206  setupface(id_extedge,faceinfo);
2207 }
2208 
2209 
2210 //____________________________________________________________________________________________________
2211 void SbPolyhedronPolygonXSect::Internals::getSurroundingValues(const std::vector<unsigned>& numbercycle, const unsigned& centralValue,
2212  unsigned& prevValue, unsigned& nextValue) const
2213 {
2214  //First, find the index of the centralValue
2215  unsigned int i;
2216  for ( i=0/*fixme, not O(n)*/;i<numbercycle.size();++i)
2217  if (numbercycle.at(i)==centralValue)
2218  break;
2219  //Sanity check that the central value was found
2220  assert(i<numbercycle.size());//fixme
2221  prevValue = ( i==0 ? numbercycle.back() : numbercycle.at(i-1) );
2222  nextValue = ( i==numbercycle.size()-1 ? numbercycle.front() : numbercycle.at(i+1) );
2223 
2224 }
2225 
2226 
2227 //_______________________________________________________________________________________________________________________
2228 unsigned SbPolyhedronPolygonXSect::Internals::top2bottomvertexid(const unsigned& topid) const {
2229  if (topid<=n)
2230  return topid+n;//original vertex
2231  return topid+nextraexternalvertices+nextrainternalvertices;//added extra vertex;
2232 }
2233 
2234 //_______________________________________________________________________________________________________________________
2235 unsigned SbPolyhedronPolygonXSect::Internals::top2bottomfaceid(const unsigned& topid) const {
2236  if (topid<=ntriangles)
2237  return topid+ntriangles;//triangle
2238  return topid;//side face
2239 }
2240 
2241 
2242 //_______________________________________________________________________________________________________________________
2244  //This is definitely the most tricky part to get right...
2245 
2246  //Make triangle_2_id map.
2247  std::map<const Triangle*,unsigned> triangle_2_id;
2248  Triangles::const_iterator itt=poly->triangles()->begin();
2249  Triangles::const_iterator ittE=poly->triangles()->end();
2250  unsigned itri = 0;
2251  for(; itt!=ittE; ++itt) {
2252  Triangle* triangle = const_cast<Triangle*>( &(*itt) );
2253  triangle_2_id[triangle]=++itri;
2254  }
2255  assert(triangle_2_id.size()==ntriangles);
2256 
2257  //Loop over the external edges, and get a list of edge face id's.
2258  std::vector<unsigned> edgeface_ids;
2259  std::map<Edge,unsigned> edge_2_firstid;
2260  std::map<Edge,unsigned> edge_2_secondid;
2261  for (unsigned i = 1; i<=n;++i) {
2262  Edge eun((i==n?1:i),(i==n?n:i+1));
2263  edgeface_ids.push_back(2*n-4+i);
2264  edge_2_firstid[eun]=edgeface_ids.back();
2265  bool hasextravertex=
2266  (edges_with_extra_vertex.find(eun)!=edges_with_extra_vertex.end());
2267  if (hasextravertex) {
2268  edgeface_ids.push_back(externaledgewithextravertex_2_2ndedgeid[eun]);
2269  edge_2_secondid[eun]=edgeface_ids.back();
2270  }
2271  }
2272 
2273  std::vector<unsigned> faceinfo_top, faceinfo_bottom(8), faceinfo_sides;
2274  itt = poly->triangles()->begin();
2275  for(; itt!=ittE; ++itt) {
2276  unsigned triangleid = triangle_2_id[&(*itt)];
2277 
2279  for (unsigned iedge=0;iedge<3;++iedge) {
2280  Edge e = GetEdge(&(*itt),iedge,true/*oriented!*/);
2281  Edge eun = GetEdge(&(*itt),iedge,false/*unoriented!*/);
2282 
2283  bool hasextravertex(edges_with_extra_vertex.find(eun)!=edges_with_extra_vertex.end());
2284  bool isinternal(edges_external.find(eun)==edges_external.end());
2285  //add for the triangle.
2286  faceinfo_top.push_back(e.first);
2287  if (isinternal) {
2288  assert(neighbourmap[e]);//fixme
2289  assert(triangle_2_id.find(neighbourmap[e])!=triangle_2_id.end());
2290  faceinfo_top.push_back(triangle_2_id[neighbourmap[e]]);
2291  if (hasextravertex) {
2292  faceinfo_top.push_back(edgewithextravertex_2_id[eun]);
2293  faceinfo_top.push_back(triangle_2_id[neighbourmap[e]]);//NB: Same neighbourface again!
2294  // assert(false&&"need to implement extravertex for internal edges");
2295  //fixme!
2296  }
2297  } else {
2298  assert(edge_2_firstid.find(eun)!=edge_2_firstid.end());
2299  unsigned firstsideid = edge_2_firstid[eun];
2300  faceinfo_top.push_back(firstsideid);
2301  if (hasextravertex) {
2302  faceinfo_top.push_back(edgewithextravertex_2_id[eun]);
2303  assert(edge_2_secondid.find(eun)!=edge_2_secondid.end());
2304  faceinfo_top.push_back(edge_2_secondid[eun]);
2305  }
2307  unsigned prev_edgeface,next_edgeface;
2308  getSurroundingValues(edgeface_ids, firstsideid,prev_edgeface, next_edgeface);
2309  unsigned firstvertexid = e.first;
2310  unsigned secondvertexid = (hasextravertex? edgewithextravertex_2_id[eun] : e.second);
2311  setupExternalFace(firstsideid,prev_edgeface, next_edgeface,firstvertexid,secondvertexid,
2312  /*firstvertexid+n,secondvertexid + (hasextravertex?externaledgewithextravertex_2_2ndedgeid.size():n),*/
2313  triangleid,n,faceinfo_sides);
2314  if (hasextravertex) {
2315  firstsideid=next_edgeface;
2316  getSurroundingValues(edgeface_ids, firstsideid,prev_edgeface, next_edgeface);
2317  firstvertexid=secondvertexid;
2318  secondvertexid=e.second;
2319  setupExternalFace(firstsideid,prev_edgeface, next_edgeface,
2320  edgewithextravertex_2_id[eun],secondvertexid,
2321  /*edgewithextravertex_2_id[eun]+externaledgewithextravertex_2_2ndedgeid.size(),secondvertexid+n,*/
2322  triangleid,n,faceinfo_sides);
2323  }
2324 
2325  }
2326  }//End edge loop
2327 
2329  assert(faceinfo_bottom.size()==8);//fixme
2330  assert(faceinfo_top.size()==8);//fixme
2331  faceinfo_bottom[0]=top2bottomvertexid(faceinfo_top[0]);
2332  faceinfo_bottom[1]=top2bottomfaceid(faceinfo_top[7]);
2333  faceinfo_bottom[2]=top2bottomvertexid(faceinfo_top[6]);
2334  faceinfo_bottom[3]=top2bottomfaceid(faceinfo_top[5]);
2335  faceinfo_bottom[4]=top2bottomvertexid(faceinfo_top[4]);
2336  faceinfo_bottom[5]=top2bottomfaceid(faceinfo_top[3]);
2337  faceinfo_bottom[6]=top2bottomvertexid(faceinfo_top[2]);
2338  faceinfo_bottom[7]=top2bottomfaceid(faceinfo_top[1]);
2339 
2340  //Actually setup triangle faces.
2341  setupface(triangleid,faceinfo_top);
2342  setupface(triangleid+ntriangles,faceinfo_bottom);
2343  faceinfo_top.clear();
2344  }//End triangle loop
2345 
2346 }
2347 
2348 // ************************************************************************************************ //
2349 // * * //
2350 // * Class SbPolyhedronArbitrary for Tessellated Solids * //
2351 // * * //
2352 // ************************************************************************************************ //
2353 
2355 {
2356  AllocateMemory(nVertices, nFacets);
2357  m_nVertexCount = 0;
2358  m_nFacetCount = 0;
2359 }
2360 
2362 {
2363 }
2364 
2365 void SbPolyhedronArbitrary::AddVertex(const double v1, const double v2, const double v3)
2366 {
2367  if(m_nVertexCount==m_nvert+1) {
2368  std::cerr <<"ERROR in SbPolyhedronArbitrary::AddVertex. Attempt to exceed maximum number of vertices: "
2369  << m_nVertexCount << std::endl;
2370  }
2371  else {
2372  m_nVertexCount++;
2373  m_pV[m_nVertexCount] = HVPoint3D(v1,v2,v3);
2374  }
2375 }
2376 
2377 void SbPolyhedronArbitrary::AddFacet(const int iv1, const int iv2, const int iv3, const int iv4)
2378 {
2379  if(m_nFacetCount==m_nface) {
2380  std::cerr <<"ERROR in SbPolyhedronArbitrary::AddFacet. Attempt to exceed maximum number of facets: "
2381  << m_nFacetCount << std::endl;
2382  }
2383  else if(iv1 < 1
2384  || iv1 > m_nvert
2385  || iv2 < 1
2386  || iv2 > m_nvert
2387  || iv3 < 1
2388  || iv3 > m_nvert
2389  || iv4 > m_nvert) {
2390  std::cerr <<"ERROR in SbPolyhedronArbitrary::AddFacet. Attempt to index vertex number which is out-of-range: ("
2391  << iv1 << "," << iv2 << "," << iv3 << "," << iv4 << ")" << std::endl;
2392  }
2393  else if(iv1 > m_nVertexCount
2394  || iv2 > m_nVertexCount
2395  || iv3 > m_nVertexCount
2396  || iv4 > m_nVertexCount) {
2397  std::cerr <<"ERROR in SbPolyhedronArbitrary::AddFacet. Vertex needs to be defined first : ("
2398  << iv1 << "," << iv2 << "," << iv3 << "," << iv4 << ")" << std::endl;
2399  }
2400  else {
2401  m_nFacetCount++;
2402  m_pF[m_nFacetCount] = SbFacet(iv1, 0, iv2, 0, iv3, 0, iv4, 0);
2403  }
2404 }
2405 
2407 {
2408  SetReferences();
2409 }
2410 
2411 SbPolyhedronGenericTrap::SbPolyhedronGenericTrap(double Dz, const std::vector<std::pair<double,double> >& Vertices)
2412 {
2413  AllocateMemory(8,6);
2414 
2415  m_pV[1] = HVPoint3D(Vertices[0].first,Vertices[0].second,-Dz);
2416  m_pV[2] = HVPoint3D(Vertices[1].first,Vertices[1].second,-Dz);
2417  m_pV[3] = HVPoint3D(Vertices[2].first,Vertices[2].second,-Dz);
2418  m_pV[4] = HVPoint3D(Vertices[3].first,Vertices[3].second,-Dz);
2423 
2424  CreatePrism();
2425 }
2426 
2428 {
2429 
2430 }
InDetVertexTruthMatchUtils::DUMMY
@ DUMMY
Definition: InDetVertexTruthMatchUtils.h:25
SbPolyhedron::m_pF
SbFacet * m_pF
Definition: SbPolyhedron.h:240
PlotCalibFromCool.il
il
Definition: PlotCalibFromCool.py:381
SbPolyhedronPolygonXSect::Internals::addExtraVertices
void addExtraVertices()
Definition: SbPolyhedron.cxx:2051
beamspotman.r
def r
Definition: beamspotman.py:676
python.SystemOfUnits.second
int second
Definition: SystemOfUnits.py:120
SbPolyhedron
Definition: SbPolyhedron.h:231
SbPolyhedronPolygonXSect::~SbPolyhedronPolygonXSect
virtual ~SbPolyhedronPolygonXSect()
Definition: SbPolyhedron.cxx:1979
python.CaloRecoConfig.f
f
Definition: CaloRecoConfig.py:127
SbPolyhedron::FindNodeNormal
HVNormal3D FindNodeNormal(int iFace, int iNode) const
Definition: SbPolyhedron.cxx:217
SbPolyhedron::GetVertex
HVPoint3D GetVertex(int index) const
Definition: SbPolyhedron.cxx:929
SbPolyhedronArbitrary::~SbPolyhedronArbitrary
virtual ~SbPolyhedronArbitrary()
Definition: SbPolyhedron.cxx:2361
SbPolyhedron::GetVertexFast
const HVPoint3D & GetVertexFast(int index) const
Definition: SbPolyhedron.cxx:950
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
SbPolyhedronPolygonXSect::Internals
Definition: SbPolyhedron.cxx:1906
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
SbPolyhedron::s_numberOfRotationSteps
static int s_numberOfRotationSteps
Definition: SbPolyhedron.h:235
NSWL1::nVertices
int nVertices(const Polygon &p)
Definition: GeoUtils.cxx:35
SbPolyhedron::GetNextVertex
bool GetNextVertex(HVPoint3D &vertex, int &edgeFlag) const
Definition: SbPolyhedron.cxx:957
SbPolyhedronPolygonXSect::Internals::externaledgewithextravertex_2_2ndedgeid
std::map< Edge, unsigned > externaledgewithextravertex_2_2ndedgeid
Definition: SbPolyhedron.cxx:1943
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
PlotCalibFromCool.yy
yy
Definition: PlotCalibFromCool.py:714
SbPolyhedron::add
SbPolyhedron add(const SbPolyhedron &p) const
Definition: SbPolyhedron.cxx:1842
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
operator<<
std::ostream & operator<<(std::ostream &ostr, const SbFacet &facet)
Definition: SbPolyhedron.cxx:118
SbPolyhedronPolygonXSect::Internals::GetEdge
Edge GetEdge(const Triangle *tr, const unsigned &i, const bool &oriented=false) const
Definition: SbPolyhedron.cxx:2178
SbPolyhedronArbitrary::SbPolyhedronArbitrary
SbPolyhedronArbitrary(const int nVertices, const int nFacets)
Definition: SbPolyhedron.cxx:2354
HVPoint3D::HVPoint3D
HVPoint3D()
Definition: SbPolyhedron.cxx:31
beamspotman.cur
def cur
Definition: beamspotman.py:671
index
Definition: index.py:1
make_hlt_rep.rep
rep
Definition: make_hlt_rep.py:32
SbPolyhedronBox::SbPolyhedronBox
SbPolyhedronBox(double Dx, double Dy, double Dz)
Definition: SbPolyhedron.cxx:1382
OP_INTERSECTION
#define OP_INTERSECTION
Definition: BooleanProcessor.h:90
SbPolyhedronTube::SbPolyhedronTube
SbPolyhedronTube(double Rmin, double Rmax, double Dz)
Definition: SbPolyhedron.cxx:1544
SbPolyhedron::GetNextVertexIndex
bool GetNextVertexIndex(int &index, int &edgeFlag) const
Definition: SbPolyhedron.cxx:902
SbPolyhedron::GetNextEdgeIndeces
bool GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag, int &iface1, int &iface2) const
Definition: SbPolyhedron.cxx:1007
skel.it
it
Definition: skel.GENtoEVGEN.py:423
SbPolyhedronCone::SbPolyhedronCone
SbPolyhedronCone(double Rmn1, double Rmx1, double Rmn2, double Rmx2, double Dz)
Definition: SbPolyhedron.cxx:1530
SbPolyhedron::GetFacet
void GetFacet(int iFace, int &n, int *iNodes, int *edgeFlags=0, int *iFaces=0) const
Definition: SbPolyhedron.cxx:1117
M_PI
#define M_PI
Definition: ActiveFraction.h:11
ffabs
float ffabs(float a)
Definition: SbPolyhedron.cxx:25
deg
#define deg
Definition: SbPolyhedron.cxx:17
PolygonTriangulator::Triangles
std::list< Triangle > Triangles
Definition: PolygonTriangulator.h:36
SbPolyhedronPolygonXSect::Internals::sbpolyhedron
SbPolyhedronPolygonXSect * sbpolyhedron
Definition: SbPolyhedron.cxx:1931
SbPolyhedronPolygonXSect::Internals::neighbourmap
std::map< Edge, const Triangle * > neighbourmap
Definition: SbPolyhedron.cxx:1938
perMillion
#define perMillion
Definition: SbPolyhedron.cxx:16
SbPolyhedronTubs::SbPolyhedronTubs
SbPolyhedronTubs(double Rmin, double Rmax, double Dz, double Phi1, double Dphi)
Definition: SbPolyhedron.cxx:1537
SbPolyhedronTrap::SbPolyhedronTrap
SbPolyhedronTrap(double Dz, double Theta, double Phi, double Dy1, double Dx1, double Dx2, double Alp1, double Dy2, double Dx3, double Dx4, double Alp2)
Definition: SbPolyhedron.cxx:1387
Phi
@ Phi
Definition: RPCdef.h:8
SbPolyhedronSphere::SbPolyhedronSphere
SbPolyhedronSphere(double rmin, double rmax, double phi, double dphi, double the, double dthe)
Definition: SbPolyhedron.cxx:1649
SbPolyhedronTorus::SbPolyhedronTorus
SbPolyhedronTorus(double rmin, double rmax, double rtor, double phi, double dphi)
Definition: SbPolyhedron.cxx:1749
DEFAULT_NUMBER_OF_STEPS
#define DEFAULT_NUMBER_OF_STEPS
Definition: SbPolyhedron.h:214
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
SbPolyhedronPolygonXSect::Internals::dz
double dz
Definition: SbPolyhedron.cxx:1926
const
bool const RAWDATA *ch2 const
Definition: LArRodBlockPhysicsV0.cxx:562
SbFacet::f
int f
Definition: SbPolyhedron.h:222
SbPolyhedron::operator=
virtual SbPolyhedron & operator=(const SbPolyhedron &from)
Definition: SbPolyhedron.cxx:161
SbPolyhedron::SetReferences
void SetReferences()
Definition: SbPolyhedron.cxx:682
x
#define x
iabs
int iabs(int a)
Definition: SbPolyhedron.cxx:22
SbPolyhedronTubs::~SbPolyhedronTubs
virtual ~SbPolyhedronTubs()
Definition: SbPolyhedron.cxx:1542
SbPolyhedron::GetNormal
HVNormal3D GetNormal(int iFace) const
Definition: SbPolyhedron.cxx:1203
SbPolyhedronPolygonXSect::Internals::allocateMemoryAndDefineVertexCoordinates
void allocateMemoryAndDefineVertexCoordinates()
Definition: SbPolyhedron.cxx:2118
SbPolyhedronCons::~SbPolyhedronCons
virtual ~SbPolyhedronCons()
Definition: SbPolyhedron.cxx:1528
SbPolyhedronPolygonXSect::Internals
friend class Internals
Definition: SbPolyhedron.h:547
SbPolyhedron::GetNextEdge
bool GetNextEdge(HVPoint3D &p1, HVPoint3D &p2, int &edgeFlag) const
Definition: SbPolyhedron.cxx:1076
SbPolyhedron::CreatePrism
void CreatePrism()
Definition: SbPolyhedron.cxx:286
SbPolyhedron::GetNextNormal
bool GetNextNormal(HVNormal3D &normal) const
Definition: SbPolyhedron.cxx:1255
SbPolyhedronPolygonXSect::Internals::Internals
Internals(SbPolyhedronPolygonXSect *sbp)
Definition: SbPolyhedron.cxx:1914
SbPolyhedronTrd2::SbPolyhedronTrd2
SbPolyhedronTrd2(double Dx1, double Dx2, double Dy1, double Dy2, double Dz)
Definition: SbPolyhedron.cxx:1342
SbPolyhedronPgon::SbPolyhedronPgon
SbPolyhedronPgon(double phi, double dphi, int npdv, int nz, const double *z, const double *rmin, const double *rmax)
Definition: SbPolyhedron.cxx:1550
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
python.changerun.kk
list kk
Definition: changerun.py:41
SbPolyhedronPolygonXSect::Internals::top2bottomfaceid
unsigned top2bottomfaceid(const unsigned &topid) const
Definition: SbPolyhedron.cxx:2235
SbPolyhedron::GetVolume
double GetVolume() const
Definition: SbPolyhedron.cxx:1314
HVPoint3D::operator=
HVPoint3D & operator=(const HVPoint3D &v)
Definition: SbPolyhedron.cxx:35
SbPolyhedron::Transform
SbPolyhedron & Transform(const HEPVis::SbRotation &rot, const SbVec3d &trans)
Definition: SbPolyhedron.cxx:878
OP_UNION
#define OP_UNION
Definition: BooleanProcessor.h:89
SbPolyhedronPolygonXSect::Internals::Triangle
PolygonTriangulator::Triangle Triangle
Definition: SbPolyhedron.cxx:1910
SbPolyhedronTrap
Definition: SbPolyhedron.h:425
SbPolyhedronPolygonXSect::Internals::setData
void setData(const std::vector< double > *xx, const std::vector< double > *yy, const double &dz)
Definition: SbPolyhedron.cxx:1982
SbPolyhedronPcon::SbPolyhedronPcon
SbPolyhedronPcon(double phi, double dphi, int nz, const double *z, const double *rmin, const double *rmax)
Definition: SbPolyhedron.cxx:1641
SbPolyhedronTube::~SbPolyhedronTube
virtual ~SbPolyhedronTube()
Definition: SbPolyhedron.cxx:1548
dqt_zlumi_pandas.err
err
Definition: dqt_zlumi_pandas.py:193
fillPileUpNoiseLumi.next
next
Definition: fillPileUpNoiseLumi.py:52
lumiFormat.i
int i
Definition: lumiFormat.py:92
PolygonTriangulator.h
z
#define z
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
beamspotman.n
n
Definition: beamspotman.py:731
SbPolyhedronPolygonXSect::Internals::nextraexternalvertices
unsigned nextraexternalvertices
Definition: SbPolyhedron.cxx:1940
Dz
double Dz
Definition: LArDetectorConstructionTBEC.cxx:57
xAOD::rotation
rotation
Definition: TrackSurface_v1.cxx:15
SbPolyhedronPara::~SbPolyhedronPara
virtual ~SbPolyhedronPara()
Definition: SbPolyhedron.cxx:1447
SbPolyhedron::AllocateMemory
void AllocateMemory(int Nvert, int Nface)
Definition: SbPolyhedron.cxx:267
SbPolyhedron::GetNumberOfRotationSteps
static int GetNumberOfRotationSteps()
Definition: SbPolyhedron.h:375
SbPolyhedronSphere::~SbPolyhedronSphere
virtual ~SbPolyhedronSphere()
Definition: SbPolyhedron.cxx:1747
SbPolyhedronTrd2
Definition: SbPolyhedron.h:396
Trk::RIGHT
@ RIGHT
the drift radius is positive (see Trk::AtaStraightLine)
Definition: DriftCircleSide.h:22
SbPolyhedron::GetUnitNormal
HVNormal3D GetUnitNormal(int iFace) const
Definition: SbPolyhedron.cxx:1228
SbPolyhedronTrap::~SbPolyhedronTrap
virtual ~SbPolyhedronTrap()
Definition: SbPolyhedron.cxx:1440
SbPolyhedronPolygonXSect::Internals::setupExternalFace
void setupExternalFace(const unsigned &id_extedge, const unsigned &prev_edgeface, const unsigned &next_edgeface, const unsigned &idfirstvertex, const unsigned &idsecondvertex, const unsigned &id_triangleface, const unsigned &n, std::vector< unsigned > &faceinfo) const
Definition: SbPolyhedron.cxx:2194
SbPolyhedronPgon::~SbPolyhedronPgon
virtual ~SbPolyhedronPgon()
Definition: SbPolyhedron.cxx:1639
SbPolyhedronPolygonXSect::Internals::Edge
std::pair< int, int > Edge
Definition: SbPolyhedron.cxx:1909
SbPolyhedron::SetSideFacets
void SetSideFacets(int ii[4], int vv[4], int *kk, double *r, double dphi, int ns, int &kface)
Definition: SbPolyhedron.cxx:368
SbFacet::m_edge
struct SbFacet::@63 m_edge[4]
SbPolyhedronPolygonXSect::Internals::nextrainternalvertices
unsigned nextrainternalvertices
Definition: SbPolyhedron.cxx:1941
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
SbPolyhedronPolygonXSect::Internals::setupface
void setupface(const unsigned &face_id, const std::vector< unsigned > &v) const
Definition: SbPolyhedron.cxx:2187
SbPolyhedronPgon
Definition: SbPolyhedron.h:488
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
SbPolyhedronArbitrary::Finalize
void Finalize()
Definition: SbPolyhedron.cxx:2406
dot.dot
def dot(G, fn, nodesToHighlight=[])
Definition: dot.py:5
SbPolyhedron::GetNextFacet
bool GetNextFacet(int &n, HVPoint3D *nodes, int *edgeFlags=0, HVNormal3D *normals=0) const
Definition: SbPolyhedron.cxx:1173
SbPolyhedronPolygonXSect::Internals::edges_with_extra_vertex
std::set< Edge > edges_with_extra_vertex
Definition: SbPolyhedron.cxx:1939
SbPolyhedron::m_pV
HVPoint3D * m_pV
Definition: SbPolyhedron.h:239
SbFacet::v
int v
Definition: SbPolyhedron.h:222
SbPolyhedronPcon::~SbPolyhedronPcon
virtual ~SbPolyhedronPcon()
Definition: SbPolyhedron.cxx:1647
SbPolyhedronPolygonXSect::Internals::edges_internal
std::set< Edge > edges_internal
Definition: SbPolyhedron.cxx:1936
SbPolyhedronCone::~SbPolyhedronCone
virtual ~SbPolyhedronCone()
Definition: SbPolyhedron.cxx:1535
SbPolyhedron::RotateAroundZ
void RotateAroundZ(int nstep, double phi, double dphi, int np1, int np2, const double *z, double *r, int nodeVis, int edgeVis)
Definition: SbPolyhedron.cxx:438
OP_SUBTRACTION
#define OP_SUBTRACTION
Definition: BooleanProcessor.h:91
PolygonTriangulator
Definition: PolygonTriangulator.h:32
SbPolyhedronArbitrary::AddFacet
void AddFacet(const int iv1, const int iv2, const int iv3, const int iv4=0)
Definition: SbPolyhedron.cxx:2377
SbPolyhedronCons::SbPolyhedronCons
SbPolyhedronCons(double Rmn1, double Rmx1, double Rmn2, double Rmx2, double Dz, double Phi1, double Dphi)
Definition: SbPolyhedron.cxx:1449
SbPolyhedron::FindNeighbour
int FindNeighbour(int iFace, int iNode, int iOrder) const
Definition: SbPolyhedron.cxx:189
SbPolyhedronPolygonXSect::Internals::edges_external
std::set< Edge > edges_external
Definition: SbPolyhedron.cxx:1937
SbPolyhedronPolygonXSect::Internals::ntriangles
unsigned ntriangles
Definition: SbPolyhedron.cxx:1930
SbPolyhedronPolygonXSect::Internals::edgewithextravertex_2_id
std::map< Edge, unsigned > edgewithextravertex_2_id
Definition: SbPolyhedron.cxx:1942
SbPolyhedronPolygonXSect::Internals::top2bottomvertexid
unsigned top2bottomvertexid(const unsigned &topid) const
Definition: SbPolyhedron.cxx:2228
HVNormal3D
HVPoint3D HVNormal3D
Definition: SbPolyhedron.h:199
ReadCellNoiseFromCoolCompare.v2
v2
Definition: ReadCellNoiseFromCoolCompare.py:364
SbPolyhedronPolygonXSect::Internals::isinternaledge
bool isinternaledge(const Edge &edge, const std::map< Edge, std::vector< Triangles::const_iterator > > &edge2triangles_map) const
Definition: SbPolyhedron.cxx:2164
python.PyAthena.v
v
Definition: PyAthena.py:157
SbPolyhedronTorus::~SbPolyhedronTorus
virtual ~SbPolyhedronTorus()
Definition: SbPolyhedron.cxx:1827
SbPolyhedron::GetSurfaceArea
double GetSurfaceArea() const
Definition: SbPolyhedron.cxx:1292
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
SbPolyhedron::SetNumberOfRotationSteps
static void SetNumberOfRotationSteps(int n)
Definition: SbPolyhedron.cxx:245
DeMoScan.index
string index
Definition: DeMoScan.py:362
SbPolyhedronArbitrary::AddVertex
void AddVertex(const double v1, const double v2, const double v3)
Definition: SbPolyhedron.cxx:2365
SbPolyhedron::m_nface
int m_nface
Definition: SbPolyhedron.h:238
SbPolyhedronPolygonXSect::Internals::initEdgeClassificationsAndNeighbours
void initEdgeClassificationsAndNeighbours()
Definition: SbPolyhedron.cxx:2000
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
a
TList * a
Definition: liststreamerinfos.cxx:10
CalibCoolCompareRT.nm
nm
Definition: CalibCoolCompareRT.py:110
y
#define y
SbPolyhedronPolygonXSect::Internals::x
const std::vector< double > * x
Definition: SbPolyhedron.cxx:1927
SbPolyhedron::InvertFacets
void InvertFacets()
Definition: SbPolyhedron.cxx:819
DeMoScan.first
bool first
Definition: DeMoScan.py:534
SbPolyhedronPolygonXSect::Internals::n
unsigned n
Definition: SbPolyhedron.cxx:1929
SbPolyhedronPolygonXSect::Internals::getSurroundingValues
void getSurroundingValues(const std::vector< unsigned > &numbercycle, const unsigned &centralValue, unsigned &prevValue, unsigned &nextValue) const
Definition: SbPolyhedron.cxx:2211
SbPolyhedronTrd1::SbPolyhedronTrd1
SbPolyhedronTrd1(double Dx1, double Dx2, double Dy, double Dz)
Definition: SbPolyhedron.cxx:1376
SbPolyhedronPolygonXSect::Internals::~Internals
~Internals()
Definition: SbPolyhedron.cxx:1917
SbPolyhedronPolygonXSect::Internals::Triangles
PolygonTriangulator::Triangles Triangles
Definition: SbPolyhedron.cxx:1911
SbPolyhedronPolygonXSect::Internals::poly
PolygonTriangulator * poly
Definition: SbPolyhedron.cxx:1934
Rmin
double Rmin
Definition: LArDetectorConstructionTBEC.cxx:54
SbPolyhedronPolygonXSect::SbPolyhedronPolygonXSect
SbPolyhedronPolygonXSect(const std::vector< double > &x, const std::vector< double > &y, const double &dz)
Definition: SbPolyhedron.cxx:1963
SbPolyhedron::RotateEdge
void RotateEdge(int k1, int k2, double r1, double r2, int v1, int v2, int vEdge, bool ifWholeCircle, int ns, int &kface)
Definition: SbPolyhedron.cxx:306
python.SystemOfUnits.ns
int ns
Definition: SystemOfUnits.py:130
SbPolyhedronPolygonXSect::Internals::y
const std::vector< double > * y
Definition: SbPolyhedron.cxx:1928
SbPolyhedron::SbPolyhedron
SbPolyhedron(int Nvert=0, int Nface=0)
Definition: SbPolyhedron.h:278
Trk::LEFT
@ LEFT
the drift radius is negative (see Trk::AtaStraightLine)
Definition: DriftCircleSide.h:20
SbFacet
Definition: SbPolyhedron.h:217
rr
const boost::regex rr(r_r)
SbPolyhedronPara::SbPolyhedronPara
SbPolyhedronPara(double Dx, double Dy, double Dz, double Alpha, double Theta, double Phi)
Definition: SbPolyhedron.cxx:1442
BooleanProcessor.h
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
SbPolyhedronGenericTrap::~SbPolyhedronGenericTrap
virtual ~SbPolyhedronGenericTrap()
Definition: SbPolyhedron.cxx:2427
SbPolyhedron::subtract
SbPolyhedron subtract(const SbPolyhedron &p) const
Definition: SbPolyhedron.cxx:1878
SbPolyhedronCons
Definition: SbPolyhedron.h:448
PlotCalibFromCool.vv
vv
Definition: PlotCalibFromCool.py:716
length
double length(const pvec &v)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:26
NSWL1::Vertices
std::vector< Vertex > Vertices
Definition: GeoUtils.h:19
SbPolyhedron::GetNextUnitNormal
bool GetNextUnitNormal(HVNormal3D &normal) const
Definition: SbPolyhedron.cxx:1276
SbPolyhedron.h
SbPolyhedronBox::~SbPolyhedronBox
virtual ~SbPolyhedronBox()
Definition: SbPolyhedron.cxx:1385
PolygonTriangulator::Triangle
std::vector< unsigned > Triangle
Definition: PolygonTriangulator.h:35
SbPolyhedronTrd1::~SbPolyhedronTrd1
virtual ~SbPolyhedronTrd1()
Definition: SbPolyhedron.cxx:1380
SbPolyhedronPolygonXSect::Internals::defineFacesTopology
void defineFacesTopology()
Definition: SbPolyhedron.cxx:2243
SbPolyhedron::intersect
SbPolyhedron intersect(const SbPolyhedron &p) const
Definition: SbPolyhedron.cxx:1860
operator+
HVPoint3D operator+(const HVPoint3D &v1, const HVPoint3D &v2)
Definition: SbPolyhedron.cxx:43
HVPoint3D
Definition: SbPolyhedron.h:188
SbPolyhedron::m_nvert
int m_nvert
Definition: SbPolyhedron.h:238
HEPVis::SbRotation
Definition: SbRotation.h:35
fitman.k
k
Definition: fitman.py:528
SbPolyhedronPolygonXSect::Internals::triangulate
void triangulate()
Definition: SbPolyhedron.cxx:1994
SbPolyhedronTrd2::~SbPolyhedronTrd2
virtual ~SbPolyhedronTrd2()
Definition: SbPolyhedron.cxx:1374
SbPolyhedronPolygonXSect
Definition: SbPolyhedron.h:534
SbPolyhedronGenericTrap::SbPolyhedronGenericTrap
SbPolyhedronGenericTrap(double Dz, const std::vector< std::pair< double, double > > &Vertices)
Definition: SbPolyhedron.cxx:2411