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