CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/DetectorDescription/Core/src/Polyhedra.cc

Go to the documentation of this file.
00001 #include "DetectorDescription/Core/src/Polyhedra.h" 
00002 #include "DetectorDescription/Base/interface/DDdebug.h"
00003 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00004 #include <cmath>
00005 
00006 using DDI::Polyhedra;
00007 
00008 using std::fabs;
00009 using std::cos;
00010 using std::sin;
00011 
00012 Polyhedra::Polyhedra( int sides, double startPhi, double deltaPhi,
00013                       const std::vector<double> & z,
00014                       const std::vector<double> & rmin,
00015                       const std::vector<double> & rmax) : Solid(ddpolyhedra_rrz)              
00016 {
00017    p_.push_back(sides);
00018    p_.push_back(startPhi);
00019    p_.push_back(deltaPhi);
00020    if((z.size()!=rmin.size()) || (z.size()!=rmax.size()) )
00021    {
00022       throw cms::Exception("DDException") << "Polyhedra(..): std::vectors z,rmin,rmax not of same length";
00023    } 
00024    else
00025    {
00026       for(unsigned int i=0;i<z.size(); ++i)
00027       {
00028          p_.push_back(z[i]);
00029          p_.push_back(rmin[i]);
00030          p_.push_back(rmax[i]);
00031       }
00032    }
00033 }             
00034 
00035 
00036 Polyhedra::Polyhedra( int sides, double startPhi, double deltaPhi,
00037                       const std::vector<double> & z,
00038                       const std::vector<double> & r) : Solid(ddpolyhedra_rz)          
00039 {
00040    p_.push_back(sides);
00041    p_.push_back(startPhi);
00042    p_.push_back(deltaPhi);
00043    if(z.size()!=r.size())
00044    {
00045       throw cms::Exception("DDException") << "Polyhedra(..): std::vectors z,rmin,rmax not of same length";
00046    } 
00047    else
00048    {
00049       for(unsigned int i=0;i<z.size(); ++i)
00050       {
00051          p_.push_back(z[i]);
00052          p_.push_back(r[i]);
00053       }
00054    }
00055 }            
00056 
00057 double Polyhedra::volume() const
00058 {
00059    double volume=0;
00060    /* the following assumption is made: there are at least 3 eaqual sides if there is a complete circle (this has to be done, otherwise you can not define a polygon in a circle */
00061    
00062    /* the calculation for the volume is similar as in the case of the polycone. However, the rotation is not defined as part of a circle, but as sides in a regular polygon (specified by parameter "sides"). The sides are defined betwee startPhi and endPhi and form triangles within the circle they are defined in. First we need to determine the aread of side. let alpha |startPhi-endPhi|. the half the angle of 1 side is beta=0.5*(alph/sides). If r is the raddius of the circle in which the regular polygon is defined, the are of such a side will be 0.5*(height side)*(base side)=0.5*(cos(beta)*r)*(2*sin(beta)*r)= cos(beta)sin(beta)*r*r. r is the radius that varies if we "walk" over the boundaries of the polygon that is described by the z and r values (this yields the same integral primitive as used with the Polycone. Read Polycone documentation in code first if you do not understand this */
00063    
00064    //FIXME: rz, rrz !!
00065    if (shape()==ddpolyhedra_rrz) 
00066    {
00067       int loop = (p_.size()-3)/3 -1;
00068       double sec=0;
00069       double a = 0.5*fabs(p_[2]/rad / p_[0]);
00070       DCOUT('V',"Polyhedra::volume(), loop=" << loop << " alph[deg]=" << a/deg);
00071       int i=3;
00072       for (int j=3; j<(loop+3); ++j) 
00073       {
00074          double dz= fabs(p_[i]-p_[i+3]);
00075          DCOUT('v', "  dz[m] =" << dz/m);
00076          /*
00077           double ai,  aii;
00078           ai  = (p_[i+2]*p_[i+2] - p_[i+1]*p_[i+1]);
00079           aii = (p_[i+5]*p_[i+5] - p_[i+4]*p_[i+4]);
00080           DCOUT('v', "  rx_i[m] =" << p_[i+2]/m << " rm_i[m] =" << p_[i+1]/m);
00081           DCOUT('v', "  rx_ii[m]=" << p_[i+5]/m << " rm_ii[m]=" << p_[i+4]/m);
00082           //double s = dz/3.*(ai*bi + 0.5*(ai*bii + bi*aii) + aii*bii);
00083           double s = dz/3.*sin(a)*cos(a)*(ai + aii + 0.5*(ai+aii));
00084           */
00085          double z=dz/2.;
00086          
00087          double H1=(p_[i+2]-p_[i+1])*cos(a);
00088          double Bl1=p_[i+1]*sin(a);
00089          double Tl1=p_[i+2]*sin(a);
00090          
00091          double H2=(p_[i+5]-p_[i+4])*cos(a);
00092          double Bl2=p_[i+4]*sin(a);
00093          double Tl2=p_[i+5]*sin(a);
00094          
00095          double s = (2*H1*Bl1+2*H1*Tl1)*z+(H1*Bl2-2*H1*Bl1+H1*Tl2-2*H1*Tl1+H2*Bl1+H2*Tl1+H2*Tl2-H2*Tl1)*z+(2/3)*(H2*Bl2-H2*Bl1-H1*Bl2+H1*Bl1-H1*Tl2+H1*Tl1)*z; 
00096          s = s*p_[0];
00097          sec += s;
00098          i+=3;
00099       }
00100       volume=sec;
00101       return volume;
00102    }  
00103    int sides=int(p_[0]);
00104    //double phiFrom=p_[1]/rad;
00105    double phiDelta=p_[2]/rad;
00106    
00107    double zBegin=0;
00108    double zEnd=0;
00109    double rBegin=0;
00110    double rEnd=0;
00111    double z=0;
00112    double alpha=0;
00113    double beta=0;
00114    unsigned int i=3;
00115    
00116    alpha=fabs(phiDelta);
00117    beta=0.5*(alpha/sides);
00118    
00119    while(i<(p_.size()-2))
00120    {
00121       zBegin=p_[i];
00122       zEnd=p_[i+2];
00123       rBegin=p_[i+1];
00124       rEnd=p_[i+3];
00125       z=zBegin-zEnd;
00126       
00127       /* volume for 1 side (we multiplie by cos(beta)sin(beta)sides later*/
00128       double volume1=(rEnd*rEnd+rBegin*rBegin+rBegin*rEnd)*z/3;
00129       
00130       volume=volume+volume1;
00131       
00132       i=i+2;
00133    }
00134    
00135    /* last line (goes from last z/r value to first */
00136    
00137    i=p_.size()-2;
00138    zBegin=p_[i];
00139    zEnd=p_[3];
00140    rBegin=p_[i+1];
00141    rEnd=p_[4];
00142    z=zBegin-zEnd;
00143    
00144    double volume2=(rEnd*rEnd+rBegin*rBegin+rBegin*rEnd)*z/3;
00145    
00146    volume=volume+volume2;
00147    
00148    volume=fabs(sides*cos(beta)*sin(beta)*volume);
00149    
00150    return volume;
00151 }