updated dndeta analysis
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / drawPlots.C
CommitLineData
1afae8ff 1/* $Id$ */
2
0ab29cfa 3#if !defined(__CINT__) || defined(__MAKECINT__)
4
5#include "AliPWG0Helper.h"
6#include "dNdEtaAnalysis.h"
7#include "AlidNdEtaCorrection.h"
8
9#include <TCanvas.h>
10#include <TFile.h>
11#include <TH1.h>
12#include <TH2F.h>
13#include <TH3F.h>
14#include <TLine.h>
15#include <TSystem.h>
16
17#endif
18
25db2d85 19Int_t gMax = 5;
20
0ab29cfa 21extern TSystem* gSystem;
22
dd367a14 23void loadlibs()
24{
69b09e3b 25 gSystem->Load("libVMC");
26 gSystem->Load("libTree");
27 gSystem->Load("libSTEERBase");
28 gSystem->Load("libESD");
29 gSystem->Load("libAOD");
dd367a14 30 gSystem->Load("libANALYSIS");
69b09e3b 31 gSystem->Load("libANALYSISalice");
dd367a14 32 gSystem->Load("libPWG0base");
33}
34
0ab29cfa 35void SetRanges(TAxis* axis)
92d2d8ad 36{
0ab29cfa 37 if (strcmp(axis->GetTitle(), "#eta") == 0)
69b09e3b 38 axis->SetRangeUser(-1.4999, 1.4999);
39 if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0 || strcmp(axis->GetTitle(), "p_{T} (GeV/c)") == 0)
40 {
c17301f3 41 axis->SetRangeUser(0, 4.9999);
69b09e3b 42 axis->SetTitle("p_{T} (GeV/c)");
43 }
70fdd197 44 if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0 || strcmp(axis->GetTitle(), "vtx z (cm)") == 0 || strcmp(axis->GetTitle(), "vtx-z [cm]") == 0 || strcmp(axis->GetTitle(), "vtx-z (cm)") == 0)
69b09e3b 45 {
0ab29cfa 46 axis->SetRangeUser(-15, 14.9999);
69b09e3b 47 axis->SetTitle("vtx-z (cm)");
48 }
0ab29cfa 49 if (strcmp(axis->GetTitle(), "Ntracks") == 0)
50 axis->SetRangeUser(0, 99.9999);
25db2d85 51}
52
0ab29cfa 53void SetRanges(TH1* hist)
25db2d85 54{
0ab29cfa 55 SetRanges(hist->GetXaxis());
56 SetRanges(hist->GetYaxis());
57 SetRanges(hist->GetZaxis());
58}
25db2d85 59
0ab29cfa 60
61void Prepare3DPlot(TH3* hist)
62{
63 hist->GetXaxis()->SetTitleOffset(1.5);
64 hist->GetYaxis()->SetTitleOffset(1.5);
65 hist->GetZaxis()->SetTitleOffset(1.5);
66
67 hist->SetStats(kFALSE);
68}
69
70void Prepare2DPlot(TH2* hist)
71{
72 hist->SetStats(kFALSE);
73 hist->GetYaxis()->SetTitleOffset(1.4);
74
75 hist->SetMinimum(0);
76 hist->SetMaximum(gMax);
77
78 SetRanges(hist);
79}
80
81void Prepare1DPlot(TH1* hist)
82{
83 hist->SetLineWidth(2);
84 hist->SetStats(kFALSE);
85
72e597d7 86 hist->GetXaxis()->SetLabelOffset(0.02);
87 hist->GetXaxis()->SetTitleOffset(1.3);
88 hist->GetYaxis()->SetTitleOffset(1.3);
89
0ab29cfa 90 SetRanges(hist);
91}
92
93void InitPad()
94{
95 gPad->Range(0, 0, 1, 1);
96 gPad->SetLeftMargin(0.15);
97 //gPad->SetRightMargin(0.05);
98 //gPad->SetTopMargin(0.13);
72e597d7 99 gPad->SetBottomMargin(0.12);
0ab29cfa 100
72e597d7 101 gPad->SetGridx();
102 gPad->SetGridy();
0ab29cfa 103}
104
105void InitPadCOLZ()
106{
107 gPad->Range(0, 0, 1, 1);
108 gPad->SetRightMargin(0.15);
109 gPad->SetLeftMargin(0.12);
c17301f3 110 gPad->SetTopMargin(0.05);
0ab29cfa 111
112 gPad->SetGridx();
113 gPad->SetGridy();
d09fb536 114}
115
dd367a14 116// --- end of helpers --- begin functions ---
117
3dfa46a4 118void DrawOverview(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
dd367a14 119{
120 loadlibs();
121 if (!TFile::Open(fileName))
122 return;
123
124 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
125 if (!dNdEtaCorrection->LoadHistograms())
126 return;
127
128 dNdEtaCorrection->Finish();
129
130 dNdEtaCorrection->DrawOverview();
131}
132
51f6de65 133void PrintInfo(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
134{
135 loadlibs();
136 if (!TFile::Open(fileName))
137 return;
138
139 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
140 if (!dNdEtaCorrection->LoadHistograms())
141 return;
142
143 dNdEtaCorrection->Finish();
144
145 for (Int_t i=AlidNdEtaCorrection::kTrack2Particle; i<=AlidNdEtaCorrection::kND; i++)
146 {
147 Printf("Correction %d", i);
70fdd197 148 dNdEtaCorrection->GetCorrection(i)->PrintInfo(0.2);
51f6de65 149 }
150}
151
152void PrintAllInfos()
153{
154 PrintInfo();
155
156 Printf("RAW ESD");
157 TFile::Open("analysis_esd_raw.root");
158 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
159 fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
160 fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
161
162 const Int_t num = 3;
163 const char* results[] = { "dndeta", "dndetaTr", "dndetaTrVtx" };
164
165 TFile::Open("analysis_esd.root");
166 for (Int_t i=0; i<num; i++)
167 {
168 Printf("CORRECTED %s", results[i]);
169 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
170 fdNdEtaAnalysis->LoadHistograms(results[i]);
171 fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
172 }
173}
174
8ca1a6d9 175void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
176{
177 gSystem->Load("libPWG0base");
178
179 TFile::Open(esdFile);
180 dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
181 fdNdEtaAnalysisESD->LoadHistograms();
182
183 TFile::Open(mcFile);
184 dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
185 fdNdEtaAnalysisMC->LoadHistograms();
0448e811 186 //fdNdEtaAnalysisMC->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
8ca1a6d9 187
188 for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
0448e811 189 fdNdEtaAnalysisESD->GetdNdEtaPtCutOffCorrectedHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(i));
8ca1a6d9 190
191 fdNdEtaAnalysisESD->DrawHistograms();
192}
193
194void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
195{
196 gSystem->Load("libPWG0base");
197
198 const char* ESDfolder = 0;
199
200 if (plot == 0) // all
201 ESDfolder = "dndeta";
202 else if (plot == 1) // mb
203 ESDfolder = "dndeta_mb";
204 else if (plot == 2) // mb vtx
205 ESDfolder = "dndeta_mbvtx";
206
207 TFile::Open("analysis_esd.root");
208 dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
209 fdNdEtaAnalysisESD->LoadHistograms();
210
211 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
212 dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
213
214 TH2F* hist = 0;
215
216 if (plot == 0) // all
217 hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetGeneratedHistogram();
218 else if (plot == 1) // mb
219 hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetMeasuredHistogram();
220 else if (plot == 2) // mb vtx
221 hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetMeasuredHistogram();
222
223 TH1* proj = hist->ProjectionX();
224
225 TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
226 for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
227 {
228 Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
229 if (value != 0)
230 {
231 printf("vtx = %f, esd = %f, corr = %f, ratio = %f\n", vertex->GetBinCenter(i), vertex->GetBinContent(i), value, vertex->GetBinContent(i) / value);
232 vertex->SetBinContent(i, vertex->GetBinContent(i) / value);
233 }
234 }
235
236 new TCanvas;
237 vertex->DrawCopy();
238}
239
240void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
241{
a7f69e56 242 loadlibs();
8ca1a6d9 243
244 TFile::Open("analysis_esd.root");
245 dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
246 fdNdEtaAnalysisESD->LoadHistograms();
247
248 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
249 dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
250
251 //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
252 //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
253
254 TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
255 TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
256
257 TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
258
259 new TCanvas;
260 histESD->Draw();
261
262 new TCanvas;
263 histCorr->Draw();
264
265 for (Int_t x=1; x<=histESD->GetNbinsX(); ++x)
266 for (Int_t y=1; y<=histESD->GetNbinsY(); ++y)
267 for (Int_t z=1; z<=histESD->GetNbinsZ(); ++z)
268 {
269 Float_t value1 = histESD->GetBinContent(x, y, z);
270 Float_t value2 = histCorr->GetBinContent(histCorr->FindBin(histESD->GetXaxis()->GetBinCenter(x), histESD->GetYaxis()->GetBinCenter(y), histESD->GetZaxis()->GetBinCenter(z)));
271
272 if (value2 > 0 && value1 > 0)
273 {
274 printf("%f %f %f\n", value1, value2, value1 / value2);
275 diff->Fill(value1 / value2);
276 }
277 }
278
279 new TCanvas;
280 diff->Draw();
281}
282
dd367a14 283Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc = "")
284{
285 Double_t avgMC = 0;
286 Double_t avgESD = 0;
287 for (Int_t bin = histMC->FindBin(-0.7999); bin <= histMC->FindBin(0.7999); bin++)
288 {
289 avgMC += histMC->GetBinContent(bin);
290 avgESD += histESD->GetBinContent(bin);
291 }
292 Int_t nBins = histMC->FindBin(0.7999) - histMC->FindBin(-0.7999) + 1;
293
294 avgMC /= nBins;
295 avgESD /= nBins;
296
297 // deviation when integrate in |eta| < 0.8 between mc and esd
298 Double_t diffFullRange = (avgMC - avgESD) / avgMC;
299
300 Printf("%s: Integrated deviation in |eta| < 0.8 is %.2f %%", desc, diffFullRange * 100);
301
302 return diffFullRange;
303}
304
1cbdb1a6 305void dNdEtaNoResolution()
306{
307 loadlibs();
308
309 TFile::Open("correction_map.root");
310
311 const char* correctionMapFolder = "dndeta_correction";
312 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
313 dNdEtaCorrection->LoadHistograms();
314 dNdEtaCorrection->GetTrack2ParticleCorrection()->PrintInfo(0.3);
315
316 TFile::Open("analysis_mc.root");
317 fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
318 fdNdEtaAnalysis->LoadHistograms();
319 fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, "ESD (no resolution effect) -> MB with vertex");
320 fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0)->SetMarkerStyle(21);
321
322 TFile::Open("analysis_mc.root");
323 dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
324 fdNdEtaAnalysisMC->LoadHistograms();
325 fdNdEtaAnalysisMC->Finish(0, 0, AlidNdEtaCorrection::kNone, "MC: MB with vertex");
326
327 DrawdNdEtaRatio(fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0), fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(0), "MB with vertex (no resolution effect)", 3);
328}
329
330TH1* GetMCHist(const char* folder, Float_t ptCut, const char* tag)
331{
69b09e3b 332 loadlibs();
333
1cbdb1a6 334 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(folder, folder);
335 fdNdEtaAnalysis->LoadHistograms();
336 fdNdEtaAnalysis->Finish(0, ptCut, AlidNdEtaCorrection::kNone, tag);
337 return fdNdEtaAnalysis->GetdNdEtaHistogram(0);
338}
339
69b09e3b 340void dNdEtaFinal(Bool_t spd = kTRUE)
341{
342 TFile* file = TFile::Open("analysis_esd.root");
343 TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
344 TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
345 Prepare1DPlot(histESD);
346 Prepare1DPlot(histESDnsd);
347
348 TCanvas* canvas = new TCanvas("dNdEtaFinal", "dNdEtaFinal", 600, 600);
349 gPad->SetTopMargin(0.05);
350 gPad->SetRightMargin(0.05);
351 gPad->SetLeftMargin(0.12);
352 gPad->SetBottomMargin(0.12);
353 gPad->SetGridx();
354 gPad->SetGridy();
355
356 Float_t etaMax = 1.9;
357 Float_t histMax = 1.39;
358 Float_t systErrorValue = 0.023;
359 Float_t systErrorNSDValue = 0.081;
360 if (!spd)
361 {
362 //etaMax = 1.5;
363 histMax = 0.99;
364 systErrorValue = 0.043;
365 systErrorNSDValue = 0.088;
366 }
367
368 dummy = new TH2F("dummy", ";#eta;dN_{ch}/d#eta", 100, -etaMax, etaMax, 100, 3, 8);
369 dummy->SetStats(0);
370 dummy->GetYaxis()->SetTitleOffset(1.3);
371
372 histESD->SetMarkerStyle(20);
373 histESDnsd->SetMarkerStyle(21);
374 histESDnsd->SetMarkerColor(4);
375 histESDnsd->SetLineColor(4);
376 histESD->SetMarkerSize(1.5);
377 histESDnsd->SetMarkerSize(1.5);
378
379 histESD->GetXaxis()->SetRangeUser(-histMax, histMax);
380 histESDnsd->GetXaxis()->SetRangeUser(-histMax, histMax);
381
382 legend = new TLegend(0.3, 0.2, 0.78, 0.4);
383 legend->SetFillColor(0);
384 legend->SetTextSize(0.04);
385 legend->AddEntry(histESD, "Inelastic events", "P");
386 legend->AddEntry(histESDnsd, "NSD events", "P");
387
388 dummy->Draw();
389
390 // syst errors.
391 TH1* systError = (TH1*) histESD->Clone("systError");
392 for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
393 systError->SetBinError(i, systError->GetBinContent(i) * systErrorValue);
394 // change error drawing style
395 systError->SetFillColor(15);
396 systError->DrawCopy("SAME E2 ][");
397
398 // syst error NSD
399 for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
400 {
401 systError->SetBinContent(i, histESDnsd->GetBinContent(i));
402 systError->SetBinError(i, systError->GetBinContent(i) * systErrorNSDValue);
403 }
404 systError->DrawCopy("SAME E2 ][");
405
406 histESD->Draw("SAME");
407 histESDnsd->Draw("SAME");
408 legend->Draw();
409
410 canvas->SaveAs(Form("%s_dndeta_final.eps", (spd) ? "spd" : "tpc"));
411}
412
413void dNdEtaPythiaPhojet()
414{
415 // evtl. deactivate acceptance maps in dNdEtaAnalysis.cxx
416
417 loadlibs();
418
419 TH1* hist[4];
420
421 TFile::Open("LHC08c11_10TeV_0.5T/mb1/spd/analysis_mc.root");
422 hist[0] = (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
423 hist[1] = (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
424
425 TFile::Open("LHC08c15_10TeV_0.5T_Phojet/mb1/spd/analysis_mc.root");
426 hist[2] = (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMCPhojet");
427 hist[3] = (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsdPhojet");
428
429 file = TFile::Open("pythia_phojet_dndeta.root", "RECREATE");
430 for (Int_t i=0; i<4; i++)
431 hist[i]->Write();
432 file->Close();
433}
434
3dfa46a4 435void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
d09fb536 436{
69b09e3b 437 loadlibs();
438
d09fb536 439 TFile* file = TFile::Open("analysis_esd.root");
69b09e3b 440
441 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
442 fdNdEtaAnalysis->LoadHistograms("dndeta");
443
74fd10b3 444 TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
81be4ee8 445 TH1* histESD1 = (TH1*) file->Get("dndeta/dNdEta_corrected_1");
446 TH1* histESD2 = (TH1*) file->Get("dndeta/dNdEta_corrected_2");
0fc41645 447 TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
69b09e3b 448 TH1* histESDnsdNoPt = (TH1*) file->Get("dndetaNSD/dNdEta");
81be4ee8 449 TH1* histESDonePart = (TH1*) file->Get("dndetaOnePart/dNdEta_corrected");
74fd10b3 450 TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
451 TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
452 TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
453 TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
454 TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
9e952c39 455 TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
d09fb536 456
457 Prepare1DPlot(histESD);
81be4ee8 458 Prepare1DPlot(histESD1);
459 Prepare1DPlot(histESD2);
0fc41645 460 Prepare1DPlot(histESDnsd);
81be4ee8 461 Prepare1DPlot(histESDonePart);
d09fb536 462 Prepare1DPlot(histESDMB);
463 Prepare1DPlot(histESDMBVtx);
464
74fd10b3 465 Prepare1DPlot(histESDNoPt);
466 Prepare1DPlot(histESDMBNoPt);
467 Prepare1DPlot(histESDMBVtxNoPt);
9e952c39 468 Prepare1DPlot(histESDMBTracksNoPt);
74fd10b3 469
470 histESD->SetLineWidth(0);
0fc41645 471 histESDnsd->SetLineWidth(0);
81be4ee8 472 histESDonePart->SetLineWidth(0);
74fd10b3 473 histESDMB->SetLineWidth(0);
474 histESDMBVtx->SetLineWidth(0);
475
476 histESDNoPt->SetLineWidth(0);
477 histESDMBNoPt->SetLineWidth(0);
478 histESDMBVtxNoPt->SetLineWidth(0);
72e597d7 479
9e952c39 480 histESD->SetMarkerColor(1);
0fc41645 481 histESDnsd->SetMarkerColor(6);
81be4ee8 482 histESDonePart->SetMarkerColor(3);
9e952c39 483 histESDMB->SetMarkerColor(2);
69b09e3b 484 histESDMBVtx->SetMarkerColor(4);
72e597d7 485
3dfa46a4 486 histESD->SetLineColor(1);
0fc41645 487 histESDnsd->SetLineColor(6);
81be4ee8 488 histESDonePart->SetLineColor(3);
3dfa46a4 489 histESDMB->SetLineColor(2);
69b09e3b 490 histESDMBVtx->SetLineColor(4);
3dfa46a4 491
9e952c39 492 histESDNoPt->SetMarkerColor(1);
493 histESDMBNoPt->SetMarkerColor(2);
69b09e3b 494 histESDMBVtxNoPt->SetMarkerColor(4);
495 histESDMBTracksNoPt->SetMarkerColor(3);
74fd10b3 496
72e597d7 497 histESD->SetMarkerStyle(20);
0fc41645 498 histESDnsd->SetMarkerStyle(29);
81be4ee8 499 histESDonePart->SetMarkerStyle(24);
72e597d7 500 histESDMB->SetMarkerStyle(21);
501 histESDMBVtx->SetMarkerStyle(22);
d09fb536 502
74fd10b3 503 histESDNoPt->SetMarkerStyle(20);
504 histESDMBNoPt->SetMarkerStyle(21);
505 histESDMBVtxNoPt->SetMarkerStyle(22);
9e952c39 506 histESDMBTracksNoPt->SetMarkerStyle(23);
3dfa46a4 507
81be4ee8 508 Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 0.89 : 1.99;
a7f69e56 509 Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 1.2 : 2.3;
69b09e3b 510 //Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 0.89 : 1.39;
511 //Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 1.2 : 1.9;
4c351225 512
513 histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
514 histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
515 histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
0fc41645 516 histESDnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
81be4ee8 517 histESDonePart->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
74fd10b3 518
4c351225 519 histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
74fd10b3 520 histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
521 histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
9e952c39 522 histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
523
0fc41645 524 Float_t max = TMath::Max(histESDMBVtx->GetMaximum(), histESDMB->GetMaximum());
525 max = TMath::Max(max, histESD->GetMaximum());
a7f69e56 526 max = TMath::Max(max, histESDMBTracksNoPt->GetMaximum());
0fc41645 527
69b09e3b 528 TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.4);
529 legend->SetFillColor(0);
530 legend->AddEntry(histESDMBVtx, "Triggered, vertex");
531 legend->AddEntry(histESDMB, "Triggered");
1d107532 532 legend->AddEntry(histESD, "All INEL events");
533 legend->AddEntry(histESDnsd, "All NSD events");
81be4ee8 534 legend->AddEntry(histESDonePart, "One Particle");
69b09e3b 535
70fdd197 536 TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 0, max * 1.1);
537 dummy->GetYaxis()->SetRangeUser(2.1, max * 1.1);
0fc41645 538 Prepare1DPlot(dummy);
539 dummy->SetStats(kFALSE);
540 dummy->SetXTitle("#eta");
541 dummy->SetYTitle("dN_{ch}/d#eta");
542 dummy->GetYaxis()->SetTitleOffset(1);
543
0f67a57c 544 TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
d09fb536 545
0ab29cfa 546 dummy->DrawCopy();
d09fb536 547 histESDMBVtx->Draw("SAME");
548 histESDMB->Draw("SAME");
549 histESD->Draw("SAME");
1d107532 550 histESDnsd->Draw("SAME");
81be4ee8 551 histESDonePart->Draw("SAME");
69b09e3b 552 legend->Draw();
81be4ee8 553
3dfa46a4 554 if (save)
555 {
81be4ee8 556 canvas->SaveAs("dNdEta1.png");
557 //canvas->SaveAs("dNdEta1.eps");
3dfa46a4 558 }
1d107532 559
81be4ee8 560 histESD->Fit("pol0", "0", "", -0.49, 0.49);
561 histESDnsd->Fit("pol0", "0", "", -0.49, 0.49);
562 histESDonePart->Fit("pol0", "0", "", -0.49, 0.49);
563 histESDonePart->Fit("pol0", "0", "", -0.99, 0.99);
564
565 canvas = new TCanvas("dNdEta1_mirrored", "dNdEta1_mirrored", 500, 500);
566 canvas->SetGridx();
567 canvas->SetGridy();
568
569 dummy->DrawCopy()->GetXaxis()->SetRangeUser(0, 100);
570 histESD->DrawCopy("SAME")->SetMarkerStyle(24);
571 histESDnsd->DrawCopy("SAME")->SetMarkerStyle(24);
572
573 graph = new TGraphErrors(histESD);
574 for (Int_t i=0; i<graph->GetN(); i++)
575 graph->GetX()[i] *= -1;
576 graph->SetMarkerStyle(5);
577 graph->Draw("P SAME");
578
579 graph = new TGraphErrors(histESDnsd);
580 for (Int_t i=0; i<graph->GetN(); i++)
581 graph->GetX()[i] *= -1;
582 graph->SetMarkerStyle(5);
583 graph->SetMarkerColor(histESDnsd->GetMarkerColor());
584 graph->Draw("P SAME");
585
586 canvas = new TCanvas("dNdEta1_ratio", "dNdEta1_ratio", 500, 500);
587 canvas->SetGridx();
588 canvas->SetGridy();
589
590 dummy_clone = dummy->DrawCopy();
591 dummy_clone->GetXaxis()->SetRangeUser(0, 100);
592 dummy_clone->GetYaxis()->SetRangeUser(0.5, 1.5);
593
594 graph = new TGraphErrors(histESD);
595 for (Int_t i=0; i<graph->GetN(); i++)
596 {
597 Int_t bin = histESD->GetXaxis()->FindBin(-graph->GetX()[i]);
598 if (histESD->GetBinContent(bin) > 0 && graph->GetY()[i] > 0)
599 {
600 graph->GetEY()[i] = sqrt(graph->GetEY()[i] * graph->GetEY()[i] / graph->GetY()[i] / graph->GetY()[i]
601 + histESD->GetBinError(bin) * histESD->GetBinError(bin) / histESD->GetBinContent(bin) / histESD->GetBinContent(bin));
602 graph->GetY()[i] /= histESD->GetBinContent(bin);
603 graph->GetEY()[i] *= graph->GetY()[i];
604 }
605 else
606 graph->GetY()[i] = 0;
607 }
608 graph->SetMarkerStyle(5);
609 graph->Draw("P SAME");
610
611 graph = new TGraphErrors(histESDnsd);
612 for (Int_t i=0; i<graph->GetN(); i++)
613 {
614 Int_t bin = histESDnsd->GetXaxis()->FindBin(-graph->GetX()[i]);
615 if (histESDnsd->GetBinContent(bin) > 0 && graph->GetY()[i] > 0)
616 {
617 graph->GetEY()[i] = sqrt(graph->GetEY()[i] * graph->GetEY()[i] / graph->GetY()[i] / graph->GetY()[i]
618 + histESDnsd->GetBinError(bin) * histESDnsd->GetBinError(bin) / histESDnsd->GetBinContent(bin) / histESDnsd->GetBinContent(bin));
619 graph->GetY()[i] /= histESDnsd->GetBinContent(bin);
620 graph->GetEY()[i] *= graph->GetY()[i];
621 graph->GetY()[i] += 0.2;
622 }
623 }
624 graph->SetMarkerStyle(5);
625 graph->SetMarkerColor(histESDnsd->GetMarkerColor());
626 graph->Draw("P SAME");
627
628 canvas = new TCanvas("dNdEta1_vertex", "dNdEta1_vertex", 500, 500);
629 dummy->DrawCopy();
630 histESD->DrawCopy("SAME");
631 histESD1->SetLineColor(2);
632 histESD1->DrawCopy("SAME");
633 histESD2->SetLineColor(4);
634 histESD2->DrawCopy("SAME");
635
0ab29cfa 636 if (onlyESD)
637 return;
d09fb536 638
dd367a14 639 loadlibs();
74fd10b3 640
d09fb536 641 TFile* file2 = TFile::Open("analysis_mc.root");
745d6088 642
81be4ee8 643 TH1* histMCTrVtx = (TH1*) GetMCHist("dndetaTrVtx", -1, "MC: MB with vertex")->Clone("histMCTrVtx");
69b09e3b 644 TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
645
1cbdb1a6 646 TH1* histMC = (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
647 TH1* histMCTr = (TH1*) GetMCHist("dndetaTr", -1, "MC: minimum bias")->Clone("histMCTr");
1cbdb1a6 648 TH1* histMCnsd = (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
81be4ee8 649 TH1* histMConePart = (TH1*) GetMCHist("dndetaOnePart", -1, "MC: OnePart")->Clone("histMConePart");
d09fb536 650
81be4ee8 651 TH1* histMCPtCut = (TH1*) GetMCHist("dndeta", 0.151, "MC: full inelastic, pt cut")->Clone("histMCPtCut");
652 TH1* histMCTrPtCut = (TH1*) GetMCHist("dndetaTr", 0.151, "MC: minimum bias, pt cut")->Clone("histMCTrPtCut");
653 TH1* histMCTrVtxPtCut = (TH1*) GetMCHist("dndetaTrVtx", 0.151, "MC: MB with vertex, pt cut")->Clone("histMCTrVtxPtCut");
654 TH1* histMCnsdNoPt = (TH1*) GetMCHist("dndetaNSD", 0.151, "MC: NSD, put cut")->Clone("histMCnsdNoPt");
655 TH1* histMCTracksPtCut = (TH1*) GetMCHist("dndetaTracks", 0.151, "MC: Tracks w/o resolution effect, pt cut")->Clone("histMCTracksPtCut");
9e952c39 656
d09fb536 657 Prepare1DPlot(histMC);
0fc41645 658 Prepare1DPlot(histMCnsd);
74fd10b3 659 Prepare1DPlot(histMCTr);
660 Prepare1DPlot(histMCTrVtx);
661
d09fb536 662 Prepare1DPlot(histMCPtCut);
74fd10b3 663 Prepare1DPlot(histMCTrPtCut);
664 Prepare1DPlot(histMCTrVtxPtCut);
745d6088 665 Prepare1DPlot(histMCTracksPtCut);
666
667 histMC->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
668 histMCnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
669 histMCTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
670 histMCTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
671
672 histMCPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
673 histMCTrPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
674 histMCTrVtxPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
675 histMCTracksPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
d09fb536 676
9e952c39 677 histMC->SetLineColor(1);
0fc41645 678 histMCnsd->SetLineColor(6);
9e952c39 679 histMCTr->SetLineColor(2);
69b09e3b 680 histMCTrVtx->SetLineColor(4);
74fd10b3 681
9e952c39 682 histMCPtCut->SetLineColor(1);
683 histMCTrPtCut->SetLineColor(2);
69b09e3b 684 histMCTrVtxPtCut->SetLineColor(4);
9e952c39 685 if (histMCTracksPtCut)
69b09e3b 686 histMCTracksPtCut->SetLineColor(3);
d09fb536 687
745d6088 688 TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
689
0ab29cfa 690 TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
0fc41645 691 dummy2->GetYaxis()->SetRangeUser(0, max * 1.1);
d09fb536 692
0ab29cfa 693 dummy2->DrawCopy();
d09fb536 694 histMC->Draw("SAME");
0fc41645 695 histMCnsd->Draw("SAME");
74fd10b3 696 histMCTr->Draw("SAME");
697 histMCTrVtx->Draw("SAME");
d09fb536 698 histESD->Draw("SAME");
0fc41645 699 histESDnsd->Draw("SAME");
74fd10b3 700 histESDMB->Draw("SAME");
701 histESDMBVtx->Draw("SAME");
d09fb536 702 histESDNoPt->Draw("SAME");
74fd10b3 703 histESDMBNoPt->Draw("SAME");
704 histESDMBVtxNoPt->Draw("SAME");
9e952c39 705 histESDMBTracksNoPt->Draw("SAME");
d09fb536 706 histMCPtCut->Draw("SAME");
74fd10b3 707 histMCTrPtCut->Draw("SAME");
708 histMCTrVtxPtCut->Draw("SAME");
9e952c39 709 if (histMCTracksPtCut)
710 histMCTracksPtCut->Draw("SAME");
d09fb536 711
3dfa46a4 712 if (save)
713 {
714 canvas2->SaveAs("dNdEta2.gif");
715 canvas2->SaveAs("dNdEta2.eps");
716 }
0ab29cfa 717
69b09e3b 718 TH1* ratio = (TH1*) DrawdNdEtaRatio(histESD, histMC, "full_inelastic", etaPlotLimit)->Clone();
719 TH1* ratioTr = (TH1*) DrawdNdEtaRatio(histESDMB, histMCTr, "triggered", etaPlotLimit)->Clone();
720 TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
721 TH1* ratioTrVtxNoPt = (TH1*) DrawdNdEtaRatio(histESDMBVtxNoPt, histMCTrVtxPtCut, "triggered_vertex_nopt", etaPlotLimit)->Clone();
722 TH1* ratioNSD = (TH1*) DrawdNdEtaRatio(histESDnsd, histMCnsd, "NSD", etaPlotLimit)->Clone();
81be4ee8 723 TH1* ratioOnePart = (TH1*) DrawdNdEtaRatio(histESDonePart, histMConePart, "OnePart", etaPlotLimit)->Clone();
69b09e3b 724
725 // draw ratios of single steps
726 c7 = new TCanvas("all_ratios", "all_ratios", 600, 600);
727 c7->SetRightMargin(0.05);
728 c7->SetTopMargin(0.05);
729 c7->SetGridx();
730 c7->SetGridy();
731
732 ratioTrVtxNoPt->SetMarkerStyle(20);
733 ratioTrVtx->SetMarkerStyle(21);
734 ratioTr->SetMarkerStyle(23);
735 ratio->SetMarkerStyle(22);
736 ratioNSD->SetMarkerStyle(26);
737
738 ratioTrVtxNoPt->SetMarkerSize(2);
739 ratioTrVtx->SetMarkerSize(2);
740 ratioTr->SetMarkerSize(2);
741 ratio->SetMarkerSize(2);
742 ratioNSD->SetMarkerSize(2);
743
744 ratioTrVtxNoPt->SetMarkerColor(1);
745 ratioTrVtx->SetMarkerColor(2);
746 ratioTr->SetMarkerColor(4);
747 ratio->SetMarkerColor(2);
748 ratioNSD->SetMarkerColor(1);
749
750 ratioTrVtxNoPt->SetLineColor(1);
751 ratioTrVtx->SetLineColor(2);
752 ratioTr->SetLineColor(4);
753 ratio->SetLineColor(2);
754 ratioNSD->SetLineColor(1);
755
756 legend7 = new TLegend(0.13, 0.7, 0.94, 0.9);
757 legend7->SetFillColor(0);
758 legend7->SetTextSize(0.035);
759 legend7->SetNColumns(2);
760
761 flat = new TF1("flat", "-1", -5, 5);
762 ratioTrVtxNoPt->Add(flat);
763 ratioTrVtx->Add(flat);
764 ratioTr->Add(flat);
765 ratio->Add(flat);
766 ratioNSD->Add(flat);
767
768 ratioTrVtxNoPt->Scale(100);
769 ratioTrVtx->Scale(100);
770 ratioTr->Scale(100);
771 ratio->Scale(100);
772 ratioNSD->Scale(100);
773
774 ratio->Add(ratioTr, -1);
775 ratioNSD->Add(ratioTr, -1);
776 ratioTr->Add(ratioTrVtx, -1);
777 ratioTrVtx->Add(ratioTrVtxNoPt, -1);
778
779 legend7->AddEntry(ratioTrVtxNoPt, "Track-to-particle", "P");
780 legend7->AddEntry(ratio, "Trigger-bias INEL", "P");
781 legend7->AddEntry(ratioTr, "Vertex-reconstruction", "P");
782 legend7->AddEntry(ratioNSD, "Trigger-bias NSD", "P");
a7f69e56 783 if ((fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kFieldOn) && (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC || fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPCITS))
69b09e3b 784 legend7->AddEntry(ratioTrVtx, "p_{T} cut-off", "P");
785
786 TH1* dummy7 = new TH2F("dummy7", ";#eta;Deviation in %", 100, -etaPlotLimit, etaPlotLimit, 100, -5, 7);
787 dummy7->SetStats(0);
788 dummy7->Draw();
789
790 ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
791 ratioTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
792 ratioTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
793 ratioTrVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
794 ratioNSD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
795
796 ratio->Draw("HIST EP SAME");
797 ratioTr->Draw("HIST EP SAME");
a7f69e56 798 if ((fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kFieldOn) && (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC || fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPCITS))
69b09e3b 799 ratioTrVtx->Draw("HIST EP SAME");
800 ratioTrVtxNoPt->Draw("HIST EP SAME");
801 ratioNSD->Draw("HIST EP SAME");
802 legend7->Draw();
803
81be4ee8 804 //c7->SaveAs("ratios.eps");
745d6088 805
0fc41645 806 new TCanvas;
807 dummy2->DrawCopy();
808 histMCnsd->Draw("SAME");
809 histESDnsd->Draw("SAME");
810
69b09e3b 811 ratio = (TH1*) histMC->Clone("ratio");
8ca1a6d9 812 TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
813
814 ratio->Divide(histESD);
815 ratioNoPt->Divide(histESDNoPt);
816
3dfa46a4 817 ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
8ca1a6d9 818
819 ratio->SetLineColor(1);
820 ratioNoPt->SetLineColor(2);
821
dd367a14 822 Double_t average = 0; // average deviation from 1 in ratio (depends on the number of bins if statistical)
823 for (Int_t bin = ratio->FindBin(-0.7999); bin <= ratio->FindBin(0.7999); bin++)
824 average += TMath::Abs(ratio->GetBinContent(bin) - 1);
825 Int_t nBins = ratio->FindBin(0.7999) - ratio->FindBin(-0.7999) + 1;
826 average /= nBins;
827 Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
828
829 PrintIntegratedDeviation(histMC, histESD, "all events");
51f6de65 830 PrintIntegratedDeviation(histMCnsd, histESDnsd, "all events (NSD)");
81be4ee8 831 PrintIntegratedDeviation(histMConePart, histESDonePart, "all events (INEL>0)");
dd367a14 832 PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
833 PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
834 PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
69b09e3b 835 PrintIntegratedDeviation(histMCnsdNoPt, histESDnsdNoPt, "all events (NSD) (no pt corr)");
dd367a14 836 PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
837 PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
a7f69e56 838 PrintIntegratedDeviation(histESD, histESDNoPt, "pt cut off correction");
dd367a14 839
69b09e3b 840 TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 600, 600);
8ca1a6d9 841 canvas3->Range(0, 0, 1, 1);
842 //canvas3->Divide(1, 2, 0, 0);
843
844 //canvas3->cd(1);
9e952c39 845 TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.5, 0.98, 0.98);
69b09e3b 846 pad1->SetTopMargin(0.05);
847 pad1->SetLeftMargin(0.13);
8ca1a6d9 848 pad1->Draw();
849
9e952c39 850 TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.5);
69b09e3b 851 pad2->SetLeftMargin(0.13);
8ca1a6d9 852 pad2->Draw();
853
69b09e3b 854 pad1->SetRightMargin(0.01);
855 pad2->SetRightMargin(0.01);
0ab29cfa 856
8ca1a6d9 857 // no border between them
858 pad1->SetBottomMargin(0);
859 pad2->SetTopMargin(0);
860
861 pad1->cd();
69b09e3b 862 pad1->SetGridx();
863 pad1->SetGridy();
8ca1a6d9 864
72e597d7 865 legend->AddEntry(histMC, "MC prediction");
0ab29cfa 866
69b09e3b 867 dummy->GetXaxis()->SetLabelSize(0.08);
868 dummy->GetYaxis()->SetLabelSize(0.08);
869 dummy->GetXaxis()->SetTitleSize(0.08);
870 dummy->GetYaxis()->SetTitleSize(0.08);
871 dummy->GetYaxis()->SetTitleOffset(0.8);
72e597d7 872 dummy->DrawCopy();
0ab29cfa 873 histESDMBVtx->Draw("SAME");
874 histESDMB->Draw("SAME");
875 histESD->Draw("SAME");
0ab29cfa 876 histMC->Draw("SAME");
72e597d7 877
69b09e3b 878 legend->SetTextSize(0.08);
72e597d7 879 legend->Draw();
0ab29cfa 880
8ca1a6d9 881 pad2->cd();
882 pad2->SetBottomMargin(0.15);
69b09e3b 883 //pad2->SetGridx();
884 //pad2->SetGridy();
72e597d7 885
69b09e3b 886 Float_t minR = 0.91; //TMath::Min(0.961, ratio->GetMinimum() * 0.95);
887 Float_t maxR = 1.09; //TMath::Max(1.049, ratio->GetMaximum() * 1.05);
9e952c39 888
69b09e3b 889 TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
8ca1a6d9 890 dummy3.SetStats(kFALSE);
745d6088 891 for (Int_t i=1; i<=100; ++i)
892 dummy3.SetBinContent(i, 1);
0fc41645 893 dummy3.GetYaxis()->SetRangeUser(minR, maxR);
8ca1a6d9 894 dummy3.SetLineWidth(2);
69b09e3b 895 dummy3.GetXaxis()->SetLabelSize(0.08);
896 dummy3.GetYaxis()->SetLabelSize(0.08);
897 dummy3.GetXaxis()->SetTitleSize(0.08);
898 dummy3.GetYaxis()->SetTitleSize(0.08);
899 dummy3.GetYaxis()->SetTitleOffset(0.8);
8ca1a6d9 900 dummy3.DrawCopy();
72e597d7 901
8ca1a6d9 902 ratio->Draw("SAME");
72e597d7 903
8ca1a6d9 904 //pad2->Draw();
72e597d7 905
8ca1a6d9 906 canvas3->Modified();
72e597d7 907
3dfa46a4 908 if (save)
909 {
910 canvas3->SaveAs("dNdEta.gif");
911 canvas3->SaveAs("dNdEta.eps");
912 }
8ca1a6d9 913
914 TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
72e597d7 915
916 ratio->Draw();
917 ratioNoPt->Draw("SAME");
918
4c351225 919 TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
920 legend->SetFillColor(0);
921 legend->AddEntry(ratio, "mc/esd");
922 legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
923 legend->Draw();
d09fb536 924}
925
81be4ee8 926void CompareTwodNdEta(const char* fileName1, const char* fileName2, Bool_t errorsCorrelated = kFALSE)
927{
928 c = new TCanvas;
929
930 c->SetGridx();
931 c->SetGridy();
932
933 hist = new TH2F("dummy", ";#eta;dN_{ch}/d#eta", 100, -2.5, 2.5, 100, 0, 8);
934 hist->SetStats(0);
935 hist->DrawCopy();//->GetYaxis()->SetRangeUser(2, 4.5);
936
937 l = new TLegend(0.2, 0.13, 0.8, 0.35);
938 l->SetNColumns(2);
939 l->SetFillColor(0);
940
941 TH1* histESD[2];
942 TH1* histESDnsd[2];
943
944 for (Int_t i=0; i<2; i++)
945 {
946 if (i == 0)
947 file = TFile::Open(fileName1);
948 if (i == 1)
949 {
950 if (fileName2 == 0)
951 break;
952 file = TFile::Open(fileName2);
953 }
954
955 histESD[i] = (TH1*) file->Get("dndeta/dNdEta_corrected");
956 histESDnsd[i] = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
957
958 histESD[i]->SetMarkerStyle(20 + i*4);
959 histESDnsd[i]->SetMarkerStyle(21 + i*4);
960
961 histESD[i]->SetMarkerColor(i+1);
962 histESD[i]->SetLineColor(i+1);
963 histESDnsd[i]->SetMarkerColor(i+1);
964 histESDnsd[i]->SetLineColor(i+1);
965
966 histESD[i]->DrawCopy("SAME");
967 histESDnsd[i]->DrawCopy("SAME");
968
969 l->AddEntry(histESD[i], Form("Data %d INEL", i), "P");
970 l->AddEntry(histESDnsd[i], Form("Data %d NSD", i), "P");
971 }
972
973 if (0)
974 {
975 TGraphErrors *gre = new TGraphErrors(16);
976 gre->SetFillColor(4);
977 gre->SetMarkerColor(4);
978 gre->SetMarkerStyle(26);
979 gre->SetPoint(0,0.125,3.14);
980 gre->SetPointError(0,0,0.07);
981 gre->SetPoint(1,0.375,3.04);
982 gre->SetPointError(1,0,0.07);
983 gre->SetPoint(2,0.625,3.17);
984 gre->SetPointError(2,0,0.07);
985 gre->SetPoint(3,0.875,3.33);
986 gre->SetPointError(3,0,0.07);
987 gre->SetPoint(4,1.125,3.33);
988 gre->SetPointError(4,0,0.07);
989 gre->SetPoint(5,1.375,3.53);
990 gre->SetPointError(5,0,0.07);
991 gre->SetPoint(6,1.625,3.46);
992 gre->SetPointError(6,0,0.07);
993 gre->SetPoint(7,1.875,3.41);
994 gre->SetPointError(7,0,0.07);
995 gre->SetPoint(8,-0.125,3.14);
996 gre->SetPointError(8,0,0.07);
997 gre->SetPoint(9,-0.375,3.04);
998 gre->SetPointError(9,0,0.07);
999 gre->SetPoint(10,-0.625,3.17);
1000 gre->SetPointError(10,0,0.07);
1001 gre->SetPoint(11,-0.875,3.33);
1002 gre->SetPointError(11,0,0.07);
1003 gre->SetPoint(12,-1.125,3.33);
1004 gre->SetPointError(12,0,0.07);
1005 gre->SetPoint(13,-1.375,3.53);
1006 gre->SetPointError(13,0,0.07);
1007 gre->SetPoint(14,-1.625,3.46);
1008 gre->SetPointError(14,0,0.07);
1009 gre->SetPoint(15,-1.875,3.41);
1010 gre->SetPointError(15,0,0.07);
1011 gre->Draw("p");
1012
1013 l->AddEntry(gre, "UA5 INEL", "P");
1014
1015 gre = new TGraphErrors(16);
1016 gre->SetMarkerColor(4);
1017 gre->SetFillColor(4);
1018 gre->SetMarkerStyle(22);
1019 gre->SetPoint(0,0.125,3.48);
1020 gre->SetPointError(0,0,0.07);
1021 gre->SetPoint(1,0.375,3.38);
1022 gre->SetPointError(1,0,0.07);
1023 gre->SetPoint(2,0.625,3.52);
1024 gre->SetPointError(2,0,0.07);
1025 gre->SetPoint(3,0.875,3.68);
1026 gre->SetPointError(3,0,0.07);
1027 gre->SetPoint(4,1.125,3.71);
1028 gre->SetPointError(4,0,0.07);
1029 gre->SetPoint(5,1.375,3.86);
1030 gre->SetPointError(5,0,0.07);
1031 gre->SetPoint(6,1.625,3.76);
1032 gre->SetPointError(6,0,0.07);
1033 gre->SetPoint(7,1.875,3.66);
1034 gre->SetPointError(7,0,0.07);
1035 gre->SetPoint(8,-0.125,3.48);
1036 gre->SetPointError(8,0,0.07);
1037 gre->SetPoint(9,-0.375,3.38);
1038 gre->SetPointError(9,0,0.07);
1039 gre->SetPoint(10,-0.625,3.52);
1040 gre->SetPointError(10,0,0.07);
1041 gre->SetPoint(11,-0.875,3.68);
1042 gre->SetPointError(11,0,0.07);
1043 gre->SetPoint(12,-1.125,3.71);
1044 gre->SetPointError(12,0,0.07);
1045 gre->SetPoint(13,-1.375,3.86);
1046 gre->SetPointError(13,0,0.07);
1047 gre->SetPoint(14,-1.625,3.76);
1048 gre->SetPointError(14,0,0.07);
1049 gre->SetPoint(15,-1.875,3.66);
1050 gre->SetPointError(15,0,0.07);
1051 gre->Draw("p");
1052
1053 l->AddEntry(gre, "UA5 NSD", "P");
1054 }
1055
1056 l->Draw();
1057
1058 if (fileName2 == 0)
1059 return;
1060
1061 new TCanvas;
1062 gPad->SetGridx();
1063 gPad->SetGridy();
1064
1065 if (errorsCorrelated)
1066 {
1067 for (Int_t i=1; i<=histESD[1]->GetNbinsX(); i++)
1068 {
1069 histESD[1]->SetBinError(i, 0);
1070 histESDnsd[1]->SetBinError(i, 0);
1071 }
1072 }
1073
1074 histESD[0]->Divide(histESD[0], histESD[1]);
1075 histESDnsd[0]->Divide(histESDnsd[0], histESDnsd[1]);
1076
1077 for (Int_t i=1; i<=histESD[1]->GetNbinsX(); i++)
1078 histESDnsd[0]->SetBinContent(i, histESDnsd[0]->GetBinContent(i) + 0.2);
1079
1080 hist->DrawCopy()->GetYaxis()->SetRangeUser(0.8, 1.4);
1081 histESD[0]->Draw("SAME");
1082 histESDnsd[0]->Draw("SAME");
1083}
1084
69b09e3b 1085TH1* DrawdNdEtaRatio(TH1* corr, TH1* mc, const char* name, Float_t etaPlotLimit)
745d6088 1086{
69b09e3b 1087 TCanvas* canvas3 = new TCanvas(name, name, 600, 600);
745d6088 1088 canvas3->Range(0, 0, 1, 1);
1089
1090 TPad* pad1 = new TPad(Form("%s_1", name), "", 0, 0.5, 0.98, 0.98);
1091 pad1->Draw();
1092
1093 TPad* pad2 = new TPad(Form("%s_2", name), "", 0, 0.02, 0.98, 0.5);
1094 pad2->Draw();
1095
69b09e3b 1096 pad1->SetRightMargin(0.01);
1097 pad2->SetRightMargin(0.01);
1098 pad1->SetTopMargin(0.05);
1099 pad1->SetLeftMargin(0.13);
1100 pad2->SetLeftMargin(0.13);
1101 pad2->SetBottomMargin(0.15);
1102
745d6088 1103 // no border between them
1104 pad1->SetBottomMargin(0);
1105 pad2->SetTopMargin(0);
1106
1107 pad1->cd();
69b09e3b 1108 pad1->SetGridx();
1109 pad1->SetGridy();
745d6088 1110
69b09e3b 1111 TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.3);
745d6088 1112 legend->SetFillColor(0);
69b09e3b 1113 legend->AddEntry(corr, "Corrected");
745d6088 1114 legend->AddEntry(mc, "MC prediction");
69b09e3b 1115 legend->SetTextSize(0.08);
745d6088 1116
a7f69e56 1117 TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 2.7, corr->GetMaximum() * 1.1);
745d6088 1118 Prepare1DPlot(dummy);
1119 dummy->SetStats(kFALSE);
1120 dummy->SetXTitle("#eta");
1121 dummy->SetYTitle("dN_{ch}/d#eta");
1122 dummy->GetYaxis()->SetTitleOffset(1);
1123
69b09e3b 1124 dummy->GetXaxis()->SetLabelSize(0.08);
1125 dummy->GetYaxis()->SetLabelSize(0.08);
1126 dummy->GetXaxis()->SetTitleSize(0.08);
1127 dummy->GetYaxis()->SetTitleSize(0.08);
1128 dummy->GetYaxis()->SetTitleOffset(0.8);
745d6088 1129 dummy->DrawCopy();
1130
1131 corr->Draw("SAME");
1132 mc->Draw("SAME");
1133
1134 legend->Draw();
1135
1136 pad2->cd();
1137 pad2->SetBottomMargin(0.15);
69b09e3b 1138 //pad2->SetGridx();
1139 //pad2->SetGridy();
745d6088 1140
1141 TH1* ratio = (TH1*) mc->Clone("ratio");
1142 ratio->Divide(corr);
1143
69b09e3b 1144 Float_t minR = TMath::Min(0.91, ratio->GetMinimum() * 0.95);
1145 Float_t maxR = TMath::Max(1.09, ratio->GetMaximum() * 1.05);
745d6088 1146
1147 TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
1148 dummy3.SetStats(kFALSE);
1149 for (Int_t i=1; i<=100; ++i)
1150 dummy3.SetBinContent(i, 1);
1151 dummy3.GetYaxis()->SetRangeUser(minR, maxR);
1152 dummy3.SetLineWidth(2);
69b09e3b 1153 dummy3.GetXaxis()->SetLabelSize(0.08);
1154 dummy3.GetYaxis()->SetLabelSize(0.08);
1155 dummy3.GetXaxis()->SetTitleSize(0.08);
1156 dummy3.GetYaxis()->SetTitleSize(0.08);
1157 dummy3.GetYaxis()->SetTitleOffset(0.8);
745d6088 1158 dummy3.DrawCopy();
1159
1160 ratio->Draw("SAME");
1161
1162 canvas3->Modified();
69b09e3b 1163
1164 return ratio;
745d6088 1165}
1166
d09fb536 1167void ptSpectrum()
1168{
1169 TFile* file = TFile::Open("analysis_esd.root");
0ab29cfa 1170 TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
d09fb536 1171
1172 TFile* file2 = TFile::Open("analysis_mc.root");
0ab29cfa 1173 TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
d09fb536 1174
1175 TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
1176 InitPad();
1177 gPad->SetLogy();
1178
1179 Prepare1DPlot(histMC);
1180 Prepare1DPlot(histESD);
1181
1182 histESD->SetTitle("");
25db2d85 1183 histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1184 histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
d09fb536 1185
1186 histMC->SetLineColor(kBlue);
1187 histESD->SetLineColor(kRed);
1188
1189 histESD->GetYaxis()->SetTitleOffset(1.5);
1190 histESD->GetXaxis()->SetRangeUser(0, 4.9999);
1191
1192 histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
1193
1194 histESD->Draw();
1195 histMC->Draw("SAME");
1196
1197 canvas->SaveAs("ptSpectrum.gif");
0ab29cfa 1198 canvas->SaveAs("ptSpectrum.eps");
92d2d8ad 1199}
1200
0bd1f8a0 1201void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
92d2d8ad 1202{
0448e811 1203 gSystem->Load("libPWG0base");
1204
1205 TFile::Open(fileName);
1206 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1207 dNdEtaCorrection->LoadHistograms();
0ab29cfa 1208
0448e811 1209 TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
1210 TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
0ab29cfa 1211
1212 Prepare2DPlot(corrTrigger);
0448e811 1213 corrTrigger->SetTitle("b) Trigger bias correction");
0ab29cfa 1214
1215 Prepare2DPlot(corrVtx);
0448e811 1216 corrVtx->SetTitle("a) Vertex reconstruction correction");
0ab29cfa 1217
5c495d37 1218 corrTrigger->GetYaxis()->SetTitle("Multiplicity");
1219 corrVtx->GetYaxis()->SetTitle("Multiplicity");
1220
0ab29cfa 1221 TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
1222 canvas->Divide(2, 1);
1223
1224 canvas->cd(1);
1225 InitPadCOLZ();
0448e811 1226 corrVtx->DrawCopy("COLZ");
0ab29cfa 1227
1228 canvas->cd(2);
1229 InitPadCOLZ();
0448e811 1230 corrTrigger->DrawCopy("COLZ");
0ab29cfa 1231
1232 canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
1233 canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
1234
1235 canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
1236 canvas->Divide(2, 1);
1237
1238 corrTrigger->GetYaxis()->SetRangeUser(0, 5);
1239 corrVtx->GetYaxis()->SetRangeUser(0, 5);
1240
1241 canvas->cd(1);
1242 InitPadCOLZ();
0448e811 1243 corrVtx->DrawCopy("COLZ");
0ab29cfa 1244
1245 canvas->cd(2);
1246 InitPadCOLZ();
0448e811 1247 corrTrigger->DrawCopy("COLZ");
0ab29cfa 1248
1249 canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
1250 canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
1251}
1252
1253void TriggerBias(const char* fileName = "correction_map.root")
1254{
1255 TFile* file = TFile::Open(fileName);
92d2d8ad 1256
5c495d37 1257 TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
92d2d8ad 1258
1259 Prepare2DPlot(corr);
1260 corr->SetTitle("Trigger bias correction");
1261
1262 TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
1263 InitPadCOLZ();
1264 corr->DrawCopy("COLZ");
1265
25db2d85 1266 canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
0ab29cfa 1267 canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
92d2d8ad 1268
1269 corr->GetYaxis()->SetRangeUser(0, 5);
1270
1271 canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
1272 InitPadCOLZ();
1273 corr->DrawCopy("COLZ");
1274
25db2d85 1275 canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
0ab29cfa 1276 canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
92d2d8ad 1277}
1278
72e597d7 1279void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
1280{
1281 gSystem->Load("libPWG0base");
1282
1283 TFile* file = TFile::Open(fileName);
1284 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
0448e811 1285 dNdEtaCorrection->LoadHistograms();
72e597d7 1286
81be4ee8 1287 TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->Get1DCorrection("x");
1288 TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->Get1DCorrection("y", -5, 5);
72e597d7 1289
0448e811 1290 TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
8ca1a6d9 1291 canvas->Divide(2, 1);
1292
1293 canvas->cd(1);
72e597d7 1294 InitPad();
1295
1296 Prepare1DPlot(hist);
1297 hist->SetTitle("");
1298 hist->GetYaxis()->SetTitle("correction factor");
1299 hist->GetYaxis()->SetRangeUser(1, 1.5);
1300 hist->GetYaxis()->SetTitleOffset(1.6);
1301 hist->Draw();
1302
8ca1a6d9 1303 canvas->cd(2);
1304 InitPad();
1305
1306 Prepare1DPlot(hist2);
1307 hist2->SetTitle("");
1308 hist2->GetYaxis()->SetTitle("correction factor");
1309 hist2->GetXaxis()->SetRangeUser(0, 5);
1310 hist2->GetYaxis()->SetTitleOffset(1.6);
1311 hist2->GetXaxis()->SetTitle("multiplicity");
1312 hist2->Draw();
1313
1314 TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
1315 pave->SetFillColor(0);
81be4ee8 1316 pave->AddText("|z| < 5 cm");
8ca1a6d9 1317 pave->Draw();
81be4ee8 1318
1319 Float_t triggerEff = 100.0 / hist2->GetBinContent(1);
1320 Printf("trigger eff in 0 bin is: %.2f %%", triggerEff);
1321
1322 return;
8ca1a6d9 1323
81be4ee8 1324 TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1325 //new TCanvas;
1326 //hist2->Draw();
1327
1328 Printf("vertex reco eff in 0 bin is: %.2f %%", 100.0 / hist2->GetBinContent(1));
1329
1330 Printf("combined efficiency is %.2f %%", triggerEff / hist2->GetBinContent(1));
72e597d7 1331}
1332
92d2d8ad 1333void VtxRecon()
1334{
1335 TFile* file = TFile::Open("correction_map.root");
1336
72e597d7 1337 TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
92d2d8ad 1338
1339 Prepare2DPlot(corr);
1340 corr->SetTitle("Vertex reconstruction correction");
1341
1342 TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
1343 InitPadCOLZ();
25db2d85 1344 corr->DrawCopy("COLZ");
1345
0ab29cfa 1346 canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
1347 canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
25db2d85 1348
1349 corr->GetYaxis()->SetRangeUser(0, 5);
1350
0ab29cfa 1351 canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
25db2d85 1352 InitPadCOLZ();
1353 corr->DrawCopy("COLZ");
92d2d8ad 1354
25db2d85 1355 canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
0ab29cfa 1356 canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
92d2d8ad 1357}
1358
72e597d7 1359void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
1360{
1361 gSystem->Load("libPWG0base");
1362
1363 TFile* file = TFile::Open(fileName);
1364 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
0448e811 1365 dNdEtaCorrection->LoadHistograms();
72e597d7 1366
0448e811 1367 TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
1368 TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
8ca1a6d9 1369
1370 TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
1371 canvas->Divide(2, 1);
72e597d7 1372
8ca1a6d9 1373 canvas->cd(1);
72e597d7 1374 InitPad();
1375
1376 Prepare1DPlot(hist);
1377 hist->SetTitle("");
1378 hist->GetYaxis()->SetTitle("correction factor");
1379 hist->GetYaxis()->SetRangeUser(1, 1.8);
1380 hist->GetYaxis()->SetTitleOffset(1.6);
8ca1a6d9 1381 hist->DrawCopy();
1382
1383 canvas->cd(2);
1384 InitPad();
1385
1386 Prepare1DPlot(hist2);
1387 hist2->SetTitle("");
1388 hist2->GetYaxis()->SetTitle("correction factor");
1389 hist2->GetXaxis()->SetRangeUser(0, 20);
1390 hist2->GetYaxis()->SetTitleOffset(1.6);
1391 hist2->GetXaxis()->SetTitle("multiplicity");
1392 hist2->Draw();
1393
1394 TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
1395 pave->SetFillColor(0);
1396 pave->AddText("|z| < 10 cm");
1397 pave->Draw();
72e597d7 1398
1399 canvas->SaveAs("VtxRecon1D.eps");
8ca1a6d9 1400
0448e811 1401 Correction1DCreatePlots(fileName, folderName, 9.9, 2);
8ca1a6d9 1402
0448e811 1403 TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
1404 TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
1405
1406 Prepare1DPlot(corrX);
1407 Prepare1DPlot(corrZ);
1408
1409 corrX->GetYaxis()->SetTitleOffset(1.5);
1410 corrZ->GetYaxis()->SetTitleOffset(1.5);
1411
1412 corrX->SetTitle("a) z projection");
1413 corrZ->SetTitle("b) p_{T} projection");
1414
69b09e3b 1415 corrX->GetYaxis()->SetTitle("Correction factor");
1416 corrZ->GetYaxis()->SetTitle("Correction factor");
0448e811 1417
1418 corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
1419
1420 TString canvasName;
1421 canvasName.Form("VtxRecon1D_Track");
1422 TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
1423 canvas->Divide(2, 1);
1424
1425 canvas->cd(1);
1426 InitPad();
1427 corrX->DrawCopy();
1428
1429 canvas->cd(2);
1430 InitPad();
1431 gPad->SetLogx();
1432 corrZ->Draw();
1433
1434 canvas->SaveAs("VtxRecon1D_Track.eps");
1435 canvas->SaveAs("VtxRecon1D_Track.gif");
72e597d7 1436}
1437
0ab29cfa 1438void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
1afae8ff 1439{
25db2d85 1440 gSystem->Load("libPWG0base");
1441
0ab29cfa 1442 TFile::Open(fileName);
8b3563f4 1443 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
0ab29cfa 1444 dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
1445
1446 TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1447 TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1448
1449 TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1450 TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1451
1452 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1453 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1454 gene->GetXaxis()->SetRangeUser(-10, 10);
1455 meas->GetXaxis()->SetRangeUser(-10, 10);
1456
1457 Float_t eff1 = gene->Integral() / meas->Integral();
1458 Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1459
1460 printf("Correction without pT cut: %f +- %f\n", eff1, error1);
1461
1462 gene->GetZaxis()->SetRangeUser(0.3, 10);
1463 meas->GetZaxis()->SetRangeUser(0.3, 10);
1464
1465 Float_t eff2 = gene->Integral() / meas->Integral();
1466 Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1467
1468 printf("Correction with pT cut: %f +- %f\n", eff2, error2);
1469
1470 gene->GetZaxis()->SetRangeUser(0.3, 1);
1471 meas->GetZaxis()->SetRangeUser(0.3, 1);
1472
1473 Float_t eff3 = gene->Integral() / meas->Integral();
1474 Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1475
1476 printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
1477}
1478
69b09e3b 1479void Correction1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType = 0, Int_t correctionType2 = -1)
0ab29cfa 1480{
69b09e3b 1481 if (correctionType2 == -1)
1482 correctionType2 = correctionType;
1483
0ab29cfa 1484 TFile::Open(fileName);
bdfe2916 1485 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
0448e811 1486 dNdEtaCorrection->LoadHistograms();
0ab29cfa 1487
0448e811 1488 TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
69b09e3b 1489 TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType2)->GetTrackCorrection()->GetMeasuredHistogram();
0ab29cfa 1490
1491 gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
1492 meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
1493 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1494 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
5a6310fe 1495 AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kFALSE);
0ab29cfa 1496 gene->GetYaxis()->SetRange(0, 0);
1497 meas->GetYaxis()->SetRange(0, 0);
1498
5a6310fe 1499 gene->GetXaxis()->SetRangeUser(-9.9, 9.9);
1500 meas->GetXaxis()->SetRangeUser(-9.9, 9.9);
1501 AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kFALSE);
0ab29cfa 1502 gene->GetZaxis()->SetRange(0, 0);
1503 meas->GetZaxis()->SetRange(0, 0);
1504
1505 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1506 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
5a6310fe 1507 AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kFALSE);
0ab29cfa 1508}
25db2d85 1509
69b09e3b 1510TCanvas* Correction1D(Int_t correctionType = 0, const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType2 = -1)
0448e811 1511{
1512 gSystem->Load("libPWG0base");
1513
69b09e3b 1514 Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType, correctionType2);
0448e811 1515
1516 TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1517 TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1518 TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1519
1520 Prepare1DPlot(corrX);
1521 Prepare1DPlot(corrY);
1522 Prepare1DPlot(corrZ);
1523
69b09e3b 1524 /*
0448e811 1525 corrX->SetTitle("a) z projection");
1526 corrY->SetTitle("b) #eta projection");
1527 corrZ->SetTitle("c) p_{T} projection");
69b09e3b 1528 */
1529
1530 corrX->SetTitle("");
1531 corrY->SetTitle("");
1532 corrZ->SetTitle("");
1533
1534 corrX->SetTitleSize(0.06, "xyz");
1535 corrX->SetLabelSize(0.06, "xyz");
1536 corrY->SetTitleSize(0.06, "xyz");
1537 corrY->SetLabelSize(0.06, "xyz");
1538 corrZ->SetTitleSize(0.06, "xyz");
1539 corrZ->SetLabelSize(0.06, "xyz");
1540
1541 corrX->GetYaxis()->SetTitle("Correction factor");
1542 corrY->GetYaxis()->SetTitle("Correction factor");
1543 corrZ->GetYaxis()->SetTitle("Correction factor");
1544 //corrX->GetYaxis()->SetTitleOffset(1.7);
1545 //corrY->GetYaxis()->SetTitleOffset(1.7);
1546 //corrZ->GetYaxis()->SetTitleOffset(1.7);
5a6310fe 1547 corrX->GetYaxis()->SetRangeUser(0.8, 1.5);
1548 corrY->GetYaxis()->SetRangeUser(0.8, 1.5);
1549 corrZ->GetYaxis()->SetRangeUser(0.8, 1.5);
0448e811 1550
5a6310fe 1551 corrZ->GetXaxis()->SetRangeUser(0.11, upperPtLimit);
0448e811 1552
1553 TString canvasName;
5a6310fe 1554 canvasName.Form(Form("Correction1D_%d_%s_%f", correctionType, fileName, upperPtLimit));
0448e811 1555 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1556 canvas->Divide(3, 1);
1557
5a6310fe 1558 TLatex* Tl = new TLatex;
69b09e3b 1559 Tl->SetTextSize(0.06);
5a6310fe 1560 Tl->SetBit(TLatex::kTextNDC);
1561
0448e811 1562 canvas->cd(1);
1563 InitPad();
69b09e3b 1564 gPad->SetTopMargin(0.05);
1565 gPad->SetBottomMargin(0.15);
0448e811 1566 corrX->DrawCopy();
69b09e3b 1567 Tl->DrawLatex(0.5, 0.88, "0.3 < p_{T} < 10");
1568 Tl->DrawLatex(0.5, 0.8, "|#eta| < 0.8");
0448e811 1569
1570 canvas->cd(2);
1571 InitPad();
69b09e3b 1572 gPad->SetTopMargin(0.05);
1573 gPad->SetBottomMargin(0.15);
0448e811 1574 corrY->Draw();
69b09e3b 1575 Tl->DrawLatex(0.5, 0.88, "0.3 < p_{T} < 10");
1576 Tl->DrawLatex(0.5, 0.8, "|vtx-z| < 10 cm");
0448e811 1577
1578 canvas->cd(3);
1579 InitPad();
69b09e3b 1580 gPad->SetTopMargin(0.05);
1581 gPad->SetBottomMargin(0.15);
5a6310fe 1582 gPad->SetLogx();
0448e811 1583 corrZ->Draw();
69b09e3b 1584 corrZ->GetXaxis()->SetLabelOffset(0.005);
1585 corrZ->GetXaxis()->SetTitleOffset(1.2);
1586 Tl->DrawLatex(0.5, 0.88, "|vtx-z| < 10 cm");
1587 Tl->DrawLatex(0.5, 0.8, "|#eta| < 0.8");
0448e811 1588
5a6310fe 1589 return canvas;
0448e811 1590}
1591
72e597d7 1592void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
0ab29cfa 1593{
1594 gSystem->Load("libPWG0base");
1595
0448e811 1596 Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
0ab29cfa 1597
0448e811 1598 TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1599 TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1600 TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
0ab29cfa 1601
1602 Prepare1DPlot(corrX);
1603 Prepare1DPlot(corrY);
1604 Prepare1DPlot(corrZ);
1605
0448e811 1606 corrX->SetTitle("a) z projection");
72e597d7 1607 corrY->SetTitle("a) #eta projection");
1608 corrZ->SetTitle("b) p_{T} projection");
1609
1610 corrY->GetYaxis()->SetTitle("correction factor");
1611 corrZ->GetYaxis()->SetTitle("correction factor");
0ab29cfa 1612
1613 corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1614
bdfe2916 1615 TString canvasName;
72e597d7 1616 canvasName.Form("Track2Particle1D_%s", folder);
bdfe2916 1617 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
0ab29cfa 1618 canvas->Divide(3, 1);
1619
1620 canvas->cd(1);
1621 InitPad();
bdfe2916 1622 corrX->DrawCopy();
0ab29cfa 1623
1624 canvas->cd(2);
1625 InitPad();
0448e811 1626 corrY->Draw();
0ab29cfa 1627
1628 canvas->cd(3);
1629 InitPad();
0448e811 1630 corrZ->Draw();
0ab29cfa 1631
72e597d7 1632 canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
1633 canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
1634
5c495d37 1635 //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
1636
72e597d7 1637 canvasName.Form("Track2Particle1D_%s_etapt", folder);
1638 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
1639 canvas->Divide(2, 1);
1640
1641 canvas->cd(1);
1642 InitPad();
1643 corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
1644 corrY->GetYaxis()->SetRangeUser(1, 1.5);
1645 corrY->GetYaxis()->SetTitleOffset(1.5);
1646 corrY->DrawCopy();
5c495d37 1647 TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
1648 pave->AddText("|z| < 10 cm");
1649 pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
1650 pave->Draw();
72e597d7 1651
1652 canvas->cd(2);
1653 InitPad();
5c495d37 1654 gPad->SetLogx();
1655 corrZ->GetYaxis()->SetRangeUser(1, 2.5);
1656 corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
72e597d7 1657 corrZ->GetYaxis()->SetTitleOffset(1.5);
1658 corrZ->DrawCopy();
5c495d37 1659 pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
1660 pave->AddText("|z| < 10 cm");
1661 pave->AddText("|#eta| < 0.8");
1662 pave->Draw();
72e597d7 1663
1664 canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
4c351225 1665 canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
0ab29cfa 1666}
1667
69b09e3b 1668/*
0ab29cfa 1669void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
1670{
1671 gSystem->Load("libPWG0base");
1672
8ca1a6d9 1673 // particle type
1674 for (Int_t particle=0; particle<4; ++particle)
1675 {
1676 TString dirName;
1677 dirName.Form("correction_%d", particle);
1678 Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
0ab29cfa 1679
8ca1a6d9 1680 TString tmpx, tmpy, tmpz;
1681 tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
1682 tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
1683 tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
0ab29cfa 1684
8ca1a6d9 1685 TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
1686 TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
1687 TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
0ab29cfa 1688
8ca1a6d9 1689 Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
0ab29cfa 1690
8ca1a6d9 1691 TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
1692 TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
1693 TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
0ab29cfa 1694
8ca1a6d9 1695 posX->Divide(negX);
1696 posY->Divide(negY);
1697 posZ->Divide(negZ);
0ab29cfa 1698
8ca1a6d9 1699 Prepare1DPlot(posX);
1700 Prepare1DPlot(posY);
1701 Prepare1DPlot(posZ);
0ab29cfa 1702
8ca1a6d9 1703 Float_t min = 0.8;
1704 Float_t max = 1.2;
0ab29cfa 1705
8ca1a6d9 1706 posX->SetMinimum(min);
1707 posX->SetMaximum(max);
1708 posY->SetMinimum(min);
1709 posY->SetMaximum(max);
1710 posZ->SetMinimum(min);
1711 posZ->SetMaximum(max);
0ab29cfa 1712
8ca1a6d9 1713 posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
0ab29cfa 1714
8ca1a6d9 1715 posX->GetYaxis()->SetTitleOffset(1.7);
1716 posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
1717 posY->GetYaxis()->SetTitleOffset(1.7);
1718 posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
1719 posZ->GetYaxis()->SetTitleOffset(1.7);
1720 posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
0ab29cfa 1721
8ca1a6d9 1722 posZ->GetXaxis()->SetRangeUser(0, 1);
0ab29cfa 1723
8ca1a6d9 1724 TString canvasName;
1725 canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
0ab29cfa 1726
8ca1a6d9 1727 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1728 canvas->Divide(3, 1);
0ab29cfa 1729
8ca1a6d9 1730 canvas->cd(1);
1731 InitPad();
1732 posX->DrawCopy();
0ab29cfa 1733
8ca1a6d9 1734 canvas->cd(2);
1735 InitPad();
1736 posY->DrawCopy();
1737
1738 canvas->cd(3);
1739 InitPad();
1740 posZ->DrawCopy();
0ab29cfa 1741
8ca1a6d9 1742 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1743 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1744 }
0ab29cfa 1745}
69b09e3b 1746*/
0ab29cfa 1747
51f6de65 1748void CompareTrack2Particle1D(const char* file1, const char* file2, Float_t upperPtLimit = 9.9)
1749{
1750 loadlibs();
1751
1752 const char* folderName = "dndeta_correction";
1753
1754 c = new TCanvas("CompareTrack2Particle1D", "CompareTrack2Particle1D", 1200, 400);
1755 c->Divide(3, 1);
1756
1757 for (Int_t fileId = 0; fileId < 2; fileId++)
1758 {
1759 const char* file = ((fileId == 0) ? file1 : file2);
1760 Correction1DCreatePlots(file, folderName, upperPtLimit, 1);
1761
1762 TH1* corr[3];
1763 corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
1764 corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y_div_measured_y"));
1765 corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
1766 /*corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x"))->Clone(Form("hist_x_%d", fileId));
1767 corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y"))->Clone(Form("hist_y_%d", fileId));
1768 corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z"))->Clone(Form("hist_z_%d", fileId));*/
1769
1770 for (Int_t i=0; i<3; i++)
1771 {
1772 c->cd(i+1);
1773 InitPad();
1774 corr[i]->GetYaxis()->SetRangeUser(0.8, 2);
1775 corr[i]->SetLineColor(fileId+1);
1776 corr[i]->DrawCopy((fileId == 0) ? "" : "SAME");
1777 }
1778 }
1779
1780 return;
1781
1782 c->SaveAs(Form("%s.gif", canvas->GetName()));
1783 c->SaveAs(Form("%s.eps", canvas->GetName()));
1784}
1785
81be4ee8 1786void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
0ab29cfa 1787{
1788 TFile::Open(fileName);
81be4ee8 1789 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folder, folder);
0448e811 1790 dNdEtaCorrection->LoadHistograms();
0ab29cfa 1791
0448e811 1792 TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1793 TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1afae8ff 1794
69b09e3b 1795 gene->GetZaxis()->SetRangeUser(0.2, 10);
1796 meas->GetZaxis()->SetRangeUser(0.2, 10);
25db2d85 1797 AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
0ab29cfa 1798 gene->GetZaxis()->SetRange(0, 0);
1799 meas->GetZaxis()->SetRange(0, 0);
25db2d85 1800
1801 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1802 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1803 AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
0ab29cfa 1804 gene->GetYaxis()->SetRange(0, 0);
1805 meas->GetYaxis()->SetRange(0, 0);
25db2d85 1806
1807 gene->GetXaxis()->SetRangeUser(-10, 10);
1808 meas->GetXaxis()->SetRangeUser(-10, 10);
1809 AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
0ab29cfa 1810 gene->GetXaxis()->SetRange(0, 0);
1811 meas->GetXaxis()->SetRange(0, 0);
1812}
1813
c17301f3 1814TCanvas* Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
0ab29cfa 1815{
1816 gSystem->Load("libPWG0base");
1817
81be4ee8 1818 Track2Particle2DCreatePlots(fileName, folder);
25db2d85 1819
0448e811 1820 TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1821 TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1822 TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
25db2d85 1823
1afae8ff 1824 Prepare2DPlot(corrYX);
1825 Prepare2DPlot(corrZX);
1826 Prepare2DPlot(corrZY);
1827
0ab29cfa 1828 const char* title = "";
92d2d8ad 1829 corrYX->SetTitle(title);
1830 corrZX->SetTitle(title);
1831 corrZY->SetTitle(title);
1afae8ff 1832
d09fb536 1833 TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1afae8ff 1834 canvas->Divide(3, 1);
1835
1836 canvas->cd(1);
1837 InitPadCOLZ();
1838 corrYX->Draw("COLZ");
1839
1840 canvas->cd(2);
1841 InitPadCOLZ();
1842 corrZX->Draw("COLZ");
1843
1844 canvas->cd(3);
1845 InitPadCOLZ();
1846 corrZY->Draw("COLZ");
92d2d8ad 1847
51f6de65 1848 canvas->SaveAs(Form("corr_track2particle_%d.gif", gMax));
1849 canvas->SaveAs(Form("corr_track2particle_%d.eps", gMax));
c17301f3 1850
1851 return canvas;
0ab29cfa 1852}
1853
1854void CompareTrack2Particle2D()
1855{
1856 gSystem->Load("libPWG0base");
1857
1858 Track2Particle2DCreatePlots("correction_maponly-positive.root");
1859
1c15d51a 1860 TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("pos_yx"));
1861 TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("pos_zx"));
1862 TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("pos_zy"));
0ab29cfa 1863
1864 Track2Particle2DCreatePlots("correction_maponly-negative.root");
1865
1c15d51a 1866 TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("neg_yx"));
1867 TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("neg_zx"));
1868 TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("neg_zy"));
0ab29cfa 1869
1870 posYX->Divide(negYX);
1871 posZX->Divide(negZX);
1872 posZY->Divide(negZY);
1873
1874 Prepare2DPlot(posYX);
1875 Prepare2DPlot(posZX);
1876 Prepare2DPlot(posZY);
1877
1878 Float_t min = 0.8;
1879 Float_t max = 1.2;
1880
1881 posYX->SetMinimum(min);
1882 posYX->SetMaximum(max);
1883 posZX->SetMinimum(min);
1884 posZX->SetMaximum(max);
1885 posZY->SetMinimum(min);
1886 posZY->SetMaximum(max);
1887
1888 TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1889 canvas->Divide(3, 1);
1890
1891 canvas->cd(1);
1892 InitPadCOLZ();
1893 posYX->Draw("COLZ");
1894
1895 canvas->cd(2);
1896 InitPadCOLZ();
1897 posZX->Draw("COLZ");
1898
1899 canvas->cd(3);
1900 InitPadCOLZ();
1901 posZY->Draw("COLZ");
1902
1903 canvas->SaveAs("CompareTrack2Particle2D.gif");
1904 canvas->SaveAs("CompareTrack2Particle2D.eps");
1afae8ff 1905}
1906
1907void Track2Particle3D()
1908{
1909 // get left margin proper
1910
1911 TFile* file = TFile::Open("correction_map.root");
1912
d09fb536 1913 TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1914
1915 corr->SetTitle("Correction Factor");
1916 SetRanges(corr->GetZaxis());
1917
1918 Prepare3DPlot(corr);
1919
1920 TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1921 canvas->SetTheta(29.428);
1922 canvas->SetPhi(16.5726);
1923
1924 corr->Draw();
1925
1926 canvas->SaveAs("Track2Particle3D.gif");
0ab29cfa 1927 canvas->SaveAs("Track2Particle3D.eps");
d09fb536 1928}
1929
1930void Track2Particle3DAll()
1931{
d09fb536 1932 TFile* file = TFile::Open("correction_map.root");
1933
1afae8ff 1934 TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1935 TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1936 TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1937
1938 gene->SetTitle("Generated Particles");
1939 meas->SetTitle("Measured Tracks");
1940 corr->SetTitle("Correction Factor");
1941
1942 Prepare3DPlot(gene);
1943 Prepare3DPlot(meas);
1944 Prepare3DPlot(corr);
1945
d09fb536 1946 TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1afae8ff 1947 canvas->Divide(3, 1);
1948
1949 canvas->cd(1);
1950 InitPad();
1951 gene->Draw();
1952
1953 canvas->cd(2);
1954 meas->Draw();
1955
1956 canvas->cd(3);
1957 corr->Draw();
d09fb536 1958
1959 canvas->SaveAs("Track2Particle3DAll.gif");
0ab29cfa 1960 canvas->SaveAs("Track2Particle3DAll.eps");
1afae8ff 1961}
1962
6b7fa615 1963void MultiplicityMC(Int_t xRangeMax = 50)
4c6b34a8 1964{
1965 TFile* file = TFile::Open("multiplicityMC.root");
1966
1967 if (!file)
1968 {
1969 printf("multiplicityMC.root could not be opened.\n");
1970 return;
1971 }
1972
1973 TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1974 TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1975 TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1976
1977 TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1978 TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1979 //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1980 for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1981 {
1982 TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1983 proj->Fit("gaus", "0");
1984 correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1985 correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1986
1987 continue;
1988
1989 // draw for debugging
1990 new TCanvas;
1991 proj->DrawCopy();
1992 proj->GetFunction("gaus")->DrawCopy("SAME");
1993 }
1994
1995 TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1996
1997 for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1998 {
1999 Float_t mean = correction->GetBinContent(i);
2000 Float_t width = correctionWidth->GetBinContent(i);
2001
2002 Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
2003 Int_t fillEnd = fMultiplicityESDCorrected->FindBin(mean + width * 3);
2004 printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
2005
2006 for (Int_t j=fillBegin; j <= fillEnd; ++j)
2007 {
2008 fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
2009 }
2010 }
2011
2012 TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
2013 fMultiplicityESDCorrectedRebinned->Rebin(10);
2014 fMultiplicityESDCorrectedRebinned->Scale(0.1);
2015
6b7fa615 2016 TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
2017 ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
2018 ratio->Divide(fMultiplicityMC);
2019
4c6b34a8 2020 TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
2021 ratio2->Divide(fMultiplicityMC);
2022
2023 TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
2024 canvas->Divide(3, 2);
2025
6b7fa615 2026 fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
2027 ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
2028 fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
2029 fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
2030 correction->GetXaxis()->SetRangeUser(0, xRangeMax);
2031 fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
2032 fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
2033
2034 canvas->cd(1); //InitPad();
4c6b34a8 2035 fMultiplicityESD->Draw();
2036 fMultiplicityMC->SetLineColor(2);
2037 fMultiplicityMC->Draw("SAME");
2038
6b7fa615 2039 TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
2040 legend->AddEntry(fMultiplicityESD, "ESD");
2041 legend->AddEntry(fMultiplicityMC, "MC");
2042 legend->Draw();
4c6b34a8 2043
6b7fa615 2044 canvas->cd(2);
4c6b34a8 2045 fCorrelation->Draw("COLZ");
2046
6b7fa615 2047 canvas->cd(3);
4c6b34a8 2048 correction->Draw();
2049 //correction->Fit("pol1");
2050 correctionWidth->SetLineColor(2);
2051 correctionWidth->Draw("SAME");
2052
6b7fa615 2053 legend = new TLegend(0.2, 0.7, 0.45, 0.85);
2054 legend->AddEntry(correction, "#bar{x}");
2055 legend->AddEntry(correctionWidth, "#sigma");
2056 legend->Draw();
2057
2058 canvas->cd(4);
2059 ratio->Draw();
2060
2061 ratio2->SetLineColor(2);
2062 ratio2->Draw("SAME");
2063
2064 legend = new TLegend(0.6, 0.7, 0.85, 0.85);
2065 legend->AddEntry(ratio, "uncorrected");
2066 legend->AddEntry(ratio2, "corrected");
2067 legend->Draw();
2068
4c6b34a8 2069 canvas->cd(5);
6b7fa615 2070 fMultiplicityESDCorrected->SetLineColor(kBlue);
4c6b34a8 2071 fMultiplicityESDCorrected->Draw();
2072 fMultiplicityMC->Draw("SAME");
2073 fMultiplicityESD->Draw("SAME");
2074
6b7fa615 2075 legend = new TLegend(0.6, 0.7, 0.85, 0.85);
2076 legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
2077 legend->AddEntry(fMultiplicityMC, "MC");
2078 legend->AddEntry(fMultiplicityESD, "ESD");
2079 legend->Draw();
2080
4c6b34a8 2081 canvas->cd(6);
6b7fa615 2082 fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
4c6b34a8 2083 fMultiplicityESDCorrectedRebinned->Draw();
2084 fMultiplicityMC->Draw("SAME");
6b7fa615 2085
2086 legend = new TLegend(0.6, 0.7, 0.85, 0.85);
2087 legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
2088 legend->AddEntry(fMultiplicityMC, "MC");
2089 legend->Draw();
2090
2091 canvas->SaveAs("MultiplicityMC.gif");
4c6b34a8 2092}
2093
2094void MultiplicityESD()
2095{
2096 TFile* file = TFile::Open("multiplicityESD.root");
2097
2098 if (!file)
2099 {
2100 printf("multiplicityESD.root could not be opened.\n");
2101 return;
2102 }
2103
2104 TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
2105
2106 TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
2107
2108 fMultiplicityESD->Draw();
2109}
2110
69b09e3b 2111void CompareCorrection2Measured(Float_t ptMin = 0.301, const char* dataInput = "analysis_esd_raw.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
770a1f1d 2112{
2113 loadlibs();
2114
2115 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
2116 TFile::Open(correctionMapFile);
2117 dNdEtaCorrection->LoadHistograms();
2118
2119 TFile* file = TFile::Open(dataInput);
2120
2121 if (!file)
2122 {
2123 cout << "Error. File not found" << endl;
2124 return;
2125 }
2126
2127 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
2128 fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
2129
2130 gROOT->cd();
2131
2132 TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
2133 hist1->SetTitle("mc");
3dfa46a4 2134 Printf("mc contains %f entries", hist1->Integral());
81be4ee8 2135 Printf("mc contains %f entries in |vtx-z| < 10, |eta| < 1, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), hist1->GetYaxis()->FindBin(-0.99), hist1->GetYaxis()->FindBin(0.99), hist1->GetZaxis()->FindBin(ptMin), hist1->GetNbinsZ()));
770a1f1d 2136
2137 TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
2138 hist2->SetTitle("esd");
3dfa46a4 2139 Printf("esd contains %f entries", hist2->Integral());
81be4ee8 2140 Printf("esd contains %f entries in |vtx-z| < 10, |eta| < 1, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), hist2->GetYaxis()->FindBin(-0.99), hist2->GetYaxis()->FindBin(0.99), hist2->GetZaxis()->FindBin(ptMin), hist2->GetNbinsZ()));
69b09e3b 2141
2142 AliPWG0Helper::CreateDividedProjections(hist1, hist2);
2143 AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
2144
2145 hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
2146 hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
2147 AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
2148
2149 new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
2150 new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
2151 new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
2152 new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
2153 new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
a7f69e56 2154
2155 TH2* hist3 = (TH2*) dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram()->Clone("mc2");
2156 hist3->SetTitle("mc2");
2157 Printf("mc event contains %f entries", hist3->Integral());
2158 Printf("mc event contains %f entries in |vtx-z| < 10", hist3->Integral(hist3->GetXaxis()->FindBin(-9.9), hist3->GetXaxis()->FindBin(9.9), 1, hist3->GetNbinsY()));
2159
2160 TH2* hist4 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("esd2");
2161 hist4->SetTitle("esd2");
2162 Printf("esd event contains %f entries", hist4->Integral());
2163 Printf("esd event contains %f entries in |vtx-z| < 10", hist4->Integral(hist4->GetXaxis()->FindBin(-9.9), hist4->GetXaxis()->FindBin(9.9), 1, hist4->GetNbinsY()));
2164
2165 ratio = (TH2*) hist3->Clone("ratio");
2166 ratio->Divide(hist4);
2167
2168 new TCanvas; ratio->Draw("COLZ");
69b09e3b 2169}
2170
2171void CompareCorrection2Generated(Float_t ptMin = 0.301, const char* dataInput = "analysis_mc.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
2172{
2173 loadlibs();
2174
2175 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
2176 TFile::Open(correctionMapFile);
2177 dNdEtaCorrection->LoadHistograms();
2178
2179 TFile* file = TFile::Open(dataInput);
2180
2181 if (!file)
2182 {
2183 cout << "Error. File not found" << endl;
2184 return;
2185 }
2186
2187 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
2188 fdNdEtaAnalysis->LoadHistograms("dndetaTrVtx");
2189
2190 gROOT->cd();
2191
2192 TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram()->Clone("mc");
2193 hist1->SetTitle("mc");
2194 Printf("mc contains %f entries", hist1->Integral());
2195 Printf("mc contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), hist1->GetYaxis()->FindBin(-0.99), hist1->GetYaxis()->FindBin(0.99), hist1->GetZaxis()->FindBin(ptMin), hist1->GetNbinsZ()));
2196
2197 TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram()->Clone("esd");
2198 hist2->SetTitle("esd");
2199 Printf("esd contains %f entries", hist2->Integral());
2200 Printf("esd contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), hist2->GetYaxis()->FindBin(-0.99), hist2->GetYaxis()->FindBin(0.99), hist2->GetZaxis()->FindBin(ptMin), hist2->GetNbinsZ()));
770a1f1d 2201
2202 AliPWG0Helper::CreateDividedProjections(hist1, hist2);
0fc41645 2203 AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
2204
2205 hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
2206 hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
2207 AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
770a1f1d 2208
2209 new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
2210 new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
2211 new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
0fc41645 2212 new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
2213 new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
3dfa46a4 2214}
2215
2216void CompareMeasured2Measured(const char* dataInput = "analysis_esd_raw.root", const char* dataInput2 = "analysis_esd_raw.root")
2217{
2218 loadlibs();
2219
2220 TFile* file = TFile::Open(dataInput);
2221
2222 if (!file)
2223 {
2224 cout << "Error. File not found" << endl;
2225 return;
2226 }
2227
2228 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
2229 fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
2230
2231 TFile* file = TFile::Open(dataInput2);
2232
2233 if (!file)
2234 {
2235 cout << "Error. File not found" << endl;
2236 return;
2237 }
2238
2239 dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta2", "dndeta2");
2240 fdNdEtaAnalysis2->LoadHistograms("fdNdEtaAnalysisESD");
2241
2242 gROOT->cd();
2243
2244 TH3* hist1 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd1");
2245 hist1->SetTitle("esd1");
2246 Printf("esd1 contains %f entries", hist1->GetEntries());
2247 Printf("esd1 contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), 1, hist1->GetNbinsY(), hist1->GetZaxis()->FindBin(0.301), hist1->GetNbinsZ()));
2248
2249 TH3* hist2 = (TH3*) fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd2");
2250 hist2->SetTitle("esd2");
2251 Printf("esd2 contains %f entries", hist2->GetEntries());
2252 Printf("esd2 contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), 1, hist2->GetNbinsY(), hist2->GetZaxis()->FindBin(0.301), hist2->GetNbinsZ()));
2253
2254 AliPWG0Helper::CreateDividedProjections(hist1, hist2);
2255
2256 new TCanvas; gROOT->FindObject("esd1_yx_div_esd2_yx")->Draw("COLZ");
2257 new TCanvas; gROOT->FindObject("esd1_zx_div_esd2_zx")->Draw("COLZ");
2258 new TCanvas; gROOT->FindObject("esd1_zy_div_esd2_zy")->Draw("COLZ");
2259
2260 TH2* event1 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event1");
2261 TH2* event2 = (TH2*) fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event2");
2262
2263 Printf("event1 contains %f entries", event1->GetEntries());
2264 Printf("event2 contains %f entries", event2->GetEntries());
2265 Printf("event1 integral is %f", event1->Integral());
2266 Printf("event2 integral is %f", event2->Integral());
2267 Printf("event1 contains %f entries in |vtx-z| < 10", event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), 1, event1->GetNbinsY()));
2268 Printf("event2 contains %f entries in |vtx-z| < 10", event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), 1, event2->GetNbinsY()));
2269
2270 projx1 = event1->ProjectionX();
2271 projx2 = event2->ProjectionX();
2272
2273 new TCanvas; projx1->DrawCopy(); projx2->SetLineColor(2); projx2->DrawCopy("SAME");
2274
2275 projx1->Divide(projx2);
2276 new TCanvas; projx1->Draw();
2277
2278 event1->Divide(event2);
2279 new TCanvas; event1->Draw("COLZ");
770a1f1d 2280
2281}
3dfa46a4 2282
69b09e3b 2283void DrawTrackletOrigin(const char* fileName = "correction_map.root", Bool_t myFile = kTRUE)
0fc41645 2284{
69b09e3b 2285 TFile::Open(fileName);
0fc41645 2286
2287 Int_t maxHists = 8;
2288 TH1* hist[8];
69b09e3b 2289
2290 const Int_t kRebin = 8;
0fc41645 2291
69b09e3b 2292 const char* titles[] = { "PP", "SS", "PP'", "PS'", "PS", "SP'", "SS'", "" };
0fc41645 2293
69b09e3b 2294 if (myFile)
2295 {
2296 for (Int_t i=0; i<maxHists; i++)
2297 {
2298 hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
2299 if (hist[i]->GetDimension() == 2)
2300 hist[i] = ((TH2*) hist[i])->ProjectionX(Form("fDeltaPhi_clone_%d", i));
2301 }
2302 }
2303 else
2304 {
2305 maxHists = 6;
2306 const char* names[] = { "DePhiPPTracklets", "DePhiSecTracklets", "DePhiPpTracklets", "DePhiPSTracklets", "DePhiPSdaugTracklets", "DePhiSPTracklets" };
2307 for (Int_t i=0; i<maxHists; i++)
2308 hist[i] = (TH1*) gFile->Get(names[i]);
2309 }
2310
2311 // clone before rebinning
2312 good = (TH1*) hist[0]->Clone("good");
2313 good->Add(hist[4]);
2314
2315 bad = (TH1*) hist[1]->Clone("bad");
2316 bad->Add(hist[2]);
2317 bad->Add(hist[3]);
2318 bad->Add(hist[5]);
2319 if (myFile)
2320 bad->Add(hist[6]);
2321
2322 c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c");
2323 TH1* ref = 0;
2324 Bool_t nw = kFALSE;
2325 if (!c)
2326 {
2327 c = new TCanvas("c", "c", 600, 600);
2328 nw = kTRUE;
2329 ref = (TH1*) c->GetListOfPrimitives()->At(1);
2330 }
2331 c->cd();
2332 c->SetRightMargin(0.05);
2333 c->SetTopMargin(0.05);
2334 c->SetLogy();
2335 c->SetGridx();
2336 c->SetGridy();
2337
2338 Int_t order[] = { 0, 4, 1, 2, 3, 5, 6, 7 };
2339 //Int_t colors[] = {1,2,4,1,2,4,1,2,4};
2340 Int_t colors[] = {1,2,3,4,6,7,8,102};
2341 Int_t markers[] = {20, 21, 22, 23, 24, 25, 26, 27, 28};
2342
2343 TLegend* legend = new TLegend(0.75, 0.6, 0.93, 0.93);
2344 legend->SetFillColor(0);
2345 legend->SetTextSize(0.04);
0fc41645 2346
2347 Int_t total = 0;
69b09e3b 2348 for (Int_t ii=0; ii<maxHists; ii++)
0fc41645 2349 {
69b09e3b 2350 i = order[ii];
2351
2352 hist[i]->Rebin(kRebin);
0fc41645 2353 hist[i]->SetStats(kFALSE);
2354 hist[i]->SetLineColor(colors[i]);
69b09e3b 2355 hist[i]->SetLineWidth(2);
2356 //hist[i]->SetMarkerStyle(markers[i]);
2357 //hist[i]->SetMarkerColor(colors[i]);
2358 //hist[i]->SetLineStyle(ii+1);
2359 hist[i]->GetXaxis()->SetRangeUser(-0.09, 0.09);
2360 hist[i]->GetYaxis()->SetRangeUser(5, hist[i]->GetMaximum() * 2);
2361 hist[i]->GetYaxis()->SetTitleOffset(1.3);
2362 hist[i]->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2363
2364 if (i == 0 && ref)
2365 hist[i]->Scale(1.0 / hist[i]->GetMaximum() * ref->GetMaximum());
2366
2367 hist[i]->DrawCopy(((i == 0 && nw) ? "" : "SAME"));
0fc41645 2368
2369 total += hist[i]->GetEntries();
2370
2371 if (i != 7)
69b09e3b 2372 legend->AddEntry(hist[i], titles[i], "L");
0fc41645 2373 }
2374
2375 legend->Draw();
69b09e3b 2376 c->SaveAs("spd_tracklets_deltaphi_detailed.eps");
0fc41645 2377
2378 Printf("Total: %d", total);
2379 for (Int_t i=0; i<maxHists; i++)
69b09e3b 2380 Printf("Histogram %d (%s) contains %.2f %% of the entries", i, titles[i], 100.0 * hist[i]->GetEntries() / total);
0fc41645 2381
2382 printf("| Delta phi | Acc. %% | ");
2383 for (Int_t i=0; i<maxHists; i++)
2384 printf("%3s %% | ", titles[i]);
2385 Printf("");
2386
2387 for (Float_t f = 0.01; f < 0.09; f += 0.01)
2388 {
2389 Int_t integralBegin = hist[0]->GetXaxis()->FindBin(-f);
2390 Int_t integralEnd = hist[0]->GetXaxis()->FindBin(f);
2391
2392 Int_t total2 = 0;
2393 for (Int_t i=0; i<maxHists; i++)
2394 total2 += (Int_t) hist[i]->Integral(integralBegin, integralEnd);
2395
2396 printf("| %.2f | %6.2f | ", f, 100.0 * total2 / total);
2397
2398 for (Int_t i=0; i<maxHists; i++)
2399 printf("%6.2f | ", (hist[i]->GetEntries() > 0) ? (100.0 * hist[i]->Integral(integralBegin, integralEnd) / hist[i]->GetEntries()) : -1.0);
2400 Printf("");
2401 }
69b09e3b 2402
2403 eff = new TH1F("eff", ";#Delta#varphi cut (rad.)", 101,-0.0005, 0.1005);
2404 cont = new TH1F("cont", "cont", 101,-0.0005, 0.1005);
2405 signalOverBg = new TH1F("signalOverBg", "signalOverBg", 101,-0.0005, 0.1005);
2406 for (Float_t cut=0.000; cut<0.10; cut += 0.001)
2407 {
2408 Float_t accGood = good->Integral(good->GetXaxis()->FindBin(-cut), good->GetXaxis()->FindBin(cut));
2409 Float_t accBad = bad->Integral(bad->GetXaxis()->FindBin(-cut), bad->GetXaxis()->FindBin(cut));
2410 Float_t sB = accGood / accBad;
2411 eff->Fill(cut, 100.0 * accGood / good->Integral());
2412 cont->Fill(cut, 100.0 * accBad / (accGood + accBad));
2413 signalOverBg->Fill(cut, sB);
2414 }
2415
2416 //new TCanvas; signalOverBg->Draw();
2417
2418 c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c2");
2419 Bool_t nw = kFALSE;
2420 if (!c)
2421 {
2422 c = new TCanvas("c2", "c2", 600, 600);
2423 nw = kTRUE;
2424 }
2425 c->cd();
2426 c->SetRightMargin(0.05);
2427 c->SetTopMargin(0.05);
2428 c->SetGridx();
2429 c->SetGridy();
2430 gPad->SetLogy();
2431 good->Rebin(kRebin);
2432 bad->Rebin(kRebin);
2433 good->GetXaxis()->SetRangeUser(-0.09, 0.09);
2434 good->GetYaxis()->SetTitleOffset(1.3);
2435 good->SetStats(0);
2436 good->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2437 good->DrawCopy((nw) ? "" : "SAME");
2438
2439 bad->SetLineColor(2);
2440 bad->SetLineStyle(2);
2441 bad->SetLineWidth(2);
2442 //bad->SetMarkerColor(2);
2443 //bad->SetMarkerStyle(7);
2444 bad->DrawCopy("SAME");
2445
2446 TLegend* legend = new TLegend(0.2, 0.13, 0.85, 0.25);
2447 legend->SetFillColor(0);
2448 legend->SetTextSize(0.04);
2449 legend->AddEntry(good, "Primaries", "L");
2450 legend->AddEntry(bad, "Secondaries + Background", "L");
2451 legend->Draw();
2452
2453 c->SaveAs("spd_tracklets_deltaphi.eps");
2454
2455 c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c3");
2456 Bool_t nw = kFALSE;
2457 if (!c)
2458 {
2459 c = new TCanvas("c3", "c3", 600, 600);
2460 nw = kTRUE;
2461 }
2462 c->cd();
2463 c->SetRightMargin(0.05);
2464 c->SetTopMargin(0.05);
2465 c->SetGridx();
2466 c->SetGridy();
2467
2468 TLegend* legend = new TLegend(0.5, 0.6, 0.93, 0.75);
2469 legend->SetFillColor(0);
2470 legend->SetTextSize(0.04);
2471 legend->AddEntry(eff, "Efficiency (%)", "L");
2472 legend->AddEntry(cont, "Contamination (%)", "L");
2473
2474 eff->SetStats(0);
2475 eff->GetXaxis()->SetRangeUser(0, 0.08);
2476 eff->GetYaxis()->SetRangeUser(1e-3, 105);
2477 eff->SetLineWidth(2);
2478 eff->DrawCopy((nw) ? "" : "SAME");
2479 cont->SetLineStyle(2);
2480 cont->SetLineWidth(2);
2481 cont->SetLineColor(2);
2482 cont->DrawCopy("SAME");
2483 legend->Draw();
2484
2485 c->SaveAs("spd_tracklets_efficiency.eps");
a7f69e56 2486}
2487
2488void DrawTrackletOrigin_Compare(const char* file1, const char* file2)
2489{
2490 DrawTrackletOrigin(file1);
2491 good1 = (TH1*) gROOT->FindObject("good")->Clone("good1");
2492 bad1 = (TH1*) gROOT->FindObject("bad")->Clone("bad1");
2493
2494 DrawTrackletOrigin(file2);
2495 good2 = (TH1*) gROOT->FindObject("good")->Clone("good2");
2496 bad2 = (TH1*) gROOT->FindObject("bad")->Clone("bad2");
2497
2498 c = new TCanvas("c4", "c4", 600, 600);
2499 c->SetRightMargin(0.05);
2500 c->SetTopMargin(0.05);
2501 c->SetGridx();
2502 c->SetGridy();
2503 gPad->SetLogy();
2504
2505 good1->Draw();
2506 bad1->SetLineColor(1);
2507 bad1->SetMarkerColor(1);
2508 bad1->Draw("SAME");
69b09e3b 2509
a7f69e56 2510 Float_t factor = (good1->Integral() + bad1->Integral()) / (good2->Integral() + bad2->Integral());
69b09e3b 2511
a7f69e56 2512 good2->Scale(factor);
2513 bad2->Scale(factor);
2514
2515 good2->SetLineColor(2);
2516 bad2->SetMarkerColor(2);
2517
2518 good2->Draw("SAME");
2519 bad2->Draw("SAME");
2520
2521 good1->GetYaxis()->SetRangeUser(1, TMath::Max(good1->GetMaximum(), good2->GetMaximum()) * 1.1);
69b09e3b 2522}
a7f69e56 2523
69b09e3b 2524void Tracklets_Asymmetry()
2525{
2526 TFile::Open("correction_map.root");
2527
2528 Int_t maxHists = 7;
2529 TH1* hist[8];
2530
2531 Int_t colors[] = {1,2,3,4,6,7,8,102};
2532 const char* titles[] = { "PP", "SS", "PP'", "PS'", "PS", "SP'", "SS'", "" };
2533
2534 TLegend* legend = new TLegend(0.75, 0.6, 0.93, 0.93);
2535
2536 for (Int_t i=0; i<maxHists; i++)
2537 {
2538 hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
2539 hist[i]->Rebin(10);
2540
2541 for (Int_t j=hist[i]->GetNbinsX()/2; j<=hist[i]->GetNbinsX(); j++)
2542 if (hist[i]->GetBinContent(j) > 0)
2543 hist[i]->SetBinContent(j, (hist[i]->GetBinContent(j) - hist[i]->GetBinContent(hist[i]->GetXaxis()->FindBin(-hist[i]->GetXaxis()->GetBinCenter(j)))) / hist[i]->GetBinContent(j));
2544
2545 hist[i]->SetStats(kFALSE);
2546 hist[i]->SetLineColor(colors[i]);
2547 hist[i]->GetXaxis()->SetRangeUser(0.001, 0.09);
2548 //hist[i]->GetYaxis()->SetRangeUser(5, hist[i]->GetMaximum() * 2);
2549 hist[i]->GetYaxis()->SetTitleOffset(1.3);
2550 hist[i]->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2551 hist[i]->Draw(((i == 0) ? "" : "SAME"));
2552
2553 legend->AddEntry(hist[i], titles[i], "L");
2554 }
2555
2556 legend->Draw();
0fc41645 2557}
567160d6 2558
1c15d51a 2559TH2* GetCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
567160d6 2560{
1c15d51a 2561 // returns the correction factor with pt integrated out
2562
567160d6 2563 loadlibs();
2564
2565 TFile::Open(fileName);
2566
2567 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
2568 if (!dNdEtaCorrection->LoadHistograms())
2569 return;
2570
ea441adf 2571 // hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram();
567160d6 2572
ea441adf 2573 gener = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
2574 measu = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
2575
2576 gener->GetZaxis()->SetRange(gener->GetZaxis()->FindBin(ptmin), gener->GetNbinsZ()+1);
2577 TH2D *gener_xy = gener->Project3D("yx");
2578
2579 measu->GetZaxis()->SetRange(measu->GetZaxis()->FindBin(ptmin), measu->GetNbinsZ()+1);
2580 TH2D *measu_xy = measu->Project3D("yx");
2581
2582 cout << measu->GetZaxis()->FindBin(ptmin) << " " << measu->GetNbinsZ()+1 << endl;
2583
2584 TCanvas *canp = new TCanvas("canp","canp",600,1000);
2585 canp->Divide(1,2,0.0001,0.0001);
2586 canp->cd(1);
2587 gener_xy->Draw("COLZ");
2588 canp->cd(2);
2589 measu_xy->Draw("COLZ");
2590
2591
2592 TCanvas *canpr = new TCanvas("canpr","canpr",700,500);
2593 canpr->cd();
2594 TH2D *proj = new TH2D(*gener_xy);
2595 proj->Divide(measu_xy);
2596
2597// proj = hist->Project3D("yx");
1cbdb1a6 2598 proj->Draw("COLZ");
2599
1c15d51a 2600 return proj;
2601}
2602
2603void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2604{
2605 TH2* proj = GetCorrection(fileName, dirName, ptmin);
2606
ea441adf 2607 const Float_t limit = 5;
1cbdb1a6 2608
2609 TString array = "{";
2610 TString arrayEnd = "}";
2611
2612 for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
2613 {
2614 Int_t begin = -1;
2615 Int_t end = -1;
2616 for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
2617 {
2618 if (begin == -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
2619 begin = x;
2620 if (begin != -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
2621 end = x;
2622 }
2623 Printf("Limits for y = %d are %d to %d", y, begin, end);
2624
2625 if (y > 1)
2626 array += ", ";
2627 array += Form("%d", begin);
2628
2629 if (y > 1)
2630 arrayEnd.Prepend(", ");
2631 arrayEnd.Prepend(Form("%d", (end == -1) ? -1 : proj->GetNbinsX() + 1 - end));
2632 }
2633 array += "}";
2634 arrayEnd.Prepend("{");
2635
2636 Printf("Begin array:");
2637 Printf("%s", array.Data());
567160d6 2638
1cbdb1a6 2639 Printf("End array (mirrored) (should be the same):");
2640 Printf("%s", arrayEnd.Data());
567160d6 2641}
eaa3702a 2642
2643void AverageMultiplicity(const char* fileName = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
2644{
2645 loadlibs();
2646
2647 TFile::Open(fileName);
2648
2649 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
2650 dNdEtaCorrection->LoadHistograms();
2651 TH2* events = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram();
2652 TH3* tracks = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram();
2653
2654 Float_t nEvents = events->Integral(events->GetXaxis()->FindBin(-1), events->GetXaxis()->FindBin(1), 0, events->GetNbinsY()+1);
2655 Float_t nTracks = tracks->Integral(tracks->GetXaxis()->FindBin(-1), tracks->GetXaxis()->FindBin(1), tracks->GetYaxis()->FindBin(-0.39), tracks->GetYaxis()->FindBin(0.59), 0, tracks->GetNbinsZ()+1);
2656
2657 Printf("%f %f --> %f", nEvents, nTracks, nTracks / nEvents);
2658}
2659
1c15d51a 2660void GetAverageCorrectionFactor(Float_t etaRange = 1.5, Float_t vertexRange = 9.9, const char* rawFile = "analysis_esd_raw.root", const char* mcFile = "analysis_mc.root")
2661{
2662 loadlibs();
2663
2664 TFile::Open(rawFile);
2665 dNdEtaAnalysis* raw = new dNdEtaAnalysis("dndeta", "dndeta");
2666 raw->LoadHistograms("fdNdEtaAnalysisESD");
2667 raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->GetXaxis()->SetRangeUser(-vertexRange, vertexRange);
2668 tracks = raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Project3D("y");
2669 events = raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("events", 0, raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->GetNbinsY() + 1);
2670 Float_t nEvents = events->Integral(events->FindBin(-vertexRange), events->FindBin(vertexRange));
2671 tracks->Scale(1.0 / nEvents / tracks->GetBinWidth(1));
2672
2673 TFile::Open(mcFile);
2674 dNdEtaAnalysis* mc = new dNdEtaAnalysis("dndeta", "dndeta");
2675 mc->LoadHistograms("dndetaTrVtx");
2676 mcH = mc->GetdNdEtaPtCutOffCorrectedHistogram(0);
2677
2678 new TCanvas;
2679 mcH->SetLineColor(2);
2680 mcH->DrawCopy();
2681 tracks->DrawCopy("SAME");
2682
2683 new TCanvas;
2684 mcH->GetYaxis()->SetRangeUser(0, 5);
2685 mcH->Divide(tracks);
2686 mcH->DrawCopy();
2687 mcH->Fit("pol0", "", "", -etaRange, etaRange);
2688}
2689
a7f69e56 2690void TrackCuts_Comparison_MC(char* histName, Int_t plotWhich = 0, const char* fileName = "correction_map.root", Bool_t mirror = kFALSE)
1c15d51a 2691{
2692 // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
2693 // --> manually disable it in the run.C
5a6310fe 2694 //
2695 // plotWhich: 0 = only before
2696 // 1 = both
2697 // 2 = only after
69b09e3b 2698 //
2699 // mirror: kTRUE --> project negative values on the positive side
2700
1c15d51a 2701
2702 file = TFile::Open(fileName);
2703
2704 Int_t count = 0;
2705 Int_t colors[] = { 1, 2, 3, 4, 5, 6 };
2706
5a6310fe 2707 TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
1c15d51a 2708 TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
5a6310fe 2709 TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
1c15d51a 2710
69b09e3b 2711 TCanvas* c1 = new TCanvas(histName, histName, 800, 1200);
5a6310fe 2712 c1->Divide(1, 2);
1c15d51a 2713 //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
5a6310fe 2714 //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
1c15d51a 2715
2716 const char* folders2[] = { "before_cuts", "after_cuts" };
5a6310fe 2717 Bool_t first = kTRUE;
2718 for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
1c15d51a 2719 {
2720 const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
5a6310fe 2721 const char* names[] = { "all", "primaries", "secondaries" };
1c15d51a 2722 TH1* base = 0;
2723 TH1* prim = 0;
2724 TH1* sec = 0;
2725 for (Int_t i = 0; i < 3; i++)
2726 {
2727 TString folder;
2728 folder.Form("%s/%s/%s", folders1[i], folders2[j], histName);
2729 TH1* hist = (TH1*) file->Get(folder);
69b09e3b 2730
2731 if (mirror)
2732 {
2733 for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
2734 {
2735 Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
2736 if (bin != newBin)
2737 {
2738 hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
2739 hist->SetBinContent(bin, 0);
2740 }
2741 }
2742 }
2743
5a6310fe 2744 legend->AddEntry(hist, Form("%s %s", names[i], folders2[j]));
1c15d51a 2745
5a6310fe 2746 c1->cd(1);
1c15d51a 2747 hist->SetLineColor(colors[count]);
2748 hist->DrawCopy((count == 0) ? "" : "SAME");
2749
2750 switch (i)
2751 {
2752 case 0: base = hist; break;
2753 case 1: prim = hist; break;
2754 case 2: sec = hist; break;
2755 }
2756
2757 count++;
2758 }
2759
2760 TH1* eff = (TH1*) prim->Clone("eff"); eff->Reset();
2761 TH1* purity = (TH1*) prim->Clone("purity"); purity->Reset();
2762
2763 for (Int_t bin = 1; bin <= prim->GetNbinsX(); bin++)
2764 {
2765 eff->SetBinContent(bin, prim->Integral(1, bin) / prim->Integral(1, prim->GetNbinsX() + 1));
5a6310fe 2766 if (prim->Integral(1, bin) + sec->Integral(1, bin) > 0)
2767 purity->SetBinContent(bin, sec->Integral(1, bin) / (prim->Integral(1, bin) + sec->Integral(1, bin)));
1c15d51a 2768 }
2769
2770 eff->GetYaxis()->SetRangeUser(0, 1);
2771 eff->SetLineColor(colors[0+j*2]);
2772 eff->SetStats(kFALSE);
2773 purity->SetLineColor(colors[1+j*2]);
2774
2775 legend3->AddEntry(eff, Form("%s: efficiency", folders2[j]));
2776 legend3->AddEntry(purity, Form("%s: contamination", folders2[j]));
2777
5a6310fe 2778 c1->cd(2);
2779 eff->DrawCopy((first) ? "" : "SAME");
2780 first = kFALSE;
1c15d51a 2781 purity->DrawCopy("SAME");
2782 }
2783
5a6310fe 2784 c1->cd(1)->SetLogy();
2785 c1->cd(1)->SetGridx();
2786 c1->cd(1)->SetGridy();
1c15d51a 2787 legend->Draw();
2788
2789 //c2->cd();
2790 // c2->SetGridx();
2791 // c2->SetGridy();
2792 //legend2->Draw();
2793
5a6310fe 2794 c1->cd(2)->SetGridx();
2795 c1->cd(2)->SetGridy();
1c15d51a 2796 legend3->Draw();
5a6310fe 2797
69b09e3b 2798 //c1->SaveAs(Form("%s.png", histName));
1c15d51a 2799}
2800
a7f69e56 2801void TrackCuts_Comparison_Data(char* histName, Int_t plotWhich, const char* fileName1, const char* fileName2, Bool_t mirror = kFALSE, const char* label1 = "file1", const char* label2 = "file2")
2802{
2803 // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
2804 // --> manually disable it in the run.C
2805 //
2806 // plotWhich: 0 = only before
2807 // 1 = both
2808 // 2 = only after
2809 //
2810 // mirror: kTRUE --> project negative values on the positive side
2811
2812
2813 Int_t count = 0;
2814 Int_t colors[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
2815
2816 TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
2817 legend->SetFillColor(0);
2818 TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
2819 TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
2820
2821 TCanvas* c1 = new TCanvas(histName, histName, 600, 600);
2822 //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
2823 //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
2824
2825 const char* folders2[] = { "before_cuts", "after_cuts" };
2826 Bool_t first = kTRUE;
2827 for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
2828 {
2829 const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
2830 const char* names[] = { "all", "primaries", "secondaries" };
2831
2832 Float_t normalize[3];
2833
2834 for (Int_t i = 0; i < 2; i++)
2835 {
2836 file = TFile::Open((i == 0) ? fileName1 : fileName2);
2837
2838 for (Int_t k = 1; k < 3; k++)
2839 {
2840 TString folder;
2841 folder.Form("%s/%s/%s", folders1[k], folders2[j], histName);
2842 Printf("%s", folder.Data());
2843 TH1* hist = (TH1*) file->Get(folder);
2844
2845 if (mirror)
2846 {
2847 for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
2848 {
2849 Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
2850 if (bin != newBin)
2851 {
2852 hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
2853 hist->SetBinContent(bin, 0);
2854 }
2855 }
2856 }
2857
2858 if (i == 0)
2859 {
2860 normalize[k] = hist->Integral();
2861 }
2862 else
2863 hist->Scale(normalize[k] / hist->Integral());
2864
2865 legend->AddEntry(hist, Form("%s %s %s", (i == 0) ? label1 : label2, (k == 1) ? "primaries" : "secondaries", folders2[j]));
2866
2867 c1->cd();
2868 hist->SetStats(0);
2869 hist->SetLineColor(colors[count]);
2870 hist->DrawCopy((count == 0) ? "" : "SAME");
2871
2872 count++;
2873 }
2874 }
2875
2876 }
2877
2878 //c1->SetLogy();
2879 c1->SetGridx();
2880 c1->SetGridy();
2881 legend->Draw();
2882}
2883
1c15d51a 2884void TrackCuts_DCA()
2885{
2886 file = TFile::Open("correction_map.root");
2887 hist = (TH2*) file->Get("esd_track_cuts/before_cuts/dXYvsDZ");
2888
2889 TCanvas* c1 = new TCanvas("c1", "c1", 600, 600);
2890 c1->SetLogz();
2891 c1->SetRightMargin(0.12);
2892 c1->SetBottomMargin(0.12);
2893
2894 hist->SetStats(kFALSE);
2895 hist->Draw("COLZ");
2896
2897 ellipse = new TEllipse(0, 0, 4);
2898 ellipse->SetLineWidth(2);
2899 ellipse->SetLineStyle(2);
2900 ellipse->SetFillStyle(0);
2901 ellipse->Draw();
2902
2903 c1->SaveAs("trackcuts_dca_2d.eps");
2904}
2905
2906void FindNSigma(TH2* hist, Int_t nSigma = 1)
2907{
2908 TH1* proj = hist->ProjectionY();
2909 proj->Reset();
2910
2911 for (Int_t bin=1; bin<=proj->GetNbinsX(); bin++)
2912 {
2913 if (hist->Integral(1, hist->GetNbinsX(), bin, bin) == 0)
2914 continue;
2915
2916 Int_t limit = -1;
2917 for (limit = 1; limit<=hist->GetNbinsX(); limit++)
2918 }
2919}
2920
2921void ShowOnlyAccepted(TH2* input, const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2922{
2923 TH2* proj = GetCorrection(fileName, dirName, ptmin);
2924
2925 for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
2926 for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
2927 if (proj->GetBinContent(x, y) > 5 || proj->GetBinContent(x, y) == 0)
2928 {
2929 proj->SetBinContent(x, y, 0);
2930 }
2931 else
2932 proj->SetBinContent(x, y, 1);
2933
2934
2935 input->Multiply(proj);
2936}
51f6de65 2937
2938void MakeGaussianProfile(const char* histName = "fVertexCorrelation", Bool_t subtractMean = kFALSE)
2939{
2940 TFile::Open("correction_map.root");
2941
2942 TH2* hist2d = (TH2*) gFile->Get(histName);
2943 hist2d->Sumw2();
2944
2945 TH1* result = hist2d->ProjectionX("result");
2946 result->GetYaxis()->SetTitle(hist2d->GetYaxis()->GetTitle());
2947 result->Reset();
2948
2949 for (Int_t x=1; x<hist2d->GetNbinsX(); ++x)
2950 {
2951 hist = hist2d->ProjectionY(Form("temp_%d", x), x, x);
2952 if (hist->GetEntries() == 0)
2953 continue;
2954 if (hist->Fit("gaus") == 0)
2955 {
2956 func = hist->GetFunction("gaus");
2957 mean = func->GetParameter(1);
2958 error = func->GetParError(1);
2959
2960 if (subtractMean)
2961 mean = hist2d->GetXaxis()->GetBinCenter(x) - mean;
2962
2963 result->SetBinContent(x, mean);
2964 result->SetBinError(x, error);
2965
2966 if (x % 10 == 0)
2967 {
2968 new TCanvas;
2969 ((TH1*) hist->Clone())->DrawCopy();
2970 }
2971 }
2972 //break;
2973 }
2974
2975 new TCanvas;
2976 result->GetYaxis()->SetRangeUser(-0.2, 0.2);
2977 result->Draw();
2978}
69b09e3b 2979
2980TH2* GetAcceptance(void* corr2d_void)
2981{
2982 corr2d = (AliCorrectionMatrix2D*) corr2d_void;
2983 corr_xy = (TH2*) corr2d->GetCorrectionHistogram()->Clone("acceptance");
2984
2985 // fold in acceptance
2986 for (Int_t x=1; x<=corr_xy->GetNbinsX(); ++x)
2987 for (Int_t y=1; y<=corr_xy->GetNbinsY(); ++y)
2988 {
2989 if (corr_xy->GetBinContent(x, y) > 1.5)
2990 corr_xy->SetBinContent(x, y, 0);
2991
2992 if (corr_xy->GetBinContent(x, y) > 0)
2993 corr_xy->SetBinContent(x, y, 1);
2994
2995 corr_xy->SetBinError(x, y, 0);
2996 }
2997
2998 return corr_xy;
2999}
3000
3001void ZeroOutsideAcceptance(TH2* acc, TH3* hist)
3002{
3003 for (Int_t x=0; x<=acc->GetNbinsX()+1; ++x)
3004 for (Int_t y=0; y<=acc->GetNbinsY()+1; ++y)
3005 {
1d107532 3006 if (acc->GetBinContent(x, y) > 2 || acc->GetBinContent(x, y) == 0)
69b09e3b 3007 {
3008 for (Int_t z=0; z<=hist->GetNbinsZ()+1; ++z)
3009 {
3010 hist->SetBinContent(x, y, z, 0);
3011 hist->SetBinError(x, y, z, 0);
3012 }
3013 }
3014 }
3015}
3016
3017void DrawPhi()
3018{
3019 loadlibs();
3020
3021 TFile::Open("correction_map.root");
3022 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
3023 if (!dNdEtaCorrection->LoadHistograms())
3024 return 0;
3025
3026 TFile::Open("analysis_esd.root");
3027 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
3028 fdNdEtaAnalysis->LoadHistograms();
3029
3030 // acc. map!
3031 //acc = GetAcceptance(dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000));
3032 acc = dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000)->GetCorrectionHistogram();
3033 //new TCanvas; acc->Draw("COLZ");
3034
3035 histG = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
3036 ZeroOutsideAcceptance(acc, histG);
3037 //new TCanvas; histG->Project3D("yx")->Draw("COLZ");
3038 //histG->GetYaxis()->SetRangeUser(-0.9, 0.9);
3039
3040 histM = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
3041 ZeroOutsideAcceptance(acc, histM);
3042 //histM->GetYaxis()->SetRangeUser(-0.9, 0.9);
3043
3044 TFile::Open("analysis_mc.root");
3045 dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndetaTrVtxMC", "dndetaTrVtxMC");
3046 fdNdEtaAnalysis2->LoadHistograms("dndetaTrVtx");
3047
3048 histMC = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
3049 ZeroOutsideAcceptance(acc, histMC);
3050 //new TCanvas; histMC->Project3D("yx2")->Draw("COLZ");
3051
3052 //histG->GetZaxis()->SetRangeUser(1,2); histMC->GetZaxis()->SetRangeUser(1,2);
3053 new TCanvas; a = histG->Project3D("yx3"); a->Add(histMC->Project3D("yx4"), -1); a->Draw("COLZ");
3054
3055 //histMC->GetYaxis()->SetRangeUser(-0.9, 0.9);
3056
3057 c = new TCanvas;
3058
3059 histG->GetXaxis()->SetRangeUser(-9.9, 9.9);
1d107532 3060 histG->Project3D("z")->DrawCopy();
69b09e3b 3061
3062 histM->GetXaxis()->SetRangeUser(-9.9, 9.9);
3063 proj = histM->Project3D("z2");
3064 proj->SetLineColor(2);
1d107532 3065 proj->DrawCopy("SAME");
69b09e3b 3066
3067 histMC->GetXaxis()->SetRangeUser(-9.9, 9.9);
3068 projMC = histMC->Project3D("z3");
1d107532 3069 projMC->SetLineColor(4);
3070 projMC->DrawCopy("SAME");
69b09e3b 3071}
3072
1d107532 3073void PrintEventStats(Int_t corrID = 3, const char* fileName = "correction_map.root", const char* dir = "dndeta_correction")
69b09e3b 3074{
3075 loadlibs();
3076
3077 /*
3078 TFile::Open("analysis_mc.root");
3079 fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD");
3080 fdNdEtaAnalysis->LoadHistograms();
3081 trackHist = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
3082 eventHist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetGeneratedHistogram();
3083 */
3084
1d107532 3085 TFile::Open(fileName);
3086 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dir, dir);
69b09e3b 3087 if (!dNdEtaCorrection->LoadHistograms())
3088 return;
3089 trackHist = dNdEtaCorrection->GetCorrection(corrID)->GetTrackCorrection()->GetGeneratedHistogram();
3090 eventHist = dNdEtaCorrection->GetCorrection(corrID)->GetEventCorrection()->GetGeneratedHistogram();
3091
3092 trackHist->GetXaxis()->SetRange(0, trackHist->GetNbinsX()+1);
3093 trackHist->GetZaxis()->SetRange(0, trackHist->GetNbinsZ()+1);
3094 eta = trackHist->Project3D("y");
3095
3096 events = eventHist->Integral(0, eventHist->GetNbinsX()+1, 0, eventHist->GetNbinsY()+1);
3097
3098 eta->Scale(1.0 / events);
3099
3100 Float_t avgN = eta->Integral(0, eta->GetNbinsX()+1);
3101 Printf("<N> = %f", avgN);
3102
3103 eta->Scale(1.0 / eta->GetXaxis()->GetBinWidth(1));
3104
3105 Printf("dndeta | eta = 0 is %f", (eta->GetBinContent(eta->FindBin(0.01)) + eta->GetBinContent(eta->FindBin(-0.01))) / 2);
a7f69e56 3106 Printf("dndeta in |eta| < 0.5 is %f", eta->Integral(eta->FindBin(-0.49), eta->FindBin(0.49)) / (eta->FindBin(0.49) - eta->FindBin(-0.49) + 1));
69b09e3b 3107 Printf("dndeta in |eta| < 1 is %f", eta->Integral(eta->FindBin(-0.99), eta->FindBin(0.99)) / (eta->FindBin(0.99) - eta->FindBin(-0.99) + 1));
3108 Printf("dndeta in |eta| < 1.5 is %f", eta->Integral(eta->FindBin(-1.49), eta->FindBin(1.49)) / (eta->FindBin(1.49) - eta->FindBin(-1.49) + 1));
3109
3110 stats = (TH2*) gFile->Get("fEventStats");
3111 proj = stats->ProjectionX();
3112 gROOT->ProcessLine(".L PrintHist.C");
3113 PrintHist2D(stats);
3114 PrintHist(proj);
3115
81be4ee8 3116 Float_t ua5_SD = 0.153;
3117 Float_t ua5_DD = 0.080;
3118 Float_t ua5_ND = 0.767;
3119
3120 Printf("+++ FRACTIONS +++");
3121
3122 Printf("ND: %f", proj->GetBinContent(3) / proj->GetBinContent(1));
3123 Printf("SD: %f", proj->GetBinContent(4) / proj->GetBinContent(1));
3124 Printf("DD: %f", proj->GetBinContent(5) / proj->GetBinContent(1));
3125
69b09e3b 3126 Printf("+++ TRIGGER EFFICIENCIES +++");
3127
a7f69e56 3128 Printf("INEL = %.1f", 100. * (proj->GetBinContent(1) - stats->GetBinContent(1, 1) - stats->GetBinContent(1, 3)) / proj->GetBinContent(1));
3129 Printf("NSD = %.1f", 100. * (proj->GetBinContent(2) - stats->GetBinContent(2, 1) - stats->GetBinContent(2, 3)) / proj->GetBinContent(2));
81be4ee8 3130
3131
3132 Float_t trigND = 100. * (proj->GetBinContent(3) - stats->GetBinContent(3, 1) - stats->GetBinContent(3, 3)) / proj->GetBinContent(3);
3133 Float_t trigSD = 100. * (proj->GetBinContent(4) - stats->GetBinContent(4, 1) - stats->GetBinContent(4, 3)) / proj->GetBinContent(4);
3134 Float_t trigDD = 100. * (proj->GetBinContent(5) - stats->GetBinContent(5, 1) - stats->GetBinContent(5, 3)) / proj->GetBinContent(5);
3135
3136 Printf("ND = %.1f", trigND);
3137 Printf("SD = %.1f", trigSD);
3138 Printf("DD = %.1f", trigDD);
3139
3140 Float_t trigINELUA5 = ua5_SD * trigSD + ua5_DD * trigDD + ua5_ND * trigND;
3141 Float_t trigNSDUA5 = (ua5_DD * trigDD + ua5_ND * trigND) / (ua5_DD + ua5_ND);
3142 Printf("INEL (UA5) = %.1f", trigINELUA5);
3143 Printf("NSD (UA5) = %.1f", trigNSDUA5);
a7f69e56 3144
1d107532 3145 Printf("+++ VERTEX EFFICIENCIES +++");
3146
3147 Printf("INEL = %.1f", 100. * (stats->GetBinContent(1, 3) + stats->GetBinContent(1, 4)) / proj->GetBinContent(1));
3148 Printf("NSD = %.1f", 100. * (stats->GetBinContent(2, 3) + stats->GetBinContent(2, 4)) / proj->GetBinContent(2));
3149
3150 Float_t vtxND = 100. * (stats->GetBinContent(3, 3) + stats->GetBinContent(3, 4)) / proj->GetBinContent(3);
3151 Float_t vtxSD = 100. * (stats->GetBinContent(4, 3) + stats->GetBinContent(4, 4)) / proj->GetBinContent(4);
3152 Float_t vtxDD = 100. * (stats->GetBinContent(5, 3) + stats->GetBinContent(5, 4)) / proj->GetBinContent(5);
3153 Printf("ND = %.1f", vtxND);
3154 Printf("SD = %.1f", vtxSD);
3155 Printf("DD = %.1f", vtxDD);
3156
1d107532 3157 Float_t vtxINELUA5 = ua5_SD * vtxSD + ua5_DD * vtxDD + ua5_ND * vtxND;
3158 Float_t vtxNSDUA5 = (ua5_DD * vtxDD + ua5_ND * vtxND) / (ua5_DD + ua5_ND);
3159 Printf("INEL (UA5) = %.1f", vtxINELUA5);
3160 Printf("NSD (UA5) = %.1f", vtxNSDUA5);
3161
a7f69e56 3162 Printf("+++ TRIGGER + VERTEX EFFICIENCIES +++");
3163
3164 Printf("INEL = %.1f", 100. * stats->GetBinContent(1, 4) / proj->GetBinContent(1));
3165 Printf("NSD = %.1f", 100. * stats->GetBinContent(2, 4) / proj->GetBinContent(2));
3166 Printf("ND = %.1f", 100. * stats->GetBinContent(3, 4) / proj->GetBinContent(3));
3167 Printf("SD = %.1f", 100. * stats->GetBinContent(4, 4) / proj->GetBinContent(4));
3168 Printf("DD = %.1f", 100. * stats->GetBinContent(5, 4) / proj->GetBinContent(5));
69b09e3b 3169
1d107532 3170
3171
69b09e3b 3172 for (Int_t i=7; i<=proj->GetNbinsX(); i++)
3173 if (proj->GetBinContent(i) > 0)
3174 Printf("bin %d (process type %d) = %.2f", i, (Int_t) proj->GetXaxis()->GetBinCenter(i), 100.0 * (proj->GetBinContent(i) - stats->GetBinContent(i, 1)) / proj->GetBinContent(i));
3175
3176 //eta->Draw();
3177}
3178
3179void TestAsymmetry()
3180{
3181 loadlibs();
3182
3183 TFile* file2 = TFile::Open("analysis_mc.root");
3184
3185 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
3186 fdNdEtaAnalysis->LoadHistograms();
3187 fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
3188
3189 fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy();
3190
3191 hist = (TH1*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
3192 hist2 = (TH1*) hist->Clone("hist2");
3193 for (Int_t x=1; x<=hist->GetNbinsX(); x++)
3194 for (Int_t y=1; y<=hist->GetNbinsY(); y++)
3195 for (Int_t z=1; z<=hist->GetNbinsZ(); z++)
3196 {
3197 Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
3198 hist->SetBinContent(x, y, z, hist2->GetBinContent(hist->GetNbinsX() / 2, TMath::Min(y, hist->GetNbinsY() + 1 - y), z));
3199 }
3200
3201 hist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
3202 for (Int_t x=1; x<=hist->GetNbinsX(); x++)
3203 for (Int_t y=1; y<=hist->GetNbinsY(); y++)
3204 {
3205 //Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
3206 hist->SetBinContent(x, y, hist->GetBinContent(hist->GetNbinsX() / 2, y));
3207 }
3208
3209 fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
3210 fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerColor(2);
3211 fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetLineColor(2);
3212 fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerStyle(5);
3213 fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy("SAMEP");
3214}
3215
3216void DeltaPhiFromPt(Float_t smearing = 0.005)
3217{
3218 loadlibs();
3219
3220 TFile::Open("analysis_mc.root");
3221 hist = (TH1*) gFile->Get("dndeta_check_pt");
3222
3223 dPhiHist = new TH1F("dPhiHist", ";#Delta phi", 400, -0.1, 0.1);
3224 dPhiHist2 = new TH1F("dPhiHist2", ";#Delta phi", 400, -0.1, 0.1);
3225
3226 for (Int_t i=1; i<=hist->GetNbinsX(); i++)
3227 {
3228 Float_t pt = hist->GetBinCenter(i);
3229 Float_t deltaPhi = (0.076 - 0.039) / 2 / (pt / 0.15);
3230
3231 if (smearing > 0)
3232 {
3233 gaus = new TF1("mygaus", "gaus(0)", -0.1, 0.1);
3234 gaus->SetParameters(1, -deltaPhi, smearing);
3235
3236 dPhiHist->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
3237
3238 dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
3239 gaus->SetParameters(1, deltaPhi, smearing);
3240 dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
3241 }
3242 else
3243{
3244dPhiHist->Fill(deltaPhi, hist->GetBinContent(i) / 2);
3245dPhiHist2->Fill(deltaPhi, hist->GetBinContent(i) / 2);
3246dPhiHist2->Fill(-deltaPhi, hist->GetBinContent(i) / 2);
3247}
3248 }
3249
3250 new TCanvas;
3251 dPhiHist->Draw();
3252 dPhiHist2->SetLineColor(2);
3253 dPhiHist2->Draw("SAME");
3254 gPad->SetLogy();
3255
3256 TFile::Open("trackletsDePhi.root");
3257 //TFile::Open("tmp/correction_maponly-positive.root");
3258 //TFile::Open("tmp/correction_map.root");
3259 //tracklets = (TH1*) gFile->Get(Form("fDeltaPhi_%d", 0));
3260 tracklets = (TH1*) gFile->Get("DePhiPPTracklets");
3261 tracklets->Scale(1.0 / tracklets->GetMaximum() * dPhiHist->GetMaximum());
3262 tracklets->SetLineColor(4);
3263 tracklets->Draw("SAME");
3264}
a7f69e56 3265
3266void VertexDistributions()
3267{
3268 loadlibs();
3269
3270 TFile::Open("correction_map.root");
3271 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
3272 if (!dNdEtaCorrection->LoadHistograms())
3273 return;
3274
3275 all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all");
3276 trigger = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("trigger");
3277 vtx = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("vtx");
3278
3279 nottriggered = (TH1*) all->Clone("nottriggered");
3280 nottriggered->Add(trigger, -1);
3281
3282 novertex = (TH1*) trigger->Clone("novertex");
3283 novertex->Add(vtx, -1);
3284
3285 temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
3286 highmult = temphist->ProjectionX("highmult", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
3287 //all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
3288
3289 for (Int_t i=1; i<=trigger->GetNbinsX(); i++)
3290 {
3291 all->SetBinContent(i, all->GetBinContent(i) / all->GetBinWidth(i));
3292 trigger->SetBinContent(i, trigger->GetBinContent(i) / trigger->GetBinWidth(i));
3293 vtx->SetBinContent(i, vtx->GetBinContent(i) / vtx->GetBinWidth(i));
3294 nottriggered->SetBinContent(i, nottriggered->GetBinContent(i) / nottriggered->GetBinWidth(i));
3295 novertex->SetBinContent(i, novertex->GetBinContent(i) / novertex->GetBinWidth(i));
3296 highmult->SetBinContent(i, highmult->GetBinContent(i) / highmult->GetBinWidth(i));
3297 }
3298
3299 new TCanvas;
3300 vtx->SetTitle("");
3301 vtx->SetStats(0);
3302 vtx->DrawCopy("HIST");
3303
3304 all->Scale(1.0 / all->Integral());
3305 nottriggered->Scale(1.0 / nottriggered->Integral());
3306 novertex->Scale(1.0 / novertex->Integral());
3307 highmult->Scale(1.0 / highmult->Integral());
3308
3309 new TCanvas;
3310 all->Draw("HIST");
3311 novertex->SetLineColor(2);
3312 novertex->Draw("HISTSAME");
3313 highmult->SetLineColor(3);
3314 highmult->Draw("HISTSAME");
3315
3316 legend = new TLegend(0.5, 0.5, 0.8, 0.8);
3317 legend->SetFillColor(0);
3318 legend->AddEntry(all, "all");
3319 legend->AddEntry(novertex, "no vertex");
3320 legend->AddEntry(highmult, "mult > 10");
3321 legend->Draw();
3322
3323 new TCanvas;
3324 trigger->Scale(1.0 / trigger->Integral());
3325 vtx->Scale(1.0 / vtx->Integral());
3326
3327 trigger->Divide(vtx);
3328
3329 trigger->Draw();
3330 //vtx->SetLineColor(2);
3331 //vtx->Draw("SAME");
3332
3333 //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
3334 temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram();
3335 //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram();
3336
70fdd197 3337 temphist = (TH2*) gFile->Get("fTemp1");
3338
a7f69e56 3339 new TCanvas;
3340 legend = new TLegend(0.7, 0.7, 0.99, 0.99);
3341 legend->SetFillColor(0);
3342
3343 Bool_t first = kTRUE;
3344 for (Int_t i=0; i<20; i+=5)
3345 {
3346 highmult = temphist->ProjectionX("highmult", i+1, i+1+4);
70fdd197 3347 highmult->Rebin(10);
a7f69e56 3348 Printf("%f", highmult->Integral());
3349 if (highmult->Integral() <= 0)
3350 continue;
3351
3352 for (Int_t j=1; j<=trigger->GetNbinsX(); j++)
3353 highmult->SetBinContent(j, highmult->GetBinContent(j) / highmult->GetBinWidth(j));
3354
3355 highmult->Scale(1.0 / highmult->Integral());
3356 highmult->SetLineColor((i/5)+1);
3357 highmult->GetYaxis()->SetRangeUser(0, 0.15);
3358 if (first)
3359 {
3360 highmult->DrawCopy();
3361 first = kFALSE;
3362 }
3363 else
3364 highmult->DrawCopy("SAME");
3365 legend->AddEntry(highmult->Clone(), Form("%d <= N <= %d", i, i+4));
3366 }
3367 legend->Draw();
3368
3369}
3370
3371void PlotPt1DCorrection()
3372{
3373 const char* files[] = { "field.root", "field_onlyprim.root", "nofield.root", "nofield_onlyprim.root" };
3374 const char* names[] = { "B: all", "B: primaries", "No B: all", "No B: primaries" };
3375 Int_t colors[] = { 1, 2, 3, 4 };
3376
3377 loadlibs();
3378
3379 dummy = new TH2F("dummy", ";p_{T};correction", 100, 0, 1.4, 100, 0.5, 3);
3380 dummy->SetStats(0);
3381 //dummy->GetYaxis()->SetTitleOffset(1.3);
3382 dummy->Draw();
3383
3384 legend = new TLegend(0.48, 0.57, 0.88, 0.88);
3385 legend->SetFillColor(0);
3386
3387 for (Int_t i=0; i<4; i++)
3388 {
3389 TFile::Open(files[i]);
3390 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
3391 if (!dNdEtaCorrection->LoadHistograms())
3392 return;
3393
3394 hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get1DCorrectionHistogram("z", -9.9, 9.9, -0.79, 0.79);