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

   1 /*
   2  *  CalcET.C
   3  */
   4
   5 void CalcET(const char* fileName, const char* fileNameMC)
   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.25,0.5,0.75}; // Choose the value of the cuts
  16         Float_t ptCut[NCUTS] = {0.25,0.5,0.75}; // Choose the value of the cuts
  17         Double_t matchCut = 0.02;
  18
  19         const Int_t fgNumOfEBins = 78; // Check the number of eta bins in the histograms
  20         const Int_t fgNumOfEtaBins = 16; // Check the number of E bins in the histograms
  21         const Int_t fgNumOfPtBins = 111; // Check the number of pT bins in the histograms
  22
  23         // declare histograms and graphs
  24         TGraphErrors *graph[NCUTS];
  25
  26         // retrieve the input list of histogram. Check the TList name in the input file.
  27         TList *list = (TList*) file->Get("out1");
  28
  29         // retrieve the histograms in the list. Check the name of the histograms
  30         THnSparseD* histAll = (THnSparseD*)list->FindObject("fHistAllRec_ETDep_EmcalRec");
  31
  32         THnSparseD* hnSparse = (THnSparseD*)list->FindObject("fHistMuon_ETDep_EmcalRec");
  33         hnSparse->GetAxis(6)->SetRangeUser(0.,matchCut);
  34         TH3D* histMuon = hnSparse->Projection(2,0,1);
  35
  36         hnSparse = (THnSparseD*)list->FindObject("fHistMuon_ETDep_EmcalRec");
  37         hnSparse->GetAxis(6)->SetRangeUser(0.,matchCut);
  38
  39         hnSparse = (THnSparseD*)list->FindObject("fHistMuon_ETDep_EmcalRec");
  40         hnSparse->GetAxis(6)->SetRangeUser(0.,matchCut);
  41         TH3D* histPion = ((THnSparseD*)list->FindObject("fHistPion_ETDep_EmcalRec"))->Projection(2,0,1);
  42
  43         hnSparse = (THnSparseD*)list->FindObject("fHistMuon_ETDep_EmcalRec");
  44         hnSparse->GetAxis(6)->SetRangeUser(0.,matchCut);
  45         TH3D* histKaon = ((THnSparseD*)list->FindObject("fHistKaon_ETDep_EmcalRec"))->Projection(2,0,1);
  46
  47         hnSparse = (THnSparseD*)list->FindObject("fHistMuon_ETDep_EmcalRec");
  48         hnSparse->GetAxis(6)->SetRangeUser(0.,matchCut);
  49         TH3D* histProton = ((THnSparseD*)list->FindObject("fHistProton_ETDep_EmcalRec"))->Projection(2,0,1);
  50
  51         TH3D* histCharged = new TH3D(*histMuon);
  52         histCharged->Add(histPion);
  53         histCharged->Add(histKaon);
  54         histCharged->Add(histProton);
  55
  56         TEfficiency *efficClu = CalcEffic(fileNameMC);
  57         TEfficiency *efficChHad = CalcEfficHadrons(fileNameMC);
  58         Double_t f_acc = CalcCorrAcc(fileNameMC);
  59         Double_t f_Ecut = CalcCorrEcut(fileNameMC);
  60         Double_t f_neutral = CalcCorrNeutral(fileNameMC);
  61
  62         // ********************************************
  63
  64         Float_t x[fgNumOfEtaBins]={0}, ex[fgNumOfEtaBins]={0};
  65         Float_t y[fgNumOfEtaBins]={0}, ey[fgNumOfEtaBins]={0};
  66
  67         // loop over different E cuts
  68         for (int iCut=0; iCut<NCUTS; iCut++)
  69         {
  70                 // loop over eta bins
  71                 for (int iy=0; iy<fgNumOfEtaBins; iy++)
  72                 {
  73                         // initialize ET variables for a new particle species
  74                         Float_t E=0, pt=0, ET=0, ET_All=0, ET_Had=0, errorSqET_All=0, errorSqET_Had=0;
  75                         x[iy]=0;
  76                         y[iy]=0;
  77                         ex[iy]=0;
  78                         ey[iy]=0;
  79
  80                         // loop over E bins
  81                         for (int ix=0; ix<fgNumOfEBins; ix++)
  82                         {
  83                                 E = histAll->GetXaxis()->GetBinCenter(ix+1);
  84                                 if (E > eCut[iCut])
  85                                 {
  86                                         ET = histAll->GetBinContent(ix+1,iy+1); // sum over all E bins in order to get total ET
  87                                         errorSqET_All += pow(E_err(E)*(ET/E),2);
  88                                         ET_All += ET/efficClu->GetEfficiency(ix);
  89
  90                                 }
  91
  92                                 // loop over pt bins
  93                                 for (int iz=0; iz<fgNumOfPtBins; iz++)
  94                                 {
  95                                         pt = histCharged->GetZaxis()->GetBinCenter(iz+1);
  96                                         if (pt > ptCut[iCut])
  97                                         {
  98                                                 ET = histCharged->GetBinContent(ix+1,iy+1,iz+1); // sum over all pt bins in order to get total ET
  99                                                 errorSqET_Had += pow(E_err(E)*(ET/E),2);
 100                                                 ET_Had += ET/efficChHad->GetEfficiency(iz);
 101                                         }
 102                                 } // end of loop over pt bins
 103
 104                         } // end of loop over E bins
 105
 106                         x[iy] = histChHad->GetYaxis()->GetBinCenter(iy+1); // x coordinate of eta bin
 107
 108                         y[iy] = (ET_All-(ET_Had/f_neutral))/(f_acc*f_Ecut); // y coordinate of eta bin (for a given particle species)
 109
 110                         ey[iy]=(sqrt(errorSqET_AllHad + errorSqET_Had/pow(f_neutral,2))/(f_acc*f_Ecut);
 111
 112                 } // end of loop over eta bins
 113
 114                 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)
 115         } // end of loop over different E cuts
 116
 117
 118         // Draw the plot
 119
 120         gStyle->SetOptTitle(0);
 121         gStyle->SetOptStat(0);
 122         gStyle->SetOptFit(0);
 123
 124         TCanvas *c = new TCanvas("c","c",500,400);
 125         //c->SetTopMargin(0.04);
 126         //c->SetRightMargin(0.04);
 127         //c->SetLeftMargin(0.181452);
 128         //c->SetBottomMargin(0.134409);
 129         c->SetBorderSize(0);
 130         c->SetFillColor(0);
 131         c->SetBorderMode(0);
 132         c->SetFrameFillColor(0);
 133         c->SetFrameBorderMode(0);
 134
 135         Int_t style[NCUTS] = {20,21,22};
 136         Int_t color[NCUTS] = {2,3,4};
 137
 138         for (int i=0; i<NCUTS; i++)
 139         {
 140                 graph[i]->SetMarkerStyle(style[i]);
 141                 graph[i]->SetMarkerColor(color[i]);
 142                 graph[i]->SetLineColor(color[i]);
 143                 graph[i]->SetFillColor(0);
 144                 if (i == 0)
 145                 {
 146                         graph[i]->GetXaxis()->SetTitle("#eta");
 147                         graph[i]->GetYaxis()->SetTitle("E_{T}^{Neutral+Charged Hadrons}/E_{T}^{Charged hadrons}");
 148                         //graph[i]->SetMaximum(1.0);
 149                         graph[i]->SetMinimum(0.0);
 150                         graph[i]->Draw("AP");
 151                 }
 152                 else
 153                         graph[i]->Draw("P");
 154         }
 155
 156         TLegend *leg = new TLegend(0.65,0.2,0.95,0.5);
 157         leg->AddEntry(graph[0],"E>250 MeV");
 158         leg->AddEntry(graph[1],"E>500 MeV");
 159         leg->AddEntry(graph[2],"E>750 MeV");
 160         leg->SetFillStyle(0);
 161         leg->SetFillColor(0);
 162         leg->SetBorderSize(0);
 163         leg->SetTextSize(0.03);
 164         leg->Draw();
 165 }
 166
 167 Float_t E_err(Float_t E)
 168 {
 169         return (E*sqrt(pow(4.35/E,2)+pow(9.07,2)/E+pow(1.63,2)))/100;
 170 }