fixing warnings
[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{
25 gSystem->Load("libANALYSIS");
26 gSystem->Load("libPWG0base");
27}
28
0ab29cfa 29void SetRanges(TAxis* axis)
92d2d8ad 30{
0ab29cfa 31 if (strcmp(axis->GetTitle(), "#eta") == 0)
32 axis->SetRangeUser(-1.7999, 1.7999);
33 if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0)
34 axis->SetRangeUser(0, 9.9999);
35 if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
36 axis->SetRangeUser(-15, 14.9999);
37 if (strcmp(axis->GetTitle(), "Ntracks") == 0)
38 axis->SetRangeUser(0, 99.9999);
25db2d85 39}
40
0ab29cfa 41void SetRanges(TH1* hist)
25db2d85 42{
0ab29cfa 43 SetRanges(hist->GetXaxis());
44 SetRanges(hist->GetYaxis());
45 SetRanges(hist->GetZaxis());
46}
25db2d85 47
0ab29cfa 48
49void Prepare3DPlot(TH3* hist)
50{
51 hist->GetXaxis()->SetTitleOffset(1.5);
52 hist->GetYaxis()->SetTitleOffset(1.5);
53 hist->GetZaxis()->SetTitleOffset(1.5);
54
55 hist->SetStats(kFALSE);
56}
57
58void Prepare2DPlot(TH2* hist)
59{
60 hist->SetStats(kFALSE);
61 hist->GetYaxis()->SetTitleOffset(1.4);
62
63 hist->SetMinimum(0);
64 hist->SetMaximum(gMax);
65
66 SetRanges(hist);
67}
68
69void Prepare1DPlot(TH1* hist)
70{
71 hist->SetLineWidth(2);
72 hist->SetStats(kFALSE);
73
72e597d7 74 hist->GetXaxis()->SetLabelOffset(0.02);
75 hist->GetXaxis()->SetTitleOffset(1.3);
76 hist->GetYaxis()->SetTitleOffset(1.3);
77
0ab29cfa 78 SetRanges(hist);
79}
80
81void InitPad()
82{
83 gPad->Range(0, 0, 1, 1);
84 gPad->SetLeftMargin(0.15);
85 //gPad->SetRightMargin(0.05);
86 //gPad->SetTopMargin(0.13);
72e597d7 87 gPad->SetBottomMargin(0.12);
0ab29cfa 88
72e597d7 89 gPad->SetGridx();
90 gPad->SetGridy();
0ab29cfa 91}
92
93void InitPadCOLZ()
94{
95 gPad->Range(0, 0, 1, 1);
96 gPad->SetRightMargin(0.15);
97 gPad->SetLeftMargin(0.12);
98
99 gPad->SetGridx();
100 gPad->SetGridy();
d09fb536 101}
102
dd367a14 103// --- end of helpers --- begin functions ---
104
3dfa46a4 105void DrawOverview(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
dd367a14 106{
107 loadlibs();
108 if (!TFile::Open(fileName))
109 return;
110
111 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
112 if (!dNdEtaCorrection->LoadHistograms())
113 return;
114
115 dNdEtaCorrection->Finish();
116
117 dNdEtaCorrection->DrawOverview();
118}
119
8ca1a6d9 120void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
121{
122 gSystem->Load("libPWG0base");
123
124 TFile::Open(esdFile);
125 dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
126 fdNdEtaAnalysisESD->LoadHistograms();
127
128 TFile::Open(mcFile);
129 dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
130 fdNdEtaAnalysisMC->LoadHistograms();
0448e811 131 //fdNdEtaAnalysisMC->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
8ca1a6d9 132
133 for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
0448e811 134 fdNdEtaAnalysisESD->GetdNdEtaPtCutOffCorrectedHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(i));
8ca1a6d9 135
136 fdNdEtaAnalysisESD->DrawHistograms();
137}
138
139void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
140{
141 gSystem->Load("libPWG0base");
142
143 const char* ESDfolder = 0;
144
145 if (plot == 0) // all
146 ESDfolder = "dndeta";
147 else if (plot == 1) // mb
148 ESDfolder = "dndeta_mb";
149 else if (plot == 2) // mb vtx
150 ESDfolder = "dndeta_mbvtx";
151
152 TFile::Open("analysis_esd.root");
153 dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
154 fdNdEtaAnalysisESD->LoadHistograms();
155
156 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
157 dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
158
159 TH2F* hist = 0;
160
161 if (plot == 0) // all
162 hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetGeneratedHistogram();
163 else if (plot == 1) // mb
164 hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetMeasuredHistogram();
165 else if (plot == 2) // mb vtx
166 hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetMeasuredHistogram();
167
168 TH1* proj = hist->ProjectionX();
169
170 TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
171 for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
172 {
173 Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
174 if (value != 0)
175 {
176 printf("vtx = %f, esd = %f, corr = %f, ratio = %f\n", vertex->GetBinCenter(i), vertex->GetBinContent(i), value, vertex->GetBinContent(i) / value);
177 vertex->SetBinContent(i, vertex->GetBinContent(i) / value);
178 }
179 }
180
181 new TCanvas;
182 vertex->DrawCopy();
183}
184
185void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
186{
187 gSystem->Load("libPWG0base");
188
189 TFile::Open("analysis_esd.root");
190 dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
191 fdNdEtaAnalysisESD->LoadHistograms();
192
193 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
194 dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
195
196 //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
197 //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
198
199 TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
200 TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
201
202 TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
203
204 new TCanvas;
205 histESD->Draw();
206
207 new TCanvas;
208 histCorr->Draw();
209
210 for (Int_t x=1; x<=histESD->GetNbinsX(); ++x)
211 for (Int_t y=1; y<=histESD->GetNbinsY(); ++y)
212 for (Int_t z=1; z<=histESD->GetNbinsZ(); ++z)
213 {
214 Float_t value1 = histESD->GetBinContent(x, y, z);
215 Float_t value2 = histCorr->GetBinContent(histCorr->FindBin(histESD->GetXaxis()->GetBinCenter(x), histESD->GetYaxis()->GetBinCenter(y), histESD->GetZaxis()->GetBinCenter(z)));
216
217 if (value2 > 0 && value1 > 0)
218 {
219 printf("%f %f %f\n", value1, value2, value1 / value2);
220 diff->Fill(value1 / value2);
221 }
222 }
223
224 new TCanvas;
225 diff->Draw();
226}
227
dd367a14 228Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc = "")
229{
230 Double_t avgMC = 0;
231 Double_t avgESD = 0;
232 for (Int_t bin = histMC->FindBin(-0.7999); bin <= histMC->FindBin(0.7999); bin++)
233 {
234 avgMC += histMC->GetBinContent(bin);
235 avgESD += histESD->GetBinContent(bin);
236 }
237 Int_t nBins = histMC->FindBin(0.7999) - histMC->FindBin(-0.7999) + 1;
238
239 avgMC /= nBins;
240 avgESD /= nBins;
241
242 // deviation when integrate in |eta| < 0.8 between mc and esd
243 Double_t diffFullRange = (avgMC - avgESD) / avgMC;
244
245 Printf("%s: Integrated deviation in |eta| < 0.8 is %.2f %%", desc, diffFullRange * 100);
246
247 return diffFullRange;
248}
249
3dfa46a4 250void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
d09fb536 251{
252 TFile* file = TFile::Open("analysis_esd.root");
74fd10b3 253 TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
0fc41645 254 TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
74fd10b3 255 TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
256 TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
257 TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
258 TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
259 TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
9e952c39 260 TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
d09fb536 261
262 Prepare1DPlot(histESD);
0fc41645 263 Prepare1DPlot(histESDnsd);
d09fb536 264 Prepare1DPlot(histESDMB);
265 Prepare1DPlot(histESDMBVtx);
266
74fd10b3 267 Prepare1DPlot(histESDNoPt);
268 Prepare1DPlot(histESDMBNoPt);
269 Prepare1DPlot(histESDMBVtxNoPt);
9e952c39 270 Prepare1DPlot(histESDMBTracksNoPt);
74fd10b3 271
272 histESD->SetLineWidth(0);
0fc41645 273 histESDnsd->SetLineWidth(0);
74fd10b3 274 histESDMB->SetLineWidth(0);
275 histESDMBVtx->SetLineWidth(0);
276
277 histESDNoPt->SetLineWidth(0);
278 histESDMBNoPt->SetLineWidth(0);
279 histESDMBVtxNoPt->SetLineWidth(0);
72e597d7 280
9e952c39 281 histESD->SetMarkerColor(1);
0fc41645 282 histESDnsd->SetMarkerColor(6);
9e952c39 283 histESDMB->SetMarkerColor(2);
284 histESDMBVtx->SetMarkerColor(3);
72e597d7 285
3dfa46a4 286 histESD->SetLineColor(1);
0fc41645 287 histESDnsd->SetLineColor(6);
3dfa46a4 288 histESDMB->SetLineColor(2);
289 histESDMBVtx->SetLineColor(3);
290
9e952c39 291 histESDNoPt->SetMarkerColor(1);
292 histESDMBNoPt->SetMarkerColor(2);
293 histESDMBVtxNoPt->SetMarkerColor(3);
294 histESDMBTracksNoPt->SetMarkerColor(4);
74fd10b3 295
72e597d7 296 histESD->SetMarkerStyle(20);
0fc41645 297 histESDnsd->SetMarkerStyle(29);
72e597d7 298 histESDMB->SetMarkerStyle(21);
299 histESDMBVtx->SetMarkerStyle(22);
d09fb536 300
74fd10b3 301 histESDNoPt->SetMarkerStyle(20);
302 histESDMBNoPt->SetMarkerStyle(21);
303 histESDMBVtxNoPt->SetMarkerStyle(22);
9e952c39 304 histESDMBTracksNoPt->SetMarkerStyle(23);
3dfa46a4 305
745d6088 306 //Float_t etaLimit = 1.2999;
307 Float_t etaLimit = 2.01;
308 Float_t etaPlotLimit = 2.2;
4c351225 309
310 histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
311 histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
312 histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
0fc41645 313 histESDnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
74fd10b3 314
4c351225 315 histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
74fd10b3 316 histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
317 histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
9e952c39 318 histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
319
0fc41645 320 Float_t max = TMath::Max(histESDMBVtx->GetMaximum(), histESDMB->GetMaximum());
321 max = TMath::Max(max, histESD->GetMaximum());
322
745d6088 323 TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 0, max * 1.1);
0fc41645 324 Prepare1DPlot(dummy);
325 dummy->SetStats(kFALSE);
326 dummy->SetXTitle("#eta");
327 dummy->SetYTitle("dN_{ch}/d#eta");
328 dummy->GetYaxis()->SetTitleOffset(1);
329
0f67a57c 330 TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
d09fb536 331
0ab29cfa 332 dummy->DrawCopy();
d09fb536 333 histESDMBVtx->Draw("SAME");
334 histESDMB->Draw("SAME");
335 histESD->Draw("SAME");
336
3dfa46a4 337 if (save)
338 {
339 canvas->SaveAs("dNdEta1.gif");
340 canvas->SaveAs("dNdEta1.eps");
341 }
0ab29cfa 342
343 if (onlyESD)
344 return;
d09fb536 345
dd367a14 346 loadlibs();
74fd10b3 347
d09fb536 348 TFile* file2 = TFile::Open("analysis_mc.root");
74fd10b3 349 TH1* histMC = (TH1*) file2->Get("dndeta/dNdEta_corrected")->Clone("cloned");
0fc41645 350 TH1* histMCnsd = (TH1*) file2->Get("dndetaNSD/dNdEta_corrected");
745d6088 351 if (histMCnsd)
352 {
353 histMCnsd = (TH1*) histMCnsd->Clone("clonedNSD");
354 }
355 else
356 histMCnsd = new TH1F;
357
74fd10b3 358 TH1* histMCTr = (TH1*) file2->Get("dndetaTr/dNdEta_corrected")->Clone("cloned2");
359 TH1* histMCTrVtx = (TH1*) file2->Get("dndetaTrVtx/dNdEta_corrected")->Clone("cloned3");
d09fb536 360
0ab29cfa 361 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
d09fb536 362 fdNdEtaAnalysis->LoadHistograms();
3dfa46a4 363 dNdEtaAnalysis* fdNdEtaAnalysis2 = (dNdEtaAnalysis*) fdNdEtaAnalysis->Clone();
364
0fc41645 365 fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, "MC: full inelastic");
3dfa46a4 366 TH1* histMCPtCut = (TH1*) fdNdEtaAnalysis->GetdNdEtaHistogram(0)->Clone("histMCPtCut");
367
74fd10b3 368 fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
369 fdNdEtaAnalysis->LoadHistograms();
0fc41645 370 fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, "MC: minimum bias");
74fd10b3 371 TH1* histMCTrPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
372
373 fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
374 fdNdEtaAnalysis->LoadHistograms();
0fc41645 375 fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, "MC: MB with trigger");
74fd10b3 376 TH1* histMCTrVtxPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
d09fb536 377
9e952c39 378 fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
379 fdNdEtaAnalysis->LoadHistograms();
0fc41645 380 fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, "MC: Tracks w/o resolution effect");
9e952c39 381 TH1* histMCTracksPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
382
d09fb536 383 Prepare1DPlot(histMC);
0fc41645 384 Prepare1DPlot(histMCnsd);
74fd10b3 385 Prepare1DPlot(histMCTr);
386 Prepare1DPlot(histMCTrVtx);
387
d09fb536 388 Prepare1DPlot(histMCPtCut);
74fd10b3 389 Prepare1DPlot(histMCTrPtCut);
390 Prepare1DPlot(histMCTrVtxPtCut);
745d6088 391 Prepare1DPlot(histMCTracksPtCut);
392
393 histMC->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
394 histMCnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
395 histMCTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
396 histMCTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
397
398 histMCPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
399 histMCTrPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
400 histMCTrVtxPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
401 histMCTracksPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
d09fb536 402
9e952c39 403 histMC->SetLineColor(1);
0fc41645 404 histMCnsd->SetLineColor(6);
9e952c39 405 histMCTr->SetLineColor(2);
406 histMCTrVtx->SetLineColor(3);
74fd10b3 407
9e952c39 408 histMCPtCut->SetLineColor(1);
409 histMCTrPtCut->SetLineColor(2);
410 histMCTrVtxPtCut->SetLineColor(3);
411 if (histMCTracksPtCut)
412 histMCTracksPtCut->SetLineColor(4);
d09fb536 413
745d6088 414 TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
415
0ab29cfa 416 TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
0fc41645 417 dummy2->GetYaxis()->SetRangeUser(0, max * 1.1);
d09fb536 418
0ab29cfa 419 dummy2->DrawCopy();
d09fb536 420 histMC->Draw("SAME");
0fc41645 421 histMCnsd->Draw("SAME");
74fd10b3 422 histMCTr->Draw("SAME");
423 histMCTrVtx->Draw("SAME");
d09fb536 424 histESD->Draw("SAME");
0fc41645 425 histESDnsd->Draw("SAME");
74fd10b3 426 histESDMB->Draw("SAME");
427 histESDMBVtx->Draw("SAME");
d09fb536 428 histESDNoPt->Draw("SAME");
74fd10b3 429 histESDMBNoPt->Draw("SAME");
430 histESDMBVtxNoPt->Draw("SAME");
9e952c39 431 histESDMBTracksNoPt->Draw("SAME");
d09fb536 432 histMCPtCut->Draw("SAME");
74fd10b3 433 histMCTrPtCut->Draw("SAME");
434 histMCTrVtxPtCut->Draw("SAME");
9e952c39 435 if (histMCTracksPtCut)
436 histMCTracksPtCut->Draw("SAME");
d09fb536 437
3dfa46a4 438 if (save)
439 {
440 canvas2->SaveAs("dNdEta2.gif");
441 canvas2->SaveAs("dNdEta2.eps");
442 }
0ab29cfa 443
745d6088 444 DrawdNdEtaRatio(histESD, histMC, "full_inelastic", etaPlotLimit);
445 DrawdNdEtaRatio(histESDMB, histMCTr, "triggered", etaPlotLimit);
446 DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit);
447
0fc41645 448 new TCanvas;
449 dummy2->DrawCopy();
450 histMCnsd->Draw("SAME");
451 histESDnsd->Draw("SAME");
452
8ca1a6d9 453 TH1* ratio = (TH1*) histMC->Clone("ratio");
454 TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
455
456 ratio->Divide(histESD);
457 ratioNoPt->Divide(histESDNoPt);
458
3dfa46a4 459 ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
8ca1a6d9 460
461 ratio->SetLineColor(1);
462 ratioNoPt->SetLineColor(2);
463
dd367a14 464 Double_t average = 0; // average deviation from 1 in ratio (depends on the number of bins if statistical)
465 for (Int_t bin = ratio->FindBin(-0.7999); bin <= ratio->FindBin(0.7999); bin++)
466 average += TMath::Abs(ratio->GetBinContent(bin) - 1);
467 Int_t nBins = ratio->FindBin(0.7999) - ratio->FindBin(-0.7999) + 1;
468 average /= nBins;
469 Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
470
471 PrintIntegratedDeviation(histMC, histESD, "all events");
472 PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
473 PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
474 PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
475 PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
476 PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
477
9e952c39 478 TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 700, 600);
8ca1a6d9 479 canvas3->Range(0, 0, 1, 1);
480 //canvas3->Divide(1, 2, 0, 0);
481
482 //canvas3->cd(1);
9e952c39 483 TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.5, 0.98, 0.98);
8ca1a6d9 484 pad1->Draw();
485
9e952c39 486 TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.5);
8ca1a6d9 487 pad2->Draw();
488
489 pad1->SetRightMargin(0.05);
490 pad2->SetRightMargin(0.05);
0ab29cfa 491
8ca1a6d9 492 // no border between them
493 pad1->SetBottomMargin(0);
494 pad2->SetTopMargin(0);
495
496 pad1->cd();
497
9e952c39 498 TLegend* legend = new TLegend(0.4, 0.05, 0.65, 0.3);
72e597d7 499 legend->SetFillColor(0);
500 legend->AddEntry(histESDMBVtx, "triggered, vertex");
501 legend->AddEntry(histESDMB, "triggered");
502 legend->AddEntry(histESD, "all events");
503 legend->AddEntry(histMC, "MC prediction");
0ab29cfa 504
9e952c39 505 dummy->GetXaxis()->SetLabelSize(0.06);
506 dummy->GetYaxis()->SetLabelSize(0.06);
507 dummy->GetXaxis()->SetTitleSize(0.06);
508 dummy->GetYaxis()->SetTitleSize(0.06);
509 dummy->GetYaxis()->SetTitleOffset(0.7);
72e597d7 510 dummy->DrawCopy();
0ab29cfa 511 histESDMBVtx->Draw("SAME");
512 histESDMB->Draw("SAME");
513 histESD->Draw("SAME");
0ab29cfa 514 histMC->Draw("SAME");
72e597d7 515
516 legend->Draw();
0ab29cfa 517
8ca1a6d9 518 pad2->cd();
519 pad2->SetBottomMargin(0.15);
72e597d7 520
0fc41645 521 Float_t minR = TMath::Min(0.961, ratio->GetMinimum() * 0.95);
522 Float_t maxR = TMath::Max(1.049, ratio->GetMaximum() * 1.05);
9e952c39 523
745d6088 524 TH1F dummy3("dummy3", ";#eta;Ratio: MC / ESD", 100, -etaPlotLimit, etaPlotLimit);
8ca1a6d9 525 dummy3.SetStats(kFALSE);
745d6088 526 for (Int_t i=1; i<=100; ++i)
527 dummy3.SetBinContent(i, 1);
0fc41645 528 dummy3.GetYaxis()->SetRangeUser(minR, maxR);
8ca1a6d9 529 dummy3.SetLineWidth(2);
530 dummy3.GetXaxis()->SetLabelSize(0.06);
531 dummy3.GetYaxis()->SetLabelSize(0.06);
532 dummy3.GetXaxis()->SetTitleSize(0.06);
533 dummy3.GetYaxis()->SetTitleSize(0.06);
534 dummy3.GetYaxis()->SetTitleOffset(0.7);
535 dummy3.DrawCopy();
72e597d7 536
8ca1a6d9 537 ratio->Draw("SAME");
72e597d7 538
8ca1a6d9 539 //pad2->Draw();
72e597d7 540
8ca1a6d9 541 canvas3->Modified();
72e597d7 542
3dfa46a4 543 if (save)
544 {
545 canvas3->SaveAs("dNdEta.gif");
546 canvas3->SaveAs("dNdEta.eps");
547 }
8ca1a6d9 548
549 TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
72e597d7 550
551 ratio->Draw();
552 ratioNoPt->Draw("SAME");
553
4c351225 554 TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
555 legend->SetFillColor(0);
556 legend->AddEntry(ratio, "mc/esd");
557 legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
558 legend->Draw();
d09fb536 559}
560
745d6088 561void DrawdNdEtaRatio(TH1* corr, TH1* mc, const char* name, Float_t etaPlotLimit)
562{
563 TCanvas* canvas3 = new TCanvas(name, name, 700, 600);
564 canvas3->Range(0, 0, 1, 1);
565
566 TPad* pad1 = new TPad(Form("%s_1", name), "", 0, 0.5, 0.98, 0.98);
567 pad1->Draw();
568
569 TPad* pad2 = new TPad(Form("%s_2", name), "", 0, 0.02, 0.98, 0.5);
570 pad2->Draw();
571
572 pad1->SetRightMargin(0.05);
573 pad2->SetRightMargin(0.05);
574
575 // no border between them
576 pad1->SetBottomMargin(0);
577 pad2->SetTopMargin(0);
578
579 pad1->cd();
580
581 TLegend* legend = new TLegend(0.4, 0.05, 0.65, 0.3);
582 legend->SetFillColor(0);
583 legend->AddEntry(corr, "corrected");
584 legend->AddEntry(mc, "MC prediction");
585
586 TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 0, corr->GetMaximum() * 1.1);
587 Prepare1DPlot(dummy);
588 dummy->SetStats(kFALSE);
589 dummy->SetXTitle("#eta");
590 dummy->SetYTitle("dN_{ch}/d#eta");
591 dummy->GetYaxis()->SetTitleOffset(1);
592
593 dummy->GetXaxis()->SetLabelSize(0.06);
594 dummy->GetYaxis()->SetLabelSize(0.06);
595 dummy->GetXaxis()->SetTitleSize(0.06);
596 dummy->GetYaxis()->SetTitleSize(0.06);
597 dummy->GetYaxis()->SetTitleOffset(0.7);
598 dummy->DrawCopy();
599
600 corr->Draw("SAME");
601 mc->Draw("SAME");
602
603 legend->Draw();
604
605 pad2->cd();
606 pad2->SetBottomMargin(0.15);
607
608 TH1* ratio = (TH1*) mc->Clone("ratio");
609 ratio->Divide(corr);
610
611 Float_t minR = TMath::Min(0.961, ratio->GetMinimum() * 0.95);
612 Float_t maxR = TMath::Max(1.049, ratio->GetMaximum() * 1.05);
613
614 TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
615 dummy3.SetStats(kFALSE);
616 for (Int_t i=1; i<=100; ++i)
617 dummy3.SetBinContent(i, 1);
618 dummy3.GetYaxis()->SetRangeUser(minR, maxR);
619 dummy3.SetLineWidth(2);
620 dummy3.GetXaxis()->SetLabelSize(0.06);
621 dummy3.GetYaxis()->SetLabelSize(0.06);
622 dummy3.GetXaxis()->SetTitleSize(0.06);
623 dummy3.GetYaxis()->SetTitleSize(0.06);
624 dummy3.GetYaxis()->SetTitleOffset(0.7);
625 dummy3.DrawCopy();
626
627 ratio->Draw("SAME");
628
629 canvas3->Modified();
630}
631
d09fb536 632void ptSpectrum()
633{
634 TFile* file = TFile::Open("analysis_esd.root");
0ab29cfa 635 TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
d09fb536 636
637 TFile* file2 = TFile::Open("analysis_mc.root");
0ab29cfa 638 TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
d09fb536 639
640 TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
641 InitPad();
642 gPad->SetLogy();
643
644 Prepare1DPlot(histMC);
645 Prepare1DPlot(histESD);
646
647 histESD->SetTitle("");
25db2d85 648 histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
649 histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
d09fb536 650
651 histMC->SetLineColor(kBlue);
652 histESD->SetLineColor(kRed);
653
654 histESD->GetYaxis()->SetTitleOffset(1.5);
655 histESD->GetXaxis()->SetRangeUser(0, 4.9999);
656
657 histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
658
659 histESD->Draw();
660 histMC->Draw("SAME");
661
662 canvas->SaveAs("ptSpectrum.gif");
0ab29cfa 663 canvas->SaveAs("ptSpectrum.eps");
92d2d8ad 664}
665
666void ptCutoff()
667{
668 gSystem->Load("libPWG0base");
669
0448e811 670 TFile::Open("correction_map.root");
8b3563f4 671 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
0448e811 672 dNdEtaCorrection->LoadHistograms();
92d2d8ad 673
0448e811 674 dNdEtaCorrection->GetMeasuredFraction(AlidNdEtaCorrection::kINEL, 0.3, -100, kTRUE);
92d2d8ad 675
0448e811 676 TH1* hist = dynamic_cast<TH1*> (gROOT->FindObject("generated_pt")->Clone("ptcutoff"));
92d2d8ad 677
678 hist->GetXaxis()->SetRangeUser(0, 0.9999);
679 hist->SetMinimum(0);
92d2d8ad 680
681 hist->SetTitle("Generated Particles");
d09fb536 682 Prepare1DPlot(hist);
92d2d8ad 683
72e597d7 684 TCanvas* canvas = new TCanvas("ptCutoff", "ptCutoff", 700, 500);
685 hist->DrawCopy();
92d2d8ad 686
d09fb536 687 TLine* line = new TLine(0.3, 0 - hist->GetMaximum() * 0, 0.3, hist->GetMaximum() * 1.1);
92d2d8ad 688 line->SetLineWidth(3);
d09fb536 689 line->SetLineColor(kRed);
92d2d8ad 690 line->Draw();
691
692 canvas->SaveAs("ptCutoff.gif");
0ab29cfa 693 canvas->SaveAs("ptCutoff.eps");
72e597d7 694
0448e811 695 TH1F* factor = new TH1F("factor", ";#eta;correction factor", 20, -1, 1.000001);
696 factor->SetLineWidth(2);
697 for (Float_t eta = -0.95; eta<1; eta += 0.1)
698 factor->Fill(eta, 1.0 / dNdEtaCorrection->GetMeasuredFraction(AlidNdEtaCorrection::kINEL, 0.3, eta, kFALSE));
72e597d7 699
700 TCanvas* canvas = new TCanvas("ptCutoff_factor", "ptCutoff_factor", 700, 500);
701 InitPad();
702
703 Prepare1DPlot(factor);
704 factor->GetYaxis()->SetRangeUser(1, 2);
705 factor->GetYaxis()->SetTitleOffset(1);
706 factor->Draw();
707
708 canvas->SaveAs("ptCutoff_factor.eps");
92d2d8ad 709}
710
0bd1f8a0 711void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
92d2d8ad 712{
0448e811 713 gSystem->Load("libPWG0base");
714
715 TFile::Open(fileName);
716 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
717 dNdEtaCorrection->LoadHistograms();
0ab29cfa 718
0448e811 719 TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
720 TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
0ab29cfa 721
722 Prepare2DPlot(corrTrigger);
0448e811 723 corrTrigger->SetTitle("b) Trigger bias correction");
0ab29cfa 724
725 Prepare2DPlot(corrVtx);
0448e811 726 corrVtx->SetTitle("a) Vertex reconstruction correction");
0ab29cfa 727
5c495d37 728 corrTrigger->GetYaxis()->SetTitle("Multiplicity");
729 corrVtx->GetYaxis()->SetTitle("Multiplicity");
730
0ab29cfa 731 TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
732 canvas->Divide(2, 1);
733
734 canvas->cd(1);
735 InitPadCOLZ();
0448e811 736 corrVtx->DrawCopy("COLZ");
0ab29cfa 737
738 canvas->cd(2);
739 InitPadCOLZ();
0448e811 740 corrTrigger->DrawCopy("COLZ");
0ab29cfa 741
742 canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
743 canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
744
745 canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
746 canvas->Divide(2, 1);
747
748 corrTrigger->GetYaxis()->SetRangeUser(0, 5);
749 corrVtx->GetYaxis()->SetRangeUser(0, 5);
750
751 canvas->cd(1);
752 InitPadCOLZ();
0448e811 753 corrVtx->DrawCopy("COLZ");
0ab29cfa 754
755 canvas->cd(2);
756 InitPadCOLZ();
0448e811 757 corrTrigger->DrawCopy("COLZ");
0ab29cfa 758
759 canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
760 canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
761}
762
763void TriggerBias(const char* fileName = "correction_map.root")
764{
765 TFile* file = TFile::Open(fileName);
92d2d8ad 766
5c495d37 767 TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
92d2d8ad 768
769 Prepare2DPlot(corr);
770 corr->SetTitle("Trigger bias correction");
771
772 TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
773 InitPadCOLZ();
774 corr->DrawCopy("COLZ");
775
25db2d85 776 canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
0ab29cfa 777 canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
92d2d8ad 778
779 corr->GetYaxis()->SetRangeUser(0, 5);
780
781 canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
782 InitPadCOLZ();
783 corr->DrawCopy("COLZ");
784
25db2d85 785 canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
0ab29cfa 786 canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
92d2d8ad 787}
788
72e597d7 789void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
790{
791 gSystem->Load("libPWG0base");
792
793 TFile* file = TFile::Open(fileName);
794 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
0448e811 795 dNdEtaCorrection->LoadHistograms();
72e597d7 796
0448e811 797 TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("x");
798 TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
72e597d7 799
0448e811 800 TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
8ca1a6d9 801 canvas->Divide(2, 1);
802
803 canvas->cd(1);
72e597d7 804 InitPad();
805
806 Prepare1DPlot(hist);
807 hist->SetTitle("");
808 hist->GetYaxis()->SetTitle("correction factor");
809 hist->GetYaxis()->SetRangeUser(1, 1.5);
810 hist->GetYaxis()->SetTitleOffset(1.6);
811 hist->Draw();
812
8ca1a6d9 813 canvas->cd(2);
814 InitPad();
815
816 Prepare1DPlot(hist2);
817 hist2->SetTitle("");
818 hist2->GetYaxis()->SetTitle("correction factor");
819 hist2->GetXaxis()->SetRangeUser(0, 5);
820 hist2->GetYaxis()->SetTitleOffset(1.6);
821 hist2->GetXaxis()->SetTitle("multiplicity");
822 hist2->Draw();
823
824 TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
825 pave->SetFillColor(0);
826 pave->AddText("|z| < 10 cm");
827 pave->Draw();
828
72e597d7 829 canvas->SaveAs("TriggerBias1D.eps");
830}
831
92d2d8ad 832void VtxRecon()
833{
834 TFile* file = TFile::Open("correction_map.root");
835
72e597d7 836 TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
92d2d8ad 837
838 Prepare2DPlot(corr);
839 corr->SetTitle("Vertex reconstruction correction");
840
841 TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
842 InitPadCOLZ();
25db2d85 843 corr->DrawCopy("COLZ");
844
0ab29cfa 845 canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
846 canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
25db2d85 847
848 corr->GetYaxis()->SetRangeUser(0, 5);
849
0ab29cfa 850 canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
25db2d85 851 InitPadCOLZ();
852 corr->DrawCopy("COLZ");
92d2d8ad 853
25db2d85 854 canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
0ab29cfa 855 canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
92d2d8ad 856}
857
72e597d7 858void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
859{
860 gSystem->Load("libPWG0base");
861
862 TFile* file = TFile::Open(fileName);
863 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
0448e811 864 dNdEtaCorrection->LoadHistograms();
72e597d7 865
0448e811 866 TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
867 TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
8ca1a6d9 868
869 TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
870 canvas->Divide(2, 1);
72e597d7 871
8ca1a6d9 872 canvas->cd(1);
72e597d7 873 InitPad();
874
875 Prepare1DPlot(hist);
876 hist->SetTitle("");
877 hist->GetYaxis()->SetTitle("correction factor");
878 hist->GetYaxis()->SetRangeUser(1, 1.8);
879 hist->GetYaxis()->SetTitleOffset(1.6);
8ca1a6d9 880 hist->DrawCopy();
881
882 canvas->cd(2);
883 InitPad();
884
885 Prepare1DPlot(hist2);
886 hist2->SetTitle("");
887 hist2->GetYaxis()->SetTitle("correction factor");
888 hist2->GetXaxis()->SetRangeUser(0, 20);
889 hist2->GetYaxis()->SetTitleOffset(1.6);
890 hist2->GetXaxis()->SetTitle("multiplicity");
891 hist2->Draw();
892
893 TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
894 pave->SetFillColor(0);
895 pave->AddText("|z| < 10 cm");
896 pave->Draw();
72e597d7 897
898 canvas->SaveAs("VtxRecon1D.eps");
8ca1a6d9 899
0448e811 900 Correction1DCreatePlots(fileName, folderName, 9.9, 2);
8ca1a6d9 901
0448e811 902 TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
903 TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
904
905 Prepare1DPlot(corrX);
906 Prepare1DPlot(corrZ);
907
908 corrX->GetYaxis()->SetTitleOffset(1.5);
909 corrZ->GetYaxis()->SetTitleOffset(1.5);
910
911 corrX->SetTitle("a) z projection");
912 corrZ->SetTitle("b) p_{T} projection");
913
914 corrX->GetYaxis()->SetTitle("correction factor");
915 corrZ->GetYaxis()->SetTitle("correction factor");
916
917 corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
918
919 TString canvasName;
920 canvasName.Form("VtxRecon1D_Track");
921 TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
922 canvas->Divide(2, 1);
923
924 canvas->cd(1);
925 InitPad();
926 corrX->DrawCopy();
927
928 canvas->cd(2);
929 InitPad();
930 gPad->SetLogx();
931 corrZ->Draw();
932
933 canvas->SaveAs("VtxRecon1D_Track.eps");
934 canvas->SaveAs("VtxRecon1D_Track.gif");
72e597d7 935}
936
0ab29cfa 937void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
1afae8ff 938{
25db2d85 939 gSystem->Load("libPWG0base");
940
0ab29cfa 941 TFile::Open(fileName);
8b3563f4 942 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
0ab29cfa 943 dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
944
945 TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
946 TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
947
948 TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
949 TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
950
951 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
952 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
953 gene->GetXaxis()->SetRangeUser(-10, 10);
954 meas->GetXaxis()->SetRangeUser(-10, 10);
955
956 Float_t eff1 = gene->Integral() / meas->Integral();
957 Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
958
959 printf("Correction without pT cut: %f +- %f\n", eff1, error1);
960
961 gene->GetZaxis()->SetRangeUser(0.3, 10);
962 meas->GetZaxis()->SetRangeUser(0.3, 10);
963
964 Float_t eff2 = gene->Integral() / meas->Integral();
965 Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
966
967 printf("Correction with pT cut: %f +- %f\n", eff2, error2);
968
969 gene->GetZaxis()->SetRangeUser(0.3, 1);
970 meas->GetZaxis()->SetRangeUser(0.3, 1);
971
972 Float_t eff3 = gene->Integral() / meas->Integral();
973 Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
974
975 printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
976}
977
0448e811 978void Correction1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType = 0)
0ab29cfa 979{
980 TFile::Open(fileName);
bdfe2916 981 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
0448e811 982 dNdEtaCorrection->LoadHistograms();
0ab29cfa 983
0448e811 984 TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
985 TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetMeasuredHistogram();
0ab29cfa 986
987 gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
988 meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
989 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
990 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
991 AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kTRUE);
992 gene->GetYaxis()->SetRange(0, 0);
993 meas->GetYaxis()->SetRange(0, 0);
994
995 gene->GetXaxis()->SetRangeUser(-10, 10);
996 meas->GetXaxis()->SetRangeUser(-10, 10);
997 AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kTRUE);
998 gene->GetZaxis()->SetRange(0, 0);
999 meas->GetZaxis()->SetRange(0, 0);
1000
1001 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1002 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1003 AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kTRUE);
1004}
25db2d85 1005
0448e811 1006void Correction1D(Int_t correctionType = 0, const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
1007{
1008 gSystem->Load("libPWG0base");
1009
1010 Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType);
1011
1012 TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1013 TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1014 TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1015
1016 Prepare1DPlot(corrX);
1017 Prepare1DPlot(corrY);
1018 Prepare1DPlot(corrZ);
1019
1020 corrX->SetTitle("a) z projection");
1021 corrY->SetTitle("b) #eta projection");
1022 corrZ->SetTitle("c) p_{T} projection");
1023
1024 corrX->GetYaxis()->SetTitle("correction factor");
1025 corrY->GetYaxis()->SetTitle("correction factor");
1026 corrZ->GetYaxis()->SetTitle("correction factor");
1027
1028 corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1029
1030 TString canvasName;
1031 canvasName.Form("Correction1D_%s", folder);
1032 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1033 canvas->Divide(3, 1);
1034
1035 canvas->cd(1);
1036 InitPad();
1037 corrX->DrawCopy();
1038
1039 canvas->cd(2);
1040 InitPad();
1041 corrY->Draw();
1042
1043 canvas->cd(3);
1044 InitPad();
1045 corrZ->Draw();
1046
1047 canvas->SaveAs(Form("Correction1D_%d_%s_%f.gif", correctionType, fileName, upperPtLimit));
1048 canvas->SaveAs(Form("Correction1D_%d_%s_%f.eps", correctionType, fileName, upperPtLimit));
1049}
1050
72e597d7 1051void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
0ab29cfa 1052{
1053 gSystem->Load("libPWG0base");
1054
0448e811 1055 Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
0ab29cfa 1056
0448e811 1057 TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1058 TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1059 TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
0ab29cfa 1060
1061 Prepare1DPlot(corrX);
1062 Prepare1DPlot(corrY);
1063 Prepare1DPlot(corrZ);
1064
0448e811 1065 corrX->SetTitle("a) z projection");
72e597d7 1066 corrY->SetTitle("a) #eta projection");
1067 corrZ->SetTitle("b) p_{T} projection");
1068
1069 corrY->GetYaxis()->SetTitle("correction factor");
1070 corrZ->GetYaxis()->SetTitle("correction factor");
0ab29cfa 1071
1072 corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1073
bdfe2916 1074 TString canvasName;
72e597d7 1075 canvasName.Form("Track2Particle1D_%s", folder);
bdfe2916 1076 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
0ab29cfa 1077 canvas->Divide(3, 1);
1078
1079 canvas->cd(1);
1080 InitPad();
bdfe2916 1081 corrX->DrawCopy();
0ab29cfa 1082
1083 canvas->cd(2);
1084 InitPad();
0448e811 1085 corrY->Draw();
0ab29cfa 1086
1087 canvas->cd(3);
1088 InitPad();
0448e811 1089 corrZ->Draw();
0ab29cfa 1090
72e597d7 1091 canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
1092 canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
1093
5c495d37 1094 //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
1095
72e597d7 1096 canvasName.Form("Track2Particle1D_%s_etapt", folder);
1097 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
1098 canvas->Divide(2, 1);
1099
1100 canvas->cd(1);
1101 InitPad();
1102 corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
1103 corrY->GetYaxis()->SetRangeUser(1, 1.5);
1104 corrY->GetYaxis()->SetTitleOffset(1.5);
1105 corrY->DrawCopy();
5c495d37 1106 TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
1107 pave->AddText("|z| < 10 cm");
1108 pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
1109 pave->Draw();
72e597d7 1110
1111 canvas->cd(2);
1112 InitPad();
5c495d37 1113 gPad->SetLogx();
1114 corrZ->GetYaxis()->SetRangeUser(1, 2.5);
1115 corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
72e597d7 1116 corrZ->GetYaxis()->SetTitleOffset(1.5);
1117 corrZ->DrawCopy();
5c495d37 1118 pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
1119 pave->AddText("|z| < 10 cm");
1120 pave->AddText("|#eta| < 0.8");
1121 pave->Draw();
72e597d7 1122
1123 canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
4c351225 1124 canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
0ab29cfa 1125}
1126
1127void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
1128{
1129 gSystem->Load("libPWG0base");
1130
8ca1a6d9 1131 // particle type
1132 for (Int_t particle=0; particle<4; ++particle)
1133 {
1134 TString dirName;
1135 dirName.Form("correction_%d", particle);
1136 Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
0ab29cfa 1137
8ca1a6d9 1138 TString tmpx, tmpy, tmpz;
1139 tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
1140 tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
1141 tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
0ab29cfa 1142
8ca1a6d9 1143 TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
1144 TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
1145 TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
0ab29cfa 1146
8ca1a6d9 1147 Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
0ab29cfa 1148
8ca1a6d9 1149 TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
1150 TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
1151 TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
0ab29cfa 1152
8ca1a6d9 1153 posX->Divide(negX);
1154 posY->Divide(negY);
1155 posZ->Divide(negZ);
0ab29cfa 1156
8ca1a6d9 1157 Prepare1DPlot(posX);
1158 Prepare1DPlot(posY);
1159 Prepare1DPlot(posZ);
0ab29cfa 1160
8ca1a6d9 1161 Float_t min = 0.8;
1162 Float_t max = 1.2;
0ab29cfa 1163
8ca1a6d9 1164 posX->SetMinimum(min);
1165 posX->SetMaximum(max);
1166 posY->SetMinimum(min);
1167 posY->SetMaximum(max);
1168 posZ->SetMinimum(min);
1169 posZ->SetMaximum(max);
0ab29cfa 1170
8ca1a6d9 1171 posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
0ab29cfa 1172
8ca1a6d9 1173 posX->GetYaxis()->SetTitleOffset(1.7);
1174 posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
1175 posY->GetYaxis()->SetTitleOffset(1.7);
1176 posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
1177 posZ->GetYaxis()->SetTitleOffset(1.7);
1178 posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
0ab29cfa 1179
8ca1a6d9 1180 posZ->GetXaxis()->SetRangeUser(0, 1);
0ab29cfa 1181
8ca1a6d9 1182 TString canvasName;
1183 canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
0ab29cfa 1184
8ca1a6d9 1185 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1186 canvas->Divide(3, 1);
0ab29cfa 1187
8ca1a6d9 1188 canvas->cd(1);
1189 InitPad();
1190 posX->DrawCopy();
0ab29cfa 1191
8ca1a6d9 1192 canvas->cd(2);
1193 InitPad();
1194 posY->DrawCopy();
1195
1196 canvas->cd(3);
1197 InitPad();
1198 posZ->DrawCopy();
0ab29cfa 1199
8ca1a6d9 1200 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1201 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1202 }
0ab29cfa 1203}
1204
1205void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
1206{
1207 TFile::Open(fileName);
1208 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
0448e811 1209 dNdEtaCorrection->LoadHistograms();
0ab29cfa 1210
0448e811 1211 TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1212 TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1afae8ff 1213
25db2d85 1214 gene->GetZaxis()->SetRangeUser(0.3, 10);
1215 meas->GetZaxis()->SetRangeUser(0.3, 10);
1216 AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
0ab29cfa 1217 gene->GetZaxis()->SetRange(0, 0);
1218 meas->GetZaxis()->SetRange(0, 0);
25db2d85 1219
1220 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1221 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1222 AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
0ab29cfa 1223 gene->GetYaxis()->SetRange(0, 0);
1224 meas->GetYaxis()->SetRange(0, 0);
25db2d85 1225
1226 gene->GetXaxis()->SetRangeUser(-10, 10);
1227 meas->GetXaxis()->SetRangeUser(-10, 10);
1228 AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
0ab29cfa 1229 gene->GetXaxis()->SetRange(0, 0);
1230 meas->GetXaxis()->SetRange(0, 0);
1231}
1232
0bd1f8a0 1233void Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
0ab29cfa 1234{
1235 gSystem->Load("libPWG0base");
1236
1237 Track2Particle2DCreatePlots(fileName);
25db2d85 1238
0448e811 1239 TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1240 TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1241 TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
25db2d85 1242
1243 /* this reads them from the file
d09fb536 1244 TH2* corrYX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx"));
1245 TH2* corrZX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx"));
25db2d85 1246 TH2* corrZY = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy"));*/
1afae8ff 1247
1248 Prepare2DPlot(corrYX);
1249 Prepare2DPlot(corrZX);
1250 Prepare2DPlot(corrZY);
1251
0ab29cfa 1252 const char* title = "";
92d2d8ad 1253 corrYX->SetTitle(title);
1254 corrZX->SetTitle(title);
1255 corrZY->SetTitle(title);
1afae8ff 1256
d09fb536 1257 TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1afae8ff 1258 canvas->Divide(3, 1);
1259
1260 canvas->cd(1);
1261 InitPadCOLZ();
1262 corrYX->Draw("COLZ");
1263
1264 canvas->cd(2);
1265 InitPadCOLZ();
1266 corrZX->Draw("COLZ");
1267
1268 canvas->cd(3);
1269 InitPadCOLZ();
1270 corrZY->Draw("COLZ");
92d2d8ad 1271
25db2d85 1272 canvas->SaveAs(Form("Track2Particle2D_%d.gif", gMax));
0ab29cfa 1273 canvas->SaveAs(Form("Track2Particle2D_%d.eps", gMax));
1274}
1275
1276void CompareTrack2Particle2D()
1277{
1278 gSystem->Load("libPWG0base");
1279
1280 Track2Particle2DCreatePlots("correction_maponly-positive.root");
1281
1282 TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("pos_yx"));
1283 TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("pos_zx"));
1284 TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("pos_zy"));
1285
1286 Track2Particle2DCreatePlots("correction_maponly-negative.root");
1287
1288 TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("neg_yx"));
1289 TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("neg_zx"));
1290 TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("neg_zy"));
1291
1292 posYX->Divide(negYX);
1293 posZX->Divide(negZX);
1294 posZY->Divide(negZY);
1295
1296 Prepare2DPlot(posYX);
1297 Prepare2DPlot(posZX);
1298 Prepare2DPlot(posZY);
1299
1300 Float_t min = 0.8;
1301 Float_t max = 1.2;
1302
1303 posYX->SetMinimum(min);
1304 posYX->SetMaximum(max);
1305 posZX->SetMinimum(min);
1306 posZX->SetMaximum(max);
1307 posZY->SetMinimum(min);
1308 posZY->SetMaximum(max);
1309
1310 TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1311 canvas->Divide(3, 1);
1312
1313 canvas->cd(1);
1314 InitPadCOLZ();
1315 posYX->Draw("COLZ");
1316
1317 canvas->cd(2);
1318 InitPadCOLZ();
1319 posZX->Draw("COLZ");
1320
1321 canvas->cd(3);
1322 InitPadCOLZ();
1323 posZY->Draw("COLZ");
1324
1325 canvas->SaveAs("CompareTrack2Particle2D.gif");
1326 canvas->SaveAs("CompareTrack2Particle2D.eps");
1afae8ff 1327}
1328
1329void Track2Particle3D()
1330{
1331 // get left margin proper
1332
1333 TFile* file = TFile::Open("correction_map.root");
1334
d09fb536 1335 TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1336
1337 corr->SetTitle("Correction Factor");
1338 SetRanges(corr->GetZaxis());
1339
1340 Prepare3DPlot(corr);
1341
1342 TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1343 canvas->SetTheta(29.428);
1344 canvas->SetPhi(16.5726);
1345
1346 corr->Draw();
1347
1348 canvas->SaveAs("Track2Particle3D.gif");
0ab29cfa 1349 canvas->SaveAs("Track2Particle3D.eps");
d09fb536 1350}
1351
1352void Track2Particle3DAll()
1353{
d09fb536 1354 TFile* file = TFile::Open("correction_map.root");
1355
1afae8ff 1356 TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1357 TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1358 TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1359
1360 gene->SetTitle("Generated Particles");
1361 meas->SetTitle("Measured Tracks");
1362 corr->SetTitle("Correction Factor");
1363
1364 Prepare3DPlot(gene);
1365 Prepare3DPlot(meas);
1366 Prepare3DPlot(corr);
1367
d09fb536 1368 TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1afae8ff 1369 canvas->Divide(3, 1);
1370
1371 canvas->cd(1);
1372 InitPad();
1373 gene->Draw();
1374
1375 canvas->cd(2);
1376 meas->Draw();
1377
1378 canvas->cd(3);
1379 corr->Draw();
d09fb536 1380
1381 canvas->SaveAs("Track2Particle3DAll.gif");
0ab29cfa 1382 canvas->SaveAs("Track2Particle3DAll.eps");
1afae8ff 1383}
1384
6b7fa615 1385void MultiplicityMC(Int_t xRangeMax = 50)
4c6b34a8 1386{
1387 TFile* file = TFile::Open("multiplicityMC.root");
1388
1389 if (!file)
1390 {
1391 printf("multiplicityMC.root could not be opened.\n");
1392 return;
1393 }
1394
1395 TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1396 TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1397 TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1398
1399 TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1400 TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1401 //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1402 for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1403 {
1404 TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1405 proj->Fit("gaus", "0");
1406 correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1407 correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1408
1409 continue;
1410
1411 // draw for debugging
1412 new TCanvas;
1413 proj->DrawCopy();
1414 proj->GetFunction("gaus")->DrawCopy("SAME");
1415 }
1416
1417 TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1418
1419 for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1420 {
1421 Float_t mean = correction->GetBinContent(i);
1422 Float_t width = correctionWidth->GetBinContent(i);
1423
1424 Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
1425 Int_t fillEnd = fMultiplicityESDCorrected->FindBin(mean + width * 3);
1426 printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1427
1428 for (Int_t j=fillBegin; j <= fillEnd; ++j)
1429 {
1430 fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
1431 }
1432 }
1433
1434 TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
1435 fMultiplicityESDCorrectedRebinned->Rebin(10);
1436 fMultiplicityESDCorrectedRebinned->Scale(0.1);
1437
6b7fa615 1438 TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
1439 ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
1440 ratio->Divide(fMultiplicityMC);
1441
4c6b34a8 1442 TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
1443 ratio2->Divide(fMultiplicityMC);
1444
1445 TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
1446 canvas->Divide(3, 2);
1447
6b7fa615 1448 fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
1449 ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
1450 fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
1451 fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
1452 correction->GetXaxis()->SetRangeUser(0, xRangeMax);
1453 fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
1454 fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
1455
1456 canvas->cd(1); //InitPad();
4c6b34a8 1457 fMultiplicityESD->Draw();
1458 fMultiplicityMC->SetLineColor(2);
1459 fMultiplicityMC->Draw("SAME");
1460
6b7fa615 1461 TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1462 legend->AddEntry(fMultiplicityESD, "ESD");
1463 legend->AddEntry(fMultiplicityMC, "MC");
1464 legend->Draw();
4c6b34a8 1465
6b7fa615 1466 canvas->cd(2);
4c6b34a8 1467 fCorrelation->Draw("COLZ");
1468
6b7fa615 1469 canvas->cd(3);
4c6b34a8 1470 correction->Draw();
1471 //correction->Fit("pol1");
1472 correctionWidth->SetLineColor(2);
1473 correctionWidth->Draw("SAME");
1474
6b7fa615 1475 legend = new TLegend(0.2, 0.7, 0.45, 0.85);
1476 legend->AddEntry(correction, "#bar{x}");
1477 legend->AddEntry(correctionWidth, "#sigma");
1478 legend->Draw();
1479
1480 canvas->cd(4);
1481 ratio->Draw();
1482
1483 ratio2->SetLineColor(2);
1484 ratio2->Draw("SAME");
1485
1486 legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1487 legend->AddEntry(ratio, "uncorrected");
1488 legend->AddEntry(ratio2, "corrected");
1489 legend->Draw();
1490
4c6b34a8 1491 canvas->cd(5);
6b7fa615 1492 fMultiplicityESDCorrected->SetLineColor(kBlue);
4c6b34a8 1493 fMultiplicityESDCorrected->Draw();
1494 fMultiplicityMC->Draw("SAME");
1495 fMultiplicityESD->Draw("SAME");
1496
6b7fa615 1497 legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1498 legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1499 legend->AddEntry(fMultiplicityMC, "MC");
1500 legend->AddEntry(fMultiplicityESD, "ESD");
1501 legend->Draw();
1502
4c6b34a8 1503 canvas->cd(6);
6b7fa615 1504 fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
4c6b34a8 1505 fMultiplicityESDCorrectedRebinned->Draw();
1506 fMultiplicityMC->Draw("SAME");
6b7fa615 1507
1508 legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1509 legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
1510 legend->AddEntry(fMultiplicityMC, "MC");
1511 legend->Draw();
1512
1513 canvas->SaveAs("MultiplicityMC.gif");
4c6b34a8 1514}
1515
1516void MultiplicityESD()
1517{
1518 TFile* file = TFile::Open("multiplicityESD.root");
1519
1520 if (!file)
1521 {
1522 printf("multiplicityESD.root could not be opened.\n");
1523 return;
1524 }
1525
1526 TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
1527
1528 TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
1529
1530 fMultiplicityESD->Draw();
1531}
1532
0ab29cfa 1533void drawPlots(Int_t max)
1afae8ff 1534{
0ab29cfa 1535 gMax = max;
1afae8ff 1536
0ab29cfa 1537 ptCutoff();
1538 TriggerBias();
1539 VtxRecon();
1540 Track2Particle2D();
1541 Track2Particle3D();
1542 ptSpectrum();
1543 dNdEta();
1afae8ff 1544}
1545
0ab29cfa 1546void drawPlots()
d09fb536 1547{
0ab29cfa 1548 drawPlots(5);
1549 drawPlots(2);
1afae8ff 1550}
770a1f1d 1551
1552void CompareCorrection2Measured(const char* dataInput = "analysis_esd_raw.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1553{
1554 loadlibs();
1555
1556 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1557 TFile::Open(correctionMapFile);
1558 dNdEtaCorrection->LoadHistograms();
1559
1560 TFile* file = TFile::Open(dataInput);
1561
1562 if (!file)
1563 {
1564 cout << "Error. File not found" << endl;
1565 return;
1566 }
1567
1568 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1569 fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1570
1571 gROOT->cd();
1572
1573 TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
1574 hist1->SetTitle("mc");
3dfa46a4 1575 Printf("mc contains %f entries", hist1->Integral());
1576 Printf("mc 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()));
770a1f1d 1577
1578 TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
1579 hist2->SetTitle("esd");
3dfa46a4 1580 Printf("esd contains %f entries", hist2->Integral());
1581 Printf("esd 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()));
770a1f1d 1582
1583 AliPWG0Helper::CreateDividedProjections(hist1, hist2);
0fc41645 1584 AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
1585
1586 hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1587 hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1588 AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
770a1f1d 1589
1590 new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
1591 new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
1592 new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
0fc41645 1593 new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
1594 new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
3dfa46a4 1595}
1596
1597void CompareMeasured2Measured(const char* dataInput = "analysis_esd_raw.root", const char* dataInput2 = "analysis_esd_raw.root")
1598{
1599 loadlibs();
1600
1601 TFile* file = TFile::Open(dataInput);
1602
1603 if (!file)
1604 {
1605 cout << "Error. File not found" << endl;
1606 return;
1607 }
1608
1609 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1610 fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1611
1612 TFile* file = TFile::Open(dataInput2);
1613
1614 if (!file)
1615 {
1616 cout << "Error. File not found" << endl;
1617 return;
1618 }
1619
1620 dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta2", "dndeta2");
1621 fdNdEtaAnalysis2->LoadHistograms("fdNdEtaAnalysisESD");
1622
1623 gROOT->cd();
1624
1625 TH3* hist1 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd1");
1626 hist1->SetTitle("esd1");
1627 Printf("esd1 contains %f entries", hist1->GetEntries());
1628 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()));
1629
1630 TH3* hist2 = (TH3*) fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd2");
1631 hist2->SetTitle("esd2");
1632 Printf("esd2 contains %f entries", hist2->GetEntries());
1633 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()));
1634
1635 AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1636
1637 new TCanvas; gROOT->FindObject("esd1_yx_div_esd2_yx")->Draw("COLZ");
1638 new TCanvas; gROOT->FindObject("esd1_zx_div_esd2_zx")->Draw("COLZ");
1639 new TCanvas; gROOT->FindObject("esd1_zy_div_esd2_zy")->Draw("COLZ");
1640
1641 TH2* event1 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event1");
1642 TH2* event2 = (TH2*) fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event2");
1643
1644 Printf("event1 contains %f entries", event1->GetEntries());
1645 Printf("event2 contains %f entries", event2->GetEntries());
1646 Printf("event1 integral is %f", event1->Integral());
1647 Printf("event2 integral is %f", event2->Integral());
1648 Printf("event1 contains %f entries in |vtx-z| < 10", event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), 1, event1->GetNbinsY()));
1649 Printf("event2 contains %f entries in |vtx-z| < 10", event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), 1, event2->GetNbinsY()));
1650
1651 projx1 = event1->ProjectionX();
1652 projx2 = event2->ProjectionX();
1653
1654 new TCanvas; projx1->DrawCopy(); projx2->SetLineColor(2); projx2->DrawCopy("SAME");
1655
1656 projx1->Divide(projx2);
1657 new TCanvas; projx1->Draw();
1658
1659 event1->Divide(event2);
1660 new TCanvas; event1->Draw("COLZ");
770a1f1d 1661
1662}
3dfa46a4 1663
0fc41645 1664void DrawTrackletOrigin()
1665{
1666 TFile::Open("correction_map.root");
1667
1668 Int_t colors[] = {1,2,3,4,6,7,8,102};
1669
1670 Int_t maxHists = 8;
1671 TH1* hist[8];
1672
1673 const char* titles[] = { "PP", "SS", "PP'", "PS", "PS*", "SP", "SS'", "" };
1674
1675 TLegend* legend = new TLegend(0.75, 0.6, 0.95, 0.95);
1676
1677 Int_t total = 0;
1678 for (Int_t i=0; i<maxHists; i++)
1679 {
1680 hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
1681 //hist[i]->Rebin(20);
1682 hist[i]->SetStats(kFALSE);
1683 hist[i]->SetLineColor(colors[i]);
1684 hist[i]->GetXaxis()->SetRangeUser(-0.2, 0.2);
1685 hist[i]->Draw(((i == 0) ? "" : "SAME"));
1686
1687 total += hist[i]->GetEntries();
1688
1689 if (i != 7)
1690 legend->AddEntry(hist[i], titles[i]);
1691 }
1692
1693 legend->Draw();
1694 gPad->SetLogy();
1695
1696 Printf("Total: %d", total);
1697 for (Int_t i=0; i<maxHists; i++)
1698 Printf("Histogram %d (%s) containts %.2f %% of the entries", i, titles[i], 100.0 * hist[i]->GetEntries() / total);
1699
1700 printf("| Delta phi | Acc. %% | ");
1701 for (Int_t i=0; i<maxHists; i++)
1702 printf("%3s %% | ", titles[i]);
1703 Printf("");
1704
1705 for (Float_t f = 0.01; f < 0.09; f += 0.01)
1706 {
1707 Int_t integralBegin = hist[0]->GetXaxis()->FindBin(-f);
1708 Int_t integralEnd = hist[0]->GetXaxis()->FindBin(f);
1709
1710 Int_t total2 = 0;
1711 for (Int_t i=0; i<maxHists; i++)
1712 total2 += (Int_t) hist[i]->Integral(integralBegin, integralEnd);
1713
1714 printf("| %.2f | %6.2f | ", f, 100.0 * total2 / total);
1715
1716 for (Int_t i=0; i<maxHists; i++)
1717 printf("%6.2f | ", (hist[i]->GetEntries() > 0) ? (100.0 * hist[i]->Integral(integralBegin, integralEnd) / hist[i]->GetEntries()) : -1.0);
1718 Printf("");
1719 }
1720}