CalorimeterProperties/src/HCALProperties.cc

0001 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0002
0003 //This class header
0004 #include "FastSimulation/CalorimeterProperties/interface/HCALProperties.h"
0005 #include <cmath>
0006 #include <iostream>
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 HCALProperties::HCALProperties(const edm::ParameterSet& fastDet) : CalorimeterProperties() {
0010   edm::ParameterSet fastDetHCAL = fastDet.getParameter<edm::ParameterSet>("HadronicCalorimeterProperties");
0011   hOPi = fastDetHCAL.getParameter<double>("HCAL_PiOverE");
0012   spotFrac = fastDetHCAL.getParameter<double>("HCAL_Sampling");
0013   HCALAeff_ = fastDetHCAL.getParameter<double>("HCALAeff");
0014   HCALZeff_ = fastDetHCAL.getParameter<double>("HCALZeff");
0015   HCALrho_ = fastDetHCAL.getParameter<double>("HCALrho");
0016   HCALradiationLengthIncm_ = fastDetHCAL.getParameter<double>("HCALradiationLengthIncm");
0017   HCALradLenIngcm2_ = fastDetHCAL.getParameter<double>("HCALradLenIngcm2");
0018   HCALmoliereRadius_ = fastDetHCAL.getParameter<double>("HCALmoliereRadius");
0019   HCALcriticalEnergy_ = fastDetHCAL.getParameter<double>("HCALcriticalEnergy");
0020   HCALinteractionLength_ = fastDetHCAL.getParameter<double>("HCALinteractionLength");
0021   etatow_ = fastDetHCAL.getParameter<std::vector<double>>("HCALetatow");
0022   hcalDepthLam_ = fastDetHCAL.getParameter<std::vector<double>>("HCALDepthLam");
0023
0024   // in principle this splitting into 42 bins may change with future detectors, but let's add a protection to make sure that differences are not typos in the configuration file:
0025   if (etatow_.size() != 42)
0026     std::cout << " HCALProperties::eta2ieta - WARNING: here we expect 42 entries instead of " << etatow_.size()
0027               << "; is the change intentional?" << std::endl;
0028   // splitting of  28-th tower is taken into account (2.65-2.853-3.0)
0029   if (hcalDepthLam_.size() != etatow_.size() - 1)
0030     std::cout << " HCALProperties::eta2ieta - WARNING: the sizes of HCALetatow and HCALDepthLam should differ by 1 "
0031                  "unit! HCALDepthLam has size "
0032               << hcalDepthLam_.size() << " and HCALetatow has size " << etatow_.size() << std::endl;
0033 }
0034
0035 double HCALProperties::getHcalDepth(double eta) const {
0036   int ieta = eta2ieta(eta);
0037
0038   /*
0039   std::cout << " HCALProperties::getHcalDepth for eta = " << eta
0040         << "  returns lam.thickness = " << hcalDepthLam_[ieta] << std::endl;
0041   */
0042
0043   return hcalDepthLam_[ieta];
0044 }
0045
0046 int HCALProperties::eta2ieta(double eta) const {
0047   // binary search in the array of towers eta edges
0048   int size = etatow_.size();
0049
0050   double x = fabs(eta);
0051   int curr = size / 2;
0052   int step = size / 4;
0053   int iter;
0054   int prevdir = 0;
0055   int actudir = 0;
0056
0057   for (iter = 0; iter < size; iter++) {
0058     if (curr >= size || curr < 1)
0059       std::cout << " HCALProperties::eta2ieta - wrong current index = " << curr << " !!!" << std::endl;
0060
0061     if ((x <= etatow_[curr]) && (x > etatow_[curr - 1]))
0062       break;
0063     prevdir = actudir;
0064     if (x > etatow_[curr]) {
0065       actudir = 1;
0066     } else {
0067       actudir = -1;
0068     }
0069     if (prevdir * actudir < 0) {
0070       if (step > 1)
0071         step /= 2;
0072     }
0073     curr += actudir * step;
0074     if (curr > size)
0075       curr = size;
0076     else {
0077       if (curr < 1) {
0078         curr = 1;
0079       }
0080     }
0081
0082     /*
0083     std::cout << " HCALProperties::eta2ieta  end of iter." << iter
0084           << " curr, etatow_[curr-1], etatow_[curr] = "
0085           << curr << " " << etatow_[curr-1] << " " << etatow_[curr] << std::endl;
0086     */
0087   }
0088
0089   /*
0090   std::cout << " HCALProperties::eta2ieta  for input x = " << x
0091         << "  found index = " << curr-1
0092         << std::endl;
0093   */
0094
0095   return curr - 1;
0096 }