CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/SimG4CMS/Calo/src/HFFibre.cc

Go to the documentation of this file.
00001 
00002 // File: HFFibre.cc
00003 // Description: Loads the table for attenuation length and calculates it
00005 
00006 #include "SimG4CMS/Calo/interface/HFFibre.h"
00007 
00008 #include "DetectorDescription/Core/interface/DDFilter.h"
00009 #include "DetectorDescription/Core/interface/DDFilteredView.h"
00010 #include "DetectorDescription/Core/interface/DDValue.h"
00011 #include "FWCore/Utilities/interface/Exception.h"
00012 
00013 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00014 #include <iostream>
00015 
00016 //#define DebugLog
00017 
00018 HFFibre::HFFibre(std::string & name, const DDCompactView & cpv, 
00019                  edm::ParameterSet const & p) {
00020 
00021   edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
00022   cFibre           = c_light*(m_HF.getParameter<double>("CFibre"));
00023   
00024   edm::LogInfo("HFShower") << "HFFibre:: Speed of light in fibre " << cFibre
00025                            << " m/ns";
00026 
00027   std::string attribute = "Volume"; 
00028   std::string value     = "HF";
00029   DDSpecificsFilter filter1;
00030   DDValue           ddv1(attribute,value,0);
00031   filter1.setCriteria(ddv1,DDSpecificsFilter::equals);
00032   DDFilteredView fv1(cpv);
00033   fv1.addFilter(filter1);
00034   bool dodet = fv1.firstChild();
00035 
00036   if (dodet) {
00037     DDsvalues_type sv(fv1.mergedSpecifics());
00038 
00039     // Attenuation length
00040     nBinAtt      = -1;
00041     attL         = getDDDArray("attl",sv,nBinAtt);
00042     edm::LogInfo("HFShower") << "HFFibre: " << nBinAtt << " attL ";
00043     for (int it=0; it<nBinAtt; it++) 
00044       edm::LogInfo("HFShower") << "HFFibre: attL[" << it << "] = " 
00045                                << attL[it]*cm << "(1/cm)";
00046 
00047     // Limits on Lambda
00048     int nb   = 2;
00049     std::vector<double>  nvec = getDDDArray("lambLim",sv,nb);
00050     lambLim[0] = static_cast<int>(nvec[0]);
00051     lambLim[1] = static_cast<int>(nvec[1]);
00052     edm::LogInfo("HFShower") << "HFFibre: Limits on lambda " << lambLim[0]
00053                              << " and " << lambLim[1];
00054 
00055     // Fibre Lengths
00056     nb       = 0;
00057     longFL   = getDDDArray("LongFL",sv,nb);
00058     edm::LogInfo("HFShower") << "HFFibre: " << nb << " Long Fibre Length";
00059     for (int it=0; it<nb; it++) 
00060       edm::LogInfo("HFShower") << "HFFibre: longFL[" << it << "] = " 
00061                                << longFL[it]/cm << " cm";
00062     nb       = 0;
00063     shortFL   = getDDDArray("ShortFL",sv,nb);
00064     edm::LogInfo("HFShower") << "HFFibre: " << nb << " Short Fibre Length";
00065     for (int it=0; it<nb; it++) 
00066       edm::LogInfo("HFShower") << "HFFibre: shortFL[" << it << "] = " 
00067                                << shortFL[it]/cm << " cm";
00068   } else {
00069     edm::LogError("HFShower") << "HFFibre: cannot get filtered "
00070                               << " view for " << attribute << " matching "
00071                               << name;
00072     throw cms::Exception("Unknown", "HFFibre")
00073       << "cannot match " << attribute << " to " << name <<"\n";
00074   }
00075 
00076   // Now geometry parameters
00077   attribute = "ReadOutName";
00078   value     = name;
00079   DDSpecificsFilter filter2;
00080   DDValue           ddv2(attribute,value,0);
00081   filter2.setCriteria(ddv2,DDSpecificsFilter::equals);
00082   DDFilteredView fv2(cpv);
00083   fv2.addFilter(filter2);
00084   dodet     = fv2.firstChild();
00085   if (dodet) {
00086     DDsvalues_type sv(fv2.mergedSpecifics());
00087 
00088     //Special Geometry parameters
00089     int nb    = -1;
00090     gpar      = getDDDArray("gparHF",sv,nb);
00091     edm::LogInfo("HFShower") << "HFFibre: " << nb <<" gpar (cm)";
00092     for (int i=0; i<nb; i++)
00093       edm::LogInfo("HFShower") << "HFFibre: gpar[" << i << "] = "
00094                                << gpar[i]/cm << " cm";
00095 
00096     nBinR     = -1;
00097     radius    = getDDDArray("rTable",sv,nBinR);
00098     edm::LogInfo("HFShower") << "HFFibre: " << nBinR <<" rTable (cm)";
00099     for (int i=0; i<nBinR; i++)
00100       edm::LogInfo("HFShower") << "HFFibre: radius[" << i << "] = "
00101                                << radius[i]/cm << " cm";
00102   } else {
00103     edm::LogError("HFShower") << "HFFibre: cannot get filtered "
00104                               << " view for " << attribute << " matching "
00105                               << name;
00106     throw cms::Exception("Unknown", "HFFibre")
00107       << "cannot match " << attribute << " to " << name <<"\n";
00108   }
00109 }
00110 
00111 HFFibre::~HFFibre() {}
00112 
00113 double HFFibre::attLength(double lambda) {
00114 
00115   int i = int(nBinAtt*(lambda - lambLim[0])/(lambLim[1]-lambLim[0]));
00116 
00117   int j =i;
00118   if (i >= nBinAtt) 
00119     j = nBinAtt-1;
00120   else if (i < 0)
00121     j = 0;
00122   double att = attL[j];
00123 #ifdef DebugLog
00124   edm::LogInfo("HFShower") << "HFFibre::attLength for Lambda " << lambda
00125                            << " index " << i  << " " << j << " Att. Length " 
00126                            << att;
00127 #endif
00128   return att;
00129 }
00130 
00131 double HFFibre::tShift(G4ThreeVector point, int depth, int fromEndAbs) {
00132 
00133   double zFibre = zShift(point, depth, fromEndAbs);
00134   double time   = zFibre/cFibre;
00135 #ifdef DebugLog
00136   edm::LogInfo("HFShower") << "HFFibre::tShift for point " << point
00137                            << " ( depth = " << depth <<", traversed length = " 
00138                            << zFibre/cm  << " cm) = " << time/ns << " ns";
00139 #endif
00140   return time;
00141 }
00142 
00143 double HFFibre::zShift(G4ThreeVector point, int depth, int fromEndAbs) {
00144 
00145   double zFibre = 0;
00146   double hR     = sqrt((point.x())*(point.x())+(point.y())*(point.y()));
00147   int    ieta   = 0;
00148   double length = 250*cm;
00149   if (fromEndAbs < 0) {
00150     zFibre      = 0.5*gpar[1] - point.z();
00151   } else {
00152     for (int i = nBinR-1; i > 0; i--)
00153       if (hR < radius[i]) ieta = nBinR - i - 1;
00154     if (depth == 2) {
00155       if ((int)(shortFL.size()) > ieta) length = shortFL[ieta];
00156     } else {
00157       if ((int)(longFL.size())  > ieta) length = longFL[ieta];
00158     }
00159     zFibre = length;
00160     if (fromEndAbs > 0) {
00161       zFibre   -= gpar[1];
00162     } else {
00163       double zz = 0.5*gpar[1] + point.z();
00164       zFibre   -= zz;
00165     }
00166     if (depth == 2) zFibre += gpar[0];
00167   }
00168 
00169 #ifdef DebugLog
00170   edm::LogInfo("HFShower") << "HFFibre::zShift for point " << point
00171                            << " (R = " << hR/cm << " cm, Index = " << ieta 
00172                            << ", depth = " << depth << ", Fibre Length = " 
00173                            << length/cm       << " cm = " << zFibre/cm  
00174                            << " cm)";
00175 #endif
00176   return zFibre;
00177 }
00178 
00179 std::vector<double> HFFibre::getDDDArray(const std::string & str, 
00180                                          const DDsvalues_type & sv, 
00181                                          int & nmin) {
00182 
00183 #ifdef DebugLog
00184   LogDebug("HFShower") << "HFFibre:getDDDArray called for " << str 
00185                        << " with nMin " << nmin;
00186 #endif
00187   DDValue value(str);
00188   if (DDfetch(&sv,value)) {
00189 #ifdef DebugLog
00190     LogDebug("HFShower") << value;
00191 #endif
00192     const std::vector<double> & fvec = value.doubles();
00193     int nval = fvec.size();
00194     if (nmin > 0) {
00195       if (nval < nmin) {
00196         edm::LogError("HFShower") << "HFFibre : # of " << str << " bins " 
00197                                   << nval << " < " << nmin << " ==> illegal";
00198         throw cms::Exception("Unknown", "HFFibre")
00199           << "nval < nmin for array " << str <<"\n";
00200       }
00201     } else {
00202       if (nval < 1 && nmin != 0) {
00203         edm::LogError("HFShower") << "HFFibre : # of " << str << " bins " 
00204                                   << nval << " < 1 ==> illegal (nmin=" 
00205                                   << nmin << ")";
00206         throw cms::Exception("Unknown", "HFFibre")
00207           << "nval < 1 for array " << str <<"\n";
00208       }
00209     }
00210     nmin = nval;
00211     return fvec;
00212   } else {
00213     if (nmin != 0) {
00214       edm::LogError("HFShower") << "HFFibre : cannot get array " << str;
00215       throw cms::Exception("Unknown", "HFFibre")
00216         << "cannot get array " << str <<"\n";
00217     } else {
00218       std::vector<double> fvec;
00219       return fvec;
00220     }
00221   }
00222 }