Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / 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 }