PWGLF/SPECTRA/IdentifiedHighPt/efficiency/createEfficiency.C

   1 #include <TFile.h>
   2 #include <TH1.h>
   3 #include <TF1.h>
   4 #include <TCanvas.h>
   5 #include <TAxis.h>
   6 #include <TStyle.h>
   7 #include <TChain.h>
   8 #include <TTree.h>
   9 #include <TMath.h>
  10 #include <TClonesArray.h>
  11 #include <TLegend.h>
  12
  13 #include <AliXRDPROOFtoolkit.h>
  14
  15 #include "DebugClasses.C"
  16 #include "my_tools.C"
  17
  18 #include <iostream>
  19
  20 using namespace std;
  21
  22 /*
  23   To run code:
  24   ============
  25
  26   Info:
  27   * I did not recheck this code. For now I would just use the efficiency values I have.
  28   * The code could be made nicer. Esepcially some plots could be put in folders.
  29
  30   Use AliRoot because of AliXRDPROOFtoolkit:
  31   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC")
  32   gSystem->AddIncludePath("-I../lib")
  33   gSystem->AddIncludePath("-I../grid")
  34   gSystem->AddIncludePath("-I../macros")
  35   gROOT->SetMacroPath(".:../macros:../grid:../lib/")
  36   .L my_tools.C+
  37   .L DebugClasses.C+
  38   .L createEfficiency.C+
  39
  40   Examples of visualization:
  41   DrawEfficiency("lhc10d_eff_pythia.root", "eff.root")
  42
  43   // This is the correction I did for the low pT guys
  44   DrawCorrection("lhc10d_eff_pythia.root", "lhc10d_eff_phojet.root")
  45
  46
  47
  48   CreateEff("lhc10d_mc_all.dat", 0, "lhc10d_eff_all.root")
  49
  50
  51 */
  52
  53 void DrawEfficiency(const Char_t* fileName, const Char_t* effFileName);
  54 void DrawCorrection(const Char_t* fileNamePYTHIA, const Char_t* fileNamePHOJET);
  55 TH1D* HistInvert(TH1D* hist);
  56
  57
  58 void CreateEff(const Char_t* mcFileName, Int_t maxEvents, const Char_t* mcOutFileName,
  59                Float_t centLow=-20, Float_t centHigh=-5)
  60 {
  61   gStyle->SetOptStat(0);
  62   //
  63   // Create output
  64   //
  65   TFile* outFile = new TFile(mcOutFileName, "RECREATE");
  66
  67   const Int_t nPid = 7;
  68   TH1D* hMcIn[nPid]     = {0, 0, 0, 0, 0, 0, 0 };
  69   TH1D* hMcOut[nPid]    = {0, 0, 0, 0, 0, 0, 0 };
  70   TH1D* hMcSec[nPid]    = {0, 0, 0, 0, 0, 0, 0 };
  71   TH1D* hMcEff[nPid]     = {0, 0, 0, 0, 0, 0, 0 };
  72   TH1D* hMcInNeg[nPid]  = {0, 0, 0, 0, 0, 0, 0 };
  73   TH1D* hMcOutNeg[nPid] = {0, 0, 0, 0, 0, 0, 0 };
  74   TH1D* hMcSecNeg[nPid] = {0, 0, 0, 0, 0, 0, 0 };
  75   TH1D* hMcEffNeg[nPid]     = {0, 0, 0, 0, 0, 0, 0 };
  76   TH1D* hMcInPos[nPid]  = {0, 0, 0, 0, 0, 0, 0 };
  77   TH1D* hMcOutPos[nPid] = {0, 0, 0, 0, 0, 0, 0 };
  78   TH1D* hMcSecPos[nPid] = {0, 0, 0, 0, 0, 0, 0 };
  79   TH1D* hMcEffPos[nPid]     = {0, 0, 0, 0, 0, 0, 0 };
  80
  81   Int_t color[nPid] = {1, 2, 3, 4, 5, 1, 1};
  82
  83   const Int_t nPtBins = 68;
  84   Double_t xBins[nPtBins+1] = {0. ,  0.05, 0.1,  0.15, 0.2,  0.25, 0.3,  0.35, 0.4,  0.45,
  85                                0.5,  0.55, 0.6,  0.65, 0.7,  0.75, 0.8,  0.85, 0.9,  0.95,
  86                                1.0,  1.1 , 1.2,  1.3 , 1.4,  1.5 , 1.6,  1.7 , 1.8,  1.9 ,
  87                                2.0,  2.2 , 2.4,  2.6 , 2.8,  3.0 , 3.2,  3.4 , 3.6,  3.8 ,
  88                                4.0,  4.5 , 5.0,  5.5 , 6.0,  6.5 , 7.0,  8.0 , 9.0,  10.0,
  89                                11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
  90                                26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
  91
  92   for(Int_t pid = 0; pid < nPid; pid++) {
  93
  94     hMcIn[pid] = new TH1D(Form("hIn%d", pid), Form("Efficiency (pid %d)", pid),
  95                           nPtBins, xBins);
  96     hMcInNeg[pid] = new TH1D(Form("hInNeg%d", pid), Form("Efficiency (pid %d, q < 0)", pid),
  97                              nPtBins, xBins);
  98     hMcInPos[pid] = new TH1D(Form("hInPos%d", pid), Form("Efficiency (pid %d, q < 0)", pid),
  99                              nPtBins, xBins);
 100
 101     hMcIn[pid]->Sumw2();
 102     hMcIn[pid]->SetMarkerStyle(29);
 103     hMcIn[pid]->SetMarkerColor(color[pid]);
 104     hMcInNeg[pid]->Sumw2();
 105     hMcInNeg[pid]->SetMarkerStyle(24);
 106     hMcInNeg[pid]->SetMarkerColor(color[pid]);
 107     hMcInPos[pid]->Sumw2();
 108     hMcInPos[pid]->SetMarkerStyle(20);
 109     hMcInPos[pid]->SetMarkerColor(color[pid]);
 110
 111     hMcOut[pid] = new TH1D(Form("hMcOut%d", pid), Form("MC out (pid %d)", pid),
 112                            nPtBins, xBins);
 113     hMcOutNeg[pid] = new TH1D(Form("hMcOutNeg%d", pid), Form("MC out (pid %d, q < 0)", pid),
 114                               nPtBins, xBins);
 115     hMcOutPos[pid] = new TH1D(Form("hMcOutPos%d", pid), Form("MC out (pid %d, q < 0)", pid),
 116                               nPtBins, xBins);
 117
 118     hMcOut[pid]->Sumw2();
 119     hMcOut[pid]->SetMarkerStyle(29);
 120     hMcOut[pid]->SetMarkerColor(color[pid]);
 121     hMcOutNeg[pid]->Sumw2();
 122     hMcOutNeg[pid]->SetMarkerStyle(24);
 123     hMcOutNeg[pid]->SetMarkerColor(color[pid]);
 124     hMcOutPos[pid]->Sumw2();
 125     hMcOutPos[pid]->SetMarkerStyle(20);
 126     hMcOutPos[pid]->SetMarkerColor(color[pid]);
 127
 128     hMcSec[pid] = new TH1D(Form("hSec%d", pid), Form("Secondaries (pid %d)", pid),
 129                           nPtBins, xBins);
 130     hMcSecNeg[pid] = new TH1D(Form("hSecNeg%d", pid), Form("Secondaries (pid %d, q < 0)", pid),
 131                              nPtBins, xBins);
 132     hMcSecPos[pid] = new TH1D(Form("hSecPos%d", pid), Form("Secondaries (pid %d, q < 0)", pid),
 133                              nPtBins, xBins);
 134
 135     hMcSec[pid]->Sumw2();
 136     hMcSec[pid]->SetMarkerStyle(29);
 137     hMcSec[pid]->SetMarkerColor(color[pid]);
 138     hMcSecNeg[pid]->Sumw2();
 139     hMcSecNeg[pid]->SetMarkerStyle(24);
 140     hMcSecNeg[pid]->SetMarkerColor(color[pid]);
 141     hMcSecPos[pid]->Sumw2();
 142     hMcSecPos[pid]->SetMarkerStyle(20);
 143     hMcSecPos[pid]->SetMarkerColor(color[pid]);
 144
 145     hMcEff[pid] = new TH1D(Form("hEff%d", pid), Form("Efficiency (pid %d)", pid),
 146                           nPtBins, xBins);
 147     hMcEffNeg[pid] = new TH1D(Form("hEffNeg%d", pid), Form("Efficiency (pid %d, q < 0)", pid),
 148                              nPtBins, xBins);
 149     hMcEffPos[pid] = new TH1D(Form("hEffPos%d", pid), Form("Efficiency (pid %d, q < 0)", pid),
 150                              nPtBins, xBins);
 151
 152     hMcEff[pid]->Sumw2();
 153     hMcEff[pid]->SetMarkerStyle(29);
 154     hMcEff[pid]->SetMarkerColor(color[pid]);
 155     hMcEffNeg[pid]->Sumw2();
 156     hMcEffNeg[pid]->SetMarkerStyle(24);
 157     hMcEffNeg[pid]->SetMarkerColor(color[pid]);
 158     hMcEffPos[pid]->Sumw2();
 159     hMcEffPos[pid]->SetMarkerStyle(20);
 160     hMcEffPos[pid]->SetMarkerColor(color[pid]);
 161   }
 162
 163   TTree* Tree   = 0;
 164   if(strstr(mcFileName, ".dat")) {
 165
 166     AliXRDPROOFtoolkit tool;
 167     TChain* chain = tool.MakeChain(mcFileName,"tree", 0, 1000);
 168     chain->Lookup();
 169     Tree = chain;
 170   } else {
 171     TFile* mcFile = FindFileFresh(mcFileName);
 172     if(!mcFile)
 173       return;
 174
 175     Tree = (TTree*)mcFile->Get("tree");
 176   }
 177   if(!Tree)
 178     return;
 179
 180
 181   DeDxEvent* event = 0;
 182   TClonesArray* trackArray = 0;
 183   TClonesArray* mcTrackArray = 0;
 184   Tree->SetBranchAddress("event", &event);
 185   Tree->SetBranchAddress("track"  , &trackArray);
 186   Tree->SetBranchAddress("trackMC"  , &mcTrackArray);
 187   Int_t nEvents = Tree->GetEntries();
 188   cout << "Number of events: " << nEvents << endl;
 189
 190   if(maxEvents>0 && maxEvents < nEvents) {
 191
 192     nEvents = maxEvents;
 193     cout << "N events was reduced to: " << maxEvents << endl;
 194   }
 195
 196   Int_t currentRun = 0;
 197
 198   for(Int_t n = 0; n < nEvents; n++) {
 199
 200     Tree->GetEntry(n);
 201
 202     if((n+1)%1000000==0)
 203       cout << "Event: " << n+1 << "/" << nEvents << endl;
 204
 205     if(event->run != currentRun) {
 206
 207       cout << "New run: " << event->run << endl;
 208       currentRun = event->run;
 209     }
 210
 211     if(event->cent < centLow || event->cent > centHigh)
 212       continue;
 213
 214     const Int_t nMcTracks = mcTrackArray->GetEntries();
 215
 216     for(Int_t i = 0; i < nMcTracks; i++) {
 217
 218       DeDxTrackMC* trackMC = (DeDxTrackMC*)mcTrackArray->At(i);
 219
 220       // if(TMath::Abs(trackMC->pdgMC)==3312 || TMath::Abs(trackMC->pdgMC)==3334)
 221       //        continue; // Xi-!
 222
 223       hMcIn[0]->Fill(trackMC->ptMC);
 224       hMcIn[trackMC->pidMC]->Fill(trackMC->ptMC);
 225
 226       if(trackMC->qMC < 0) {
 227
 228         hMcInNeg[0]->Fill(trackMC->ptMC);
 229         hMcInNeg[trackMC->pidMC]->Fill(trackMC->ptMC);
 230       } else {
 231
 232         hMcInPos[0]->Fill(trackMC->ptMC);
 233         hMcInPos[trackMC->pidMC]->Fill(trackMC->ptMC);
 234       }
 235     }
 236
 237     const Int_t nTracks = trackArray->GetEntries();
 238
 239     for(Int_t i = 0; i < nTracks; i++) {
 240
 241       DeDxTrack* track = (DeDxTrack*)trackArray->At(i);
 242
 243       if(!(track->filter&1))
 244         continue;
 245
 246       // if(TMath::Abs(track->mother)==3312 || TMath::Abs(track->mother)==3334)
 247       //        continue; // Xi+- or Omega+-!
 248
 249       hMcOut[0]->Fill(track->pt);
 250       hMcOut[track->pid]->Fill(track->pt);
 251
 252       if(track->q < 0) {
 253
 254         hMcOutNeg[0]->Fill(track->pt);
 255         hMcOutNeg[track->pid]->Fill(track->pt);
 256       } else {
 257
 258         hMcOutPos[0]->Fill(track->pt);
 259         hMcOutPos[track->pid]->Fill(track->pt);
 260       }
 261
 262       if(track->primary==0) {
 263         hMcSec[0]->Fill(track->pt);
 264         hMcSec[track->pid]->Fill(track->pt);
 265
 266         if(track->q < 0) {
 267
 268           hMcSecNeg[0]->Fill(track->pt);
 269           hMcSecNeg[track->pid]->Fill(track->pt);
 270         } else {
 271
 272           hMcSecPos[0]->Fill(track->pt);
 273           hMcSecPos[track->pid]->Fill(track->pt);
 274         }
 275       }
 276     }
 277   }
 278
 279   TH1D* hMcInPiKP = (TH1D*)hMcIn[1]->Clone("hMcInPiKP");
 280   hMcInPiKP->Add(hMcIn[2]);
 281   hMcInPiKP->Add(hMcIn[3]);
 282   hMcInPiKP->Divide(hMcIn[0]);
 283
 284   for(Int_t pid = 0; pid < nPid; pid++) {
 285
 286     hMcSec[pid]->Divide(hMcSec[pid], hMcOut[pid]);
 287     hMcSecNeg[pid]->Divide(hMcSecNeg[pid], hMcOutNeg[pid]);
 288     hMcSecPos[pid]->Divide(hMcSecPos[pid], hMcOutPos[pid]);
 289
 290     hMcEff[pid]   ->Divide(hMcOut[pid], hMcIn[pid]);
 291     hMcEffNeg[pid]->Divide(hMcOutNeg[pid], hMcInNeg[pid]);
 292     hMcEffPos[pid]->Divide(hMcOutPos[pid], hMcInPos[pid]);
 293   }
 294
 295   outFile->Write();
 296   outFile->Close();
 297 }
 298
 299 //_______________________________________________________________________
 300 void DrawEfficiency(const Char_t* fileName, const Char_t* effFileName)
 301 {
 302   TFile* file = FindFileFresh(fileName);
 303   if(!file)
 304     return;
 305
 306   gStyle->SetOptStat(0);
 307
 308   const Double_t ptmin  =  3.01;
 309   const Double_t ptmax  = 19.999;
 310   //  const Double_t ptmax  = 49.999;
 311
 312   const Double_t effmin = 0.6;
 313   const Double_t effmax = 0.9;
 314
 315   // const Double_t effmin = 0.0;
 316   // const Double_t effmax = 1.5;
 317
 318   const Double_t secmin = 0.0;
 319   const Double_t secmax = 0.1;
 320
 321   const Int_t nPid = 7;
 322
 323   TH1D* hEff[nPid]     = {0, 0, 0, 0, 0, 0, 0 };
 324   TH1D* hEffNeg[nPid]  = {0, 0, 0, 0, 0, 0, 0 };
 325   TH1D* hEffPos[nPid]  = {0, 0, 0, 0, 0, 0, 0 };
 326   TH1D* hSec[nPid]    = {0, 0, 0, 0, 0, 0, 0 };
 327   TH1D* hSecNeg[nPid] = {0, 0, 0, 0, 0, 0, 0 };
 328   TH1D* hSecPos[nPid] = {0, 0, 0, 0, 0, 0, 0 };
 329   TH1D* hOut[nPid]    = {0, 0, 0, 0, 0, 0, 0 };
 330
 331   for(Int_t pid = 0; pid < nPid; pid++) {
 332
 333     hEff[pid] = (TH1D*)file->Get(Form("hEff%d", pid));
 334
 335     hEffNeg[pid] = (TH1D*)file->Get(Form("hEffNeg%d", pid));
 336
 337     hEffPos[pid] = (TH1D*)file->Get(Form("hEffPos%d", pid));
 338
 339     hSec[pid] = (TH1D*)file->Get(Form("hSec%d", pid));
 340
 341     hSecNeg[pid] = (TH1D*)file->Get(Form("hSecNeg%d", pid));
 342
 343     hSecPos[pid] = (TH1D*)file->Get(Form("hSecPos%d", pid));
 344
 345     hOut[pid] = (TH1D*)file->Get(Form("hMcOut%d", pid));
 346
 347     // hEff[pid]->Rebin(4);
 348     // hEffNeg[pid]->Rebin(4);
 349     // hEffPos[pid]->Rebin(4);
 350
 351     // hSec[pid]->Rebin(4);
 352     // hSecNeg[pid]->Rebin(4);
 353     // hSecPos[pid]->Rebin(4);
 354
 355     // hOut[pid]->Rebin(4);
 356
 357     // hEff[pid]->Scale(0.25);
 358     // hEffNeg[pid]->Scale(0.25);
 359     // hEffPos[pid]->Scale(0.25);
 360
 361     // hSec[pid]->Scale(0.25);
 362     // hSecNeg[pid]->Scale(0.25);
 363     // hSecPos[pid]->Scale(0.25);
 364
 365     // hOut[pid]->Scale(0.25);
 366
 367     hEff[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax);
 368     hEff[pid]->GetYaxis()->SetRangeUser(effmin, effmax);
 369     hEffNeg[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax);
 370     hEffNeg[pid]->GetYaxis()->SetRangeUser(effmin, effmax);
 371     hEffPos[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax);
 372     hEffPos[pid]->GetYaxis()->SetRangeUser(effmin, effmax);
 373     hSec[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax);
 374     hSec[pid]->GetYaxis()->SetRangeUser(secmin, secmax);
 375     hSecNeg[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax);
 376     hSecNeg[pid]->GetYaxis()->SetRangeUser(secmin, secmax);
 377     hSecPos[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax);
 378     hSecPos[pid]->GetYaxis()->SetRangeUser(secmin, secmax);
 379     hOut[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax);
 380   }
 381
 382   TCanvas* cChEff = new TCanvas("cChEff", "All efficiency", 800, 300);
 383   cChEff->Clear();
 384   cChEff->Divide(2, 1);
 385   cChEff->cd(1);
 386   //  TF1* pionEff = new TF1("pionEff", "pol0", 0.0, 50.0);
 387   TF1* chEff = new TF1("chEff", "[0]*(1-[1]/x)", 0.0, 50.0);
 388   chEff->SetLineColor(1);
 389   chEff->SetParameters(0.7, 0.5);
 390   hEff[0]->Fit(chEff, "", "", ptmin, ptmax);
 391
 392   cChEff->cd(2);
 393   hEffPos[0]->Draw();
 394   hEffNeg[0]->Draw("SAME");
 395   chEff->DrawCopy("SAME");
 396   cChEff->SaveAs("eff_charged.gif");
 397
 398   TCanvas* cPionEff = new TCanvas("cPionEff", "Pion efficiency", 800, 300);
 399   cPionEff->Clear();
 400   cPionEff->Divide(2, 1);
 401   cPionEff->cd(1);
 402   //  TF1* pionEff = new TF1("pionEff", "pol0", 0.0, 50.0);
 403   TF1* pionEff = new TF1("pionEff", "[0]*(1-[1]/x)", 0.0, 50.0);
 404   pionEff->SetLineColor(2);
 405   pionEff->SetParameters(0.7, 0.5);
 406   hEff[1]->Fit(pionEff, "", "", ptmin, ptmax);
 407
 408   cPionEff->cd(2);
 409   hEffPos[1]->Draw();
 410   hEffNeg[1]->Draw("SAME");
 411   pionEff->DrawCopy("SAME");
 412   cPionEff->SaveAs("eff_pion.gif");
 413
 414   TCanvas* cKaonEff = new TCanvas("cKaonEff", "Kaon efficiency", 800, 300);
 415   cKaonEff->Clear();
 416   cKaonEff->Divide(2, 1);
 417   cKaonEff->cd(1);
 418   //  TF1* kaonEff = new TF1("kaonEff", "pol0", 0.0, 50.0);
 419   TF1* kaonEff = new TF1("kaonEff", "[0]*(1-[1]/x)", 0.0, 50.0);
 420   kaonEff->SetLineColor(3);
 421   kaonEff->SetParameters(0.7, 0.5);
 422   hEff[2]->Fit(kaonEff, "", "", ptmin, ptmax);
 423
 424   cKaonEff->cd(2);
 425   hEffPos[2]->Draw();
 426   hEffNeg[2]->Draw("SAME");
 427   kaonEff->DrawCopy("SAME");
 428   cKaonEff->SaveAs("eff_kaon.gif");
 429
 430   TCanvas* cProtonEff = new TCanvas("cProtonEff", "Proton efficiency", 800, 300);
 431   cProtonEff->Clear();
 432   cProtonEff->Divide(2, 1);
 433   cProtonEff->cd(1);
 434   //  TF1* protonEff = new TF1("protonEff", "pol0", 0.0, 50.0);
 435   TF1* protonEff = new TF1("protonEff", "[0]*(1-[1]/x)", 0.0, 50.0);
 436   protonEff->SetLineColor(4);
 437   protonEff->SetParameters(0.7, 0.5);
 438   hEff[3]->Fit(protonEff, "", "", ptmin, ptmax);
 439
 440   cProtonEff->cd(2);
 441   hEffPos[3]->Draw();
 442   hEffNeg[3]->Draw("SAME");
 443   protonEff->DrawCopy("SAME");
 444   cProtonEff->SaveAs("eff_proton.gif");
 445
 446   TCanvas* cAllEff = new TCanvas("cAllEff", "All efficiency", 400, 300);
 447   cAllEff->Clear();
 448
 449   TH1D* hDummy = (TH1D*)hEff[1]->Clone("hDummy");
 450   hDummy->Reset();
 451   hDummy->SetTitle("Efficiency vs p_{T}; p_{T} [GeV/c]; Efficiency");
 452   hDummy->Draw();
 453
 454   chEff->DrawCopy("SAME");
 455   pionEff->DrawCopy("SAME");
 456   kaonEff->DrawCopy("SAME");
 457   protonEff->DrawCopy("SAME");
 458   cAllEff->SaveAs("eff_all.gif");
 459
 460   // TCanvas* cPionSec = new TCanvas("cPionSec", "Pion sec", 800, 300);
 461   // cPionSec->Clear();
 462   // cPionSec->Divide(2, 1);
 463   // cPionSec->cd(1);
 464   // TF1* pionSec = new TF1("pionSec", "pol0", 0.0, 50.0);
 465   // hSec[1]->Fit(pionSec, "", "", ptmin, ptmax);
 466
 467   // cPionSec->cd(2);
 468   // hSecPos[1]->Draw();
 469   // hSecNeg[1]->Draw("SAME");
 470   // pionSec->Draw("SAME");
 471   // cPionSec->SaveAs("sec_pion.gif");
 472
 473   // TCanvas* cKaonSec = new TCanvas("cKaonSec", "Kaon sec", 800, 300);
 474   // cKaonSec->Clear();
 475   // cKaonSec->Divide(2, 1);
 476   // cKaonSec->cd(1);
 477   // TF1* kaonSec = new TF1("kaonSec", "pol0", 0.0, 50.0);
 478   // hSec[2]->Fit(kaonSec, "", "", ptmin, ptmax);
 479
 480   // cKaonSec->cd(2);
 481   // hSecPos[2]->Draw();
 482   // hSecNeg[2]->Draw("SAME");
 483   // kaonSec->Draw("SAME");
 484   // cKaonSec->SaveAs("sec_kaon.gif");
 485
 486   // TCanvas* cProtonSec = new TCanvas("cProtonSec", "Proton sec", 800, 300);
 487   // cProtonSec->Clear();
 488   // cProtonSec->Divide(2, 1);
 489   // cProtonSec->cd(1);
 490   // TF1* protonSec = new TF1("protonSec", "pol0", 0.0, 50.0);
 491   // hSec[3]->Fit(protonSec, "", "", ptmin, ptmax);
 492
 493   // cProtonSec->cd(2);
 494   // hSecPos[3]->Draw();
 495   // hSecNeg[3]->Draw("SAME");
 496   // protonSec->Draw("SAME");
 497   // cProtonSec->SaveAs("sec_proton.gif");
 498
 499
 500   TCanvas* cEffRatioPi = new TCanvas("cEffRatioPi", "eff pi / eff all", 400, 300);
 501   cEffRatioPi->Clear();
 502   //  TF1* pionEff = new TF1("pionEff", "pol0", 0.0, 50.0);
 503   TF1* effRatioPi = new TF1("effRatioPi", "pol0", 0.0, 50.0);
 504   effRatioPi->SetLineColor(6);
 505   effRatioPi->SetParameters(0.7, 0.5);
 506   TH1D* hEffRatioPi = (TH1D*)hEff[1]->Clone("hEffRatioPi");
 507   hEffRatioPi->SetTitle("; p_{T} [GeV/c]; #epsilon_{ch}/#epsilon_{pi}");
 508   hEffRatioPi->Divide(hEff[0], hEff[1], 1, 1, "B");
 509   //  hEffRatioPi->Divide(hEff[1], hEff[0]);
 510   hEffRatioPi->GetXaxis()->SetRangeUser(0.0, 19.99);
 511   hEffRatioPi->GetYaxis()->SetRangeUser(0.8, 1.1);
 512   hEffRatioPi->SetStats(kTRUE);
 513   hEffRatioPi->Fit(effRatioPi, "", "", ptmin, ptmax);
 514   cEffRatioPi->SaveAs("eff_ratio_pi.gif");
 515   cEffRatioPi->SaveAs("eff_ratio_pi.pdf");
 516
 517   TCanvas* cEffRatioK = new TCanvas("cEffRatioK", "eff K / eff all", 400, 300);
 518   cEffRatioK->Clear();
 519   TF1* effRatioK = new TF1("effRatioK", "exp(-[1]*x)+[0]", 0.0, 50.0);
 520   effRatioK->SetParameters(1.0, 1.0);
 521   effRatioK->SetLineColor(6);
 522   TH1D* hEffChargedRebinned = (TH1D*)hEff[0]->Clone("hEffChargedRebinned");
 523   hEffChargedRebinned->Rebin(2);
 524   TH1D* hEffRatioK = (TH1D*)hEff[2]->Clone("hEffRatioK");
 525   hEffRatioK->Rebin(2);
 526   hEffRatioK->SetTitle("; p_{T} [GeV/c]; #epsilon_{ch}/#epsilon_{K}");
 527   hEffRatioK->Divide(hEffChargedRebinned, hEffRatioK, 1, 1, "B");
 528   //  hEffRatioK->Divide(hEff[1], hEff[0]);
 529   hEffRatioK->GetXaxis()->SetRangeUser(0.0, 19.99);
 530   hEffRatioK->GetYaxis()->SetRangeUser(0.8, 1.1);
 531   hEffRatioK->SetStats(kTRUE);
 532   hEffRatioK->Fit(effRatioK, "", "", ptmin, ptmax);
 533   cEffRatioK->SaveAs("eff_ratio_K.gif");
 534   cEffRatioK->SaveAs("eff_ratio_K.pdf");
 535
 536   TCanvas* cEffRatioP = new TCanvas("cEffRatioP", "eff p / eff all", 400, 300);
 537   cEffRatioP->Clear();
 538   TF1* effRatioP = new TF1("effRatioP", "pol0", 0.0, 50.0);
 539   effRatioP->SetLineColor(6);
 540   effRatioP->SetParameters(0.7, 0.5);
 541   // TH1D* hEffChargedRebinned = (TH1D*)hEff[0]->Clone("hEffChargedRebinned");
 542   // hEffChargedRebinned->Rebin(2);
 543   TH1D* hEffRatioP = (TH1D*)hEff[3]->Clone("hEffRatioP");
 544   hEffRatioP->Rebin(2);
 545   hEffRatioP->SetTitle("; p_{T} [GeV/c]; #epsilon_{ch}/#epsilon_{p}");
 546   hEffRatioP->Divide(hEffChargedRebinned, hEffRatioP, 1, 1, "B");
 547   //  hEffRatioP->Divide(hEff[1], hEff[0]);
 548   hEffRatioP->GetXaxis()->SetRangeUser(0.0, 19.99);
 549   hEffRatioP->GetYaxis()->SetRangeUser(0.8, 1.1);
 550   hEffRatioP->SetStats(kTRUE);
 551   hEffRatioP->Fit(effRatioP, "", "", ptmin, ptmax);
 552   cEffRatioP->SaveAs("eff_ratio_p.gif");
 553   cEffRatioP->SaveAs("eff_ratio_p.pdf");
 554
 555
 556   // TCanvas* cMuonRatio = new TCanvas("cMuonRatio", "muon / all", 400, 300);
 557   // cMuonRatio->Clear();
 558   // //  TF1* pionMuon = new TF1("pionMuon", "pol0", 0.0, 50.0);
 559   // TF1* muonRatio = new TF1("muonRatio", "pol0", 0.0, 50.0);
 560   // muonRatio->SetLineColor(1);
 561   // muonRatio->SetParameter(0, 0.01);
 562   // TH1D* hMuonRatio = (TH1D*)hOut[5]->Clone("hMuonRatio");
 563   // hMuonRatio->SetTitle("; p_{T} [GeV/c]; muon/all");
 564   // hMuonRatio->Divide(hOut[5], hOut[0], 1, 1, "B");
 565   // hMuonRatio->GetYaxis()->SetRangeUser(0.0, 0.05);
 566   // hMuonRatio->SetStats(kTRUE);
 567   // hMuonRatio->Fit(muonRatio, "", "", ptmin, ptmax);
 568   // cMuonRatio->SaveAs("muon_ratio.gif");
 569   // cMuonRatio->SaveAs("muon_ratio.pdf");
 570
 571   // TCanvas* cElectronRatio = new TCanvas("cElectronRatio", "electron / all", 400, 300);
 572   // cElectronRatio->Clear();
 573   // //  TF1* pionElectron = new TF1("pionElectron", "pol0", 0.0, 50.0);
 574   // TF1* electronRatio = new TF1("electronRatio", "pol0", 0.0, 50.0);
 575   // electronRatio->SetLineColor(1);
 576   // electronRatio->SetParameter(0, 0.01);
 577   // TH1D* hElectronRatio = (TH1D*)hOut[4]->Clone("hElectronRatio");
 578   // hElectronRatio->SetTitle("; p_{T} [GeV/c]; electron/all");
 579   // hElectronRatio->Divide(hOut[4], hOut[0], 1, 1, "B");
 580   // hElectronRatio->GetYaxis()->SetRangeUser(0.0, 0.05);
 581   // hElectronRatio->SetStats(kTRUE);
 582   // hElectronRatio->SetMarkerColor(6);
 583   // hElectronRatio->Fit(electronRatio, "", "", ptmin, ptmax);
 584   // cElectronRatio->SaveAs("electron_ratio.gif");
 585   // cElectronRatio->SaveAs("electron_ratio.pdf");
 586
 587
 588   TFile* effFile = new TFile(effFileName, "RECREATE");
 589   chEff->Write();
 590   pionEff->Write();
 591   kaonEff->Write();
 592   protonEff->Write();
 593   effFile->Close();
 594 }
 595
 596
 597 //_______________________________________________________________________
 598 void DrawCorrection(const Char_t* fileNamePYTHIA, const Char_t* fileNamePHOJET)
 599 {
 600   TFile* filePYTHIA = FindFileFresh(fileNamePYTHIA);
 601   if(!filePYTHIA)
 602     return;
 603
 604   TFile* filePHOJET = FindFileFresh(fileNamePHOJET);
 605   if(!filePHOJET)
 606     return;
 607
 608   gStyle->SetOptStat(0);
 609
 610   const Double_t ptmin  =  0.0;
 611   const Double_t ptmax  = 4.999;
 612   //  const Double_t ptmax  = 49.999;
 613
 614   const Double_t effmin = 0.95;
 615   const Double_t effmax = 1.12;
 616
 617   const Int_t nPid = 7;
 618
 619   TH1D* hInPYTHIA[nPid]     = {0, 0, 0, 0, 0, 0, 0 };
 620   TH1D* hInPHOJET[nPid]     = {0, 0, 0, 0, 0, 0, 0 };
 621
 622   for(Int_t pid = 0; pid < nPid; pid++) {
 623
 624     hInPYTHIA[pid] = (TH1D*)filePYTHIA->Get(Form("hIn%d", pid));
 625     hInPYTHIA[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax);
 626     hInPYTHIA[pid]->GetYaxis()->SetRangeUser(effmin, effmax);
 627
 628     hInPHOJET[pid] = (TH1D*)filePHOJET->Get(Form("hIn%d", pid));
 629     hInPHOJET[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax);
 630     hInPHOJET[pid]->GetYaxis()->SetRangeUser(effmin, effmax);
 631   }
 632
 633   hInPYTHIA[1]->Add(hInPYTHIA[2]);
 634   hInPYTHIA[1]->Add(hInPYTHIA[3]);
 635   hInPYTHIA[0]->Divide(hInPYTHIA[1], hInPYTHIA[0], 1, 1, "B");
 636
 637   TH1D* histPYTHIA = HistInvert(hInPYTHIA[0]);
 638   histPYTHIA->SetLineColor(2);
 639   histPYTHIA->SetMarkerColor(2);
 640
 641   hInPHOJET[1]->Add(hInPHOJET[2]);
 642   hInPHOJET[1]->Add(hInPHOJET[3]);
 643   hInPHOJET[0]->Divide(hInPHOJET[1], hInPHOJET[0], 1, 1, "B");
 644
 645   TH1D* histPHOJET = HistInvert(hInPHOJET[0]);
 646   histPHOJET->SetLineColor(4);
 647   histPHOJET->SetMarkerColor(4);
 648
 649   TCanvas* cChCorrection = new TCanvas("cChCorrection", "All efficiency", 400, 300);
 650   cChCorrection->Clear();
 651
 652   cChCorrection->SetGridy();
 653
 654   histPYTHIA->GetYaxis()->SetRangeUser(effmin, effmax);
 655   histPYTHIA->SetTitle("N_{ch}/(#pi + K + p) for primaries at generator level; p_{T} [GeV/c]; Ratio");
 656   histPYTHIA->Draw();
 657
 658   histPHOJET->Draw("SAME");
 659
 660   TLegend* legend = new TLegend(0.55, 0.22, 0.79, 0.42);
 661   legend->SetBorderSize(0);
 662   legend->SetFillColor(0);
 663   legend->AddEntry(histPYTHIA, "PYTHIA", "P");
 664   legend->AddEntry(histPHOJET, "PHOJET", "P");
 665   legend->Draw();
 666
 667   cChCorrection->SaveAs("charged_over_piKp.gif");
 668
 669   TFile* file = new TFile("correction.root", "RECREATE");
 670   histPYTHIA->SetName("histPYTHIA");
 671   histPHOJET->SetName("histPHOJET");
 672   histPYTHIA->Write();
 673   histPHOJET->Write();
 674   file->Close();
 675 }
 676
 677 //_______________________________________________________________________
 678 TH1D* HistInvert(TH1D* hist)
 679 {
 680   TH1D* histNew = new TH1D(*hist);
 681   histNew->Reset();
 682   histNew->SetName(Form("%s_inv", hist->GetName()));
 683   histNew->SetTitle(Form("%s (inverted)", hist->GetTitle()));
 684
 685   const Int_t nBins = hist->GetNbinsX();
 686
 687   for(Int_t i = 1; i <= nBins; i++) {
 688
 689     if(hist->GetBinContent(i) == 0)
 690       continue;
 691
 692     histNew->SetBinContent(i, 1.0/hist->GetBinContent(i));
 693     histNew->SetBinError(i, hist->GetBinError(i)/hist->GetBinContent(i)/hist->GetBinContent(i));
 694   }
 695
 696   return histNew;
 697 }