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