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