]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG0/dNdEta/plotsMultiplicity.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bugfix: duplicated line
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / plotsMultiplicity.C
CommitLineData
0b4bfd98 1//
2// plots for the note about multplicity measurements
3//
4
6d81c2de 5#if !defined(__CINT__) || defined(__MAKECINT__)
6
7#include <TCanvas.h>
8#include <TPad.h>
9#include <TH1F.h>
10#include <TH2F.h>
11#include <TH3F.h>
12#include <TLine.h>
13#include <TF1.h>
14#include <TSystem.h>
15#include <TFile.h>
16#include <TLegend.h>
17#include <TStopwatch.h>
18#include <TROOT.h>
19#include <TGraph.h>
20#include <TMath.h>
21
22#include "AliMultiplicityCorrection.h"
23#include "AliCorrection.h"
24#include "AliCorrectionMatrix3D.h"
25
26#endif
27
0b4bfd98 28const char* correctionFile = "multiplicityMC_2M.root";
29const char* measuredFile = "multiplicityMC_1M_3.root";
30Int_t etaRange = 3;
6d81c2de 31Int_t drawRatioRange = 150;
0b4bfd98 32
33const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
34const char* measuredFileTPC = "multiplicityMC_TPC_0.6M.root";
35Int_t etaRangeTPC = 1;
36
37void SetTPC()
38{
39 correctionFile = correctionFileTPC;
40 measuredFile = measuredFileTPC;
41 etaRange = etaRangeTPC;
6d81c2de 42 drawRatioRange = 100;
43}
44
45void Smooth(TH1* hist, Int_t windowWidth = 20)
46{
47 TH1* clone = (TH1*) hist->Clone("clone");
48 for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin)
49 {
50 Int_t min = TMath::Max(2, bin-windowWidth);
51 Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
52 Float_t average = clone->Integral(min, max) / (max - min + 1);
53
54 hist->SetBinContent(bin, average);
55 hist->SetBinError(bin, 0);
56 }
57
58 delete clone;
0b4bfd98 59}
60
61void responseMatrixPlot()
62{
63 gSystem->Load("libPWG0base");
64
65 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
66
67 TFile::Open(correctionFile);
68 mult->LoadHistograms("Multiplicity");
69
70 TH1* hist = mult->GetCorrelation(etaRange)->Project3D("zy");
71 hist->SetStats(kFALSE);
72
73 hist->SetTitle(";true multiplicity;measured multiplicity;Entries");
74 hist->GetXaxis()->SetRangeUser(0, 200);
75 hist->GetYaxis()->SetRangeUser(0, 200);
76
77 TCanvas* canvas = new TCanvas("c1", "c1", 800, 600);
78 canvas->SetRightMargin(0.15);
79 canvas->SetTopMargin(0.05);
80
81 gPad->SetLogz();
82 hist->Draw("COLZ");
83
84 canvas->SaveAs("responsematrix.eps");
85}
86
87TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
88{
89 // normalize unfolded result to mc hist
90 result->Scale(1.0 / result->Integral(2, 200));
91 result->Scale(mcHist->Integral(2, 200));
92
93 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
94 canvas->Range(0, 0, 1, 1);
95
6d81c2de 96 TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
0b4bfd98 97 pad1->Draw();
98
6d81c2de 99 TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
0b4bfd98 100 pad2->Draw();
101
102 pad1->SetRightMargin(0.05);
103 pad2->SetRightMargin(0.05);
104
105 // no border between them
106 pad1->SetBottomMargin(0);
107 pad2->SetTopMargin(0);
108
109 pad1->cd();
110
111 mcHist->GetXaxis()->SetLabelSize(0.06);
112 mcHist->GetYaxis()->SetLabelSize(0.06);
113 mcHist->GetXaxis()->SetTitleSize(0.06);
114 mcHist->GetYaxis()->SetTitleSize(0.06);
115 mcHist->GetYaxis()->SetTitleOffset(0.6);
116
117 mcHist->GetXaxis()->SetRangeUser(0, 200);
118
119 mcHist->SetTitle(";true multiplicity;Entries");
120 mcHist->SetStats(kFALSE);
121
122 mcHist->DrawCopy("HIST E");
123 gPad->SetLogy();
124
125 result->SetLineColor(2);
126 result->DrawCopy("SAME HISTE");
127
128 TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
129 legend->AddEntry(mcHist, "true distribution");
130 legend->AddEntry(result, "unfolded distribution");
131 legend->SetFillColor(0);
132 legend->Draw();
133
134 pad2->cd();
135 pad2->SetBottomMargin(0.15);
136
137 // calculate ratio
138 mcHist->Sumw2();
139 TH1* ratio = (TH1*) mcHist->Clone("ratio");
140 result->Sumw2();
141 ratio->Divide(ratio, result, 1, 1, "");
142 ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
143 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
144
145 ratio->DrawCopy();
146
147 // get average of ratio
148 Float_t sum = 0;
149 for (Int_t i=2; i<=150; ++i)
150 {
151 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
152 }
153 sum /= 149;
154
155 printf("Average (2..150) of |ratio - 1| is %f\n", sum);
156
157 TLine* line = new TLine(0, 1, 200, 1);
158 line->SetLineWidth(2);
159 line->Draw();
160
161 line = new TLine(0, 1.1, 200, 1.1);
162 line->SetLineWidth(2);
163 line->SetLineStyle(2);
164 line->Draw();
165 line = new TLine(0, 0.9, 200, 0.9);
166 line->SetLineWidth(2);
167 line->SetLineStyle(2);
168 line->Draw();
169
170 canvas->Modified();
171
172 canvas->SaveAs(epsName);
173
174 return canvas;
175}
176
177TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
178{
179 // draws the 3 plots in the upper plot
180 // draws the ratio between result1 and result2 in the lower plot
181
182 // normalize unfolded result to mc hist
183 result1->Scale(1.0 / result1->Integral(2, 200));
184 result1->Scale(mcHist->Integral(2, 200));
185 result2->Scale(1.0 / result2->Integral(2, 200));
186 result2->Scale(mcHist->Integral(2, 200));
187
188 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
189 canvas->Range(0, 0, 1, 1);
190
6d81c2de 191 TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
0b4bfd98 192 pad1->Draw();
193
6d81c2de 194 TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
0b4bfd98 195 pad2->Draw();
196
197 pad1->SetRightMargin(0.05);
198 pad2->SetRightMargin(0.05);
199
200 // no border between them
201 pad1->SetBottomMargin(0);
202 pad2->SetTopMargin(0);
203
204 pad1->cd();
205
206 mcHist->GetXaxis()->SetLabelSize(0.06);
207 mcHist->GetYaxis()->SetLabelSize(0.06);
208 mcHist->GetXaxis()->SetTitleSize(0.06);
209 mcHist->GetYaxis()->SetTitleSize(0.06);
210 mcHist->GetYaxis()->SetTitleOffset(0.6);
211
6d81c2de 212 mcHist->GetXaxis()->SetRangeUser(0, drawRatioRange);
0b4bfd98 213
214 mcHist->SetTitle(";true multiplicity;Entries");
215 mcHist->SetStats(kFALSE);
216
217 mcHist->DrawCopy("HIST E");
218 gPad->SetLogy();
219
220 result1->SetLineColor(2);
221 result1->DrawCopy("SAME HISTE");
222
223 result2->SetLineColor(4);
224 result2->DrawCopy("SAME HISTE");
225
226 TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
227 legend->AddEntry(mcHist, "true distribution");
228 legend->AddEntry(result1, "unfolded distribution (syst)");
229 legend->AddEntry(result2, "unfolded distribution (normal)");
230 legend->SetFillColor(0);
231 legend->Draw();
232
233 pad2->cd();
234 pad2->SetBottomMargin(0.15);
235
236 result1->GetXaxis()->SetLabelSize(0.06);
237 result1->GetYaxis()->SetLabelSize(0.06);
238 result1->GetXaxis()->SetTitleSize(0.06);
239 result1->GetYaxis()->SetTitleSize(0.06);
240 result1->GetYaxis()->SetTitleOffset(0.6);
241
6d81c2de 242 result1->GetXaxis()->SetRangeUser(0, drawRatioRange);
0b4bfd98 243
244 result1->SetTitle(";true multiplicity;Entries");
245 result1->SetStats(kFALSE);
246
247 // calculate ratio
248 result1->Sumw2();
249 TH1* ratio = (TH1*) result1->Clone("ratio");
250 result2->Sumw2();
251 ratio->Divide(ratio, result2, 1, 1, "");
252 ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
253 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
254
255 ratio->DrawCopy();
256
257 // get average of ratio
258 Float_t sum = 0;
6d81c2de 259 for (Int_t i=2; i<=drawRatioRange; ++i)
0b4bfd98 260 {
261 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
262 }
6d81c2de 263 sum /= drawRatioRange-1;
0b4bfd98 264
6d81c2de 265 printf("Average (2..%d) of |ratio - 1| is %f\n", drawRatioRange, sum);
0b4bfd98 266
6d81c2de 267 TLine* line = new TLine(0, 1, drawRatioRange, 1);
0b4bfd98 268 line->SetLineWidth(2);
269 line->Draw();
270
6d81c2de 271 line = new TLine(0, 1.1, drawRatioRange, 1.1);
0b4bfd98 272 line->SetLineWidth(2);
273 line->SetLineStyle(2);
274 line->Draw();
6d81c2de 275 line = new TLine(0, 0.9, drawRatioRange, 0.9);
0b4bfd98 276 line->SetLineWidth(2);
277 line->SetLineStyle(2);
278 line->Draw();
279
280 canvas->Modified();
281
282 canvas->SaveAs(epsName);
283
284 return canvas;
285}
286
6d81c2de 287TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0)
0b4bfd98 288{
289 // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
290 // a systematic effect
291
292 // normalize results
293 result->Scale(1.0 / result->Integral(2, 200));
294
295 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
296
6d81c2de 297 result->GetXaxis()->SetRangeUser(1, drawRatioRange);
0b4bfd98 298 result->SetStats(kFALSE);
299
6d81c2de 300 Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
301
302 TLegend* legend = new TLegend(0.2, 0.75, 0.35, 0.95);
303 legend->SetFillColor(0);
304
0b4bfd98 305 for (Int_t n=0; n<nResultSyst; ++n)
306 {
307 resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
308
309 // calculate ratio
310 TH1* ratio = (TH1*) result->Clone("ratio");
311 ratio->Divide(ratio, resultSyst[n], 1, 1, "B");
312 ratio->SetTitle(";true multiplicity;Ratio");
313 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
314
315 if (firstMarker)
316 ratio->SetMarkerStyle(5);
317
6d81c2de 318 ratio->SetLineColor(colors[n / 2]);
319 if ((n % 2))
320 ratio->SetLineStyle(2);
0b4bfd98 321
322 ratio->DrawCopy((n == 0) ? ((firstMarker) ? "P" : "HIST") : "SAME HIST");
323
6d81c2de 324 if (legendStrings && legendStrings[n])
325 legend->AddEntry(ratio, legendStrings[n]);
326
0b4bfd98 327 // get average of ratio
328 Float_t sum = 0;
6d81c2de 329 for (Int_t i=2; i<=drawRatioRange; ++i)
0b4bfd98 330 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
6d81c2de 331 sum /= drawRatioRange-1;
0b4bfd98 332
6d81c2de 333 printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, drawRatioRange, sum);
0b4bfd98 334 }
335
6d81c2de 336 if (legendStrings)
337 legend->Draw();
338
339 TLine* line = new TLine(1, 1, drawRatioRange, 1);
0b4bfd98 340 line->SetLineWidth(2);
341 line->Draw();
342
6d81c2de 343 line = new TLine(1, 1.1, drawRatioRange, 1.1);
0b4bfd98 344 line->SetLineWidth(2);
345 line->SetLineStyle(2);
346 line->Draw();
6d81c2de 347 line = new TLine(1, 0.9, drawRatioRange, 0.9);
0b4bfd98 348 line->SetLineWidth(2);
349 line->SetLineStyle(2);
350 line->Draw();
351
352 canvas->Modified();
353
354 canvas->SaveAs(epsName);
355 canvas->SaveAs(Form("%s.gif", epsName.Data()));
356
357 return canvas;
358}
359
6d81c2de 360TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE)
0b4bfd98 361{
362 // draws the ratios of each mc to the corresponding result
363
364 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
6d81c2de 365 canvas->SetRightMargin(0.05);
366 canvas->SetTopMargin(0.05);
0b4bfd98 367
368 for (Int_t n=0; n<nResultSyst; ++n)
369 {
370 // normalize
371 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
372 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
373
6d81c2de 374 result[n]->GetXaxis()->SetRangeUser(0, drawRatioRange);
0b4bfd98 375 result[n]->SetStats(kFALSE);
376
377 // calculate ratio
378 TH1* ratio = (TH1*) result[n]->Clone("ratio");
379 ratio->Divide(mc[n], ratio, 1, 1, "B");
6d81c2de 380
381 // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis...
382 ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0);
383
384 if (smooth)
385 Smooth(ratio);
386
387 ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : "")));
0b4bfd98 388 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
389
6d81c2de 390 if (dashed)
391 {
392 ratio->SetLineColor((n/2)+1);
393 ratio->SetLineStyle((n%2)+1);
394 }
395 else
396 ratio->SetLineColor(n+1);
0b4bfd98 397
398 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
399
400 // get average of ratio
401 Float_t sum = 0;
6d81c2de 402 for (Int_t i=2; i<=drawRatioRange; ++i)
0b4bfd98 403 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
6d81c2de 404 sum /= drawRatioRange-1;
0b4bfd98 405
6d81c2de 406 printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, drawRatioRange, sum);
0b4bfd98 407 }
408
6d81c2de 409 TLine* line = new TLine(0, 1, drawRatioRange, 1);
0b4bfd98 410 line->SetLineWidth(2);
411 line->Draw();
412
6d81c2de 413 line = new TLine(0, 1.1, drawRatioRange, 1.1);
0b4bfd98 414 line->SetLineWidth(2);
415 line->SetLineStyle(2);
416 line->Draw();
6d81c2de 417 line = new TLine(0, 0.9, drawRatioRange, 0.9);
0b4bfd98 418 line->SetLineWidth(2);
419 line->SetLineStyle(2);
420 line->Draw();
421
422 canvas->Modified();
423
424 canvas->SaveAs(epsName);
425 canvas->SaveAs(Form("%s.gif", epsName.Data()));
426
427 return canvas;
428}
429
430TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
431{
432 // draws the ratios of each mc to the corresponding result
433 // deducts from each ratio the ratio of mcBase / resultBase
434
435 // normalize
436 resultBase->Scale(1.0 / resultBase->Integral(2, 200));
437 mcBase->Scale(1.0 / mcBase->Integral(2, 200));
438
439 // calculate ratio
440 TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
441 ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
442
443 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
444 canvas->SetRightMargin(0.05);
445 canvas->SetTopMargin(0.05);
446
447 for (Int_t n=0; n<nResultSyst; ++n)
448 {
449 // normalize
450 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
451 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
452
6d81c2de 453 result[n]->GetXaxis()->SetRangeUser(0, drawRatioRange);
0b4bfd98 454 result[n]->SetStats(kFALSE);
455
456 // calculate ratio
457 TH1* ratio = (TH1*) result[n]->Clone("ratio");
458 ratio->Divide(mc[n], ratio, 1, 1, "B");
459 ratio->Add(ratioBase, -1);
460
461 ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
462 ratio->GetYaxis()->SetRangeUser(-1, 1);
463 ratio->SetLineColor(n+1);
464 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
465
466 // get average of ratio
467 Float_t sum = 0;
6d81c2de 468 for (Int_t i=2; i<=drawRatioRange; ++i)
0b4bfd98 469 sum += TMath::Abs(ratio->GetBinContent(i));
6d81c2de 470 sum /= drawRatioRange-1;
0b4bfd98 471
6d81c2de 472 printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, drawRatioRange, sum);
0b4bfd98 473 }
474
6d81c2de 475 TLine* line = new TLine(0, 0, drawRatioRange, 0);
0b4bfd98 476 line->SetLineWidth(2);
477 line->Draw();
478
6d81c2de 479 line = new TLine(0, 0.1, drawRatioRange, 0.1);
0b4bfd98 480 line->SetLineWidth(2);
481 line->SetLineStyle(2);
482 line->Draw();
6d81c2de 483 line = new TLine(0, -0.1, drawRatioRange, -0.1);
0b4bfd98 484 line->SetLineWidth(2);
485 line->SetLineStyle(2);
486 line->Draw();
487
488 canvas->Modified();
489
490 canvas->SaveAs(epsName);
491 canvas->SaveAs(Form("%s.gif", epsName.Data()));
492
493 return canvas;
494}
495
0b4bfd98 496TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
497{
498 // draws the ratios of each mc to the corresponding result
499 // deducts from each ratio the ratio of mcBase / resultBase
500 // smoothens the ratios by a sliding window
501
502 // normalize
503 resultBase->Scale(1.0 / resultBase->Integral(2, 200));
504 mcBase->Scale(1.0 / mcBase->Integral(2, 200));
505
506 // calculate ratio
507 TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
508 ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
509
510 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
511 canvas->SetRightMargin(0.05);
512 canvas->SetTopMargin(0.05);
513
514 for (Int_t n=0; n<nResultSyst; ++n)
515 {
516 // normalize
517 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
518 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
519
6d81c2de 520 result[n]->GetXaxis()->SetRangeUser(0, drawRatioRange);
0b4bfd98 521 result[n]->SetStats(kFALSE);
522
523 // calculate ratio
524 TH1* ratio = (TH1*) result[n]->Clone("ratio");
525 ratio->Divide(mc[n], ratio, 1, 1, "B");
526 ratio->Add(ratioBase, -1);
527
6d81c2de 528 //new TCanvas; ratio->DrawCopy();
529 // clear 0 bin
530 ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0);
531
0b4bfd98 532 Smooth(ratio);
6d81c2de 533
534 //ratio->SetLineColor(1); ratio->DrawCopy("SAME");
0b4bfd98 535
536 canvas->cd();
6d81c2de 537 ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)");
538 ratio->GetYaxis()->SetRangeUser(-0.3, 0.3);
539 ratio->SetLineColor((n / 2)+1);
540 ratio->SetLineStyle((n % 2)+1);
0b4bfd98 541 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
542
543 // get average of ratio
544 Float_t sum = 0;
545 for (Int_t i=2; i<=150; ++i)
546 sum += TMath::Abs(ratio->GetBinContent(i));
547 sum /= 149;
548
549 printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
550 }
551
6d81c2de 552 TLine* line = new TLine(0, 0, drawRatioRange, 0);
0b4bfd98 553 line->SetLineWidth(2);
554 line->Draw();
555
6d81c2de 556 line = new TLine(0, 0.1, drawRatioRange, 0.1);
0b4bfd98 557 line->SetLineWidth(2);
558 line->SetLineStyle(2);
559 line->Draw();
6d81c2de 560 line = new TLine(0, -0.1, drawRatioRange, -0.1);
0b4bfd98 561 line->SetLineWidth(2);
562 line->SetLineStyle(2);
563 line->Draw();
564
565 canvas->Modified();
566
567 canvas->SaveAs(epsName);
568 canvas->SaveAs(Form("%s.gif", epsName.Data()));
569
570 return canvas;
571}
572
573void DrawResiduals(TH1* measured, TH1* unfoldedFolded, const char* epsName)
574{
575 // normalize
576 unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, 200));
577 unfoldedFolded->Scale(measured->Integral(2, 200));
578
579 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
580 canvas->Range(0, 0, 1, 1);
581
582 TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
583 pad1->Draw();
584 pad1->SetGridx();
585 pad1->SetGridy();
586
587 TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
588 pad2->Draw();
589 pad2->SetGridx();
590 pad2->SetGridy();
591
592 TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
593 pad3->SetGridx();
594 pad3->SetGridy();
595 pad3->SetRightMargin(0.05);
596 pad3->SetTopMargin(0.05);
597 pad3->Draw();
598
599 pad1->SetRightMargin(0.05);
600 pad2->SetRightMargin(0.05);
601
602 // no border between them
603 pad1->SetBottomMargin(0);
604 pad2->SetTopMargin(0);
605
606 pad1->cd();
607
608 measured->GetXaxis()->SetLabelSize(0.06);
609 measured->GetYaxis()->SetLabelSize(0.06);
610 measured->GetXaxis()->SetTitleSize(0.06);
611 measured->GetYaxis()->SetTitleSize(0.06);
612 measured->GetYaxis()->SetTitleOffset(0.6);
613
614 measured->GetXaxis()->SetRangeUser(0, 150);
615
616 measured->SetTitle(";measured multiplicity;Entries");
617 measured->SetStats(kFALSE);
618
619 measured->DrawCopy("HIST");
620 gPad->SetLogy();
621
622 unfoldedFolded->SetMarkerStyle(5);
623 unfoldedFolded->SetMarkerColor(2);
624 unfoldedFolded->SetLineColor(0);
625 unfoldedFolded->DrawCopy("SAME P");
626
627 TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
628 legend->AddEntry(measured, "measured distribution");
629 legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution");
630 legend->SetFillColor(0);
631 legend->Draw();
632
633 pad2->cd();
634 pad2->SetBottomMargin(0.15);
635
636 // calculate ratio
637 measured->Sumw2();
638 TH1* residual = (TH1*) measured->Clone("residual");
639 unfoldedFolded->Sumw2();
640
641 residual->Add(unfoldedFolded, -1);
642
643 // projection
644 TH1* residualHist = new TH1F("residualHist", ";", 15, -3, 3);
645
646 for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
647 {
648 if (measured->GetBinError(i) > 0)
649 {
650 residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
651 residual->SetBinError(i, 1);
652
653 residualHist->Fill(residual->GetBinContent(i));
654 }
655 else
656 {
657 residual->SetBinContent(i, 0);
658 residual->SetBinError(i, 0);
659 }
660 }
661
662 residual->GetYaxis()->SetTitle("Residuals 1/e (M - R #otimes U)");
663 residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
664 residual->DrawCopy();
665
666 TLine* line = new TLine(-0.5, 0, 150.5, 0);
667 line->SetLineWidth(2);
668 line->Draw();
669
670 pad3->cd();
671 residualHist->SetStats(kFALSE);
672 residualHist->GetXaxis()->SetLabelSize(0.08);
673 residualHist->Fit("gaus");
674 residualHist->Draw();
675
676 canvas->Modified();
677 canvas->SaveAs(canvas->GetName());
678
679 //const char* epsName2 = "proj.eps";
680 //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
681 //canvas->SetGridx();
682 //canvas->SetGridy();
683
684 //canvas->SaveAs(canvas->GetName());
685}
686
687void bayesianExample()
688{
689 TStopwatch watch;
690 watch.Start();
691
692 gSystem->Load("libPWG0base");
693
694 TFile::Open(correctionFile);
695 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
696 mult->LoadHistograms("Multiplicity");
697
698 TFile::Open(measuredFile);
699 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
700 mult2->LoadHistograms("Multiplicity");
701
702 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
703
704 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
705
706 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
707 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
708
709 //mult->DrawComparison("bayesianExample", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
710 DrawResultRatio(mcHist, result, "bayesianExample.eps");
711
712 // draw residual plot
713
714 // TODO take out efficiency correction if other than AliMultiplicityCorrection::kTrVtx
715 TH2* convoluted = mult->CalculateMultiplicityESD(result, etaRange);
716 TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
717
718 TH1* measured = mult2->GetMultiplicityESD(etaRange)->ProjectionY("measured");
719
720 DrawResiduals(measured, convolutedProj, "bayesianResiduals.eps");
721
722 watch.Stop();
723 watch.Print();
724}
725
726void chi2FluctuationResult()
727{
728 gSystem->Load("libPWG0base");
729
730 TFile::Open(correctionFile);
731 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
732 mult->LoadHistograms("Multiplicity");
733
734 TFile::Open(measuredFile);
735 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
736 mult2->LoadHistograms("Multiplicity");
737
738 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
739 mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
740 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
741
742 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
6d81c2de 743 //TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
0b4bfd98 744
745 mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
746
747 TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_3");
748 canvas->SaveAs("chi2FluctuationResult.eps");
749}
750
751void chi2Example()
752{
753 TStopwatch watch;
754 watch.Start();
755
756 gSystem->Load("libPWG0base");
757
758 TFile::Open(correctionFile);
759 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
760 mult->LoadHistograms("Multiplicity");
761
762 TFile::Open(measuredFile);
763 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
764 mult2->LoadHistograms("Multiplicity");
765
766 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
767 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
768 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
769
770 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
771 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
772
773 DrawResultRatio(mcHist, result, "chi2Example.eps");
774
775 watch.Stop();
776 watch.Print();
777}
778
779void chi2ExampleTPC()
780{
781 TStopwatch watch;
782 watch.Start();
783
784 gSystem->Load("libPWG0base");
785
786 TFile::Open(correctionFileTPC);
787 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
788 mult->LoadHistograms("Multiplicity");
789
790 TFile::Open(measuredFileTPC);
791 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
792 mult2->LoadHistograms("Multiplicity");
793
794 mult->SetMultiplicityESD(etaRangeTPC, mult2->GetMultiplicityESD(etaRangeTPC));
795 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
796 mult->ApplyMinuitFit(etaRangeTPC, kFALSE, AliMultiplicityCorrection::kTrVtx);
797
798 TH1* mcHist = mult2->GetMultiplicityVtx(etaRangeTPC)->ProjectionY("mymc");
799 TH1* result = mult->GetMultiplicityESDCorrected(etaRangeTPC);
800
801 DrawResultRatio(mcHist, result, "chi2ExampleTPC.eps");
802
803 watch.Stop();
804 watch.Print();
805}
806
807void bayesianNBD()
808{
809 gSystem->Load("libPWG0base");
810 TFile::Open("multiplicityMC_3M.root");
811
812 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
813 mult->LoadHistograms("Multiplicity");
814
815 TFile::Open("multiplicityMC_3M_NBD.root");
816 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
817 mult2->LoadHistograms("Multiplicity");
818
819 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
820 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1.0, 100);
821
822 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
823
824 mcHist->Sumw2();
825 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
826
827 //mult->DrawComparison("bayesianNBD", etaRange, kFALSE, kTRUE, mcHist);
828 DrawResultRatio(mcHist, result, "bayesianNBD.eps");
829}
830
831void minimizationNBD()
832{
833 gSystem->Load("libPWG0base");
834 TFile::Open("multiplicityMC_3M.root");
835
836 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
837 mult->LoadHistograms("Multiplicity");
838
839 TFile::Open("multiplicityMC_3M_NBD.root");
840 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
841 mult2->LoadHistograms("Multiplicity");
842
843 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
844 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
845 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
846
847 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
848
849 mcHist->Sumw2();
850 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
851
852 //mult->DrawComparison("minimizationNBD", etaRange, kFALSE, kTRUE, mcHist);
853 DrawResultRatio(mcHist, result, "minimizationNBD.eps");
854}
855
856void minimizationInfluenceAlpha()
857{
858 gSystem->Load("libPWG0base");
859
860 TFile::Open(measuredFile);
861 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
862 mult2->LoadHistograms("Multiplicity");
863
864 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
865 mcHist->Scale(1.0 / mcHist->Integral());
866 mcHist->GetXaxis()->SetRangeUser(0, 200);
867 mcHist->SetStats(kFALSE);
868 mcHist->SetTitle(";true multiplicity;P_{N}");
869
870 TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 300);
871 canvas->Divide(3, 1);
872
873 TFile::Open("eval-2M-1M/EvaluateChi2MethodDetail.root");
874
6d81c2de 875 TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_100.000000");
876 TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_03_2_100000.000000");
877 TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_06_2_100000000.000000");
0b4bfd98 878
879 mcHist->Rebin(2); mcHist->Scale(0.5);
880 hist1->Rebin(2); hist1->Scale(0.5);
881 hist2->Rebin(2); hist2->Scale(0.5);
882 hist3->Rebin(2); hist3->Scale(0.5);
883
884 mcHist->GetXaxis()->SetRangeUser(0, 200);
885
886 canvas->cd(1);
887 gPad->SetLogy();
888 mcHist->Draw();
889 hist1->SetMarkerStyle(5);
890 hist1->SetMarkerColor(2);
891 hist1->Draw("SAME PE");
892
893 canvas->cd(2);
894 gPad->SetLogy();
895 mcHist->Draw();
896 hist2->SetMarkerStyle(5);
897 hist2->SetMarkerColor(2);
898 hist2->Draw("SAME PE");
899
900 canvas->cd(3);
901 gPad->SetLogy();
902 mcHist->Draw();
903 hist3->SetMarkerStyle(5);
904 hist3->SetMarkerColor(2);
905 hist3->Draw("SAME PE");
906
907 canvas->SaveAs("minimizationInfluenceAlpha.eps");
908}
909
910void NBDFit()
911{
912 gSystem->Load("libPWG0base");
913
914 TFile::Open(correctionFile);
915 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
916 mult->LoadHistograms("Multiplicity");
917
918 TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
919 fCurrentESD->Sumw2();
920 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
921
922 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])");
923 func->SetParNames("scaling", "averagen", "k");
924 func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
925 func->SetParLimits(1, 0.001, 1000);
926 func->SetParLimits(2, 0.001, 1000);
927 func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
928
929 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
930 lognormal->SetParNames("scaling", "mean", "sigma");
931 lognormal->SetParameters(1, 1, 1);
932 lognormal->SetParLimits(0, 0, 10);
933 lognormal->SetParLimits(1, 0, 100);
934 lognormal->SetParLimits(2, 1e-3, 10);
935
936 TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
937 fCurrentESD->SetStats(kFALSE);
938 fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
939 fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
940 fCurrentESD->Draw("HIST");
941 fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
942 fCurrentESD->Fit(func, "W0", "", 0, 50);
943 func->SetRange(0, 100);
944 func->Draw("SAME");
945 printf("chi2 = %f\n", func->GetChisquare());
946
947 fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
948 lognormal->SetLineColor(2);
949 lognormal->SetLineStyle(2);
950 lognormal->SetRange(0, 100);
951 lognormal->Draw("SAME");
952
953 canvas->SaveAs("NBDFit.eps");
954}
955
956void DifferentSamples()
957{
958 // data generated by runMultiplicitySelector.C DifferentSamples
959
960 const char* name = "DifferentSamples";
961
962 TFile* file = TFile::Open(Form("%s.root", name));
963
964 TCanvas* canvas = new TCanvas(name, name, 800, 600);
965 canvas->Divide(2, 2);
966
967 for (Int_t i=0; i<4; ++i)
968 {
969 canvas->cd(i+1);
970 gPad->SetTopMargin(0.05);
971 gPad->SetRightMargin(0.05);
972 TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", i));
973 TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", i));
974 TH1* mc = (TH1*) file->Get(Form("mc_%d", i));
975 mc->Sumw2();
976
977 chi2Result->Divide(chi2Result, mc, 1, 1, "");
978 bayesResult->Divide(bayesResult, mc, 1, 1, "");
979
980 chi2Result->SetTitle(";true multiplicity;unfolded measured/MC");
981 chi2Result->GetXaxis()->SetRangeUser(0, 150);
982 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
983 chi2Result->GetYaxis()->SetTitleOffset(1.2);
984 chi2Result->SetLineColor(1);
985 chi2Result->SetStats(kFALSE);
986
987 bayesResult->SetStats(kFALSE);
988 bayesResult->SetLineColor(2);
989
990 chi2Result->DrawCopy("HIST");
991 bayesResult->DrawCopy("SAME HIST");
992
993 TLine* line = new TLine(0, 1, 150, 1);
994 line->Draw();
995 }
996
997 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
998}
999
1000void StartingConditions()
1001{
1002 // data generated by runMultiplicitySelector.C StartingConditions
1003
1004 const char* name = "StartingConditions";
1005
1006 TFile* file = TFile::Open(Form("%s.root", name));
1007
1008 TCanvas* canvas = new TCanvas(name, name, 800, 400);
1009 canvas->Divide(2, 1);
1010
1011 TH1* mc = (TH1*) file->Get("mc");
1012 mc->Sumw2();
1013 mc->Scale(1.0 / mc->Integral());
1014
6d81c2de 1015 //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
0b4bfd98 1016
1017 for (Int_t i=0; i<6; ++i)
1018 {
1019 Int_t id = i;
1020 if (id > 2)
1021 id += 2;
1022
1023 TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
1024 TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
1025
1026 chi2Result->Divide(chi2Result, mc, 1, 1, "");
1027 bayesResult->Divide(bayesResult, mc, 1, 1, "");
1028
1029 chi2Result->SetTitle("a) #chi^{2} minimization;true multiplicity;unfolded / MC");
1030 chi2Result->GetXaxis()->SetRangeUser(0, 150);
1031 chi2Result->GetYaxis()->SetRangeUser(0.8, 1.2);
1032 chi2Result->GetYaxis()->SetTitleOffset(1.5);
1033 //chi2Result->SetMarkerStyle(marker[i]);
1034 chi2Result->SetLineColor(i+1);
1035 chi2Result->SetStats(kFALSE);
1036
1037 bayesResult->SetTitle("b) Bayesian method;true multiplicity;unfolded / MC");
1038 bayesResult->GetXaxis()->SetRangeUser(0, 150);
1039 bayesResult->GetYaxis()->SetRangeUser(0.8, 1.2);
1040 bayesResult->GetYaxis()->SetTitleOffset(1.5);
1041 bayesResult->SetStats(kFALSE);
1042 bayesResult->SetLineColor(2);
1043 bayesResult->SetLineColor(i+1);
1044
1045 canvas->cd(1);
1046 gPad->SetLeftMargin(0.12);
1047 chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1048
1049 canvas->cd(2);
1050 gPad->SetLeftMargin(0.12);
1051 bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1052
1053 //TLine* line = new TLine(0, 1, 150, 1);
1054 //line->Draw();
1055 }
1056
1057 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1058}
1059
1060void StatisticsPlot()
1061{
1062 const char* name = "StatisticsPlot";
1063
1064 TFile* file = TFile::Open(Form("%s.root", name));
1065
1066 TCanvas* canvas = new TCanvas(name, name, 600, 400);
1067
6d81c2de 1068 TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
0b4bfd98 1069 fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1070 fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1071 fitResultsChi2->Draw("AP");
1072
1073 TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1074 fitResultsChi2->Fit(f);
1075
1076 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1077
1078 TH1* mc[5];
1079 TH1* result[5];
1080
1081 const char* plotname = "chi2Result";
1082
6d81c2de 1083 name = "StatisticsPlotRatios";
1084 canvas = new TCanvas(name, name, 600, 400);
0b4bfd98 1085
1086 for (Int_t i=0; i<5; ++i)
1087 {
1088 mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1089 result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1090
1091 result[i]->SetLineColor(i+1);
1092 result[i]->Draw(((i == 0) ? "" : "SAME"));
1093 }
1094}
1095
1096void SystematicLowEfficiency()
1097{
1098 gSystem->Load("libPWG0base");
1099
1100 TFile::Open(correctionFile);
1101 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1102 mult->LoadHistograms("Multiplicity");
1103
1104 // calculate result with systematic effect
1105 TFile::Open("multiplicityMC_100k_1_loweff.root");
1106 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1107 mult2->LoadHistograms("Multiplicity");
1108
1109 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1110 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1111 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1112
1113 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
6d81c2de 1114 TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
0b4bfd98 1115
1116 DrawResultRatio(mcHist, result1, "SystematicLowEfficiencyLow.eps");
1117
1118 // calculate normal result
1119 TFile::Open("multiplicityMC_100k_1.root");
1120 mult2->LoadHistograms("Multiplicity");
1121
1122 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1123 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1124
6d81c2de 1125 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
0b4bfd98 1126
1127 DrawResultRatio(mcHist, result2, "SystematicLowEfficiencyOK.eps");
1128
1129 Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
1130}
1131
1132void SystematicMisalignment()
1133{
1134 gSystem->Load("libPWG0base");
1135
1136 TFile::Open(correctionFile);
1137 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1138 mult->LoadHistograms("Multiplicity");
1139
1140 TFile::Open("multiplicityMC_100k_fullmis.root");
1141 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1142 mult2->LoadHistograms("Multiplicity");
1143
1144 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1145 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1146 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1147
1148 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1149 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1150
1151 DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1152}
1153
6d81c2de 1154void SystematicMisalignmentTPC()
0b4bfd98 1155{
1156 gSystem->Load("libPWG0base");
1157
6d81c2de 1158 SetTPC();
0b4bfd98 1159
6d81c2de 1160 TFile::Open(correctionFile);
1161 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1162 mult->LoadHistograms("Multiplicity");
0b4bfd98 1163
6d81c2de 1164 TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
1165 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1166 mult2->LoadHistograms("Multiplicity");
1167
1168 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1169 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1170 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1171
1172 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1173 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1174
1175 DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
1176}
1177
1178void EfficiencySpecies()
1179{
1180 gSystem->Load("libPWG0base");
0b4bfd98 1181
1182 Int_t marker[] = {24, 25, 26};
1183 Int_t color[] = {1, 2, 4};
1184
6d81c2de 1185 // SPD TPC
1186 const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
1187 Float_t etaRange[] = {0.49, 0.9};
1188 const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1189
1190 TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1191 canvas->Divide(2, 1);
0b4bfd98 1192
6d81c2de 1193 for (Int_t loop=0; loop<2; ++loop)
0b4bfd98 1194 {
6d81c2de 1195 Printf("%s", fileName[loop]);
0b4bfd98 1196
6d81c2de 1197 TFile::Open(fileName[loop]);
1198 AliCorrection* correction[4];
0b4bfd98 1199
6d81c2de 1200 canvas->cd(loop+1);
0b4bfd98 1201
6d81c2de 1202 gPad->SetGridx();
1203 gPad->SetGridy();
1204 gPad->SetRightMargin(0.05);
1205 //gPad->SetTopMargin(0.05);
0b4bfd98 1206
6d81c2de 1207 TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
1208 legend->SetFillColor(0);
1209 legend->SetEntrySeparation(0.2);
1210
1211 Float_t below = 0;
1212 Float_t total = 0;
0b4bfd98 1213
6d81c2de 1214 for (Int_t i=0; i<3; ++i)
1215 {
1216 Printf("correction %d", i);
0b4bfd98 1217
6d81c2de 1218 TString name; name.Form("correction_%d", i);
1219 correction[i] = new AliCorrection(name, name);
1220 correction[i]->LoadHistograms();
0b4bfd98 1221
6d81c2de 1222 TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1223 TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
0b4bfd98 1224
6d81c2de 1225 // limit vtx axis
1226 gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
1227 meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
0b4bfd98 1228
6d81c2de 1229 // 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)
1230 /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1231 for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1232 {
1233 gene->SetBinContent(x, 0, z, 0);
1234 gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1235 meas->SetBinContent(x, 0, z, 0);
1236 meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1237 }*/
0b4bfd98 1238
6d81c2de 1239 // limit eta axis
1240 gene->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1241 meas->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
0b4bfd98 1242
6d81c2de 1243 TH1* genePt = gene->Project3D(Form("z_%d", i));
1244 TH1* measPt = meas->Project3D(Form("z_%d", i));
0b4bfd98 1245
6d81c2de 1246 genePt->Sumw2();
1247 measPt->Sumw2();
0b4bfd98 1248
6d81c2de 1249 TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1250 effPt->Reset();
1251 effPt->Divide(measPt, genePt, 1, 1, "B");
0b4bfd98 1252
6d81c2de 1253 Int_t bin = 0;
1254 for (bin=20; bin>=1; bin--)
1255 {
1256 if (effPt->GetBinContent(bin) < 0.5)
1257 break;
1258 }
0b4bfd98 1259
6d81c2de 1260 Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
0b4bfd98 1261
6d81c2de 1262 Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1263 Printf("%.4f of the particles are below that momentum", fraction);
0b4bfd98 1264
6d81c2de 1265 below += genePt->Integral(1, bin);
1266 total += genePt->Integral();
0b4bfd98 1267
6d81c2de 1268 effPt->SetLineColor(color[i]);
1269 effPt->SetMarkerColor(color[i]);
1270 effPt->SetMarkerStyle(marker[i]);
0b4bfd98 1271
6d81c2de 1272 effPt->GetXaxis()->SetRangeUser(0.06, 1);
1273 effPt->GetYaxis()->SetRangeUser(0, 1);
1274
1275 effPt->GetYaxis()->SetTitleOffset(1.2);
1276
1277 effPt->SetStats(kFALSE);
1278 effPt->SetTitle(titles[loop]);
1279 effPt->GetYaxis()->SetTitle("Efficiency");
1280
1281 effPt->DrawCopy((i == 0) ? "" : "SAME");
1282
1283 legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")));
1284 }
1285
1286 Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1287
1288 legend->Draw();
1289 }
0b4bfd98 1290
1291 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1292}
1293
6d81c2de 1294void ParticleSpeciesComparison1(const char* fileNameMC = "multiplicityMC_400k_syst_species.root", const char* fileNameESD = "multiplicityMC_100k_syst.root")
0b4bfd98 1295{
1296 gSystem->Load("libPWG0base");
1297
6d81c2de 1298 Bool_t chi2 = 1;
0b4bfd98 1299
1300 TFile::Open(fileNameESD);
1301 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
1302 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", etaRange));
1303
1304 TH1* results[10];
1305
1306 // loop over cases (normal, enhanced/reduced ratios)
1307 Int_t nMax = 7;
1308 for (Int_t i = 0; i<nMax; ++i)
1309 {
1310 TString folder;
1311 folder.Form("Multiplicity_%d", i);
1312
1313 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
1314
1315 TFile::Open(fileNameMC);
1316 mult->LoadHistograms();
1317
1318 mult->SetMultiplicityESD(etaRange, hist);
1319
1320 if (chi2)
1321 {
1322 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1323 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1324 //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1325 }
1326 else
1327 {
1328 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1329 //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1330 }
1331
1332 //Float_t averageRatio = 0;
1333 //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1334
1335 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1336
1337 //Printf("Case %d. Average ratio is %f", i, averageRatio);
1338 }
1339
1340 DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[0], "ParticleSpeciesComparison1_1.eps");
1341
1342 TH1* mc = hist2->ProjectionY("mymchist2", -1, -1, "e");
1343
1344 for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1345 {
1346 results[0]->SetBinError(i, 0);
1347 mc->SetBinError(i, 0);
1348 }
1349
6d81c2de 1350 const char* legendStrings[] = { "#pi^{#pm}", 0, "K^{#pm}", 0, "p,#bar{p}", 0 };
1351
1352 DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
0b4bfd98 1353
6d81c2de 1354 //not valid: draw chi2 uncertainty on top!
1355 /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
0b4bfd98 1356 TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1357 errorHist->SetLineColor(1);
1358 errorHist->SetLineWidth(2);
1359 TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1360 for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1361 {
1362 errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1363 errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1364 }
1365
1366 errorHist->DrawCopy("SAME");
6d81c2de 1367 errorHist2->DrawCopy("SAME");*/
0b4bfd98 1368
6d81c2de 1369 //canvas->SaveAs(canvas->GetName());
0b4bfd98 1370
6d81c2de 1371 DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
0b4bfd98 1372
6d81c2de 1373 //errorHist->DrawCopy("SAME");
1374 //errorHist2->DrawCopy("SAME");
0b4bfd98 1375
6d81c2de 1376 //canvas2->SaveAs(canvas2->GetName());
0b4bfd98 1377}
1378
1379void ParticleSpeciesComparison2()
1380{
1381 gSystem->Load("libPWG0base");
1382
1383 const char* fileNameMC = "multiplicityMC_400k_syst.root";
1384 const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1385 Bool_t chi2 = 0;
1386
1387 TFile::Open(fileNameMC);
1388 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1389 mult->LoadHistograms();
1390
1391 TH1* mc[10];
1392 TH1* results[10];
1393
1394 // loop over cases (normal, enhanced/reduced ratios)
1395 Int_t nMax = 7;
1396 for (Int_t i = 0; i<nMax; ++i)
1397 {
1398 TString folder;
1399 folder.Form("Multiplicity_%d", i);
1400
1401 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1402
1403 TFile::Open(fileNameESD);
1404 mult2->LoadHistograms();
1405
1406 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1407
1408 if (chi2)
1409 {
1410 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1411 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1412 //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1413 }
1414 else
1415 {
1416 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1417 //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1418 }
1419
1420 //Float_t averageRatio = 0;
1421 //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1422
1423 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1424
1425 TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1426 mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1427
1428 //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1429 //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1430
1431 //Printf("Case %d. Average ratio is %f", i, averageRatio);
1432 }
1433
1434 DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1435}
1436
1437TH1* Invert(TH1* eff)
1438{
1439 // calculate corr = 1 / eff
1440
1441 TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1442 corr->Reset();
1443
1444 for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1445 {
1446 if (eff->GetBinContent(i) > 0)
1447 {
1448 corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1449 corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1450 }
1451 }
1452
1453 return corr;
1454}
1455
1456void TriggerVertexCorrection()
1457{
1458 //
1459 // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1460 //
1461
1462 gSystem->Load("libPWG0base");
1463
1464 TFile::Open(correctionFile);
1465 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1466 mult->LoadHistograms("Multiplicity");
1467
1468 TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1469 TH1* corrMB = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1470
1471 TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 600);
1472
1473 corrINEL->SetStats(kFALSE);
1474 corrINEL->GetXaxis()->SetRangeUser(0, 30);
1475 corrINEL->GetYaxis()->SetRangeUser(0, 5);
1476 corrINEL->SetTitle(";true multiplicity;correction factor");
1477 corrINEL->SetMarkerStyle(22);
1478 corrINEL->Draw("PE");
1479
1480 corrMB->SetStats(kFALSE);
1481 corrMB->SetLineColor(2);
1482 corrMB->SetMarkerStyle(25);
1483 corrMB->SetMarkerColor(2);
1484 corrMB->Draw("SAME PE");
1485
1486 TLegend* legend = new TLegend(0.3, 0.5, 0.85, 0.65);
1487 legend->SetFillColor(0);
1488 legend->AddEntry(corrINEL, "correction to inelastic sample");
1489 legend->AddEntry(corrMB, "correction to minimum bias sample");
1490
1491 legend->Draw();
1492
1493 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1494}
1495
6d81c2de 1496void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
0b4bfd98 1497{
1498 gSystem->Load("libPWG0base");
1499
1500 TFile::Open(correctionFile);
1501 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1502 mult->LoadHistograms("Multiplicity");
1503
1504 TFile::Open(measuredFile);
1505 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1506 mult2->LoadHistograms("Multiplicity");
1507
1508 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1509
1510 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1511
6d81c2de 1512 TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1513
1514 TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1515 TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
0b4bfd98 1516
1517 if (!mc)
1518 {
1519 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
6d81c2de 1520 DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
0b4bfd98 1521 }
1522
6d81c2de 1523 TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
0b4bfd98 1524 canvas->SetGridx();
1525 canvas->SetGridy();
1526 canvas->SetRightMargin(0.05);
1527 canvas->SetTopMargin(0.05);
1528
1529 errorResponse->SetLineColor(1);
1530 errorResponse->GetXaxis()->SetRangeUser(0, 200);
1531 errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1532 errorResponse->SetStats(kFALSE);
1533 errorResponse->SetTitle(";true multiplicity;Uncertainty");
1534
1535 errorResponse->Draw();
1536
1537 errorMeasured->SetLineColor(2);
1538 errorMeasured->Draw("SAME");
1539
6d81c2de 1540 errorBoth->SetLineColor(4);
0b4bfd98 1541 errorBoth->Draw("SAME");
1542
1543 Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
1544 Printf("Average errorMeasured: %f", errorMeasured->Integral(2, 150) / 149);
1545 Printf("Average errorBoth: %f", errorBoth->Integral(2, 150) / 149);
1546
1547 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1548
1549 TFile* file = new TFile(Form("%s.root", canvas->GetName()), "RECREATE");
1550 errorResponse->Write();
1551 errorMeasured->Write();
1552 errorBoth->Write();
1553 file->Close();
1554}
1555
6d81c2de 1556void StatisticalUncertaintyCompare(const char* det = "SPD")
1557{
1558 TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1559 TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1560 TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1561 TH1* errorBoth = (TH1*) file1->Get("errorBoth");
1562
1563 TString str;
1564 str.Form("StatisticalUncertaintyCompare%s", det);
1565
1566 TCanvas* canvas = new TCanvas(str, str, 600, 400);
1567 canvas->SetGridx();
1568 canvas->SetGridy();
1569 canvas->SetRightMargin(0.05);
1570 canvas->SetTopMargin(0.05);
1571
1572 errorResponse->SetLineColor(1);
1573 errorResponse->GetXaxis()->SetRangeUser(1, (strcmp(det, "TPC") ? 200 : 100));
1574 errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1575 errorResponse->SetStats(kFALSE);
1576 errorResponse->SetTitle(";true multiplicity;Uncertainty");
1577
1578 errorResponse->Draw();
1579
1580 errorMeasured->SetLineColor(2);
1581 errorMeasured->Draw("SAME");
1582
1583 errorBoth->SetLineColor(4);
1584 errorBoth->Draw("SAME");
1585
1586 TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
1587 TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
1588
1589 errorBoth2->SetLineColor(4);
1590 errorBoth2->SetLineStyle(2);
1591 errorBoth2->Draw("SAME");
1592
1593 TLegend* legend = new TLegend(0.2, 0.6, 0.6, 0.9);
1594 legend->SetFillColor(0);
1595 legend->AddEntry(errorResponse, "response matrix (Bayesian)");
1596 legend->AddEntry(errorMeasured, "measured (Bayesian)");
1597 legend->AddEntry(errorBoth, "both (Bayesian)");
1598 legend->AddEntry(errorBoth2, "both (#chi^{2} minimization)");
1599 legend->Draw();
1600
1601 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1602}
1603
0b4bfd98 1604void EfficiencyComparison(Int_t eventType = 2)
1605{
1606 const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
1607
1608 gSystem->Load("libPWG0base");
1609
1610 TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
1611 canvas->SetGridx();
1612 canvas->SetGridy();
1613 canvas->SetRightMargin(0.05);
1614 canvas->SetTopMargin(0.05);
1615
1616 AliMultiplicityCorrection* data[4];
1617 TH1* effArray[4];
1618
1619 Int_t markers[] = { 24, 25, 26, 5 };
1620 Int_t colors[] = { 1, 2, 3, 4 };
1621
1622 TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
1623 legend->SetFillColor(0);
1624
1625 TH1* effError = 0;
1626
1627 for (Int_t i=0; i<4; ++i)
1628 {
1629 TString name;
1630 name.Form("Multiplicity_%d", i);
1631
1632 TFile::Open(files[i]);
1633 data[i] = new AliMultiplicityCorrection(name, name);
1634
1635 if (i < 3)
1636 {
1637 data[i]->LoadHistograms("Multiplicity");
1638 }
1639 else
1640 data[i]->LoadHistograms("Multiplicity_0");
1641
6d81c2de 1642 TH1* eff = (TH1*) data[i]->GetEfficiency(3, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
0b4bfd98 1643 effArray[i] = eff;
1644
1645 eff->GetXaxis()->SetRangeUser(0, 15);
1646 eff->GetYaxis()->SetRangeUser(0, 1.1);
1647 eff->SetStats(kFALSE);
1648 eff->SetTitle(";true multiplicity;Efficiency");
1649 eff->SetLineColor(colors[i]);
1650 eff->SetMarkerColor(colors[i]);
1651 eff->SetMarkerStyle(markers[i]);
1652
1653 if (i == 3)
1654 {
1655 for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1656 eff->SetBinError(bin, 0);
1657
1658 // loop over cross section combinations
1659 for (Int_t j=1; j<7; ++j)
1660 {
1661 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
1662 mult->LoadHistograms(Form("Multiplicity_%d", j));
1663
6d81c2de 1664 TH1* eff2 = mult->GetEfficiency(3, (AliMultiplicityCorrection::EventType) eventType);
0b4bfd98 1665
1666 for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1667 {
1668 // TODO we could also do assymetric errors here
1669 Float_t deviation = 10.0 * TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
1670
6d81c2de 1671 eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), (Double_t) deviation));
0b4bfd98 1672 }
1673 }
1674
1675 for (Int_t bin=1; bin<=20; bin++)
1676 if (eff->GetBinContent(bin) > 0)
1677 Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1678
6d81c2de 1679 effError = (TH1*) eff->Clone("effError");
0b4bfd98 1680 effError->Reset();
1681
1682 for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
1683 if (eff->GetBinContent(bin) > 0)
1684 effError->SetBinContent(bin, eff->GetBinError(bin) / eff->GetBinContent(bin));
1685
1686 effError->DrawCopy("SAME HIST");
1687 }
1688
1689 eff->SetBinContent(1, 0);
1690 eff->SetBinError(1, 0);
1691
1692 canvas->cd();
1693 if (i == 0)
1694 {
1695 eff->DrawCopy("P");
1696 }
1697 else
1698 eff->DrawCopy("SAME P");
1699
1700 legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
1701 }
1702
1703 legend->AddEntry(effError, "relative syst. uncertainty #times 10");
1704
1705 legend->Draw();
1706
1707 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1708}
1709
1710void ModelDependencyPlot()
1711{
1712 gSystem->Load("libPWG0base");
1713
1714 TFile::Open("multiplicityMC_3M.root");
1715 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1716 mult->LoadHistograms("Multiplicity");
1717
1718 TH2* proj = (TH2*) mult->GetCorrelation(3)->Project3D("zy");
1719
1720 TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 800, 400);
1721 canvas->SetGridx();
1722 canvas->SetGridy();
1723 //canvas->SetRightMargin(0.05);
1724 //canvas->SetTopMargin(0.05);
1725
1726 canvas->Divide(2, 1);
1727
1728 canvas->cd(2);
1729 gPad->SetLogy();
1730
1731 Int_t selectedMult = 30;
1732 Int_t yMax = 200000;
1733
1734 TH1* full = proj->ProjectionX("full");
1735 TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
1736
1737 full->SetStats(kFALSE);
1738 full->GetXaxis()->SetRangeUser(0, 200);
1739 full->GetYaxis()->SetRangeUser(5, yMax);
1740 full->SetTitle(";multiplicity");
1741
1742 selected->SetLineColor(0);
1743 selected->SetMarkerColor(2);
1744 selected->SetMarkerStyle(7);
1745
1746 full->Draw();
1747 selected->Draw("SAME P");
1748
1749 TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
1750 legend->SetFillColor(0);
1751 legend->AddEntry(full, "true");
1752 legend->AddEntry(selected, "measured");
1753 legend->Draw();
1754
1755 TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
1756 line->SetLineWidth(2);
1757 line->Draw();
1758
1759 canvas->cd(1);
1760 gPad->SetLogy();
1761
6d81c2de 1762 full = proj->ProjectionY("full2");
1763 selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
0b4bfd98 1764
1765 full->SetStats(kFALSE);
1766 full->GetXaxis()->SetRangeUser(0, 200);
1767 full->GetYaxis()->SetRangeUser(5, yMax);
1768 full->SetTitle(";multiplicity");
1769
1770 full->SetLineColor(0);
1771 full->SetMarkerColor(2);
1772 full->SetMarkerStyle(7);
1773
1774 full->Draw("P");
1775 selected->Draw("SAME");
1776
1777 legend->Draw();
1778
6d81c2de 1779 line = new TLine(selectedMult, 5, selectedMult, yMax);
0b4bfd98 1780 line->SetLineWidth(2);
1781 line->Draw();
1782
1783 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1784}
1785
1786void SystematicpTSpectrum()
1787{
1788 gSystem->Load("libPWG0base");
1789
1790 TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
1791 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1792 mult->LoadHistograms("Multiplicity");
1793
1794 TFile::Open("multiplicityMC_100k_syst.root");
1795 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1796 mult2->LoadHistograms("Multiplicity");
1797
1798 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1799 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1800 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1801
1802 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1803 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1804
1805 DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
1806}
1807
1808// to be deleted
1809/*void covMatrix(Bool_t mc = kTRUE)
1810{
1811 gSystem->Load("libPWG0base");
1812
1813 TFile::Open(correctionFile);
1814 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1815 mult->LoadHistograms("Multiplicity");
1816
1817 TFile::Open(measuredFile);
1818 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1819 mult2->LoadHistograms("Multiplicity");
1820
1821 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1822
1823 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1824
1825 mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0));
1826}*/
1827
1828Double_t FitPtFunc(Double_t *x, Double_t *par)
1829{
1830 Double_t xx = x[0];
1831
1832 Float_t val1 = par[1] + par[2] * xx + par[3] * xx * xx;
1833 Float_t val2 = TMath::Exp(par[4] + par[5] * xx);
1834
1835 const Float_t kTransitionWidth = 0;
1836
1837 // power law part
1838 if (xx < par[0] - kTransitionWidth)
1839 {
1840 return val1;
1841 }
1842 /*else if (xx < par[0] + kTransitionWidth)
1843 {
1844 // smooth transition
1845 Float_t factor = (xx - par[0] + kTransitionWidth) / kTransitionWidth / 2;
1846 return (1 - factor) * val1 + factor * val2;
1847 }*/
1848 else
1849 {
1850 return val2;
1851 }
1852}
1853
1854void FitPt(const char* fileName = "firstplots100k_truept.root")
1855{
1856 gSystem->Load("libPWG0base");
1857
1858 TFile::Open(fileName);
1859
1860 /*
1861 // merge corrections
1862 AliCorrection* correction[4];
1863 TList list;
1864
1865 for (Int_t i=0; i<4; ++i)
1866 {
1867 Printf("correction %d", i);
1868
1869 TString name; name.Form("correction_%d", i);
1870 correction[i] = new AliCorrection(name, name);
1871 correction[i]->LoadHistograms();
1872
1873 if (i > 0)
1874 list.Add(correction[i]);
1875 }
1876
1877 correction[0]->Merge(&list);
1878
1879 TH3* gene = correction[0]->GetTrackCorrection()->GetGeneratedHistogram();
1880
1881 // limit vtx, eta axis
1882 gene->GetXaxis()->SetRangeUser(-5.9, 5.9);
1883 gene->GetYaxis()->SetRangeUser(-1.99, 0.99);
1884
1885 TH1* genePt = gene->Project3D("z");*/
1886 TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
1887 genePt->Sumw2();
1888
1889 //genePt->Scale(1.0 / genePt->Integral());
1890
1891 // normalize by bin width
1892 for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1893 genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1894
1895 /// genePt->GetXaxis()->GetBinCenter(x));
1896
1897 genePt->GetXaxis()->SetRangeUser(0, 7.9);
6d81c2de 1898 //genePt->GetYaxis()->SetTitle("a.u.");
0b4bfd98 1899
1900 TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
1901 //func->SetLineColor(2);
1902 func->SetParameters(1, -1, 1, 1, 1);
1903 func->SetParLimits(3, 1, 10);
1904 func->SetParLimits(4, 0, 10);
1905
1906 //TF1* func = new TF1("func", "[1]*x**[0]", 0.001, 100);
1907
1908 //TF1* func = new TF1("func", &FitPtFunc, 0, 2, 6);
1909 //func->SetParameters(0.3, -2.34909e-01, 1.54394e+01, -3.04134e+01, 1.41912e+00, -2.79284e+00);
1910 //func->FixParameter(0, 0.314);
1911 //func->SetParLimits(0, 0.1, 0.3);
1912
1913 genePt->SetMarkerStyle(25);
1914 genePt->SetTitle("");
1915 genePt->SetStats(kFALSE);
1916 genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
1917 //func->Draw("SAME");
1918
1919 // fit only exp. part
1920 func->SetParameters(1, -1);
1921 func->FixParameter(2, 0);
1922 func->FixParameter(3, 1);
1923 func->FixParameter(4, 1);
1924 genePt->Fit(func, "0", "", 0.2, 1);
1925
1926 new TCanvas;
1927 genePt->DrawCopy("P");
1928 func->SetRange(0.02, 8);
1929 func->DrawCopy("SAME");
1930 gPad->SetLogy();
1931
1932 // now fix exp. parameters and fit second part
1933 Double_t param0 = func->GetParameter(0);
1934 Double_t param1 = func->GetParameter(1);
1935 func->SetParameters(0, -1, 1, 1, 1);
1936 func->FixParameter(0, 0);
1937 func->FixParameter(1, -1);
1938 func->ReleaseParameter(2);
1939 func->ReleaseParameter(3);
1940 func->ReleaseParameter(4);
1941 func->SetParLimits(3, 1, 10);
1942 func->SetParLimits(4, 0, 10);
1943
1944 genePt->Fit(func, "0", "", 1.5, 4);
1945
1946 new TCanvas;
1947 genePt->DrawCopy("P");
1948 func->SetRange(0.02, 8);
1949 func->DrawCopy("SAME");
1950 gPad->SetLogy();
1951
1952 // fit both
1953 func->SetParameter(0, param0);
1954 func->SetParameter(1, param1);
1955 func->ReleaseParameter(0);
1956 func->ReleaseParameter(1);
1957
1958 new TCanvas;
1959 genePt->DrawCopy("P");
6d81c2de 1960 func->SetRange(0.02, 5);
0b4bfd98 1961 func->DrawCopy("SAME");
1962 gPad->SetLogy();
1963
1964 genePt->Fit(func, "0", "", 0.2, 4);
1965
6d81c2de 1966 TCanvas* canvas = new TCanvas("FitPt", "FitPt", 800, 400);
1967 canvas->Divide(2, 1);
1968 canvas->cd(1);
1969
1970 gPad->SetGridx();
1971 gPad->SetGridy();
1972 gPad->SetLeftMargin(0.13);
1973 gPad->SetRightMargin(0.05);
1974 gPad->SetTopMargin(0.05);
0b4bfd98 1975
6d81c2de 1976 genePt->GetXaxis()->SetRangeUser(0, 4.9);
1977 genePt->GetYaxis()->SetRangeUser(1e-2, 1e4);
1978 genePt->GetYaxis()->SetTitleOffset(1.4);
1979 genePt->GetXaxis()->SetTitleOffset(1.1);
0b4bfd98 1980 genePt->DrawCopy("P");
6d81c2de 1981 func->SetRange(0.02, 5);
1982 func->DrawCopy("SAME");
1983 gPad->SetLogy();
1984
1985 canvas->cd(2);
1986
1987 TH1* genePtClone = (TH1*) genePt->Clone("genePtClone");
1988 genePtClone->Reset();
1989 genePtClone->DrawCopy("P");
1990
1991 gPad->SetGridx();
1992 gPad->SetGridy();
1993 gPad->SetLeftMargin(0.13);
1994 gPad->SetRightMargin(0.05);
1995 gPad->SetTopMargin(0.05);
1996
0b4bfd98 1997 func->DrawCopy("SAME");
1998 gPad->SetLogy();
1999
6d81c2de 2000 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2001
0b4bfd98 2002 TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2003
2004 TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2005
2006 TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2007
2008 for (Int_t param=0; param<5; param++)
2009 {
2010 for (Int_t sign=0; sign<2; sign++)
2011 {
2012 TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign)); // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2013 func2->SetParameters(func->GetParameters());
2014 //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2015
2016 Float_t factor = ((sign == 0) ? 0.9 : 1.1);
2017 func2->SetParameter(param, func2->GetParameter(param) * factor);
2018 //func2->Print();
2019
6d81c2de 2020 canvas->cd(2);
0b4bfd98 2021 func2->SetLineWidth(1);
6d81c2de 2022 func2->SetLineColor(2);
0b4bfd98 2023 func2->DrawCopy("SAME");
2024
2025 canvas2->cd();
2026 TH1* second = func2->GetHistogram();
2027 second->Divide(first);
2028 second->SetLineColor(param + 1);
2029 second->GetYaxis()->SetRangeUser(0, 2);
2030 second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2031 second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2032 }
2033 }
2034
2035 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
6d81c2de 2036 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
0b4bfd98 2037
2038 file->Close();
2039}
2040
6d81c2de 2041void DrawSystematicpT()
2042{
2043 TFile* file = TFile::Open("SystematicpT.root");
2044
2045 TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2046 TH1* result2 = (TH1*) file->Get("result_unity");
2047
2048 TH1* mcHist[12];
2049 TH1* result[12];
2050
2051 Int_t nParams = 5;
2052
2053 for (Int_t id=0; id<nParams*2; ++id)
2054 {
2055 mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2056 result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2057 }
2058
2059 DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2060
2061 //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2062
2063 DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2064
2065 //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2066
2067 // does not make sense: mc is different
2068 //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2069}
2070
2071void SystematicpT(Bool_t chi2 = 1)
0b4bfd98 2072{
2073 gSystem->Load("libPWG0base");
2074
6d81c2de 2075 TFile::Open("ptspectrum900.root");
0b4bfd98 2076 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2077 mult->LoadHistograms("Multiplicity");
2078
2079 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2080
2081 TH1* mcHist[12];
2082 TH1* result[12];
2083
6d81c2de 2084 Int_t nParams = 5;
0b4bfd98 2085
2086 for (Int_t param=0; param<nParams; param++)
2087 {
2088 for (Int_t sign=0; sign<2; sign++)
2089 {
2090 // calculate result with systematic effect
6d81c2de 2091 TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
0b4bfd98 2092 mult2->LoadHistograms("Multiplicity");
2093
2094 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2095
6d81c2de 2096 if (chi2)
2097 {
2098 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2099 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2100 }
2101 else
2102 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
0b4bfd98 2103
2104 Int_t id = param * 2 + sign;
2105
2106 mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2107 result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2108
2109 TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2110 DrawResultRatio(mcHist[id], result[id], tmp);
2111 }
2112 }
2113
2114 // calculate normal result
2115 TFile::Open("ptspectrum100_1.root");
2116 mult2->LoadHistograms("Multiplicity");
6d81c2de 2117 TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
0b4bfd98 2118
2119 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2120
6d81c2de 2121 if (chi2)
2122 {
2123 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2124 }
2125 else
2126 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
0b4bfd98 2127
6d81c2de 2128 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
0b4bfd98 2129
6d81c2de 2130 TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2131 mcHist2->Write();
2132 result2->Write();
2133 for (Int_t id=0; id<nParams*2; ++id)
2134 {
2135 mcHist[id]->Write();
2136 result[id]->Write();
2137 }
2138 file->Close();
2139
2140 DrawSystematicpT();
2141}
0b4bfd98 2142
6d81c2de 2143void SystematicpTCutOff()
2144{
2145 gSystem->Load("libPWG0base");
2146 SetTPC();
0b4bfd98 2147
6d81c2de 2148 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2149 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2150 mult->LoadHistograms("Multiplicity");
0b4bfd98 2151
6d81c2de 2152 TFile::Open(measuredFile);
2153 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2154 mult2->LoadHistograms("Multiplicity");
0b4bfd98 2155
6d81c2de 2156 // "normal" result
2157 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2158 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2159 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
0b4bfd98 2160
6d81c2de 2161 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2162 TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
0b4bfd98 2163
6d81c2de 2164 // change of pt spectrum (down)
2165 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2166 AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2167 mult3->LoadHistograms("Multiplicity");
0b4bfd98 2168
6d81c2de 2169 mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2170 mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2171 mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2172
2173 TH1* result2 = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2174
2175 Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
0b4bfd98 2176}