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

   1 /*
   2  *  CalcCorrAcc.C
   3  */
   4
   5 void CalcCorrAcc(const char* fileName, Bool_t makeDraw=kFALSE, Float_t* yy=0, Float_t* eyy=0)
   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 NCUTS = 3; // Choose the number of cuts to check
  15         Float_t eCut[NCUTS] = {0.1,0.25,0.5}; // Choose the value of the cuts
  16         //Float_t eCut[NCUTS] = {0.5}; // Choose the value of the cuts
  17
  18         Int_t style[NCUTS] = {20,21,22};
  19         Int_t color[NCUTS] = {1,2,4};
  20         //Int_t style[NCUTS] = {20};
  21         //Int_t color[NCUTS] = {1};
  22
  23         const Int_t fgNumOfEBins = 78; // 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
  26         // declare histograms and graphs
  27         TH2F *hist;
  28         TH2F *histAllEM;
  29         TH2F *histAccEM;
  30
  31         TGraphErrors *graph[NCUTS];
  32
  33         // retrieve the input list of histogram. Check the TList name in the input file.
  34         TList *list = (TList*) file->Get("out1");
  35
  36         // retrieve the histograms in the list. Check the name of the histograms
  37         hist = (TH2F*)list->FindObject("fHistElectron_EtaE_ET_EmcalMC");
  38         histAllEM = new TH2F(*hist);
  39         hist = (TH2F*)list->FindObject("fHistConvElectron_EtaE_ET_EmcalMC");
  40         histAllEM->Add(hist);
  41         hist = (TH2F*)list->FindObject("fHistScatElectron_EtaE_ET_EmcalMC");
  42         histAllEM->Add(hist);
  43         hist = (TH2F*)list->FindObject("fHistGamma_EtaE_ET_EmcalMC");
  44         histAllEM->Add(hist);
  45         hist = (TH2F*)list->FindObject("fHistAnnihGamma_EtaE_ET_EmcalMC");
  46         histAllEM->Add(hist);
  47         hist = (TH2F*)list->FindObject("fHistScatGamma_EtaE_ET_EmcalMC");
  48         histAllEM->Add(hist);
  49
  50         hist = (TH2F*)list->FindObject("fHistElectronAcc_EtaE_ET_EmcalMC");
  51         histAccEM = new TH2F(*hist);
  52         hist = (TH2F*)list->FindObject("fHistConvElectronAcc_EtaE_ET_EmcalMC");
  53         histAccEM->Add(hist);
  54         hist = (TH2F*)list->FindObject("fHistScatElectronAcc_EtaE_ET_EmcalMC");
  55         histAccEM->Add(hist);
  56         hist = (TH2F*)list->FindObject("fHistGammaAcc_EtaE_ET_EmcalMC");
  57         histAccEM->Add(hist);
  58         hist = (TH2F*)list->FindObject("fHistScatGammaAcc_EtaE_ET_EmcalMC");
  59         histAccEM->Add(hist);
  60         hist = (TH2F*)list->FindObject("fHistAnnihGammaAcc_EtaE_ET_EmcalMC");
  61         histAccEM->Add(hist);
  62
  63
  64         // ********************************************
  65
  66         Float_t x[fgNumOfEtaBins]={0}, ex[fgNumOfEtaBins]={0};
  67         Float_t y[fgNumOfEtaBins]={0}, ey[fgNumOfEtaBins]={0};
  68
  69         // returning the output that is in memory
  70         yy=y;
  71         eyy = ey;
  72
  73         // loop over different E cuts
  74         for (int iCut=0; iCut<NCUTS; iCut++)
  75         {
  76                 // loop over eta bins
  77                 for (int iy=0; iy<fgNumOfEtaBins; iy++)
  78                 {
  79                         // initialize ET variables for a new particle species
  80                         Float_t E=0, ET=0, ET_AllEM=0, ET_AccEM=0, errorSqET_AllEM=0, errorSqET_AccEM=0;
  81                         x[iy]=0;
  82                         y[iy]=0;
  83                         ex[iy]=0;
  84                         ey[iy]=0;
  85
  86                         // loop over E bins
  87                         for (int ix=0; ix<fgNumOfEBins; ix++)
  88                         {
  89                                 E = histAccEM->GetXaxis()->GetBinLowEdge(ix+1);
  90
  91                                 if (E >= eCut[iCut])
  92                                 {
  93                                         ET = histAccEM->GetBinContent(ix+1,iy+1); // sum over all E bins in order to get total ET
  94                                         errorSqET_AccEM += pow(E_err(E)*(ET/E),2);
  95                                         ET_AccEM += ET;
  96                                 }
  97
  98                                 ET = histAllEM->GetBinContent(ix+1,iy+1); // sum over all E bins in order to get total ET
  99                                 if (E > 0)
 100                                         errorSqET_AllEM += pow(E_err(E)*(ET/E),2);
 101                                 ET_AllEM += ET;
 102                         } // end of loop over E bins
 103
 104                         x[iy] = histAccEM->GetYaxis()->GetBinCenter(iy+1); // x coordinate of eta bin
 105
 106                         if ( (ET_AllEM > 0) && (ET_AccEM > 0) )
 107                         {
 108                                 y[iy] = ET_AccEM/ET_AllEM; // y coordinate of eta bin (for a given particle species)
 109                                 ey[iy]=y[iy]*sqrt(errorSqET_AllEM/pow(ET_AllEM,2) + errorSqET_AccEM/pow(ET_AccEM,2));
 110                         }
 111                         else
 112                         {
 113                                 y[iy] = 0;
 114                                 ey[iy] = 0;
 115                         }
 116                 } // end of loop over eta bins
 117
 118                 graph[iCut] = new TGraphErrors(fgNumOfEtaBins,x,y,ex,ey); // graphic of ET(>E_cut)/ET(total) for a given E cut (all particle species combined)
 119         } // end of loop over different E cuts
 120
 121
 122         // Draw the plot
 123         if (makeDraw)
 124         {
 125                 gStyle->SetOptTitle(0);
 126                 gStyle->SetOptStat(0);
 127                 gStyle->SetOptFit(0);
 128
 129                 TCanvas *c = new TCanvas("c","c",500,400);
 130                 //c->SetTopMargin(0.04);
 131                 //c->SetRightMargin(0.04);
 132                 //c->SetLeftMargin(0.181452);
 133                 //c->SetBottomMargin(0.134409);
 134                 c->SetBorderSize(0);
 135                 c->SetFillColor(0);
 136                 c->SetBorderMode(0);
 137                 c->SetFrameFillColor(0);
 138                 c->SetFrameBorderMode(0);
 139
 140                 for (int i=0; i<NCUTS; i++)
 141                 {
 142                         graph[i]->SetMarkerStyle(style[i]);
 143                         graph[i]->SetMarkerColor(color[i]);
 144                         graph[i]->SetLineColor(color[i]);
 145                         graph[i]->SetFillColor(0);
 146                         if (i == 0)
 147                         {
 148                                 graph[i]->GetXaxis()->SetTitle("#eta");
 149                                 graph[i]->GetYaxis()->SetTitle("E_{T}^{EM Acceptance}/E_{T}^{EM All}");
 150                                 graph[i]->SetMaximum(0.1);
 151                                 graph[i]->SetMinimum(0.0);
 152                                 graph[i]->Draw("AP");
 153                         }
 154                         else
 155                                 graph[i]->Draw("P");
 156                 }
 157
 158                 TLegend *leg = new TLegend(0.65,0.2,0.95,0.5);
 159                 leg->AddEntry(graph[0],"E>100 MeV");
 160                 leg->AddEntry(graph[1],"E>250 MeV");
 161                 leg->AddEntry(graph[2],"E>500 MeV");
 162                 leg->SetFillStyle(0);
 163                 leg->SetFillColor(0);
 164                 leg->SetBorderSize(0);
 165                 leg->SetTextSize(0.03);
 166                 leg->Draw();
 167         }
 168 }
 169
 170 Float_t E_err(Float_t E)
 171 {
 172         return (E*sqrt(pow(4.35/E,2)+pow(9.07,2)/E+pow(1.63,2)))/100;
 173 }