Core/src/Cons.cc

0001 #include "DetectorDescription/Core/interface/Cons.h"
0002 #include "DataFormats/Math/interface/GeantUnits.h"
0003
0004 #include <cmath>
0005 #include <ostream>
0006 #include <vector>
0007
0008 #include "DetectorDescription/Core/interface/DDSolidShapes.h"
0009 #include "DetectorDescription/Core/interface/Solid.h"
0010
0011 using namespace geant_units;
0012 using namespace geant_units::operators;
0013
0014 DDI::Cons::Cons(double zhalf,
0015                 double rInMinusZ,
0016                 double rOutMinusZ,
0017                 double rInPlusZ,
0018                 double rOutPlusZ,
0019                 double startPhi,
0020                 double deltaPhi)
0021     : Solid(DDSolidShape::ddcons) {
0022   p_.emplace_back(zhalf);
0023   p_.emplace_back(rInMinusZ);
0024   p_.emplace_back(rOutMinusZ);
0025   p_.emplace_back(rInPlusZ);
0026   p_.emplace_back(rOutPlusZ);
0027   p_.emplace_back(startPhi);
0028   p_.emplace_back(deltaPhi);
0029 }
0030
0031 void DDI::Cons::stream(std::ostream& os) const {
0032   os << " zhalf=" << convertMmToCm(p_[0]) << " rIn-Z=" << convertMmToCm(p_[1]) << " rOut-Z=" << convertMmToCm(p_[2])
0033      << " rIn+Z=" << convertMmToCm(p_[3]) << " rOut+Z=" << convertMmToCm(p_[4])
0034      << " startPhi=" << convertRadToDeg(p_[5]) << " deltaPhi=" << convertRadToDeg(p_[6]);
0035 }
0036
0037 double DDI::Cons::volume() const {
0038   /* zhalf is the half length of the cone,
0039      phiTo is always clockwise rotated from phiFrom
0040      rInMinusZ is always smaller than rOutMinusZ (same for r*PlusZ)
0041      They are the distances relative to the rotation axes */
0042
0043   /* calculation normalize from 0 to z */
0044
0045   /* The function f=rInMinusZ+((rInPlusZ-rInMinusZ)/z)*x defines
0046      the radius of the the rotation from 0 to z. Raised to the power
0047      of 2 integrated on x from 0 to z. Multiplied by pi, gives the
0048      volume that needs to substracted from the other volume */
0049
0050   /* f^2=rInMinusZ*rInMinusZ+2*rInMinusZ*((rInPlusZ-rInMinusZ)/z)*x+((rInPlusZ-rInMinusZ)*(rInPlusZ-rInMinusZ)*x*x)/(z*z) */
0051
0052   /* primitive of f^2 is: rInMinusZ*rInMinusZ*x+rInMinusZ*((rInPlusZ-rInMinusZ)/z)*(x*x)+(rInPlusZ-rInMinusZ)*(rInPlusZ-rInMinusZ)*(x*x*x)/(3*z*z) */
0053
0054   /*integration from 0 to z yields: pi*( rInMinusZ*rInMinusZ*z+rInMinusZ*(rInPlusZ-rInMinusZ)*z+((rInPlusZ-rInMinusZ)*(rInPlusZ-rInMinusZ)*z)/(3) ) */
0055
0056   double zhalf = p_[0];
0057   double rInMinusZ = p_[1];
0058   double rOutMinusZ = p_[2];
0059   double rInPlusZ = p_[3];
0060   double rOutPlusZ = p_[4];
0061   //double phiFrom=p_[5]/rad;
0062   double deltaPhi = std::abs(p_[6]);
0063   double z = 2 * zhalf;
0064
0065   double volume1 = 1_pi * (rInPlusZ * rInPlusZ + rInMinusZ * rInMinusZ + rInMinusZ * rInPlusZ) * z / 3;
0066
0067   double volume2 = 1_pi * (rOutPlusZ * rOutPlusZ + rOutMinusZ * rOutMinusZ + rOutMinusZ * rOutPlusZ) * z / 3;
0068
0069   double slice = deltaPhi / (2_pi);
0070   double volume = slice * (volume2 - volume1);
0071
0072   return volume;
0073 }