Calo/src/HFFibre.cc

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/SystemOfUnits.h>
0010 #include <CLHEP/Units/PhysicalConstants.h>
0011 #include <iostream>
0012 #include <sstream>
0013
0014 //#define EDM_ML_DEBUG
0015 HFFibre::Params::Params(double iFractionOfSpeedOfLightInFibre,
0016                         const HcalDDDSimConstants* hcons,
0017                         const HcalSimulationParameters* hps)
0018     : fractionOfSpeedOfLightInFibre_{iFractionOfSpeedOfLightInFibre},
0019       gParHF_{hcons->getGparHF()},
0020       rTableHF_{hcons->getRTableHF()},
0021       shortFibreLength_{hps->shortFiberLength_},
0022       longFibreLength_{hps->longFiberLength_},
0023       attenuationLength_{hps->attenuationLength_},
0024       lambdaLimits_{{static_cast<double>(hps->lambdaLimits_[0]), static_cast<double>(hps->lambdaLimits_[1])}} {}
0025
0026 HFFibre::HFFibre(const HcalDDDSimConstants* hcons, const HcalSimulationParameters* hps, edm::ParameterSet const& p)
0027     : HFFibre(Params(p.getParameter<edm::ParameterSet>("HFShower")
0028                          .getParameter<edm::ParameterSet>("HFShowerBlock")
0029                          .getParameter<double>("CFibre"),
0030                      hcons,
0031                      hps)) {}
0032
0033 HFFibre::HFFibre(Params iP)
0034     : cFibre_(CLHEP::c_light * iP.fractionOfSpeedOfLightInFibre_),
0035       gpar_(std::move(iP.gParHF_)),
0036       radius_(std::move(iP.rTableHF_)),
0037       shortFL_(std::move(iP.shortFibreLength_)),
0038       longFL_(std::move(iP.longFibreLength_)),
0039       attL_(std::move(iP.attenuationLength_)),
0040       lambLim_(iP.lambdaLimits_) {
0041   edm::LogVerbatim("HFShower") << "HFFibre:: Speed of light in fibre " << cFibre_ / (CLHEP::m / CLHEP::ns) << " m/ns";
0042   // Attenuation length
0043   nBinAtt_ = static_cast<int>(attL_.size());
0044
0045   edm::LogVerbatim("HFShower").log([&](auto& logger) {
0046     logger << "HFFibre: " << nBinAtt_ << " attL(1/cm): ";
0047     for (int it = 0; it < nBinAtt_; it++) {
0048       if (it / 10 * 10 == it) {
0049         logger << "\n";
0050       }
0051       logger << "  " << attL_[it] * CLHEP::cm;
0052     }
0053   });
0054   edm::LogVerbatim("HFShower") << "HFFibre: Limits on lambda " << lambLim_[0] << " and " << lambLim_[1];
0055   // Fibre Lengths
0056   edm::LogVerbatim("HFShower").log([&](auto& logger) {
0057     logger << "HFFibre: " << longFL_.size() << " Long Fibre Length(cm):";
0058     for (unsigned int it = 0; it < longFL_.size(); it++) {
0059       if (it / 10 * 10 == it) {
0060         logger << "\n";
0061       }
0062       logger << "  " << longFL_[it] / CLHEP::cm;
0063     }
0064   });
0065   edm::LogVerbatim("HFShower").log([&](auto& logger) {
0066     logger << "HFFibre: " << shortFL_.size() << " Short Fibre Length(cm):";
0067     for (unsigned int it = 0; it < shortFL_.size(); it++) {
0068       if (it / 10 * 10 == it) {
0069         logger << "\n";
0070       }
0071       logger << "  " << shortFL_[it] / CLHEP::cm;
0072     }
0073   });
0074   // Now geometry parameters
0075
0076   nBinR_ = static_cast<int>(radius_.size());
0077   edm::LogVerbatim("HFShower").log([&](auto& logger) {
0078     logger << "HFFibre: " << radius_.size() << " rTable(cm):";
0079     for (int i = 0; i < nBinR_; ++i) {
0080       if (i / 10 * 10 == i) {
0081         logger << "\n";
0082       }
0083       logger << "  " << radius_[i] / CLHEP::cm;
0084     }
0085   });
0086
0087   edm::LogVerbatim("HFShower").log([&](auto& logger) {
0088     logger << "HFFibre: " << gpar_.size() << " gParHF:";
0089     for (std::size_t i = 0; i < gpar_.size(); ++i) {
0090       if (i / 10 * 10 == i) {
0091         logger << "\n";
0092       }
0093       logger << "  " << gpar_[i];
0094     }
0095   });
0096 }
0097
0098 double HFFibre::attLength(double lambda) const {
0099   int i = int(nBinAtt_ * (lambda - lambLim_[0]) / (lambLim_[1] - lambLim_[0]));
0100
0101   int j = i;
0102   if (i >= nBinAtt_)
0103     j = nBinAtt_ - 1;
0104   else if (i < 0)
0105     j = 0;
0106   double att = attL_[j];
0107 #ifdef EDM_ML_DEBUG
0108   edm::LogVerbatim("HFShower") << "HFFibre::attLength for Lambda " << lambda << " index " << i << " " << j
0109                                << " Att. Length " << att;
0110 #endif
0111   return att;
0112 }
0113
0114 double HFFibre::tShift(const G4ThreeVector& point, int depth, int fromEndAbs) const {
0115   double zFibre = zShift(point, depth, fromEndAbs);
0116   double time = zFibre / cFibre_;
0117 #ifdef EDM_ML_DEBUG
0118   edm::LogVerbatim("HFShower") << "HFFibre::tShift for point " << point << " ( depth = " << depth
0119                                << ", traversed length = " << zFibre / CLHEP::cm << " cm) = " << time / CLHEP::ns
0120                                << " ns";
0121 #endif
0122   return time;
0123 }
0124
0125 double HFFibre::zShift(const G4ThreeVector& point, int depth, int fromEndAbs) const {  // point is z-local
0126
0127   double zFibre = 0;
0128   int ieta = 0;
0129   double length = 250 * CLHEP::cm;
0130   double hR = std::sqrt((point.x()) * (point.x()) + (point.y()) * (point.y()));
0131
0132   // Defines the Radius bin by radial subdivision
0133   if (fromEndAbs >= 0) {
0134     for (int i = nBinR_ - 1; i > 0; --i)
0135       if (hR < radius_[i])
0136         ieta = nBinR_ - i - 1;
0137   }
0138
0139   // Defines the full length of the fibre
0140   if (depth == 2) {
0141     if (static_cast<int>(shortFL_.size()) > ieta)
0142       length = shortFL_[ieta] + gpar_[0];
0143   } else {
0144     if (static_cast<int>(longFL_.size()) > ieta)
0145       length = longFL_[ieta];
0146   }
0147   zFibre = length;  // from beginning of abs (full length)
0148
0149   if (fromEndAbs > 0) {
0150     zFibre -= gpar_[1];  // length from end of HF
0151   } else {
0152     double zz = 0.5 * gpar_[1] + point.z();  // depth of point of photon emission (from beginning of HF)
0153     zFibre -= zz;                            // length of fiber from point of photon emission
0154   }
0155
0156 #ifdef EDM_ML_DEBUG
0157   edm::LogVerbatim("HFShower") << "HFFibre::zShift for point " << point << " (R = " << hR / CLHEP::cm
0158                                << " cm, Index = " << ieta << ", depth = " << depth
0159                                << ", Fibre Length = " << length / CLHEP::cm << " cm = " << zFibre / CLHEP::cm << " cm)";
0160 #endif
0161   return zFibre;
0162 }