PWG3/vertexingHF/macros/ComputeEfficiencyInputFromTwoSteps.C

   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.;
  27   for (Int_t i=1; i<=nbins1; i++) {
  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.;
  37   for (Int_t ibin=1; ibin<=nbins1; ibin++) {
  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 }