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