Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-15 08:48:10

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HFFibre.cc
0003 // Description: Loads the table for attenuation length and calculates it
0004 ///////////////////////////////////////////////////////////////////////////////
0005 
0006 #include "SimG4CMS/Calo/interface/HFFibre.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008 
0009 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0010 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0011 #include <iostream>
0012 #include <sstream>
0013 
0014 //#define EDM_ML_DEBUG
0015 
0016 HFFibre::HFFibre(const std::string& name,
0017                  const HcalDDDSimConstants* hcons,
0018                  const HcalSimulationParameters* hps,
0019                  edm::ParameterSet const& p)
0020     : hcalConstant_(hcons), hcalsimpar_(hps) {
0021   edm::ParameterSet m_HF =
0022       (p.getParameter<edm::ParameterSet>("HFShower")).getParameter<edm::ParameterSet>("HFShowerBlock");
0023   cFibre = c_light * (m_HF.getParameter<double>("CFibre"));
0024   edm::LogVerbatim("HFShower") << "HFFibre:: Speed of light in fibre " << cFibre / (CLHEP::m / CLHEP::ns) << " m/ns";
0025   // Attenuation length
0026   attL = hcalsimpar_->attenuationLength_;
0027   nBinAtt = static_cast<int>(attL.size());
0028   std::stringstream ss1;
0029   for (int it = 0; it < nBinAtt; it++) {
0030     if (it / 10 * 10 == it) {
0031       ss1 << "\n";
0032     }
0033     ss1 << "  " << attL[it] * CLHEP::cm;
0034   }
0035   edm::LogVerbatim("HFShower") << "HFFibre: " << nBinAtt << " attL(1/cm): " << ss1.str();
0036   // Limits on Lambda
0037   std::vector<int> nvec = hcalsimpar_->lambdaLimits_;
0038   lambLim[0] = nvec[0];
0039   lambLim[1] = nvec[1];
0040   edm::LogVerbatim("HFShower") << "HFFibre: Limits on lambda " << lambLim[0] << " and " << lambLim[1];
0041   // Fibre Lengths
0042   longFL = hcalsimpar_->longFiberLength_;
0043   std::stringstream ss2;
0044   for (unsigned int it = 0; it < longFL.size(); it++) {
0045     if (it / 10 * 10 == it) {
0046       ss2 << "\n";
0047     }
0048     ss2 << "  " << longFL[it] / CLHEP::cm;
0049   }
0050   edm::LogVerbatim("HFShower") << "HFFibre: " << longFL.size() << " Long Fibre Length(cm):" << ss2.str();
0051   shortFL = hcalsimpar_->shortFiberLength_;
0052   std::stringstream ss3;
0053   for (unsigned int it = 0; it < shortFL.size(); it++) {
0054     if (it / 10 * 10 == it) {
0055       ss3 << "\n";
0056     }
0057     ss3 << "  " << shortFL[it] / CLHEP::cm;
0058   }
0059   edm::LogVerbatim("HFShower") << "HFFibre: " << shortFL.size() << " Short Fibre Length(cm):" << ss3.str();
0060 
0061   // Now geometry parameters
0062   gpar = hcalConstant_->getGparHF();
0063   radius = hcalConstant_->getRTableHF();
0064 
0065   nBinR = static_cast<int>(radius.size());
0066   std::stringstream sss;
0067   for (int i = 0; i < nBinR; ++i) {
0068     if (i / 10 * 10 == i) {
0069       sss << "\n";
0070     }
0071     sss << "  " << radius[i] / CLHEP::cm;
0072   }
0073   edm::LogVerbatim("HFShower") << "HFFibre: " << radius.size() << " rTable(cm):" << sss.str();
0074 }
0075 
0076 double HFFibre::attLength(double lambda) {
0077   int i = int(nBinAtt * (lambda - lambLim[0]) / (lambLim[1] - lambLim[0]));
0078 
0079   int j = i;
0080   if (i >= nBinAtt)
0081     j = nBinAtt - 1;
0082   else if (i < 0)
0083     j = 0;
0084   double att = attL[j];
0085 #ifdef EDM_ML_DEBUG
0086   edm::LogVerbatim("HFShower") << "HFFibre::attLength for Lambda " << lambda << " index " << i << " " << j
0087                                << " Att. Length " << att;
0088 #endif
0089   return att;
0090 }
0091 
0092 double HFFibre::tShift(const G4ThreeVector& point, int depth, int fromEndAbs) {
0093   double zFibre = zShift(point, depth, fromEndAbs);
0094   double time = zFibre / cFibre;
0095 #ifdef EDM_ML_DEBUG
0096   edm::LogVerbatim("HFShower") << "HFFibre::tShift for point " << point << " ( depth = " << depth
0097                                << ", traversed length = " << zFibre / CLHEP::cm << " cm) = " << time / CLHEP::ns
0098                                << " ns";
0099 #endif
0100   return time;
0101 }
0102 
0103 double HFFibre::zShift(const G4ThreeVector& point, int depth, int fromEndAbs) {  // point is z-local
0104 
0105   double zFibre = 0;
0106   int ieta = 0;
0107   double length = 250 * CLHEP::cm;
0108   double hR = sqrt((point.x()) * (point.x()) + (point.y()) * (point.y()));
0109 
0110   // Defines the Radius bin by radial subdivision
0111   if (fromEndAbs >= 0) {
0112     for (int i = nBinR - 1; i > 0; --i)
0113       if (hR < radius[i])
0114         ieta = nBinR - i - 1;
0115   }
0116 
0117   // Defines the full length of the fibre
0118   if (depth == 2) {
0119     if (static_cast<int>(shortFL.size()) > ieta)
0120       length = shortFL[ieta] + gpar[0];
0121   } else {
0122     if (static_cast<int>(longFL.size()) > ieta)
0123       length = longFL[ieta];
0124   }
0125   zFibre = length;  // from beginning of abs (full length)
0126 
0127   if (fromEndAbs > 0) {
0128     zFibre -= gpar[1];  // length from end of HF
0129   } else {
0130     double zz = 0.5 * gpar[1] + point.z();  // depth of point of photon emission (from beginning of HF)
0131     zFibre -= zz;                           // length of fiber from point of photon emission
0132   }
0133 
0134 #ifdef EDM_ML_DEBUG
0135   edm::LogVerbatim("HFShower") << "HFFibre::zShift for point " << point << " (R = " << hR / CLHEP::cm
0136                                << " cm, Index = " << ieta << ", depth = " << depth
0137                                << ", Fibre Length = " << length / CLHEP::cm << " cm = " << zFibre / CLHEP::cm << " cm)";
0138 #endif
0139   return zFibre;
0140 }