TkMSParametrization/src/MultipleScatteringX0Data.cc

0001 #include "MultipleScatteringX0Data.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include "FWCore/ParameterSet/interface/FileInPath.h"
0004
0005 #include "TH2.h"
0006 #include "TFile.h"
0007 #include "TKey.h"
0008
0009 #include <iostream>
0010 #include <string>
0011 #include <cmath>
0012
0013 using namespace std;
0014
0015 MultipleScatteringX0Data::MultipleScatteringX0Data() : theData(nullptr) {
0016   string filename = fileName();
0017   TFile theFile(filename.c_str(), "READ");
0018   if (not theFile.IsZombie()) {
0019     theData.reset(dynamic_cast<TH2F*>(theFile.GetKey("h100")->ReadObj()));
0020     theData->SetDirectory(nullptr);
0021   }
0022   if (!theData) {
0023     throw cms::Exception("Data not found")
0024         << " ** MultipleScatteringX0Data ** file: " << filename << " <-- no data found!!!";
0025   }
0026 }
0027
0028 MultipleScatteringX0Data::~MultipleScatteringX0Data() {}
0029
0030 string MultipleScatteringX0Data::fileName() {
0031   string defName = "RecoTracker/TkMSParametrization/data/MultipleScatteringX0Data.root";
0032   edm::FileInPath f(defName);
0033   return f.fullPath();
0034 }
0035
0036 int MultipleScatteringX0Data::nBinsEta() const {
0037   if (theData)
0038     return theData->GetNbinsX();
0039   else
0040     return 0;
0041 }
0042
0043 float MultipleScatteringX0Data::minEta() const {
0044   if (theData)
0045     return theData->GetXaxis()->GetXmin();
0046   else
0047     return 0;
0048 }
0049
0050 float MultipleScatteringX0Data::maxEta() const {
0051   if (theData)
0052     return theData->GetXaxis()->GetXmax();
0053   else
0054     return 0;
0055 }
0056
0057 float MultipleScatteringX0Data::sumX0atEta(float eta, float r) const {
0058   if (!theData)
0059     return 0.;
0060   eta = fabs(eta);
0061
0062   int ieta = theData->GetXaxis()->FindBin(eta);
0063   int irad = theData->GetYaxis()->FindBin(r);
0064
0065   if (irad < theData->GetNbinsY()) {
0066     return theData->GetBinContent(ieta, irad);
0067   } else {
0068     float sumX0 = 0;
0069     for (int ir = theData->GetNbinsY(); ir > 0; ir--) {
0070       float val = theData->GetBinContent(ieta, ir);
0071       if (val > sumX0)
0072         sumX0 = val;
0073     }
0074     return sumX0;
0075   }
0076 }