PWGLF/totEt/macros/emEt/EMCal/CalcEfficHadrons.C

   1 /*
   2  *  CalcEffic.C
   3  */
   4
   5 TEfficiency* CalcEfficHadrons(const char* fileName, Bool_t makeDraw=kFALSE)
   6 {
   7         // open the input file
   8         TFile *file = new TFile(fileName);
   9
  10         // ********************************************
  11         // parameters to check before running the macro
  12         // ********************************************
  13
  14         const Int_t NHISTS = 4; // Check the number of histograms for different particle species
  15         const Int_t NOUTPUTS = 3;
  16         const Int_t NHISTOUT[NOUTPUTS] = {1,1,1};
  17         const Int_t IHISTOUT[NOUTPUTS][NHISTS] = {{0,-1,-1,-1},{1,-1,-1,-1},{2,-1,-1,-1}};
  18         const Float_t CUT_RES = 0.02;
  19
  20         Int_t style[NOUTPUTS] = {20,21,22};
  21         Int_t color[NOUTPUTS] = {1,2,4};
  22
  23         const Int_t fgNumOfPtBins = 111; // Check the number of eta bins in the histograms
  24         const Int_t fgNumOfEtaBins = 16; // Check the number of E bins in the histograms
  25         const Int_t fgNumOfRBins = 45;
  26
  27         Double_t fgPtAxis[117]= {0.0,0.01,0.02,0.03,0.04, 0.05, 0.06,0.07,0.08,0.09, 0.10,0.11, .12,0.13, .14,0.15, .16,0.17, .18,0.19,
  28                 0.2, .22, .24, .26, .28, 0.30, 0.32, .34, .36, .38, 0.40, .42, .44, .46, .48,
  29                 0.5, .52, .54, .56, .58, 0.60, 0.62, .64, .66, .68, 0.70, .72, .74, .76, .78,
  30                 .80, .82, .84, .86, .88, 0.90, 0.92, .94, .96, .98, 1.00,1.05, 1.1,1.15, 1.2,
  31                 1.25, 1.3,1.35,1.40,1.45, 1.50, 1.55, 1.6,1.65, 1.7, 1.75, 1.8,1.85, 1.9,1.95,
  32                 2.0, 2.2, 2.4, 2.6, 2.8, 3.00, 3.20, 3.4, 3.6, 3.8, 4.00, 4.2, 4.4, 4.6, 4.8,
  33                 5.0, 5.5, 6.0, 6.5, 7.0, 7.50, 8.00, 8.5, 9.0, 9.5, 10.0,12.0,14.0,16.0,18.0,
  34                 20.0,25.0,30.0,35.0,40.0, 45.0, 50.0};
  35
  36         // declare histograms and graphs
  37         TH2F *histNum[NHISTS];
  38         TH2F *histDen[NHISTS];
  39         TGraphErrors *graph[NOUTPUTS];
  40         TH1D* projYNum;
  41         TEfficiency *effic[NOUTPUTS];
  42         char efficName[50];
  43
  44         // retrieve the input list of histogram. Check the TList name in the input file.
  45         TList *list = (TList*) file->Get("out1");
  46
  47         // retrieve the histograms in the list. Check the name of the histograms
  48         histNum[0] = (TH2F*)list->FindObject("fHistPionRec_ResPt_EmcalMC");
  49         histNum[1] = (TH2F*)list->FindObject("fHistKaonRec_ResPt_EmcalMC");
  50         histNum[2] = (TH2F*)list->FindObject("fHistProtonRec_ResPt_EmcalMC");
  51         histNum[3] = (TH2F*)list->FindObject("fHistMuonRec_ResPt_EmcalMC");
  52
  53         // retrieve the histograms in the list. Check the name of the histograms
  54         histDen[0] = (TH2F*)list->FindObject("fHistPionAcc_EtaPt_EmcalMC");
  55         histDen[1] = (TH2F*)list->FindObject("fHistKaonAcc_EtaPt_EmcalMC");
  56         histDen[2] = (TH2F*)list->FindObject("fHistProtonAcc_EtaPt_EmcalMC");
  57         histDen[3] = (TH2F*)list->FindObject("fHistMuonAcc_EtaPt_EmcalMC");
  58
  59         // ********************************************
  60
  61         Float_t x[fgNumOfPtBins]={0}, ex[fgNumOfPtBins]={0};
  62         Float_t y[fgNumOfPtBins]={0}, ey[fgNumOfPtBins]={0};
  63         Float_t num=0, den=0;
  64         //Int_t num=0, den=0;
  65         Float_t Res=0;
  66
  67         // loop over different desired outputs
  68         for (int iOut=0; iOut<NOUTPUTS; iOut++)
  69         {
  70                 sprintf(efficName,"effic_%d",iOut);
  71                 effic[iOut] = new TEfficiency(efficName,efficName,fgNumOfPtBins,fgPtAxis);
  72
  73                 // loop over E bins
  74                 for (int ix=0; ix<fgNumOfPtBins; ix++)
  75                 {
  76                         // initialize ET variables for a new particle species
  77                         x[ix]=histNum[0]->GetXaxis()->GetBinCenter(ix+1);
  78                         y[ix]=0;
  79                         ex[ix]=0;
  80                         ey[ix]=0;
  81                         num = 0;
  82                         den = 0;
  83
  84                         // loop over eta bins
  85                         for (int iy=0; iy<fgNumOfEtaBins; iy++)
  86                         {
  87                                 for (int iHist=0; iHist<NHISTOUT[iOut]; iHist++)
  88                                 {
  89                                         den += histDen[IHISTOUT[iOut][iHist]]->GetBinContent(ix+1,iy+1); // sum over all E bins in order to get total ET
  90                                 }
  91                         }
  92
  93                         // loop over residual bins
  94                         for (int iHist=0; iHist<NHISTOUT[iOut]; iHist++)
  95                         {
  96                                 projYNum = histNum[IHISTOUT[iOut][iHist]]->ProjectionY();
  97                                 for (int iy=0; iy<fgNumOfRBins; iy++)
  98                                 {
  99                                         Res = projYNum->GetBinCenter(iy+1);
 100                                         if (Res<CUT_RES)
 101                                                 num += histNum[IHISTOUT[iOut][iHist]]->GetBinContent(ix+1,iy+1); // sum over all E bins in order to get total ET
 102                                 }
 103                         }
 104
 105                         if ((num>0) && (den>0))
 106                         {
 107                                 effic[iOut]->SetTotalEvents(ix,den);
 108                                 effic[iOut]->SetPassedEvents(ix,num);
 109                                 y[ix] = num/den;
 110                                 ey[ix] = y[ix]*sqrt(1/num+1/den);
 111                                 //ey[ix] = ((num+1)*(num+2))/((den+2)*(den+3))-((num+1)*(num+1))/((den+2)*(den+2));
 112                         }
 113                         else
 114                         {
 115                                 y[ix] = 0;
 116                                 ey[ix] = 0;
 117                         }
 118
 119                 } // end of loop over E bins
 120
 121                 graph[iOut] = new TGraphErrors(fgNumOfPtBins,x,y,ex,ey); // graphic of ET(>E_cut)/ET(total) for a given particle species and E cut
 122
 123         } // end of loop over different outputs
 124
 125
 126         // Draw the plot
 127
 128         if (makeDraw)
 129         {
 130                 gStyle->SetOptTitle(0);
 131                 gStyle->SetOptStat(0);
 132                 gStyle->SetOptFit(0);
 133
 134                 TCanvas *c = new TCanvas("c","c",500,400);
 135                 //c->SetTopMargin(0.04);
 136                 //c->SetRightMargin(0.04);
 137                 //c->SetLeftMargin(0.181452);
 138                 //c->SetBottomMargin(0.134409);
 139                 c->SetBorderSize(0);
 140                 c->SetFillColor(0);
 141                 c->SetBorderMode(0);
 142                 c->SetFrameFillColor(0);
 143                 c->SetFrameBorderMode(0);
 144
 145                 /*
 146                  for (int i=0; i<NOUTPUTS; i++)
 147                  {
 148                  graph[i]->SetMarkerStyle(style[i]);
 149                  graph[i]->SetMarkerColor(color[i]);
 150                  graph[i]->SetLineColor(color[i]);
 151                  graph[i]->SetFillColor(0);
 152                  if (i == 0)
 153                  {
 154                  graph[i]->GetXaxis()->SetTitle("E (GeV)");
 155                  graph[i]->GetYaxis()->SetTitle("effic");
 156                  graph[i]->SetMaximum(1.0);
 157                  graph[i]->SetMinimum(0.0);
 158                  graph[i]->Draw("AP");
 159                  }
 160                  else
 161                  graph[i]->Draw("P");
 162                  }
 163                  */
 164                 for (int i=0; i<NOUTPUTS; i++)
 165                 {
 166                         effic[i]->SetMarkerStyle(style[i]);
 167                         effic[i]->SetMarkerColor(color[i]);
 168                         effic[i]->SetLineColor(color[i]);
 169                         effic[i]->SetFillColor(0);
 170                         effic[i]->SetTitle("efficiency; p_{T} (GeV/c); #epsilon");
 171                         if (i == 0)
 172                         {
 173                                 effic[i]->Draw();
 174                         }
 175                         else
 176                                 effic[i]->Draw("Psame");
 177                 }
 178
 179                 TLegend *leg = new TLegend(0.65,0.2,0.95,0.5);
 180                 leg->AddEntry(effic[0],"pions");
 181                 leg->AddEntry(effic[1],"kaons");
 182                 leg->AddEntry(effic[2],"protons");
 183                 //leg->AddEntry(effic[3],"muons");
 184                 leg->SetFillStyle(0);
 185                 leg->SetFillColor(0);
 186                 leg->SetBorderSize(0);
 187                 leg->SetTextSize(0.03);
 188                 leg->Draw();
 189         }
 190
 191         return effic[0];
 192 }