]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG0/multiplicity/plots.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
skeleton for documentation
[u/mrichter/AliRoot.git] / PWG0 / multiplicity / plots.C
CommitLineData
f3eb27f6 1/* $Id$ */
2
3//
4// plots for the note about multplicity measurements
5//
6
7#if !defined(__CINT__) || defined(__MAKECINT__)
8
9#include <TCanvas.h>
10#include <TPad.h>
11#include <TH1F.h>
12#include <TH2F.h>
13#include <TH3F.h>
14#include <TLine.h>
15#include <TF1.h>
16#include <TSystem.h>
17#include <TFile.h>
18#include <TLegend.h>
19#include <TStopwatch.h>
20#include <TROOT.h>
21#include <TGraph.h>
22#include <TMath.h>
23#include <TPaveText.h>
24#include <TImage.h>
25#include <TLatex.h>
26
27#include "AliMultiplicityCorrection.h"
28#include "AliCorrection.h"
29#include "AliCorrectionMatrix3D.h"
30
31#endif
32
51f6de65 33const char* correctionFile = "multiplicity.root";
f3eb27f6 34const char* measuredFile = "multiplicityMC_1M_3.root";
51f6de65 35Int_t etaRange = 2;
f3eb27f6 36Int_t displayRange = 200; // axis range
37Int_t ratioRange = 151; // range to calculate difference
38Int_t longDisplayRange = 200;
39
40const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
41const char* measuredFileTPC = "multiplicityMC_TPC_0.6M.root";
42Int_t etaRangeTPC = 1;
43
5a6310fe 44void loadlibs()
45{
46 gSystem->Load("libANALYSIS");
47 gSystem->Load("libPWG0base");
48}
49
f3eb27f6 50void SetTPC()
51{
52 correctionFile = correctionFileTPC;
53 measuredFile = measuredFileTPC;
54 etaRange = etaRangeTPC;
55 displayRange = 100;
56 ratioRange = 76;
57 longDisplayRange = 100;
58}
59
60void Smooth(TH1* hist, Int_t windowWidth = 20)
61{
62 TH1* clone = (TH1*) hist->Clone("clone");
63 for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin)
64 {
65 Int_t min = TMath::Max(2, bin-windowWidth);
66 Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
67 Float_t average = clone->Integral(min, max) / (max - min + 1);
68
69 hist->SetBinContent(bin, average);
70 hist->SetBinError(bin, 0);
71 }
72
73 delete clone;
74}
75
76void responseMatrixPlot()
77{
78 gSystem->Load("libPWG0base");
79
80 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
81
82 TFile::Open(correctionFile);
83 mult->LoadHistograms("Multiplicity");
84
85 TH1* hist = mult->GetCorrelation(etaRange)->Project3D("zy");
86 hist->SetStats(kFALSE);
87
88 hist->SetTitle(";true multiplicity;measured multiplicity;Entries");
89 hist->GetXaxis()->SetRangeUser(0, longDisplayRange);
90 hist->GetYaxis()->SetRangeUser(0, longDisplayRange);
91
92 TCanvas* canvas = new TCanvas("c1", "c1", 800, 600);
93 canvas->SetRightMargin(0.15);
94 canvas->SetTopMargin(0.05);
95
96 gPad->SetLogz();
97 hist->Draw("COLZ");
98
99 canvas->SaveAs("responsematrix.eps");
100}
101
102TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
103{
104 // normalize unfolded result to mc hist
105 result->Scale(1.0 / result->Integral(2, 200));
106 result->Scale(mcHist->Integral(2, 200));
107
108 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
109 canvas->Range(0, 0, 1, 1);
110
111 TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
112 pad1->Draw();
113
114 TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
115 pad2->Draw();
116
117 pad1->SetRightMargin(0.05);
118 pad2->SetRightMargin(0.05);
119
120 // no border between them
121 pad1->SetBottomMargin(0);
122 pad2->SetTopMargin(0);
123
124 pad1->cd();
125
126 mcHist->GetXaxis()->SetLabelSize(0.06);
127 mcHist->GetYaxis()->SetLabelSize(0.06);
128 mcHist->GetXaxis()->SetTitleSize(0.06);
129 mcHist->GetYaxis()->SetTitleSize(0.06);
130 mcHist->GetYaxis()->SetTitleOffset(0.6);
131
132 mcHist->GetXaxis()->SetRangeUser(0, displayRange);
133
134 mcHist->SetTitle(";true multiplicity;Entries");
135 mcHist->SetStats(kFALSE);
136
137 mcHist->DrawCopy("HIST E");
138 gPad->SetLogy();
139
140 result->SetLineColor(2);
141 result->DrawCopy("SAME HISTE");
142
143 TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
144 legend->AddEntry(mcHist, "true distribution");
145 legend->AddEntry(result, "unfolded distribution");
146 legend->SetFillColor(0);
147 legend->Draw();
148
149 pad2->cd();
150 pad2->SetBottomMargin(0.15);
151
152 // calculate ratio
153 mcHist->Sumw2();
154 TH1* ratio = (TH1*) mcHist->Clone("ratio");
155 result->Sumw2();
156 ratio->Divide(ratio, result, 1, 1, "");
157 ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
158 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
159
160 ratio->DrawCopy();
161
162 // get average of ratio
163 Float_t sum = 0;
164 for (Int_t i=2; i<=ratioRange; ++i)
165 {
166 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
167 }
168 sum /= ratioRange-1;
169
170 printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
171
172 TLine* line = new TLine(0, 1, displayRange, 1);
173 line->SetLineWidth(2);
174 line->Draw();
175
176 line = new TLine(0, 1.1, displayRange, 1.1);
177 line->SetLineWidth(2);
178 line->SetLineStyle(2);
179 line->Draw();
180 line = new TLine(0, 0.9, displayRange, 0.9);
181 line->SetLineWidth(2);
182 line->SetLineStyle(2);
183 line->Draw();
184
185 canvas->Modified();
186
187 canvas->SaveAs(epsName);
188
189 return canvas;
190}
191
192TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
193{
194 // draws the 3 plots in the upper plot
195 // draws the ratio between result1 and result2 in the lower plot
196
197 // normalize unfolded result to mc hist
198 result1->Scale(1.0 / result1->Integral(2, 200));
199 result1->Scale(mcHist->Integral(2, 200));
200 result2->Scale(1.0 / result2->Integral(2, 200));
201 result2->Scale(mcHist->Integral(2, 200));
202
203 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
204 canvas->Range(0, 0, 1, 1);
205
206 TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
207 pad1->Draw();
208
209 TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
210 pad2->Draw();
211
212 pad1->SetRightMargin(0.05);
213 pad2->SetRightMargin(0.05);
214
215 // no border between them
216 pad1->SetBottomMargin(0);
217 pad2->SetTopMargin(0);
218
219 pad1->cd();
220
221 mcHist->GetXaxis()->SetLabelSize(0.06);
222 mcHist->GetYaxis()->SetLabelSize(0.06);
223 mcHist->GetXaxis()->SetTitleSize(0.06);
224 mcHist->GetYaxis()->SetTitleSize(0.06);
225 mcHist->GetYaxis()->SetTitleOffset(0.6);
226
227 mcHist->GetXaxis()->SetRangeUser(0, displayRange);
228
229 mcHist->SetTitle(";true multiplicity;Entries");
230 mcHist->SetStats(kFALSE);
231
232 mcHist->DrawCopy("HIST E");
233 gPad->SetLogy();
234
235 result1->SetLineColor(2);
236 result1->DrawCopy("SAME HISTE");
237
238 result2->SetLineColor(4);
239 result2->DrawCopy("SAME HISTE");
240
241 TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
242 legend->AddEntry(mcHist, "true distribution");
243 legend->AddEntry(result1, "unfolded distribution (syst)");
244 legend->AddEntry(result2, "unfolded distribution (normal)");
245 legend->SetFillColor(0);
246 legend->Draw();
247
248 pad2->cd();
249 pad2->SetBottomMargin(0.15);
250
251 result1->GetXaxis()->SetLabelSize(0.06);
252 result1->GetYaxis()->SetLabelSize(0.06);
253 result1->GetXaxis()->SetTitleSize(0.06);
254 result1->GetYaxis()->SetTitleSize(0.06);
255 result1->GetYaxis()->SetTitleOffset(0.6);
256
257 result1->GetXaxis()->SetRangeUser(0, displayRange);
258
259 result1->SetTitle(";true multiplicity;Entries");
260 result1->SetStats(kFALSE);
261
262 // calculate ratio
263 result1->Sumw2();
264 TH1* ratio = (TH1*) result1->Clone("ratio");
265 result2->Sumw2();
266 ratio->Divide(ratio, result2, 1, 1, "");
267 ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
268 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
269
270 ratio->DrawCopy();
271
272 // get average of ratio
273 Float_t sum = 0;
274 for (Int_t i=2; i<=ratioRange; ++i)
275 {
276 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
277 }
278 sum /= ratioRange-1;
279
280 printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
281
282 TLine* line = new TLine(0, 1, displayRange, 1);
283 line->SetLineWidth(2);
284 line->Draw();
285
286 line = new TLine(0, 1.1, displayRange, 1.1);
287 line->SetLineWidth(2);
288 line->SetLineStyle(2);
289 line->Draw();
290 line = new TLine(0, 0.9, displayRange, 0.9);
291 line->SetLineWidth(2);
292 line->SetLineStyle(2);
293 line->Draw();
294
295 canvas->Modified();
296
297 canvas->SaveAs(epsName);
298
299 return canvas;
300}
301
302TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE)
303{
304 // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
305 // a systematic effect
306
307 // normalize results
308 result->Scale(1.0 / result->Integral(2, 200));
309
310 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
311 canvas->SetTopMargin(0.05);
312 canvas->SetRightMargin(0.05);
313
314 result->GetXaxis()->SetRangeUser(0, displayRange);
315 result->GetYaxis()->SetRangeUser(0.55, 1.45);
316 result->SetStats(kFALSE);
317
318 // to get the axis how we want it
319 TH1* dummy = (TH1*) result->Clone("dummy");
320 dummy->Reset();
321 dummy->SetTitle(";true multiplicity;Ratio");
322 dummy->DrawCopy();
323 delete dummy;
324
325 Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
326
327 TLegend* legend = new TLegend(0.2, 0.75, 0.35, 0.95);
328 legend->SetFillColor(0);
329
330 for (Int_t n=0; n<nResultSyst; ++n)
331 {
332 resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
333
334 // calculate ratio
335 TH1* ratio = (TH1*) result->Clone("ratio");
336 ratio->Divide(ratio, resultSyst[n], 1, 1, "");
337 ratio->GetXaxis()->SetRangeUser(1, displayRange);
338
339 if (firstMarker)
340 ratio->SetMarkerStyle(5);
341
342 ratio->SetLineColor(colors[n / 2]);
343 if ((n % 2))
344 ratio->SetLineStyle(2);
345
346 TString drawStr("SAME HIST");
347 if (n == 0 && firstMarker)
348 drawStr = "SAME P";
349 if (errors)
350 drawStr += " E";
351
352 ratio->DrawCopy(drawStr);
353
354 if (legendStrings && legendStrings[n])
355 legend->AddEntry(ratio, legendStrings[n]);
356
357 // get average of ratio
358 Float_t sum = 0;
359 for (Int_t i=2; i<=ratioRange; ++i)
360 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
361 sum /= ratioRange-1;
362
363 printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
364 }
365
366 if (legendStrings)
367 legend->Draw();
368
369 TLine* line = new TLine(0, 1, displayRange, 1);
370 line->SetLineWidth(2);
371 line->Draw();
372
373 line = new TLine(0, 1.1, displayRange, 1.1);
374 line->SetLineWidth(2);
375 line->SetLineStyle(2);
376 line->Draw();
377 line = new TLine(0, 0.9, displayRange, 0.9);
378 line->SetLineWidth(2);
379 line->SetLineStyle(2);
380 line->Draw();
381
382 canvas->SaveAs(epsName);
383 canvas->SaveAs(Form("%s.gif", epsName.Data()));
384
385 return canvas;
386}
387
388TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE)
389{
390 // draws the ratios of each mc to the corresponding result
391
392 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
393 canvas->SetRightMargin(0.05);
394 canvas->SetTopMargin(0.05);
395
396 for (Int_t n=0; n<nResultSyst; ++n)
397 {
398 // normalize
399 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
400 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
401
402 result[n]->GetXaxis()->SetRangeUser(0, displayRange);
403 result[n]->SetStats(kFALSE);
404
405 // calculate ratio
406 TH1* ratio = (TH1*) result[n]->Clone("ratio");
407 ratio->Divide(mc[n], ratio, 1, 1, "B");
408
409 // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis...
410 ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0);
411
412 if (smooth)
413 Smooth(ratio);
414
415 ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : "")));
416 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
417
418 if (dashed)
419 {
420 ratio->SetLineColor((n/2)+1);
421 ratio->SetLineStyle((n%2)+1);
422 }
423 else
424 ratio->SetLineColor(n+1);
425
426 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
427
428 // get average of ratio
429 Float_t sum = 0;
430 for (Int_t i=2; i<=ratioRange; ++i)
431 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
432 sum /= ratioRange-1;
433
434 printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
435 }
436
437 TLine* line = new TLine(0, 1, displayRange, 1);
438 line->SetLineWidth(2);
439 line->Draw();
440
441 line = new TLine(0, 1.1, displayRange, 1.1);
442 line->SetLineWidth(2);
443 line->SetLineStyle(2);
444 line->Draw();
445 line = new TLine(0, 0.9, displayRange, 0.9);
446 line->SetLineWidth(2);
447 line->SetLineStyle(2);
448 line->Draw();
449
450 canvas->Modified();
451
452 canvas->SaveAs(epsName);
453 canvas->SaveAs(Form("%s.gif", epsName.Data()));
454
455 return canvas;
456}
457
458TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
459{
460 // draws the ratios of each mc to the corresponding result
461 // deducts from each ratio the ratio of mcBase / resultBase
462
463 // normalize
464 resultBase->Scale(1.0 / resultBase->Integral(2, 200));
465 mcBase->Scale(1.0 / mcBase->Integral(2, 200));
466
467 // calculate ratio
468 TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
469 ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
470
471 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
472 canvas->SetRightMargin(0.05);
473 canvas->SetTopMargin(0.05);
474
475 for (Int_t n=0; n<nResultSyst; ++n)
476 {
477 // normalize
478 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
479 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
480
481 result[n]->GetXaxis()->SetRangeUser(0, displayRange);
482 result[n]->SetStats(kFALSE);
483
484 // calculate ratio
485 TH1* ratio = (TH1*) result[n]->Clone("ratio");
486 ratio->Divide(mc[n], ratio, 1, 1, "B");
487 ratio->Add(ratioBase, -1);
488
489 ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
490 ratio->GetYaxis()->SetRangeUser(-1, 1);
491 ratio->SetLineColor(n+1);
492 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
493
494 // get average of ratio
495 Float_t sum = 0;
496 for (Int_t i=2; i<=ratioRange; ++i)
497 sum += TMath::Abs(ratio->GetBinContent(i));
498 sum /= ratioRange-1;
499
500 printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, ratioRange, sum);
501 }
502
503 TLine* line = new TLine(0, 0, displayRange, 0);
504 line->SetLineWidth(2);
505 line->Draw();
506
507 line = new TLine(0, 0.1, displayRange, 0.1);
508 line->SetLineWidth(2);
509 line->SetLineStyle(2);
510 line->Draw();
511 line = new TLine(0, -0.1, displayRange, -0.1);
512 line->SetLineWidth(2);
513 line->SetLineStyle(2);
514 line->Draw();
515
516 canvas->Modified();
517
518 canvas->SaveAs(epsName);
519 canvas->SaveAs(Form("%s.gif", epsName.Data()));
520
521 return canvas;
522}
523
524TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
525{
526 // draws the ratios of each mc to the corresponding result
527 // deducts from each ratio the ratio of mcBase / resultBase
528 // smoothens the ratios by a sliding window
529
530 // normalize
531 resultBase->Scale(1.0 / resultBase->Integral(2, 200));
532 mcBase->Scale(1.0 / mcBase->Integral(2, 200));
533
534 // calculate ratio
535 TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
536 ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
537
538 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
539 canvas->SetRightMargin(0.05);
540 canvas->SetTopMargin(0.05);
541
542 for (Int_t n=0; n<nResultSyst; ++n)
543 {
544 // normalize
545 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
546 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
547
548 result[n]->GetXaxis()->SetRangeUser(0, displayRange);
549 result[n]->SetStats(kFALSE);
550
551 // calculate ratio
552 TH1* ratio = (TH1*) result[n]->Clone("ratio");
553 ratio->Divide(mc[n], ratio, 1, 1, "B");
554 ratio->Add(ratioBase, -1);
555
556 //new TCanvas; ratio->DrawCopy();
557 // clear 0 bin
558 ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0);
559
560 Smooth(ratio);
561
562 //ratio->SetLineColor(1); ratio->DrawCopy("SAME");
563
564 canvas->cd();
565 ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)");
566 ratio->GetYaxis()->SetRangeUser(-0.3, 0.3);
567 ratio->SetLineColor((n / 2)+1);
568 ratio->SetLineStyle((n % 2)+1);
569 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
570
571 // get average of ratio
572 Float_t sum = 0;
573 for (Int_t i=2; i<=150; ++i)
574 sum += TMath::Abs(ratio->GetBinContent(i));
575 sum /= 149;
576
577 printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
578 }
579
580 TLine* line = new TLine(0, 0, displayRange, 0);
581 line->SetLineWidth(2);
582 line->Draw();
583
584 line = new TLine(0, 0.1, displayRange, 0.1);
585 line->SetLineWidth(2);
586 line->SetLineStyle(2);
587 line->Draw();
588 line = new TLine(0, -0.1, displayRange, -0.1);
589 line->SetLineWidth(2);
590 line->SetLineStyle(2);
591 line->Draw();
592
593 canvas->Modified();
594
595 canvas->SaveAs(epsName);
596 canvas->SaveAs(Form("%s.gif", epsName.Data()));
597
598 return canvas;
599}
600
601void DrawResiduals(TH1* measured, TH1* unfoldedFolded, const char* epsName)
602{
603 // normalize
604 unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, 200));
605 unfoldedFolded->Scale(measured->Integral(2, 200));
606
607 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
608 canvas->Range(0, 0, 1, 1);
609
610 TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
611 pad1->Draw();
612 pad1->SetGridx();
613 pad1->SetGridy();
614
615 TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
616 pad2->Draw();
617 pad2->SetGridx();
618 pad2->SetGridy();
619
620 TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
621 pad3->SetGridx();
622 pad3->SetGridy();
623 pad3->SetRightMargin(0.05);
624 pad3->SetTopMargin(0.05);
625 pad3->Draw();
626
627 pad1->SetRightMargin(0.05);
628 pad2->SetRightMargin(0.05);
629
630 // no border between them
631 pad1->SetBottomMargin(0);
632 pad2->SetTopMargin(0);
633
634 pad1->cd();
635
636 measured->GetXaxis()->SetLabelSize(0.06);
637 measured->GetYaxis()->SetLabelSize(0.06);
638 measured->GetXaxis()->SetTitleSize(0.06);
639 measured->GetYaxis()->SetTitleSize(0.06);
640 measured->GetYaxis()->SetTitleOffset(0.6);
641
642 measured->GetXaxis()->SetRangeUser(0, 150);
643
644 measured->SetTitle(";measured multiplicity;Entries");
645 measured->SetStats(kFALSE);
646
647 measured->DrawCopy("HIST");
648 gPad->SetLogy();
649
650 unfoldedFolded->SetMarkerStyle(5);
651 unfoldedFolded->SetMarkerColor(2);
652 unfoldedFolded->SetLineColor(0);
653 unfoldedFolded->DrawCopy("SAME P");
654
655 TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
656 legend->AddEntry(measured, "measured distribution");
657 legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution");
658 legend->SetFillColor(0);
659 legend->Draw();
660
661 pad2->cd();
662 pad2->SetBottomMargin(0.15);
663
664 // calculate ratio
665 measured->Sumw2();
666 TH1* residual = (TH1*) measured->Clone("residual");
667 unfoldedFolded->Sumw2();
668
669 residual->Add(unfoldedFolded, -1);
670
671 // projection
672 TH1* residualHist = new TH1F("residualHist", ";", 15, -3, 3);
673
674 for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
675 {
676 if (measured->GetBinError(i) > 0)
677 {
678 residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
679 residual->SetBinError(i, 1);
680
681 residualHist->Fill(residual->GetBinContent(i));
682 }
683 else
684 {
685 residual->SetBinContent(i, 0);
686 residual->SetBinError(i, 0);
687 }
688 }
689
690 residual->GetYaxis()->SetTitle("Residuals 1/e (M - R #otimes U)");
691 residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
692 residual->DrawCopy();
693
694 TLine* line = new TLine(-0.5, 0, 150.5, 0);
695 line->SetLineWidth(2);
696 line->Draw();
697
698 pad3->cd();
699 residualHist->SetStats(kFALSE);
700 residualHist->GetXaxis()->SetLabelSize(0.08);
701 residualHist->Fit("gaus");
702 residualHist->Draw();
703
704 canvas->Modified();
705 canvas->SaveAs(canvas->GetName());
706
707 //const char* epsName2 = "proj.eps";
708 //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
709 //canvas->SetGridx();
710 //canvas->SetGridy();
711
712 //canvas->SaveAs(canvas->GetName());
713}
714
715void bayesianExample()
716{
717 TStopwatch watch;
718 watch.Start();
719
720 gSystem->Load("libPWG0base");
721
722 TFile::Open(correctionFile);
723 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
724 mult->LoadHistograms("Multiplicity");
725
726 TFile::Open(measuredFile);
727 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
728 mult2->LoadHistograms("Multiplicity");
729
730 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
731
732 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
733
734 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
735 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
736
737 //mult->DrawComparison("bayesianExample", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
738 DrawResultRatio(mcHist, result, "bayesianExample.eps");
739
740 //Printf("KolmogorovTest says PROB = %f", mcHist->KolmogorovTest(result, "D"));
741 //Printf("Chi2Test says PROB = %f", mcHist->Chi2Test(result));
742
743 // draw residual plot
744
745 // TODO take out efficiency correction if other than AliMultiplicityCorrection::kTrVtx
746 TH2* convoluted = mult->CalculateMultiplicityESD(result, etaRange);
747 TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
748
749 TH1* measured = mult2->GetMultiplicityESD(etaRange)->ProjectionY("measured");
750
751 DrawResiduals(measured, convolutedProj, "bayesianResiduals.eps");
752
753 watch.Stop();
754 watch.Print();
755}
756
757void chi2FluctuationResult()
758{
759 gSystem->Load("libPWG0base");
760
761 TFile::Open(correctionFile);
762 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
763 mult->LoadHistograms("Multiplicity");
764
765 TFile::Open(measuredFile);
766 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
767 mult2->LoadHistograms("Multiplicity");
768
769 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
770 mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
771 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
772
773 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
774 //TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
775
776 mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
777
778 TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_3");
779 canvas->SaveAs("chi2FluctuationResult.eps");
780}
781
782void chi2Example()
783{
784 TStopwatch watch;
785 watch.Start();
786
787 gSystem->Load("libPWG0base");
788
789 TFile::Open(correctionFile);
790 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
791 mult->LoadHistograms("Multiplicity");
792
793 TFile::Open(measuredFile);
794 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
795 mult2->LoadHistograms("Multiplicity");
796
797 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
798 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
799 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
800
801 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
802 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
803
804 DrawResultRatio(mcHist, result, "chi2Example.eps");
805
806 watch.Stop();
807 watch.Print();
808}
809
810void chi2ExampleTPC()
811{
812 TStopwatch watch;
813 watch.Start();
814
815 gSystem->Load("libPWG0base");
816
817 TFile::Open(correctionFileTPC);
818 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
819 mult->LoadHistograms("Multiplicity");
820
821 TFile::Open(measuredFileTPC);
822 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
823 mult2->LoadHistograms("Multiplicity");
824
825 mult->SetMultiplicityESD(etaRangeTPC, mult2->GetMultiplicityESD(etaRangeTPC));
826 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
827 mult->ApplyMinuitFit(etaRangeTPC, kFALSE, AliMultiplicityCorrection::kTrVtx);
828
829 TH1* mcHist = mult2->GetMultiplicityVtx(etaRangeTPC)->ProjectionY("mymc");
830 TH1* result = mult->GetMultiplicityESDCorrected(etaRangeTPC);
831
832 DrawResultRatio(mcHist, result, "chi2ExampleTPC.eps");
833
834 watch.Stop();
835 watch.Print();
836}
837
838void bayesianNBD()
839{
840 gSystem->Load("libPWG0base");
841 TFile::Open("multiplicityMC_3M.root");
842
843 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
844 mult->LoadHistograms("Multiplicity");
845
846 TFile::Open("multiplicityMC_3M_NBD.root");
847 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
848 mult2->LoadHistograms("Multiplicity");
849
850 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
851 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1.0, 100);
852
853 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
854
855 mcHist->Sumw2();
856 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
857
858 //mult->DrawComparison("bayesianNBD", etaRange, kFALSE, kTRUE, mcHist);
859 DrawResultRatio(mcHist, result, "bayesianNBD.eps");
860}
861
862void minimizationNBD()
863{
864 gSystem->Load("libPWG0base");
865 TFile::Open("multiplicityMC_3M.root");
866
867 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
868 mult->LoadHistograms("Multiplicity");
869
870 TFile::Open("multiplicityMC_3M_NBD.root");
871 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
872 mult2->LoadHistograms("Multiplicity");
873
874 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
875 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
876 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
877
878 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
879
880 mcHist->Sumw2();
881 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
882
883 //mult->DrawComparison("minimizationNBD", etaRange, kFALSE, kTRUE, mcHist);
884 DrawResultRatio(mcHist, result, "minimizationNBD.eps");
885}
886
887void minimizationInfluenceAlpha()
888{
889 gSystem->Load("libPWG0base");
890
891 TFile::Open(measuredFile);
892 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
893 mult2->LoadHistograms("Multiplicity");
894
895 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
896 mcHist->Scale(1.0 / mcHist->Integral());
897 mcHist->GetXaxis()->SetRangeUser(0, 200);
898 mcHist->SetStats(kFALSE);
899 mcHist->SetTitle(";true multiplicity;P_{N}");
900
901 TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 300);
902 canvas->Divide(3, 1);
903
904 TFile::Open("eval-2M-1M/EvaluateChi2MethodDetail.root");
905
906 TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_100.000000");
907 TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_03_2_100000.000000");
908 TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_06_2_100000000.000000");
909
910 mcHist->Rebin(2); mcHist->Scale(0.5);
911 hist1->Rebin(2); hist1->Scale(0.5);
912 hist2->Rebin(2); hist2->Scale(0.5);
913 hist3->Rebin(2); hist3->Scale(0.5);
914
915 mcHist->GetXaxis()->SetRangeUser(0, 200);
916
917 canvas->cd(1);
918 gPad->SetLogy();
919 mcHist->Draw();
920 hist1->SetMarkerStyle(5);
921 hist1->SetMarkerColor(2);
922 hist1->Draw("SAME PE");
923
924 canvas->cd(2);
925 gPad->SetLogy();
926 mcHist->Draw();
927 hist2->SetMarkerStyle(5);
928 hist2->SetMarkerColor(2);
929 hist2->Draw("SAME PE");
930
931 canvas->cd(3);
932 gPad->SetLogy();
933 mcHist->Draw();
934 hist3->SetMarkerStyle(5);
935 hist3->SetMarkerColor(2);
936 hist3->Draw("SAME PE");
937
938 canvas->SaveAs("minimizationInfluenceAlpha.eps");
939}
940
941void NBDFit()
942{
943 gSystem->Load("libPWG0base");
944
945 TFile::Open(correctionFile);
946 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
947 mult->LoadHistograms("Multiplicity");
948
949 TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
950 fCurrentESD->Sumw2();
951 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
952
953 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])");
954 func->SetParNames("scaling", "averagen", "k");
955 func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
956 func->SetParLimits(1, 0.001, 1000);
957 func->SetParLimits(2, 0.001, 1000);
958 func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
959
960 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
961 lognormal->SetParNames("scaling", "mean", "sigma");
962 lognormal->SetParameters(1, 1, 1);
963 lognormal->SetParLimits(0, 0, 10);
964 lognormal->SetParLimits(1, 0, 100);
965 lognormal->SetParLimits(2, 1e-3, 10);
966
967 TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
968 fCurrentESD->SetStats(kFALSE);
969 fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
970 fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
971 fCurrentESD->Draw("HIST");
972 fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
973 fCurrentESD->Fit(func, "W0", "", 0, 50);
974 func->SetRange(0, 100);
975 func->Draw("SAME");
976 printf("chi2 = %f\n", func->GetChisquare());
977
978 fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
979 lognormal->SetLineColor(2);
980 lognormal->SetLineStyle(2);
981 lognormal->SetRange(0, 100);
982 lognormal->Draw("SAME");
983
984 canvas->SaveAs("NBDFit.eps");
985}
986
987void DifferentSamples()
988{
989 // data generated by runMultiplicitySelector.C DifferentSamples
990
991 const char* name = "DifferentSamples";
992
993 TFile* file = TFile::Open(Form("%s.root", name));
994
995 TCanvas* canvas = new TCanvas(name, name, 800, 600);
996 canvas->Divide(2, 2);
997
998 for (Int_t i=0; i<4; ++i)
999 {
1000 canvas->cd(i+1);
1001 gPad->SetTopMargin(0.05);
1002 gPad->SetRightMargin(0.05);
1003 TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", i));
1004 TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", i));
1005 TH1* mc = (TH1*) file->Get(Form("mc_%d", i));
1006 mc->Sumw2();
1007
1008 chi2Result->Divide(chi2Result, mc, 1, 1, "");
1009 bayesResult->Divide(bayesResult, mc, 1, 1, "");
1010
1011 chi2Result->SetTitle(";true multiplicity;unfolded measured/MC");
1012 chi2Result->GetXaxis()->SetRangeUser(0, 150);
1013 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
1014 chi2Result->GetYaxis()->SetTitleOffset(1.2);
1015 chi2Result->SetLineColor(1);
1016 chi2Result->SetStats(kFALSE);
1017
1018 bayesResult->SetStats(kFALSE);
1019 bayesResult->SetLineColor(2);
1020
1021 chi2Result->DrawCopy("HIST");
1022 bayesResult->DrawCopy("SAME HIST");
1023
1024 TLine* line = new TLine(0, 1, 150, 1);
1025 line->Draw();
1026 }
1027
1028 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1029}
1030
1031void StartingConditions()
1032{
1033 // data generated by runMultiplicitySelector.C StartingConditions
1034
1035 const char* name = "StartingConditions";
1036
1037 TFile* file = TFile::Open(Form("%s.root", name));
1038
1039 TCanvas* canvas = new TCanvas(name, name, 800, 400);
1040 canvas->Divide(2, 1);
1041
1042 TH1* mc = (TH1*) file->Get("mc");
1043 mc->Sumw2();
1044 mc->Scale(1.0 / mc->Integral());
1045
1046 //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
1047
1048 TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.95);
1049 legend->SetFillColor(0);
1050
1051 const char* names[] = { "True", "Measured 1", "Measured 2", "Measured 3", "NBD", "Flat" };
1052
1053 for (Int_t i=0; i<6; ++i)
1054 {
1055 Int_t id = i;
1056 if (id > 2)
1057 id += 2;
1058
1059 TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
1060 TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
1061
1062 chi2Result->Divide(chi2Result, mc, 1, 1, "");
1063 bayesResult->Divide(bayesResult, mc, 1, 1, "");
1064
1065 chi2Result->SetTitle("a) #chi^{2} minimization;true multiplicity;unfolded / MC");
1066 chi2Result->GetXaxis()->SetRangeUser(0, 150);
1067 chi2Result->GetYaxis()->SetRangeUser(0.8, 1.2);
1068 chi2Result->GetYaxis()->SetTitleOffset(1.5);
1069 //chi2Result->SetMarkerStyle(marker[i]);
1070 chi2Result->SetLineColor(i+1);
1071 chi2Result->SetMarkerColor(i+1);
1072 chi2Result->SetStats(kFALSE);
1073
1074 bayesResult->SetTitle("b) Bayesian method;true multiplicity;unfolded / MC");
1075 bayesResult->GetXaxis()->SetRangeUser(0, 150);
1076 bayesResult->GetYaxis()->SetRangeUser(0.8, 1.2);
1077 bayesResult->GetYaxis()->SetTitleOffset(1.5);
1078 bayesResult->SetStats(kFALSE);
1079 //bayesResult->SetLineColor(2);
1080 bayesResult->SetLineColor(i+1);
1081
1082 canvas->cd(1);
1083 gPad->SetLeftMargin(0.12);
1084 chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1085
1086 canvas->cd(2);
1087 gPad->SetLeftMargin(0.12);
1088 bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1089
1090 //TLine* line = new TLine(0, 1, 150, 1);
1091 //line->Draw();
1092
1093 legend->AddEntry(chi2Result, names[i]);
1094 }
1095
1096 canvas->cd(1);
1097 legend->Draw();
1098 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1099}
1100
1101void StatisticsPlot()
1102{
1103 const char* name = "StatisticsPlot";
1104
1105 TFile* file = TFile::Open(Form("%s.root", name));
1106
1107 TCanvas* canvas = new TCanvas(name, name, 600, 400);
1108
1109 TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
1110 fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1111 fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1112 fitResultsChi2->Draw("AP");
1113
1114 TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1115 fitResultsChi2->Fit(f);
1116
1117 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1118
1119 TH1* mc[5];
1120 TH1* result[5];
1121
1122 const char* plotname = "chi2Result";
1123
1124 name = "StatisticsPlotRatios";
1125 canvas = new TCanvas(name, name, 600, 400);
1126
1127 for (Int_t i=0; i<5; ++i)
1128 {
1129 mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1130 result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1131
1132 result[i]->SetLineColor(i+1);
1133 result[i]->Draw(((i == 0) ? "" : "SAME"));
1134 }
1135}
1136
1137void SystematicLowEfficiency()
1138{
1139 gSystem->Load("libPWG0base");
1140
1141 TFile::Open(correctionFile);
1142 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1143 mult->LoadHistograms("Multiplicity");
1144
1145 // calculate result with systematic effect
1146 TFile::Open("multiplicityMC_100k_1_loweff.root");
1147 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1148 mult2->LoadHistograms("Multiplicity");
1149
1150 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1151 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1152 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1153
1154 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1155 TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
1156
1157 DrawResultRatio(mcHist, result1, "SystematicLowEfficiencyLow.eps");
1158
1159 // calculate normal result
1160 TFile::Open("multiplicityMC_100k_1.root");
1161 mult2->LoadHistograms("Multiplicity");
1162
1163 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1164 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1165
1166 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
1167
1168 DrawResultRatio(mcHist, result2, "SystematicLowEfficiencyOK.eps");
1169
1170 Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
1171}
1172
1173void SystematicMisalignment()
1174{
1175 gSystem->Load("libPWG0base");
1176
1177 TFile::Open(correctionFile);
1178 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1179 mult->LoadHistograms("Multiplicity");
1180
1181 TFile::Open("multiplicityMC_100k_fullmis.root");
1182 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1183 mult2->LoadHistograms("Multiplicity");
1184
1185 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1186 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1187 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1188
1189 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1190 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1191
1192 DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1193}
1194
1195void SystematicMisalignmentTPC()
1196{
1197 gSystem->Load("libPWG0base");
1198
1199 SetTPC();
1200
1201 TFile::Open(correctionFile);
1202 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1203 mult->LoadHistograms("Multiplicity");
1204
1205 TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
1206 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1207 mult2->LoadHistograms("Multiplicity");
1208
1209 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1210 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1211 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1212
1213 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1214 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1215
1216 DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
1217}
1218
1219void EfficiencySpecies()
1220{
5a6310fe 1221 loadlibs();
f3eb27f6 1222
1223 Int_t marker[] = {24, 25, 26};
1224 Int_t color[] = {1, 2, 4};
1225
1226 // SPD TPC
5a6310fe 1227 //const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
1228 const char* fileName[] = { "spd/multiplicity.root", "tpc/multiplicity.root" };
f3eb27f6 1229 Float_t etaRange[] = {0.49, 0.9};
1230 const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1231
1232 TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1233 canvas->Divide(2, 1);
1234
5a6310fe 1235 for (Int_t loop=0; loop<2; ++loop)
f3eb27f6 1236 {
1237 Printf("%s", fileName[loop]);
1238
1239 AliCorrection* correction[4];
1240
1241 canvas->cd(loop+1);
1242
1243 gPad->SetGridx();
1244 gPad->SetGridy();
1245 gPad->SetRightMargin(0.05);
1246 //gPad->SetTopMargin(0.05);
1247
1248 TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
1249 legend->SetFillColor(0);
1250 legend->SetEntrySeparation(0.2);
1251
1252 Float_t below = 0;
1253 Float_t total = 0;
1254
1255 TFile* file = TFile::Open(fileName[loop]);
1256 if (!file)
1257 {
1258 Printf("Could not open %s", fileName[loop]);
1259 return;
1260 }
1261
1262 Float_t sumGen = 0;
1263 Float_t sumMeas = 0;
1264
1265 for (Int_t i=0; i<3; ++i)
1266 {
1267 Printf("correction %d", i);
1268
1269 TString name; name.Form("correction_%d", i);
1270 correction[i] = new AliCorrection(name, name);
1271 correction[i]->LoadHistograms();
1272
1273 TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1274 TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1275
1276 // limit vtx axis
1277 gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
1278 meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
1279
1280 // 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)
1281 /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1282 for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1283 {
1284 gene->SetBinContent(x, 0, z, 0);
1285 gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1286 meas->SetBinContent(x, 0, z, 0);
1287 meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1288 }*/
1289
1290 // limit eta axis
1291 gene->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1292 meas->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1293
1294 TH1* genePt = gene->Project3D(Form("z_%d", i));
1295 TH1* measPt = meas->Project3D(Form("z_%d", i));
1296
1297 genePt->Sumw2();
1298 measPt->Sumw2();
1299
1300 sumGen += genePt->Integral();
1301 sumMeas += measPt->Integral();
1302
1303 TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1304 effPt->Reset();
1305 effPt->Divide(measPt, genePt, 1, 1, "B");
1306
1307 Int_t bin = 0;
1308 for (bin=20; bin>=1; bin--)
1309 {
1310 if (effPt->GetBinContent(bin) < 0.5)
1311 break;
1312 }
1313
1314 Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1315
1316 Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1317 Printf("%.4f of the particles are below that momentum", fraction);
1318
1319 below += genePt->Integral(1, bin);
1320 total += genePt->Integral();
1321
1322 effPt->SetLineColor(color[i]);
1323 effPt->SetMarkerColor(color[i]);
1324 effPt->SetMarkerStyle(marker[i]);
1325
1326 effPt->GetXaxis()->SetRangeUser(0.06, 1);
1327 effPt->GetYaxis()->SetRangeUser(0, 1);
1328
5a6310fe 1329 effPt->GetXaxis()->SetTitleOffset(1.1);
f3eb27f6 1330 effPt->GetYaxis()->SetTitleOffset(1.2);
1331
1332 effPt->SetStats(kFALSE);
1333 effPt->SetTitle(titles[loop]);
1334 effPt->GetYaxis()->SetTitle("Efficiency");
1335
1336 effPt->DrawCopy((i == 0) ? "" : "SAME");
1337
1338 legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")));
1339 }
1340
1341 Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1342
1343 Printf("%f measured, %f generated, effiency: %f", sumGen, sumMeas, sumMeas / sumGen);
1344
1345 legend->Draw();
1346 }
1347
1348 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1349}
1350
1351void ParticleSpeciesComparison1(Bool_t chi2 = kTRUE, const char* fileNameMC = "multiplicityMC_400k_syst_species.root", const char* fileNameESD = "multiplicityMC_100k_syst.root")
1352{
1353 gSystem->Load("libPWG0base");
1354
1355 TFile::Open(fileNameESD);
1356 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
1357 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", etaRange));
1358
1359 TH1* results[10];
1360
1361 // loop over cases (normal, enhanced/reduced ratios)
1362 Int_t nMax = 7;
1363 for (Int_t i = 0; i<nMax; ++i)
1364 {
1365 TString folder;
1366 folder.Form("Multiplicity_%d", i);
1367
1368 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
1369
1370 TFile::Open(fileNameMC);
1371 mult->LoadHistograms();
1372
1373 mult->SetMultiplicityESD(etaRange, hist);
1374
1375 if (chi2)
1376 {
1377 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1378 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1379 //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1380 }
1381 else
1382 {
1383 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1384 //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1385 }
1386
1387 //Float_t averageRatio = 0;
1388 //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1389
1390 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1391
1392 //Printf("Case %d. Average ratio is %f", i, averageRatio);
1393 }
1394
1395 DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[0], "ParticleSpeciesComparison1_1.eps");
1396
1397 TH1* mc = hist2->ProjectionY("mymchist2", -1, -1, "e");
1398
1399 for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1400 {
1401 results[0]->SetBinError(i, 0);
1402 mc->SetBinError(i, 0);
1403 }
1404
1405 const char* legendStrings[] = { "#pi^{#pm}", 0, "K^{#pm}", 0, "p,#bar{p}", 0 };
1406
1407 DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
1408
1409 //not valid: draw chi2 uncertainty on top!
1410 /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1411 TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1412 errorHist->SetLineColor(1);
1413 errorHist->SetLineWidth(2);
1414 TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1415 for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1416 {
1417 errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1418 errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1419 }
1420
1421 errorHist->DrawCopy("SAME");
1422 errorHist2->DrawCopy("SAME");*/
1423
1424 //canvas->SaveAs(canvas->GetName());
1425
1426 DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
1427
1428 //errorHist->DrawCopy("SAME");
1429 //errorHist2->DrawCopy("SAME");
1430
1431 //canvas2->SaveAs(canvas2->GetName());
1432}
1433
1434/*void ParticleSpeciesComparison2()
1435{
1436 gSystem->Load("libPWG0base");
1437
1438 const char* fileNameMC = "multiplicityMC_400k_syst.root";
1439 const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1440 Bool_t chi2 = 0;
1441
1442 TFile::Open(fileNameMC);
1443 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1444 mult->LoadHistograms();
1445
1446 TH1* mc[10];
1447 TH1* results[10];
1448
1449 // loop over cases (normal, enhanced/reduced ratios)
1450 Int_t nMax = 7;
1451 for (Int_t i = 0; i<nMax; ++i)
1452 {
1453 TString folder;
1454 folder.Form("Multiplicity_%d", i);
1455
1456 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1457
1458 TFile::Open(fileNameESD);
1459 mult2->LoadHistograms();
1460
1461 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1462
1463 if (chi2)
1464 {
1465 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1466 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1467 //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1468 }
1469 else
1470 {
1471 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1472 //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1473 }
1474
1475 //Float_t averageRatio = 0;
1476 //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1477
1478 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1479
1480 TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1481 mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1482
1483 //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1484 //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1485
1486 //Printf("Case %d. Average ratio is %f", i, averageRatio);
1487 }
1488
1489 DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1490}*/
1491
1492TH1* Invert(TH1* eff)
1493{
1494 // calculate corr = 1 / eff
1495
1496 TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1497 corr->Reset();
1498
1499 for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1500 {
1501 if (eff->GetBinContent(i) > 0)
1502 {
1503 corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1504 corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1505 }
1506 }
1507
1508 return corr;
1509}
1510
1511void TriggerVertexCorrection()
1512{
1513 //
1514 // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1515 //
1516
1517 gSystem->Load("libPWG0base");
1518
1519 TFile::Open(correctionFile);
1520 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1521 mult->LoadHistograms("Multiplicity");
1522
1523 TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1524 TH1* corrMB = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1525
1526 TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 600);
1527
1528 corrINEL->SetStats(kFALSE);
1529 corrINEL->GetXaxis()->SetRangeUser(0, 20);
1530 corrINEL->GetYaxis()->SetRangeUser(0.5, 2.5);
1531 corrINEL->SetTitle(";true multiplicity;correction factor");
1532 corrINEL->SetMarkerStyle(22);
1533 corrINEL->Draw("PE");
1534
1535 corrMB->SetStats(kFALSE);
1536 corrMB->SetLineColor(2);
1537 corrMB->SetMarkerStyle(25);
1538 corrMB->SetMarkerColor(2);
1539 corrMB->Draw("SAME PE");
1540
1541 TLegend* legend = new TLegend(0.3, 0.5, 0.85, 0.65);
1542 legend->SetFillColor(0);
1543 legend->AddEntry(corrINEL, "correction to inelastic sample");
1544 legend->AddEntry(corrMB, "correction to minimum bias sample");
1545
1546 legend->Draw();
1547
1548 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1549}
1550
1551void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
1552{
1553 gSystem->Load("libPWG0base");
1554
1555 TFile::Open(correctionFile);
1556 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1557 mult->LoadHistograms("Multiplicity");
1558
1559 TFile::Open(measuredFile);
1560 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1561 mult2->LoadHistograms("Multiplicity");
1562
1563 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1564
1565 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1566
1567 TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1568
1569 TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1570 TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
1571
1572 if (!mc)
1573 {
1574 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1575 DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
1576 }
1577
1578 TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
1579 canvas->SetGridx();
1580 canvas->SetGridy();
1581 canvas->SetRightMargin(0.05);
1582 canvas->SetTopMargin(0.05);
1583
1584 errorResponse->SetLineColor(1);
1585 errorResponse->GetXaxis()->SetRangeUser(0, 200);
1586 errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1587 errorResponse->SetStats(kFALSE);
1588 errorResponse->SetTitle(";true multiplicity;Uncertainty");
1589
1590 errorResponse->Draw();
1591
1592 errorMeasured->SetLineColor(2);
1593 errorMeasured->Draw("SAME");
1594
1595 errorBoth->SetLineColor(4);
1596 errorBoth->Draw("SAME");
1597
1598 Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
1599 Printf("Average errorMeasured: %f", errorMeasured->Integral(2, 150) / 149);
1600 Printf("Average errorBoth: %f", errorBoth->Integral(2, 150) / 149);
1601
1602 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1603
1604 TFile* file = new TFile(Form("%s.root", canvas->GetName()), "RECREATE");
1605 errorResponse->Write();
1606 errorMeasured->Write();
1607 errorBoth->Write();
1608 file->Close();
1609}
1610
1611void StatisticalUncertaintyCompare(const char* det = "SPD")
1612{
1613 TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1614 TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1615 TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1616 TH1* errorBoth = (TH1*) file1->Get("errorBoth");
1617
1618 TString str;
1619 str.Form("StatisticalUncertaintyCompare%s", det);
1620
1621 TCanvas* canvas = new TCanvas(str, str, 600, 400);
1622 canvas->SetGridx();
1623 canvas->SetGridy();
1624 canvas->SetRightMargin(0.05);
1625 canvas->SetTopMargin(0.05);
1626
1627 errorResponse->SetLineColor(1);
1628 errorResponse->GetXaxis()->SetRangeUser(1, (strcmp(det, "TPC") ? 200 : 100));
1629 errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1630 errorResponse->SetStats(kFALSE);
1631 errorResponse->GetYaxis()->SetTitleOffset(1.2);
1632 errorResponse->SetTitle(";true multiplicity;#sigma(U-T)/T");
1633
1634 errorResponse->Draw();
1635
1636 errorMeasured->SetLineColor(2);
1637 errorMeasured->Draw("SAME");
1638
1639 errorBoth->SetLineColor(4);
1640 errorBoth->Draw("SAME");
1641
1642 TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
1643 TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
1644
1645 errorBoth2->SetLineColor(4);
1646 errorBoth2->SetLineStyle(2);
1647 errorBoth2->Draw("SAME");
1648
1649 TLegend* legend = new TLegend(0.2, 0.6, 0.6, 0.9);
1650 legend->SetFillColor(0);
1651 legend->AddEntry(errorResponse, "response matrix (Bayesian)");
1652 legend->AddEntry(errorMeasured, "measured (Bayesian)");
1653 legend->AddEntry(errorBoth, "both (Bayesian)");
1654 legend->AddEntry(errorBoth2, "both (#chi^{2} minimization)");
1655 legend->Draw();
1656
1657 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1658}
1659
1660void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
1661{
1662 const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
1663
1664 gSystem->Load("libPWG0base");
1665
1666 TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
1667 canvas->SetGridx();
1668 canvas->SetGridy();
1669 canvas->SetRightMargin(0.05);
1670 canvas->SetTopMargin(0.05);
1671
1672 AliMultiplicityCorrection* data[4];
1673 TH1* effArray[4];
1674
1675 Int_t markers[] = { 24, 25, 26, 5 };
1676 Int_t colors[] = { 1, 2, 3, 4 };
1677
1678 TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
1679 legend->SetFillColor(0);
1680
1681 TH1* effError = 0;
1682
1683 for (Int_t i=0; i<4; ++i)
1684 {
1685 TString name;
1686 name.Form("Multiplicity_%d", i);
1687
1688 TFile::Open(files[i]);
1689 data[i] = new AliMultiplicityCorrection(name, name);
1690
1691 if (i < 3)
1692 {
1693 data[i]->LoadHistograms("Multiplicity");
1694 }
1695 else
1696 data[i]->LoadHistograms("Multiplicity_0");
1697
1698 TH1* eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
1699 effArray[i] = eff;
1700
1701 eff->GetXaxis()->SetRangeUser(0, 15);
1702 eff->GetYaxis()->SetRangeUser(0, 1.1);
1703 eff->SetStats(kFALSE);
1704 eff->SetTitle(";true multiplicity;Efficiency");
1705 eff->SetLineColor(colors[i]);
1706 eff->SetMarkerColor(colors[i]);
1707 eff->SetMarkerStyle(markers[i]);
1708
1709 if (i == 3)
1710 {
1711 for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1712 eff->SetBinError(bin, 0);
1713
1714 // loop over cross section combinations
1715 for (Int_t j=1; j<7; ++j)
1716 {
1717 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
1718 mult->LoadHistograms(Form("Multiplicity_%d", j));
1719
1720 TH1* eff2 = mult->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType);
1721
1722 for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1723 {
1724 // TODO we could also do asymmetric errors here
1725 Float_t deviation = TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
1726
1727 eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), (Double_t) deviation));
1728 }
1729 }
1730
1731 for (Int_t bin=1; bin<=20; bin++)
1732 if (eff->GetBinContent(bin) > 0)
1733 Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1734
1735 if (uncertainty) {
1736 effError = (TH1*) eff->Clone("effError");
1737 effError->Reset();
1738
1739 for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
1740 if (eff->GetBinContent(bin) > 0)
1741 effError->SetBinContent(bin, 10.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1742
1743 effError->SetLineColor(1);
1744 effError->SetMarkerStyle(1);
1745 effError->DrawCopy("SAME HIST");
1746 }
1747 }
1748
1749 eff->SetBinContent(1, 0);
1750 eff->SetBinError(1, 0);
1751
1752 canvas->cd();
1753 if (i == 0)
1754 {
1755 eff->DrawCopy("P");
1756 }
1757 else
1758 eff->DrawCopy("SAME P");
1759
1760 legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
1761 }
1762
1763 if (uncertainty)
1764 legend->AddEntry(effError, "relative syst. uncertainty #times 10");
1765
1766 legend->Draw();
1767
1768 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1769}
1770
1771void ModelDependencyPlot()
1772{
1773 gSystem->Load("libPWG0base");
1774
1775 TFile::Open("multiplicityMC_3M.root");
1776 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1777 mult->LoadHistograms("Multiplicity");
1778
1779 TH2* proj = (TH2*) mult->GetCorrelation(3)->Project3D("zy");
1780
1781 TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 800, 400);
1782 canvas->SetGridx();
1783 canvas->SetGridy();
1784 //canvas->SetRightMargin(0.05);
1785 //canvas->SetTopMargin(0.05);
1786
1787 canvas->Divide(2, 1);
1788
1789 canvas->cd(2);
1790 gPad->SetLogy();
1791
1792 Int_t selectedMult = 30;
1793 Int_t yMax = 200000;
1794
1795 TH1* full = proj->ProjectionX("full");
1796 TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
1797
1798 full->SetStats(kFALSE);
1799 full->GetXaxis()->SetRangeUser(0, 200);
1800 full->GetYaxis()->SetRangeUser(5, yMax);
1801 full->SetTitle(";multiplicity");
1802
1803 selected->SetLineColor(0);
1804 selected->SetMarkerColor(2);
1805 selected->SetMarkerStyle(7);
1806
1807 full->Draw();
1808 selected->Draw("SAME P");
1809
1810 TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
1811 legend->SetFillColor(0);
1812 legend->AddEntry(full, "true");
1813 legend->AddEntry(selected, "measured");
1814 legend->Draw();
1815
1816 TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
1817 line->SetLineWidth(2);
1818 line->Draw();
1819
1820 canvas->cd(1);
1821 gPad->SetLogy();
1822
1823 full = proj->ProjectionY("full2");
1824 selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
1825
1826 full->SetStats(kFALSE);
1827 full->GetXaxis()->SetRangeUser(0, 200);
1828 full->GetYaxis()->SetRangeUser(5, yMax);
1829 full->SetTitle(";multiplicity");
1830
1831 full->SetLineColor(0);
1832 full->SetMarkerColor(2);
1833 full->SetMarkerStyle(7);
1834
1835 full->Draw("P");
1836 selected->Draw("SAME");
1837
1838 legend->Draw();
1839
1840 line = new TLine(selectedMult, 5, selectedMult, yMax);
1841 line->SetLineWidth(2);
1842 line->Draw();
1843
1844 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1845}
1846
1847void SystematicpTSpectrum()
1848{
1849 gSystem->Load("libPWG0base");
1850
1851 TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
1852 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1853 mult->LoadHistograms("Multiplicity");
1854
1855 TFile::Open("multiplicityMC_100k_syst.root");
1856 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1857 mult2->LoadHistograms("Multiplicity");
1858
1859 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1860 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1861 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1862
1863 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1864 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1865
1866 DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
1867}
1868
1869// to be deleted
1870/*void covMatrix(Bool_t mc = kTRUE)
1871{
1872 gSystem->Load("libPWG0base");
1873
1874 TFile::Open(correctionFile);
1875 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1876 mult->LoadHistograms("Multiplicity");
1877
1878 TFile::Open(measuredFile);
1879 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1880 mult2->LoadHistograms("Multiplicity");
1881
1882 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1883
1884 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1885
1886 mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0));
1887}*/
1888
1889Double_t FitPtFunc(Double_t *x, Double_t *par)
1890{
1891 Double_t xx = x[0];
1892
1893 Float_t val1 = par[1] + par[2] * xx + par[3] * xx * xx;
1894 Float_t val2 = TMath::Exp(par[4] + par[5] * xx);
1895
1896 const Float_t kTransitionWidth = 0;
1897
1898 // power law part
1899 if (xx < par[0] - kTransitionWidth)
1900 {
1901 return val1;
1902 }
1903 /*else if (xx < par[0] + kTransitionWidth)
1904 {
1905 // smooth transition
1906 Float_t factor = (xx - par[0] + kTransitionWidth) / kTransitionWidth / 2;
1907 return (1 - factor) * val1 + factor * val2;
1908 }*/
1909 else
1910 {
1911 return val2;
1912 }
1913}
1914
1915void FitPtNew(const char* fileName = "TruePt14TeV.root")
1916{
1917 gSystem->Load("libANALYSIS");
1918 gSystem->Load("libPWG0base");
1919
1920 TFile::Open(fileName);
1921
1922 TH1* genePt = (TH1*) gFile->Get("fHistPt");
1923 genePt->Sumw2();
1924
1925 // normalize by bin width
1926 for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1927 genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1928
1929 genePt->GetXaxis()->SetRangeUser(0.05, 2.0);
1930
1931 genePt->Scale(1.0 / genePt->Integral());
1932
1933 TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x*x)", 0.001, 100);
1934 //func->SetLineColor(2);
1935 func->SetParameters(1, -1);
1936
1937 genePt->SetMarkerStyle(25);
1938 genePt->SetTitle("");
1939 genePt->SetStats(kFALSE);
1940 genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
1941 //func->Draw("SAME");
1942
1943 genePt->Fit(func, "0", "", 0.05, 1);
1944
1945 new TCanvas;
1946 genePt->DrawCopy("P");
1947 func->SetRange(0.02, 8);
1948 func->DrawCopy("SAME");
1949 gPad->SetLogy();
1950}
1951
1952void FitPt(const char* fileName = "firstplots100k_truept.root")
1953{
1954 gSystem->Load("libPWG0base");
1955
1956 TFile::Open(fileName);
1957
1958 /*
1959 // merge corrections
1960 AliCorrection* correction[4];
1961 TList list;
1962
1963 for (Int_t i=0; i<4; ++i)
1964 {
1965 Printf("correction %d", i);
1966
1967 TString name; name.Form("correction_%d", i);
1968 correction[i] = new AliCorrection(name, name);
1969 correction[i]->LoadHistograms();
1970
1971 if (i > 0)
1972 list.Add(correction[i]);
1973 }
1974
1975 correction[0]->Merge(&list);
1976
1977 TH3* gene = correction[0]->GetTrackCorrection()->GetGeneratedHistogram();
1978
1979 // limit vtx, eta axis
1980 gene->GetXaxis()->SetRangeUser(-5.9, 5.9);
1981 gene->GetYaxis()->SetRangeUser(-1.99, 0.99);
1982
1983 TH1* genePt = gene->Project3D("z");*/
1984 TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
1985 if (!genePt)
1986 genePt = (TH1*) gFile->Get("fHistPt");
1987
1988 genePt->Sumw2();
1989
1990 //genePt->Scale(1.0 / genePt->Integral());
1991
1992 // normalize by bin width
1993 for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1994 genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1995
1996 /// genePt->GetXaxis()->GetBinCenter(x));
1997
1998 genePt->GetXaxis()->SetRangeUser(0, 7.9);
1999 //genePt->GetYaxis()->SetTitle("a.u.");
2000
2001 //TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x*x)", 0.001, 100);
2002 TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
2003 //func->SetLineColor(2);
2004 func->SetParameters(1, -1, 1, 1, 1);
2005 func->SetParLimits(3, 1, 10);
2006 func->SetParLimits(4, 0, 10);
2007
2008 //TF1* func = new TF1("func", "[1]*x**[0]", 0.001, 100);
2009
2010 //TF1* func = new TF1("func", &FitPtFunc, 0, 2, 6);
2011 //func->SetParameters(0.3, -2.34909e-01, 1.54394e+01, -3.04134e+01, 1.41912e+00, -2.79284e+00);
2012 //func->FixParameter(0, 0.314);
2013 //func->SetParLimits(0, 0.1, 0.3);
2014
2015 genePt->SetMarkerStyle(25);
2016 genePt->SetTitle("");
2017 genePt->SetStats(kFALSE);
2018 genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
2019 //func->Draw("SAME");
2020
2021 // fit only exp. part
2022 func->SetParameters(1, -1);
2023 func->FixParameter(2, 0);
2024 func->FixParameter(3, 1);
2025 func->FixParameter(4, 1);
2026 genePt->Fit(func, "0", "", 0.2, 1);
2027
2028 new TCanvas;
2029 genePt->DrawCopy("P");
2030 func->SetRange(0.02, 8);
2031 func->DrawCopy("SAME");
2032 gPad->SetLogy();
2033
2034 // now fix exp. parameters and fit second part
2035 Double_t param0 = func->GetParameter(0);
2036 Double_t param1 = func->GetParameter(1);
2037 func->SetParameters(0, -1, 1, 1, 1);
2038 func->FixParameter(0, 0);
2039 func->FixParameter(1, -1);
2040 func->ReleaseParameter(2);
2041 func->ReleaseParameter(3);
2042 func->ReleaseParameter(4);
2043 func->SetParLimits(3, 1, 10);
2044 func->SetParLimits(4, 0, 10);
2045
2046 genePt->Fit(func, "0", "", 1.5, 4);
2047
2048 new TCanvas;
2049 genePt->DrawCopy("P");
2050 func->SetRange(0.02, 8);
2051 func->DrawCopy("SAME");
2052 gPad->SetLogy();
2053
2054 // fit both
2055 func->SetParameter(0, param0);
2056 func->SetParameter(1, param1);
2057 func->ReleaseParameter(0);
2058 func->ReleaseParameter(1);
2059
2060 new TCanvas;
2061 genePt->DrawCopy("P");
2062 func->SetRange(0.02, 5);
2063 func->DrawCopy("SAME");
2064 gPad->SetLogy();
2065
2066 genePt->Fit(func, "0", "", 0.2, 4);
2067
2068 TCanvas* canvas = new TCanvas("FitPt", "FitPt", 800, 400);
2069 canvas->Divide(2, 1);
2070 canvas->cd(1);
2071
2072 gPad->SetGridx();
2073 gPad->SetGridy();
2074 gPad->SetLeftMargin(0.13);
2075 gPad->SetRightMargin(0.05);
2076 gPad->SetTopMargin(0.05);
2077
2078 genePt->GetXaxis()->SetRangeUser(0, 4.9);
2079 genePt->GetYaxis()->SetRangeUser(1e-2, 1e4);
2080 genePt->GetYaxis()->SetTitleOffset(1.4);
2081 genePt->GetXaxis()->SetTitleOffset(1.1);
2082 genePt->DrawCopy("P");
2083 func->SetRange(0.02, 5);
2084 func->DrawCopy("SAME");
2085 gPad->SetLogy();
2086
2087 canvas->cd(2);
2088
2089 TH1* genePtClone = (TH1*) genePt->Clone("genePtClone");
2090 genePtClone->Reset();
2091 genePtClone->DrawCopy("P");
2092
2093 gPad->SetGridx();
2094 gPad->SetGridy();
2095 gPad->SetLeftMargin(0.13);
2096 gPad->SetRightMargin(0.05);
2097 gPad->SetTopMargin(0.05);
2098
2099 func->DrawCopy("SAME");
2100 gPad->SetLogy();
2101
2102 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2103
2104 TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2105
2106 TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2107
2108 TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2109
2110 for (Int_t param=0; param<5; param++)
2111 {
2112 for (Int_t sign=0; sign<2; sign++)
2113 {
2114 TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign)); // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2115 func2->SetParameters(func->GetParameters());
2116 //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2117
2118 Float_t factor = ((sign == 0) ? 0.9 : 1.1);
2119 func2->SetParameter(param, func2->GetParameter(param) * factor);
2120 //func2->Print();
2121
2122 canvas->cd(2);
2123 func2->SetLineWidth(1);
2124 func2->SetLineColor(2);
2125 func2->DrawCopy("SAME");
2126
2127 canvas2->cd();
2128 TH1* second = func2->GetHistogram();
2129 second->Divide(first);
2130 second->SetLineColor(param + 1);
2131 second->GetYaxis()->SetRangeUser(0, 2);
2132 second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2133 second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2134 }
2135 }
2136
2137 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2138 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2139
2140 file->Close();
2141}
2142
2143void DrawSystematicpT()
2144{
2145 TFile* file = TFile::Open("SystematicpT.root");
2146
2147 TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2148 TH1* result2 = (TH1*) file->Get("result_unity");
2149
2150 TH1* mcHist[12];
2151 TH1* result[12];
2152
2153 Int_t nParams = 5;
2154
2155 for (Int_t id=0; id<nParams*2; ++id)
2156 {
2157 mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2158 result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2159 }
2160
2161 DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2162
2163 //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2164
2165 DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2166
2167 //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2168
2169 // does not make sense: mc is different
2170 //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2171}
2172
2173void SystematicpT(Bool_t chi2 = 1)
2174{
2175 gSystem->Load("libPWG0base");
2176
2177 TFile::Open("ptspectrum900.root");
2178 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2179 mult->LoadHistograms("Multiplicity");
2180
2181 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2182
2183 TH1* mcHist[12];
2184 TH1* result[12];
2185
2186 Int_t nParams = 5;
2187
2188 for (Int_t param=0; param<nParams; param++)
2189 {
2190 for (Int_t sign=0; sign<2; sign++)
2191 {
2192 // calculate result with systematic effect
2193 TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
2194 mult2->LoadHistograms("Multiplicity");
2195
2196 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2197
2198 if (chi2)
2199 {
2200 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2201 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2202 }
2203 else
2204 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2205
2206 Int_t id = param * 2 + sign;
2207
2208 mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2209 result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2210
2211 TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2212 DrawResultRatio(mcHist[id], result[id], tmp);
2213 }
2214 }
2215
2216 // calculate normal result
2217 TFile::Open("ptspectrum100_1.root");
2218 mult2->LoadHistograms("Multiplicity");
2219 TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
2220
2221 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2222
2223 if (chi2)
2224 {
2225 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2226 }
2227 else
2228 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2229
2230 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
2231
2232 TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2233 mcHist2->Write();
2234 result2->Write();
2235 for (Int_t id=0; id<nParams*2; ++id)
2236 {
2237 mcHist[id]->Write();
2238 result[id]->Write();
2239 }
2240 file->Close();
2241
2242 DrawSystematicpT();
2243}
2244
2245void DrawSystematicpT2()
2246{
2247 //displayRange = 200;
2248
2249 // read from file
2250 TFile* file = TFile::Open("SystematicpT2.root");
2251 TH1* mcHist = (TH1*) file->Get("mymc");
2252 TH1* result[12];
2253 result[0] = (TH1*) file->Get("result_unity");
2254 Int_t nParams = 5;
2255 for (Int_t id=0; id<nParams*2; ++id)
2256 result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
2257
2258 DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
2259 DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
2260 DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
2261}
2262
2263void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
2264{
2265 gSystem->Load("libPWG0base");
2266
2267 if (tpc)
2268 {
2269 SetTPC();
2270 TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2271 }
2272 else
2273 TFile::Open("ptspectrum100_1.root");
2274
2275 AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2276 measured->LoadHistograms("Multiplicity");
2277 TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2278
2279 TH1* result[12];
2280
2281 Int_t nParams = 5;
2282
2283 // -1 = unity change, 0...4 parameters
2284 for (Int_t id=-1; id<nParams*2; id++)
2285 {
2286 Int_t param = id / 2;
2287 Int_t sign = id % 2;
2288
2289 TString idStr;
2290 if (id == -1)
2291 {
2292 idStr = "unity";
2293 }
2294 else
2295 idStr.Form("%d_%d", param, sign);
2296
2297 // calculate result with systematic effect
2298 if (tpc)
2299 {
2300 TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
2301 }
2302 else
2303 TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
2304
2305 AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2306 response->LoadHistograms("Multiplicity");
2307
2308 response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
2309
2310 if (chi2)
2311 {
2312 response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2313 response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2314 }
2315 else
2316 response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2317
2318 result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
2319
2320 TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
2321 DrawResultRatio(mcHist, result[id+1], tmp);
2322 }
2323
2324 TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
2325 mcHist->Write();
2326 for (Int_t id=0; id<nParams*2+1; ++id)
2327 result[id]->Write();
2328 file->Close();
2329
2330 DrawSystematicpT2();
2331}
2332
2333void SystematicpTCutOff(Bool_t chi2 = kTRUE)
2334{
2335 // only needed for TPC
2336 SetTPC();
2337
2338 gSystem->Load("libPWG0base");
2339
2340 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2341 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2342 mult->LoadHistograms("Multiplicity");
2343
2344 TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2345 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2346 mult2->LoadHistograms("Multiplicity");
2347
2348 // "normal" result
2349 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2350
2351 if (chi2)
2352 {
2353 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2354 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2355 }
2356 else
2357 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2358
2359 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2360 TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
2361
2362 TH1* syst[2];
2363
2364 // change of pt spectrum (down)
2365 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2366 AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2367 mult3->LoadHistograms("Multiplicity");
2368 mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2369 if (chi2)
2370 {
2371 mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2372 mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2373 }
2374 else
2375 mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2376 syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2377
2378 // change of pt spectrum (up)
2379 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
2380 AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
2381 mult4->LoadHistograms("Multiplicity");
2382 mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2383 if (chi2)
2384 {
2385 mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2386 mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2387 }
2388 else
2389 mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2390 syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
2391
2392 DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
2393
2394 Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
2395 Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
2396}
2397
2398TH1* SystematicsSummary(Bool_t tpc = 1)
2399{
2400 Int_t nEffects = 7;
2401
2402 TH1* effects[10];
2403 const char** names = 0;
2404 Int_t colors[] = { 1, 2, 3, 4, 6, 7, 8 };
2405 Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
2406
2407 for (Int_t i=0; i<nEffects; ++i)
2408 effects[i] = new TH1F("SystematicsSummary", ";true multiplicity;Effect", 201, -0.5, 200.5);
2409
2410 if (tpc)
2411 {
2412 SetTPC();
2413
2414 const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
2415 names = namesTPC;
2416
2417 // method
2418 TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
2419 TH1* hist = (TH1*) file->Get("errorBoth");
2420
2421 // smooth a bit, but skip 0 bin
2422 effects[0]->SetBinContent(2, hist->GetBinContent(2));
2423 for (Int_t i=3; i<=200; ++i)
2424 effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2425
2426 // relative x-section
2427 effects[1]->SetBinContent(2, 0.005);
2428 effects[1]->SetBinContent(3, 0.0025);
2429 effects[1]->SetBinContent(4, 0.0025);
2430
2431 // particle composition
2432 for (Int_t i=2; i<=101; ++i)
2433 {
2434 if (i < 41)
2435 {
2436 effects[2]->SetBinContent(i, 0.01);
2437 }
2438 else if (i < 76)
2439 {
2440 effects[2]->SetBinContent(i, 0.02);
2441 }
2442 else
2443 effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
2444 }
2445
2446 // pt cut off (only tpc)
2447 for (Int_t i=2; i<=101; ++i)
2448 {
2449 if (i < 11)
2450 {
2451 effects[3]->SetBinContent(i, 0.05);
2452 }
2453 else if (i < 51)
2454 {
2455 effects[3]->SetBinContent(i, 0.01);
2456 }
2457 else
2458 effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
2459 }
2460
2461 // track selection (only tpc)
2462 for (Int_t i=2; i<=101; ++i)
2463 effects[4]->SetBinContent(i, 0.03);
2464
2465 // secondaries
2466 for (Int_t i=2; i<=101; ++i)
2467 effects[5]->SetBinContent(i, 0.01);
2468
2469 // pt spectrum
2470 for (Int_t i=2; i<=101; ++i)
2471 {
2472 if (i < 21)
2473 {
2474 effects[6]->SetBinContent(i, 0.05);
2475 }
2476 else if (i < 51)
2477 {
2478 effects[6]->SetBinContent(i, 0.02);
2479 }
2480 else
2481 effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
2482 }
2483
2484 }
2485 else
2486 {
2487 displayRange = 200;
2488 nEffects = 5;
2489
2490 const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Rel. cross-section", "Particle composition", "Secondaries", "p_{t} spectrum"};
2491 names = namesSPD;
2492
2493 // method
2494 TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
2495 TH1* hist = (TH1*) file->Get("errorBoth");
2496
2497 // smooth a bit, but skip 0 bin
2498 effects[0]->SetBinContent(2, hist->GetBinContent(2));
2499 for (Int_t i=3; i<=201; ++i)
2500 effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2501
2502 // relative x-section
2503 effects[1]->SetBinContent(2, 0.01);
2504 effects[1]->SetBinContent(3, 0.005);
2505
2506 // particle composition
2507 for (Int_t i=2; i<=201; ++i)
2508 {
2509 if (i < 6)
2510 {
2511 effects[2]->SetBinContent(i, 0.3);
2512 }
2513 else if (i < 11)
2514 {
2515 effects[2]->SetBinContent(i, 0.05);
2516 }
2517 else if (i < 121)
2518 {
2519 effects[2]->SetBinContent(i, 0.02);
2520 }
2521 else if (i < 151)
2522 {
2523 effects[2]->SetBinContent(i, 0.02 + 0.04 / 30 * (i - 121));
2524 }
2525 else
2526 effects[2]->SetBinContent(i, 0.06 + 0.1 / 30 * (i - 151));
2527 }
2528
2529 // secondaries
2530 for (Int_t i=2; i<=201; ++i)
2531 effects[3]->SetBinContent(i, 0.01);
2532
2533 // pt spectrum
2534 for (Int_t i=2; i<=201; ++i)
2535 {
2536 if (i < 6)
2537 {
2538 effects[4]->SetBinContent(i, 1);
2539 }
2540 else if (i < 121)
2541 {
2542 effects[4]->SetBinContent(i, 0.03);
2543 }
2544 else if (i < 151)
2545 {
2546 effects[4]->SetBinContent(i, 0.03 + 0.07 / 30 * (i - 121));
2547 }
2548 else
2549 effects[4]->SetBinContent(i, 0.1);
2550 }
2551 }
2552
2553 TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 400);
2554 canvas->SetRightMargin(0.25);
2555 canvas->SetTopMargin(0.05);
2556 TLegend* legend = new TLegend(0.2, 0.4, 0.5, 0.4 + 0.5 * nEffects / 7);
2557 legend->SetFillColor(0);
2558
2559 for (Int_t i=0; i<nEffects; ++i)
2560 {
2561 TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
2562 /*current->Reset();
2563 for (Int_t j=0; j<nEffects-i; ++j)
2564 current->Add(effects[j]);*/
2565
2566 current->SetLineColor(colors[i]);
2567 //current->SetFillColor(colors[i]);
2568 current->SetMarkerColor(colors[i]);
2569 //current->SetMarkerStyle(markers[i]);
2570
2571 current->SetStats(kFALSE);
2572 current->GetYaxis()->SetRangeUser(0, 0.4);
2573 current->GetXaxis()->SetRangeUser(0, displayRange);
2574 current->DrawCopy(((i == 0) ? "" : "SAME"));
2575 legend->AddEntry(current, names[i]);
2576
2577 TLatex* text = new TLatex(displayRange+5, current->GetBinContent(displayRange+1), names[i]);
2578 text->SetTextColor(colors[i]);
2579 text->Draw();
2580 }
2581
2582 // add total in square
2583 TH1* total = (TH1*) effects[0]->Clone("total");
2584 total->Reset();
2585
2586 for (Int_t i=0; i<nEffects; ++i)
2587 {
2588 //Printf("%d %f", i, effects[i]->GetBinContent(20));
2589 effects[i]->Multiply(effects[i]);
2590 total->Add(effects[i]);
2591 }
2592
2593 for (Int_t i=1; i<=total->GetNbinsX(); ++i)
2594 if (total->GetBinContent(i) > 0)
2595 total->SetBinContent(i, TMath::Min(sqrt(total->GetBinContent(i)), 1.0));
2596
2597 //Printf("%f", total->GetBinContent(20));
2598
2599 total->SetMarkerStyle(3);
2600 total->SetMarkerColor(1);
2601 legend->AddEntry(total, "total");
2602 total->DrawCopy("SAME P");
2603
2604 legend->Draw();
2605
2606 canvas->SaveAs(canvas->GetName());
2607
2608 return total;
2609}
2610
2611void finalPlot(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE, Bool_t small = kFALSE)
2612{
2613 gSystem->Load("libPWG0base");
2614
2615 if (tpc)
2616 SetTPC();
2617
2618 if (!chi2)
2619 Printf("WARNING: Bayesian set. This is only for test!");
2620
2621 // systematic error
2622 TH1* error = SystematicsSummary(tpc);
2623
2624 TFile::Open(correctionFile);
2625 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2626 mult->LoadHistograms("Multiplicity");
2627
2628 TFile::Open(measuredFile);
2629 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2630 mult2->LoadHistograms("Multiplicity");
2631
2632 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2633
2634 if (chi2)
2635 {
2636 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2637 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kINEL);
2638 }
2639 else
2640 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kINEL, 1, 100, 0, kFALSE);
2641
2642 TH1* mcHist = mult2->GetMultiplicityINEL(etaRange)->ProjectionY("mymc");
2643 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2644
2645 DrawResultRatio(mcHist, result, "finalPlotCheck.eps");
2646
2647 // normalize result
2648 result->Scale(1.0 / result->Integral(2, 200));
2649
2650 result->GetXaxis()->SetRangeUser(0, ((tpc) ? displayRange : 200));
2651 result->SetBinContent(1, 0); result->SetBinError(1, 0);
2652 result->SetTitle(";true multiplicity;Probability");
2653 result->SetLineColor(1);
2654 result->SetStats(kFALSE);
2655
2656 TH1* systError = (TH1*) result->Clone("systError");
2657 for (Int_t i=2; i<=systError->GetNbinsX(); ++i)
2658 systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
2659
2660 // change error drawing style
2661 systError->SetFillColor(15);
2662
2663 TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", (small) ? 600 : 800, 400);
2664 canvas->SetRightMargin(0.05);
2665 canvas->SetTopMargin(0.05);
2666
2667 systError->Draw("E2 ][");
2668 result->DrawCopy("SAME E ][");
2669 canvas->SetLogy();
2670
2671 //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
2672 TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
2673 text->SetFillColor(0);
2674 text->SetTextAlign(12);
2675 text->AddText("Systematic errors summed quadratically");
2676 text->AddText("0.6 million minimum bias events");
2677 text->AddText("corrected to inelastic events");
2678 text->Draw("B");
2679
2680 TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
2681 text2->SetFillColor(0);
2682 text2->SetTextAlign(12);
2683 text2->AddText("#sqrt{s} = 14 TeV");
2684 if (tpc)
2685 {
2686 text2->AddText("|#eta| < 0.9");
2687 }
2688 else
2689 text2->AddText("|#eta| < 2.0");
2690 text2->AddText("simulated data (PYTHIA)");
2691 text2->Draw("B");
2692
2693 if (tpc)
2694 {
2695 TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
2696 text3->SetNDC();
2697 text3->Draw();
2698 }
2699 else
2700 {
2701 TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
2702 text3->SetNDC();
2703 text3->Draw();
2704 }
2705
2706 // alice logo
2707 TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
2708 pad->Draw();
2709 pad->cd();
2710 TImage* img = TImage::Open("$HOME/alice.png");
2711 img->SetImageQuality(TAttImage::kImgBest);
2712 img->Draw();
2713
2714 canvas->Modified();
2715
2716/* TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
2717 text->SetTextSize(0.04);
2718 text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
2719 text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
2720 text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
2721 text->Draw();*/
2722
2723
2724 canvas->SaveAs(canvas->GetName());
2725}
2726
2727void BlobelUnfoldingExample()
2728{
2729 const Int_t kSize = 20;
2730
2731 TMatrixD matrix(kSize, kSize);
2732 for (Int_t x=0; x<kSize; x++)
2733 {
2734 for (Int_t y=0; y<kSize; y++)
2735 {
2736 if (x == y)
2737 {
2738 if (x == 0 || x == kSize -1)
2739 {
2740 matrix(x, y) = 0.75;
2741 }
2742 else
2743 matrix(x, y) = 0.5;
2744 }
2745 else if (TMath::Abs(x - y) == 1)
2746 {
2747 matrix(x, y) = 0.25;
2748 }
2749 }
2750 }
2751
2752 //matrix.Print();
2753
2754 TMatrixD inverted(matrix);
2755 inverted.Invert();
2756
2757 //inverted.Print();
2758
2759 TH1F* inputDist = new TH1F("inputDist", ";t;#tilde{T}(t)", kSize, -0.5, (Float_t) kSize - 0.5);
2760 TVectorD inputDistVector(kSize);
2761 TH1F* unfolded = inputDist->Clone("unfolded");
2762 TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
2763 measuredIdealDist->SetTitle(";m;#tilde{M}(m)");
2764 TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
2765
2766 TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
2767 // norm: 1/(sqrt(2pi)sigma)
2768 gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
2769 //gaus->Print();
2770
2771 for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
2772 {
2773 Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
2774 inputDist->SetBinContent(x, value);
2775 inputDistVector(x-1) = value;
2776 }
2777
2778 TVectorD measuredDistIdealVector = matrix * inputDistVector;
2779
2780 for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
2781 measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
2782
2783 measuredDist->FillRandom(measuredIdealDist, 10000);
2784
2785 // fill error matrix before scaling
2786 TMatrixD covarianceMatrixMeasured(kSize, kSize);
2787 for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
2788 covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
2789
2790 TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
2791 //covarianceMatrix.Print();
2792
2793 TVectorD measuredDistVector(kSize);
2794 for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
2795 measuredDistVector(x-1) = measuredDist->GetBinContent(x);
2796
2797 TVectorD unfoldedVector = inverted * measuredDistVector;
2798 for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
2799 unfolded->SetBinContent(x, unfoldedVector(x-1));
2800
2801 TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1000, 500);
2802 canvas->SetTopMargin(0.05);
2803 canvas->Divide(2, 1);
2804
2805 canvas->cd(1);
2806 canvas->cd(1)->SetLeftMargin(0.15);
2807 canvas->cd(1)->SetRightMargin(0.05);
2808 measuredDist->GetYaxis()->SetTitleOffset(1.7);
2809 measuredDist->SetStats(0);
2810 measuredDist->DrawCopy();
2811 gaus->Draw("SAME");
2812
2813 canvas->cd(2);
2814 canvas->cd(2)->SetLeftMargin(0.15);
2815 canvas->cd(2)->SetRightMargin(0.05);
2816 unfolded->GetYaxis()->SetTitleOffset(1.7);
2817 unfolded->SetStats(0);
2818 unfolded->DrawCopy();
2819 gaus->Draw("SAME");
2820
2821 canvas->SaveAs("BlobelUnfoldingExample.eps");
2822}
2823
2824void E735Fit()
2825{
2826 TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
2827 fCurrentESD->Sumw2();
2828
2829 // Open the input stream
2830 ifstream in;
2831 in.open("e735data.txt");
2832
2833 while(in.good())
2834 {
2835 Float_t x, y, ye;
2836 in >> x >> y >> ye;
2837
2838 //Printf("%f %f %f", x, y, ye);
2839 fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
2840 fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
2841 }
2842
2843 in.close();
2844
2845 //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
2846
2847 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2848
2849 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])");
2850 func->SetParNames("scaling", "averagen", "k");
2851 func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
2852 func->SetParLimits(1, 0.001, 1000);
2853 func->SetParLimits(2, 0.001, 1000);
2854 func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
2855
2856 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
2857 lognormal->SetParNames("scaling", "mean", "sigma");
2858 lognormal->SetParameters(1, 1, 1);
2859 lognormal->SetParLimits(0, 0, 10);
2860 lognormal->SetParLimits(1, 0, 100);
2861 lognormal->SetParLimits(2, 1e-3, 10);
2862
2863 TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
2864 fCurrentESD->SetStats(kFALSE);
2865 fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
2866 fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
2867 fCurrentESD->Draw("");
2868 fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
2869 fCurrentESD->Fit(func, "0", "", 0, 150);
2870 func->SetRange(0, 250);
2871 func->Draw("SAME");
2872 printf("chi2 = %f\n", func->GetChisquare());
2873
2874 fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
2875 lognormal->SetLineColor(2);
2876 lognormal->SetLineStyle(2);
2877 lognormal->SetRange(0, 250);
2878 lognormal->Draw("SAME");
2879
2880 gPad->SetLogy();
2881
2882 canvas->SaveAs("E735Fit.eps");
2883}