]>
Commit | Line | Data |
---|---|---|
5f3c1b97 | 1 | #include "TFile.h" |
2 | #include "TH1.h" | |
3 | #include "TH1D.h" | |
4 | #include "TMath.h" | |
5 | ||
6 | #include <Riostream.h> | |
7 | ||
8 | // | |
9 | // Multiply two histos of variable bin width | |
10 | // | |
11 | TH1 * MultiplyHistos(TH1 *h1, TH1 *h2) { | |
12 | ||
13 | Int_t nbins1 = h1->GetNbinsX(); | |
14 | Int_t nbins2 = h2->GetNbinsX(); | |
15 | Double_t binwidth1 = h1->GetBinWidth(1); | |
16 | Double_t binwidth2 = h2->GetBinWidth(1); | |
17 | ||
18 | if ( nbins1!=nbins2 || binwidth1!=binwidth2 ) { | |
19 | cout << "Histos do not have the same binning, do not seem compatible" << endl; | |
20 | return NULL; | |
21 | } | |
22 | ||
23 | // | |
24 | // Get the bins & limits | |
25 | Double_t *limits = new Double_t[nbins1+1]; | |
26 | Double_t xlow=0., binwidth=0.; | |
aceff457 | 27 | for (Int_t i=1; i<=nbins1; i++) { |
5f3c1b97 | 28 | binwidth = h1->GetBinWidth(i); |
29 | xlow = h1->GetBinLowEdge(i); | |
30 | limits[i-1] = xlow; | |
31 | } | |
32 | limits[nbins1] = xlow + binwidth; | |
33 | ||
34 | TH1D *hMultiply = new TH1D("hMultiply","hMultiply",nbins1,limits); | |
35 | ||
36 | Double_t value=0., err=0.; | |
aceff457 | 37 | for (Int_t ibin=1; ibin<=nbins1; ibin++) { |
5f3c1b97 | 38 | value = h1->GetBinContent(ibin) * h2->GetBinContent(ibin); |
39 | err = value * TMath::Sqrt( (h1->GetBinError(ibin)/h1->GetBinContent(ibin)) * (h1->GetBinError(ibin)/h1->GetBinContent(ibin)) + | |
40 | (h2->GetBinError(ibin)/h2->GetBinContent(ibin)) * (h2->GetBinError(ibin)/h2->GetBinContent(ibin)) ); | |
41 | hMultiply->SetBinContent(ibin,value); | |
42 | hMultiply->SetBinError(ibin,err); | |
43 | } | |
44 | ||
45 | return (TH1*)hMultiply; | |
46 | } | |
47 | ||
48 | // | |
49 | // Main function | |
50 | // | |
51 | void ComputeEfficiencyInputFromTwoSteps (const char* recolhc10d3filename="Distributions.root", | |
52 | const char* recolhc10d3histoname="RECPIDpt", | |
53 | const char* simuAcclhc10d3filename="Distributions.root", | |
54 | const char* simuAcclhc10d3histoname="MCAccpt", | |
55 | const char* simuLimAcclhc10d4filename="", | |
56 | const char* simuLimAcclhc10d4histoname="MCLimAccpt", | |
57 | const char* simuAcclhc10d4filename="", | |
58 | const char* simuAcclhc10d4histoname="MCAccpt", | |
59 | const char* outfilename="ComputeEfficiencyInputFromTwoSteps.root") | |
60 | { | |
61 | ||
62 | TFile *recolhc10d3file = new TFile(recolhc10d3filename,"read"); | |
63 | TH1D *hrecolhc10d3 = (TH1D*)recolhc10d3file->Get(recolhc10d3histoname); | |
64 | ||
65 | TFile *simuAcclhc10d3file = new TFile(simuAcclhc10d3filename,"read"); | |
66 | TH1D *hsimuAcclhc10d3 = (TH1D*)simuAcclhc10d3file->Get(simuAcclhc10d3histoname); | |
67 | ||
68 | TFile *simuLimAcclhc10d4file = new TFile(simuLimAcclhc10d4filename,"read"); | |
69 | TH1D *hsimuLimAcclhc10d4 = (TH1D*)simuLimAcclhc10d4file->Get(simuLimAcclhc10d4histoname); | |
70 | ||
71 | TFile *simuAcclhc10d4file = new TFile(simuAcclhc10d4filename,"read"); | |
72 | TH1D *hsimuAcclhc10d4 = (TH1D*)simuAcclhc10d4file->Get(simuAcclhc10d4histoname); | |
73 | ||
74 | ||
75 | TFile *out = new TFile(outfilename,"recreate"); | |
76 | TH1D *hRecoPIDCorr = (TH1D*)MultiplyHistos(hrecolhc10d3,hsimuAcclhc10d4); | |
77 | hRecoPIDCorr->SetNameTitle("hRecoPIDCorr","hRecoPIDCorr"); | |
78 | TH1D *hSimuCorr = (TH1D*)MultiplyHistos(hsimuAcclhc10d3,hsimuLimAcclhc10d4); | |
79 | hSimuCorr->SetNameTitle("hSimuCorr","hSimuCorr"); | |
80 | ||
81 | out->cd(); | |
82 | hRecoPIDCorr->Write(); | |
83 | hSimuCorr->Write(); | |
84 | out->Close(); | |
85 | ||
86 | } |