CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_2_7_hltpatch1/src/DetectorDescription/Core/src/Ellipsoid.cc

Go to the documentation of this file.
00001 #include "DetectorDescription/Core/src/Ellipsoid.h"
00002 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00003 #include <iostream>
00004 #include <cmath>
00005 
00006 void DDI::Ellipsoid::stream(std::ostream & os) const
00007 {
00008    os << " xSemiAxis[cm]=" << p_[0]/cm
00009       << " ySemiAxis[cm]=" << p_[1]/cm
00010       << " zSemiAxis" << p_[2]/cm
00011       << " zBottomCut" << p_[3]/cm
00012       << " zTopCut" << p_[4]/cm;
00013 }
00014 
00015 double DDI::Ellipsoid::halfVol(double dz, double maxz)  const {
00016   double volume(0.);
00017   double z(0.);
00018   double x(p_[0]), y(p_[1]);
00019   double c2 = p_[2] * p_[2];
00020   for (;z <= maxz; z=z+dz) {
00021     // what is x here? what is y?  This assumes a TRIANGLE approximation 
00022     // This assumes that I can estimate the decrease in x and y at current z
00023     // at current z, x is the projection down to the x axis.
00024     // x = a * sqrt(1-z^2 / c^2 )
00025     // unneccesary paranoia?
00026     if ( z*z / c2 >  1.0 ) {
00027       std::cout << "error!" << std::endl;
00028       exit(0);
00029     }
00030     x = p_[0] * sqrt( 1 - z*z/c2);
00031     y = p_[1] * sqrt( 1 - z*z/c2);
00032 //     if (dispcount < 100)
00033 //     std::cout << "x = " << x << " y = " << y << " z = " << z << std::endl;
00034 //     ++dispcount;
00035     volume = volume + Geom::pi() * dz * x  * y;
00036   }
00037 //   std::cout << " x = " << x;
00038 //   std::cout << " y = " << y;
00039 //   std::cout << " z = " << z;
00040 //   std::cout << " vol = " << volume << std::endl;
00041   return volume;
00042 }
00043 
00044 double DDI::Ellipsoid::volume() const { 
00045   double volume(0.);
00046   // default if both p_[3] and p_[4] are 0
00047   volume = 4./3. * Geom::pi() * p_[0] * p_[1] * p_[2];
00048   if ( p_[3] > 0.0 ) {
00049     //fail
00050     std::cout << "FAIL: p_[3] > 0.0" <<std::endl;
00051   } else if ( p_[4] < 0.0 ) {
00052     //fail
00053     std::cout << "FAIL: p_[4] <  0.0" <<std::endl;
00054   } else if ( p_[3] < 0. && p_[4] > 0. ) {
00055     volume = halfVol (p_[4]/100000., p_[4]) + halfVol (std::fabs(p_[3]/100000.), std::fabs(p_[3]));
00056   } else if ( p_[3] < 0. ) {
00057     volume = volume / 2 + halfVol(std::fabs(p_[3]/100000.), std::fabs(p_[3]));
00058   } else if ( p_[4] > 0. ) {
00059     volume = volume / 2 + halfVol (p_[4]/100000., p_[4]);
00060   }
00061   return volume; 
00062 
00063 }