3 #if !defined(__CINT__) || defined(__MAKECINT__)
5 #include "AliPWG0Helper.h"
6 #include "dNdEtaAnalysis.h"
7 #include "AlidNdEtaCorrection.h"
21 extern TSystem* gSystem;
25 gSystem->Load("libANALYSIS");
26 gSystem->Load("libPWG0base");
29 void SetRanges(TAxis* axis)
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);
41 void SetRanges(TH1* hist)
43 SetRanges(hist->GetXaxis());
44 SetRanges(hist->GetYaxis());
45 SetRanges(hist->GetZaxis());
49 void Prepare3DPlot(TH3* hist)
51 hist->GetXaxis()->SetTitleOffset(1.5);
52 hist->GetYaxis()->SetTitleOffset(1.5);
53 hist->GetZaxis()->SetTitleOffset(1.5);
55 hist->SetStats(kFALSE);
58 void Prepare2DPlot(TH2* hist)
60 hist->SetStats(kFALSE);
61 hist->GetYaxis()->SetTitleOffset(1.4);
64 hist->SetMaximum(gMax);
69 void Prepare1DPlot(TH1* hist)
71 hist->SetLineWidth(2);
72 hist->SetStats(kFALSE);
74 hist->GetXaxis()->SetLabelOffset(0.02);
75 hist->GetXaxis()->SetTitleOffset(1.3);
76 hist->GetYaxis()->SetTitleOffset(1.3);
83 gPad->Range(0, 0, 1, 1);
84 gPad->SetLeftMargin(0.15);
85 //gPad->SetRightMargin(0.05);
86 //gPad->SetTopMargin(0.13);
87 gPad->SetBottomMargin(0.12);
95 gPad->Range(0, 0, 1, 1);
96 gPad->SetRightMargin(0.15);
97 gPad->SetLeftMargin(0.12);
103 // --- end of helpers --- begin functions ---
105 void DrawOverview(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
108 if (!TFile::Open(fileName))
111 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
112 if (!dNdEtaCorrection->LoadHistograms())
115 dNdEtaCorrection->Finish();
117 dNdEtaCorrection->DrawOverview();
120 void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
122 gSystem->Load("libPWG0base");
124 TFile::Open(esdFile);
125 dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
126 fdNdEtaAnalysisESD->LoadHistograms();
129 dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
130 fdNdEtaAnalysisMC->LoadHistograms();
131 //fdNdEtaAnalysisMC->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
133 for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
134 fdNdEtaAnalysisESD->GetdNdEtaPtCutOffCorrectedHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(i));
136 fdNdEtaAnalysisESD->DrawHistograms();
139 void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
141 gSystem->Load("libPWG0base");
143 const char* ESDfolder = 0;
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";
152 TFile::Open("analysis_esd.root");
153 dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
154 fdNdEtaAnalysisESD->LoadHistograms();
156 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
157 dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
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();
168 TH1* proj = hist->ProjectionX();
170 TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
171 for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
173 Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
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);
185 void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
187 gSystem->Load("libPWG0base");
189 TFile::Open("analysis_esd.root");
190 dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
191 fdNdEtaAnalysisESD->LoadHistograms();
193 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
194 dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
196 //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
197 //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
199 TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
200 TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
202 TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
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)
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)));
217 if (value2 > 0 && value1 > 0)
219 printf("%f %f %f\n", value1, value2, value1 / value2);
220 diff->Fill(value1 / value2);
228 Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc = "")
232 for (Int_t bin = histMC->FindBin(-0.7999); bin <= histMC->FindBin(0.7999); bin++)
234 avgMC += histMC->GetBinContent(bin);
235 avgESD += histESD->GetBinContent(bin);
237 Int_t nBins = histMC->FindBin(0.7999) - histMC->FindBin(-0.7999) + 1;
242 // deviation when integrate in |eta| < 0.8 between mc and esd
243 Double_t diffFullRange = (avgMC - avgESD) / avgMC;
245 Printf("%s: Integrated deviation in |eta| < 0.8 is %.2f %%", desc, diffFullRange * 100);
247 return diffFullRange;
250 void dNdEtaNoResolution()
254 TFile::Open("correction_map.root");
256 const char* correctionMapFolder = "dndeta_correction";
257 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
258 dNdEtaCorrection->LoadHistograms();
259 dNdEtaCorrection->GetTrack2ParticleCorrection()->PrintInfo(0.3);
261 TFile::Open("analysis_mc.root");
262 fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
263 fdNdEtaAnalysis->LoadHistograms();
264 fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, "ESD (no resolution effect) -> MB with vertex");
265 fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0)->SetMarkerStyle(21);
267 TFile::Open("analysis_mc.root");
268 dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
269 fdNdEtaAnalysisMC->LoadHistograms();
270 fdNdEtaAnalysisMC->Finish(0, 0, AlidNdEtaCorrection::kNone, "MC: MB with vertex");
272 DrawdNdEtaRatio(fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0), fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(0), "MB with vertex (no resolution effect)", 3);
275 TH1* GetMCHist(const char* folder, Float_t ptCut, const char* tag)
277 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(folder, folder);
278 fdNdEtaAnalysis->LoadHistograms();
279 fdNdEtaAnalysis->Finish(0, ptCut, AlidNdEtaCorrection::kNone, tag);
280 return fdNdEtaAnalysis->GetdNdEtaHistogram(0);
283 void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
285 TFile* file = TFile::Open("analysis_esd.root");
286 TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
287 TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
288 TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
289 TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
290 TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
291 TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
292 TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
293 TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
295 Prepare1DPlot(histESD);
296 Prepare1DPlot(histESDnsd);
297 Prepare1DPlot(histESDMB);
298 Prepare1DPlot(histESDMBVtx);
300 Prepare1DPlot(histESDNoPt);
301 Prepare1DPlot(histESDMBNoPt);
302 Prepare1DPlot(histESDMBVtxNoPt);
303 Prepare1DPlot(histESDMBTracksNoPt);
305 histESD->SetLineWidth(0);
306 histESDnsd->SetLineWidth(0);
307 histESDMB->SetLineWidth(0);
308 histESDMBVtx->SetLineWidth(0);
310 histESDNoPt->SetLineWidth(0);
311 histESDMBNoPt->SetLineWidth(0);
312 histESDMBVtxNoPt->SetLineWidth(0);
314 histESD->SetMarkerColor(1);
315 histESDnsd->SetMarkerColor(6);
316 histESDMB->SetMarkerColor(2);
317 histESDMBVtx->SetMarkerColor(3);
319 histESD->SetLineColor(1);
320 histESDnsd->SetLineColor(6);
321 histESDMB->SetLineColor(2);
322 histESDMBVtx->SetLineColor(3);
324 histESDNoPt->SetMarkerColor(1);
325 histESDMBNoPt->SetMarkerColor(2);
326 histESDMBVtxNoPt->SetMarkerColor(3);
327 histESDMBTracksNoPt->SetMarkerColor(4);
329 histESD->SetMarkerStyle(20);
330 histESDnsd->SetMarkerStyle(29);
331 histESDMB->SetMarkerStyle(21);
332 histESDMBVtx->SetMarkerStyle(22);
334 histESDNoPt->SetMarkerStyle(20);
335 histESDMBNoPt->SetMarkerStyle(21);
336 histESDMBVtxNoPt->SetMarkerStyle(22);
337 histESDMBTracksNoPt->SetMarkerStyle(23);
339 //Float_t etaLimit = 1.2999;
340 Float_t etaLimit = 2.41;
341 Float_t etaPlotLimit = 2.6;
343 histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
344 histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
345 histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
346 histESDnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
348 histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
349 histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
350 histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
351 histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
353 Float_t max = TMath::Max(histESDMBVtx->GetMaximum(), histESDMB->GetMaximum());
354 max = TMath::Max(max, histESD->GetMaximum());
356 TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 0, max * 1.1);
357 Prepare1DPlot(dummy);
358 dummy->SetStats(kFALSE);
359 dummy->SetXTitle("#eta");
360 dummy->SetYTitle("dN_{ch}/d#eta");
361 dummy->GetYaxis()->SetTitleOffset(1);
363 TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
366 histESDMBVtx->Draw("SAME");
367 histESDMB->Draw("SAME");
368 histESD->Draw("SAME");
372 canvas->SaveAs("dNdEta1.gif");
373 canvas->SaveAs("dNdEta1.eps");
381 TFile* file2 = TFile::Open("analysis_mc.root");
383 TH1* histMC = (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
384 TH1* histMCTr = (TH1*) GetMCHist("dndetaTr", -1, "MC: minimum bias")->Clone("histMCTr");
385 TH1* histMCTrVtx = (TH1*) GetMCHist("dndetaTrVtx", -1, "MC: MB with trigger")->Clone("histMCTrVtx");
386 TH1* histMCnsd = (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
388 TH1* histMCPtCut = (TH1*) GetMCHist("dndeta", 0.3, "MC: full inelastic, pt cut")->Clone("histMCPtCut");
389 TH1* histMCTrPtCut = (TH1*) GetMCHist("dndetaTr", 0.3, "MC: minimum bias, pt cut")->Clone("histMCTrPtCut");
390 TH1* histMCTrVtxPtCut = (TH1*) GetMCHist("dndetaTrVtx", 0.3, "MC: MB with trigger, pt cut")->Clone("histMCTrVtxPtCut");
391 TH1* histMCTracksPtCut = (TH1*) GetMCHist("dndetaTracks", 0.3, "MC: Tracks w/o resolution effect, pt cut")->Clone("histMCTracksPtCut");
393 Prepare1DPlot(histMC);
394 Prepare1DPlot(histMCnsd);
395 Prepare1DPlot(histMCTr);
396 Prepare1DPlot(histMCTrVtx);
398 Prepare1DPlot(histMCPtCut);
399 Prepare1DPlot(histMCTrPtCut);
400 Prepare1DPlot(histMCTrVtxPtCut);
401 Prepare1DPlot(histMCTracksPtCut);
403 histMC->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
404 histMCnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
405 histMCTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
406 histMCTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
408 histMCPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
409 histMCTrPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
410 histMCTrVtxPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
411 histMCTracksPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
413 histMC->SetLineColor(1);
414 histMCnsd->SetLineColor(6);
415 histMCTr->SetLineColor(2);
416 histMCTrVtx->SetLineColor(3);
418 histMCPtCut->SetLineColor(1);
419 histMCTrPtCut->SetLineColor(2);
420 histMCTrVtxPtCut->SetLineColor(3);
421 if (histMCTracksPtCut)
422 histMCTracksPtCut->SetLineColor(4);
424 TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
426 TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
427 dummy2->GetYaxis()->SetRangeUser(0, max * 1.1);
430 histMC->Draw("SAME");
431 histMCnsd->Draw("SAME");
432 histMCTr->Draw("SAME");
433 histMCTrVtx->Draw("SAME");
434 histESD->Draw("SAME");
435 histESDnsd->Draw("SAME");
436 histESDMB->Draw("SAME");
437 histESDMBVtx->Draw("SAME");
438 histESDNoPt->Draw("SAME");
439 histESDMBNoPt->Draw("SAME");
440 histESDMBVtxNoPt->Draw("SAME");
441 histESDMBTracksNoPt->Draw("SAME");
442 histMCPtCut->Draw("SAME");
443 histMCTrPtCut->Draw("SAME");
444 histMCTrVtxPtCut->Draw("SAME");
445 if (histMCTracksPtCut)
446 histMCTracksPtCut->Draw("SAME");
450 canvas2->SaveAs("dNdEta2.gif");
451 canvas2->SaveAs("dNdEta2.eps");
454 DrawdNdEtaRatio(histESD, histMC, "full_inelastic", etaPlotLimit);
455 DrawdNdEtaRatio(histESDMB, histMCTr, "triggered", etaPlotLimit);
456 DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit);
460 histMCnsd->Draw("SAME");
461 histESDnsd->Draw("SAME");
463 TH1* ratio = (TH1*) histMC->Clone("ratio");
464 TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
466 ratio->Divide(histESD);
467 ratioNoPt->Divide(histESDNoPt);
469 ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
471 ratio->SetLineColor(1);
472 ratioNoPt->SetLineColor(2);
474 Double_t average = 0; // average deviation from 1 in ratio (depends on the number of bins if statistical)
475 for (Int_t bin = ratio->FindBin(-0.7999); bin <= ratio->FindBin(0.7999); bin++)
476 average += TMath::Abs(ratio->GetBinContent(bin) - 1);
477 Int_t nBins = ratio->FindBin(0.7999) - ratio->FindBin(-0.7999) + 1;
479 Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
481 PrintIntegratedDeviation(histMC, histESD, "all events");
482 PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
483 PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
484 PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
485 PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
486 PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
488 TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 700, 600);
489 canvas3->Range(0, 0, 1, 1);
490 //canvas3->Divide(1, 2, 0, 0);
493 TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.5, 0.98, 0.98);
496 TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.5);
499 pad1->SetRightMargin(0.05);
500 pad2->SetRightMargin(0.05);
502 // no border between them
503 pad1->SetBottomMargin(0);
504 pad2->SetTopMargin(0);
508 TLegend* legend = new TLegend(0.4, 0.05, 0.65, 0.3);
509 legend->SetFillColor(0);
510 legend->AddEntry(histESDMBVtx, "triggered, vertex");
511 legend->AddEntry(histESDMB, "triggered");
512 legend->AddEntry(histESD, "all events");
513 legend->AddEntry(histMC, "MC prediction");
515 dummy->GetXaxis()->SetLabelSize(0.06);
516 dummy->GetYaxis()->SetLabelSize(0.06);
517 dummy->GetXaxis()->SetTitleSize(0.06);
518 dummy->GetYaxis()->SetTitleSize(0.06);
519 dummy->GetYaxis()->SetTitleOffset(0.7);
521 histESDMBVtx->Draw("SAME");
522 histESDMB->Draw("SAME");
523 histESD->Draw("SAME");
524 histMC->Draw("SAME");
529 pad2->SetBottomMargin(0.15);
531 Float_t minR = 0.9; //TMath::Min(0.961, ratio->GetMinimum() * 0.95);
532 Float_t maxR = 1.1; //TMath::Max(1.049, ratio->GetMaximum() * 1.05);
534 TH1F dummy3("dummy3", ";#eta;Ratio: MC / ESD", 100, -etaPlotLimit, etaPlotLimit);
535 dummy3.SetStats(kFALSE);
536 for (Int_t i=1; i<=100; ++i)
537 dummy3.SetBinContent(i, 1);
538 dummy3.GetYaxis()->SetRangeUser(minR, maxR);
539 dummy3.SetLineWidth(2);
540 dummy3.GetXaxis()->SetLabelSize(0.06);
541 dummy3.GetYaxis()->SetLabelSize(0.06);
542 dummy3.GetXaxis()->SetTitleSize(0.06);
543 dummy3.GetYaxis()->SetTitleSize(0.06);
544 dummy3.GetYaxis()->SetTitleOffset(0.7);
555 canvas3->SaveAs("dNdEta.gif");
556 canvas3->SaveAs("dNdEta.eps");
559 TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
562 ratioNoPt->Draw("SAME");
564 TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
565 legend->SetFillColor(0);
566 legend->AddEntry(ratio, "mc/esd");
567 legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
571 void DrawdNdEtaRatio(TH1* corr, TH1* mc, const char* name, Float_t etaPlotLimit)
573 TCanvas* canvas3 = new TCanvas(name, name, 700, 600);
574 canvas3->Range(0, 0, 1, 1);
576 TPad* pad1 = new TPad(Form("%s_1", name), "", 0, 0.5, 0.98, 0.98);
579 TPad* pad2 = new TPad(Form("%s_2", name), "", 0, 0.02, 0.98, 0.5);
582 pad1->SetRightMargin(0.05);
583 pad2->SetRightMargin(0.05);
585 // no border between them
586 pad1->SetBottomMargin(0);
587 pad2->SetTopMargin(0);
591 TLegend* legend = new TLegend(0.4, 0.05, 0.65, 0.3);
592 legend->SetFillColor(0);
593 legend->AddEntry(corr, "corrected");
594 legend->AddEntry(mc, "MC prediction");
596 TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 0, corr->GetMaximum() * 1.1);
597 Prepare1DPlot(dummy);
598 dummy->SetStats(kFALSE);
599 dummy->SetXTitle("#eta");
600 dummy->SetYTitle("dN_{ch}/d#eta");
601 dummy->GetYaxis()->SetTitleOffset(1);
603 dummy->GetXaxis()->SetLabelSize(0.06);
604 dummy->GetYaxis()->SetLabelSize(0.06);
605 dummy->GetXaxis()->SetTitleSize(0.06);
606 dummy->GetYaxis()->SetTitleSize(0.06);
607 dummy->GetYaxis()->SetTitleOffset(0.7);
616 pad2->SetBottomMargin(0.15);
618 TH1* ratio = (TH1*) mc->Clone("ratio");
621 Float_t minR = TMath::Min(0.961, ratio->GetMinimum() * 0.95);
622 Float_t maxR = TMath::Max(1.049, ratio->GetMaximum() * 1.05);
624 TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
625 dummy3.SetStats(kFALSE);
626 for (Int_t i=1; i<=100; ++i)
627 dummy3.SetBinContent(i, 1);
628 dummy3.GetYaxis()->SetRangeUser(minR, maxR);
629 dummy3.SetLineWidth(2);
630 dummy3.GetXaxis()->SetLabelSize(0.06);
631 dummy3.GetYaxis()->SetLabelSize(0.06);
632 dummy3.GetXaxis()->SetTitleSize(0.06);
633 dummy3.GetYaxis()->SetTitleSize(0.06);
634 dummy3.GetYaxis()->SetTitleOffset(0.7);
644 TFile* file = TFile::Open("analysis_esd.root");
645 TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
647 TFile* file2 = TFile::Open("analysis_mc.root");
648 TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
650 TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
654 Prepare1DPlot(histMC);
655 Prepare1DPlot(histESD);
657 histESD->SetTitle("");
658 histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
659 histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
661 histMC->SetLineColor(kBlue);
662 histESD->SetLineColor(kRed);
664 histESD->GetYaxis()->SetTitleOffset(1.5);
665 histESD->GetXaxis()->SetRangeUser(0, 4.9999);
667 histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
670 histMC->Draw("SAME");
672 canvas->SaveAs("ptSpectrum.gif");
673 canvas->SaveAs("ptSpectrum.eps");
678 gSystem->Load("libPWG0base");
680 TFile::Open("correction_map.root");
681 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
682 dNdEtaCorrection->LoadHistograms();
684 dNdEtaCorrection->GetMeasuredFraction(AlidNdEtaCorrection::kINEL, 0.3, -100, kTRUE);
686 TH1* hist = dynamic_cast<TH1*> (gROOT->FindObject("generated_pt")->Clone("ptcutoff"));
688 hist->GetXaxis()->SetRangeUser(0, 0.9999);
691 hist->SetTitle("Generated Particles");
694 TCanvas* canvas = new TCanvas("ptCutoff", "ptCutoff", 700, 500);
697 TLine* line = new TLine(0.3, 0 - hist->GetMaximum() * 0, 0.3, hist->GetMaximum() * 1.1);
698 line->SetLineWidth(3);
699 line->SetLineColor(kRed);
702 canvas->SaveAs("ptCutoff.gif");
703 canvas->SaveAs("ptCutoff.eps");
705 TH1F* factor = new TH1F("factor", ";#eta;correction factor", 20, -1, 1.000001);
706 factor->SetLineWidth(2);
707 for (Float_t eta = -0.95; eta<1; eta += 0.1)
708 factor->Fill(eta, 1.0 / dNdEtaCorrection->GetMeasuredFraction(AlidNdEtaCorrection::kINEL, 0.3, eta, kFALSE));
710 TCanvas* canvas = new TCanvas("ptCutoff_factor", "ptCutoff_factor", 700, 500);
713 Prepare1DPlot(factor);
714 factor->GetYaxis()->SetRangeUser(1, 2);
715 factor->GetYaxis()->SetTitleOffset(1);
718 canvas->SaveAs("ptCutoff_factor.eps");
721 void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
723 gSystem->Load("libPWG0base");
725 TFile::Open(fileName);
726 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
727 dNdEtaCorrection->LoadHistograms();
729 TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
730 TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
732 Prepare2DPlot(corrTrigger);
733 corrTrigger->SetTitle("b) Trigger bias correction");
735 Prepare2DPlot(corrVtx);
736 corrVtx->SetTitle("a) Vertex reconstruction correction");
738 corrTrigger->GetYaxis()->SetTitle("Multiplicity");
739 corrVtx->GetYaxis()->SetTitle("Multiplicity");
741 TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
742 canvas->Divide(2, 1);
746 corrVtx->DrawCopy("COLZ");
750 corrTrigger->DrawCopy("COLZ");
752 canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
753 canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
755 canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
756 canvas->Divide(2, 1);
758 corrTrigger->GetYaxis()->SetRangeUser(0, 5);
759 corrVtx->GetYaxis()->SetRangeUser(0, 5);
763 corrVtx->DrawCopy("COLZ");
767 corrTrigger->DrawCopy("COLZ");
769 canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
770 canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
773 void TriggerBias(const char* fileName = "correction_map.root")
775 TFile* file = TFile::Open(fileName);
777 TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
780 corr->SetTitle("Trigger bias correction");
782 TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
784 corr->DrawCopy("COLZ");
786 canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
787 canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
789 corr->GetYaxis()->SetRangeUser(0, 5);
791 canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
793 corr->DrawCopy("COLZ");
795 canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
796 canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
799 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
801 gSystem->Load("libPWG0base");
803 TFile* file = TFile::Open(fileName);
804 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
805 dNdEtaCorrection->LoadHistograms();
807 TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("x");
808 TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
810 TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
811 canvas->Divide(2, 1);
818 hist->GetYaxis()->SetTitle("correction factor");
819 hist->GetYaxis()->SetRangeUser(1, 1.5);
820 hist->GetYaxis()->SetTitleOffset(1.6);
826 Prepare1DPlot(hist2);
828 hist2->GetYaxis()->SetTitle("correction factor");
829 hist2->GetXaxis()->SetRangeUser(0, 5);
830 hist2->GetYaxis()->SetTitleOffset(1.6);
831 hist2->GetXaxis()->SetTitle("multiplicity");
834 TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
835 pave->SetFillColor(0);
836 pave->AddText("|z| < 10 cm");
839 canvas->SaveAs("TriggerBias1D.eps");
844 TFile* file = TFile::Open("correction_map.root");
846 TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
849 corr->SetTitle("Vertex reconstruction correction");
851 TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
853 corr->DrawCopy("COLZ");
855 canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
856 canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
858 corr->GetYaxis()->SetRangeUser(0, 5);
860 canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
862 corr->DrawCopy("COLZ");
864 canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
865 canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
868 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
870 gSystem->Load("libPWG0base");
872 TFile* file = TFile::Open(fileName);
873 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
874 dNdEtaCorrection->LoadHistograms();
876 TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
877 TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
879 TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
880 canvas->Divide(2, 1);
887 hist->GetYaxis()->SetTitle("correction factor");
888 hist->GetYaxis()->SetRangeUser(1, 1.8);
889 hist->GetYaxis()->SetTitleOffset(1.6);
895 Prepare1DPlot(hist2);
897 hist2->GetYaxis()->SetTitle("correction factor");
898 hist2->GetXaxis()->SetRangeUser(0, 20);
899 hist2->GetYaxis()->SetTitleOffset(1.6);
900 hist2->GetXaxis()->SetTitle("multiplicity");
903 TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
904 pave->SetFillColor(0);
905 pave->AddText("|z| < 10 cm");
908 canvas->SaveAs("VtxRecon1D.eps");
910 Correction1DCreatePlots(fileName, folderName, 9.9, 2);
912 TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
913 TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
915 Prepare1DPlot(corrX);
916 Prepare1DPlot(corrZ);
918 corrX->GetYaxis()->SetTitleOffset(1.5);
919 corrZ->GetYaxis()->SetTitleOffset(1.5);
921 corrX->SetTitle("a) z projection");
922 corrZ->SetTitle("b) p_{T} projection");
924 corrX->GetYaxis()->SetTitle("correction factor");
925 corrZ->GetYaxis()->SetTitle("correction factor");
927 corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
930 canvasName.Form("VtxRecon1D_Track");
931 TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
932 canvas->Divide(2, 1);
943 canvas->SaveAs("VtxRecon1D_Track.eps");
944 canvas->SaveAs("VtxRecon1D_Track.gif");
947 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
949 gSystem->Load("libPWG0base");
951 TFile::Open(fileName);
952 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
953 dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
955 TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
956 TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
958 TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
959 TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
961 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
962 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
963 gene->GetXaxis()->SetRangeUser(-10, 10);
964 meas->GetXaxis()->SetRangeUser(-10, 10);
966 Float_t eff1 = gene->Integral() / meas->Integral();
967 Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
969 printf("Correction without pT cut: %f +- %f\n", eff1, error1);
971 gene->GetZaxis()->SetRangeUser(0.3, 10);
972 meas->GetZaxis()->SetRangeUser(0.3, 10);
974 Float_t eff2 = gene->Integral() / meas->Integral();
975 Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
977 printf("Correction with pT cut: %f +- %f\n", eff2, error2);
979 gene->GetZaxis()->SetRangeUser(0.3, 1);
980 meas->GetZaxis()->SetRangeUser(0.3, 1);
982 Float_t eff3 = gene->Integral() / meas->Integral();
983 Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
985 printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
988 void Correction1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType = 0)
990 TFile::Open(fileName);
991 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
992 dNdEtaCorrection->LoadHistograms();
994 TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
995 TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetMeasuredHistogram();
997 gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
998 meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
999 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1000 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1001 AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kTRUE);
1002 gene->GetYaxis()->SetRange(0, 0);
1003 meas->GetYaxis()->SetRange(0, 0);
1005 gene->GetXaxis()->SetRangeUser(-10, 10);
1006 meas->GetXaxis()->SetRangeUser(-10, 10);
1007 AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kTRUE);
1008 gene->GetZaxis()->SetRange(0, 0);
1009 meas->GetZaxis()->SetRange(0, 0);
1011 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1012 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1013 AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kTRUE);
1016 void Correction1D(Int_t correctionType = 0, const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
1018 gSystem->Load("libPWG0base");
1020 Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType);
1022 TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1023 TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1024 TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1026 Prepare1DPlot(corrX);
1027 Prepare1DPlot(corrY);
1028 Prepare1DPlot(corrZ);
1030 corrX->SetTitle("a) z projection");
1031 corrY->SetTitle("b) #eta projection");
1032 corrZ->SetTitle("c) p_{T} projection");
1034 corrX->GetYaxis()->SetTitle("correction factor");
1035 corrY->GetYaxis()->SetTitle("correction factor");
1036 corrZ->GetYaxis()->SetTitle("correction factor");
1038 corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1041 canvasName.Form("Correction1D_%s", folder);
1042 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1043 canvas->Divide(3, 1);
1057 canvas->SaveAs(Form("Correction1D_%d_%s_%f.gif", correctionType, fileName, upperPtLimit));
1058 canvas->SaveAs(Form("Correction1D_%d_%s_%f.eps", correctionType, fileName, upperPtLimit));
1061 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
1063 gSystem->Load("libPWG0base");
1065 Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
1067 TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1068 TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1069 TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1071 Prepare1DPlot(corrX);
1072 Prepare1DPlot(corrY);
1073 Prepare1DPlot(corrZ);
1075 corrX->SetTitle("a) z projection");
1076 corrY->SetTitle("a) #eta projection");
1077 corrZ->SetTitle("b) p_{T} projection");
1079 corrY->GetYaxis()->SetTitle("correction factor");
1080 corrZ->GetYaxis()->SetTitle("correction factor");
1082 corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1085 canvasName.Form("Track2Particle1D_%s", folder);
1086 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1087 canvas->Divide(3, 1);
1101 canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
1102 canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
1104 //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
1106 canvasName.Form("Track2Particle1D_%s_etapt", folder);
1107 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
1108 canvas->Divide(2, 1);
1112 corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
1113 corrY->GetYaxis()->SetRangeUser(1, 1.5);
1114 corrY->GetYaxis()->SetTitleOffset(1.5);
1116 TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
1117 pave->AddText("|z| < 10 cm");
1118 pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
1124 corrZ->GetYaxis()->SetRangeUser(1, 2.5);
1125 corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
1126 corrZ->GetYaxis()->SetTitleOffset(1.5);
1128 pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
1129 pave->AddText("|z| < 10 cm");
1130 pave->AddText("|#eta| < 0.8");
1133 canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
1134 canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
1137 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
1139 gSystem->Load("libPWG0base");
1142 for (Int_t particle=0; particle<4; ++particle)
1145 dirName.Form("correction_%d", particle);
1146 Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
1148 TString tmpx, tmpy, tmpz;
1149 tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
1150 tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
1151 tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
1153 TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
1154 TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
1155 TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
1157 Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
1159 TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
1160 TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
1161 TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
1167 Prepare1DPlot(posX);
1168 Prepare1DPlot(posY);
1169 Prepare1DPlot(posZ);
1174 posX->SetMinimum(min);
1175 posX->SetMaximum(max);
1176 posY->SetMinimum(min);
1177 posY->SetMaximum(max);
1178 posZ->SetMinimum(min);
1179 posZ->SetMaximum(max);
1181 posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1183 posX->GetYaxis()->SetTitleOffset(1.7);
1184 posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
1185 posY->GetYaxis()->SetTitleOffset(1.7);
1186 posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
1187 posZ->GetYaxis()->SetTitleOffset(1.7);
1188 posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
1190 posZ->GetXaxis()->SetRangeUser(0, 1);
1193 canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
1195 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1196 canvas->Divide(3, 1);
1210 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1211 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1215 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
1217 TFile::Open(fileName);
1218 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1219 dNdEtaCorrection->LoadHistograms();
1221 TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1222 TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1224 gene->GetZaxis()->SetRangeUser(0.3, 10);
1225 meas->GetZaxis()->SetRangeUser(0.3, 10);
1226 AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
1227 gene->GetZaxis()->SetRange(0, 0);
1228 meas->GetZaxis()->SetRange(0, 0);
1230 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1231 meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1232 AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
1233 gene->GetYaxis()->SetRange(0, 0);
1234 meas->GetYaxis()->SetRange(0, 0);
1236 gene->GetXaxis()->SetRangeUser(-10, 10);
1237 meas->GetXaxis()->SetRangeUser(-10, 10);
1238 AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
1239 gene->GetXaxis()->SetRange(0, 0);
1240 meas->GetXaxis()->SetRange(0, 0);
1243 void Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1245 gSystem->Load("libPWG0base");
1247 Track2Particle2DCreatePlots(fileName);
1249 TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1250 TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1251 TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
1253 /* this reads them from the file
1254 TH2* corrYX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx"));
1255 TH2* corrZX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx"));
1256 TH2* corrZY = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy"));*/
1258 Prepare2DPlot(corrYX);
1259 Prepare2DPlot(corrZX);
1260 Prepare2DPlot(corrZY);
1262 const char* title = "";
1263 corrYX->SetTitle(title);
1264 corrZX->SetTitle(title);
1265 corrZY->SetTitle(title);
1267 TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1268 canvas->Divide(3, 1);
1272 corrYX->Draw("COLZ");
1276 corrZX->Draw("COLZ");
1280 corrZY->Draw("COLZ");
1282 canvas->SaveAs(Form("Track2Particle2D_%d.gif", gMax));
1283 canvas->SaveAs(Form("Track2Particle2D_%d.eps", gMax));
1286 void CompareTrack2Particle2D()
1288 gSystem->Load("libPWG0base");
1290 Track2Particle2DCreatePlots("correction_maponly-positive.root");
1292 TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("pos_yx"));
1293 TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("pos_zx"));
1294 TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("pos_zy"));
1296 Track2Particle2DCreatePlots("correction_maponly-negative.root");
1298 TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("neg_yx"));
1299 TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("neg_zx"));
1300 TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("neg_zy"));
1302 posYX->Divide(negYX);
1303 posZX->Divide(negZX);
1304 posZY->Divide(negZY);
1306 Prepare2DPlot(posYX);
1307 Prepare2DPlot(posZX);
1308 Prepare2DPlot(posZY);
1313 posYX->SetMinimum(min);
1314 posYX->SetMaximum(max);
1315 posZX->SetMinimum(min);
1316 posZX->SetMaximum(max);
1317 posZY->SetMinimum(min);
1318 posZY->SetMaximum(max);
1320 TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1321 canvas->Divide(3, 1);
1325 posYX->Draw("COLZ");
1329 posZX->Draw("COLZ");
1333 posZY->Draw("COLZ");
1335 canvas->SaveAs("CompareTrack2Particle2D.gif");
1336 canvas->SaveAs("CompareTrack2Particle2D.eps");
1339 void Track2Particle3D()
1341 // get left margin proper
1343 TFile* file = TFile::Open("correction_map.root");
1345 TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1347 corr->SetTitle("Correction Factor");
1348 SetRanges(corr->GetZaxis());
1350 Prepare3DPlot(corr);
1352 TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1353 canvas->SetTheta(29.428);
1354 canvas->SetPhi(16.5726);
1358 canvas->SaveAs("Track2Particle3D.gif");
1359 canvas->SaveAs("Track2Particle3D.eps");
1362 void Track2Particle3DAll()
1364 TFile* file = TFile::Open("correction_map.root");
1366 TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1367 TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1368 TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1370 gene->SetTitle("Generated Particles");
1371 meas->SetTitle("Measured Tracks");
1372 corr->SetTitle("Correction Factor");
1374 Prepare3DPlot(gene);
1375 Prepare3DPlot(meas);
1376 Prepare3DPlot(corr);
1378 TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1379 canvas->Divide(3, 1);
1391 canvas->SaveAs("Track2Particle3DAll.gif");
1392 canvas->SaveAs("Track2Particle3DAll.eps");
1395 void MultiplicityMC(Int_t xRangeMax = 50)
1397 TFile* file = TFile::Open("multiplicityMC.root");
1401 printf("multiplicityMC.root could not be opened.\n");
1405 TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1406 TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1407 TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1409 TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1410 TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1411 //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1412 for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1414 TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1415 proj->Fit("gaus", "0");
1416 correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1417 correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1421 // draw for debugging
1424 proj->GetFunction("gaus")->DrawCopy("SAME");
1427 TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1429 for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1431 Float_t mean = correction->GetBinContent(i);
1432 Float_t width = correctionWidth->GetBinContent(i);
1434 Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
1435 Int_t fillEnd = fMultiplicityESDCorrected->FindBin(mean + width * 3);
1436 printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1438 for (Int_t j=fillBegin; j <= fillEnd; ++j)
1440 fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
1444 TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
1445 fMultiplicityESDCorrectedRebinned->Rebin(10);
1446 fMultiplicityESDCorrectedRebinned->Scale(0.1);
1448 TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
1449 ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
1450 ratio->Divide(fMultiplicityMC);
1452 TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
1453 ratio2->Divide(fMultiplicityMC);
1455 TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
1456 canvas->Divide(3, 2);
1458 fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
1459 ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
1460 fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
1461 fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
1462 correction->GetXaxis()->SetRangeUser(0, xRangeMax);
1463 fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
1464 fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
1466 canvas->cd(1); //InitPad();
1467 fMultiplicityESD->Draw();
1468 fMultiplicityMC->SetLineColor(2);
1469 fMultiplicityMC->Draw("SAME");
1471 TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1472 legend->AddEntry(fMultiplicityESD, "ESD");
1473 legend->AddEntry(fMultiplicityMC, "MC");
1477 fCorrelation->Draw("COLZ");
1481 //correction->Fit("pol1");
1482 correctionWidth->SetLineColor(2);
1483 correctionWidth->Draw("SAME");
1485 legend = new TLegend(0.2, 0.7, 0.45, 0.85);
1486 legend->AddEntry(correction, "#bar{x}");
1487 legend->AddEntry(correctionWidth, "#sigma");
1493 ratio2->SetLineColor(2);
1494 ratio2->Draw("SAME");
1496 legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1497 legend->AddEntry(ratio, "uncorrected");
1498 legend->AddEntry(ratio2, "corrected");
1502 fMultiplicityESDCorrected->SetLineColor(kBlue);
1503 fMultiplicityESDCorrected->Draw();
1504 fMultiplicityMC->Draw("SAME");
1505 fMultiplicityESD->Draw("SAME");
1507 legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1508 legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1509 legend->AddEntry(fMultiplicityMC, "MC");
1510 legend->AddEntry(fMultiplicityESD, "ESD");
1514 fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
1515 fMultiplicityESDCorrectedRebinned->Draw();
1516 fMultiplicityMC->Draw("SAME");
1518 legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1519 legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
1520 legend->AddEntry(fMultiplicityMC, "MC");
1523 canvas->SaveAs("MultiplicityMC.gif");
1526 void MultiplicityESD()
1528 TFile* file = TFile::Open("multiplicityESD.root");
1532 printf("multiplicityESD.root could not be opened.\n");
1536 TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
1538 TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
1540 fMultiplicityESD->Draw();
1543 void drawPlots(Int_t max)
1562 void CompareCorrection2Measured(const char* dataInput = "analysis_esd_raw.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1566 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1567 TFile::Open(correctionMapFile);
1568 dNdEtaCorrection->LoadHistograms();
1570 TFile* file = TFile::Open(dataInput);
1574 cout << "Error. File not found" << endl;
1578 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1579 fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1583 TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
1584 hist1->SetTitle("mc");
1585 Printf("mc contains %f entries", hist1->Integral());
1586 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()));
1588 TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
1589 hist2->SetTitle("esd");
1590 Printf("esd contains %f entries", hist2->Integral());
1591 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()));
1593 AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1594 AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
1596 hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1597 hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1598 AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
1600 new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
1601 new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
1602 new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
1603 new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
1604 new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
1607 void CompareMeasured2Measured(const char* dataInput = "analysis_esd_raw.root", const char* dataInput2 = "analysis_esd_raw.root")
1611 TFile* file = TFile::Open(dataInput);
1615 cout << "Error. File not found" << endl;
1619 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1620 fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1622 TFile* file = TFile::Open(dataInput2);
1626 cout << "Error. File not found" << endl;
1630 dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta2", "dndeta2");
1631 fdNdEtaAnalysis2->LoadHistograms("fdNdEtaAnalysisESD");
1635 TH3* hist1 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd1");
1636 hist1->SetTitle("esd1");
1637 Printf("esd1 contains %f entries", hist1->GetEntries());
1638 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()));
1640 TH3* hist2 = (TH3*) fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd2");
1641 hist2->SetTitle("esd2");
1642 Printf("esd2 contains %f entries", hist2->GetEntries());
1643 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()));
1645 AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1647 new TCanvas; gROOT->FindObject("esd1_yx_div_esd2_yx")->Draw("COLZ");
1648 new TCanvas; gROOT->FindObject("esd1_zx_div_esd2_zx")->Draw("COLZ");
1649 new TCanvas; gROOT->FindObject("esd1_zy_div_esd2_zy")->Draw("COLZ");
1651 TH2* event1 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event1");
1652 TH2* event2 = (TH2*) fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event2");
1654 Printf("event1 contains %f entries", event1->GetEntries());
1655 Printf("event2 contains %f entries", event2->GetEntries());
1656 Printf("event1 integral is %f", event1->Integral());
1657 Printf("event2 integral is %f", event2->Integral());
1658 Printf("event1 contains %f entries in |vtx-z| < 10", event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), 1, event1->GetNbinsY()));
1659 Printf("event2 contains %f entries in |vtx-z| < 10", event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), 1, event2->GetNbinsY()));
1661 projx1 = event1->ProjectionX();
1662 projx2 = event2->ProjectionX();
1664 new TCanvas; projx1->DrawCopy(); projx2->SetLineColor(2); projx2->DrawCopy("SAME");
1666 projx1->Divide(projx2);
1667 new TCanvas; projx1->Draw();
1669 event1->Divide(event2);
1670 new TCanvas; event1->Draw("COLZ");
1674 void DrawTrackletOrigin()
1676 TFile::Open("correction_map.root");
1678 Int_t colors[] = {1,2,3,4,6,7,8,102};
1683 const char* titles[] = { "PP", "SS", "PP'", "PS", "PS*", "SP", "SS'", "" };
1685 TLegend* legend = new TLegend(0.75, 0.6, 0.95, 0.95);
1688 for (Int_t i=0; i<maxHists; i++)
1690 hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
1691 //hist[i]->Rebin(20);
1692 hist[i]->SetStats(kFALSE);
1693 hist[i]->SetLineColor(colors[i]);
1694 hist[i]->GetXaxis()->SetRangeUser(-0.2, 0.2);
1695 hist[i]->Draw(((i == 0) ? "" : "SAME"));
1697 total += hist[i]->GetEntries();
1700 legend->AddEntry(hist[i], titles[i]);
1706 Printf("Total: %d", total);
1707 for (Int_t i=0; i<maxHists; i++)
1708 Printf("Histogram %d (%s) containts %.2f %% of the entries", i, titles[i], 100.0 * hist[i]->GetEntries() / total);
1710 printf("| Delta phi | Acc. %% | ");
1711 for (Int_t i=0; i<maxHists; i++)
1712 printf("%3s %% | ", titles[i]);
1715 for (Float_t f = 0.01; f < 0.09; f += 0.01)
1717 Int_t integralBegin = hist[0]->GetXaxis()->FindBin(-f);
1718 Int_t integralEnd = hist[0]->GetXaxis()->FindBin(f);
1721 for (Int_t i=0; i<maxHists; i++)
1722 total2 += (Int_t) hist[i]->Integral(integralBegin, integralEnd);
1724 printf("| %.2f | %6.2f | ", f, 100.0 * total2 / total);
1726 for (Int_t i=0; i<maxHists; i++)
1727 printf("%6.2f | ", (hist[i]->GetEntries() > 0) ? (100.0 * hist[i]->Integral(integralBegin, integralEnd) / hist[i]->GetEntries()) : -1.0);
1732 TH2* GetCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
1734 // returns the correction factor with pt integrated out
1738 TFile::Open(fileName);
1740 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
1741 if (!dNdEtaCorrection->LoadHistograms())
1744 // hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram();
1746 gener = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1747 measu = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1749 gener->GetZaxis()->SetRange(gener->GetZaxis()->FindBin(ptmin), gener->GetNbinsZ()+1);
1750 TH2D *gener_xy = gener->Project3D("yx");
1752 measu->GetZaxis()->SetRange(measu->GetZaxis()->FindBin(ptmin), measu->GetNbinsZ()+1);
1753 TH2D *measu_xy = measu->Project3D("yx");
1755 cout << measu->GetZaxis()->FindBin(ptmin) << " " << measu->GetNbinsZ()+1 << endl;
1757 TCanvas *canp = new TCanvas("canp","canp",600,1000);
1758 canp->Divide(1,2,0.0001,0.0001);
1760 gener_xy->Draw("COLZ");
1762 measu_xy->Draw("COLZ");
1765 TCanvas *canpr = new TCanvas("canpr","canpr",700,500);
1767 TH2D *proj = new TH2D(*gener_xy);
1768 proj->Divide(measu_xy);
1770 // proj = hist->Project3D("yx");
1776 void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
1778 TH2* proj = GetCorrection(fileName, dirName, ptmin);
1780 const Float_t limit = 5;
1782 TString array = "{";
1783 TString arrayEnd = "}";
1785 for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
1789 for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
1791 if (begin == -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
1793 if (begin != -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
1796 Printf("Limits for y = %d are %d to %d", y, begin, end);
1800 array += Form("%d", begin);
1803 arrayEnd.Prepend(", ");
1804 arrayEnd.Prepend(Form("%d", (end == -1) ? -1 : proj->GetNbinsX() + 1 - end));
1807 arrayEnd.Prepend("{");
1809 Printf("Begin array:");
1810 Printf("%s", array.Data());
1812 Printf("End array (mirrored) (should be the same):");
1813 Printf("%s", arrayEnd.Data());
1816 void AverageMultiplicity(const char* fileName = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1820 TFile::Open(fileName);
1822 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1823 dNdEtaCorrection->LoadHistograms();
1824 TH2* events = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram();
1825 TH3* tracks = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram();
1827 Float_t nEvents = events->Integral(events->GetXaxis()->FindBin(-1), events->GetXaxis()->FindBin(1), 0, events->GetNbinsY()+1);
1828 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);
1830 Printf("%f %f --> %f", nEvents, nTracks, nTracks / nEvents);
1833 void GetAverageCorrectionFactor(Float_t etaRange = 1.5, Float_t vertexRange = 9.9, const char* rawFile = "analysis_esd_raw.root", const char* mcFile = "analysis_mc.root")
1837 TFile::Open(rawFile);
1838 dNdEtaAnalysis* raw = new dNdEtaAnalysis("dndeta", "dndeta");
1839 raw->LoadHistograms("fdNdEtaAnalysisESD");
1840 raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->GetXaxis()->SetRangeUser(-vertexRange, vertexRange);
1841 tracks = raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Project3D("y");
1842 events = raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("events", 0, raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->GetNbinsY() + 1);
1843 Float_t nEvents = events->Integral(events->FindBin(-vertexRange), events->FindBin(vertexRange));
1844 tracks->Scale(1.0 / nEvents / tracks->GetBinWidth(1));
1846 TFile::Open(mcFile);
1847 dNdEtaAnalysis* mc = new dNdEtaAnalysis("dndeta", "dndeta");
1848 mc->LoadHistograms("dndetaTrVtx");
1849 mcH = mc->GetdNdEtaPtCutOffCorrectedHistogram(0);
1852 mcH->SetLineColor(2);
1854 tracks->DrawCopy("SAME");
1857 mcH->GetYaxis()->SetRangeUser(0, 5);
1858 mcH->Divide(tracks);
1860 mcH->Fit("pol0", "", "", -etaRange, etaRange);
1863 void TrackCuts_Comparison(char* histName, Bool_t after = kFALSE, const char* fileName = "correction_map.root")
1865 // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
1866 // --> manually disable it in the run.C
1868 file = TFile::Open(fileName);
1871 Int_t colors[] = { 1, 2, 3, 4, 5, 6 };
1873 TLegend* legend = new TLegend(0.4, 0.6, 1, 1);
1874 TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
1875 TLegend* legend3 = new TLegend(0.7, 0.5, 1, 0.7);
1877 TCanvas* c1 = new TCanvas("c1", "c1", 800, 600);
1878 //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
1879 TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
1881 const char* folders2[] = { "before_cuts", "after_cuts" };
1882 for (Int_t j = 0; j < ((after) ? 2 : 1); j++)
1884 const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
1888 for (Int_t i = 0; i < 3; i++)
1891 folder.Form("%s/%s/%s", folders1[i], folders2[j], histName);
1892 TH1* hist = (TH1*) file->Get(folder);
1893 legend->AddEntry(hist, folder);
1896 hist->SetLineColor(colors[count]);
1897 hist->DrawCopy((count == 0) ? "" : "SAME");
1901 case 0: base = hist; break;
1902 case 1: prim = hist; break;
1903 case 2: sec = hist; break;
1909 TH1* eff = (TH1*) prim->Clone("eff"); eff->Reset();
1910 TH1* purity = (TH1*) prim->Clone("purity"); purity->Reset();
1912 for (Int_t bin = 1; bin <= prim->GetNbinsX(); bin++)
1914 eff->SetBinContent(bin, prim->Integral(1, bin) / prim->Integral(1, prim->GetNbinsX() + 1));
1915 purity->SetBinContent(bin, sec->Integral(1, bin) / (prim->Integral(1, bin) + sec->Integral(1, bin)));
1918 eff->GetYaxis()->SetRangeUser(0, 1);
1919 eff->SetLineColor(colors[0+j*2]);
1920 eff->SetStats(kFALSE);
1921 purity->SetLineColor(colors[1+j*2]);
1923 legend3->AddEntry(eff, Form("%s: efficiency", folders2[j]));
1924 legend3->AddEntry(purity, Form("%s: contamination", folders2[j]));
1927 eff->DrawCopy((j == 0) ? "" : "SAME");
1928 purity->DrawCopy("SAME");
1948 void TrackCuts_DCA()
1950 file = TFile::Open("correction_map.root");
1951 hist = (TH2*) file->Get("esd_track_cuts/before_cuts/dXYvsDZ");
1953 TCanvas* c1 = new TCanvas("c1", "c1", 600, 600);
1955 c1->SetRightMargin(0.12);
1956 c1->SetBottomMargin(0.12);
1958 hist->SetStats(kFALSE);
1961 ellipse = new TEllipse(0, 0, 4);
1962 ellipse->SetLineWidth(2);
1963 ellipse->SetLineStyle(2);
1964 ellipse->SetFillStyle(0);
1967 c1->SaveAs("trackcuts_dca_2d.eps");
1970 void FindNSigma(TH2* hist, Int_t nSigma = 1)
1972 TH1* proj = hist->ProjectionY();
1975 for (Int_t bin=1; bin<=proj->GetNbinsX(); bin++)
1977 if (hist->Integral(1, hist->GetNbinsX(), bin, bin) == 0)
1981 for (limit = 1; limit<=hist->GetNbinsX(); limit++)
1985 void ShowOnlyAccepted(TH2* input, const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
1987 TH2* proj = GetCorrection(fileName, dirName, ptmin);
1989 for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
1990 for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
1991 if (proj->GetBinContent(x, y) > 5 || proj->GetBinContent(x, y) == 0)
1993 proj->SetBinContent(x, y, 0);
1996 proj->SetBinContent(x, y, 1);
1999 input->Multiply(proj);