Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-04 02:02:27

0001 #include "DetectorDescription/Core/interface/Polyhedra.h"
0002 #include "DataFormats/Math/interface/GeantUnits.h"
0003 
0004 #include <cmath>
0005 
0006 #include "DetectorDescription/Core/interface/DDSolidShapes.h"
0007 #include "DetectorDescription/Core/interface/Solid.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010 
0011 using DDI::Polyhedra;
0012 
0013 using std::cos;
0014 using std::fabs;
0015 using std::sin;
0016 using namespace geant_units;
0017 using namespace geant_units::operators;
0018 
0019 Polyhedra::Polyhedra(int sides,
0020                      double startPhi,
0021                      double deltaPhi,
0022                      const std::vector<double>& z,
0023                      const std::vector<double>& rmin,
0024                      const std::vector<double>& rmax)
0025     : Solid(DDSolidShape::ddpolyhedra_rrz) {
0026   p_.emplace_back(sides);
0027   p_.emplace_back(startPhi);
0028   p_.emplace_back(deltaPhi);
0029   if ((z.size() != rmin.size()) || (z.size() != rmax.size())) {
0030     throw cms::Exception("DDException") << "Polyhedra(..): std::vectors z,rmin,rmax not of same length";
0031   } else {
0032     for (unsigned int i = 0; i < z.size(); ++i) {
0033       p_.emplace_back(z[i]);
0034       p_.emplace_back(rmin[i]);
0035       p_.emplace_back(rmax[i]);
0036     }
0037   }
0038 }
0039 
0040 Polyhedra::Polyhedra(
0041     int sides, double startPhi, double deltaPhi, const std::vector<double>& z, const std::vector<double>& r)
0042     : Solid(DDSolidShape::ddpolyhedra_rz) {
0043   p_.emplace_back(sides);
0044   p_.emplace_back(startPhi);
0045   p_.emplace_back(deltaPhi);
0046   if (z.size() != r.size()) {
0047     throw cms::Exception("DDException") << "Polyhedra(..): std::vectors z,rmin,rmax not of same length";
0048   } else {
0049     for (unsigned int i = 0; i < z.size(); ++i) {
0050       p_.emplace_back(z[i]);
0051       p_.emplace_back(r[i]);
0052     }
0053   }
0054 }
0055 
0056 double Polyhedra::volume() const {
0057   double volume = 0;
0058   /* the following assumption is made:
0059      there are at least 3 eaqual sides
0060      if there is a complete circle (this has to be done,
0061      otherwise you can not define a polygon in a circle */
0062 
0063   /* the calculation for the volume is similar as in
0064      the case of the polycone. However, the rotation
0065      is not defined as part of a circle, but as sides
0066      in a regular polygon (specified by parameter "sides").
0067      The sides are defined betwee startPhi and endPhi and
0068      form triangles within the circle they are defined in.
0069      First we need to determine the aread of side.
0070      let alpha |startPhi-endPhi|.
0071      the half the angle of 1 side is beta=0.5*(alph/sides).
0072      If r is the raddius of the circle in which the regular polygon is defined,
0073      the are of such a side will be
0074      0.5*(height side)*(base side)=0.5*(cos(beta)*r)*(2*sin(beta)*r)= cos(beta)sin(beta)*r*r.
0075      r is the radius that varies if we "walk" over the boundaries of
0076      the polygon that is described by the z and r values
0077      (this yields the same integral primitive as used with the Polycone.
0078      Read Polycone documentation in code first if you do not understand this */
0079 
0080   //FIXME: rz, rrz !!
0081   if (shape() == DDSolidShape::ddpolyhedra_rrz) {
0082     int loop = (p_.size() - 3) / 3 - 1;
0083     double sec = 0;
0084     double a = 0.5 * fabs(p_[2] / p_[0]);
0085     int i = 3;
0086     for (int j = 3; j < (loop + 3); ++j) {
0087       double dz = fabs(p_[i] - p_[i + 3]);
0088       /*
0089     double ai,  aii;
0090     ai  = (p_[i+2]*p_[i+2] - p_[i+1]*p_[i+1]);
0091     aii = (p_[i+5]*p_[i+5] - p_[i+4]*p_[i+4]);
0092     //double s = dz/3.*(ai*bi + 0.5*(ai*bii + bi*aii) + aii*bii);
0093     double s = dz/3.*sin(a)*cos(a)*(ai + aii + 0.5*(ai+aii));
0094       */
0095       double z = dz / 2.;
0096 
0097       double H1 = (p_[i + 2] - p_[i + 1]) * cos(a);
0098       double Bl1 = p_[i + 1] * sin(a);
0099       double Tl1 = p_[i + 2] * sin(a);
0100 
0101       double H2 = (p_[i + 5] - p_[i + 4]) * cos(a);
0102       double Bl2 = p_[i + 4] * sin(a);
0103       double Tl2 = p_[i + 5] * sin(a);
0104 
0105       double s = (2 * H1 * Bl1 + 2 * H1 * Tl1) * z +
0106                  (H1 * Bl2 - 2 * H1 * Bl1 + H1 * Tl2 - 2 * H1 * Tl1 + H2 * Bl1 + H2 * Tl1 + H2 * Tl2 - H2 * Tl1) * z +
0107                  (2 / 3) * (H2 * Bl2 - H2 * Bl1 - H1 * Bl2 + H1 * Bl1 - H1 * Tl2 + H1 * Tl1) * z;
0108       s = s * p_[0];
0109       sec += s;
0110       i += 3;
0111     }
0112     volume = sec;
0113     return volume;
0114   }
0115   int sides = int(p_[0]);
0116   //double phiFrom=p_[1]/rad;
0117   double phiDelta = p_[2];
0118 
0119   double zBegin = 0;
0120   double zEnd = 0;
0121   double rBegin = 0;
0122   double rEnd = 0;
0123   double z = 0;
0124   double alpha = 0;
0125   double beta = 0;
0126   unsigned int i = 3;
0127 
0128   alpha = fabs(phiDelta);
0129   beta = 0.5 * (alpha / sides);
0130 
0131   while (i < (p_.size() - 2)) {
0132     zBegin = p_[i];
0133     zEnd = p_[i + 2];
0134     rBegin = p_[i + 1];
0135     rEnd = p_[i + 3];
0136     z = zBegin - zEnd;
0137 
0138     /* volume for 1 side (we multiplie by cos(beta)sin(beta)sides later*/
0139     double volume1 = (rEnd * rEnd + rBegin * rBegin + rBegin * rEnd) * z / 3;
0140 
0141     volume = volume + volume1;
0142 
0143     i = i + 2;
0144   }
0145 
0146   /* last line (goes from last z/r value to first */
0147 
0148   i = p_.size() - 2;
0149   zBegin = p_[i];
0150   zEnd = p_[3];
0151   rBegin = p_[i + 1];
0152   rEnd = p_[4];
0153   z = zBegin - zEnd;
0154 
0155   double volume2 = (rEnd * rEnd + rBegin * rBegin + rBegin * rEnd) * z / 3;
0156 
0157   volume = volume + volume2;
0158 
0159   volume = fabs(sides * cos(beta) * sin(beta) * volume);
0160 
0161   return volume;
0162 }
0163 
0164 void DDI::Polyhedra::stream(std::ostream& os) const {
0165   os << " sides=" << p_[0] << " startPhi[deg]=" << convertRadToDeg(p_[1]) << " dPhi[deg]=" << convertRadToDeg(p_[2])
0166      << " Sizes[cm]=";
0167   for (unsigned k = 3; k < p_.size(); ++k)
0168     os << convertMmToCm(p_[k]) << " ";
0169 }