Core/src/Polycone.cc

0001 #include "DetectorDescription/Core/interface/Polycone.h"
0002 #include "DataFormats/Math/interface/GeantUnits.h"
0003
0004 #include <cassert>
0005 #include <cmath>
0006
0007 #include "DetectorDescription/Core/interface/DDSolidShapes.h"
0008 #include "DetectorDescription/Core/interface/Solid.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011
0012 using DDI::Polycone;
0013 using namespace geant_units::operators;
0014
0015 Polycone::Polycone(double startPhi,
0016                    double deltaPhi,
0017                    const std::vector<double>& z,
0018                    const std::vector<double>& rmin,
0019                    const std::vector<double>& rmax)
0020     : Solid(DDSolidShape::ddpolycone_rrz) {
0021   p_.emplace_back(startPhi);
0022   p_.emplace_back(deltaPhi);
0023   if ((z.size() != rmin.size()) || (z.size() != rmax.size())) {
0024     throw cms::Exception("DDException") << "Polycone(..): std::vectors z,rmin,rmax not of same length";
0025   } else {
0026     for (unsigned int i = 0; i < z.size(); ++i) {
0027       p_.emplace_back(z[i]);
0028       p_.emplace_back(rmin[i]);
0029       p_.emplace_back(rmax[i]);
0030     }
0031   }
0032 }
0033
0034 Polycone::Polycone(double startPhi, double deltaPhi, const std::vector<double>& z, const std::vector<double>& r)
0035     : Solid(DDSolidShape::ddpolycone_rz) {
0036   p_.emplace_back(startPhi);
0037   p_.emplace_back(deltaPhi);
0038   if (z.size() != r.size()) {
0039     throw cms::Exception("DDException") << "Polycone(..): std::vectors z,rmin,rmax not of same length";
0040   } else {
0041     for (unsigned int i = 0; i < z.size(); ++i) {
0042       p_.emplace_back(z[i]);
0043       p_.emplace_back(r[i]);
0044     }
0045   }
0046 }
0047
0048 double Polycone::volume() const {
0049   double result = -1.;
0050   if (shape_ == DDSolidShape::ddpolycone_rrz) {
0051     unsigned int loop = (p_.size() - 2) / 3 - 1;
0052     assert(loop > 0);
0053     double sec = 0;
0054     int i = 2;
0055     for (unsigned int j = 2; j < (loop + 2); ++j) {
0056       double dz = std::fabs(p_[i] - p_[i + 3]);
0057       double v_min = dz * 1_pi / 3. * (p_[i + 1] * p_[i + 1] + p_[i + 4] * p_[i + 4] + p_[i + 1] * p_[i + 4]);
0058       double v_max = dz * 1_pi / 3. * (p_[i + 2] * p_[i + 2] + p_[i + 5] * p_[i + 5] + p_[i + 2] * p_[i + 5]);
0059       double s = v_max - v_min;
0060       //assert(s>=0);
0061       sec += s;
0062       i += 3;
0063     }
0064     result = sec * std::fabs(p_[1]) / 2._pi;
0065   }
0066
0067   if (shape_ == DDSolidShape::ddpolycone_rz) {
0068     double volume = 0;
0069     double phiFrom = p_[0];
0070     double phiTo = p_[0] + p_[1];
0071     double slice = (std::fabs(phiFrom - phiTo)) / 2_pi;
0072     double zBegin = 0;
0073     double zEnd = 0;
0074     double rBegin = 0;
0075     double rEnd = 0;
0076     double z = 0;
0077     unsigned int i = 2;
0078
0079     while (i < (p_.size() - 2)) {
0080       zBegin = p_[i];
0081       zEnd = p_[i + 2];
0082       rBegin = p_[i + 1];
0083       rEnd = p_[i + 3];
0084       z = zBegin - zEnd;
0085
0086       /* for calculation of volume1 look at calculation of DDConsImpl of a volume1. Furthermore z can be smaller than zero. This makes sense since volumes we have to substract */
0087       double volume1 = (rEnd * rEnd + rBegin * rBegin + rBegin * rEnd) * z / 3;
0088       volume = volume + volume1;
0089       i = i + 2;
0090     }
0091
0092     /* last line (goes from last z/r value to first */
0093     i = p_.size() - 2;
0094     zBegin = p_[i];
0095     zEnd = p_[2];
0096     rBegin = p_[i + 1];
0097     rEnd = p_[3];
0098     z = zBegin - zEnd;
0099
0100     double volume2 = (rEnd * rEnd + rBegin * rBegin + rBegin * rEnd) * z / 3;
0101     volume = volume + volume2;
0102     volume = std::fabs(slice * 1_pi * volume);
0103     result = volume;
0104   }
0105   return result;
0106 }
0107
0108 void DDI::Polycone::stream(std::ostream& os) const {
0109   os << " startPhi[deg]=" << convertRadToDeg(p_[0]) << " dPhi[deg]=" << convertRadToDeg(p_[1]) << " Sizes[cm]=";
0110   for (unsigned k = 2; k < p_.size(); ++k)
0111     os << convertMmToCm(p_[k]) << " ";
0112 }