Reference Manual

Go to the documentation of this file.
00001 #include "DetectorDescription/Core/src/Trap.h"
00002 #include "DetectorDescription/Base/interface/DDdebug.h"
00003 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00004 #include <cmath>
00005 
00006 using std::sqrt;
00007 
00008 
00009 DDI::Trap::Trap( double pDz, 
00010                  double pTheta,
00011                  double pPhi,
00012                  double pDy1, double pDx1,double pDx2,
00013                  double pAlp1,
00014                  double pDy2, double pDx3, double pDx4,
00015                  double pAlp2 )
00016  : Solid(ddtrap) 
00017 {                
00018   p_.push_back(pDz); // ......... 0
00019   p_.push_back(pTheta); // .. 1
00020   p_.push_back(pPhi); // ....... 2
00021   p_.push_back(pDy1); // ........ 3
00022   p_.push_back(pDx1); // ........ 4
00023   p_.push_back(pDx2); // ........ 5
00024   p_.push_back(pAlp1); // ....... 6
00025   p_.push_back(pDy2); // ........ 7
00026   p_.push_back(pDx3); // ......... 8
00027   p_.push_back(pDx4); // ........ 9
00028   p_.push_back(pAlp2);
00029 }
00030 
00031 
00032 void DDI::Trap::stream(std::ostream & os) const
00033 {
00034   os << " dz=" << p_[0]/cm
00035      << " theta=" << p_[1]/deg
00036      << " phi=" << p_[2]/deg
00037      << " dy1=" << p_[3]/cm
00038      << " dx1=" << p_[4]/cm
00039      << " dx2=" << p_[5]/cm
00040      << " alpha1=" << p_[6]/deg
00041      << " dy2=" << p_[7]/cm
00042      << " dx3=" << p_[8]/cm
00043      << " dx4=" << p_[9]/cm
00044      << " alpha2=" << p_[10]/deg;
00045 }
00046 
00047 double DDI::Trap::volume() const
00048 {
00049  double volume=0;
00050 
00051   /* use notation as described in documentation about geant 3 shapes */
00052   /* we do not need all the parameters.*/
00053 
00054   double Dz=p_[0];
00055   double H1=p_[3];
00056   double Bl1=p_[4];
00057   double Tl1=p_[5];
00058   double H2=p_[7];
00059   double Bl2=p_[8];
00060   double Tl2=p_[9];
00061 
00062   double z=2*Dz;
00063 
00064   /* the area of a trapezoid with one side of length 2*Bl1 and other side 2*Tl1,     and height 2*H1 is 0.5*(2*Bl1+2*Tl1)*2H1=2*H1(Bl1+Tl1) */
00065 
00066   /* the volume of a geometry defined by 2 2D parallel trapezoids is (in this case the integral over the area of a trapezoid that is defined as function x between these two trapezoids */
00067 
00068   /* the following formula describes this parmeterized area in x. x=0: trapezoid defined by H1,Bl1,Tl1,  x=z: trapezoid defined by H2,Bl2,Tl2 */
00069 
00070   /* area(x)=2*(H1+x/z*(H2-H1))*(Bl1+x/z*(Bl2-Bl1)+Tl1+x/z*(Tl2-Tl1)) */
00071  
00072  /* primitive of area(x):
00073     (2*H1*Bl1+2*H1*Tl1)*x+(H1*Bl2-2*H1*Bl1+H1*Tl2-2*H1*Tl1+H2*Bl1+H2*Tl1+H2*Tl2-H2*Tl1)*x*x/z+(2/3)*(H2*Bl2-H2*Bl1-H1*Bl2+H1*Bl1-H1*Tl2+H1*Tl1)*x*x*x/(z*z)   */
00074 
00075 // volume=(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*z+(2/3)*(H2*Bl2-H2*Bl1-H1*Bl2+H1*Bl1-H1*Tl2+H1*Tl1)*z*z*z; 
00076   volume=(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; 
00077 
00078 
00079  /* 
00080     Alternative:
00081     A ... height of bottom trapez, B ... middle line perpendicular to A
00082     a ... height of top trapez,    b ... middle line perpendicular to a
00083     H ... heigt of the solid
00084     
00085     V = H/3. * ( A*B + 0.5 * ( A*b + B*a ) + a*b ) <-- this is wrong ..
00086     V = H/3 * ( A*B + sqrt( A*B*a*b ) + a*b )
00087  */ 
00088   double A = 2.*p_[3];
00089   double B = p_[4] + p_[5];
00090   double a = 2.*p_[7];
00091   double b = p_[8] + p_[9];
00092   
00093 
00094   double volu_alt = 2.*p_[0]/3. * ( A*B + sqrt ( A*b*B*a ) + a*b );
00095   DCOUT('V', "alternative-volume=" << volu_alt/m3 << std::endl);
00096   
00097   //DCOUT_V('C',"DC: solid=" << this->ddname() << " vol=" << volume << " vol_a=" << volu_alt << " d=" << (volume-volu_alt)/volume*100. << "%");
00098   return volume;
00099 }
/data/refman/pasoursint/CMSSW_5_3_3/src/DetectorDescription/Core/src/Trap.cc
