PWGLF/totEt/macros/emEt/hadCorr/PlotTrackMatchFractions.C

   1 void SetStyles(TH1 *histo,int marker, int color,char *xtitle, char *ytitle);
   2 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) ;
   3 void PlotTrackMatchFractions(TString filename="rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root", TString det = "EMCal"){
   4   gStyle->SetOptTitle(0);
   5   gStyle->SetOptStat(0);
   6   gStyle->SetOptFit(0);
   7   TString outname = "/tmp/TrackMatchFractionsGammas"+det+".png";
   8   TString outname2 = "/tmp/TrackMatchFractionsHadrons"+det+".png";
   9   //TString outnamebin = Form("%iTo%i",bin,binLast);
  10
  11   TFile *f = TFile::Open(filename, "READ");
  12
  13   TList *l = (TList*)f->Get("out1");
  14   TH1F *fHistGammasCut = l->FindObject("fHistGammasTracksCut");
  15   TH1F *fHistGammasAccepted = l->FindObject("fHistGammasTracksAccepted");
  16   TH1F *fHistChargedTracksCut = l->FindObject("fHistChargedTracksCut");
  17   TH1F *fHistChargedTracksAccepted = l->FindObject("fHistChargedTracksAccepted");
  18   TH1F *fHistGammasTotal = fHistGammasCut->Clone("fHistGammasTotal");
  19   fHistGammasTotal->Add(fHistGammasAccepted);
  20   TH1F *fHistGammasFracCut = bayneseffdiv(fHistGammasCut,fHistGammasTotal,"fHistGammasFracCut");
  21   TH1F *fHistGammasFracAcc = bayneseffdiv(fHistGammasAccepted,fHistGammasTotal,"fHistGammasFracCut");
  22   SetStyles(fHistGammasFracCut,20,2,"E_{T}","Fraction");
  23   SetStyles(fHistGammasFracAcc,25,4,"E_{T}","Fraction");
  24   TH1F *fHistChargedTracksTotal = fHistChargedTracksCut->Clone("fHistChargedTracksTotal");
  25   fHistChargedTracksTotal->Add(fHistChargedTracksAccepted);
  26   TH1F *fHistChargedTracksFracCut = bayneseffdiv(fHistChargedTracksCut,fHistChargedTracksTotal,"fHistChargedTracksFracCut");
  27   TH1F *fHistChargedTracksFracAcc = bayneseffdiv(fHistChargedTracksAccepted,fHistChargedTracksTotal,"fHistChargedTracksFracCut");
  28   SetStyles(fHistChargedTracksFracCut,20,2,"E_{T}","Fraction");
  29   SetStyles(fHistChargedTracksFracAcc,25,4,"E_{T}","Fraction");
  30   TLegend *leg = new TLegend(0.357383,0.47043,0.458054,0.75);
  31   leg->SetFillStyle(0);
  32   leg->SetFillColor(0);
  33   leg->SetBorderSize(0);
  34   leg->SetTextSize(0.03);
  35   leg->AddEntry(fHistGammasFracCut,"Fraction of #gamma cut");
  36   leg->AddEntry(fHistGammasFracAcc,"Fraction of #gamma accepted");
  37   leg->SetTextSize(0.061828);
  38   TCanvas *c1 = new TCanvas("c1","c1",600,400);
  39   c1->SetTopMargin(0.02);
  40   c1->SetRightMargin(0.02);
  41   c1->SetBorderSize(0);
  42   c1->SetFillColor(0);
  43   c1->SetFillColor(0);
  44   c1->SetBorderMode(0);
  45   c1->SetFrameFillColor(0);
  46   c1->SetFrameBorderMode(0);
  47   fHistGammasFracCut->SetMaximum(1.0);
  48   fHistGammasFracCut->SetMinimum(0.0);
  49   fHistGammasFracCut->Draw();
  50   fHistGammasFracAcc->Draw("same");
  51   leg->Draw();
  52   c1->SaveAs(outname.Data());
  53
  54   TLegend *leg2 = new TLegend(0.357383,0.47043,0.458054,0.75);
  55   leg2->SetFillStyle(0);
  56   leg2->SetFillColor(0);
  57   leg2->SetBorderSize(0);
  58   leg2->SetTextSize(0.03);
  59   leg2->AddEntry(fHistGammasFracCut,"Fraction of hadrons cut");
  60   leg2->AddEntry(fHistGammasFracAcc,"Fraction of hadrons accepted");
  61   leg2->SetTextSize(0.061828);
  62   TCanvas *c2 = new TCanvas("c2","c2",600,400);
  63   c2->SetTopMargin(0.02);
  64   c2->SetRightMargin(0.02);
  65   c2->SetBorderSize(0);
  66   c2->SetFillColor(0);
  67   c2->SetFillColor(0);
  68   c2->SetBorderMode(0);
  69   c2->SetFrameFillColor(0);
  70   c2->SetFrameBorderMode(0);
  71   fHistChargedTracksFracCut->SetMaximum(1.0);
  72   fHistChargedTracksFracCut->SetMinimum(0.0);
  73   fHistChargedTracksFracCut->Draw();
  74   fHistChargedTracksFracAcc->Draw("same");
  75   leg2->Draw();
  76
  77   c2->SaveAs(outname2.Data());
  78
  79 }
  80 void SetStyles(TH1 *histo,int marker, int color,char *xtitle, char *ytitle){
  81   histo->SetMarkerStyle(marker);
  82   histo->SetMarkerColor(color);
  83   histo->SetLineColor(color);
  84   histo->GetXaxis()->SetTitle(xtitle);
  85   histo->GetYaxis()->SetTitle(ytitle);
  86 }
  87
  88
  89 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name)
  90 {
  91     if(!numerator){
  92       cerr<<"Error:  numerator does not exist!"<<endl;
  93       return NULL;
  94     }
  95     if(!denominator){
  96       cerr<<"Error:  denominator does not exist!"<<endl;
  97       return NULL;
  98     }
  99     TH1* result = (TH1*)numerator->Clone(name);
 100     Int_t nbins = numerator->GetNbinsX();
 101     for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
 102       Double_t numeratorVal = numerator->GetBinContent(ibin);
 103       Double_t denominatorVal = denominator->GetBinContent(ibin);
 104       // Check if the errors are right or the thing is scaled
 105       Double_t numeratorValErr = numerator->GetBinError(ibin);
 106       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
 107         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
 108         numeratorVal /= rescale;
 109       }
 110       Double_t denominatorValErr = denominator->GetBinError(ibin);
 111       if (!(denominatorValErr==0. || denominatorVal==0. )) {
 112         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
 113         denominatorVal /= rescale;
 114       }
 115       Double_t quotient = 0.;
 116       if (denominatorVal!=0.) {
 117         quotient = numeratorVal/denominatorVal;
 118       }
 119       Double_t quotientErr=0;
 120       quotientErr = TMath::Sqrt(
 121                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
 122                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
 123       result->SetBinContent(ibin,quotient);
 124       result->SetBinError(ibin,quotientErr);
 125       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
 126     }
 127     return result;
 128 }