PWGLF/totEt/macros/emEt/PlotGammaEfficiency.C

   1 void SetStyles(TH1 *histo,int marker, int color,char *xtitle, char *ytitle, Bool_t scale = kFALSE);
   2
   3 void PlotGammaEfficiency(Bool_t isPhos = kFALSE){
   4
   5   gStyle->SetOptTitle(0);
   6   gStyle->SetOptStat(0);
   7   gStyle->SetOptFit(0);
   8   TString det = "";
   9     if(isPhos){
  10       det = "PHOS";
  11       filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.PHOS.LHC11a10a_bis.Run139465.root";
  12     }
  13     else{
  14       det = "EMCal";
  15       filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root";
  16     }
  17
  18
  19
  20   TFile *f = TFile::Open(filename, "READ");
  21
  22   TList *l = (TList*)f->Get("out1");
  23   TH1F  *fHistGammasGenerated = l->FindObject("fHistGammasGenerated");
  24   TH1F  *fHistGammasFound = l->FindObject("fHistGammasFound");
  25   TH1F *hEfficiency = bayneseffdiv(fHistGammasFound,fHistGammasGenerated,"Efficiency");
  26   TH2F  *fHistGammasGeneratedMult = l->FindObject("fHistGammasGeneratedCent");
  27   //TH2F  *fHistGammasFoundMult = l->FindObject("fHistGammasFoundCent");
  28   //TH2F  *fHistGammasFoundOutOfAccMult = l->FindObject("fHistGammasFoundOutOfAccCent");
  29   //    TH2F *fHistGammasFoundRecoEnergyCent;
  30   //  TH2F *fHistGammasFoundOutOfAccRecoEnergyCent;
  31   TH2F  *fHistGammasFoundMult = l->FindObject("fHistGammasFoundRecoEnergyCent");
  32   TH2F  *fHistGammasFoundOutOfAccMult = l->FindObject("fHistGammasFoundOutOfAccRecoEnergyCent");
  33    fHistGammasFoundMult->Add(fHistGammasFoundOutOfAccMult);
  34   hEfficiency->GetXaxis()->SetTitle("Energy");
  35   hEfficiency->GetYaxis()->SetTitle("efficiency");
  36   hEfficiency->GetXaxis()->SetRange(1,hEfficiency->FindBin(3));
  37
  38   TCanvas *c = new TCanvas("c","c",600,400);
  39   c->SetTopMargin(0.0187166);
  40   c->SetRightMargin(0.0201342);
  41   c->SetLeftMargin(0.112416);
  42   c->SetBottomMargin(0.15508);
  43   c->SetBorderSize(0);
  44   c->SetFillColor(0);
  45   c->SetFillColor(0);
  46   c->SetBorderMode(0);
  47   c->SetFrameFillColor(0);
  48   c->SetFrameBorderMode(0);
  49
  50   hEfficiency->Draw();
  51   c->SaveAs("/tmp/GammaEfficiency.png");
  52
  53   TCanvas *c1 = new TCanvas("c1","c1",600,400);
  54   c1->SetTopMargin(0.0187166);
  55   c1->SetRightMargin(0.0201342);
  56   c1->SetLeftMargin(0.112416);
  57   c1->SetBottomMargin(0.15508);
  58   c1->SetBorderSize(0);
  59   c1->SetFillColor(0);
  60   c1->SetFillColor(0);
  61   c1->SetBorderMode(0);
  62   c1->SetFrameFillColor(0);
  63   c1->SetFrameBorderMode(0);
  64   TH1D *fHistGammasGeneratedCent[22] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL, NULL};
  65   TH1D *fHistGammasFoundCent[22] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL, NULL};
  66   TH1D *fEff[22] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL, NULL};
  67   int colors[] = {TColor::kRed,TColor::kRed, TColor::kOrange-3, TColor::kOrange-3, TColor::kGreen+3,
  68                   TColor::kGreen+3, TColor::kBlue, TColor::kBlue, TColor::kViolet, TColor::kViolet,
  69                   TColor::kMagenta+3,TColor::kRed,TColor::kRed, TColor::kOrange-3, TColor::kOrange-3, TColor::kGreen+3,
  70                   TColor::kGreen+3, TColor::kBlue, TColor::kBlue, TColor::kViolet, TColor::kViolet,
  71                   TColor::kMagenta+3};
  72   int markers[] = {20,24,21,25,22, 26,23,32,33,27, 29,20,24,21,25,22, 26,23,32,33,27, 29};
  73   TLegend *leg1 = new TLegend(0.129195,0.692513,0.29698,0.970588);
  74   leg1->SetFillStyle(0);
  75   leg1->SetFillColor(0);
  76   leg1->SetBorderSize(0);
  77   leg1->SetTextSize(0.03);
  78   leg1->SetTextSize(0.0534759);
  79   Int_t binwidth = 4;
  80   for(int i=0;i<20;i+=binwidth){
  81     Int_t maxbin = i+1+binwidth;
  82     fHistGammasGeneratedCent[i] = fHistGammasGeneratedMult->ProjectionX(Form("All%i",i),i+1,maxbin);
  83     fHistGammasFoundCent[i] = fHistGammasFoundMult->ProjectionX(Form("Reco%i",i),i+1,maxbin);
  84     fEff[i] = (TH1D*) bayneseffdiv(fHistGammasFoundCent[i],fHistGammasGeneratedCent[i],Form("eff%i",i));
  85     SetStyles(fEff[i],markers[i],colors[i],"energy","efficiency",1.0);
  86     TString legentry = Form("%i-%i",i*5,(maxbin-1)*5);
  87     legentry +="%";
  88     leg1->AddEntry(fEff[i],legentry.Data());
  89     if(i==0){
  90       fEff[i]->SetMaximum(1.5);
  91       fEff[i]->SetMinimum(0.0);
  92       fEff[i]->GetXaxis()->SetRange(1,fEff[i]->FindBin(3.0));
  93       fEff[i]->GetXaxis()->SetTitleSize(0.07);
  94       fEff[i]->GetYaxis()->SetTitleSize(0.07);
  95       fEff[i]->GetYaxis()->SetTitleOffset(0.7);
  96       fEff[i]->GetXaxis()->SetLabelSize(0.07);
  97       fEff[i]->GetYaxis()->SetLabelSize(0.07);
  98       fEff[i]->Draw();
  99     }
 100     else{ fEff[i]->Draw("same");}
 101
 102   }
 103   leg1->Draw();
 104   TString outname = "/tmp/"+det+"GammaEfficiencyCent.eps";
 105     c1->SaveAs(outname.Data());
 106 //   fHistGammasFound->Draw();
 107
 108 }
 109
 110
 111 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name)
 112 {
 113     if(!numerator){
 114       cerr<<"Error:  numerator does not exist!"<<endl;
 115       return NULL;
 116     }
 117     if(!denominator){
 118       cerr<<"Error:  denominator does not exist!"<<endl;
 119       return NULL;
 120     }
 121     TH1F* result = (TH1F*)numerator->Clone(name);
 122     Int_t nbins = numerator->GetNbinsX();
 123     for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
 124       Double_t numeratorVal = numerator->GetBinContent(ibin);
 125       Double_t denominatorVal = denominator->GetBinContent(ibin);
 126       // Check if the errors are right or the thing is scaled
 127       Double_t numeratorValErr = numerator->GetBinError(ibin);
 128       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
 129         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
 130         numeratorVal /= rescale;
 131       }
 132       Double_t denominatorValErr = denominator->GetBinError(ibin);
 133       if (!(denominatorValErr==0. || denominatorVal==0. )) {
 134         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
 135         denominatorVal /= rescale;
 136       }
 137       Double_t quotient = 0.;
 138       if (denominatorVal!=0.) {
 139         quotient = numeratorVal/denominatorVal;
 140       }
 141       Double_t quotientErr=0;
 142       quotientErr = TMath::Sqrt(
 143                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
 144                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
 145       if(quotient>1) quotientErr = 0;
 146       result->SetBinContent(ibin,quotient);
 147       result->SetBinError(ibin,quotientErr);
 148       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
 149     }
 150     return result;
 151 }
 152 void SetStyles(TH1 *histo,int marker, int color,char *xtitle, char *ytitle, Bool_t scale){
 153   histo->Sumw2();
 154   histo->SetMarkerStyle(marker);
 155   histo->SetMarkerColor(color);
 156   histo->SetLineColor(color);
 157   histo->GetXaxis()->SetTitle(xtitle);
 158   histo->GetYaxis()->SetTitle(ytitle);
 159 }