4 // plots for the note about multplicity measurements
7 #if !defined(__CINT__) || defined(__MAKECINT__)
19 #include <TStopwatch.h>
23 #include <TPaveText.h>
27 #include "AliMultiplicityCorrection.h"
28 #include "AliCorrection.h"
29 #include "AliCorrectionMatrix3D.h"
33 const char* correctionFile = "multiplicityMC.root";
34 const char* measuredFile = "multiplicityESD.root";
36 Int_t displayRange = 80; // axis range
37 Int_t ratioRange = 151; // range to calculate difference
38 Int_t longDisplayRange = 120;
40 const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
41 const char* measuredFileTPC = "multiplicityMC_TPC_0.6M.root";
42 Int_t etaRangeTPC = 1;
46 gSystem->Load("libANALYSIS");
47 gSystem->Load("libPWG0base");
52 correctionFile = correctionFileTPC;
53 measuredFile = measuredFileTPC;
54 etaRange = etaRangeTPC;
57 longDisplayRange = 100;
60 const char* GetMultLabel(Int_t etaR = -1, Bool_t trueM = kTRUE)
65 TString tmpStr((trueM) ? "True " : "Measured ");
67 tmpStr += "multiplicity";
68 //return Form("%s", tmpStr.Data());
72 tmpStr += " (full phase space)";
75 tmpStr += Form(" in |#eta| < %.1f", (etaR+1)* 0.5);
76 return Form("%s", tmpStr.Data());
79 void Smooth(TH1* hist, Int_t windowWidth = 20)
81 TH1* clone = (TH1*) hist->Clone("clone");
82 for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin)
84 Int_t min = TMath::Max(2, bin-windowWidth);
85 Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
86 Float_t average = clone->Integral(min, max) / (max - min + 1);
88 hist->SetBinContent(bin, average);
89 hist->SetBinError(bin, 0);
95 void responseMatrixPlot(const char* fileName = 0)
99 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
102 fileName = correctionFile;
104 TFile::Open(fileName);
105 mult->LoadHistograms("Multiplicity");
107 // empty under/overflow bins in x, otherwise Project3D takes them into account
108 TH1* hist = mult->GetCorrelation(etaRange);
109 for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
111 for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
113 hist->SetBinContent(0, y, z, 0);
114 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
117 hist = ((TH3*) mult->GetCorrelation(etaRange))->Project3D("zy");
118 hist->SetStats(kFALSE);
120 hist->SetTitle(Form(";%s;%s;Entries", GetMultLabel(), GetMultLabel(etaRange, kFALSE)));
121 hist->GetXaxis()->SetRangeUser(0, longDisplayRange);
122 hist->GetYaxis()->SetRangeUser(0, longDisplayRange);
124 hist->GetYaxis()->SetTitleOffset(1.3);
125 hist->GetZaxis()->SetTitleOffset(1.2);
127 TCanvas* canvas = new TCanvas("c1", "c1", 600, 600);
128 canvas->SetRightMargin(0.15);
129 canvas->SetTopMargin(0.05);
134 canvas->SaveAs("responsematrix.eps");
137 void multPythiaPhojet()
141 TFile::Open("LHC08c11_10TeV_0.5T/mb1/spd/multiplicity.root");
142 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
143 mult->LoadHistograms("Multiplicity");
144 hist1 = mult->GetMultiplicityINEL(1)->ProjectionY();
147 TFile::Open("LHC08c15_10TeV_0.5T_Phojet/mb1/spd/multiplicity.root");
148 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
149 mult->LoadHistograms("Multiplicity");
150 hist2 = mult->GetMultiplicityINEL(1)->ProjectionY();
153 legend = new TLegend(0.6, 0.7, 0.9, 0.9);
154 legend->SetFillColor(0);
155 legend->SetTextSize(0.04);
156 legend->AddEntry(hist1, "Pythia", "L");
157 legend->AddEntry(hist2, "Phojet", "L");
159 c1 = new TCanvas("c", "c", 600, 600);
160 c1->SetTopMargin(0.05);
161 c1->SetRightMargin(0.05);
162 c1->SetLeftMargin(0.12);
167 //hist1->SetMarkerStyle(20);
168 //hist2->SetMarkerStyle(24);
169 //hist2->SetMarkerColor(2);
170 hist1->SetLineWidth(2);
171 hist2->SetLineWidth(2);
172 hist2->SetLineStyle(2);
173 hist2->SetLineColor(2);
175 hist1->Scale(1.0 / hist1->Integral());
176 hist2->Scale(1.0 / hist2->Integral());
179 hist1->GetYaxis()->SetTitleOffset(1.3);
180 hist1->GetXaxis()->SetRangeUser(0, 100);
181 hist1->SetTitle(";N_{ch};P(N_{ch})");
187 c1->SaveAs("mult_pythia_phojet.eps");
190 TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
192 // normalize unfolded result to mc hist
193 result->Scale(1.0 / result->Integral(2, displayRange));
194 result->Scale(mcHist->Integral(2, displayRange));
196 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
197 canvas->Range(0, 0, 1, 1);
199 TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
202 TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
205 pad1->SetRightMargin(0.05);
206 pad2->SetRightMargin(0.05);
208 // no border between them
209 pad1->SetBottomMargin(0);
210 pad2->SetTopMargin(0);
216 mcHist->GetXaxis()->SetLabelSize(0.06);
217 mcHist->GetYaxis()->SetLabelSize(0.06);
218 mcHist->GetXaxis()->SetTitleSize(0.06);
219 mcHist->GetYaxis()->SetTitleSize(0.06);
220 mcHist->GetYaxis()->SetTitleOffset(0.6);
222 mcHist->GetXaxis()->SetRangeUser(0, displayRange);
224 mcHist->SetTitle(Form(";%s;Entries", GetMultLabel()));
225 mcHist->SetStats(kFALSE);
227 mcHist->DrawCopy("HIST E");
230 result->SetLineColor(2);
231 result->SetMarkerColor(2);
232 result->SetMarkerStyle(5);
233 result->DrawCopy("SAME PE");
235 TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
236 legend->AddEntry(mcHist, "True distribution");
237 legend->AddEntry(result, "Unfolded distribution", "P");
238 legend->SetFillColor(0);
242 pad2->SetBottomMargin(0.15);
246 TH1* ratio = (TH1*) mcHist->Clone("ratio");
248 ratio->Divide(ratio, result, 1, 1, "");
249 ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
250 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
254 // get average of ratio
256 for (Int_t i=2; i<=ratioRange; ++i)
258 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
262 printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
264 TLine* line = new TLine(0, 1, displayRange, 1);
265 line->SetLineWidth(2);
268 line = new TLine(0, 1.1, displayRange, 1.1);
269 line->SetLineWidth(2);
270 line->SetLineStyle(2);
272 line = new TLine(0, 0.9, displayRange, 0.9);
273 line->SetLineWidth(2);
274 line->SetLineStyle(2);
279 canvas->SaveAs(epsName);
284 TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
286 // draws the 3 plots in the upper plot
287 // draws the ratio between result1 and result2 in the lower plot
289 // normalize unfolded result to mc hist
290 result1->Scale(1.0 / result1->Integral(2, displayRange));
291 result1->Scale(mcHist->Integral(2, displayRange));
292 result2->Scale(1.0 / result2->Integral(2, displayRange));
293 result2->Scale(mcHist->Integral(2, displayRange));
295 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
296 canvas->Range(0, 0, 1, 1);
298 TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
301 TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
304 pad1->SetRightMargin(0.05);
305 pad2->SetRightMargin(0.05);
307 // no border between them
308 pad1->SetBottomMargin(0);
309 pad2->SetTopMargin(0);
315 mcHist->GetXaxis()->SetLabelSize(0.06);
316 mcHist->GetYaxis()->SetLabelSize(0.06);
317 mcHist->GetXaxis()->SetTitleSize(0.06);
318 mcHist->GetYaxis()->SetTitleSize(0.06);
319 mcHist->GetYaxis()->SetTitleOffset(0.6);
321 mcHist->GetXaxis()->SetRangeUser(0, displayRange);
323 mcHist->SetTitle(";true multiplicity;Entries");
324 mcHist->SetStats(kFALSE);
326 mcHist->DrawCopy("HIST E");
329 result1->SetLineColor(2);
330 result1->SetMarkerStyle(24);
331 result1->DrawCopy("SAME HISTE");
333 result2->SetLineColor(4);
334 result2->SetMarkerColor(4);
335 result2->DrawCopy("SAME HISTE");
337 TLegend* legend = new TLegend(0.5, 0.6, 0.95, 0.9);
338 legend->AddEntry(mcHist, "True distribution");
339 legend->AddEntry(result1, "Unfolded distribution (syst)", "P");
340 legend->AddEntry(result2, "Unfolded distribution (normal)", "P");
341 legend->SetFillColor(0);
342 legend->SetTextSize(0.06);
346 pad2->SetBottomMargin(0.15);
350 result1->GetXaxis()->SetLabelSize(0.06);
351 result1->GetYaxis()->SetLabelSize(0.06);
352 result1->GetXaxis()->SetTitleSize(0.06);
353 result1->GetYaxis()->SetTitleSize(0.06);
354 result1->GetYaxis()->SetTitleOffset(0.6);
356 result1->GetXaxis()->SetRangeUser(0, displayRange);
358 result1->SetTitle(Form(";%s;Entries", GetMultLabel()));
359 result1->SetStats(kFALSE);
363 TH1* ratio = (TH1*) result1->Clone("ratio");
365 ratio->Divide(ratio, result2, 1, 1, "");
366 ratio->SetLineColor(1);
367 ratio->SetMarkerColor(1);
368 ratio->SetMarkerStyle(0);
369 ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
370 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
374 // get average of ratio
376 for (Int_t i=2; i<=ratioRange; ++i)
378 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
382 printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
384 TLine* line = new TLine(0, 1, displayRange, 1);
385 line->SetLineWidth(2);
388 line = new TLine(0, 1.1, displayRange, 1.1);
389 line->SetLineWidth(2);
390 line->SetLineStyle(2);
392 line = new TLine(0, 0.9, displayRange, 0.9);
393 line->SetLineWidth(2);
394 line->SetLineStyle(2);
399 canvas->SaveAs(epsName);
404 TCanvas* Draw2ResultRatios(TH1* mcHist, TH1* result1, TH1* result2, TH1* ratio2, TH1* ratio3, TString epsName)
406 // draws the 3 plots in the upper plot
407 // draws the ratio between result1 and result2 in the lower plot
408 // also draws ratio2 and ratio3 in the lower plot, uses their name for the legend
410 // normalize unfolded result to mc hist
411 result1->Scale(1.0 / result1->Integral(2, 200));
412 result1->Scale(mcHist->Integral(2, 200));
413 result2->Scale(1.0 / result2->Integral(2, 200));
414 result2->Scale(mcHist->Integral(2, 200));
416 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
417 canvas->Range(0, 0, 1, 1);
419 TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
422 TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
425 pad1->SetRightMargin(0.05);
426 pad2->SetRightMargin(0.05);
428 // no border between them
429 pad1->SetBottomMargin(0);
430 pad2->SetTopMargin(0);
436 mcHist->GetXaxis()->SetLabelSize(0.06);
437 mcHist->GetYaxis()->SetLabelSize(0.06);
438 mcHist->GetXaxis()->SetTitleSize(0.06);
439 mcHist->GetYaxis()->SetTitleSize(0.06);
440 mcHist->GetYaxis()->SetTitleOffset(0.6);
442 mcHist->GetXaxis()->SetRangeUser(0, displayRange);
444 mcHist->SetTitle(";True multiplicity;Entries");
445 mcHist->SetStats(kFALSE);
447 mcHist->DrawCopy("HIST E");
450 result1->SetLineColor(2);
451 result1->SetMarkerColor(2);
452 result1->SetMarkerStyle(24);
453 result1->DrawCopy("SAME HISTE");
455 result2->SetLineColor(4);
456 result2->SetMarkerColor(4);
457 result2->SetMarkerStyle(5);
458 result2->DrawCopy("SAME HISTE");
460 TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
461 legend->AddEntry(mcHist, "True distribution");
462 legend->AddEntry(result1, "Unfolded distribution (5%)", "P");
463 legend->AddEntry(result2, "Unfolded distribution (normal)", "P");
464 legend->SetFillColor(0);
465 legend->SetTextSize(0.06);
469 pad2->SetBottomMargin(0.15);
473 result1->GetXaxis()->SetLabelSize(0.06);
474 result1->GetYaxis()->SetLabelSize(0.06);
475 result1->GetXaxis()->SetTitleSize(0.06);
476 result1->GetYaxis()->SetTitleSize(0.06);
477 result1->GetYaxis()->SetTitleOffset(0.6);
479 result1->GetXaxis()->SetRangeUser(0, displayRange);
481 result1->SetTitle(Form(";%s;Entries", GetMultLabel()));
482 result1->SetStats(kFALSE);
486 TH1* ratio = (TH1*) result1->Clone("ratio");
488 ratio->Divide(ratio, result2, 1, 1, "");
489 ratio->SetLineColor(1);
490 ratio->SetMarkerColor(1);
491 ratio->SetMarkerStyle(24);
492 ratio->GetYaxis()->SetTitle("Ratio (change / normal)");
493 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
495 ratio2->SetLineColor(2);
496 ratio2->SetMarkerColor(2);
497 ratio2->SetMarkerStyle(25);
499 ratio3->SetLineColor(4);
500 ratio3->SetMarkerColor(4);
501 ratio3->SetMarkerStyle(26);
504 ratio2->DrawCopy("SAME");
505 ratio3->DrawCopy("SAME");
507 legend2 = new TLegend(0.3, 0.8, 0.8, 0.95);
508 legend2->SetNColumns(3);
509 legend2->SetFillColor(0);
510 legend2->SetTextSize(0.06);
511 legend2->AddEntry(ratio, "5% change", "P");
512 legend2->AddEntry(ratio2, "2% change", "P");
513 legend2->AddEntry(ratio3, "1% change", "P");
516 // get average of ratio
518 for (Int_t i=2; i<=ratioRange; ++i)
520 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
524 printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
526 TLine* line = new TLine(0, 1, displayRange, 1);
527 line->SetLineWidth(2);
530 line = new TLine(0, 1.1, displayRange, 1.1);
531 line->SetLineWidth(2);
532 line->SetLineStyle(2);
534 line = new TLine(0, 0.9, displayRange, 0.9);
535 line->SetLineWidth(2);
536 line->SetLineStyle(2);
541 canvas->SaveAs(epsName);
546 TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE)
548 // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
549 // a systematic effect
552 result->Scale(1.0 / result->Integral(2, 200));
554 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 500);
555 canvas->SetTopMargin(0.05);
556 canvas->SetRightMargin(0.05);
560 result->GetXaxis()->SetRangeUser(0, displayRange);
561 result->GetYaxis()->SetRangeUser(0.55, 1.45);
562 result->SetStats(kFALSE);
564 // to get the axis how we want it
565 TH1* dummy = (TH1*) result->Clone("dummy");
567 dummy->SetTitle(Form(";%s;Ratio", GetMultLabel()));
571 Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
573 TLegend* legend = new TLegend(0.2, 0.7, 0.7, 0.93);
574 legend->SetFillColor(0);
575 legend->SetTextSize(0.04);
577 legend->SetNColumns(2);
579 for (Int_t n=0; n<nResultSyst; ++n)
581 resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
584 TH1* ratio = (TH1*) result->Clone("ratio");
585 ratio->Divide(ratio, resultSyst[n], 1, 1, "");
586 ratio->GetXaxis()->SetRangeUser(0, displayRange);
589 ratio->SetMarkerStyle(5);
591 ratio->SetLineColor(colors[n / 2]);
593 ratio->SetLineStyle(2);
595 TString drawStr("SAME HIST");
596 if (n == 0 && firstMarker)
601 ratio->DrawCopy(drawStr);
603 if (legendStrings && legendStrings[n])
604 legend->AddEntry(ratio, legendStrings[n], "L");
606 // get average of ratio
608 for (Int_t i=2; i<=ratioRange; ++i)
609 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
612 printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
618 TLine* line = new TLine(-0.5, 1, displayRange, 1);
619 line->SetLineWidth(2);
622 line = new TLine(-0.5, 1.1, displayRange, 1.1);
623 line->SetLineWidth(2);
624 line->SetLineStyle(2);
626 line = new TLine(-0.5, 0.9, displayRange, 0.9);
627 line->SetLineWidth(2);
628 line->SetLineStyle(2);
631 canvas->SaveAs(epsName);
632 canvas->SaveAs(Form("%s.gif", epsName.Data()));
637 TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE)
639 // draws the ratios of each mc to the corresponding result
641 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
642 canvas->SetRightMargin(0.05);
643 canvas->SetTopMargin(0.05);
645 for (Int_t n=0; n<nResultSyst; ++n)
648 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
649 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
651 result[n]->GetXaxis()->SetRangeUser(0, displayRange);
652 result[n]->SetStats(kFALSE);
655 TH1* ratio = (TH1*) result[n]->Clone("ratio");
656 ratio->Divide(mc[n], ratio, 1, 1, "B");
658 // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis...
659 ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0);
664 ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : "")));
665 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
669 ratio->SetLineColor((n/2)+1);
670 ratio->SetLineStyle((n%2)+1);
673 ratio->SetLineColor(n+1);
675 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
677 // get average of ratio
679 for (Int_t i=2; i<=ratioRange; ++i)
680 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
683 printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
686 TLine* line = new TLine(0, 1, displayRange, 1);
687 line->SetLineWidth(2);
690 line = new TLine(0, 1.1, displayRange, 1.1);
691 line->SetLineWidth(2);
692 line->SetLineStyle(2);
694 line = new TLine(0, 0.9, displayRange, 0.9);
695 line->SetLineWidth(2);
696 line->SetLineStyle(2);
701 canvas->SaveAs(epsName);
702 canvas->SaveAs(Form("%s.gif", epsName.Data()));
707 TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
709 // draws the ratios of each mc to the corresponding result
710 // deducts from each ratio the ratio of mcBase / resultBase
713 resultBase->Scale(1.0 / resultBase->Integral(2, 200));
714 mcBase->Scale(1.0 / mcBase->Integral(2, 200));
717 TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
718 ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
720 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
721 canvas->SetRightMargin(0.05);
722 canvas->SetTopMargin(0.05);
724 for (Int_t n=0; n<nResultSyst; ++n)
727 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
728 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
730 result[n]->GetXaxis()->SetRangeUser(0, displayRange);
731 result[n]->SetStats(kFALSE);
734 TH1* ratio = (TH1*) result[n]->Clone("ratio");
735 ratio->Divide(mc[n], ratio, 1, 1, "B");
736 ratio->Add(ratioBase, -1);
738 ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
739 ratio->GetYaxis()->SetRangeUser(-1, 1);
740 ratio->SetLineColor(n+1);
741 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
743 // get average of ratio
745 for (Int_t i=2; i<=ratioRange; ++i)
746 sum += TMath::Abs(ratio->GetBinContent(i));
749 printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, ratioRange, sum);
752 TLine* line = new TLine(0, 0, displayRange, 0);
753 line->SetLineWidth(2);
756 line = new TLine(0, 0.1, displayRange, 0.1);
757 line->SetLineWidth(2);
758 line->SetLineStyle(2);
760 line = new TLine(0, -0.1, displayRange, -0.1);
761 line->SetLineWidth(2);
762 line->SetLineStyle(2);
767 canvas->SaveAs(epsName);
768 canvas->SaveAs(Form("%s.gif", epsName.Data()));
773 TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
775 // draws the ratios of each mc to the corresponding result
776 // deducts from each ratio the ratio of mcBase / resultBase
777 // smoothens the ratios by a sliding window
780 resultBase->Scale(1.0 / resultBase->Integral(2, 200));
781 mcBase->Scale(1.0 / mcBase->Integral(2, 200));
784 TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
785 ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
787 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
788 canvas->SetRightMargin(0.05);
789 canvas->SetTopMargin(0.05);
791 for (Int_t n=0; n<nResultSyst; ++n)
794 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
795 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
797 result[n]->GetXaxis()->SetRangeUser(0, displayRange);
798 result[n]->SetStats(kFALSE);
801 TH1* ratio = (TH1*) result[n]->Clone("ratio");
802 ratio->Divide(mc[n], ratio, 1, 1, "B");
803 ratio->Add(ratioBase, -1);
805 //new TCanvas; ratio->DrawCopy();
807 ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0);
811 //ratio->SetLineColor(1); ratio->DrawCopy("SAME");
814 ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)");
815 ratio->GetYaxis()->SetRangeUser(-0.3, 0.3);
816 ratio->SetLineColor((n / 2)+1);
817 ratio->SetLineStyle((n % 2)+1);
818 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
820 // get average of ratio
822 for (Int_t i=2; i<=150; ++i)
823 sum += TMath::Abs(ratio->GetBinContent(i));
826 printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
829 TLine* line = new TLine(0, 0, displayRange, 0);
830 line->SetLineWidth(2);
833 line = new TLine(0, 0.1, displayRange, 0.1);
834 line->SetLineWidth(2);
835 line->SetLineStyle(2);
837 line = new TLine(0, -0.1, displayRange, -0.1);
838 line->SetLineWidth(2);
839 line->SetLineStyle(2);
844 canvas->SaveAs(epsName);
845 canvas->SaveAs(Form("%s.gif", epsName.Data()));
850 void DrawResiduals(const char* fileName, const char* epsName)
854 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
855 TFile::Open(fileName);
856 mult->LoadHistograms("Multiplicity");
858 TH1* measured = mult->GetMultiplicityESD(etaRange)->ProjectionY("myesd", 1, 1);
859 TH1* unfoldedFolded = mult->CalculateMultiplicityESD(mult->GetMultiplicityESDCorrected(etaRange), etaRange)->ProjectionY("myfolded", 1, 1);
862 unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, displayRange+1));
863 unfoldedFolded->Scale(measured->Integral(2, displayRange+1));
865 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
866 canvas->Range(0, 0, 1, 1);
868 TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 1, 1);
873 TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 1, 0.5);
878 TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
881 pad3->SetRightMargin(0.05);
882 pad3->SetTopMargin(0.05);
885 pad1->SetRightMargin(0.05);
886 pad2->SetRightMargin(0.05);
888 // no border between them
889 pad1->SetBottomMargin(0);
890 pad2->SetTopMargin(0);
894 measured->GetXaxis()->SetLabelSize(0.06);
895 measured->GetYaxis()->SetLabelSize(0.06);
896 measured->GetXaxis()->SetTitleSize(0.06);
897 measured->GetYaxis()->SetTitleSize(0.06);
898 measured->GetYaxis()->SetTitleOffset(0.6);
900 measured->GetXaxis()->SetRangeUser(0, displayRange);
902 measured->SetTitle(Form(";%s;Entries", GetMultLabel(etaRange, kFALSE)));
903 measured->SetStats(kFALSE);
905 measured->DrawCopy("HIST");
908 unfoldedFolded->SetMarkerStyle(5);
909 unfoldedFolded->SetMarkerColor(2);
910 unfoldedFolded->SetLineColor(2);
911 unfoldedFolded->DrawCopy("SAME PHIST");
913 TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
914 legend->AddEntry(measured, "Measured distribution", "L");
915 legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution", "P");
916 legend->SetFillColor(0);
917 legend->SetTextSize(0.06);
921 pad2->SetBottomMargin(0.15);
925 TH1* residual = (TH1*) measured->Clone("residual");
926 unfoldedFolded->Sumw2();
928 residual->Add(unfoldedFolded, -1);
931 TH1* residualHist = new TH1F("residualHist", ";", 11, -3, 3);
932 residualHist->Sumw2();
935 for (Int_t i=1; i<=displayRange+1; ++i)
937 if (measured->GetBinError(i) > 0)
939 residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
940 residual->SetBinError(i, 1);
942 residualHist->Fill(residual->GetBinContent(i));
943 chi2 += residual->GetBinContent(i) * residual->GetBinContent(i);
947 residual->SetBinContent(i, 0);
948 residual->SetBinError(i, 0);
952 Printf("chi2 / ndf = %f / %d = %f", chi2, displayRange+1, chi2 / (displayRange+1));
954 residual->GetYaxis()->SetTitle("Residuals: (1/e) (M - R #otimes U)");
955 residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
956 residual->DrawCopy();
958 TLine* line = new TLine(-0.5, 0, displayRange + 0.5, 0);
959 line->SetLineWidth(2);
963 residualHist->SetStats(kFALSE);
964 residualHist->SetLabelSize(0.08, "xy");
965 residualHist->Fit("gaus");
966 residualHist->Draw("HIST");
967 residualHist->FindObject("gaus")->Draw("SAME");
970 canvas->SaveAs(canvas->GetName());
972 //const char* epsName2 = "proj.eps";
973 //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
974 //canvas->SetGridx();
975 //canvas->SetGridy();
977 //canvas->SaveAs(canvas->GetName());
980 void chi2FluctuationResult()
984 TFile::Open("chi2_noregularization.root");
985 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
986 mult->LoadHistograms("Multiplicity");
988 TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
990 mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
992 TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_1");
993 ((TH1*) canvas->FindObject("proj"))->GetXaxis()->SetRangeUser(0, displayRange);
994 ((TH1*) canvas->FindObject("fCurrentESD"))->GetXaxis()->SetRangeUser(0, displayRange);
995 //((TH1*) canvas->FindObject("proj"))->GetXaxis()->SetTitle(GetMultTitle());
996 //((TH1*) canvas->FindObject("fCurrentESD"))->GetXaxis()->->SetTitle(GetMultTitle(etaRange, kFALSE));
997 canvas->SaveAs("chi2FluctuationResult.eps");
1000 void DrawUnfolded(const char* fileName, const char* eps)
1004 TFile::Open(fileName);
1006 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1007 mult->LoadHistograms("Multiplicity");
1009 TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, mult->GetMultiplicityVtx(etaRange)->GetNbinsX());
1010 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1012 DrawResultRatio(mcHist, result, eps);
1015 void minimizationInfluenceAlpha()
1019 TFile::Open(measuredFile);
1020 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1021 mult2->LoadHistograms("Multiplicity");
1023 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, mult2->GetMultiplicityVtx(etaRange)->GetNbinsX());
1024 mcHist->Scale(1.0 / mcHist->Integral());
1025 mcHist->SetStats(kFALSE);
1026 mcHist->SetTitle(";True multiplicity n in |#eta| < 1.0;P(n)");
1028 TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 350);
1029 canvas->Divide(3, 1, 0.005);
1031 TFile::Open("chi2compare-influencealpha/EvaluateChi2Method.root");
1033 TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_3.162278");
1034 TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_07_2_10000.000000");
1035 TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_13_2_10000000.000000");
1037 /*mcHist->Rebin(2); mcHist->Scale(0.5);
1038 hist1->Rebin(2); hist1->Scale(0.5);
1039 hist2->Rebin(2); hist2->Scale(0.5);
1040 hist3->Rebin(2); hist3->Scale(0.5);*/
1042 mcHist->GetXaxis()->SetRangeUser(0, displayRange);
1043 mcHist->SetLabelSize(0.06, "xy");
1044 mcHist->SetTitleSize(0.06, "xy");
1045 mcHist->GetYaxis()->SetTitleOffset(1.5);
1050 gPad->SetRightMargin(0.03);
1051 gPad->SetLeftMargin(0.19);
1052 gPad->SetTopMargin(0.05);
1053 gPad->SetBottomMargin(0.13);
1057 hist1->SetMarkerStyle(5);
1058 hist1->SetMarkerColor(2);
1059 hist1->Draw("SAME PE");
1062 gPad->SetRightMargin(0.03);
1063 gPad->SetLeftMargin(0.19);
1064 gPad->SetTopMargin(0.05);
1065 gPad->SetBottomMargin(0.13);
1070 hist2->SetMarkerStyle(5);
1071 hist2->SetMarkerColor(2);
1072 hist2->Draw("SAME PE");
1075 gPad->SetRightMargin(0.03);
1076 gPad->SetLeftMargin(0.19);
1077 gPad->SetTopMargin(0.05);
1078 gPad->SetBottomMargin(0.13);
1083 hist3->SetMarkerStyle(5);
1084 hist3->SetMarkerColor(2);
1085 hist3->Draw("SAME PE");
1087 canvas->SaveAs("minimizationInfluenceAlpha.eps");
1092 gSystem->Load("libPWG0base");
1094 TFile::Open(correctionFile);
1095 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1096 mult->LoadHistograms("Multiplicity");
1098 TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
1099 fCurrentESD->Sumw2();
1100 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1102 TF1* func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])");
1103 func->SetParNames("scaling", "averagen", "k");
1104 func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
1105 func->SetParLimits(1, 0.001, 1000);
1106 func->SetParLimits(2, 0.001, 1000);
1107 func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
1109 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
1110 lognormal->SetParNames("scaling", "mean", "sigma");
1111 lognormal->SetParameters(1, 1, 1);
1112 lognormal->SetParLimits(0, 0, 10);
1113 lognormal->SetParLimits(1, 0, 100);
1114 lognormal->SetParLimits(2, 1e-3, 10);
1116 TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
1117 fCurrentESD->SetStats(kFALSE);
1118 fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
1119 fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
1120 fCurrentESD->Draw("HIST");
1121 fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1122 fCurrentESD->Fit(func, "W0", "", 0, 50);
1123 func->SetRange(0, 100);
1125 printf("chi2 = %f\n", func->GetChisquare());
1127 fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
1128 lognormal->SetLineColor(2);
1129 lognormal->SetLineStyle(2);
1130 lognormal->SetRange(0, 100);
1131 lognormal->Draw("SAME");
1133 canvas->SaveAs("NBDFit.eps");
1136 void StartingConditions()
1138 // data generated by runMultiplicitySelector.C StartingConditions
1140 const char* name = "StartingConditions";
1142 TFile* file = TFile::Open(Form("%s.root", name));
1144 TCanvas* canvas = new TCanvas(name, name, 1200, 600);
1145 canvas->Divide(2, 1);
1147 TH1* mc = (TH1*) file->Get("mc");
1149 mc->Scale(1.0 / mc->Integral(2, displayRange));
1151 //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
1153 TLegend* legend = new TLegend(0.6, 0.7, 0.99, 0.99);
1154 legend->SetFillColor(0);
1155 legend->SetTextSize(0.04);
1157 const char* names[] = { "True", "Measured 1", "Measured 2", "Measured 3", "NBD", "Flat" };
1159 for (Int_t i=0; i<6; ++i)
1165 TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
1166 TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
1168 chi2Result->Scale(1.0 / chi2Result->Integral(2, displayRange));
1169 bayesResult->Scale(1.0 / bayesResult->Integral(2, displayRange));
1171 chi2Result->Divide(mc, chi2Result, 1, 1, "");
1172 bayesResult->Divide(mc, bayesResult, 1, 1, "");
1174 chi2Result->SetTitle(Form("a) #chi^{2}-minimization;%s;MC / unfolded", GetMultLabel()));
1175 chi2Result->GetXaxis()->SetRangeUser(0, displayRange);
1176 chi2Result->GetYaxis()->SetRangeUser(0.7, 1.3);
1177 chi2Result->GetYaxis()->SetTitleOffset(1.7);
1178 //chi2Result->SetMarkerStyle(marker[i]);
1179 chi2Result->SetLineColor(i+1);
1180 chi2Result->SetMarkerColor(i+1);
1181 chi2Result->SetStats(kFALSE);
1183 bayesResult->SetTitle(Form("b) Bayesian unfolding;%s;MC / unfolded", GetMultLabel()));
1184 bayesResult->GetXaxis()->SetRangeUser(0, displayRange);
1185 bayesResult->GetYaxis()->SetRangeUser(0.7, 1.3);
1186 bayesResult->GetYaxis()->SetTitleOffset(1.7);
1187 bayesResult->SetStats(kFALSE);
1188 //bayesResult->SetLineColor(2);
1189 bayesResult->SetLineColor(i+1);
1192 gPad->SetRightMargin(0.05);
1193 gPad->SetLeftMargin(0.12);
1196 chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1199 gPad->SetRightMargin(0.05);
1200 gPad->SetLeftMargin(0.12);
1203 bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1205 //TLine* line = new TLine(0, 1, 150, 1);
1208 legend->AddEntry(chi2Result, names[i]);
1215 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1218 void StatisticsPlot()
1220 const char* name = "StatisticsPlot";
1222 TFile* file = TFile::Open(Form("%s.root", name));
1224 TCanvas* canvas = new TCanvas(name, name, 600, 400);
1226 TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
1227 fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1228 fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1229 fitResultsChi2->Draw("AP");
1231 TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1232 fitResultsChi2->Fit(f);
1234 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1239 const char* plotname = "chi2Result";
1241 name = "StatisticsPlotRatios";
1242 canvas = new TCanvas(name, name, 600, 400);
1244 for (Int_t i=0; i<5; ++i)
1246 mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1247 result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1249 result[i]->SetLineColor(i+1);
1250 result[i]->Draw(((i == 0) ? "" : "SAME"));
1254 void Draw2Unfolded(const char* file1, const char* file2, const char* output)
1259 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1260 mult->LoadHistograms("Multiplicity");
1262 // result with systematic effect
1264 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1265 mult2->LoadHistograms("Multiplicity");
1267 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1268 TH1* result1 = (TH1*) mult2->GetMultiplicityESDCorrected(etaRange)->Clone("result1"); // from file2 (with syst)
1269 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2"); // from file1 (without syst)
1271 DrawResultRatio(mcHist, result1, "tmp1.eps");
1272 DrawResultRatio(mcHist, result2, "tmp2.eps");
1273 Draw2ResultRatio(mcHist, result1, result2, output);
1281 Draw2Unfolded("self.root", "pythia.root", "test.eps");
1283 canvas = (TCanvas*) gROOT->GetListOfCanvases()->Last();
1284 pad1 = (TCanvas*)canvas->GetListOfPrimitives()->First();
1285 pad2 = (TCanvas*)canvas->GetListOfPrimitives()->Last();
1286 legend = (TLegend*)pad1->GetListOfPrimitives()->Last();
1288 ((TH1*)pad2->GetListOfPrimitives()->At(1))->GetYaxis()->SetTitle("Ratio (Pythia / Phojet)");
1289 ((TLegendEntry*)legend->GetListOfPrimitives()->At(1))->SetLabel("Unfolded distribution (Pythia)");
1290 ((TLegendEntry*)legend->GetListOfPrimitives()->At(2))->SetLabel("Unfolded distribution (Phojet)");
1291 canvas->SaveAs("PythiaPhojet.eps");
1298 Draw2Unfolded("chi2_ideal.root", "chi2_misaligned.root", "test.eps");
1300 canvas = (TCanvas*) gROOT->GetListOfCanvases()->Last();
1301 pad1 = (TCanvas*)canvas->GetListOfPrimitives()->First();
1302 pad2 = (TCanvas*)canvas->GetListOfPrimitives()->Last();
1303 legend = (TLegend*)pad1->GetListOfPrimitives()->Last();
1305 ((TH1*)pad2->GetListOfPrimitives()->At(1))->GetYaxis()->SetTitle("Ratio (misaligned / realigned)");
1306 ((TLegendEntry*)legend->GetListOfPrimitives()->At(1))->SetLabel("Unfolded distribution (misaligned)");
1307 ((TLegendEntry*)legend->GetListOfPrimitives()->At(2))->SetLabel("Unfolded distribution (realigned)");
1308 canvas->SaveAs("SystematicMisalignment.eps");
1311 void SystematicLowEfficiency()
1313 Draw2Unfolded("chi2.root", "chi2_loweff_5.root", "SystematicLowEfficiency.eps");
1316 void SystematicLowEfficiency2()
1320 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("chi2.root");
1321 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
1322 TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1323 result2->Scale(1.0 / result2->Integral(2, displayRange));
1324 result2->Scale(mcHist->Integral(2, displayRange));
1326 mult = AliMultiplicityCorrection::Open("chi2_loweff_5.root");
1327 TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
1329 mult = AliMultiplicityCorrection::Open("chi2_loweff_2.root");
1330 TH1* ratio2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("ratio2");
1331 ratio2->Scale(1.0 / ratio2->Integral(2, displayRange));
1332 ratio2->Scale(mcHist->Integral(2, displayRange));
1333 ratio2->Divide(result2);
1335 mult = AliMultiplicityCorrection::Open("chi2_loweff_1.root");
1336 TH1* ratio3 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("ratio3");
1337 ratio3->Scale(1.0 / ratio3->Integral(2, displayRange));
1338 ratio3->Scale(mcHist->Integral(2, displayRange));
1339 ratio3->Divide(result2);
1341 Draw2ResultRatios(mcHist, result1, result2, ratio2, ratio3, "SystematicLowEfficiency2.eps");
1344 void SystematicMisalignment()
1346 gSystem->Load("libPWG0base");
1348 TFile::Open(correctionFile);
1349 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1350 mult->LoadHistograms("Multiplicity");
1352 TFile::Open("multiplicityMC_100k_fullmis.root");
1353 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1354 mult2->LoadHistograms("Multiplicity");
1356 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1357 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1358 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1360 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1361 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1363 DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1366 void SystematicMisalignmentTPC()
1368 gSystem->Load("libPWG0base");
1372 TFile::Open(correctionFile);
1373 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1374 mult->LoadHistograms("Multiplicity");
1376 TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
1377 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1378 mult2->LoadHistograms("Multiplicity");
1380 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1381 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1382 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1384 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1385 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1387 DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
1390 void LowMomentumEffectSPD()
1392 // this function increases/reduces the correction as function of pt between 0 and 0.2 gev by +-50% to 0% and checks the effect on the overall correction factor
1393 // only a normal acceptance region is considered to not get a bias by the edges
1396 TFile::Open("multiplicity.root");
1399 AliCorrection* correction[8];
1402 for (Int_t loop=0; loop<3; loop++)
1405 Float_t sumMeas = 0;
1407 Printf("loop %d", loop);
1408 for (Int_t i=0; i<4; ++i)
1410 Printf("correction %d", i);
1412 TString name; name.Form("correction_%d", i);
1413 correction[i] = new AliCorrection(name, name);
1414 correction[i]->LoadHistograms();
1416 TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1417 TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1419 Float_t vtxRange = 5.9;
1420 gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1421 meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1423 Float_t etaRange = 0.99;
1424 gene->GetYaxis()->SetRangeUser(-etaRange, etaRange);
1425 meas->GetYaxis()->SetRangeUser(-etaRange, etaRange);
1427 TH1* genePt = gene->Project3D(Form("z_%d", i));
1428 TH1* measPt = meas->Project3D(Form("z_%d", i));
1432 for (Int_t x=1; x<=genePt->GetNbinsX(); x++)
1434 Float_t pt = genePt->GetXaxis()->GetBinCenter(x);
1440 factor = 1.5 - pt / 0.2 * 0.5;
1442 factor = 0.5 + pt / 0.2 * 0.5;
1443 //Printf("%f", factor);
1444 genePt->SetBinContent(x, genePt->GetBinContent(x) * factor);
1445 measPt->SetBinContent(x, measPt->GetBinContent(x) * factor);
1450 //new TCanvas; genePt->DrawCopy(); measPt->DrawCopy("SAME");
1452 sumGen += genePt->Integral();
1453 sumMeas += measPt->Integral();
1455 Float_t average = measPt->Integral() / genePt->Integral();
1457 Printf("The average efficiency of this correction is %f", average);
1460 Float_t average = sumMeas / sumGen;
1462 Printf("The average efficiency of all corrections is %f", average);
1463 values[loop] = average;
1466 Printf("relative is %f and %f", values[1] / values[0], values[2] / values[0]);
1470 void EfficiencySpecies(Bool_t addDecayStopped = kFALSE)
1474 Int_t marker[] = {24, 25, 26, 27};
1475 Int_t color[] = {1, 2, 4, 3};
1478 //const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
1479 //const char* fileName[] = { "spd/multiplicity.root", "tpc/multiplicity.root" };
1480 const char* fileName[] = { "multiplicity.root", "multiplicity.root" };
1481 Float_t etaRangeArr[] = {0.49, 0.9};
1482 const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1484 TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1485 canvas->Divide(2, 1);
1487 TCanvas* canvas3 = new TCanvas("EfficiencySpecies_comb", "EfficiencySpecies_comb", 600, 600);
1490 gPad->SetRightMargin(0.05);
1491 gPad->SetTopMargin(0.05);
1493 TLegend* legends[2];
1495 for (Int_t loop=1; loop<2; ++loop)
1497 Printf("%s", fileName[loop]);
1499 TCanvas* canvas2 = new TCanvas(Form("EfficiencySpecies_%d", loop), Form("EfficiencySpecies_%d", loop), 600, 600);
1502 gPad->SetRightMargin(0.05);
1503 gPad->SetTopMargin(0.05);
1505 AliCorrection* correction[8];
1511 gPad->SetRightMargin(0.05);
1513 TLegend* legend = new TLegend(0.6, 0.4, 0.85, 0.6);
1514 legend->SetFillColor(0);
1515 legend->SetEntrySeparation(0.2);
1516 legend->SetTextSize(gStyle->GetTextSize());
1518 legends[loop] = new TLegend(0.4+loop*0.3, 0.2, 0.6+loop*0.3, 0.5);
1519 legends[loop]->SetFillColor(0);
1520 legends[loop]->SetEntrySeparation(0.2);
1521 legends[loop]->SetTextSize(gStyle->GetTextSize());
1522 legends[loop]->SetHeader((loop == 0) ? "SPD" : "TPC");
1527 TFile* file = TFile::Open(fileName[loop]);
1530 Printf("Could not open %s", fileName[loop]);
1535 Float_t sumMeas = 0;
1537 for (Int_t i=0; i<3; ++i)
1539 Printf("correction %d", i);
1541 TString name; name.Form("correction_%d", i);
1542 correction[i] = new AliCorrection(name, name);
1543 correction[i]->LoadHistograms();
1545 TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1546 TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1549 Float_t vtxRange = 3.9;
1550 gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1551 meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1553 // empty over/underflow bin in eta, setting range to +-2 is not enough because this is the maximum range, Project3D takes them into account then (might be a bug)
1554 /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1555 for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1557 gene->SetBinContent(x, 0, z, 0);
1558 gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1559 meas->SetBinContent(x, 0, z, 0);
1560 meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1564 Float_t etaBegin = -etaRangeArr[loop];
1565 Float_t etaEnd = etaRangeArr[loop];
1568 gene->GetYaxis()->SetRangeUser(etaBegin, etaEnd);
1569 meas->GetYaxis()->SetRangeUser(etaBegin, etaEnd);
1571 TH1* genePt = gene->Project3D(Form("z_%d", i));
1572 TH1* measPt = meas->Project3D(Form("z_%d", i));
1577 for (Int_t x=0; x<=genePt->GetNbinsX()+1; x++)
1579 genePt->SetBinError(x, TMath::Sqrt(genePt->GetBinContent(x)));
1580 measPt->SetBinError(x, TMath::Sqrt(measPt->GetBinContent(x)));
1583 sumGen += genePt->Integral();
1584 sumMeas += measPt->Integral();
1586 TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1588 effPt->Divide(measPt, genePt, 1, 1, "B");
1591 for (bin=20; bin>=1; bin--)
1593 if (effPt->GetBinContent(bin) < 0.5)
1597 Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1599 Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1600 Printf("%.4f of the particles are below that momentum", fraction);
1602 below += genePt->Integral(1, bin);
1603 total += genePt->Integral();
1605 effPt->SetLineColor(color[i]);
1606 effPt->SetMarkerColor(color[i]);
1607 effPt->SetMarkerStyle(marker[i]);
1608 effPt->SetMarkerSize(2);
1610 effPt->GetXaxis()->SetRangeUser(0, 1);
1611 effPt->GetYaxis()->SetRangeUser(0.001, 1);
1613 effPt->GetXaxis()->SetTitleOffset(1.1);
1614 effPt->GetYaxis()->SetTitleOffset(1.2);
1616 effPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1618 effPt->SetStats(kFALSE);
1619 effPt->SetTitle(titles[loop]);
1620 effPt->GetYaxis()->SetTitle("Efficiency");
1623 effPt->DrawCopy((i == 0) ? "" : "SAME");
1626 effPt->SetTitle("");
1627 effPt->DrawCopy((i == 0) ? "" : "SAME");
1630 effPtClone = (TH1*) effPt->Clone("effPtClone");
1631 effPtClone->SetMarkerStyle(marker[i]-4*loop);
1632 effPtClone->DrawCopy((i == 0 && loop == 0) ? "" : "SAME");
1634 legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P");
1635 legends[loop]->AddEntry(effPtClone, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P");
1636 //legend2->AddEntry(effPt, Form("%s %s", (loop == 0) ? "SPD" : "TPC", ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}"))), "P");
1638 if (addDecayStopped)
1640 name.Form("correction_%d", i+4);
1641 corr = new AliCorrection(name, name);
1642 corr->LoadHistograms();
1644 TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
1645 TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
1648 gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1649 meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1650 gene->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]);
1651 meas->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]);
1653 TH1* decayed = gene->Project3D(Form("z_%d", i+4));
1654 TH1* stopped = meas->Project3D(Form("z_%d", i+4));
1656 Printf("%d: %d decayed, %d stopped, out of %d", i, (Int_t) decayed->Integral(), (Int_t) stopped->Integral(), (Int_t) genePt->Integral());
1658 decayed->Divide(decayed, genePt, 1, 1, "B");
1659 stopped->Divide(stopped, genePt, 1, 1, "B");
1661 decayed->SetMarkerStyle(20);
1662 stopped->SetMarkerStyle(21);
1663 stopped->SetMarkerColor(2);
1665 new TCanvas(Form("all_%d_%d", loop, i), Form("all_%d_%d", loop, i), 600, 600);
1667 decayed->DrawCopy("SAME");
1668 stopped->DrawCopy("SAME");
1670 decayed->Add(stopped);
1671 decayed->Add(effPt);
1672 decayed->SetMarkerStyle(22);
1673 decayed->SetMarkerColor(4);
1674 decayed->DrawCopy("SAME");
1679 Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1681 Printf("%f measured, %f generated, effiency: %f", sumGen, sumMeas, sumMeas / sumGen);
1688 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
1692 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1697 canvas3->SaveAs(Form("%s.eps", canvas3->GetName()));
1703 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt_ref.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1704 mv unfolded.root chi2_ptref.root
1705 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1706 mv unfolded.root chi2_pt0.root
1707 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1708 mv unfolded.root chi2_pt1.root
1709 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1710 mv unfolded.root chi2_pt0_25.root
1711 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1712 mv unfolded.root chi2_pt1_25.root
1718 const char* files[] = { "chi2_ptref.root", "chi2_pt0.root", "chi2_pt1.root", "chi2_pt0_25.root", "chi2_pt1_25.root"};
1721 for (Int_t i = 0; i<nMax; ++i)
1723 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(files[i]);
1724 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1727 const char* legendStrings[] = { "Reduced 50%", "Enhanced 50%", "Reduced 25%", "Enhanced 25%" };
1728 DrawRatio(results[0], nMax-1, results+1, "LowMomentumSyst.eps", kFALSE, legendStrings);
1731 void ParticleSpeciesComparison()
1738 // loop over cases (normal, enhanced/reduced ratios)
1740 for (Int_t i = 0; i<nMax; ++i)
1742 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(Form("chi2_species_%d.root", i), Form("Multiplicity_%d", i));
1744 mc = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymchist", 1, 1);
1745 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1748 DrawResultRatio(mc, results[0], "ParticleSpeciesComparison1_1.eps");
1750 for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1752 results[0]->SetBinError(i, 0);
1753 mc->SetBinError(i, 0);
1756 const char* legendStrings[] = { "K #times 0.5", "K #times 1.5", "p #times 0.5", "p #times 1.5", "K #times 0.5, p #times 0.5", "K #times 1.5, p #times 1.5", "K #times 0.5, p #times 1.5", "K #times 1.5, p #times 0.5" };
1758 DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
1760 //not valid: draw chi2 uncertainty on top!
1761 /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1762 TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1763 errorHist->SetLineColor(1);
1764 errorHist->SetLineWidth(2);
1765 TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1766 for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1768 errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1769 errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1772 errorHist->DrawCopy("SAME");
1773 errorHist2->DrawCopy("SAME");*/
1775 //canvas->SaveAs(canvas->GetName());
1777 DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
1779 //errorHist->DrawCopy("SAME");
1780 //errorHist2->DrawCopy("SAME");
1782 //canvas2->SaveAs(canvas2->GetName());
1785 /*void ParticleSpeciesComparison2()
1787 gSystem->Load("libPWG0base");
1789 const char* fileNameMC = "multiplicityMC_400k_syst.root";
1790 const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1793 TFile::Open(fileNameMC);
1794 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1795 mult->LoadHistograms();
1800 // loop over cases (normal, enhanced/reduced ratios)
1802 for (Int_t i = 0; i<nMax; ++i)
1805 folder.Form("Multiplicity_%d", i);
1807 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1809 TFile::Open(fileNameESD);
1810 mult2->LoadHistograms();
1812 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1816 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1817 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1818 //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1822 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1823 //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1826 //Float_t averageRatio = 0;
1827 //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1829 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1831 TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1832 mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1834 //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1835 //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1837 //Printf("Case %d. Average ratio is %f", i, averageRatio);
1840 DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1843 TH1* Invert(TH1* eff)
1845 // calculate corr = 1 / eff
1847 TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1850 for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1852 if (eff->GetBinContent(i) > 0)
1854 corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1855 corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1862 void TriggerVertexCorrection()
1865 // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1870 TFile::Open(correctionFile);
1871 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1872 mult->LoadHistograms("Multiplicity");
1874 TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1875 TH1* corrNSD = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kNSD));
1876 TH1* corrMB = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1878 TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 500);
1881 gPad->SetTopMargin(0.05);
1882 gPad->SetRightMargin(0.05);
1884 corrINEL->SetStats(kFALSE);
1885 corrINEL->GetXaxis()->SetRangeUser(0, 12);
1886 corrINEL->GetYaxis()->SetRangeUser(0, 8);
1887 corrINEL->SetTitle(Form(";%s;Correction factor", GetMultLabel()));
1888 corrINEL->SetMarkerStyle(22);
1889 corrINEL->Draw("PE");
1891 corrMB->SetStats(kFALSE);
1892 corrMB->SetLineColor(2);
1893 corrMB->SetMarkerStyle(25);
1894 corrMB->SetMarkerColor(2);
1895 corrMB->Draw("SAME PE");
1897 corrNSD->SetLineColor(4);
1898 corrNSD->SetMarkerStyle(24);
1899 corrNSD->SetMarkerColor(4);
1900 corrNSD->Draw("SAME PE");
1902 Printf(" MB INEL NSD");
1903 Printf("bin 0: %f %f %f", corrMB->GetBinContent(1), corrINEL->GetBinContent(1), corrNSD->GetBinContent(1));
1904 Printf("bin 1: %f %f %f", corrMB->GetBinContent(2), corrINEL->GetBinContent(2), corrNSD->GetBinContent(2));
1906 TLegend* legend = new TLegend(0.3, 0.6, 0.85, 0.85);
1907 legend->SetFillColor(0);
1908 legend->AddEntry(corrINEL, "Correction to inelastic sample");
1909 legend->AddEntry(corrNSD, "Correction to NSD sample");
1910 legend->AddEntry(corrMB, "Correction to triggered sample");
1911 legend->SetTextSize(0.04);
1915 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1918 void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
1922 TFile::Open(correctionFile);
1923 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1924 mult->LoadHistograms("Multiplicity");
1926 TFile::Open(measuredFile);
1927 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1928 mult2->LoadHistograms("Multiplicity");
1930 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1932 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1934 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 1e5);
1936 TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1940 TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1942 TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
1946 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1947 DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
1950 TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
1951 errorResponse->Write();
1952 errorMeasured->Write();
1957 void DrawStatisticalUncertainty()
1959 TFile::Open("StatisticalUncertainty.root");
1961 errorResponse = (TH1*) gFile->Get("errorResponse");
1962 errorMeasured = (TH1*) gFile->Get("errorMeasured");
1963 errorBoth = (TH1*) gFile->Get("errorBoth");
1965 TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
1968 canvas->SetRightMargin(0.05);
1969 canvas->SetTopMargin(0.05);
1971 errorResponse->SetLineColor(1);
1972 errorResponse->GetXaxis()->SetRangeUser(0, longDisplayRange);
1973 errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1974 errorResponse->SetStats(kFALSE);
1975 errorResponse->SetTitle(";true multiplicity;Uncertainty");
1977 errorResponse->Draw();
1979 errorMeasured->SetLineColor(2);
1980 errorMeasured->Draw("SAME");
1982 errorBoth->SetLineColor(4);
1983 errorBoth->Draw("SAME");
1985 Printf("Average errorResponse: %f", errorResponse->Integral(2, displayRange) / (displayRange - 1));
1986 Printf("Average errorMeasured: %f", errorMeasured->Integral(2, displayRange) / (displayRange - 1));
1987 Printf("Average errorBoth: %f", errorBoth->Integral(2, displayRange) / (displayRange - 1));
1989 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1992 void StatisticalUncertaintyCompare(const char* det = "SPD")
1994 TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1995 TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1996 TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1997 TH1* errorBoth = (TH1*) file1->Get("errorBoth");
2000 str.Form("StatisticalUncertaintyCompare%s", det);
2002 TCanvas* canvas = new TCanvas(str, str, 800, 500);
2005 canvas->SetRightMargin(0.05);
2006 canvas->SetTopMargin(0.05);
2008 errorResponse->Scale(1.0 / sqrt(2));
2009 errorMeasured->Scale(1.0 / sqrt(2));
2010 errorBoth->Scale(1.0 / sqrt(2));
2012 errorResponse->SetLineColor(1);
2013 errorResponse->GetXaxis()->SetRangeUser(0, displayRange);
2014 errorResponse->GetYaxis()->SetRangeUser(0, 0.18);
2015 errorResponse->SetStats(kFALSE);
2016 errorResponse->GetYaxis()->SetTitleOffset(1.2);
2017 errorResponse->SetTitle(Form(";%s;#sqrt{2}^{-1} #sigma(unfolded - unfolded_{0}) / unfolded_{0}", GetMultLabel()));
2019 errorResponse->Draw();
2021 errorMeasured->SetLineColor(2);
2022 errorMeasured->Draw("SAME");
2024 errorBoth->SetLineColor(4);
2025 errorBoth->Draw("SAME");
2027 TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
2028 TH1* errorResponse2 = (TH1*) file2->Get("errorResponse");
2029 TH1* errorMeasured2 = (TH1*) file2->Get("errorMeasured");
2030 TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
2032 errorResponse2->Scale(1.0 / sqrt(2));
2033 errorMeasured2->Scale(1.0 / sqrt(2));
2034 errorBoth2->Scale(1.0 / sqrt(2));
2036 errorResponse2->SetLineStyle(2);
2037 errorResponse2->Draw("SAME");
2039 errorMeasured2->SetLineColor(2);
2040 errorMeasured2->SetLineStyle(2);
2041 errorMeasured2->Draw("SAME");
2043 errorBoth2->SetLineColor(4);
2044 errorBoth2->SetLineStyle(2);
2045 errorBoth2->Draw("SAME");
2047 TLegend* legend = new TLegend(0.2, 0.5, 0.8, 0.9);
2048 legend->SetFillColor(0);
2049 legend->SetTextSize(0.04);
2050 legend->AddEntry(errorBoth, "Both (Bayesian unfolding)");
2051 legend->AddEntry(errorMeasured, "Measured (Bayesian unfolding)");
2052 legend->AddEntry(errorResponse, "Response matrix (Bayesian unfolding)");
2053 legend->AddEntry(errorBoth2, "Both (#chi^{2}-minimization)");
2054 legend->AddEntry(errorMeasured2, "Measured (#chi^{2}-minimization)");
2055 legend->AddEntry(errorResponse2, "Response matrix (#chi^{2}-minimization)");
2058 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2061 void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
2063 const char* files[] = { "multiplicityMC_nd.root", "multiplicityMC_sd.root", "multiplicityMC_dd.root", "multiplicityMC_xsection.root" };
2067 TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 600);
2070 canvas->SetRightMargin(0.05);
2071 canvas->SetTopMargin(0.05);
2073 AliMultiplicityCorrection* data[4];
2075 TH1* effErrorArray[2];
2077 Int_t markers[] = { 24, 25, 26, 5 };
2078 //Int_t markers[] = { 2, 25, 24, 5 };
2079 Int_t colors[] = { 1, 2, 4, 6 };
2080 //Int_t colors[] = { 1, 1, 1, 1 };
2082 //TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
2083 TLegend* legend = new TLegend(0.3, 0.3, 0.9, 0.6);
2084 legend->SetTextSize(0.04);
2085 legend->SetFillColor(0);
2087 for (Int_t i=0; i<4; ++i)
2090 name.Form("Multiplicity_%d", i);
2092 TFile::Open(files[i]);
2093 data[i] = new AliMultiplicityCorrection(name, name);
2097 data[i]->LoadHistograms("Multiplicity");
2100 data[i]->LoadHistograms("Multiplicity_0");
2103 if (eventType == -1)
2105 eff = (TH1*) data[i]->GetTriggerEfficiency(etaRange)->Clone(Form("eff_%d", i));
2108 eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
2111 eff->GetXaxis()->SetRangeUser(0, 15);
2112 eff->GetYaxis()->SetRangeUser(0, 1.19);
2113 eff->SetStats(kFALSE);
2114 eff->GetXaxis()->SetTitle(GetMultLabel());
2115 eff->GetYaxis()->SetTitle("Efficiency");
2117 eff->SetLineColor(colors[i]);
2118 eff->SetMarkerColor(colors[i]);
2119 eff->SetMarkerStyle(markers[i]);
2123 // once for INEL, once for NSD
2124 for (AliMultiplicityCorrection::EventType eventType2 = AliMultiplicityCorrection::kINEL; eventType2 <= AliMultiplicityCorrection::kNSD; eventType2++)
2126 effDiff = (TH1*) data[i]->GetEfficiency(etaRange, eventType2)->Clone(Form("effDiff_%d", i));
2128 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2129 effDiff->SetBinError(bin, 0);
2131 // loop over cross section combinations
2132 for (Int_t j=1; j<7; ++j)
2134 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
2135 mult->LoadHistograms(Form("Multiplicity_%d", j));
2137 TH1* eff2 = mult->GetEfficiency(etaRange, eventType2);
2139 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2141 // TODO we could also do asymmetric errors here
2142 Float_t deviation = TMath::Abs(effDiff->GetBinContent(bin) - eff2->GetBinContent(bin));
2144 effDiff->SetBinError(bin, TMath::Max(effDiff->GetBinError(bin), (Double_t) deviation));
2148 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2150 //if (eventType2 == AliMultiplicityCorrection::kINEL)
2151 //eff->SetBinError(bin, 0);
2152 //eff->SetBinError(bin, effDiff->GetBinError(bin));
2153 if (bin < 20 && effDiff->GetBinContent(bin) > 0)
2154 Printf("Bin %d: Error: %.2f", bin, 100.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin));
2158 TH1* effError = (TH1*) effDiff->Clone(Form("effError_%s", (eventType2 == AliMultiplicityCorrection::kINEL) ? "INEL" : "NSD"));
2159 effErrorArray[eventType2 - AliMultiplicityCorrection::kINEL] = effError;
2162 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2163 if (effDiff->GetBinContent(bin) > 0)
2164 effError->SetBinContent(bin, 1.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin));
2166 effError->SetLineColor(1);
2167 effError->SetMarkerStyle(1);
2169 if (eventType2 == AliMultiplicityCorrection::kNSD)
2170 effError->SetLineStyle(2);
2172 effError->DrawCopy("SAME HIST");
2183 eff->DrawCopy("SAME P");
2185 legend->AddEntry(eff, (((i == 0) ? "Non-diffractive" : ((i == 1) ? "Single-diffractive" : ((i == 2) ? "Double-diffractive" : "Pythia combined")))));
2190 legend->AddEntry(effErrorArray[0], "Relative syst. uncertainty: inelastic");
2191 legend->AddEntry(effErrorArray[1], "Relative syst. uncertainty: NSD");
2193 file = TFile::Open("uncertainty_xsection.root", "RECREATE");
2194 effErrorArray[0]->Write();
2195 effErrorArray[1]->Write();
2201 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2204 void ModelDependencyPlot()
2208 TFile::Open("multiplicityMC.root");
2209 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2210 mult->LoadHistograms("Multiplicity");
2212 hist = mult->GetCorrelation(etaRange);
2214 for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
2216 for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
2218 hist->SetBinContent(0, y, z, 0);
2219 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2223 TH2* proj = (TH2*) hist->Project3D("zy");
2225 TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 1200, 600);
2227 canvas->Divide(2, 1);
2233 gPad->SetTopMargin(0.05);
2234 gPad->SetRightMargin(0.05);
2236 Int_t selectedMult = 30;
2239 TH1* full = proj->ProjectionX("full");
2240 TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
2242 full->SetStats(kFALSE);
2243 full->GetXaxis()->SetRangeUser(0, displayRange);
2244 full->GetYaxis()->SetRangeUser(5, yMax);
2245 full->SetTitle(";Multiplicity;Entries");
2247 selected->SetLineColor(0);
2248 selected->SetMarkerColor(2);
2249 selected->SetMarkerStyle(5);
2252 selected->Draw("SAME P");
2254 TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
2255 legend->SetFillColor(0);
2256 legend->AddEntry(full, "True");
2257 legend->AddEntry(selected, "Measured");
2260 TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
2261 line->SetLineWidth(2);
2268 gPad->SetTopMargin(0.05);
2269 gPad->SetRightMargin(0.05);
2271 full = proj->ProjectionY("full2");
2272 selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
2274 full->SetStats(kFALSE);
2275 full->GetXaxis()->SetRangeUser(0, displayRange);
2276 full->GetYaxis()->SetRangeUser(5, yMax);
2277 full->SetTitle(";Multiplicity;Entries");
2279 full->SetLineColor(0);
2280 full->SetMarkerColor(2);
2281 full->SetMarkerStyle(5);
2284 selected->Draw("SAME");
2288 line = new TLine(selectedMult, 5, selectedMult, yMax);
2289 line->SetLineWidth(2);
2292 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2295 void SystematicpTSpectrum()
2297 gSystem->Load("libPWG0base");
2299 TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
2300 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2301 mult->LoadHistograms("Multiplicity");
2303 TFile::Open("multiplicityMC_100k_syst.root");
2304 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2305 mult2->LoadHistograms("Multiplicity");
2307 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2308 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2309 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2311 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2312 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2314 DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
2317 void FitPt(const char* fileName)
2319 // needs a MC file from the dNdEta analysis
2321 TFile::Open(fileName);
2323 TH1* genePt = (TH1*) gFile->Get("fdNdpT");
2325 genePt->SetTitle(";p_{T} (GeV/c);1/p_{T} dN_{ch}/dp_{T} (GeV/c)^{-2}");
2327 genePt->Scale(1.0 / 287800);
2329 genePt->Scale(1.0 / genePt->GetXaxis()->GetBinWidth(1));
2331 genePt->GetXaxis()->SetRangeUser(0, 0.4);
2333 TF1* func = new TF1("func", "[1]*x*exp(x*[0])");
2334 func->SetParameters(-1, 1);
2336 genePt->SetMarkerStyle(25);
2337 genePt->SetTitle("");
2338 genePt->SetStats(kFALSE);
2341 genePt->DrawCopy("P");
2342 func->DrawCopy("SAME");
2345 genePt->Fit(func, "0", "", 0, 0.25);
2346 genePt->Fit(func, "0", "", 0, 0.25);
2348 TCanvas* canvas = new TCanvas("FitPt", "FitPt", 600, 600);
2352 gPad->SetLeftMargin(0.13);
2353 gPad->SetRightMargin(0.05);
2354 gPad->SetTopMargin(0.05);
2356 //genePt->GetXaxis()->SetRangeUser(0, 1);
2357 genePt->GetYaxis()->SetRangeUser(2, 200);
2358 genePt->GetYaxis()->SetTitleOffset(1.4);
2359 genePt->GetXaxis()->SetTitleOffset(1.1);
2360 genePt->DrawCopy("P");
2361 //func->SetRange(0, 0.3);
2362 func->DrawCopy("SAME");
2365 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2367 TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2369 TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2371 TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2373 for (Int_t param=0; param<2; param++)
2375 for (Int_t sign=0; sign<2; sign++)
2377 TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign)); // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2378 func2->SetParameters(func->GetParameters());
2379 //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2381 Float_t factor = ((sign == 0) ? 0.75 : 1.25);
2382 func2->SetParameter(param, func2->GetParameter(param) * factor);
2386 func2->SetLineWidth(2);
2387 func2->SetLineColor(2);
2388 func2->SetLineStyle(2);
2389 func2->DrawCopy("SAME");
2392 TH1* second = func2->GetHistogram();
2393 second->Divide(first);
2394 second->SetLineColor(param + 1);
2395 // set to 1 above 0.2 GeV
2396 for (Int_t bin=second->GetXaxis()->FindBin(0.20001); bin<=second->GetNbinsX(); bin++)
2397 second->SetBinContent(bin, 1);
2398 second->GetYaxis()->SetRangeUser(0, 2);
2399 second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2400 second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2404 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2405 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2410 void DrawSystematicpT()
2412 TFile* file = TFile::Open("SystematicpT.root");
2414 TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2415 TH1* result2 = (TH1*) file->Get("result_unity");
2422 for (Int_t id=0; id<nParams*2; ++id)
2424 mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2425 result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2428 DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2430 //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2432 DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2434 //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2436 // does not make sense: mc is different
2437 //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2440 void SystematicpT(Bool_t chi2 = 1)
2442 gSystem->Load("libPWG0base");
2444 TFile::Open("ptspectrum900.root");
2445 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2446 mult->LoadHistograms("Multiplicity");
2448 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2455 for (Int_t param=0; param<nParams; param++)
2457 for (Int_t sign=0; sign<2; sign++)
2459 // calculate result with systematic effect
2460 TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
2461 mult2->LoadHistograms("Multiplicity");
2463 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2467 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2468 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2471 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2473 Int_t id = param * 2 + sign;
2475 mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2476 result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2478 TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2479 DrawResultRatio(mcHist[id], result[id], tmp);
2483 // calculate normal result
2484 TFile::Open("ptspectrum100_1.root");
2485 mult2->LoadHistograms("Multiplicity");
2486 TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
2488 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2492 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2495 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2497 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
2499 TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2502 for (Int_t id=0; id<nParams*2; ++id)
2504 mcHist[id]->Write();
2505 result[id]->Write();
2512 void DrawSystematicpT2()
2514 //displayRange = 200;
2517 TFile* file = TFile::Open("SystematicpT2.root");
2518 TH1* mcHist = (TH1*) file->Get("mymc");
2520 result[0] = (TH1*) file->Get("result_unity");
2522 for (Int_t id=0; id<nParams*2; ++id)
2523 result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
2525 DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
2526 DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
2527 DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
2530 void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
2532 gSystem->Load("libPWG0base");
2537 TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2540 TFile::Open("ptspectrum100_1.root");
2542 AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2543 measured->LoadHistograms("Multiplicity");
2544 TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2550 // -1 = unity change, 0...4 parameters
2551 for (Int_t id=-1; id<nParams*2; id++)
2553 Int_t param = id / 2;
2554 Int_t sign = id % 2;
2562 idStr.Form("%d_%d", param, sign);
2564 // calculate result with systematic effect
2567 TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
2570 TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
2572 AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2573 response->LoadHistograms("Multiplicity");
2575 response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
2579 response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2580 response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2583 response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2585 result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
2587 TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
2588 DrawResultRatio(mcHist, result[id+1], tmp);
2591 TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
2593 for (Int_t id=0; id<nParams*2+1; ++id)
2594 result[id]->Write();
2597 DrawSystematicpT2();
2600 void SystematicpTCutOff(Bool_t chi2 = kTRUE)
2602 // only needed for TPC
2605 gSystem->Load("libPWG0base");
2607 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2608 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2609 mult->LoadHistograms("Multiplicity");
2611 TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2612 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2613 mult2->LoadHistograms("Multiplicity");
2616 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2620 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2621 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2624 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2626 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2627 TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
2631 // change of pt spectrum (down)
2632 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2633 AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2634 mult3->LoadHistograms("Multiplicity");
2635 mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2638 mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2639 mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2642 mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2643 syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2645 // change of pt spectrum (up)
2646 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
2647 AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
2648 mult4->LoadHistograms("Multiplicity");
2649 mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2652 mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2653 mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2656 mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2657 syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
2659 DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
2661 Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
2662 Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
2665 TH1* SystematicsSummary(Bool_t tpc = 0, Bool_t nsd = kTRUE)
2670 const char** names = 0;
2671 Int_t colors[] = { 1, 2, 4, 1, 2, 4 };
2672 Int_t styles[] = { 1, 2, 3, 1, 2, 3 };
2673 Int_t widths[] = { 1, 1, 1, 2, 2, 2 };
2674 Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
2676 TH1* dummy = new TH2F("dummy", Form(";%s;Uncertainty", GetMultLabel()), 202, -1.5, 200.5, 100, 0, 0.4);
2679 for (Int_t i=0; i<nEffects; ++i)
2680 effects[i] = new TH1F("SystematicsSummary", Form(";%s;Uncertainty", GetMultLabel()), 201, -0.5, 200.5);
2686 const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
2690 TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
2691 TH1* hist = (TH1*) file->Get("errorBoth");
2693 // smooth a bit, but skip 0 bin
2694 effects[0]->SetBinContent(2, hist->GetBinContent(2));
2695 for (Int_t i=3; i<=200; ++i)
2696 effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2698 // relative x-section
2699 effects[1]->SetBinContent(2, 0.005);
2700 effects[1]->SetBinContent(3, 0.0025);
2701 effects[1]->SetBinContent(4, 0.0025);
2703 // particle composition
2704 for (Int_t i=2; i<=101; ++i)
2708 effects[2]->SetBinContent(i, 0.01);
2712 effects[2]->SetBinContent(i, 0.02);
2715 effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
2718 // pt cut off (only tpc)
2719 for (Int_t i=2; i<=101; ++i)
2723 effects[3]->SetBinContent(i, 0.05);
2727 effects[3]->SetBinContent(i, 0.01);
2730 effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
2733 // track selection (only tpc)
2734 for (Int_t i=2; i<=101; ++i)
2735 effects[4]->SetBinContent(i, 0.03);
2738 for (Int_t i=2; i<=101; ++i)
2739 effects[5]->SetBinContent(i, 0.01);
2742 for (Int_t i=2; i<=101; ++i)
2746 effects[6]->SetBinContent(i, 0.05);
2750 effects[6]->SetBinContent(i, 0.02);
2753 effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
2761 //const char* namesSPD[] = { "Particle composition", "p_{t} cut-off", "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)"};
2762 const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)", "Particle composition", "p_{t} cut-off"};
2768 TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
2769 TH1* hist = (TH1*) file->Get("errorBoth");
2770 hist->Scale(1.0 / sqrt(2));
2772 // smooth a bit, but skip 0 bin
2773 /*effects[currentEffect]->SetBinContent(1, hist->GetBinContent(1));
2774 for (Int_t i=2; i<=201; ++i)
2775 effects[currentEffect]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);*/
2776 effects[currentEffect] = hist;
2780 // relative x-section
2781 file = TFile::Open("uncertainty_xsection.root");
2782 effects[currentEffect++] = (TH1*) file->Get("effError_INEL");
2783 effects[currentEffect] = (TH1*) file->Get("effError_NSD");
2784 effects[currentEffect]->SetLineStyle(1);
2785 //effects[2]->SetBinContent(1, 0.20);
2786 //effects[2]->SetBinContent(2, 0.01);
2787 //effects[2]->SetBinContent(3, 0.002);
2791 // particle composition
2792 effects[currentEffect]->SetBinContent(1, 0.16);
2793 for (Int_t i=2; i<=81; ++i)
2795 effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * i / 81);
2801 effects[currentEffect]->SetBinContent(1, 0.06);
2802 effects[currentEffect]->SetBinContent(2, 0.03);
2803 for (Int_t i=3; i<=81; ++i)
2807 effects[currentEffect]->SetBinContent(i, 0.01);
2811 effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * (i - 61) / 20);
2817 // // material budget
2818 // for (Int_t i=1; i<=81; ++i)
2821 // effects[currentEffect]->SetBinContent(i, 0.05 - 0.01 * i);
2823 // effects[currentEffect]->SetBinContent(i, 0.05 * (i - 50) / 30);
2830 TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 500);
2831 canvas->SetRightMargin(0.05);
2832 canvas->SetTopMargin(0.05);
2833 //canvas->SetGridx();
2835 TLegend* legend = new TLegend(0.2, 0.4, 0.7, 0.4 + 0.5 * nEffects / 7);
2836 legend->SetFillColor(0);
2837 legend->SetTextSize(0.04);
2839 dummy->GetXaxis()->SetRangeUser(0, displayRange);
2841 for (Int_t i=0; i<nEffects; ++i)
2843 TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
2845 for (Int_t j=0; j<nEffects-i; ++j)
2846 current->Add(effects[j]);*/
2848 current->SetLineColor(colors[i]);
2849 current->SetLineStyle(styles[i]);
2850 current->SetLineWidth(widths[i]);
2851 //current->SetFillColor(colors[i]);
2852 current->SetMarkerColor(colors[i]);
2853 //current->SetMarkerStyle(markers[i]);
2855 current->SetStats(kFALSE);
2856 current->GetYaxis()->SetRangeUser(0, 0.4);
2857 current->DrawCopy("SAME");
2858 legend->AddEntry(current, names[i]);
2860 //TLatex* text = new TLatex(displayRange+2, current->GetBinContent(displayRange+1), names[i]);
2861 TLatex* text = new TLatex(displayRange+2, 0.1 - i * 0.02, names[i]);
2862 text->SetTextSize(0.04);
2863 text->SetTextColor(colors[i]);
2867 // add total in square
2868 TH1* totalINEL = (TH1*) effects[0]->Clone("totalINEL");
2870 TH1* totalNSD = (TH1*) totalINEL->Clone("totalNSD");
2872 for (Int_t i=0; i<nEffects; ++i)
2874 //Printf("%d %f", i, effects[i]->GetBinContent(20));
2875 effects[i]->Multiply(effects[i]);
2878 totalINEL->Add(effects[i]);
2880 totalNSD->Add(effects[i]);
2883 for (Int_t i=1; i<=totalINEL->GetNbinsX(); ++i)
2885 totalINEL->SetBinError(i, 0);
2886 if (totalINEL->GetBinContent(i) > 0)
2887 totalINEL->SetBinContent(i, TMath::Min(sqrt(totalINEL->GetBinContent(i)), 1.0));
2888 totalNSD->SetBinError(i, 0);
2889 if (totalNSD->GetBinContent(i) > 0)
2890 totalNSD->SetBinContent(i, TMath::Min(sqrt(totalNSD->GetBinContent(i)), 1.0));
2893 //Printf("%f", total->GetBinContent(20));
2895 totalINEL->SetMarkerStyle(5);
2896 totalINEL->SetMarkerColor(1);
2897 legend->AddEntry(totalINEL, "Total (INEL)", "P");
2899 totalNSD->SetMarkerStyle(24);
2900 totalNSD->SetMarkerColor(2);
2901 legend->AddEntry(totalNSD, "Total (NSD)", "P");
2903 Printf("total in bin 0 is INEL: %f NSD: %f", totalINEL->GetBinContent(1), totalNSD->GetBinContent(1));
2904 totalINEL->DrawCopy("SAME P"); //->SetBinContent(1, 0);
2905 totalNSD->DrawCopy("SAME P"); //->SetBinContent(1, 0);
2909 canvas->SaveAs(canvas->GetName());
2911 return (nsd) ? totalNSD : totalINEL;
2914 void finalPlot(Bool_t tpc = 0, Bool_t small = kFALSE)
2921 //TH1* errorNSD = SystematicsSummary(tpc, 1);
2923 TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", 800, 500);
2924 canvas->SetRightMargin(0.05);
2925 canvas->SetTopMargin(0.05);
2929 legend = new TLegend(0.5, 0.6, 0.9, 0.8);
2930 legend->SetFillColor(0);
2931 legend->SetTextSize(0.04);
2933 for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++)
2935 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root");
2936 TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc");
2937 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2939 DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType));
2942 //result->Scale(1.0 / result->Integral(2, displayRange));
2944 result->GetXaxis()->SetRangeUser(0, displayRange);
2945 //result->SetBinContent(1, 0); result->SetBinError(1, 0);
2946 result->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (etaRange+1) * 0.5));
2947 result->SetMarkerStyle(0);
2948 result->SetLineColor(1);
2949 result->SetStats(kFALSE);
2952 TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD));
2954 TH1* systError = (TH1*) result->Clone("systError");
2955 for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
2956 systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
2958 // change error drawing style
2959 systError->SetFillColor(15);
2961 if (eventType == AliMultiplicityCorrection::kNSD)
2963 result->SetLineColor(2);
2964 result->SetMarkerColor(2);
2965 result->SetMarkerStyle(5);
2969 systError->DrawCopy(Form("E2 ][ %s", (eventType == AliMultiplicityCorrection::kINEL) ? "" : "SAME"));
2970 result->DrawCopy("SAME E ][");
2973 legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic cross-section" : "NSD cross-section", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P");
2978 //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
2979 TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
2980 text->SetFillColor(0);
2981 text->SetTextAlign(12);
2982 text->AddText("Systematic errors summed quadratically");
2983 text->AddText("0.6 million minimum bias events");
2984 text->AddText("corrected to inelastic events");
2987 TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
2988 text2->SetFillColor(0);
2989 text2->SetTextAlign(12);
2990 text2->AddText("#sqrt{s} = 14 TeV");
2993 text2->AddText("|#eta| < 0.9");
2996 text2->AddText("|#eta| < 2.0");
2997 text2->AddText("simulated data (PYTHIA)");
3002 TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
3008 TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
3014 TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
3017 TImage* img = TImage::Open("$HOME/alice.png");
3018 img->SetImageQuality(TAttImage::kImgBest);
3024 /* TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
3025 text->SetTextSize(0.04);
3026 text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
3027 text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
3028 text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
3032 canvas->SaveAs(canvas->GetName());
3035 void finalPlot2(Bool_t tpc = 0)
3042 //TH1* errorNSD = SystematicsSummary(tpc, 1);
3044 TCanvas* canvas = new TCanvas("finalPlot2.eps", "finalPlot2.eps", 800, 600);
3045 canvas->SetRightMargin(0.05);
3046 canvas->SetTopMargin(0.05);
3051 legend = new TLegend(0.5, 0.7, 0.9, 0.9);
3052 legend->SetFillColor(0);
3053 legend->SetTextSize(0.04);
3055 Int_t displayRanges[] = { 50, 80, 120 };
3057 TH2* dummy = new TH2F("dummy", ";True multiplicity N_{ch};P(N_{ch})", 100, -0.5, displayRanges[2]+10, 1000, 5e-5, 5);
3064 for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++)
3066 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root");
3067 //TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc");
3068 for (Int_t etaR = 0; etaR <= 2; etaR++)
3070 TH1* result = mult->GetMultiplicityESDCorrected(etaR);
3072 //DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType));
3075 result->Scale(TMath::Power(5, etaR) / result->Integral(1, displayRanges[etaR]));
3077 result->GetXaxis()->SetRangeUser(0, displayRanges[etaR]);
3078 //result->SetBinContent(1, 0); result->SetBinError(1, 0);
3079 //result->SetTitle(Form(""));
3080 result->SetMarkerStyle(0);
3081 result->SetLineColor(1);
3082 result->SetLineWidth(2);
3083 //result->SetMarkerStyle(4);
3084 //result->SetStats(kFALSE);
3087 TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD));
3089 TH1* systError = (TH1*) result->Clone("systError");
3090 // small hack until we have syst. errors for all eta ranges
3091 Float_t factor = 80.0 / displayRanges[etaR];
3092 for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
3094 systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(factor * i));
3097 // change error drawing style
3098 systError->SetFillColor(15);
3100 if (eventType == AliMultiplicityCorrection::kNSD)
3102 result->SetLineColor(2);
3103 result->SetMarkerColor(2);
3104 result->SetMarkerStyle(5);
3108 systError->GetXaxis()->SetRangeUser(0, displayRanges[etaR]);
3109 systError->DrawCopy("E2 ][ SAME");
3112 if (eventType == AliMultiplicityCorrection::kINEL)
3114 TLatex* Tl = new TLatex;
3115 Tl->SetTextSize(0.04);
3116 //Tl->SetBit(TLatex::kTextNDC);
3117 Tl->SetTextAlign(32);
3119 Float_t etaRangeArr[] = { 0.5, 1.0, 1.4 };
3121 tmpStr.Form("|#eta| < %.1f (x %d)", etaRangeArr[etaR], (Int_t) TMath::Power(5, etaR));
3123 Tl->DrawLatex(36, result->GetBinContent(41), tmpStr);
3126 Tl->SetTextAlign(12);
3127 Tl->DrawLatex(82, result->GetBinContent(81), tmpStr);
3131 Tl->SetTextAlign(12);
3132 Tl->DrawLatex(106, result->GetBinContent(101), tmpStr);
3137 legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic events" : "NSD events", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P");
3143 for (Int_t i=0; i<list.GetEntries(); i++)
3144 ((TH1*) list.At(i))->DrawCopy("SAME E ][");
3148 canvas->SaveAs(canvas->GetName());
3151 TMatrixD* NonInvertable()
3153 const Int_t kSize = 5;
3155 TMatrixD matrix(kSize, kSize);
3156 for (Int_t x=0; x<kSize; x++)
3158 for (Int_t y=0; y<kSize; y++)
3162 if (x == 0 || x == kSize -1)
3164 matrix(x, y) = 0.75;
3169 else if (TMath::Abs(x - y) == 1)
3171 matrix(x, y) = 0.25;
3178 //TMatrixD inverted(matrix);
3179 //inverted.Invert();
3183 return new TMatrixD(matrix);
3186 void BlobelUnfoldingExample()
3188 const Int_t kSize = 20;
3190 TMatrixD matrix(kSize, kSize);
3191 for (Int_t x=0; x<kSize; x++)
3193 for (Int_t y=0; y<kSize; y++)
3197 if (x == 0 || x == kSize -1)
3199 matrix(x, y) = 0.75;
3204 else if (TMath::Abs(x - y) == 1)
3206 matrix(x, y) = 0.25;
3213 TMatrixD inverted(matrix);
3218 TH1F* inputDist = new TH1F("inputDist", ";t;Entries", kSize, -0.5, (Float_t) kSize - 0.5);
3219 TVectorD inputDistVector(kSize);
3220 TH1F* unfolded = inputDist->Clone("unfolded");
3221 TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
3222 measuredIdealDist->SetTitle(";m;Entries");
3223 TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
3225 TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
3226 // norm: 1/(sqrt(2pi)sigma)
3227 gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
3230 for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
3232 Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
3233 inputDist->SetBinContent(x, value);
3234 inputDistVector(x-1) = value;
3237 TVectorD measuredDistIdealVector = matrix * inputDistVector;
3239 for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
3240 measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
3242 measuredDist->FillRandom(measuredIdealDist, 10000);
3244 // fill error matrix before scaling
3245 TMatrixD covarianceMatrixMeasured(kSize, kSize);
3246 for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
3247 covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
3249 TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
3250 //covarianceMatrix.Print();
3252 TVectorD measuredDistVector(kSize);
3253 for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
3254 measuredDistVector(x-1) = measuredDist->GetBinContent(x);
3256 TVectorD unfoldedVector = inverted * measuredDistVector;
3257 for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
3258 unfolded->SetBinContent(x, unfoldedVector(x-1));
3260 TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1200, 600);
3261 canvas->SetTopMargin(0.05);
3262 canvas->Divide(2, 1);
3265 canvas->cd(1)->SetLeftMargin(0.15);
3266 canvas->cd(1)->SetRightMargin(0.05);
3267 canvas->cd(1)->SetTopMargin(0.05);
3270 measuredDist->GetYaxis()->SetRangeUser(-600, 2799);
3271 measuredDist->GetYaxis()->SetTitleOffset(1.9);
3272 measuredDist->SetStats(0);
3273 measuredDist->DrawCopy();
3277 canvas->cd(2)->SetLeftMargin(0.15);
3278 canvas->cd(2)->SetRightMargin(0.05);
3279 canvas->cd(2)->SetTopMargin(0.05);
3282 unfolded->GetYaxis()->SetRangeUser(-600, 2799);
3283 unfolded->GetYaxis()->SetTitleOffset(1.9);
3284 unfolded->SetStats(0);
3285 unfolded->DrawCopy();
3288 canvas->SaveAs("BlobelUnfoldingExample.eps");
3292 // now unfold this with Bayesian
3295 // fill a multiplicity object
3296 mult = new AliMultiplicityCorrection("mult", "mult");
3297 for (Int_t x=0; x<kSize; x++)
3299 mult->GetMultiplicityVtx(0)->SetBinContent(1, x+1, inputDistVector(x));
3300 mult->GetMultiplicityESD(0)->SetBinContent(1, x+1, measuredDistVector(x)*10000);
3301 for (Int_t y=0; y<kSize; y++)
3302 mult->GetCorrelation(0)->SetBinContent(1, x+1, y+1, matrix(x, y));
3305 //mult->DrawHistograms();
3307 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 0);
3308 //mult->SetCreateBigBin(kFALSE);
3309 mult->ApplyMinuitFit(0, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); //hist2->ProjectionY("mymchist"));
3311 //mult->ApplyBayesianMethod(0, kFALSE, AliMultiplicityCorrection::kTrVtx, 0, -1, 0, kFALSE);
3313 mult->DrawComparison("BlobelExample", 0, kFALSE, kTRUE, mult->GetMultiplicityVtx(0)->ProjectionY("mcmchist", 1, mult->GetMultiplicityVtx(0)->GetNbinsX()));
3320 TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
3321 fCurrentESD->Sumw2();
3323 // Open the input stream
3325 in.open("e735data.txt");
3332 //Printf("%f %f %f", x, y, ye);
3333 fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
3334 fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
3339 //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
3341 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
3343 TF1* func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])");
3344 func->SetParNames("scaling", "averagen", "k");
3345 func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
3346 func->SetParLimits(1, 0.001, 1000);
3347 func->SetParLimits(2, 0.001, 1000);
3348 func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
3350 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
3351 lognormal->SetParNames("scaling", "mean", "sigma");
3352 lognormal->SetParameters(1, 1, 1);
3353 lognormal->SetParLimits(0, 0, 10);
3354 lognormal->SetParLimits(1, 0, 100);
3355 lognormal->SetParLimits(2, 1e-3, 10);
3357 TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
3358 fCurrentESD->SetStats(kFALSE);
3359 fCurrentESD->SetMarkerStyle(0);
3360 fCurrentESD->SetLineColor(1);
3361 fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
3362 fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
3363 fCurrentESD->Draw("");
3364 fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
3365 fCurrentESD->Fit(func, "0", "", 0, 150);
3366 func->SetRange(0, 250);
3368 printf("chi2 = %f\n", func->GetChisquare());
3370 fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
3371 lognormal->SetLineColor(2);
3372 lognormal->SetLineStyle(2);
3373 lognormal->SetRange(0, 250);
3374 lognormal->Draw("SAME");
3378 canvas->SaveAs("E735Fit.eps");
3381 void DifferentSamples()
3386 const char* filesChi2[] = { "chi2_100k_1.root", "chi2_100k_2.root" };
3387 const char* filesBayesian[] = { "bayesian_100k_1.root", "bayesian_100k_2.root" };
3389 TCanvas* canvas = new TCanvas("DifferentSamples", "DifferentSamples", 1200, 600);
3390 canvas->Divide(2, 1);
3392 legend = new TLegend(0.15, 0.7, 0.65, 0.9);
3393 legend->SetFillColor(0);
3394 legend->SetTextSize(0.04);
3396 for (Int_t i=0; i<n; i++)
3398 AliMultiplicityCorrection* chi2 = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3399 TFile::Open(filesChi2[i]);
3400 chi2->LoadHistograms("Multiplicity");
3402 AliMultiplicityCorrection* bayesian = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3403 TFile::Open(filesBayesian[i]);
3404 bayesian->LoadHistograms("Multiplicity");
3406 chi2Hist = chi2->GetMultiplicityESDCorrected(etaRange);
3407 bayesianHist = bayesian->GetMultiplicityESDCorrected(etaRange);
3409 mc = chi2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
3411 // normalize and divide
3412 chi2Hist->Scale(1.0 / chi2Hist->Integral(2, displayRange+1) * mc->Integral(2, displayRange));
3413 bayesianHist->Scale(1.0 / bayesianHist->Integral(2, displayRange+1) * mc->Integral(2, displayRange));
3415 chi2Hist->Divide(mc, chi2Hist);
3416 bayesianHist->Divide(mc, bayesianHist);
3419 gPad->SetTopMargin(0.05);
3420 gPad->SetRightMargin(0.05);
3421 //gPad->SetLeftMargin(0.12);
3425 chi2Hist->GetXaxis()->SetRangeUser(0, displayRange);
3426 chi2Hist->GetYaxis()->SetTitleOffset(1.3);
3427 chi2Hist->SetStats(0);
3428 chi2Hist->SetTitle(Form(";%s;MC / unfolded", GetMultLabel()));
3429 chi2Hist->GetYaxis()->SetRangeUser(0.2, 1.8);
3430 chi2Hist->Draw("HIST");
3432 for (Int_t x=1; x<=bayesianHist->GetNbinsX(); x++)
3433 bayesianHist->SetBinError(x, 1e-6);
3435 bayesianHist->SetLineColor(2);
3436 bayesianHist->SetMarkerColor(2);
3437 bayesianHist->SetMarkerStyle(5);
3438 bayesianHist->Draw("HIST E SAME");
3442 legend->AddEntry(chi2Hist, "#chi^{2}-minimization", "L");
3443 legend->AddEntry(bayesianHist, "Bayesian unfolding", "LP");
3448 canvas->SaveAs("DifferentSamples.eps");
3455 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root");
3456 hist2d = mult->GetMultiplicityMC(etaRange, AliMultiplicityCorrection::kINEL);
3457 mult1 = hist2d->ProjectionY("mult1", 1, hist2d->GetNbinsX());
3459 conv = (TH1*) mult1->Clone("conv");
3462 mult1->Scale(1.0 / mult1->Integral());
3464 for (Int_t i=1; i<=mult1->GetNbinsX(); i++)
3465 for (Int_t j=1; j<=mult1->GetNbinsX(); j++)
3466 conv->Fill(mult1->GetBinCenter(i)+mult1->GetBinCenter(j), mult1->GetBinContent(i) * mult1->GetBinContent(j));
3468 conv->Scale(1.0 / conv->Integral());
3470 c = new TCanvas("c", "c", 800, 500);
3472 gPad->SetTopMargin(0.05);
3473 gPad->SetRightMargin(0.05);
3474 mult1->SetTitle(Form(";%s;Probability", GetMultLabel()));
3479 mult1->GetYaxis()->SetRangeUser(1e-7, 2 * mult1->GetMaximum());
3480 mult1->GetXaxis()->SetRangeUser(0, displayRange);
3481 mult1->GetXaxis()->SetTitleOffset(1.15);
3482 conv->SetLineColor(2);
3483 conv->SetMarkerColor(2);
3484 conv->SetMarkerStyle(5);
3485 conv->DrawCopy("SAME P");
3487 conv->Scale(0.00058);
3488 conv->DrawCopy("SAME P");
3490 legend = new TLegend(0.73, 0.73, 0.93, 0.93);
3491 legend->SetFillColor(0);
3492 legend->SetTextSize(0.04);
3493 legend->AddEntry(mult1, "1 collision");
3494 legend->AddEntry(conv, "2 collisions", "P");
3497 c->SaveAs("pileup.eps");
3500 conv->Divide(mult1);
3504 void TestErrorDetermination(Int_t nRandomizations)
3506 TF1* func = new TF1("nbd", "exp(log([0]) + TMath::LnGamma([2]+x) - TMath::LnGamma([2]) - TMath::LnGamma(x+1) + log([1] / ([1]+[2])) * x + log(1.0 + [1]/[2]) * -[2])", 0, 100);
3507 func->SetParNames("scaling", "averagen", "k");
3508 func->SetParameters(1, 15, 2);
3510 TF1* func2 = new TF1("nbd2", "exp(log([0]) + TMath::LnGamma([2]+x) - TMath::LnGamma([2]) - TMath::LnGamma(x+1) + log([1] / ([1]+[2])) * x + log(1.0 + [1]/[2]) * -[2])", 0, 100);
3511 func2->SetParNames("scaling", "averagen", "k");
3512 func2->SetParLimits(0, 0.5, 2);
3513 func2->SetParLimits(1, 1, 50);
3514 func2->SetParLimits(2, 1, 10);
3515 func2->SetParameters(1, 15, 2);
3516 //func2->FixParameter(0, 1);
3518 //new TCanvas; func->Draw("L");
3520 hist1 = new TH1F("hist1", "", 100, 0.5, 100.5);
3521 hist2 = new TH1F("hist2", "", 100, 0.5, 100.5);
3525 params[0] = new TH1F("param_0", Form("param_%d", 0), 100, 0.95, 1.05);
3526 params[1] = new TH1F("param_1", Form("param_%d", 1), 100, 14, 16);
3527 params[2] = new TH1F("param_2", Form("param_%d", 2), 100, 1.8, 2.2);
3529 const Int_t nTries = 1000;
3530 for (Int_t i=0; i<nTries; i++)
3534 if (nRandomizations == 1)
3536 hist1->FillRandom("nbd", 10000);
3538 else if (nRandomizations == 2)
3541 hist2->FillRandom("nbd", 10000);
3542 hist1->FillRandom(hist2, 10000);
3544 else if (nRandomizations == 3)
3547 hist1->FillRandom("nbd", 10000);
3548 hist2->FillRandom(hist1, 10000);
3550 hist1->FillRandom(hist2, 10000);
3555 //new TCanvas; hist1->Draw();
3557 hist1->Scale(1.0 / hist1->Integral());
3558 hist1->Fit(func2, "NQ");
3559 hist1->Fit(func2, "NQ");
3560 for (Int_t j=0; j<3; j++)
3561 params[j]->Fill(func2->GetParameter(j));
3564 for (Int_t j=0; j<3; j++)
3566 new TCanvas; params[j]->Draw();
3567 params[j]->Fit("gaus");
3568 Printf("sigma of param %d if %f", j, ((TF1*) params[j]->FindObject("gaus"))->GetParameter(2));
3572 void DrawRawDistributions(const char* fileName = "multiplicityESD.root")
3576 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileName);
3578 c = new TCanvas("c", "c", 600, 600);
3580 dummy = new TH2F("dummy", ";measured multiplicity", 100, -0.5, 149.5, 100, 0.5, 4e4);
3587 Int_t colors[] = { 1, 2, 4 };
3589 for (Int_t i=2; i>=0; i--)
3591 hist = mult->GetMultiplicityESD(i)->ProjectionY();
3593 hist->SetLineColor(colors[i]);
3594 hist->DrawCopy("SAME");