Geant4: HepPolyhedron.cc Source File

00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id: HepPolyhedron.cc 69794 2013-05-15 13:17:48Z gcosmo $
00028 //
00029 // 
00030 //
00031 // G4 Polyhedron library
00032 //
00033 // History:
00034 // 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
00035 //
00036 // 30.09.96 E.Chernyaev
00037 // - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
00038 // - added GetNextUnitNormal, GetNextEdgeIndeces, GetNextEdge
00039 //
00040 // 15.12.96 E.Chernyaev
00041 // - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
00042 // - rewritten G4PolyhedronCons;
00043 // - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
00044 //
00045 // 01.06.97 E.Chernyaev
00046 // - modified RotateAroundZ, added SetSideFacets
00047 //
00048 // 19.03.00 E.Chernyaev
00049 // - implemented boolean operations (add, subtract, intersect) on polyhedra;
00050 //
00051 // 25.05.01 E.Chernyaev
00052 // - added GetSurfaceArea() and GetVolume();
00053 //
00054 // 05.11.02 E.Chernyaev
00055 // - added createTwistedTrap() and createPolyhedron();
00056 //
00057 // 20.06.05 G.Cosmo
00058 // - added HepPolyhedronEllipsoid;
00059 //
00060 // 18.07.07 T.Nikitin
00061 // - added HepParaboloid;
00062   
00063 #include "HepPolyhedron.h"
00064 #include "G4PhysicalConstants.hh"
00065 #include "G4Vector3D.hh"
00066 
00067 #include <cstdlib>  // Required on some compilers for std::abs(int) ...
00068 #include <cmath>
00069 
00070 using CLHEP::perMillion;
00071 using CLHEP::deg;
00072 using CLHEP::pi;
00073 using CLHEP::twopi;
00074 using CLHEP::nm;
00075 const G4double spatialTolerance = 0.01*nm;
00076 
00077 /***********************************************************************
00078  *                                                                     *
00079  * Name: HepPolyhedron operator <<                   Date:    09.05.96 *
00080  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
00081  *                                                                     *
00082  * Function: Print contents of G4 polyhedron                           *
00083  *                                                                     *
00084  ***********************************************************************/
00085 std::ostream & operator<<(std::ostream & ostr, const G4Facet & facet) {
00086   for (G4int k=0; k<4; k++) {
00087     ostr << " " << facet.edge[k].v << "/" << facet.edge[k].f;
00088   }
00089   return ostr;
00090 }
00091 
00092 std::ostream & operator<<(std::ostream & ostr, const HepPolyhedron & ph) {
00093   ostr << std::endl;
00094   ostr << "Nverteces=" << ph.nvert << ", Nfacets=" << ph.nface << std::endl;
00095   G4int i;
00096   for (i=1; i<=ph.nvert; i++) {
00097      ostr << "xyz(" << i << ")="
00098           << ph.pV[i].x() << ' ' << ph.pV[i].y() << ' ' << ph.pV[i].z()
00099           << std::endl;
00100   }
00101   for (i=1; i<=ph.nface; i++) {
00102     ostr << "face(" << i << ")=" << ph.pF[i] << std::endl;
00103   }
00104   return ostr;
00105 }
00106 
00107 HepPolyhedron::HepPolyhedron(const HepPolyhedron &from)
00108 /***********************************************************************
00109  *                                                                     *
00110  * Name: HepPolyhedron copy constructor             Date:    23.07.96  *
00111  * Author: E.Chernyaev (IHEP/Protvino)              Revised:           *
00112  *                                                                     *
00113  ***********************************************************************/
00114 : nvert(0), nface(0), pV(0), pF(0)
00115 {
00116   AllocateMemory(from.nvert, from.nface);
00117   for (G4int i=1; i<=nvert; i++) pV[i] = from.pV[i];
00118   for (G4int k=1; k<=nface; k++) pF[k] = from.pF[k];
00119 }
00120 
00121 HepPolyhedron & HepPolyhedron::operator=(const HepPolyhedron &from)
00122 /***********************************************************************
00123  *                                                                     *
00124  * Name: HepPolyhedron operator =                   Date:    23.07.96  *
00125  * Author: E.Chernyaev (IHEP/Protvino)              Revised:           *
00126  *                                                                     *
00127  * Function: Copy contents of one polyhedron to another                *
00128  *                                                                     *
00129  ***********************************************************************/
00130 {
00131   if (this != &from) {
00132     AllocateMemory(from.nvert, from.nface);
00133     for (G4int i=1; i<=nvert; i++) pV[i] = from.pV[i];
00134     for (G4int k=1; k<=nface; k++) pF[k] = from.pF[k];
00135   }
00136   return *this;
00137 }
00138 
00139 G4int
00140 HepPolyhedron::FindNeighbour(G4int iFace, G4int iNode, G4int iOrder) const
00141 /***********************************************************************
00142  *                                                                     *
00143  * Name: HepPolyhedron::FindNeighbour                Date:    22.11.99 *
00144  * Author: E.Chernyaev                               Revised:          *
00145  *                                                                     *
00146  * Function: Find neighbouring face                                    *
00147  *                                                                     *
00148  ***********************************************************************/
00149 {
00150   G4int i;
00151   for (i=0; i<4; i++) {
00152     if (iNode == std::abs(pF[iFace].edge[i].v)) break;
00153   }
00154   if (i == 4) {
00155     std::cerr
00156       << "HepPolyhedron::FindNeighbour: face " << iFace
00157       << " has no node " << iNode
00158       << std::endl; 
00159     return 0;
00160   }
00161   if (iOrder < 0) {
00162     if ( --i < 0) i = 3;
00163     if (pF[iFace].edge[i].v == 0) i = 2;
00164   }
00165   return (pF[iFace].edge[i].v > 0) ? 0 : pF[iFace].edge[i].f;
00166 }
00167 
00168 G4Normal3D HepPolyhedron::FindNodeNormal(G4int iFace, G4int iNode) const
00169 /***********************************************************************
00170  *                                                                     *
00171  * Name: HepPolyhedron::FindNodeNormal               Date:    22.11.99 *
00172  * Author: E.Chernyaev                               Revised:          *
00173  *                                                                     *
00174  * Function: Find normal at given node                                 *
00175  *                                                                     *
00176  ***********************************************************************/
00177 {
00178   G4Normal3D   normal = GetUnitNormal(iFace);
00179   G4int          k = iFace, iOrder = 1, n = 1;
00180 
00181   for(;;) {
00182     k = FindNeighbour(k, iNode, iOrder);
00183     if (k == iFace) break; 
00184     if (k > 0) {
00185       n++;
00186       normal += GetUnitNormal(k);
00187     }else{
00188       if (iOrder < 0) break;
00189       k = iFace;
00190       iOrder = -iOrder;
00191     }
00192   }
00193   return normal.unit();
00194 }
00195 
00196 G4int HepPolyhedron::GetNumberOfRotationSteps()
00197 /***********************************************************************
00198  *                                                                     *
00199  * Name: HepPolyhedron::GetNumberOfRotationSteps     Date:    24.06.97 *
00200  * Author: J.Allison (Manchester University)         Revised:          *
00201  *                                                                     *
00202  * Function: Get number of steps for whole circle                      *
00203  *                                                                     *
00204  ***********************************************************************/
00205 {
00206   return fNumberOfRotationSteps;
00207 }
00208 
00209 void HepPolyhedron::SetNumberOfRotationSteps(G4int n)
00210 /***********************************************************************
00211  *                                                                     *
00212  * Name: HepPolyhedron::SetNumberOfRotationSteps     Date:    24.06.97 *
00213  * Author: J.Allison (Manchester University)         Revised:          *
00214  *                                                                     *
00215  * Function: Set number of steps for whole circle                      *
00216  *                                                                     *
00217  ***********************************************************************/
00218 {
00219   const G4int nMin = 3;
00220   if (n < nMin) {
00221     std::cerr 
00222       << "HepPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
00223       << "number of steps per circle < " << nMin << "; forced to " << nMin
00224       << std::endl;
00225     fNumberOfRotationSteps = nMin;
00226   }else{
00227     fNumberOfRotationSteps = n;
00228   }    
00229 }
00230 
00231 void HepPolyhedron::ResetNumberOfRotationSteps()
00232 /***********************************************************************
00233  *                                                                     *
00234  * Name: HepPolyhedron::GetNumberOfRotationSteps     Date:    24.06.97 *
00235  * Author: J.Allison (Manchester University)         Revised:          *
00236  *                                                                     *
00237  * Function: Reset number of steps for whole circle to default value   *
00238  *                                                                     *
00239  ***********************************************************************/
00240 {
00241   fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
00242 }
00243 
00244 void HepPolyhedron::AllocateMemory(G4int Nvert, G4int Nface)
00245 /***********************************************************************
00246  *                                                                     *
00247  * Name: HepPolyhedron::AllocateMemory               Date:    19.06.96 *
00248  * Author: E.Chernyaev (IHEP/Protvino)               Revised: 05.11.02 *
00249  *                                                                     *
00250  * Function: Allocate memory for GEANT4 polyhedron                     *
00251  *                                                                     *
00252  * Input: Nvert - number of nodes                                      *
00253  *        Nface - number of faces                                      *
00254  *                                                                     *
00255  ***********************************************************************/
00256 {
00257   if (nvert == Nvert && nface == Nface) return;
00258   if (pV != 0) delete [] pV;
00259   if (pF != 0) delete [] pF;
00260   if (Nvert > 0 && Nface > 0) {
00261     nvert = Nvert;
00262     nface = Nface;
00263     pV    = new G4Point3D[nvert+1];
00264     pF    = new G4Facet[nface+1];
00265   }else{
00266     nvert = 0; nface = 0; pV = 0; pF = 0;
00267   }
00268 }
00269 
00270 void HepPolyhedron::CreatePrism()
00271 /***********************************************************************
00272  *                                                                     *
00273  * Name: HepPolyhedron::CreatePrism                  Date:    15.07.96 *
00274  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
00275  *                                                                     *
00276  * Function: Set facets for a prism                                    *
00277  *                                                                     *
00278  ***********************************************************************/
00279 {
00280   enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
00281 
00282   pF[1] = G4Facet(1,LEFT,  4,BACK,  3,RIGHT,  2,FRONT);
00283   pF[2] = G4Facet(5,TOP,   8,BACK,  4,BOTTOM, 1,FRONT);
00284   pF[3] = G4Facet(8,TOP,   7,RIGHT, 3,BOTTOM, 4,LEFT);
00285   pF[4] = G4Facet(7,TOP,   6,FRONT, 2,BOTTOM, 3,BACK);
00286   pF[5] = G4Facet(6,TOP,   5,LEFT,  1,BOTTOM, 2,RIGHT);
00287   pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK,   8,LEFT);
00288 }
00289 
00290 void HepPolyhedron::RotateEdge(G4int k1, G4int k2, G4double r1, G4double r2,
00291                               G4int v1, G4int v2, G4int vEdge,
00292                               G4bool ifWholeCircle, G4int nds, G4int &kface)
00293 /***********************************************************************
00294  *                                                                     *
00295  * Name: HepPolyhedron::RotateEdge                   Date:    05.12.96 *
00296  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
00297  *                                                                     *
00298  * Function: Create set of facets by rotation of an edge around Z-axis *
00299  *                                                                     *
00300  * Input: k1, k2 - end vertices of the edge                            *
00301  *        r1, r2 - radiuses of the end vertices                        *
00302  *        v1, v2 - visibility of edges produced by rotation of the end *
00303  *                 vertices                                            *
00304  *        vEdge  - visibility of the edge                              *
00305  *        ifWholeCircle - is true in case of whole circle rotation     *
00306  *        nds    - number of discrete steps                            *
00307  *        r[]    - r-coordinates                                       *
00308  *        kface  - current free cell in the pF array                   *
00309  *                                                                     *
00310  ***********************************************************************/
00311 {
00312   if (r1 == 0. && r2 == 0) return;
00313 
00314   G4int i;
00315   G4int i1  = k1;
00316   G4int i2  = k2;
00317   G4int ii1 = ifWholeCircle ? i1 : i1+nds;
00318   G4int ii2 = ifWholeCircle ? i2 : i2+nds;
00319   G4int vv  = ifWholeCircle ? vEdge : 1;
00320 
00321   if (nds == 1) {
00322     if (r1 == 0.) {
00323       pF[kface++]   = G4Facet(i1,0,    v2*i2,0, (i2+1),0);
00324     }else if (r2 == 0.) {
00325       pF[kface++]   = G4Facet(i1,0,    i2,0,    v1*(i1+1),0);
00326     }else{
00327       pF[kface++]   = G4Facet(i1,0,    v2*i2,0, (i2+1),0, v1*(i1+1),0);
00328     }
00329   }else{
00330     if (r1 == 0.) {
00331       pF[kface++]   = G4Facet(vv*i1,0,    v2*i2,0, vEdge*(i2+1),0);
00332       for (i2++,i=1; i<nds-1; i2++,i++) {
00333         pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
00334       }
00335       pF[kface++]   = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
00336     }else if (r2 == 0.) {
00337       pF[kface++]   = G4Facet(vv*i1,0,    vEdge*i2,0, v1*(i1+1),0);
00338       for (i1++,i=1; i<nds-1; i1++,i++) {
00339         pF[kface++] = G4Facet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
00340       }
00341       pF[kface++]   = G4Facet(vEdge*i1,0, vv*i2,0,    v1*ii1,0);
00342     }else{
00343       pF[kface++]   = G4Facet(vv*i1,0,    v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
00344       for (i1++,i2++,i=1; i<nds-1; i1++,i2++,i++) {
00345         pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
00346       }  
00347       pF[kface++]   = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0,      v1*ii1,0);
00348     }
00349   }
00350 }
00351 
00352 void HepPolyhedron::SetSideFacets(G4int ii[4], G4int vv[4],
00353                                  G4int *kk, G4double *r,
00354                                  G4double dphi, G4int nds, G4int &kface)
00355 /***********************************************************************
00356  *                                                                     *
00357  * Name: HepPolyhedron::SetSideFacets                Date:    20.05.97 *
00358  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
00359  *                                                                     *
00360  * Function: Set side facets for the case of incomplete rotation       *
00361  *                                                                     *
00362  * Input: ii[4] - indeces of original verteces                         *
00363  *        vv[4] - visibility of edges                                  *
00364  *        kk[]  - indeces of nodes                                     *
00365  *        r[]   - radiuses                                             *
00366  *        dphi  - delta phi                                            *
00367  *        nds    - number of discrete steps                            *
00368  *        kface  - current free cell in the pF array                   *
00369  *                                                                     *
00370  ***********************************************************************/
00371 {
00372   G4int k1, k2, k3, k4;
00373   
00374   if (std::abs((G4double)(dphi-pi)) < perMillion) {          // half a circle
00375     for (G4int i=0; i<4; i++) {
00376       k1 = ii[i];
00377       k2 = (i == 3) ? ii[0] : ii[i+1];
00378       if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;      
00379     }
00380   }
00381 
00382   if (ii[1] == ii[2]) {
00383     k1 = kk[ii[0]];
00384     k2 = kk[ii[2]];
00385     k3 = kk[ii[3]];
00386     pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
00387     if (r[ii[0]] != 0.) k1 += nds;
00388     if (r[ii[2]] != 0.) k2 += nds;
00389     if (r[ii[3]] != 0.) k3 += nds;
00390     pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
00391   }else if (kk[ii[0]] == kk[ii[1]]) {
00392     k1 = kk[ii[0]];
00393     k2 = kk[ii[2]];
00394     k3 = kk[ii[3]];
00395     pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
00396     if (r[ii[0]] != 0.) k1 += nds;
00397     if (r[ii[2]] != 0.) k2 += nds;
00398     if (r[ii[3]] != 0.) k3 += nds;
00399     pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
00400   }else if (kk[ii[2]] == kk[ii[3]]) {
00401     k1 = kk[ii[0]];
00402     k2 = kk[ii[1]];
00403     k3 = kk[ii[2]];
00404     pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
00405     if (r[ii[0]] != 0.) k1 += nds;
00406     if (r[ii[1]] != 0.) k2 += nds;
00407     if (r[ii[2]] != 0.) k3 += nds;
00408     pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
00409   }else{
00410     k1 = kk[ii[0]];
00411     k2 = kk[ii[1]];
00412     k3 = kk[ii[2]];
00413     k4 = kk[ii[3]];
00414     pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
00415     if (r[ii[0]] != 0.) k1 += nds;
00416     if (r[ii[1]] != 0.) k2 += nds;
00417     if (r[ii[2]] != 0.) k3 += nds;
00418     if (r[ii[3]] != 0.) k4 += nds;
00419     pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
00420   }
00421 }
00422 
00423 void HepPolyhedron::RotateAroundZ(G4int nstep, G4double phi, G4double dphi,
00424                                  G4int np1, G4int np2,
00425                                  const G4double *z, G4double *r,
00426                                  G4int nodeVis, G4int edgeVis)
00427 /***********************************************************************
00428  *                                                                     *
00429  * Name: HepPolyhedron::RotateAroundZ                Date:    27.11.96 *
00430  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
00431  *                                                                     *
00432  * Function: Create HepPolyhedron for a solid produced by rotation of  *
00433  *           two polylines around Z-axis                               *
00434  *                                                                     *
00435  * Input: nstep - number of discrete steps, if 0 then default          *
00436  *        phi   - starting phi angle                                   *
00437  *        dphi  - delta phi                                            *
00438  *        np1   - number of points in external polyline                *
00439  *                (must be negative in case of closed polyline)        *
00440  *        np2   - number of points in internal polyline (may be 1)     *
00441  *        z[]   - z-coordinates (+z >>> -z for both polylines)         *
00442  *        r[]   - r-coordinates                                        *
00443  *        nodeVis - how to Draw edges joing consecutive positions of   *
00444  *                  node during rotation                               *
00445  *        edgeVis - how to Draw edges                                  *
00446  *                                                                     *
00447  ***********************************************************************/
00448 {
00449   static const G4double wholeCircle = twopi;
00450     
00451   //   S E T   R O T A T I O N   P A R A M E T E R S
00452 
00453   G4bool ifWholeCircle = (std::abs(dphi-wholeCircle) < perMillion) ? true : false;
00454   G4double   delPhi  = ifWholeCircle ? wholeCircle : dphi;  
00455   G4int        nSphi    = (nstep > 0) ?
00456     nstep : G4int(delPhi*GetNumberOfRotationSteps()/wholeCircle+.5);
00457   if (nSphi == 0) nSphi = 1;
00458   G4int        nVphi    = ifWholeCircle ? nSphi : nSphi+1;
00459   G4bool ifClosed = np1 > 0 ? false : true;
00460   
00461   //   C O U N T   V E R T E C E S
00462 
00463   G4int absNp1 = std::abs(np1);
00464   G4int absNp2 = std::abs(np2);
00465   G4int i1beg = 0;
00466   G4int i1end = absNp1-1;
00467   G4int i2beg = absNp1;
00468   G4int i2end = absNp1+absNp2-1; 
00469   G4int i, j, k;
00470 
00471   for(i=i1beg; i<=i2end; i++) {
00472     if (std::abs(r[i]) < spatialTolerance) r[i] = 0.;
00473   }
00474 
00475   j = 0;                                                // external nodes
00476   for (i=i1beg; i<=i1end; i++) {
00477     j += (r[i] == 0.) ? 1 : nVphi;
00478   }
00479 
00480   G4bool ifSide1 = false;                           // internal nodes
00481   G4bool ifSide2 = false;
00482 
00483   if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
00484     j += (r[i2beg] == 0.) ? 1 : nVphi;
00485     ifSide1 = true;
00486   }
00487 
00488   for(i=i2beg+1; i<i2end; i++) {
00489     j += (r[i] == 0.) ? 1 : nVphi;
00490   }
00491   
00492   if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
00493     if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
00494     ifSide2 = true;
00495   }
00496 
00497   //   C O U N T   F A C E S
00498 
00499   k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi;       // external faces
00500 
00501   if (absNp2 > 1) {                                     // internal faces
00502     for(i=i2beg; i<i2end; i++) {
00503       if (r[i] > 0. || r[i+1] > 0.)       k += nSphi;
00504     }
00505 
00506     if (ifClosed) {
00507       if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
00508     }
00509   }
00510 
00511   if (!ifClosed) {                                      // side faces
00512     if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
00513     if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) k += nSphi;
00514   }
00515 
00516   if (!ifWholeCircle) {                                 // phi_side faces
00517     k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
00518   }
00519 
00520   //   A L L O C A T E   M E M O R Y
00521 
00522   AllocateMemory(j, k);
00523 
00524   //   G E N E R A T E   V E R T E C E S
00525 
00526   G4int *kk;
00527   kk = new G4int[absNp1+absNp2];
00528 
00529   k = 1;
00530   for(i=i1beg; i<=i1end; i++) {
00531     kk[i] = k;
00532     if (r[i] == 0.)
00533     { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
00534   }
00535 
00536   i = i2beg;
00537   if (ifSide1) {
00538     kk[i] = k;
00539     if (r[i] == 0.)
00540     { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
00541   }else{
00542     kk[i] = kk[i1beg];
00543   }
00544 
00545   for(i=i2beg+1; i<i2end; i++) {
00546     kk[i] = k;
00547     if (r[i] == 0.)
00548     { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
00549   }
00550 
00551   if (absNp2 > 1) {
00552     i = i2end;
00553     if (ifSide2) {
00554       kk[i] = k;
00555       if (r[i] == 0.) pV[k] = G4Point3D(0, 0, z[i]);
00556     }else{
00557       kk[i] = kk[i1end];
00558     }
00559   }
00560 
00561   G4double cosPhi, sinPhi;
00562 
00563   for(j=0; j<nVphi; j++) {
00564     cosPhi = std::cos(phi+j*delPhi/nSphi);
00565     sinPhi = std::sin(phi+j*delPhi/nSphi);
00566     for(i=i1beg; i<=i2end; i++) {
00567       if (r[i] != 0.)
00568         pV[kk[i]+j] = G4Point3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
00569     }
00570   }
00571 
00572   //   G E N E R A T E   E X T E R N A L   F A C E S
00573 
00574   G4int v1,v2;
00575 
00576   k = 1;
00577   v2 = ifClosed ? nodeVis : 1;
00578   for(i=i1beg; i<i1end; i++) {
00579     v1 = v2;
00580     if (!ifClosed && i == i1end-1) {
00581       v2 = 1;
00582     }else{
00583       v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
00584     }
00585     RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
00586                edgeVis, ifWholeCircle, nSphi, k);
00587   }
00588   if (ifClosed) {
00589     RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
00590                edgeVis, ifWholeCircle, nSphi, k);
00591   }
00592 
00593   //   G E N E R A T E   I N T E R N A L   F A C E S
00594 
00595   if (absNp2 > 1) {
00596     v2 = ifClosed ? nodeVis : 1;
00597     for(i=i2beg; i<i2end; i++) {
00598       v1 = v2;
00599       if (!ifClosed && i==i2end-1) {
00600         v2 = 1;
00601       }else{
00602         v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 :  nodeVis;
00603       }
00604       RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
00605                  edgeVis, ifWholeCircle, nSphi, k);
00606     }
00607     if (ifClosed) {
00608       RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
00609                  edgeVis, ifWholeCircle, nSphi, k);
00610     }
00611   }
00612 
00613   //   G E N E R A T E   S I D E   F A C E S
00614 
00615   if (!ifClosed) {
00616     if (ifSide1) {
00617       RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
00618                  -1, ifWholeCircle, nSphi, k);
00619     }
00620     if (ifSide2) {
00621       RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
00622                  -1, ifWholeCircle, nSphi, k);
00623     }
00624   }
00625 
00626   //   G E N E R A T E   S I D E   F A C E S  for the case of incomplete circle
00627 
00628   if (!ifWholeCircle) {
00629 
00630     G4int  ii[4], vv[4];
00631 
00632     if (ifClosed) {
00633       for (i=i1beg; i<=i1end; i++) {
00634         ii[0] = i;
00635         ii[3] = (i == i1end) ? i1beg : i+1;
00636         ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
00637         ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
00638         vv[0] = -1;
00639         vv[1] = 1;
00640         vv[2] = -1;
00641         vv[3] = 1;
00642         SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
00643       }
00644     }else{
00645       for (i=i1beg; i<i1end; i++) {
00646         ii[0] = i;
00647         ii[3] = i+1;
00648         ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
00649         ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
00650         vv[0] = (i == i1beg)   ? 1 : -1;
00651         vv[1] = 1;
00652         vv[2] = (i == i1end-1) ? 1 : -1;
00653         vv[3] = 1;
00654         SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
00655       }
00656     }      
00657   }
00658 
00659   delete [] kk;
00660 
00661   if (k-1 != nface) {
00662     std::cerr
00663       << "Polyhedron::RotateAroundZ: number of generated faces ("
00664       << k-1 << ") is not equal to the number of allocated faces ("
00665       << nface << ")"
00666       << std::endl;
00667   }
00668 }
00669 
00670 void HepPolyhedron::SetReferences()
00671 /***********************************************************************
00672  *                                                                     *
00673  * Name: HepPolyhedron::SetReferences                Date:    04.12.96 *
00674  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
00675  *                                                                     *
00676  * Function: For each edge set reference to neighbouring facet         *
00677  *                                                                     *
00678  ***********************************************************************/
00679 {
00680   if (nface <= 0) return;
00681 
00682   struct edgeListMember {
00683     edgeListMember *next;
00684     G4int v2;
00685     G4int iface;
00686     G4int iedge;
00687   } *edgeList, *freeList, **headList;
00688 
00689   
00690   //   A L L O C A T E   A N D   I N I T I A T E   L I S T S
00691 
00692   edgeList = new edgeListMember[2*nface];
00693   headList = new edgeListMember*[nvert];
00694   
00695   G4int i;
00696   for (i=0; i<nvert; i++) {
00697     headList[i] = 0;
00698   }
00699   freeList = edgeList;
00700   for (i=0; i<2*nface-1; i++) {
00701     edgeList[i].next = &edgeList[i+1];
00702   }
00703   edgeList[2*nface-1].next = 0;
00704 
00705   //   L O O P   A L O N G   E D G E S
00706 
00707   G4int iface, iedge, nedge, i1, i2, k1, k2;
00708   edgeListMember *prev, *cur;
00709   
00710   for(iface=1; iface<=nface; iface++) {
00711     nedge = (pF[iface].edge[3].v == 0) ? 3 : 4;
00712     for (iedge=0; iedge<nedge; iedge++) {
00713       i1 = iedge;
00714       i2 = (iedge < nedge-1) ? iedge+1 : 0;
00715       i1 = std::abs(pF[iface].edge[i1].v);
00716       i2 = std::abs(pF[iface].edge[i2].v);
00717       k1 = (i1 < i2) ? i1 : i2;          // k1 = ::min(i1,i2);
00718       k2 = (i1 > i2) ? i1 : i2;          // k2 = ::max(i1,i2);
00719       
00720       // check head of the List corresponding to k1
00721       cur = headList[k1];
00722       if (cur == 0) {
00723         headList[k1] = freeList;
00724         freeList = freeList->next;
00725         cur = headList[k1];
00726         cur->next = 0;
00727         cur->v2 = k2;
00728         cur->iface = iface;
00729         cur->iedge = iedge;
00730         continue;
00731       }
00732 
00733       if (cur->v2 == k2) {
00734         headList[k1] = cur->next;
00735         cur->next = freeList;
00736         freeList = cur;      
00737         pF[iface].edge[iedge].f = cur->iface;
00738         pF[cur->iface].edge[cur->iedge].f = iface;
00739         i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
00740         i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
00741         if (i1 != i2) {
00742           std::cerr
00743             << "Polyhedron::SetReferences: different edge visibility "
00744             << iface << "/" << iedge << "/"
00745             << pF[iface].edge[iedge].v << " and "
00746             << cur->iface << "/" << cur->iedge << "/"
00747             << pF[cur->iface].edge[cur->iedge].v
00748             << std::endl;
00749         }
00750         continue;
00751       }
00752 
00753       // check List itself
00754       for (;;) {
00755         prev = cur;
00756         cur = prev->next;
00757         if (cur == 0) {
00758           prev->next = freeList;
00759           freeList = freeList->next;
00760           cur = prev->next;
00761           cur->next = 0;
00762           cur->v2 = k2;
00763           cur->iface = iface;
00764           cur->iedge = iedge;
00765           break;
00766         }
00767 
00768         if (cur->v2 == k2) {
00769           prev->next = cur->next;
00770           cur->next = freeList;
00771           freeList = cur;      
00772           pF[iface].edge[iedge].f = cur->iface;
00773           pF[cur->iface].edge[cur->iedge].f = iface;
00774           i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
00775           i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
00776             if (i1 != i2) {
00777               std::cerr
00778                 << "Polyhedron::SetReferences: different edge visibility "
00779                 << iface << "/" << iedge << "/"
00780                 << pF[iface].edge[iedge].v << " and "
00781                 << cur->iface << "/" << cur->iedge << "/"
00782                 << pF[cur->iface].edge[cur->iedge].v
00783                 << std::endl;
00784             }
00785           break;
00786         }
00787       }
00788     }
00789   }
00790 
00791   //  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
00792 
00793   for (i=0; i<nvert; i++) {
00794     if (headList[i] != 0) {
00795       std::cerr
00796         << "Polyhedron::SetReferences: List " << i << " is not empty"
00797         << std::endl;
00798     }
00799   }
00800 
00801   //   F R E E   M E M O R Y
00802 
00803   delete [] edgeList;
00804   delete [] headList;
00805 }
00806 
00807 void HepPolyhedron::InvertFacets()
00808 /***********************************************************************
00809  *                                                                     *
00810  * Name: HepPolyhedron::InvertFacets                Date:    01.12.99  *
00811  * Author: E.Chernyaev                              Revised:           *
00812  *                                                                     *
00813  * Function: Invert the order of the nodes in the facets               *
00814  *                                                                     *
00815  ***********************************************************************/
00816 {
00817   if (nface <= 0) return;
00818   G4int i, k, nnode, v[4],f[4];
00819   for (i=1; i<=nface; i++) {
00820     nnode =  (pF[i].edge[3].v == 0) ? 3 : 4;
00821     for (k=0; k<nnode; k++) {
00822       v[k] = (k+1 == nnode) ? pF[i].edge[0].v : pF[i].edge[k+1].v;
00823       if (v[k] * pF[i].edge[k].v < 0) v[k] = -v[k];
00824       f[k] = pF[i].edge[k].f;
00825     }
00826     for (k=0; k<nnode; k++) {
00827       pF[i].edge[nnode-1-k].v = v[k];
00828       pF[i].edge[nnode-1-k].f = f[k];
00829     }
00830   }
00831 }
00832 
00833 HepPolyhedron & HepPolyhedron::Transform(const G4Transform3D &t)
00834 /***********************************************************************
00835  *                                                                     *
00836  * Name: HepPolyhedron::Transform                    Date:    01.12.99  *
00837  * Author: E.Chernyaev                              Revised:           *
00838  *                                                                     *
00839  * Function: Make transformation of the polyhedron                     *
00840  *                                                                     *
00841  ***********************************************************************/
00842 {
00843   if (nvert > 0) {
00844     for (G4int i=1; i<=nvert; i++) { pV[i] = t * pV[i]; }
00845 
00846     //  C H E C K   D E T E R M I N A N T   A N D
00847     //  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
00848 
00849     G4Vector3D d = t * G4Vector3D(0,0,0);
00850     G4Vector3D x = t * G4Vector3D(1,0,0) - d;
00851     G4Vector3D y = t * G4Vector3D(0,1,0) - d;
00852     G4Vector3D z = t * G4Vector3D(0,0,1) - d;
00853     if ((x.cross(y))*z < 0) InvertFacets();
00854   }
00855   return *this;
00856 }
00857 
00858 G4bool HepPolyhedron::GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
00859 /***********************************************************************
00860  *                                                                     *
00861  * Name: HepPolyhedron::GetNextVertexIndex          Date:    03.09.96  *
00862  * Author: Yasuhide Sawada                          Revised:           *
00863  *                                                                     *
00864  * Function:                                                           *
00865  *                                                                     *
00866  ***********************************************************************/
00867 {
00868   static G4int iFace = 1;
00869   static G4int iQVertex = 0;
00870   G4int vIndex = pF[iFace].edge[iQVertex].v;
00871 
00872   edgeFlag = (vIndex > 0) ? 1 : 0;
00873   index = std::abs(vIndex);
00874 
00875   if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
00876     iQVertex = 0;
00877     if (++iFace > nface) iFace = 1;
00878     return false;  // Last Edge
00879   }else{
00880     ++iQVertex;
00881     return true;  // not Last Edge
00882   }
00883 }
00884 
00885 G4Point3D HepPolyhedron::GetVertex(G4int index) const
00886 /***********************************************************************
00887  *                                                                     *
00888  * Name: HepPolyhedron::GetVertex                   Date:    03.09.96  *
00889  * Author: Yasuhide Sawada                          Revised: 17.11.99  *
00890  *                                                                     *
00891  * Function: Get vertex of the index.                                  *
00892  *                                                                     *
00893  ***********************************************************************/
00894 {
00895   if (index <= 0 || index > nvert) {
00896     std::cerr
00897       << "HepPolyhedron::GetVertex: irrelevant index " << index
00898       << std::endl;
00899     return G4Point3D();
00900   }
00901   return pV[index];
00902 }
00903 
00904 G4bool
00905 HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag) const
00906 /***********************************************************************
00907  *                                                                     *
00908  * Name: HepPolyhedron::GetNextVertex               Date:    22.07.96  *
00909  * Author: John Allison                             Revised:           *
00910  *                                                                     *
00911  * Function: Get vertices of the quadrilaterals in order for each      *
00912  *           face in face order.  Returns false when finished each     *
00913  *           face.                                                     *
00914  *                                                                     *
00915  ***********************************************************************/
00916 {
00917   G4int index;
00918   G4bool rep = GetNextVertexIndex(index, edgeFlag);
00919   vertex = pV[index];
00920   return rep;
00921 }
00922 
00923 G4bool HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag,
00924                                   G4Normal3D &normal) const
00925 /***********************************************************************
00926  *                                                                     *
00927  * Name: HepPolyhedron::GetNextVertex               Date:    26.11.99  *
00928  * Author: E.Chernyaev                              Revised:           *
00929  *                                                                     *
00930  * Function: Get vertices with normals of the quadrilaterals in order  *
00931  *           for each face in face order.                              *
00932  *           Returns false when finished each face.                    *
00933  *                                                                     *
00934  ***********************************************************************/
00935 {
00936   static G4int iFace = 1;
00937   static G4int iNode = 0;
00938 
00939   if (nface == 0) return false;  // empty polyhedron
00940 
00941   G4int k = pF[iFace].edge[iNode].v;
00942   if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
00943   vertex = pV[k];
00944   normal = FindNodeNormal(iFace,k);
00945   if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
00946     iNode = 0;
00947     if (++iFace > nface) iFace = 1;
00948     return false;                // last node
00949   }else{
00950     ++iNode;
00951     return true;                 // not last node
00952   }
00953 }
00954 
00955 G4bool HepPolyhedron::GetNextEdgeIndeces(G4int &i1, G4int &i2, G4int &edgeFlag,
00956                                        G4int &iface1, G4int &iface2) const
00957 /***********************************************************************
00958  *                                                                     *
00959  * Name: HepPolyhedron::GetNextEdgeIndeces          Date:    30.09.96  *
00960  * Author: E.Chernyaev                              Revised: 17.11.99  *
00961  *                                                                     *
00962  * Function: Get indeces of the next edge together with indeces of     *
00963  *           of the faces which share the edge.                        *
00964  *           Returns false when the last edge.                         *
00965  *                                                                     *
00966  ***********************************************************************/
00967 {
00968   static G4int iFace    = 1;
00969   static G4int iQVertex = 0;
00970   static G4int iOrder   = 1;
00971   G4int  k1, k2, kflag, kface1, kface2;
00972 
00973   if (iFace == 1 && iQVertex == 0) {
00974     k2 = pF[nface].edge[0].v;
00975     k1 = pF[nface].edge[3].v;
00976     if (k1 == 0) k1 = pF[nface].edge[2].v;
00977     if (std::abs(k1) > std::abs(k2)) iOrder = -1;
00978   }
00979 
00980   do {
00981     k1     = pF[iFace].edge[iQVertex].v;
00982     kflag  = k1;
00983     k1     = std::abs(k1);
00984     kface1 = iFace; 
00985     kface2 = pF[iFace].edge[iQVertex].f;
00986     if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
00987       iQVertex = 0;
00988       k2 = std::abs(pF[iFace].edge[iQVertex].v);
00989       iFace++;
00990     }else{
00991       iQVertex++;
00992       k2 = std::abs(pF[iFace].edge[iQVertex].v);
00993     }
00994   } while (iOrder*k1 > iOrder*k2);
00995 
00996   i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
00997   iface1 = kface1; iface2 = kface2; 
00998 
00999   if (iFace > nface) {
01000     iFace  = 1; iOrder = 1;
01001     return false;
01002   }else{
01003     return true;
01004   }
01005 }
01006 
01007 G4bool
01008 HepPolyhedron::GetNextEdgeIndeces(G4int &i1, G4int &i2, G4int &edgeFlag) const
01009 /***********************************************************************
01010  *                                                                     *
01011  * Name: HepPolyhedron::GetNextEdgeIndeces          Date:    17.11.99  *
01012  * Author: E.Chernyaev                              Revised:           *
01013  *                                                                     *
01014  * Function: Get indeces of the next edge.                             *
01015  *           Returns false when the last edge.                         *
01016  *                                                                     *
01017  ***********************************************************************/
01018 {
01019   G4int kface1, kface2;
01020   return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
01021 }
01022 
01023 G4bool
01024 HepPolyhedron::GetNextEdge(G4Point3D &p1,
01025                            G4Point3D &p2,
01026                            G4int &edgeFlag) const
01027 /***********************************************************************
01028  *                                                                     *
01029  * Name: HepPolyhedron::GetNextEdge                 Date:    30.09.96  *
01030  * Author: E.Chernyaev                              Revised:           *
01031  *                                                                     *
01032  * Function: Get next edge.                                            *
01033  *           Returns false when the last edge.                         *
01034  *                                                                     *
01035  ***********************************************************************/
01036 {
01037   G4int i1,i2;
01038   G4bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
01039   p1 = pV[i1];
01040   p2 = pV[i2];
01041   return rep;
01042 }
01043 
01044 G4bool
01045 HepPolyhedron::GetNextEdge(G4Point3D &p1, G4Point3D &p2,
01046                           G4int &edgeFlag, G4int &iface1, G4int &iface2) const
01047 /***********************************************************************
01048  *                                                                     *
01049  * Name: HepPolyhedron::GetNextEdge                 Date:    17.11.99  *
01050  * Author: E.Chernyaev                              Revised:           *
01051  *                                                                     *
01052  * Function: Get next edge with indeces of the faces which share       *
01053  *           the edge.                                                 *
01054  *           Returns false when the last edge.                         *
01055  *                                                                     *
01056  ***********************************************************************/
01057 {
01058   G4int i1,i2;
01059   G4bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
01060   p1 = pV[i1];
01061   p2 = pV[i2];
01062   return rep;
01063 }
01064 
01065 void HepPolyhedron::GetFacet(G4int iFace, G4int &n, G4int *iNodes,
01066                             G4int *edgeFlags, G4int *iFaces) const
01067 /***********************************************************************
01068  *                                                                     *
01069  * Name: HepPolyhedron::GetFacet                    Date:    15.12.99  *
01070  * Author: E.Chernyaev                              Revised:           *
01071  *                                                                     *
01072  * Function: Get face by index                                         *
01073  *                                                                     *
01074  ***********************************************************************/
01075 {
01076   if (iFace < 1 || iFace > nface) {
01077     std::cerr 
01078       << "HepPolyhedron::GetFacet: irrelevant index " << iFace
01079       << std::endl;
01080     n = 0;
01081   }else{
01082     G4int i, k;
01083     for (i=0; i<4; i++) { 
01084       k = pF[iFace].edge[i].v;
01085       if (k == 0) break;
01086       if (iFaces != 0) iFaces[i] = pF[iFace].edge[i].f;
01087       if (k > 0) { 
01088         iNodes[i] = k;
01089         if (edgeFlags != 0) edgeFlags[i] = 1;
01090       }else{
01091         iNodes[i] = -k;
01092         if (edgeFlags != 0) edgeFlags[i] = -1;
01093       }
01094     }
01095     n = i;
01096   }
01097 }
01098 
01099 void HepPolyhedron::GetFacet(G4int index, G4int &n, G4Point3D *nodes,
01100                              G4int *edgeFlags, G4Normal3D *normals) const
01101 /***********************************************************************
01102  *                                                                     *
01103  * Name: HepPolyhedron::GetFacet                    Date:    17.11.99  *
01104  * Author: E.Chernyaev                              Revised:           *
01105  *                                                                     *
01106  * Function: Get face by index                                         *
01107  *                                                                     *
01108  ***********************************************************************/
01109 {
01110   G4int iNodes[4];
01111   GetFacet(index, n, iNodes, edgeFlags);
01112   if (n != 0) {
01113     for (G4int i=0; i<n; i++) {
01114       nodes[i] = pV[iNodes[i]];
01115       if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
01116     }
01117   }
01118 }
01119 
01120 G4bool
01121 HepPolyhedron::GetNextFacet(G4int &n, G4Point3D *nodes,
01122                            G4int *edgeFlags, G4Normal3D *normals) const
01123 /***********************************************************************
01124  *                                                                     *
01125  * Name: HepPolyhedron::GetNextFacet                Date:    19.11.99  *
01126  * Author: E.Chernyaev                              Revised:           *
01127  *                                                                     *
01128  * Function: Get next face with normals of unit length at the nodes.   *
01129  *           Returns false when finished all faces.                    *
01130  *                                                                     *
01131  ***********************************************************************/
01132 {
01133   static G4int iFace = 1;
01134 
01135   if (edgeFlags == 0) {
01136     GetFacet(iFace, n, nodes);
01137   }else if (normals == 0) {
01138     GetFacet(iFace, n, nodes, edgeFlags);
01139   }else{
01140     GetFacet(iFace, n, nodes, edgeFlags, normals);
01141   }
01142 
01143   if (++iFace > nface) {
01144     iFace  = 1;
01145     return false;
01146   }else{
01147     return true;
01148   }
01149 }
01150 
01151 G4Normal3D HepPolyhedron::GetNormal(G4int iFace) const
01152 /***********************************************************************
01153  *                                                                     *
01154  * Name: HepPolyhedron::GetNormal                    Date:    19.11.99 *
01155  * Author: E.Chernyaev                               Revised:          *
01156  *                                                                     *
01157  * Function: Get normal of the face given by index                     *
01158  *                                                                     *
01159  ***********************************************************************/
01160 {
01161   if (iFace < 1 || iFace > nface) {
01162     std::cerr 
01163       << "HepPolyhedron::GetNormal: irrelevant index " << iFace 
01164       << std::endl;
01165     return G4Normal3D();
01166   }
01167 
01168   G4int i0  = std::abs(pF[iFace].edge[0].v);
01169   G4int i1  = std::abs(pF[iFace].edge[1].v);
01170   G4int i2  = std::abs(pF[iFace].edge[2].v);
01171   G4int i3  = std::abs(pF[iFace].edge[3].v);
01172   if (i3 == 0) i3 = i0;
01173   return (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1]);
01174 }
01175 
01176 G4Normal3D HepPolyhedron::GetUnitNormal(G4int iFace) const
01177 /***********************************************************************
01178  *                                                                     *
01179  * Name: HepPolyhedron::GetNormal                    Date:    19.11.99 *
01180  * Author: E.Chernyaev                               Revised:          *
01181  *                                                                     *
01182  * Function: Get unit normal of the face given by index                *
01183  *                                                                     *
01184  ***********************************************************************/
01185 {
01186   if (iFace < 1 || iFace > nface) {
01187     std::cerr 
01188       << "HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
01189       << std::endl;
01190     return G4Normal3D();
01191   }
01192 
01193   G4int i0  = std::abs(pF[iFace].edge[0].v);
01194   G4int i1  = std::abs(pF[iFace].edge[1].v);
01195   G4int i2  = std::abs(pF[iFace].edge[2].v);
01196   G4int i3  = std::abs(pF[iFace].edge[3].v);
01197   if (i3 == 0) i3 = i0;
01198   return ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).unit();
01199 }
01200 
01201 G4bool HepPolyhedron::GetNextNormal(G4Normal3D &normal) const
01202 /***********************************************************************
01203  *                                                                     *
01204  * Name: HepPolyhedron::GetNextNormal               Date:    22.07.96  *
01205  * Author: John Allison                             Revised: 19.11.99  *
01206  *                                                                     *
01207  * Function: Get normals of each face in face order.  Returns false    *
01208  *           when finished all faces.                                  *
01209  *                                                                     *
01210  ***********************************************************************/
01211 {
01212   static G4int iFace = 1;
01213   normal = GetNormal(iFace);
01214   if (++iFace > nface) {
01215     iFace = 1;
01216     return false;
01217   }else{
01218     return true;
01219   }
01220 }
01221 
01222 G4bool HepPolyhedron::GetNextUnitNormal(G4Normal3D &normal) const
01223 /***********************************************************************
01224  *                                                                     *
01225  * Name: HepPolyhedron::GetNextUnitNormal           Date:    16.09.96  *
01226  * Author: E.Chernyaev                              Revised:           *
01227  *                                                                     *
01228  * Function: Get normals of unit length of each face in face order.    *
01229  *           Returns false when finished all faces.                    *
01230  *                                                                     *
01231  ***********************************************************************/
01232 {
01233   G4bool rep = GetNextNormal(normal);
01234   normal = normal.unit();
01235   return rep;
01236 }
01237 
01238 G4double HepPolyhedron::GetSurfaceArea() const
01239 /***********************************************************************
01240  *                                                                     *
01241  * Name: HepPolyhedron::GetSurfaceArea              Date:    25.05.01  *
01242  * Author: E.Chernyaev                              Revised:           *
01243  *                                                                     *
01244  * Function: Returns area of the surface of the polyhedron.            *
01245  *                                                                     *
01246  ***********************************************************************/
01247 {
01248   G4double srf = 0.;
01249   for (G4int iFace=1; iFace<=nface; iFace++) {
01250     G4int i0 = std::abs(pF[iFace].edge[0].v);
01251     G4int i1 = std::abs(pF[iFace].edge[1].v);
01252     G4int i2 = std::abs(pF[iFace].edge[2].v);
01253     G4int i3 = std::abs(pF[iFace].edge[3].v);
01254     if (i3 == 0) i3 = i0;
01255     srf += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).mag();
01256   }
01257   return srf/2.;
01258 }
01259 
01260 G4double HepPolyhedron::GetVolume() const
01261 /***********************************************************************
01262  *                                                                     *
01263  * Name: HepPolyhedron::GetVolume                   Date:    25.05.01  *
01264  * Author: E.Chernyaev                              Revised:           *
01265  *                                                                     *
01266  * Function: Returns volume of the polyhedron.                         *
01267  *                                                                     *
01268  ***********************************************************************/
01269 {
01270   G4double v = 0.;
01271   for (G4int iFace=1; iFace<=nface; iFace++) {
01272     G4int i0 = std::abs(pF[iFace].edge[0].v);
01273     G4int i1 = std::abs(pF[iFace].edge[1].v);
01274     G4int i2 = std::abs(pF[iFace].edge[2].v);
01275     G4int i3 = std::abs(pF[iFace].edge[3].v);
01276     G4Point3D pt;
01277     if (i3 == 0) {
01278       i3 = i0;
01279       pt = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);
01280     }else{
01281       pt = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0.25;
01282     }
01283     v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).dot(pt);
01284   }
01285   return v/6.;
01286 }
01287 
01288 G4int
01289 HepPolyhedron::createTwistedTrap(G4double Dz,
01290                                  const G4double xy1[][2],
01291                                  const G4double xy2[][2])
01292 /***********************************************************************
01293  *                                                                     *
01294  * Name: createTwistedTrap                           Date:    05.11.02 *
01295  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
01296  *                                                                     *
01297  * Function: Creates polyhedron for twisted trapezoid                  *
01298  *                                                                     *
01299  * Input: Dz       - half-length along Z             8----7            *
01300  *        xy1[2,4] - quadrilateral at Z=-Dz       5----6  !            *
01301  *        xy2[2,4] - quadrilateral at Z=+Dz       !  4-!--3            *
01302  *                                                1----2               *
01303  *                                                                     *
01304  ***********************************************************************/
01305 {
01306   AllocateMemory(12,18);
01307 
01308   pV[ 1] = G4Point3D(xy1[0][0],xy1[0][1],-Dz);
01309   pV[ 2] = G4Point3D(xy1[1][0],xy1[1][1],-Dz);
01310   pV[ 3] = G4Point3D(xy1[2][0],xy1[2][1],-Dz);
01311   pV[ 4] = G4Point3D(xy1[3][0],xy1[3][1],-Dz);
01312 
01313   pV[ 5] = G4Point3D(xy2[0][0],xy2[0][1], Dz);
01314   pV[ 6] = G4Point3D(xy2[1][0],xy2[1][1], Dz);
01315   pV[ 7] = G4Point3D(xy2[2][0],xy2[2][1], Dz);
01316   pV[ 8] = G4Point3D(xy2[3][0],xy2[3][1], Dz);
01317 
01318   pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;
01319   pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;
01320   pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;
01321   pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;
01322 
01323   enum {DUMMY, BOTTOM,
01324         LEFT_BOTTOM,  LEFT_FRONT,   LEFT_TOP,  LEFT_BACK,
01325         BACK_BOTTOM,  BACK_LEFT,    BACK_TOP,  BACK_RIGHT,
01326         RIGHT_BOTTOM, RIGHT_BACK,   RIGHT_TOP, RIGHT_FRONT,
01327         FRONT_BOTTOM, FRONT_RIGHT,  FRONT_TOP, FRONT_LEFT,
01328         TOP};
01329 
01330   pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
01331 
01332   pF[ 2]=G4Facet(4,BOTTOM,     -1,LEFT_FRONT,  -12,LEFT_BACK,    0,0);
01333   pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP,    -12,LEFT_BOTTOM,  0,0);
01334   pF[ 4]=G4Facet(5,TOP,        -8,LEFT_BACK,   -12,LEFT_FRONT,   0,0);
01335   pF[ 5]=G4Facet(8,BACK_LEFT,  -4,LEFT_BOTTOM, -12,LEFT_TOP,     0,0);
01336 
01337   pF[ 6]=G4Facet(3,BOTTOM,     -4,BACK_LEFT,   -11,BACK_RIGHT,   0,0);
01338   pF[ 7]=G4Facet(4,LEFT_BACK,  -8,BACK_TOP,    -11,BACK_BOTTOM,  0,0);
01339   pF[ 8]=G4Facet(8,TOP,        -7,BACK_RIGHT,  -11,BACK_LEFT,    0,0);
01340   pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP,     0,0);
01341 
01342   pF[10]=G4Facet(2,BOTTOM,     -3,RIGHT_BACK,  -10,RIGHT_FRONT,  0,0);
01343   pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP,   -10,RIGHT_BOTTOM, 0,0);
01344   pF[12]=G4Facet(7,TOP,        -6,RIGHT_FRONT, -10,RIGHT_BACK,   0,0);
01345   pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP,    0,0);
01346 
01347   pF[14]=G4Facet(1,BOTTOM,     -2,FRONT_RIGHT,  -9,FRONT_LEFT,   0,0);
01348   pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP,    -9,FRONT_BOTTOM, 0,0);
01349   pF[16]=G4Facet(6,TOP,        -5,FRONT_LEFT,   -9,FRONT_RIGHT,  0,0);
01350   pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP,    0,0);
01351  
01352   pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
01353 
01354   return 0;
01355 }
01356 
01357 G4int
01358 HepPolyhedron::createPolyhedron(G4int Nnodes, G4int Nfaces,
01359                                 const G4double xyz[][3],
01360                                 const G4int  faces[][4])
01361 /***********************************************************************
01362  *                                                                     *
01363  * Name: createPolyhedron                            Date:    05.11.02 *
01364  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
01365  *                                                                     *
01366  * Function: Creates user defined polyhedron                           *
01367  *                                                                     *
01368  * Input: Nnodes  - number of nodes                                    *
01369  *        Nfaces  - number of faces                                    *
01370  *        nodes[][3] - node coordinates                                *
01371  *        faces[][4] - faces                                           *
01372  *                                                                     *
01373  ***********************************************************************/
01374 {
01375   AllocateMemory(Nnodes, Nfaces);
01376   if (nvert == 0) return 1;
01377 
01378   for (G4int i=0; i<Nnodes; i++) {
01379     pV[i+1] = G4Point3D(xyz[i][0], xyz[i][1], xyz[i][2]);
01380   }
01381   for (G4int k=0; k<Nfaces; k++) {
01382     pF[k+1] = G4Facet(faces[k][0],0,faces[k][1],0,faces[k][2],0,faces[k][3],0);
01383   }
01384   SetReferences();
01385   return 0;
01386 }
01387 
01388 HepPolyhedronTrd2::HepPolyhedronTrd2(G4double Dx1, G4double Dx2,
01389                                      G4double Dy1, G4double Dy2,
01390                                      G4double Dz)
01391 /***********************************************************************
01392  *                                                                     *
01393  * Name: HepPolyhedronTrd2                           Date:    22.07.96 *
01394  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
01395  *                                                                     *
01396  * Function: Create GEANT4 TRD2-trapezoid                              *
01397  *                                                                     *
01398  * Input: Dx1 - half-length along X at -Dz           8----7            *
01399  *        Dx2 - half-length along X ay +Dz        5----6  !            *
01400  *        Dy1 - half-length along Y ay -Dz        !  4-!--3            *
01401  *        Dy2 - half-length along Y ay +Dz        1----2               *
01402  *        Dz  - half-length along Z                                    *
01403  *                                                                     *
01404  ***********************************************************************/
01405 {
01406   AllocateMemory(8,6);
01407 
01408   pV[1] = G4Point3D(-Dx1,-Dy1,-Dz);
01409   pV[2] = G4Point3D( Dx1,-Dy1,-Dz);
01410   pV[3] = G4Point3D( Dx1, Dy1,-Dz);
01411   pV[4] = G4Point3D(-Dx1, Dy1,-Dz);
01412   pV[5] = G4Point3D(-Dx2,-Dy2, Dz);
01413   pV[6] = G4Point3D( Dx2,-Dy2, Dz);
01414   pV[7] = G4Point3D( Dx2, Dy2, Dz);
01415   pV[8] = G4Point3D(-Dx2, Dy2, Dz);
01416 
01417   CreatePrism();
01418 }
01419 
01420 HepPolyhedronTrd2::~HepPolyhedronTrd2() {}
01421 
01422 HepPolyhedronTrd1::HepPolyhedronTrd1(G4double Dx1, G4double Dx2,
01423                                      G4double Dy, G4double Dz)
01424   : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
01425 
01426 HepPolyhedronTrd1::~HepPolyhedronTrd1() {}
01427 
01428 HepPolyhedronBox::HepPolyhedronBox(G4double Dx, G4double Dy, G4double Dz)
01429   : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
01430 
01431 HepPolyhedronBox::~HepPolyhedronBox() {}
01432 
01433 HepPolyhedronTrap::HepPolyhedronTrap(G4double Dz,
01434                                      G4double Theta,
01435                                      G4double Phi,
01436                                      G4double Dy1,
01437                                      G4double Dx1,
01438                                      G4double Dx2,
01439                                      G4double Alp1,
01440                                      G4double Dy2,
01441                                      G4double Dx3,
01442                                      G4double Dx4,
01443                                      G4double Alp2)
01444 /***********************************************************************
01445  *                                                                     *
01446  * Name: HepPolyhedronTrap                           Date:    20.11.96 *
01447  * Author: E.Chernyaev                               Revised:          *
01448  *                                                                     *
01449  * Function: Create GEANT4 TRAP-trapezoid                              *
01450  *                                                                     *
01451  * Input: DZ   - half-length in Z                                      *
01452  *        Theta,Phi - polar angles of the line joining centres of the  *
01453  *                    faces at Z=-Dz and Z=+Dz                         *
01454  *        Dy1  - half-length in Y of the face at Z=-Dz                 *
01455  *        Dx1  - half-length in X of low edge of the face at Z=-Dz     *
01456  *        Dx2  - half-length in X of top edge of the face at Z=-Dz     *
01457  *        Alp1 - angle between Y-axis and the median joining top and   *
01458  *               low edges of the face at Z=-Dz                        *
01459  *        Dy2  - half-length in Y of the face at Z=+Dz                 *
01460  *        Dx3  - half-length in X of low edge of the face at Z=+Dz     *
01461  *        Dx4  - half-length in X of top edge of the face at Z=+Dz     *
01462  *        Alp2 - angle between Y-axis and the median joining top and   *
01463  *               low edges of the face at Z=+Dz                        *
01464  *                                                                     *
01465  ***********************************************************************/
01466 {
01467   G4double DzTthetaCphi = Dz*std::tan(Theta)*std::cos(Phi);
01468   G4double DzTthetaSphi = Dz*std::tan(Theta)*std::sin(Phi);
01469   G4double Dy1Talp1 = Dy1*std::tan(Alp1);
01470   G4double Dy2Talp2 = Dy2*std::tan(Alp2);
01471   
01472   AllocateMemory(8,6);
01473 
01474   pV[1] = G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
01475   pV[2] = G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
01476   pV[3] = G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
01477   pV[4] = G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
01478   pV[5] = G4Point3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
01479   pV[6] = G4Point3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
01480   pV[7] = G4Point3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
01481   pV[8] = G4Point3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
01482 
01483   CreatePrism();
01484 }
01485 
01486 HepPolyhedronTrap::~HepPolyhedronTrap() {}
01487 
01488 HepPolyhedronPara::HepPolyhedronPara(G4double Dx, G4double Dy, G4double Dz,
01489                                      G4double Alpha, G4double Theta,
01490                                      G4double Phi)
01491   : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
01492 
01493 HepPolyhedronPara::~HepPolyhedronPara() {}
01494 
01495 HepPolyhedronParaboloid::HepPolyhedronParaboloid(G4double r1,
01496                                                  G4double r2,
01497                                                  G4double dz,
01498                                                  G4double sPhi,
01499                                                  G4double dPhi) 
01500 /***********************************************************************
01501  *                                                                     *
01502  * Name: HepPolyhedronParaboloid                     Date:    28.06.07 *
01503  * Author: L.Lindroos, T.Nikitina (CERN), July 2007  Revised: 28.06.07 *
01504  *                                                                     *
01505  * Function: Constructor for paraboloid                                *
01506  *                                                                     *
01507  * Input: r1    - inside and outside radiuses at -Dz                   *
01508  *        r2    - inside and outside radiuses at +Dz                   *
01509  *        dz    - half length in Z                                     *
01510  *        sPhi  - starting angle of the segment                        *
01511  *        dPhi  - segment range                                        *
01512  *                                                                     *
01513  ***********************************************************************/
01514 {
01515   static const G4double wholeCircle=twopi;
01516 
01517   //   C H E C K   I N P U T   P A R A M E T E R S
01518 
01519   G4int k = 0;
01520   if (r1 < 0. || r2 <= 0.)        k = 1;
01521 
01522   if (dz <= 0.) k += 2;
01523 
01524   G4double phi1, phi2, dphi;
01525 
01526   if(dPhi < 0.)
01527   {
01528     phi2 = sPhi; phi1 = phi2 + dPhi;
01529   }
01530   else if(dPhi == 0.) 
01531   {
01532     phi1 = sPhi; phi2 = phi1 + wholeCircle;
01533   }
01534   else
01535   {
01536     phi1 = sPhi; phi2 = phi1 + dPhi;
01537   }
01538   dphi  = phi2 - phi1;
01539 
01540   if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
01541   if (dphi > wholeCircle) k += 4; 
01542 
01543   if (k != 0) {
01544     std::cerr << "HepPolyhedronParaboloid: error in input parameters";
01545     if ((k & 1) != 0) std::cerr << " (radiuses)";
01546     if ((k & 2) != 0) std::cerr << " (half-length)";
01547     if ((k & 4) != 0) std::cerr << " (angles)";
01548     std::cerr << std::endl;
01549     std::cerr << " r1=" << r1;
01550     std::cerr << " r2=" << r2;
01551     std::cerr << " dz=" << dz << " sPhi=" << sPhi << " dPhi=" << dPhi
01552               << std::endl;
01553     return;
01554   }
01555   
01556   //   P R E P A R E   T W O   P O L Y L I N E S
01557 
01558   G4int n = GetNumberOfRotationSteps();
01559   G4double dl = (r2 - r1) / n;
01560   G4double k1 = (r2*r2 - r1*r1) / 2 / dz;
01561   G4double k2 = (r2*r2 + r1*r1) / 2;
01562 
01563   G4double *zz = new G4double[n + 2], *rr = new G4double[n + 2];
01564 
01565   zz[0] = dz;
01566   rr[0] = r2;
01567 
01568   for(G4int i = 1; i < n - 1; i++)
01569   {
01570     rr[i] = rr[i-1] - dl;
01571     zz[i] = (rr[i]*rr[i] - k2) / k1;
01572     if(rr[i] < 0)
01573     {
01574       rr[i] = 0;
01575       zz[i] = 0;
01576     }
01577   }
01578 
01579   zz[n-1] = -dz;
01580   rr[n-1] = r1;
01581 
01582   zz[n] = dz;
01583   rr[n] = 0;
01584 
01585   zz[n+1] = -dz;
01586   rr[n+1] = 0;
01587 
01588   //   R O T A T E    P O L Y L I N E S
01589 
01590   RotateAroundZ(0, phi1, dphi, n, 2, zz, rr, -1, -1); 
01591   SetReferences();
01592 
01593   delete [] zz;
01594   delete [] rr;
01595 }
01596 
01597 HepPolyhedronParaboloid::~HepPolyhedronParaboloid() {}
01598 
01599 HepPolyhedronHype::HepPolyhedronHype(G4double r1,
01600                                      G4double r2,
01601                                      G4double sqrtan1,
01602                                      G4double sqrtan2,
01603                                      G4double halfZ) 
01604 /***********************************************************************
01605  *                                                                     *
01606  * Name: HepPolyhedronHype                           Date:    14.04.08 *
01607  * Author: Tatiana Nikitina (CERN)                   Revised: 14.04.08 *
01608  *                                                                     *
01609  * Function: Constructor for Hype                                      *
01610  *                                                                     *
01611  * Input: r1       - inside radius at z=0                              *
01612  *        r2       - outside radiuses at z=0                           *
01613  *        sqrtan1  - sqr of tan of Inner Stereo Angle                  *
01614  *        sqrtan2  - sqr of tan of Outer Stereo Angle                  *
01615  *        halfZ    - half length in Z                                  *
01616  *                                                                     *
01617  ***********************************************************************/
01618 {
01619   static const G4double wholeCircle=twopi;
01620 
01621   //   C H E C K   I N P U T   P A R A M E T E R S
01622 
01623   G4int k = 0;
01624   if (r2 < 0. || r1 < 0. )        k = 1;
01625   if (r1 > r2 )                   k = 1;
01626   if (r1 == r2)                   k = 1;
01627 
01628   if (halfZ <= 0.) k += 2;
01629  
01630   if (sqrtan1<0.||sqrtan2<0.) k += 4;  
01631  
01632   if (k != 0)
01633   {
01634     std::cerr << "HepPolyhedronHype: error in input parameters";
01635     if ((k & 1) != 0) std::cerr << " (radiuses)";
01636     if ((k & 2) != 0) std::cerr << " (half-length)";
01637     if ((k & 4) != 0) std::cerr << " (angles)";
01638     std::cerr << std::endl;
01639     std::cerr << " r1=" << r1 << " r2=" << r2;
01640     std::cerr << " halfZ=" << halfZ << " sqrTan1=" << sqrtan1
01641               << " sqrTan2=" << sqrtan2
01642               << std::endl;
01643     return;
01644   }
01645   
01646   //   P R E P A R E   T W O   P O L Y L I N E S
01647 
01648   G4int n = GetNumberOfRotationSteps();
01649   G4double dz = 2.*halfZ / n;
01650   G4double k1 = r1*r1;
01651   G4double k2 = r2*r2;
01652 
01653   G4double *zz = new G4double[n+n+1], *rr = new G4double[n+n+1];
01654 
01655   zz[0] = halfZ;
01656   rr[0] = std::sqrt(sqrtan2*halfZ*halfZ+k2);
01657 
01658   for(G4int i = 1; i < n-1; i++)
01659   {
01660     zz[i] = zz[i-1] - dz;
01661     rr[i] =std::sqrt(sqrtan2*zz[i]*zz[i]+k2);
01662   }
01663 
01664   zz[n-1] = -halfZ;
01665   rr[n-1] = rr[0];
01666 
01667   zz[n] = halfZ;
01668   rr[n] =  std::sqrt(sqrtan1*halfZ*halfZ+k1);
01669 
01670   for(G4int i = n+1; i < n+n; i++)
01671   {
01672     zz[i] = zz[i-1] - dz;
01673     rr[i] =std::sqrt(sqrtan1*zz[i]*zz[i]+k1);
01674   }
01675   zz[n+n] = -halfZ;
01676   rr[n+n] = rr[n];
01677 
01678   //   R O T A T E    P O L Y L I N E S
01679 
01680   RotateAroundZ(0, 0., wholeCircle, n, n, zz, rr, -1, -1); 
01681   SetReferences();
01682 
01683   delete [] zz;
01684   delete [] rr;
01685 }
01686 
01687 HepPolyhedronHype::~HepPolyhedronHype() {}
01688 
01689 HepPolyhedronCons::HepPolyhedronCons(G4double Rmn1,
01690                                      G4double Rmx1,
01691                                      G4double Rmn2,
01692                                      G4double Rmx2, 
01693                                      G4double Dz,
01694                                      G4double Phi1,
01695                                      G4double Dphi) 
01696 /***********************************************************************
01697  *                                                                     *
01698  * Name: HepPolyhedronCons::HepPolyhedronCons        Date:    15.12.96 *
01699  * Author: E.Chernyaev (IHEP/Protvino)               Revised: 15.12.96 *
01700  *                                                                     *
01701  * Function: Constructor for CONS, TUBS, CONE, TUBE                    *
01702  *                                                                     *
01703  * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz              *
01704  *        Rmn2, Rmx2 - inside and outside radiuses at +Dz              *
01705  *        Dz         - half length in Z                                *
01706  *        Phi1       - starting angle of the segment                   *
01707  *        Dphi       - segment range                                   *
01708  *                                                                     *
01709  ***********************************************************************/
01710 {
01711   static const G4double wholeCircle=twopi;
01712 
01713   //   C H E C K   I N P U T   P A R A M E T E R S
01714 
01715   G4int k = 0;
01716   if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.)        k = 1;
01717   if (Rmn1 > Rmx1 || Rmn2 > Rmx2)                              k = 1;
01718   if (Rmn1 == Rmx1 && Rmn2 == Rmx2)                            k = 1;
01719 
01720   if (Dz <= 0.) k += 2;
01721  
01722   G4double phi1, phi2, dphi;
01723   if (Dphi < 0.) {
01724     phi2 = Phi1; phi1 = phi2 - Dphi;
01725   }else if (Dphi == 0.) {
01726     phi1 = Phi1; phi2 = phi1 + wholeCircle;
01727   }else{
01728     phi1 = Phi1; phi2 = phi1 + Dphi;
01729   }
01730   dphi  = phi2 - phi1;
01731   if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
01732   if (dphi > wholeCircle) k += 4; 
01733 
01734   if (k != 0) {
01735     std::cerr << "HepPolyhedronCone(s)/Tube(s): error in input parameters";
01736     if ((k & 1) != 0) std::cerr << " (radiuses)";
01737     if ((k & 2) != 0) std::cerr << " (half-length)";
01738     if ((k & 4) != 0) std::cerr << " (angles)";
01739     std::cerr << std::endl;
01740     std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
01741     std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
01742     std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
01743               << std::endl;
01744     return;
01745   }
01746   
01747   //   P R E P A R E   T W O   P O L Y L I N E S
01748 
01749   G4double zz[4], rr[4];
01750   zz[0] =  Dz; 
01751   zz[1] = -Dz; 
01752   zz[2] =  Dz; 
01753   zz[3] = -Dz; 
01754   rr[0] =  Rmx2;
01755   rr[1] =  Rmx1;
01756   rr[2] =  Rmn2;
01757   rr[3] =  Rmn1;
01758 
01759   //   R O T A T E    P O L Y L I N E S
01760 
01761   RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1); 
01762   SetReferences();
01763 }
01764 
01765 HepPolyhedronCons::~HepPolyhedronCons() {}
01766 
01767 HepPolyhedronCone::HepPolyhedronCone(G4double Rmn1, G4double Rmx1, 
01768                                      G4double Rmn2, G4double Rmx2,
01769                                      G4double Dz) :
01770   HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
01771 
01772 HepPolyhedronCone::~HepPolyhedronCone() {}
01773 
01774 HepPolyhedronTubs::HepPolyhedronTubs(G4double Rmin, G4double Rmax,
01775                                      G4double Dz, 
01776                                      G4double Phi1, G4double Dphi)
01777   :   HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
01778 
01779 HepPolyhedronTubs::~HepPolyhedronTubs() {}
01780 
01781 HepPolyhedronTube::HepPolyhedronTube (G4double Rmin, G4double Rmax,
01782                                       G4double Dz)
01783   : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
01784 
01785 HepPolyhedronTube::~HepPolyhedronTube () {}
01786 
01787 HepPolyhedronPgon::HepPolyhedronPgon(G4double phi,
01788                                      G4double dphi,
01789                                      G4int    npdv,
01790                                      G4int    nz,
01791                                      const G4double *z,
01792                                      const G4double *rmin,
01793                                      const G4double *rmax)
01794 /***********************************************************************
01795  *                                                                     *
01796  * Name: HepPolyhedronPgon                           Date:    09.12.96 *
01797  * Author: E.Chernyaev                               Revised:          *
01798  *                                                                     *
01799  * Function: Constructor of polyhedron for PGON, PCON                  *
01800  *                                                                     *
01801  * Input: phi  - initial phi                                           *
01802  *        dphi - delta phi                                             *
01803  *        npdv - number of steps along phi                             *
01804  *        nz   - number of z-planes (at least two)                     *
01805  *        z[]  - z coordinates of the slices                           *
01806  *        rmin[] - smaller r at the slices                             *
01807  *        rmax[] - bigger  r at the slices                             *
01808  *                                                                     *
01809  ***********************************************************************/
01810 {
01811   //   C H E C K   I N P U T   P A R A M E T E R S
01812 
01813   if (dphi <= 0. || dphi > twopi) {
01814     std::cerr
01815       << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
01816       << std::endl;
01817     return;
01818   }    
01819     
01820   if (nz < 2) {
01821     std::cerr
01822       << "HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
01823       << std::endl;
01824     return;
01825   }
01826 
01827   if (npdv < 0) {
01828     std::cerr
01829       << "HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
01830       << std::endl;
01831     return;
01832   }
01833 
01834   G4int i;
01835   for (i=0; i<nz; i++) {
01836     if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
01837       std::cerr
01838         << "HepPolyhedronPgon: error in radiuses rmin[" << i << "]="
01839         << rmin[i] << " rmax[" << i << "]=" << rmax[i]
01840         << std::endl;
01841       return;
01842     }
01843   }
01844 
01845   //   P R E P A R E   T W O   P O L Y L I N E S
01846 
01847   G4double *zz, *rr;
01848   zz = new G4double[2*nz];
01849   rr = new G4double[2*nz];
01850 
01851   if (z[0] > z[nz-1]) {
01852     for (i=0; i<nz; i++) {
01853       zz[i]    = z[i];
01854       rr[i]    = rmax[i];
01855       zz[i+nz] = z[i];
01856       rr[i+nz] = rmin[i];
01857     }
01858   }else{
01859     for (i=0; i<nz; i++) {
01860       zz[i]    = z[nz-i-1];
01861       rr[i]    = rmax[nz-i-1];
01862       zz[i+nz] = z[nz-i-1];
01863       rr[i+nz] = rmin[nz-i-1];
01864     }
01865   }
01866 
01867   //   R O T A T E    P O L Y L I N E S
01868 
01869   RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1); 
01870   SetReferences();
01871   
01872   delete [] zz;
01873   delete [] rr;
01874 }
01875 
01876 HepPolyhedronPgon::~HepPolyhedronPgon() {}
01877 
01878 HepPolyhedronPcon::HepPolyhedronPcon(G4double phi, G4double dphi, G4int nz,
01879                                      const G4double *z,
01880                                      const G4double *rmin,
01881                                      const G4double *rmax)
01882   : HepPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
01883 
01884 HepPolyhedronPcon::~HepPolyhedronPcon() {}
01885 
01886 HepPolyhedronSphere::HepPolyhedronSphere(G4double rmin, G4double rmax,
01887                                          G4double phi, G4double dphi,
01888                                          G4double the, G4double dthe)
01889 /***********************************************************************
01890  *                                                                     *
01891  * Name: HepPolyhedronSphere                         Date:    11.12.96 *
01892  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
01893  *                                                                     *
01894  * Function: Constructor of polyhedron for SPHERE                      *
01895  *                                                                     *
01896  * Input: rmin - internal radius                                       *
01897  *        rmax - external radius                                       *
01898  *        phi  - initial phi                                           *
01899  *        dphi - delta phi                                             *
01900  *        the  - initial theta                                         *
01901  *        dthe - delta theta                                           *
01902  *                                                                     *
01903  ***********************************************************************/
01904 {
01905   //   C H E C K   I N P U T   P A R A M E T E R S
01906 
01907   if (dphi <= 0. || dphi > twopi) {
01908     std::cerr
01909       << "HepPolyhedronSphere: wrong delta phi = " << dphi
01910       << std::endl;
01911     return;
01912   }    
01913 
01914   if (the < 0. || the > pi) {
01915     std::cerr
01916       << "HepPolyhedronSphere: wrong theta = " << the
01917       << std::endl;
01918     return;
01919   }    
01920   
01921   if (dthe <= 0. || dthe > pi) {
01922     std::cerr
01923       << "HepPolyhedronSphere: wrong delta theta = " << dthe
01924       << std::endl;
01925     return;
01926   }    
01927 
01928   if (the+dthe > pi) {
01929     std::cerr
01930       << "HepPolyhedronSphere: wrong theta + delta theta = "
01931       << the << " " << dthe
01932       << std::endl;
01933     return;
01934   }    
01935   
01936   if (rmin < 0. || rmin >= rmax) {
01937     std::cerr
01938       << "HepPolyhedronSphere: error in radiuses"
01939       << " rmin=" << rmin << " rmax=" << rmax
01940       << std::endl;
01941     return;
01942   }
01943 
01944   //   P R E P A R E   T W O   P O L Y L I N E S
01945 
01946   G4int nds = (GetNumberOfRotationSteps() + 1) / 2;
01947   G4int np1 = G4int(dthe*nds/pi+.5) + 1;
01948   if (np1 <= 1) np1 = 2;
01949   G4int np2 = rmin < spatialTolerance ? 1 : np1;
01950 
01951   G4double *zz, *rr;
01952   zz = new G4double[np1+np2];
01953   rr = new G4double[np1+np2];
01954 
01955   G4double a = dthe/(np1-1);
01956   G4double cosa, sina;
01957   for (G4int i=0; i<np1; i++) {
01958     cosa  = std::cos(the+i*a);
01959     sina  = std::sin(the+i*a);
01960     zz[i] = rmax*cosa;
01961     rr[i] = rmax*sina;
01962     if (np2 > 1) {
01963       zz[i+np1] = rmin*cosa;
01964       rr[i+np1] = rmin*sina;
01965     }
01966   }
01967   if (np2 == 1) {
01968     zz[np1] = 0.;
01969     rr[np1] = 0.;
01970   }
01971 
01972   //   R O T A T E    P O L Y L I N E S
01973 
01974   RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1); 
01975   SetReferences();
01976   
01977   delete [] zz;
01978   delete [] rr;
01979 }
01980 
01981 HepPolyhedronSphere::~HepPolyhedronSphere() {}
01982 
01983 HepPolyhedronTorus::HepPolyhedronTorus(G4double rmin,
01984                                        G4double rmax,
01985                                        G4double rtor,
01986                                        G4double phi,
01987                                        G4double dphi)
01988 /***********************************************************************
01989  *                                                                     *
01990  * Name: HepPolyhedronTorus                          Date:    11.12.96 *
01991  * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
01992  *                                                                     *
01993  * Function: Constructor of polyhedron for TORUS                       *
01994  *                                                                     *
01995  * Input: rmin - internal radius                                       *
01996  *        rmax - external radius                                       *
01997  *        rtor - radius of torus                                       *
01998  *        phi  - initial phi                                           *
01999  *        dphi - delta phi                                             *
02000  *                                                                     *
02001  ***********************************************************************/
02002 {
02003   //   C H E C K   I N P U T   P A R A M E T E R S
02004 
02005   if (dphi <= 0. || dphi > twopi) {
02006     std::cerr
02007       << "HepPolyhedronTorus: wrong delta phi = " << dphi
02008       << std::endl;
02009     return;
02010   }
02011 
02012   if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
02013     std::cerr
02014       << "HepPolyhedronTorus: error in radiuses"
02015       << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
02016       << std::endl;
02017     return;
02018   }
02019 
02020   //   P R E P A R E   T W O   P O L Y L I N E S
02021 
02022   G4int np1 = GetNumberOfRotationSteps();
02023   G4int np2 = rmin < spatialTolerance ? 1 : np1;
02024 
02025   G4double *zz, *rr;
02026   zz = new G4double[np1+np2];
02027   rr = new G4double[np1+np2];
02028 
02029   G4double a = twopi/np1;
02030   G4double cosa, sina;
02031   for (G4int i=0; i<np1; i++) {
02032     cosa  = std::cos(i*a);
02033     sina  = std::sin(i*a);
02034     zz[i] = rmax*cosa;
02035     rr[i] = rtor+rmax*sina;
02036     if (np2 > 1) {
02037       zz[i+np1] = rmin*cosa;
02038       rr[i+np1] = rtor+rmin*sina;
02039     }
02040   }
02041   if (np2 == 1) {
02042     zz[np1] = 0.;
02043     rr[np1] = rtor;
02044     np2 = -1;
02045   }
02046 
02047   //   R O T A T E    P O L Y L I N E S
02048 
02049   RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1); 
02050   SetReferences();
02051   
02052   delete [] zz;
02053   delete [] rr;
02054 }
02055 
02056 HepPolyhedronTorus::~HepPolyhedronTorus() {}
02057 
02058 HepPolyhedronEllipsoid::HepPolyhedronEllipsoid(G4double ax, G4double by,
02059                                                G4double cz, G4double zCut1,
02060                                                G4double zCut2)
02061 /***********************************************************************
02062  *                                                                     *
02063  * Name: HepPolyhedronEllipsoid                      Date:    25.02.05 *
02064  * Author: G.Guerrieri                               Revised:          *
02065  *                                                                     *
02066  * Function: Constructor of polyhedron for ELLIPSOID                   *
02067  *                                                                     *
02068  * Input: ax - semiaxis x                                              *
02069  *        by - semiaxis y                                              *
02070  *        cz - semiaxis z                                              *
02071  *        zCut1 - lower cut plane level (solid lies above this plane)  *
02072  *        zCut2 - upper cut plane level (solid lies below this plane)  *
02073  *                                                                     *
02074  ***********************************************************************/
02075 {
02076   //   C H E C K   I N P U T   P A R A M E T E R S
02077 
02078   if (zCut1 >= cz || zCut2 <= -cz || zCut1 > zCut2) {
02079     std::cerr << "HepPolyhedronEllipsoid: wrong zCut1 = " << zCut1
02080            << " zCut2 = " << zCut2
02081            << " for given cz = " << cz << std::endl;
02082     return;
02083   }
02084   if (cz <= 0.0) {
02085     std::cerr << "HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
02086       << std::endl;
02087     return;
02088   }
02089 
02090   G4double dthe;
02091   G4double sthe;
02092   G4int cutflag;
02093   cutflag= 0;
02094   if (zCut2 >= cz)
02095     {
02096       sthe= 0.0;
02097     }
02098   else
02099     {
02100       sthe= std::acos(zCut2/cz);
02101       cutflag++;
02102     }
02103   if (zCut1 <= -cz)
02104     {
02105       dthe= pi - sthe;
02106     }
02107   else
02108     {
02109       dthe= std::acos(zCut1/cz)-sthe;
02110       cutflag++;
02111     }
02112 
02113   //   P R E P A R E   T W O   P O L Y L I N E S
02114   //   generate sphere of radius cz first, then rescale x and y later
02115 
02116   G4int nds = (GetNumberOfRotationSteps() + 1) / 2;
02117   G4int np1 = G4int(dthe*nds/pi) + 2 + cutflag;
02118 
02119   G4double *zz, *rr;
02120   zz = new G4double[np1+1];
02121   rr = new G4double[np1+1];
02122   if (!zz || !rr)
02123     {
02124       G4Exception("HepPolyhedronEllipsoid::HepPolyhedronEllipsoid",
02125                   "greps1002", FatalException, "Out of memory");
02126     }
02127 
02128   G4double a = dthe/(np1-cutflag-1);
02129   G4double cosa, sina;
02130   G4int j=0;
02131   if (sthe > 0.0)
02132     {
02133       zz[j]= zCut2;
02134       rr[j]= 0.;
02135       j++;
02136     }
02137   for (G4int i=0; i<np1-cutflag; i++) {
02138     cosa  = std::cos(sthe+i*a);
02139     sina  = std::sin(sthe+i*a);
02140     zz[j] = cz*cosa;
02141     rr[j] = cz*sina;
02142     j++;
02143   }
02144   if (j < np1)
02145     {
02146       zz[j]= zCut1;
02147       rr[j]= 0.;
02148       j++;
02149     }
02150   if (j > np1)
02151     {
02152       std::cerr << "Logic error in HepPolyhedronEllipsoid, memory corrupted!"
02153                 << std::endl;
02154     }
02155   if (j < np1)
02156     {
02157       std::cerr << "Warning: logic error in HepPolyhedronEllipsoid."
02158                 << std::endl;
02159       np1= j;
02160     }
02161   zz[j] = 0.;
02162   rr[j] = 0.;
02163 
02164   
02165   //   R O T A T E    P O L Y L I N E S
02166 
02167   RotateAroundZ(0, 0.0, twopi, np1, 1, zz, rr, -1, 1); 
02168   SetReferences();
02169 
02170   delete [] zz;
02171   delete [] rr;
02172 
02173   // rescale x and y vertex coordinates
02174   {
02175     G4Point3D * p= pV;
02176     for (G4int i=0; i<nvert; i++, p++) {
02177       p->setX( p->x() * ax/cz );
02178       p->setY( p->y() * by/cz );
02179     }
02180   }
02181 }
02182 
02183 HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() {}
02184 
02185 HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone(G4double ax,
02186                                                          G4double ay,
02187                                                          G4double h,
02188                                                          G4double zTopCut) 
02189 /***********************************************************************
02190  *                                                                     *
02191  * Name: HepPolyhedronEllipticalCone                 Date:    8.9.2005 *
02192  * Author: D.Anninos                                 Revised: 9.9.2005 *
02193  *                                                                     *
02194  * Function: Constructor for EllipticalCone                            *
02195  *                                                                     *
02196  * Input: ax, ay     - X & Y semi axes at z = 0                        *
02197  *        h          - height of full cone                             *
02198  *        zTopCut    - Top Cut in Z Axis                               *
02199  *                                                                     *
02200  ***********************************************************************/
02201 {
02202   //   C H E C K   I N P U T   P A R A M E T E R S
02203 
02204   G4int k = 0;
02205   if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
02206 
02207   if (k != 0) {
02208     std::cerr << "HepPolyhedronCone: error in input parameters";
02209     std::cerr << std::endl;
02210     return;
02211   }
02212   
02213   //   P R E P A R E   T W O   P O L Y L I N E S
02214 
02215   zTopCut = (h >= zTopCut ? zTopCut : h);
02216 
02217   G4double *zz, *rr;
02218   zz = new G4double[4];
02219   rr = new G4double[4];
02220   zz[0] =   zTopCut; 
02221   zz[1] =  -zTopCut; 
02222   zz[2] =   zTopCut; 
02223   zz[3] =  -zTopCut; 
02224   rr[0] =  (h-zTopCut);
02225   rr[1] =  (h+zTopCut);
02226   rr[2] =  0.;
02227   rr[3] =  0.;
02228 
02229   //   R O T A T E    P O L Y L I N E S
02230 
02231   RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, -1, -1); 
02232   SetReferences();
02233 
02234   delete [] zz;
02235   delete [] rr;
02236 
02237   // rescale x and y vertex coordinates
02238  {
02239    G4Point3D * p= pV;
02240    for (G4int i=0; i<nvert; i++, p++) {
02241      p->setX( p->x() * ax );
02242      p->setY( p->y() * ay );
02243    }
02244  }
02245 }
02246 
02247 HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() {}
02248 
02249 G4int HepPolyhedron::fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
02250 /***********************************************************************
02251  *                                                                     *
02252  * Name: HepPolyhedron::fNumberOfRotationSteps       Date:    24.06.97 *
02253  * Author: J.Allison (Manchester University)         Revised:          *
02254  *                                                                     *
02255  * Function: Number of steps for whole circle                          *
02256  *                                                                     *
02257  ***********************************************************************/
02258 
02259 #include "BooleanProcessor.src"
02260 
02261 HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const 
02262 /***********************************************************************
02263  *                                                                     *
02264  * Name: HepPolyhedron::add                          Date:    19.03.00 *
02265  * Author: E.Chernyaev                               Revised:          *
02266  *                                                                     *
02267  * Function: Boolean "union" of two polyhedra                          *
02268  *                                                                     *
02269  ***********************************************************************/
02270 {
02271   G4int ierr;
02272   BooleanProcessor processor;
02273   return processor.execute(OP_UNION, *this, p,ierr);
02274 }
02275 
02276 HepPolyhedron HepPolyhedron::intersect(const HepPolyhedron & p) const 
02277 /***********************************************************************
02278  *                                                                     *
02279  * Name: HepPolyhedron::intersect                    Date:    19.03.00 *
02280  * Author: E.Chernyaev                               Revised:          *
02281  *                                                                     *
02282  * Function: Boolean "intersection" of two polyhedra                   *
02283  *                                                                     *
02284  ***********************************************************************/
02285 {
02286   G4int ierr;
02287   BooleanProcessor processor;
02288   return processor.execute(OP_INTERSECTION, *this, p,ierr);
02289 }
02290 
02291 HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const 
02292 /***********************************************************************
02293  *                                                                     *
02294  * Name: HepPolyhedron::add                          Date:    19.03.00 *
02295  * Author: E.Chernyaev                               Revised:          *
02296  *                                                                     *
02297  * Function: Boolean "subtraction" of "p" from "this"                  *
02298  *                                                                     *
02299  ***********************************************************************/
02300 {
02301   G4int ierr;
02302   BooleanProcessor processor;
02303   return processor.execute(OP_SUBTRACTION, *this, p,ierr);
02304 }
02305 
02306 //NOTE : include the code of HepPolyhedronProcessor here
02307 //       since there is no BooleanProcessor.h
02308 
02309 #undef INTERSECTION
02310 
02311 #include "HepPolyhedronProcessor.src"
02312