Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / efficiency / createEfficiency.C
CommitLineData
4ebdd20e 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
20using 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
53void DrawEfficiency(const Char_t* fileName, const Char_t* effFileName);
54void DrawCorrection(const Char_t* fileNamePYTHIA, const Char_t* fileNamePHOJET);
55TH1D* HistInvert(TH1D* hist);
56
57
58void 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//_______________________________________________________________________
300void 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//_______________________________________________________________________
598void 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//_______________________________________________________________________
678TH1D* 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}