Core/src/Polyhedra.cc

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 }