]> 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
adding macro that applies the corrections and visualizes the result
[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
6d81c2de 1228 for (Int_t loop=0; loop<2; ++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
6d81c2de 1255 for (Int_t i=0; i<3; ++i)
1256 {
1257 Printf("correction %d", i);
0b4bfd98 1258
6d81c2de 1259 TString name; name.Form("correction_%d", i);
1260 correction[i] = new AliCorrection(name, name);
1261 correction[i]->LoadHistograms();
0b4bfd98 1262
6d81c2de 1263 TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1264 TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
0b4bfd98 1265
6d81c2de 1266 // limit vtx axis
1267 gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
1268 meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
0b4bfd98 1269
6d81c2de 1270 // 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)
1271 /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1272 for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1273 {
1274 gene->SetBinContent(x, 0, z, 0);
1275 gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1276 meas->SetBinContent(x, 0, z, 0);
1277 meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1278 }*/
0b4bfd98 1279
6d81c2de 1280 // limit eta axis
1281 gene->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1282 meas->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
0b4bfd98 1283
6d81c2de 1284 TH1* genePt = gene->Project3D(Form("z_%d", i));
1285 TH1* measPt = meas->Project3D(Form("z_%d", i));
0b4bfd98 1286
6d81c2de 1287 genePt->Sumw2();
1288 measPt->Sumw2();
0b4bfd98 1289
6d81c2de 1290 TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1291 effPt->Reset();
1292 effPt->Divide(measPt, genePt, 1, 1, "B");
0b4bfd98 1293
6d81c2de 1294 Int_t bin = 0;
1295 for (bin=20; bin>=1; bin--)
1296 {
1297 if (effPt->GetBinContent(bin) < 0.5)
1298 break;
1299 }
0b4bfd98 1300
6d81c2de 1301 Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
0b4bfd98 1302
6d81c2de 1303 Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1304 Printf("%.4f of the particles are below that momentum", fraction);
0b4bfd98 1305
6d81c2de 1306 below += genePt->Integral(1, bin);
1307 total += genePt->Integral();
0b4bfd98 1308
6d81c2de 1309 effPt->SetLineColor(color[i]);
1310 effPt->SetMarkerColor(color[i]);
1311 effPt->SetMarkerStyle(marker[i]);
0b4bfd98 1312
6d81c2de 1313 effPt->GetXaxis()->SetRangeUser(0.06, 1);
1314 effPt->GetYaxis()->SetRangeUser(0, 1);
1315
1316 effPt->GetYaxis()->SetTitleOffset(1.2);
1317
1318 effPt->SetStats(kFALSE);
1319 effPt->SetTitle(titles[loop]);
1320 effPt->GetYaxis()->SetTitle("Efficiency");
1321
1322 effPt->DrawCopy((i == 0) ? "" : "SAME");
1323
1324 legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")));
1325 }
1326
1327 Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1328
1329 legend->Draw();
1330 }
0b4bfd98 1331
1332 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1333}
1334
44466df2 1335void ParticleSpeciesComparison1(Bool_t chi2 = kTRUE, const char* fileNameMC = "multiplicityMC_400k_syst_species.root", const char* fileNameESD = "multiplicityMC_100k_syst.root")
0b4bfd98 1336{
1337 gSystem->Load("libPWG0base");
1338
0b4bfd98 1339 TFile::Open(fileNameESD);
1340 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
1341 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", etaRange));
1342
1343 TH1* results[10];
1344
1345 // loop over cases (normal, enhanced/reduced ratios)
1346 Int_t nMax = 7;
1347 for (Int_t i = 0; i<nMax; ++i)
1348 {
1349 TString folder;
1350 folder.Form("Multiplicity_%d", i);
1351
1352 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
1353
1354 TFile::Open(fileNameMC);
1355 mult->LoadHistograms();
1356
1357 mult->SetMultiplicityESD(etaRange, hist);
1358
1359 if (chi2)
1360 {
1361 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1362 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1363 //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1364 }
1365 else
1366 {
1367 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1368 //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1369 }
1370
1371 //Float_t averageRatio = 0;
1372 //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1373
1374 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1375
1376 //Printf("Case %d. Average ratio is %f", i, averageRatio);
1377 }
1378
1379 DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[0], "ParticleSpeciesComparison1_1.eps");
1380
1381 TH1* mc = hist2->ProjectionY("mymchist2", -1, -1, "e");
1382
1383 for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1384 {
1385 results[0]->SetBinError(i, 0);
1386 mc->SetBinError(i, 0);
1387 }
1388
6d81c2de 1389 const char* legendStrings[] = { "#pi^{#pm}", 0, "K^{#pm}", 0, "p,#bar{p}", 0 };
1390
1391 DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
0b4bfd98 1392
6d81c2de 1393 //not valid: draw chi2 uncertainty on top!
1394 /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
0b4bfd98 1395 TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1396 errorHist->SetLineColor(1);
1397 errorHist->SetLineWidth(2);
1398 TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1399 for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1400 {
1401 errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1402 errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1403 }
1404
1405 errorHist->DrawCopy("SAME");
6d81c2de 1406 errorHist2->DrawCopy("SAME");*/
0b4bfd98 1407
6d81c2de 1408 //canvas->SaveAs(canvas->GetName());
0b4bfd98 1409
6d81c2de 1410 DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
0b4bfd98 1411
6d81c2de 1412 //errorHist->DrawCopy("SAME");
1413 //errorHist2->DrawCopy("SAME");
0b4bfd98 1414
6d81c2de 1415 //canvas2->SaveAs(canvas2->GetName());
0b4bfd98 1416}
1417
44466df2 1418/*void ParticleSpeciesComparison2()
0b4bfd98 1419{
1420 gSystem->Load("libPWG0base");
1421
1422 const char* fileNameMC = "multiplicityMC_400k_syst.root";
1423 const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1424 Bool_t chi2 = 0;
1425
1426 TFile::Open(fileNameMC);
1427 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1428 mult->LoadHistograms();
1429
1430 TH1* mc[10];
1431 TH1* results[10];
1432
1433 // loop over cases (normal, enhanced/reduced ratios)
1434 Int_t nMax = 7;
1435 for (Int_t i = 0; i<nMax; ++i)
1436 {
1437 TString folder;
1438 folder.Form("Multiplicity_%d", i);
1439
1440 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1441
1442 TFile::Open(fileNameESD);
1443 mult2->LoadHistograms();
1444
1445 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1446
1447 if (chi2)
1448 {
1449 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1450 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1451 //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1452 }
1453 else
1454 {
1455 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1456 //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1457 }
1458
1459 //Float_t averageRatio = 0;
1460 //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1461
1462 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1463
1464 TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1465 mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1466
1467 //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1468 //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1469
1470 //Printf("Case %d. Average ratio is %f", i, averageRatio);
1471 }
1472
1473 DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
44466df2 1474}*/
0b4bfd98 1475
1476TH1* Invert(TH1* eff)
1477{
1478 // calculate corr = 1 / eff
1479
1480 TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1481 corr->Reset();
1482
1483 for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1484 {
1485 if (eff->GetBinContent(i) > 0)
1486 {
1487 corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1488 corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1489 }
1490 }
1491
1492 return corr;
1493}
1494
1495void TriggerVertexCorrection()
1496{
1497 //
1498 // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1499 //
1500
1501 gSystem->Load("libPWG0base");
1502
1503 TFile::Open(correctionFile);
1504 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1505 mult->LoadHistograms("Multiplicity");
1506
1507 TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1508 TH1* corrMB = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1509
1510 TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 600);
1511
1512 corrINEL->SetStats(kFALSE);
1513 corrINEL->GetXaxis()->SetRangeUser(0, 30);
1514 corrINEL->GetYaxis()->SetRangeUser(0, 5);
1515 corrINEL->SetTitle(";true multiplicity;correction factor");
1516 corrINEL->SetMarkerStyle(22);
1517 corrINEL->Draw("PE");
1518
1519 corrMB->SetStats(kFALSE);
1520 corrMB->SetLineColor(2);
1521 corrMB->SetMarkerStyle(25);
1522 corrMB->SetMarkerColor(2);
1523 corrMB->Draw("SAME PE");
1524
1525 TLegend* legend = new TLegend(0.3, 0.5, 0.85, 0.65);
1526 legend->SetFillColor(0);
1527 legend->AddEntry(corrINEL, "correction to inelastic sample");
1528 legend->AddEntry(corrMB, "correction to minimum bias sample");
1529
1530 legend->Draw();
1531
1532 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1533}
1534
6d81c2de 1535void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
0b4bfd98 1536{
1537 gSystem->Load("libPWG0base");
1538
1539 TFile::Open(correctionFile);
1540 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1541 mult->LoadHistograms("Multiplicity");
1542
1543 TFile::Open(measuredFile);
1544 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1545 mult2->LoadHistograms("Multiplicity");
1546
1547 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1548
1549 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1550
6d81c2de 1551 TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1552
1553 TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1554 TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
0b4bfd98 1555
1556 if (!mc)
1557 {
1558 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
6d81c2de 1559 DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
0b4bfd98 1560 }
1561
6d81c2de 1562 TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
0b4bfd98 1563 canvas->SetGridx();
1564 canvas->SetGridy();
1565 canvas->SetRightMargin(0.05);
1566 canvas->SetTopMargin(0.05);
1567
1568 errorResponse->SetLineColor(1);
1569 errorResponse->GetXaxis()->SetRangeUser(0, 200);
1570 errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1571 errorResponse->SetStats(kFALSE);
1572 errorResponse->SetTitle(";true multiplicity;Uncertainty");
1573
1574 errorResponse->Draw();
1575
1576 errorMeasured->SetLineColor(2);
1577 errorMeasured->Draw("SAME");
1578
6d81c2de 1579 errorBoth->SetLineColor(4);
0b4bfd98 1580 errorBoth->Draw("SAME");
1581
1582 Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
1583 Printf("Average errorMeasured: %f", errorMeasured->Integral(2, 150) / 149);
1584 Printf("Average errorBoth: %f", errorBoth->Integral(2, 150) / 149);
1585
1586 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1587
1588 TFile* file = new TFile(Form("%s.root", canvas->GetName()), "RECREATE");
1589 errorResponse->Write();
1590 errorMeasured->Write();
1591 errorBoth->Write();
1592 file->Close();
1593}
1594
6d81c2de 1595void StatisticalUncertaintyCompare(const char* det = "SPD")
1596{
1597 TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1598 TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1599 TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1600 TH1* errorBoth = (TH1*) file1->Get("errorBoth");
1601
1602 TString str;
1603 str.Form("StatisticalUncertaintyCompare%s", det);
1604
1605 TCanvas* canvas = new TCanvas(str, str, 600, 400);
1606 canvas->SetGridx();
1607 canvas->SetGridy();
1608 canvas->SetRightMargin(0.05);
1609 canvas->SetTopMargin(0.05);
1610
1611 errorResponse->SetLineColor(1);
1612 errorResponse->GetXaxis()->SetRangeUser(1, (strcmp(det, "TPC") ? 200 : 100));
1613 errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1614 errorResponse->SetStats(kFALSE);
0f67a57c 1615 errorResponse->GetYaxis()->SetTitleOffset(1.2);
1616 errorResponse->SetTitle(";true multiplicity;#sigma(U-T)/T");
6d81c2de 1617
1618 errorResponse->Draw();
1619
1620 errorMeasured->SetLineColor(2);
1621 errorMeasured->Draw("SAME");
1622
1623 errorBoth->SetLineColor(4);
1624 errorBoth->Draw("SAME");
1625
1626 TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
1627 TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
1628
1629 errorBoth2->SetLineColor(4);
1630 errorBoth2->SetLineStyle(2);
1631 errorBoth2->Draw("SAME");
1632
1633 TLegend* legend = new TLegend(0.2, 0.6, 0.6, 0.9);
1634 legend->SetFillColor(0);
1635 legend->AddEntry(errorResponse, "response matrix (Bayesian)");
1636 legend->AddEntry(errorMeasured, "measured (Bayesian)");
1637 legend->AddEntry(errorBoth, "both (Bayesian)");
1638 legend->AddEntry(errorBoth2, "both (#chi^{2} minimization)");
1639 legend->Draw();
1640
1641 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1642}
1643
0f67a57c 1644void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
0b4bfd98 1645{
1646 const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
1647
1648 gSystem->Load("libPWG0base");
1649
1650 TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
1651 canvas->SetGridx();
1652 canvas->SetGridy();
1653 canvas->SetRightMargin(0.05);
1654 canvas->SetTopMargin(0.05);
1655
1656 AliMultiplicityCorrection* data[4];
1657 TH1* effArray[4];
1658
1659 Int_t markers[] = { 24, 25, 26, 5 };
1660 Int_t colors[] = { 1, 2, 3, 4 };
1661
1662 TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
1663 legend->SetFillColor(0);
1664
1665 TH1* effError = 0;
1666
1667 for (Int_t i=0; i<4; ++i)
1668 {
1669 TString name;
1670 name.Form("Multiplicity_%d", i);
1671
1672 TFile::Open(files[i]);
1673 data[i] = new AliMultiplicityCorrection(name, name);
1674
1675 if (i < 3)
1676 {
1677 data[i]->LoadHistograms("Multiplicity");
1678 }
1679 else
1680 data[i]->LoadHistograms("Multiplicity_0");
1681
44466df2 1682 TH1* eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
0b4bfd98 1683 effArray[i] = eff;
1684
1685 eff->GetXaxis()->SetRangeUser(0, 15);
1686 eff->GetYaxis()->SetRangeUser(0, 1.1);
1687 eff->SetStats(kFALSE);
1688 eff->SetTitle(";true multiplicity;Efficiency");
1689 eff->SetLineColor(colors[i]);
1690 eff->SetMarkerColor(colors[i]);
1691 eff->SetMarkerStyle(markers[i]);
1692
1693 if (i == 3)
1694 {
1695 for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1696 eff->SetBinError(bin, 0);
1697
1698 // loop over cross section combinations
1699 for (Int_t j=1; j<7; ++j)
1700 {
1701 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
1702 mult->LoadHistograms(Form("Multiplicity_%d", j));
1703
44466df2 1704 TH1* eff2 = mult->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType);
0b4bfd98 1705
1706 for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1707 {
44466df2 1708 // TODO we could also do asymmetric errors here
1709 Float_t deviation = TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
0b4bfd98 1710
6d81c2de 1711 eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), (Double_t) deviation));
0b4bfd98 1712 }
1713 }
1714
1715 for (Int_t bin=1; bin<=20; bin++)
1716 if (eff->GetBinContent(bin) > 0)
1717 Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
0f67a57c 1718
1719 if (uncertainty) {
1720 effError = (TH1*) eff->Clone("effError");
1721 effError->Reset();
1722
1723 for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
1724 if (eff->GetBinContent(bin) > 0)
1725 effError->SetBinContent(bin, 10.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1726
1727 effError->SetLineColor(1);
1728 effError->SetMarkerStyle(1);
1729 effError->DrawCopy("SAME HIST");
1730 }
0b4bfd98 1731 }
1732
1733 eff->SetBinContent(1, 0);
1734 eff->SetBinError(1, 0);
1735
1736 canvas->cd();
1737 if (i == 0)
1738 {
1739 eff->DrawCopy("P");
1740 }
1741 else
1742 eff->DrawCopy("SAME P");
1743
1744 legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
1745 }
1746
0f67a57c 1747 if (uncertainty)
1748 legend->AddEntry(effError, "relative syst. uncertainty #times 10");
0b4bfd98 1749
1750 legend->Draw();
1751
1752 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1753}
1754
1755void ModelDependencyPlot()
1756{
1757 gSystem->Load("libPWG0base");
1758
1759 TFile::Open("multiplicityMC_3M.root");
1760 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1761 mult->LoadHistograms("Multiplicity");
1762
1763 TH2* proj = (TH2*) mult->GetCorrelation(3)->Project3D("zy");
1764
1765 TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 800, 400);
1766 canvas->SetGridx();
1767 canvas->SetGridy();
1768 //canvas->SetRightMargin(0.05);
1769 //canvas->SetTopMargin(0.05);
1770
1771 canvas->Divide(2, 1);
1772
1773 canvas->cd(2);
1774 gPad->SetLogy();
1775
1776 Int_t selectedMult = 30;
1777 Int_t yMax = 200000;
1778
1779 TH1* full = proj->ProjectionX("full");
1780 TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
1781
1782 full->SetStats(kFALSE);
1783 full->GetXaxis()->SetRangeUser(0, 200);
1784 full->GetYaxis()->SetRangeUser(5, yMax);
1785 full->SetTitle(";multiplicity");
1786
1787 selected->SetLineColor(0);
1788 selected->SetMarkerColor(2);
1789 selected->SetMarkerStyle(7);
1790
1791 full->Draw();
1792 selected->Draw("SAME P");
1793
1794 TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
1795 legend->SetFillColor(0);
1796 legend->AddEntry(full, "true");
1797 legend->AddEntry(selected, "measured");
1798 legend->Draw();
1799
1800 TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
1801 line->SetLineWidth(2);
1802 line->Draw();
1803
1804 canvas->cd(1);
1805 gPad->SetLogy();
1806
6d81c2de 1807 full = proj->ProjectionY("full2");
1808 selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
0b4bfd98 1809
1810 full->SetStats(kFALSE);
1811 full->GetXaxis()->SetRangeUser(0, 200);
1812 full->GetYaxis()->SetRangeUser(5, yMax);
1813 full->SetTitle(";multiplicity");
1814
1815 full->SetLineColor(0);
1816 full->SetMarkerColor(2);
1817 full->SetMarkerStyle(7);
1818
1819 full->Draw("P");
1820 selected->Draw("SAME");
1821
1822 legend->Draw();
1823
6d81c2de 1824 line = new TLine(selectedMult, 5, selectedMult, yMax);
0b4bfd98 1825 line->SetLineWidth(2);
1826 line->Draw();
1827
1828 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1829}
1830
1831void SystematicpTSpectrum()
1832{
1833 gSystem->Load("libPWG0base");
1834
1835 TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
1836 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1837 mult->LoadHistograms("Multiplicity");
1838
1839 TFile::Open("multiplicityMC_100k_syst.root");
1840 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1841 mult2->LoadHistograms("Multiplicity");
1842
1843 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1844 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1845 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1846
1847 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1848 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1849
1850 DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
1851}
1852
1853// to be deleted
1854/*void covMatrix(Bool_t mc = kTRUE)
1855{
1856 gSystem->Load("libPWG0base");
1857
1858 TFile::Open(correctionFile);
1859 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1860 mult->LoadHistograms("Multiplicity");
1861
1862 TFile::Open(measuredFile);
1863 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1864 mult2->LoadHistograms("Multiplicity");
1865
1866 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1867
1868 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1869
1870 mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0));
1871}*/
1872
1873Double_t FitPtFunc(Double_t *x, Double_t *par)
1874{
1875 Double_t xx = x[0];
1876
1877 Float_t val1 = par[1] + par[2] * xx + par[3] * xx * xx;
1878 Float_t val2 = TMath::Exp(par[4] + par[5] * xx);
1879
1880 const Float_t kTransitionWidth = 0;
1881
1882 // power law part
1883 if (xx < par[0] - kTransitionWidth)
1884 {
1885 return val1;
1886 }
1887 /*else if (xx < par[0] + kTransitionWidth)
1888 {
1889 // smooth transition
1890 Float_t factor = (xx - par[0] + kTransitionWidth) / kTransitionWidth / 2;
1891 return (1 - factor) * val1 + factor * val2;
1892 }*/
1893 else
1894 {
1895 return val2;
1896 }
1897}
1898
dd367a14 1899void FitPtNew(const char* fileName = "TruePt14TeV.root")
1900{
1901 gSystem->Load("libANALYSIS");
1902 gSystem->Load("libPWG0base");
1903
1904 TFile::Open(fileName);
1905
1906 TH1* genePt = (TH1*) gFile->Get("fHistPt");
1907 genePt->Sumw2();
1908
1909 // normalize by bin width
1910 for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1911 genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1912
1913 genePt->GetXaxis()->SetRangeUser(0.05, 2.0);
1914
1915 genePt->Scale(1.0 / genePt->Integral());
1916
1917 TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x*x)", 0.001, 100);
1918 //func->SetLineColor(2);
1919 func->SetParameters(1, -1);
1920
1921 genePt->SetMarkerStyle(25);
1922 genePt->SetTitle("");
1923 genePt->SetStats(kFALSE);
1924 genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
1925 //func->Draw("SAME");
1926
1927 genePt->Fit(func, "0", "", 0.05, 1);
1928
1929 new TCanvas;
1930 genePt->DrawCopy("P");
1931 func->SetRange(0.02, 8);
1932 func->DrawCopy("SAME");
1933 gPad->SetLogy();
1934}
1935
0b4bfd98 1936void FitPt(const char* fileName = "firstplots100k_truept.root")
1937{
1938 gSystem->Load("libPWG0base");
1939
1940 TFile::Open(fileName);
1941
1942 /*
1943 // merge corrections
1944 AliCorrection* correction[4];
1945 TList list;
1946
1947 for (Int_t i=0; i<4; ++i)
1948 {
1949 Printf("correction %d", i);
1950
1951 TString name; name.Form("correction_%d", i);
1952 correction[i] = new AliCorrection(name, name);
1953 correction[i]->LoadHistograms();
1954
1955 if (i > 0)
1956 list.Add(correction[i]);
1957 }
1958
1959 correction[0]->Merge(&list);
1960
1961 TH3* gene = correction[0]->GetTrackCorrection()->GetGeneratedHistogram();
1962
1963 // limit vtx, eta axis
1964 gene->GetXaxis()->SetRangeUser(-5.9, 5.9);
1965 gene->GetYaxis()->SetRangeUser(-1.99, 0.99);
1966
1967 TH1* genePt = gene->Project3D("z");*/
1968 TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
dd367a14 1969 if (!genePt)
1970 genePt = (TH1*) gFile->Get("fHistPt");
1971
0b4bfd98 1972 genePt->Sumw2();
1973
1974 //genePt->Scale(1.0 / genePt->Integral());
1975
1976 // normalize by bin width
1977 for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1978 genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1979
1980 /// genePt->GetXaxis()->GetBinCenter(x));
1981
1982 genePt->GetXaxis()->SetRangeUser(0, 7.9);
6d81c2de 1983 //genePt->GetYaxis()->SetTitle("a.u.");
0b4bfd98 1984
dd367a14 1985 //TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x*x)", 0.001, 100);
0b4bfd98 1986 TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
1987 //func->SetLineColor(2);
1988 func->SetParameters(1, -1, 1, 1, 1);
1989 func->SetParLimits(3, 1, 10);
1990 func->SetParLimits(4, 0, 10);
1991
1992 //TF1* func = new TF1("func", "[1]*x**[0]", 0.001, 100);
1993
1994 //TF1* func = new TF1("func", &FitPtFunc, 0, 2, 6);
1995 //func->SetParameters(0.3, -2.34909e-01, 1.54394e+01, -3.04134e+01, 1.41912e+00, -2.79284e+00);
1996 //func->FixParameter(0, 0.314);
1997 //func->SetParLimits(0, 0.1, 0.3);
1998
1999 genePt->SetMarkerStyle(25);
2000 genePt->SetTitle("");
2001 genePt->SetStats(kFALSE);
2002 genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
2003 //func->Draw("SAME");
2004
2005 // fit only exp. part
2006 func->SetParameters(1, -1);
2007 func->FixParameter(2, 0);
2008 func->FixParameter(3, 1);
2009 func->FixParameter(4, 1);
2010 genePt->Fit(func, "0", "", 0.2, 1);
2011
2012 new TCanvas;
2013 genePt->DrawCopy("P");
2014 func->SetRange(0.02, 8);
2015 func->DrawCopy("SAME");
2016 gPad->SetLogy();
2017
2018 // now fix exp. parameters and fit second part
2019 Double_t param0 = func->GetParameter(0);
2020 Double_t param1 = func->GetParameter(1);
2021 func->SetParameters(0, -1, 1, 1, 1);
2022 func->FixParameter(0, 0);
2023 func->FixParameter(1, -1);
2024 func->ReleaseParameter(2);
2025 func->ReleaseParameter(3);
2026 func->ReleaseParameter(4);
2027 func->SetParLimits(3, 1, 10);
2028 func->SetParLimits(4, 0, 10);
2029
2030 genePt->Fit(func, "0", "", 1.5, 4);
2031
2032 new TCanvas;
2033 genePt->DrawCopy("P");
2034 func->SetRange(0.02, 8);
2035 func->DrawCopy("SAME");
2036 gPad->SetLogy();
2037
2038 // fit both
2039 func->SetParameter(0, param0);
2040 func->SetParameter(1, param1);
2041 func->ReleaseParameter(0);
2042 func->ReleaseParameter(1);
2043
2044 new TCanvas;
2045 genePt->DrawCopy("P");
6d81c2de 2046 func->SetRange(0.02, 5);
0b4bfd98 2047 func->DrawCopy("SAME");
2048 gPad->SetLogy();
2049
2050 genePt->Fit(func, "0", "", 0.2, 4);
2051
6d81c2de 2052 TCanvas* canvas = new TCanvas("FitPt", "FitPt", 800, 400);
2053 canvas->Divide(2, 1);
2054 canvas->cd(1);
2055
2056 gPad->SetGridx();
2057 gPad->SetGridy();
2058 gPad->SetLeftMargin(0.13);
2059 gPad->SetRightMargin(0.05);
2060 gPad->SetTopMargin(0.05);
0b4bfd98 2061
6d81c2de 2062 genePt->GetXaxis()->SetRangeUser(0, 4.9);
2063 genePt->GetYaxis()->SetRangeUser(1e-2, 1e4);
2064 genePt->GetYaxis()->SetTitleOffset(1.4);
2065 genePt->GetXaxis()->SetTitleOffset(1.1);
0b4bfd98 2066 genePt->DrawCopy("P");
6d81c2de 2067 func->SetRange(0.02, 5);
2068 func->DrawCopy("SAME");
2069 gPad->SetLogy();
2070
2071 canvas->cd(2);
2072
2073 TH1* genePtClone = (TH1*) genePt->Clone("genePtClone");
2074 genePtClone->Reset();
2075 genePtClone->DrawCopy("P");
2076
2077 gPad->SetGridx();
2078 gPad->SetGridy();
2079 gPad->SetLeftMargin(0.13);
2080 gPad->SetRightMargin(0.05);
2081 gPad->SetTopMargin(0.05);
2082
0b4bfd98 2083 func->DrawCopy("SAME");
2084 gPad->SetLogy();
2085
6d81c2de 2086 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2087
0b4bfd98 2088 TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2089
2090 TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2091
2092 TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2093
2094 for (Int_t param=0; param<5; param++)
2095 {
2096 for (Int_t sign=0; sign<2; sign++)
2097 {
2098 TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign)); // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2099 func2->SetParameters(func->GetParameters());
2100 //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2101
2102 Float_t factor = ((sign == 0) ? 0.9 : 1.1);
2103 func2->SetParameter(param, func2->GetParameter(param) * factor);
2104 //func2->Print();
2105
6d81c2de 2106 canvas->cd(2);
0b4bfd98 2107 func2->SetLineWidth(1);
6d81c2de 2108 func2->SetLineColor(2);
0b4bfd98 2109 func2->DrawCopy("SAME");
2110
2111 canvas2->cd();
2112 TH1* second = func2->GetHistogram();
2113 second->Divide(first);
2114 second->SetLineColor(param + 1);
2115 second->GetYaxis()->SetRangeUser(0, 2);
2116 second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2117 second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2118 }
2119 }
2120
2121 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
6d81c2de 2122 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
0b4bfd98 2123
2124 file->Close();
2125}
2126
6d81c2de 2127void DrawSystematicpT()
2128{
2129 TFile* file = TFile::Open("SystematicpT.root");
2130
2131 TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2132 TH1* result2 = (TH1*) file->Get("result_unity");
2133
2134 TH1* mcHist[12];
2135 TH1* result[12];
2136
2137 Int_t nParams = 5;
2138
2139 for (Int_t id=0; id<nParams*2; ++id)
2140 {
2141 mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2142 result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2143 }
2144
2145 DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2146
2147 //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2148
2149 DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2150
2151 //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2152
2153 // does not make sense: mc is different
2154 //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2155}
2156
2157void SystematicpT(Bool_t chi2 = 1)
0b4bfd98 2158{
2159 gSystem->Load("libPWG0base");
2160
6d81c2de 2161 TFile::Open("ptspectrum900.root");
0b4bfd98 2162 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2163 mult->LoadHistograms("Multiplicity");
2164
2165 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2166
2167 TH1* mcHist[12];
2168 TH1* result[12];
2169
6d81c2de 2170 Int_t nParams = 5;
0b4bfd98 2171
2172 for (Int_t param=0; param<nParams; param++)
2173 {
2174 for (Int_t sign=0; sign<2; sign++)
2175 {
2176 // calculate result with systematic effect
6d81c2de 2177 TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
0b4bfd98 2178 mult2->LoadHistograms("Multiplicity");
2179
2180 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2181
6d81c2de 2182 if (chi2)
2183 {
2184 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2185 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2186 }
2187 else
2188 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
0b4bfd98 2189
2190 Int_t id = param * 2 + sign;
2191
2192 mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2193 result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2194
2195 TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2196 DrawResultRatio(mcHist[id], result[id], tmp);
2197 }
2198 }
2199
2200 // calculate normal result
2201 TFile::Open("ptspectrum100_1.root");
2202 mult2->LoadHistograms("Multiplicity");
6d81c2de 2203 TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
0b4bfd98 2204
2205 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2206
6d81c2de 2207 if (chi2)
2208 {
2209 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2210 }
2211 else
2212 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
0b4bfd98 2213
6d81c2de 2214 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
0b4bfd98 2215
6d81c2de 2216 TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2217 mcHist2->Write();
2218 result2->Write();
2219 for (Int_t id=0; id<nParams*2; ++id)
2220 {
2221 mcHist[id]->Write();
2222 result[id]->Write();
2223 }
2224 file->Close();
2225
2226 DrawSystematicpT();
2227}
0b4bfd98 2228
44466df2 2229void DrawSystematicpT2()
2230{
2231 //displayRange = 200;
2232
2233 // read from file
2234 TFile* file = TFile::Open("SystematicpT2.root");
2235 TH1* mcHist = (TH1*) file->Get("mymc");
2236 TH1* result[12];
2237 result[0] = (TH1*) file->Get("result_unity");
2238 Int_t nParams = 5;
2239 for (Int_t id=0; id<nParams*2; ++id)
2240 result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
2241
2242 DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
2243 DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
2244 DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
2245}
2246
2247void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
6d81c2de 2248{
2249 gSystem->Load("libPWG0base");
44466df2 2250
2251 if (tpc)
2252 {
2253 SetTPC();
2254 TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2255 }
2256 else
2257 TFile::Open("ptspectrum100_1.root");
2258
2259 AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2260 measured->LoadHistograms("Multiplicity");
2261 TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2262
2263 TH1* result[12];
2264
2265 Int_t nParams = 5;
2266
2267 // -1 = unity change, 0...4 parameters
2268 for (Int_t id=-1; id<nParams*2; id++)
2269 {
2270 Int_t param = id / 2;
2271 Int_t sign = id % 2;
2272
2273 TString idStr;
2274 if (id == -1)
2275 {
2276 idStr = "unity";
2277 }
2278 else
2279 idStr.Form("%d_%d", param, sign);
2280
2281 // calculate result with systematic effect
2282 if (tpc)
2283 {
2284 TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
2285 }
2286 else
2287 TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
2288
2289 AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2290 response->LoadHistograms("Multiplicity");
2291
2292 response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
2293
2294 if (chi2)
2295 {
2296 response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2297 response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2298 }
2299 else
2300 response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2301
2302 result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
2303
2304 TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
2305 DrawResultRatio(mcHist, result[id+1], tmp);
2306 }
2307
2308 TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
2309 mcHist->Write();
2310 for (Int_t id=0; id<nParams*2+1; ++id)
2311 result[id]->Write();
2312 file->Close();
2313
2314 DrawSystematicpT2();
2315}
2316
2317void SystematicpTCutOff(Bool_t chi2 = kTRUE)
2318{
2319 // only needed for TPC
6d81c2de 2320 SetTPC();
0b4bfd98 2321
44466df2 2322 gSystem->Load("libPWG0base");
2323
6d81c2de 2324 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2325 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2326 mult->LoadHistograms("Multiplicity");
0b4bfd98 2327
44466df2 2328 TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
6d81c2de 2329 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2330 mult2->LoadHistograms("Multiplicity");
0b4bfd98 2331
6d81c2de 2332 // "normal" result
2333 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
44466df2 2334
2335 if (chi2)
2336 {
2337 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2338 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2339 }
2340 else
2341 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
0b4bfd98 2342
6d81c2de 2343 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2344 TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
0b4bfd98 2345
44466df2 2346 TH1* syst[2];
2347
6d81c2de 2348 // change of pt spectrum (down)
2349 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2350 AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2351 mult3->LoadHistograms("Multiplicity");
6d81c2de 2352 mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
44466df2 2353 if (chi2)
2354 {
2355 mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2356 mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2357 }
2358 else
2359 mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2360 syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2361
2362 // change of pt spectrum (up)
2363 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
2364 AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
2365 mult4->LoadHistograms("Multiplicity");
2366 mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2367 if (chi2)
2368 {
2369 mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2370 mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2371 }
2372 else
2373 mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2374 syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
6d81c2de 2375
44466df2 2376 DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
6d81c2de 2377
44466df2 2378 Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
2379 Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
2380}
2381
2382TH1* SystematicsSummary(Bool_t tpc = 1)
2383{
2384 Int_t nEffects = 7;
2385
2386 TH1* effects[10];
2387 const char** names = 0;
2388 Int_t colors[] = { 1, 2, 3, 4, 6, 7, 8 };
2389 Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
2390
2391 for (Int_t i=0; i<nEffects; ++i)
2392 effects[i] = new TH1F("SystematicsSummary", ";true multiplicity;Effect", 201, -0.5, 200.5);
2393
2394 if (tpc)
2395 {
2396 SetTPC();
2397
2398 const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
2399 names = namesTPC;
2400
2401 // method
2402 TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
2403 TH1* hist = (TH1*) file->Get("errorBoth");
2404
2405 // smooth a bit, but skip 0 bin
2406 effects[0]->SetBinContent(2, hist->GetBinContent(2));
2407 for (Int_t i=3; i<=200; ++i)
2408 effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2409
2410 // relative x-section
2411 effects[1]->SetBinContent(2, 0.005);
2412 effects[1]->SetBinContent(3, 0.0025);
2413 effects[1]->SetBinContent(4, 0.0025);
2414
2415 // particle composition
2416 for (Int_t i=2; i<=101; ++i)
2417 {
2418 if (i < 41)
2419 {
2420 effects[2]->SetBinContent(i, 0.01);
2421 }
2422 else if (i < 76)
2423 {
2424 effects[2]->SetBinContent(i, 0.02);
2425 }
2426 else
2427 effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
2428 }
2429
2430 // pt cut off (only tpc)
2431 for (Int_t i=2; i<=101; ++i)
2432 {
2433 if (i < 11)
2434 {
2435 effects[3]->SetBinContent(i, 0.05);
2436 }
2437 else if (i < 51)
2438 {
2439 effects[3]->SetBinContent(i, 0.01);
2440 }
2441 else
2442 effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
2443 }
2444
2445 // track selection (only tpc)
2446 for (Int_t i=2; i<=101; ++i)
2447 effects[4]->SetBinContent(i, 0.03);
2448
2449 // secondaries
2450 for (Int_t i=2; i<=101; ++i)
2451 effects[5]->SetBinContent(i, 0.01);
2452
2453 // pt spectrum
2454 for (Int_t i=2; i<=101; ++i)
2455 {
2456 if (i < 21)
2457 {
2458 effects[6]->SetBinContent(i, 0.05);
2459 }
2460 else if (i < 51)
2461 {
2462 effects[6]->SetBinContent(i, 0.02);
2463 }
2464 else
2465 effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
2466 }
2467
2468 }
2469 else
2470 {
2471 displayRange = 200;
2472 nEffects = 5;
2473
2474 const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Rel. cross-section", "Particle composition", "Secondaries", "p_{t} spectrum"};
2475 names = namesSPD;
2476
2477 // method
2478 TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
2479 TH1* hist = (TH1*) file->Get("errorBoth");
2480
2481 // smooth a bit, but skip 0 bin
2482 effects[0]->SetBinContent(2, hist->GetBinContent(2));
2483 for (Int_t i=3; i<=201; ++i)
2484 effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2485
2486 // relative x-section
2487 effects[1]->SetBinContent(2, 0.01);
2488 effects[1]->SetBinContent(3, 0.005);
2489
2490 // particle composition
2491 for (Int_t i=2; i<=201; ++i)
2492 {
2493 if (i < 6)
2494 {
2495 effects[2]->SetBinContent(i, 0.3);
2496 }
2497 else if (i < 11)
2498 {
2499 effects[2]->SetBinContent(i, 0.05);
2500 }
2501 else if (i < 121)
2502 {
2503 effects[2]->SetBinContent(i, 0.02);
2504 }
2505 else if (i < 151)
2506 {
2507 effects[2]->SetBinContent(i, 0.02 + 0.04 / 30 * (i - 121));
2508 }
2509 else
2510 effects[2]->SetBinContent(i, 0.06 + 0.1 / 30 * (i - 151));
2511 }
2512
2513 // secondaries
2514 for (Int_t i=2; i<=201; ++i)
2515 effects[3]->SetBinContent(i, 0.01);
2516
2517 // pt spectrum
2518 for (Int_t i=2; i<=201; ++i)
2519 {
2520 if (i < 6)
2521 {
2522 effects[4]->SetBinContent(i, 1);
2523 }
2524 else if (i < 121)
2525 {
2526 effects[4]->SetBinContent(i, 0.03);
2527 }
2528 else if (i < 151)
2529 {
2530 effects[4]->SetBinContent(i, 0.03 + 0.07 / 30 * (i - 121));
2531 }
2532 else
2533 effects[4]->SetBinContent(i, 0.1);
2534 }
2535 }
2536
2537 TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 400);
2538 canvas->SetRightMargin(0.25);
2539 canvas->SetTopMargin(0.05);
2540 TLegend* legend = new TLegend(0.2, 0.4, 0.5, 0.4 + 0.5 * nEffects / 7);
2541 legend->SetFillColor(0);
2542
2543 for (Int_t i=0; i<nEffects; ++i)
2544 {
2545 TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
2546 /*current->Reset();
2547 for (Int_t j=0; j<nEffects-i; ++j)
2548 current->Add(effects[j]);*/
2549
2550 current->SetLineColor(colors[i]);
2551 //current->SetFillColor(colors[i]);
2552 current->SetMarkerColor(colors[i]);
2553 //current->SetMarkerStyle(markers[i]);
2554
2555 current->SetStats(kFALSE);
2556 current->GetYaxis()->SetRangeUser(0, 0.4);
2557 current->GetXaxis()->SetRangeUser(0, displayRange);
2558 current->DrawCopy(((i == 0) ? "" : "SAME"));
2559 legend->AddEntry(current, names[i]);
2560
2561 TLatex* text = new TLatex(displayRange+5, current->GetBinContent(displayRange+1), names[i]);
2562 text->SetTextColor(colors[i]);
2563 text->Draw();
2564 }
2565
2566 // add total in square
2567 TH1* total = (TH1*) effects[0]->Clone("total");
2568 total->Reset();
2569
2570 for (Int_t i=0; i<nEffects; ++i)
2571 {
2572 //Printf("%d %f", i, effects[i]->GetBinContent(20));
2573 effects[i]->Multiply(effects[i]);
2574 total->Add(effects[i]);
2575 }
2576
2577 for (Int_t i=1; i<=total->GetNbinsX(); ++i)
2578 if (total->GetBinContent(i) > 0)
2579 total->SetBinContent(i, TMath::Min(sqrt(total->GetBinContent(i)), 1.0));
2580
2581 //Printf("%f", total->GetBinContent(20));
2582
2583 total->SetMarkerStyle(3);
2584 total->SetMarkerColor(1);
2585 legend->AddEntry(total, "total");
2586 total->DrawCopy("SAME P");
2587
2588 legend->Draw();
2589
2590 canvas->SaveAs(canvas->GetName());
2591
2592 return total;
2593}
2594
0f67a57c 2595void finalPlot(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE, Bool_t small = kFALSE)
44466df2 2596{
2597 gSystem->Load("libPWG0base");
2598
2599 if (tpc)
2600 SetTPC();
2601
2602 if (!chi2)
2603 Printf("WARNING: Bayesian set. This is only for test!");
2604
2605 // systematic error
2606 TH1* error = SystematicsSummary(tpc);
2607
2608 TFile::Open(correctionFile);
2609 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2610 mult->LoadHistograms("Multiplicity");
2611
2612 TFile::Open(measuredFile);
2613 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2614 mult2->LoadHistograms("Multiplicity");
2615
2616 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2617
2618 if (chi2)
2619 {
2620 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2621 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kINEL);
2622 }
2623 else
2624 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kINEL, 1, 100, 0, kFALSE);
2625
2626 TH1* mcHist = mult2->GetMultiplicityINEL(etaRange)->ProjectionY("mymc");
2627 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2628
2629 DrawResultRatio(mcHist, result, "finalPlotCheck.eps");
2630
2631 // normalize result
2632 result->Scale(1.0 / result->Integral(2, 200));
2633
2634 result->GetXaxis()->SetRangeUser(0, ((tpc) ? displayRange : 200));
2635 result->SetBinContent(1, 0); result->SetBinError(1, 0);
2636 result->SetTitle(";true multiplicity;Probability");
2637 result->SetLineColor(1);
2638 result->SetStats(kFALSE);
2639
2640 TH1* systError = (TH1*) result->Clone("systError");
2641 for (Int_t i=2; i<=systError->GetNbinsX(); ++i)
2642 systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
2643
2644 // change error drawing style
2645 systError->SetFillColor(15);
2646
0f67a57c 2647 TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", (small) ? 600 : 800, 400);
44466df2 2648 canvas->SetRightMargin(0.05);
2649 canvas->SetTopMargin(0.05);
2650
2651 systError->Draw("E2 ][");
2652 result->DrawCopy("SAME E ][");
2653 canvas->SetLogy();
2654
2655 //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
2656 TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
2657 text->SetFillColor(0);
2658 text->SetTextAlign(12);
2659 text->AddText("Systematic errors summed quadratically");
2660 text->AddText("0.6 million minimum bias events");
2661 text->AddText("corrected to inelastic events");
2662 text->Draw("B");
2663
2664 TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
2665 text2->SetFillColor(0);
2666 text2->SetTextAlign(12);
2667 text2->AddText("#sqrt{s} = 14 TeV");
2668 if (tpc)
2669 {
2670 text2->AddText("|#eta| < 0.9");
2671 }
2672 else
2673 text2->AddText("|#eta| < 2.0");
2674 text2->AddText("simulated data (PYTHIA)");
2675 text2->Draw("B");
2676
2677 if (tpc)
2678 {
2679 TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
2680 text3->SetNDC();
2681 text3->Draw();
2682 }
2683 else
2684 {
2685 TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
2686 text3->SetNDC();
2687 text3->Draw();
2688 }
2689
2690 // alice logo
2691 TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
2692 pad->Draw();
2693 pad->cd();
2694 TImage* img = TImage::Open("$HOME/alice.png");
2695 img->SetImageQuality(TAttImage::kImgBest);
2696 img->Draw();
2697
2698 canvas->Modified();
2699
2700/* TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
2701 text->SetTextSize(0.04);
2702 text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
2703 text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
2704 text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
2705 text->Draw();*/
2706
2707
2708 canvas->SaveAs(canvas->GetName());
0b4bfd98 2709}
0f67a57c 2710
2711void BlobelUnfoldingExample()
2712{
2713 const Int_t kSize = 20;
2714
2715 TMatrixD matrix(kSize, kSize);
2716 for (Int_t x=0; x<kSize; x++)
2717 {
2718 for (Int_t y=0; y<kSize; y++)
2719 {
2720 if (x == y)
2721 {
2722 if (x == 0 || x == kSize -1)
2723 {
2724 matrix(x, y) = 0.75;
2725 }
2726 else
2727 matrix(x, y) = 0.5;
2728 }
2729 else if (TMath::Abs(x - y) == 1)
2730 {
2731 matrix(x, y) = 0.25;
2732 }
2733 }
2734 }
2735
2736 //matrix.Print();
2737
2738 TMatrixD inverted(matrix);
2739 inverted.Invert();
2740
2741 //inverted.Print();
2742
2743 TH1F* inputDist = new TH1F("inputDist", ";t;#tilde{T}(t)", kSize, -0.5, (Float_t) kSize - 0.5);
2744 TVectorD inputDistVector(kSize);
2745 TH1F* unfolded = inputDist->Clone("unfolded");
2746 TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
2747 measuredIdealDist->SetTitle(";m;#tilde{M}(m)");
2748 TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
2749
2750 TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
2751 // norm: 1/(sqrt(2pi)sigma)
2752 gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
2753 //gaus->Print();
2754
2755 for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
2756 {
2757 Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
2758 inputDist->SetBinContent(x, value);
2759 inputDistVector(x-1) = value;
2760 }
2761
2762 TVectorD measuredDistIdealVector = matrix * inputDistVector;
2763
2764 for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
2765 measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
2766
2767 measuredDist->FillRandom(measuredIdealDist, 10000);
2768
2769 // fill error matrix before scaling
2770 TMatrixD covarianceMatrixMeasured(kSize, kSize);
2771 for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
2772 covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
2773
2774 TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
2775 //covarianceMatrix.Print();
2776
2777 TVectorD measuredDistVector(kSize);
2778 for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
2779 measuredDistVector(x-1) = measuredDist->GetBinContent(x);
2780
2781 TVectorD unfoldedVector = inverted * measuredDistVector;
2782 for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
2783 unfolded->SetBinContent(x, unfoldedVector(x-1));
2784
2785 TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1000, 500);
2786 canvas->SetTopMargin(0.05);
2787 canvas->Divide(2, 1);
2788
2789 canvas->cd(1);
2790 canvas->cd(1)->SetLeftMargin(0.15);
2791 canvas->cd(1)->SetRightMargin(0.05);
2792 measuredDist->GetYaxis()->SetTitleOffset(1.7);
2793 measuredDist->SetStats(0);
2794 measuredDist->DrawCopy();
2795 gaus->Draw("SAME");
2796
2797 canvas->cd(2);
2798 canvas->cd(2)->SetLeftMargin(0.15);
2799 canvas->cd(2)->SetRightMargin(0.05);
2800 unfolded->GetYaxis()->SetTitleOffset(1.7);
2801 unfolded->SetStats(0);
2802 unfolded->DrawCopy();
2803 gaus->Draw("SAME");
2804
2805 canvas->SaveAs("BlobelUnfoldingExample.eps");
2806}
2807
2808void E735Fit()
2809{
2810 TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
2811 fCurrentESD->Sumw2();
2812
2813 // Open the input stream
2814 ifstream in;
2815 in.open("e735data.txt");
2816
2817 while(in.good())
2818 {
2819 Float_t x, y, ye;
2820 in >> x >> y >> ye;
2821
2822 //Printf("%f %f %f", x, y, ye);
2823 fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
2824 fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
2825 }
2826
2827 in.close();
2828
2829 //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
2830
2831 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2832
2833 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])");
2834 func->SetParNames("scaling", "averagen", "k");
2835 func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
2836 func->SetParLimits(1, 0.001, 1000);
2837 func->SetParLimits(2, 0.001, 1000);
2838 func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
2839
2840 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
2841 lognormal->SetParNames("scaling", "mean", "sigma");
2842 lognormal->SetParameters(1, 1, 1);
2843 lognormal->SetParLimits(0, 0, 10);
2844 lognormal->SetParLimits(1, 0, 100);
2845 lognormal->SetParLimits(2, 1e-3, 10);
2846
2847 TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
2848 fCurrentESD->SetStats(kFALSE);
2849 fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
2850 fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
2851 fCurrentESD->Draw("");
2852 fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
2853 fCurrentESD->Fit(func, "0", "", 0, 150);
2854 func->SetRange(0, 250);
2855 func->Draw("SAME");
2856 printf("chi2 = %f\n", func->GetChisquare());
2857
2858 fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
2859 lognormal->SetLineColor(2);
2860 lognormal->SetLineStyle(2);
2861 lognormal->SetRange(0, 250);
2862 lognormal->Draw("SAME");
2863
2864 gPad->SetLogy();
2865
2866 canvas->SaveAs("E735Fit.eps");
2867}