00001 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringX0Data.h" 00002 #include "FWCore/Utilities/interface/Exception.h" 00003 #include "FWCore/ParameterSet/interface/FileInPath.h" 00004 00005 #include "TH2.h" 00006 #include "TFile.h" 00007 #include "TKey.h" 00008 00009 #include <iostream> 00010 #include <string> 00011 #include <math.h> 00012 00013 using namespace std; 00014 00015 00016 MultipleScatteringX0Data::MultipleScatteringX0Data() 00017 : theFile(0), theData(0) 00018 { 00019 string filename = fileName(); 00020 theFile = new TFile(filename.c_str(),"READ"); 00021 if (theFile) { 00022 theData = dynamic_cast<TH2F*> (theFile->GetKey("h100")->ReadObj()); 00023 } 00024 if (!theData) { 00025 throw cms::Exception("Data not found") 00026 << " ** MultipleScatteringX0Data ** file: " 00027 << filename 00028 <<" <-- no data found!!!"; 00029 } 00030 } 00031 00032 MultipleScatteringX0Data::~MultipleScatteringX0Data() 00033 { 00034 if(theFile) { 00035 theFile->Close(); 00036 delete theFile; 00037 } 00038 } 00039 00040 string MultipleScatteringX0Data::fileName() 00041 { 00042 string defName="RecoTracker/TkMSParametrization/data/MultipleScatteringX0Data.root"; 00043 edm::FileInPath f(defName); 00044 return f.fullPath(); 00045 } 00046 00047 int MultipleScatteringX0Data::nBinsEta() const 00048 { 00049 if (theData) return theData->GetNbinsX(); 00050 else return 0; 00051 } 00052 00053 float MultipleScatteringX0Data::minEta() const 00054 { 00055 if (theData) return theData->GetXaxis()->GetXmin(); 00056 else return 0; 00057 } 00058 00059 float MultipleScatteringX0Data::maxEta() const 00060 { 00061 if (theData) return theData->GetXaxis()->GetXmax(); 00062 else return 0; 00063 } 00064 00065 float MultipleScatteringX0Data::sumX0atEta(float eta, float r) const 00066 { 00067 if(!theData) return 0.; 00068 eta = fabs(eta); 00069 00070 int ieta = theData->GetXaxis()->FindBin(eta); 00071 int irad = theData->GetYaxis()->FindBin(r); 00072 00073 if (irad < theData->GetNbinsY()) { 00074 return theData->GetBinContent(ieta,irad); 00075 } 00076 else { 00077 float sumX0 = 0; 00078 for(int ir = theData->GetNbinsY(); ir > 0; ir--) { 00079 float val = theData->GetBinContent(ieta, ir); 00080 if (val > sumX0) sumX0 = val; 00081 } 00082 return sumX0; 00083 } 00084 }