Calorimetry/src/KKCorrectionFactors.cc

0001 #include "FastSimulation/Calorimetry/interface/KKCorrectionFactors.h"
0002
0003 KKCorrectionFactors::KKCorrectionFactors(const edm::ParameterSet& pset) {
0004   interpolate3D_ = pset.exists("interpolate3D") && pset.getUntrackedParameter<bool>("interpolate3D");
0005
0006   // Get the filename and histogram name
0007   std::string fileName = pset.getUntrackedParameter<std::string>("fileName");
0008   std::string histogramName = pset.getUntrackedParameter<std::string>("histogramName");
0009
0010   // Read histo
0011   edm::FileInPath myDataFile(fileName);
0012   TFile* myFile = TFile::Open(myDataFile.fullPath().c_str());
0013
0014   gROOT->cd();  // create histogram in memory
0015   auto obj = myFile->Get(histogramName.c_str());
0016
0017   // Check if histogram exists in file
0018   if (!obj) {
0019     throw cms::Exception("FileReadError")
0020         << "Histogram \"" << histogramName << "\" not found in file \"" << fileName << "\".\n";
0021   }
0022   h3_ = new TH3F(*((TH3F*)obj));
0023
0024   delete myFile;
0025 }
0026
0027 float KKCorrectionFactors::getScale(float genE, float genEta, float simE) const {
0028   float r = simE / genE;
0029   float scale = 1.;
0030
0031   if (interpolate3D_
0032       // TH3::Interpolate can only interpolate inside the bondaries of the histogram
0033       && genE > h3_->GetXaxis()->GetXmin() && genE < h3_->GetXaxis()->GetXmax() &&
0034       genEta > h3_->GetYaxis()->GetXmin() && genEta < h3_->GetYaxis()->GetXmax() && r < h3_->GetZaxis()->GetXmax() &&
0035       r > h3_->GetZaxis()->GetXmax()) {
0036     scale = h3_->Interpolate(genE, genEta, r);
0037
0038   } else {  // interpolation in r is mandatory
0039
0040     int binE = h3_->GetXaxis()->FindFixBin(genE);
0041     int binEta = h3_->GetYaxis()->FindFixBin(genEta);
0042
0043     // find the two bins which are closest to the actual value
0044     auto binWidthR = h3_->GetZaxis()->GetBinWidth(0);
0045     int binRup = h3_->GetZaxis()->FindFixBin(r + binWidthR / 2.);
0046     int binRdn = h3_->GetZaxis()->FindFixBin(r - binWidthR / 2.);
0047
0048     auto scaleUp = h3_->GetBinContent(binE, binEta, binRup);
0049     auto scaleDn = h3_->GetBinContent(binE, binEta, binRdn);
0050
0051     // make a linear extrapolation between neighbour bins if they are not zero
0052     auto Rdn = h3_->GetZaxis()->GetBinCenter(binRdn);
0053     scale = scaleDn + (scaleUp - scaleDn) * (r - Rdn) / binWidthR;
0054   }
0055   return scale;
0056 }