]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG0/multiplicity/plots.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
correction for slices in AliCFContainer ; simplification of several
[u/mrichter/AliRoot.git] / PWG0 / multiplicity / plots.C
CommitLineData
f3eb27f6 1/* $Id$ */
2
3//
4// plots for the note about multplicity measurements
5//
6
7#if !defined(__CINT__) || defined(__MAKECINT__)
8
9#include <TCanvas.h>
10#include <TPad.h>
11#include <TH1F.h>
12#include <TH2F.h>
13#include <TH3F.h>
14#include <TLine.h>
15#include <TF1.h>
16#include <TSystem.h>
17#include <TFile.h>
18#include <TLegend.h>
19#include <TStopwatch.h>
20#include <TROOT.h>
21#include <TGraph.h>
22#include <TMath.h>
23#include <TPaveText.h>
24#include <TImage.h>
25#include <TLatex.h>
26
27#include "AliMultiplicityCorrection.h"
28#include "AliCorrection.h"
29#include "AliCorrectionMatrix3D.h"
30
31#endif
32
69b09e3b 33const char* correctionFile = "multiplicityMC.root";
34const char* measuredFile = "multiplicityESD.root";
35Int_t etaRange = 1;
36Int_t displayRange = 80; // axis range
f3eb27f6 37Int_t ratioRange = 151; // range to calculate difference
69b09e3b 38Int_t longDisplayRange = 120;
f3eb27f6 39
40const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
41const char* measuredFileTPC = "multiplicityMC_TPC_0.6M.root";
42Int_t etaRangeTPC = 1;
43
5a6310fe 44void loadlibs()
45{
46 gSystem->Load("libANALYSIS");
47 gSystem->Load("libPWG0base");
48}
49
f3eb27f6 50void SetTPC()
51{
52 correctionFile = correctionFileTPC;
53 measuredFile = measuredFileTPC;
54 etaRange = etaRangeTPC;
55 displayRange = 100;
56 ratioRange = 76;
57 longDisplayRange = 100;
58}
59
69b09e3b 60const char* GetMultLabel(Int_t etaR = -1, Bool_t trueM = kTRUE)
61{
62 if (etaR == -1)
63 etaR = etaRange;
64
65 TString tmpStr((trueM) ? "True " : "Measured ");
66
67 tmpStr += "multiplicity";
68 //return Form("%s", tmpStr.Data());
69
70 if (etaR == 4)
71 {
72 tmpStr += " (full phase space)";
73 }
74 else
75 tmpStr += Form(" in |#eta| < %.1f", (etaR+1)* 0.5);
76 return Form("%s", tmpStr.Data());
77}
78
f3eb27f6 79void Smooth(TH1* hist, Int_t windowWidth = 20)
80{
81 TH1* clone = (TH1*) hist->Clone("clone");
82 for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin)
83 {
84 Int_t min = TMath::Max(2, bin-windowWidth);
85 Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
86 Float_t average = clone->Integral(min, max) / (max - min + 1);
87
88 hist->SetBinContent(bin, average);
89 hist->SetBinError(bin, 0);
90 }
91
92 delete clone;
93}
94
69b09e3b 95void responseMatrixPlot(const char* fileName = 0)
f3eb27f6 96{
69b09e3b 97 loadlibs();
f3eb27f6 98
99 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
100
69b09e3b 101 if (fileName == 0)
102 fileName = correctionFile;
103
104 TFile::Open(fileName);
f3eb27f6 105 mult->LoadHistograms("Multiplicity");
106
69b09e3b 107 // empty under/overflow bins in x, otherwise Project3D takes them into account
108 TH1* hist = mult->GetCorrelation(etaRange);
109 for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
110 {
111 for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
112 {
113 hist->SetBinContent(0, y, z, 0);
114 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
115 }
116 }
117 hist = ((TH3*) mult->GetCorrelation(etaRange))->Project3D("zy");
f3eb27f6 118 hist->SetStats(kFALSE);
119
69b09e3b 120 hist->SetTitle(Form(";%s;%s;Entries", GetMultLabel(), GetMultLabel(etaRange, kFALSE)));
f3eb27f6 121 hist->GetXaxis()->SetRangeUser(0, longDisplayRange);
122 hist->GetYaxis()->SetRangeUser(0, longDisplayRange);
69b09e3b 123
124 hist->GetYaxis()->SetTitleOffset(1.3);
125 hist->GetZaxis()->SetTitleOffset(1.2);
f3eb27f6 126
69b09e3b 127 TCanvas* canvas = new TCanvas("c1", "c1", 600, 600);
f3eb27f6 128 canvas->SetRightMargin(0.15);
129 canvas->SetTopMargin(0.05);
130
131 gPad->SetLogz();
132 hist->Draw("COLZ");
133
134 canvas->SaveAs("responsematrix.eps");
135}
136
69b09e3b 137void multPythiaPhojet()
138{
139 loadlibs();
140
141 TFile::Open("LHC08c11_10TeV_0.5T/mb1/spd/multiplicity.root");
142 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
143 mult->LoadHistograms("Multiplicity");
144 hist1 = mult->GetMultiplicityINEL(1)->ProjectionY();
145 hist1->Sumw2();
146
147 TFile::Open("LHC08c15_10TeV_0.5T_Phojet/mb1/spd/multiplicity.root");
148 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
149 mult->LoadHistograms("Multiplicity");
150 hist2 = mult->GetMultiplicityINEL(1)->ProjectionY();
151 hist2->Sumw2();
152
153 legend = new TLegend(0.6, 0.7, 0.9, 0.9);
154 legend->SetFillColor(0);
155 legend->SetTextSize(0.04);
156 legend->AddEntry(hist1, "Pythia", "L");
157 legend->AddEntry(hist2, "Phojet", "L");
158
159 c1 = new TCanvas("c", "c", 600, 600);
160 c1->SetTopMargin(0.05);
161 c1->SetRightMargin(0.05);
162 c1->SetLeftMargin(0.12);
163 c1->SetGridx();
164 c1->SetGridy();
165 c1->SetLogy();
166
167 //hist1->SetMarkerStyle(20);
168 //hist2->SetMarkerStyle(24);
169 //hist2->SetMarkerColor(2);
170 hist1->SetLineWidth(2);
171 hist2->SetLineWidth(2);
172 hist2->SetLineStyle(2);
173 hist2->SetLineColor(2);
174
175 hist1->Scale(1.0 / hist1->Integral());
176 hist2->Scale(1.0 / hist2->Integral());
177
178 hist1->SetStats(0);
179 hist1->GetYaxis()->SetTitleOffset(1.3);
180 hist1->GetXaxis()->SetRangeUser(0, 100);
181 hist1->SetTitle(";N_{ch};P(N_{ch})");
182
183 hist1->Draw("");
184 hist2->Draw("SAME");
185 legend->Draw();
186
187 c1->SaveAs("mult_pythia_phojet.eps");
188}
189
f3eb27f6 190TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
191{
192 // normalize unfolded result to mc hist
69b09e3b 193 result->Scale(1.0 / result->Integral(2, displayRange));
194 result->Scale(mcHist->Integral(2, displayRange));
f3eb27f6 195
196 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
197 canvas->Range(0, 0, 1, 1);
198
199 TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
200 pad1->Draw();
201
202 TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
203 pad2->Draw();
204
205 pad1->SetRightMargin(0.05);
206 pad2->SetRightMargin(0.05);
207
208 // no border between them
209 pad1->SetBottomMargin(0);
210 pad2->SetTopMargin(0);
211
212 pad1->cd();
69b09e3b 213 pad1->SetGridx();
214 pad1->SetGridy();
f3eb27f6 215
216 mcHist->GetXaxis()->SetLabelSize(0.06);
217 mcHist->GetYaxis()->SetLabelSize(0.06);
218 mcHist->GetXaxis()->SetTitleSize(0.06);
219 mcHist->GetYaxis()->SetTitleSize(0.06);
220 mcHist->GetYaxis()->SetTitleOffset(0.6);
221
222 mcHist->GetXaxis()->SetRangeUser(0, displayRange);
223
69b09e3b 224 mcHist->SetTitle(Form(";%s;Entries", GetMultLabel()));
f3eb27f6 225 mcHist->SetStats(kFALSE);
226
227 mcHist->DrawCopy("HIST E");
228 gPad->SetLogy();
229
230 result->SetLineColor(2);
69b09e3b 231 result->SetMarkerColor(2);
232 result->SetMarkerStyle(5);
233 result->DrawCopy("SAME PE");
f3eb27f6 234
235 TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
69b09e3b 236 legend->AddEntry(mcHist, "True distribution");
237 legend->AddEntry(result, "Unfolded distribution", "P");
f3eb27f6 238 legend->SetFillColor(0);
239 legend->Draw();
240
241 pad2->cd();
242 pad2->SetBottomMargin(0.15);
243
244 // calculate ratio
245 mcHist->Sumw2();
246 TH1* ratio = (TH1*) mcHist->Clone("ratio");
247 result->Sumw2();
248 ratio->Divide(ratio, result, 1, 1, "");
249 ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
250 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
251
252 ratio->DrawCopy();
253
254 // get average of ratio
255 Float_t sum = 0;
256 for (Int_t i=2; i<=ratioRange; ++i)
257 {
258 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
259 }
260 sum /= ratioRange-1;
261
262 printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
263
264 TLine* line = new TLine(0, 1, displayRange, 1);
265 line->SetLineWidth(2);
266 line->Draw();
267
268 line = new TLine(0, 1.1, displayRange, 1.1);
269 line->SetLineWidth(2);
270 line->SetLineStyle(2);
271 line->Draw();
272 line = new TLine(0, 0.9, displayRange, 0.9);
273 line->SetLineWidth(2);
274 line->SetLineStyle(2);
275 line->Draw();
276
277 canvas->Modified();
278
279 canvas->SaveAs(epsName);
280
281 return canvas;
282}
283
284TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
285{
286 // draws the 3 plots in the upper plot
287 // draws the ratio between result1 and result2 in the lower plot
288
69b09e3b 289 // normalize unfolded result to mc hist
290 result1->Scale(1.0 / result1->Integral(2, displayRange));
291 result1->Scale(mcHist->Integral(2, displayRange));
292 result2->Scale(1.0 / result2->Integral(2, displayRange));
293 result2->Scale(mcHist->Integral(2, displayRange));
294
295 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
296 canvas->Range(0, 0, 1, 1);
297
298 TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
299 pad1->Draw();
300
301 TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
302 pad2->Draw();
303
304 pad1->SetRightMargin(0.05);
305 pad2->SetRightMargin(0.05);
306
307 // no border between them
308 pad1->SetBottomMargin(0);
309 pad2->SetTopMargin(0);
310
311 pad1->cd();
312 gPad->SetGridx();
313 gPad->SetGridy();
314
315 mcHist->GetXaxis()->SetLabelSize(0.06);
316 mcHist->GetYaxis()->SetLabelSize(0.06);
317 mcHist->GetXaxis()->SetTitleSize(0.06);
318 mcHist->GetYaxis()->SetTitleSize(0.06);
319 mcHist->GetYaxis()->SetTitleOffset(0.6);
320
321 mcHist->GetXaxis()->SetRangeUser(0, displayRange);
322
323 mcHist->SetTitle(";true multiplicity;Entries");
324 mcHist->SetStats(kFALSE);
325
326 mcHist->DrawCopy("HIST E");
327 gPad->SetLogy();
328
329 result1->SetLineColor(2);
330 result1->SetMarkerStyle(24);
331 result1->DrawCopy("SAME HISTE");
332
333 result2->SetLineColor(4);
334 result2->SetMarkerColor(4);
335 result2->DrawCopy("SAME HISTE");
336
337 TLegend* legend = new TLegend(0.5, 0.6, 0.95, 0.9);
338 legend->AddEntry(mcHist, "True distribution");
339 legend->AddEntry(result1, "Unfolded distribution (syst)", "P");
340 legend->AddEntry(result2, "Unfolded distribution (normal)", "P");
341 legend->SetFillColor(0);
342 legend->SetTextSize(0.06);
343 legend->Draw();
344
345 pad2->cd();
346 pad2->SetBottomMargin(0.15);
347 //gPad->SetGridx();
348 //gPad->SetGridy();
349
350 result1->GetXaxis()->SetLabelSize(0.06);
351 result1->GetYaxis()->SetLabelSize(0.06);
352 result1->GetXaxis()->SetTitleSize(0.06);
353 result1->GetYaxis()->SetTitleSize(0.06);
354 result1->GetYaxis()->SetTitleOffset(0.6);
355
356 result1->GetXaxis()->SetRangeUser(0, displayRange);
357
358 result1->SetTitle(Form(";%s;Entries", GetMultLabel()));
359 result1->SetStats(kFALSE);
360
361 // calculate ratio
362 result1->Sumw2();
363 TH1* ratio = (TH1*) result1->Clone("ratio");
364 result2->Sumw2();
365 ratio->Divide(ratio, result2, 1, 1, "");
366 ratio->SetLineColor(1);
367 ratio->SetMarkerColor(1);
368 ratio->SetMarkerStyle(0);
369 ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
370 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
371
372 ratio->DrawCopy();
373
374 // get average of ratio
375 Float_t sum = 0;
376 for (Int_t i=2; i<=ratioRange; ++i)
377 {
378 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
379 }
380 sum /= ratioRange-1;
381
382 printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
383
384 TLine* line = new TLine(0, 1, displayRange, 1);
385 line->SetLineWidth(2);
386 line->Draw();
387
388 line = new TLine(0, 1.1, displayRange, 1.1);
389 line->SetLineWidth(2);
390 line->SetLineStyle(2);
391 line->Draw();
392 line = new TLine(0, 0.9, displayRange, 0.9);
393 line->SetLineWidth(2);
394 line->SetLineStyle(2);
395 line->Draw();
396
397 canvas->Modified();
398
399 canvas->SaveAs(epsName);
400
401 return canvas;
402}
403
404TCanvas* Draw2ResultRatios(TH1* mcHist, TH1* result1, TH1* result2, TH1* ratio2, TH1* ratio3, TString epsName)
405{
406 // draws the 3 plots in the upper plot
407 // draws the ratio between result1 and result2 in the lower plot
408 // also draws ratio2 and ratio3 in the lower plot, uses their name for the legend
409
f3eb27f6 410 // normalize unfolded result to mc hist
411 result1->Scale(1.0 / result1->Integral(2, 200));
412 result1->Scale(mcHist->Integral(2, 200));
413 result2->Scale(1.0 / result2->Integral(2, 200));
414 result2->Scale(mcHist->Integral(2, 200));
415
416 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
417 canvas->Range(0, 0, 1, 1);
418
419 TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
420 pad1->Draw();
421
422 TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
423 pad2->Draw();
424
425 pad1->SetRightMargin(0.05);
426 pad2->SetRightMargin(0.05);
427
428 // no border between them
429 pad1->SetBottomMargin(0);
430 pad2->SetTopMargin(0);
431
432 pad1->cd();
69b09e3b 433 gPad->SetGridx();
434 gPad->SetGridy();
f3eb27f6 435
436 mcHist->GetXaxis()->SetLabelSize(0.06);
437 mcHist->GetYaxis()->SetLabelSize(0.06);
438 mcHist->GetXaxis()->SetTitleSize(0.06);
439 mcHist->GetYaxis()->SetTitleSize(0.06);
440 mcHist->GetYaxis()->SetTitleOffset(0.6);
441
442 mcHist->GetXaxis()->SetRangeUser(0, displayRange);
443
69b09e3b 444 mcHist->SetTitle(";True multiplicity;Entries");
f3eb27f6 445 mcHist->SetStats(kFALSE);
446
447 mcHist->DrawCopy("HIST E");
448 gPad->SetLogy();
449
450 result1->SetLineColor(2);
69b09e3b 451 result1->SetMarkerColor(2);
452 result1->SetMarkerStyle(24);
f3eb27f6 453 result1->DrawCopy("SAME HISTE");
454
455 result2->SetLineColor(4);
69b09e3b 456 result2->SetMarkerColor(4);
457 result2->SetMarkerStyle(5);
f3eb27f6 458 result2->DrawCopy("SAME HISTE");
459
460 TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
69b09e3b 461 legend->AddEntry(mcHist, "True distribution");
462 legend->AddEntry(result1, "Unfolded distribution (5%)", "P");
463 legend->AddEntry(result2, "Unfolded distribution (normal)", "P");
f3eb27f6 464 legend->SetFillColor(0);
69b09e3b 465 legend->SetTextSize(0.06);
f3eb27f6 466 legend->Draw();
467
468 pad2->cd();
469 pad2->SetBottomMargin(0.15);
69b09e3b 470 //gPad->SetGridx();
471 //gPad->SetGridy();
f3eb27f6 472
473 result1->GetXaxis()->SetLabelSize(0.06);
474 result1->GetYaxis()->SetLabelSize(0.06);
475 result1->GetXaxis()->SetTitleSize(0.06);
476 result1->GetYaxis()->SetTitleSize(0.06);
477 result1->GetYaxis()->SetTitleOffset(0.6);
478
479 result1->GetXaxis()->SetRangeUser(0, displayRange);
480
69b09e3b 481 result1->SetTitle(Form(";%s;Entries", GetMultLabel()));
f3eb27f6 482 result1->SetStats(kFALSE);
483
484 // calculate ratio
485 result1->Sumw2();
486 TH1* ratio = (TH1*) result1->Clone("ratio");
487 result2->Sumw2();
488 ratio->Divide(ratio, result2, 1, 1, "");
69b09e3b 489 ratio->SetLineColor(1);
490 ratio->SetMarkerColor(1);
491 ratio->SetMarkerStyle(24);
492 ratio->GetYaxis()->SetTitle("Ratio (change / normal)");
f3eb27f6 493 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
494
69b09e3b 495 ratio2->SetLineColor(2);
496 ratio2->SetMarkerColor(2);
497 ratio2->SetMarkerStyle(25);
498
499 ratio3->SetLineColor(4);
500 ratio3->SetMarkerColor(4);
501 ratio3->SetMarkerStyle(26);
502
f3eb27f6 503 ratio->DrawCopy();
69b09e3b 504 ratio2->DrawCopy("SAME");
505 ratio3->DrawCopy("SAME");
506
507 legend2 = new TLegend(0.3, 0.8, 0.8, 0.95);
508 legend2->SetNColumns(3);
509 legend2->SetFillColor(0);
510 legend2->SetTextSize(0.06);
511 legend2->AddEntry(ratio, "5% change", "P");
512 legend2->AddEntry(ratio2, "2% change", "P");
513 legend2->AddEntry(ratio3, "1% change", "P");
514 legend2->Draw();
f3eb27f6 515
516 // get average of ratio
517 Float_t sum = 0;
518 for (Int_t i=2; i<=ratioRange; ++i)
519 {
520 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
521 }
522 sum /= ratioRange-1;
523
524 printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
525
526 TLine* line = new TLine(0, 1, displayRange, 1);
527 line->SetLineWidth(2);
528 line->Draw();
529
530 line = new TLine(0, 1.1, displayRange, 1.1);
531 line->SetLineWidth(2);
532 line->SetLineStyle(2);
533 line->Draw();
534 line = new TLine(0, 0.9, displayRange, 0.9);
535 line->SetLineWidth(2);
536 line->SetLineStyle(2);
537 line->Draw();
538
539 canvas->Modified();
540
541 canvas->SaveAs(epsName);
542
543 return canvas;
544}
545
546TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE)
547{
548 // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
549 // a systematic effect
550
551 // normalize results
552 result->Scale(1.0 / result->Integral(2, 200));
553
69b09e3b 554 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 500);
f3eb27f6 555 canvas->SetTopMargin(0.05);
556 canvas->SetRightMargin(0.05);
69b09e3b 557 canvas->SetGridx();
558 canvas->SetGridy();
f3eb27f6 559
560 result->GetXaxis()->SetRangeUser(0, displayRange);
561 result->GetYaxis()->SetRangeUser(0.55, 1.45);
562 result->SetStats(kFALSE);
563
564 // to get the axis how we want it
565 TH1* dummy = (TH1*) result->Clone("dummy");
566 dummy->Reset();
69b09e3b 567 dummy->SetTitle(Form(";%s;Ratio", GetMultLabel()));
f3eb27f6 568 dummy->DrawCopy();
569 delete dummy;
570
571 Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
572
69b09e3b 573 TLegend* legend = new TLegend(0.2, 0.7, 0.7, 0.93);
f3eb27f6 574 legend->SetFillColor(0);
69b09e3b 575 legend->SetTextSize(0.04);
576 if (nResultSyst > 6)
577 legend->SetNColumns(2);
f3eb27f6 578
579 for (Int_t n=0; n<nResultSyst; ++n)
580 {
581 resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
582
583 // calculate ratio
584 TH1* ratio = (TH1*) result->Clone("ratio");
585 ratio->Divide(ratio, resultSyst[n], 1, 1, "");
69b09e3b 586 ratio->GetXaxis()->SetRangeUser(0, displayRange);
f3eb27f6 587
588 if (firstMarker)
589 ratio->SetMarkerStyle(5);
590
591 ratio->SetLineColor(colors[n / 2]);
592 if ((n % 2))
593 ratio->SetLineStyle(2);
594
595 TString drawStr("SAME HIST");
596 if (n == 0 && firstMarker)
597 drawStr = "SAME P";
598 if (errors)
599 drawStr += " E";
600
601 ratio->DrawCopy(drawStr);
602
603 if (legendStrings && legendStrings[n])
69b09e3b 604 legend->AddEntry(ratio, legendStrings[n], "L");
f3eb27f6 605
606 // get average of ratio
607 Float_t sum = 0;
608 for (Int_t i=2; i<=ratioRange; ++i)
609 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
610 sum /= ratioRange-1;
611
612 printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
613 }
614
615 if (legendStrings)
616 legend->Draw();
617
69b09e3b 618 TLine* line = new TLine(-0.5, 1, displayRange, 1);
f3eb27f6 619 line->SetLineWidth(2);
620 line->Draw();
621
69b09e3b 622 line = new TLine(-0.5, 1.1, displayRange, 1.1);
f3eb27f6 623 line->SetLineWidth(2);
624 line->SetLineStyle(2);
625 line->Draw();
69b09e3b 626 line = new TLine(-0.5, 0.9, displayRange, 0.9);
f3eb27f6 627 line->SetLineWidth(2);
628 line->SetLineStyle(2);
629 line->Draw();
630
631 canvas->SaveAs(epsName);
632 canvas->SaveAs(Form("%s.gif", epsName.Data()));
633
634 return canvas;
635}
636
637TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE)
638{
639 // draws the ratios of each mc to the corresponding result
640
641 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
642 canvas->SetRightMargin(0.05);
643 canvas->SetTopMargin(0.05);
644
645 for (Int_t n=0; n<nResultSyst; ++n)
646 {
647 // normalize
648 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
649 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
650
651 result[n]->GetXaxis()->SetRangeUser(0, displayRange);
652 result[n]->SetStats(kFALSE);
653
654 // calculate ratio
655 TH1* ratio = (TH1*) result[n]->Clone("ratio");
656 ratio->Divide(mc[n], ratio, 1, 1, "B");
657
658 // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis...
659 ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0);
660
661 if (smooth)
662 Smooth(ratio);
663
664 ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : "")));
665 ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
666
667 if (dashed)
668 {
669 ratio->SetLineColor((n/2)+1);
670 ratio->SetLineStyle((n%2)+1);
671 }
672 else
673 ratio->SetLineColor(n+1);
674
675 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
676
677 // get average of ratio
678 Float_t sum = 0;
679 for (Int_t i=2; i<=ratioRange; ++i)
680 sum += TMath::Abs(ratio->GetBinContent(i) - 1);
681 sum /= ratioRange-1;
682
683 printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
684 }
685
686 TLine* line = new TLine(0, 1, displayRange, 1);
687 line->SetLineWidth(2);
688 line->Draw();
689
690 line = new TLine(0, 1.1, displayRange, 1.1);
691 line->SetLineWidth(2);
692 line->SetLineStyle(2);
693 line->Draw();
694 line = new TLine(0, 0.9, displayRange, 0.9);
695 line->SetLineWidth(2);
696 line->SetLineStyle(2);
697 line->Draw();
698
699 canvas->Modified();
700
701 canvas->SaveAs(epsName);
702 canvas->SaveAs(Form("%s.gif", epsName.Data()));
703
704 return canvas;
705}
706
707TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
708{
709 // draws the ratios of each mc to the corresponding result
710 // deducts from each ratio the ratio of mcBase / resultBase
711
712 // normalize
713 resultBase->Scale(1.0 / resultBase->Integral(2, 200));
714 mcBase->Scale(1.0 / mcBase->Integral(2, 200));
715
716 // calculate ratio
717 TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
718 ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
719
720 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
721 canvas->SetRightMargin(0.05);
722 canvas->SetTopMargin(0.05);
723
724 for (Int_t n=0; n<nResultSyst; ++n)
725 {
726 // normalize
727 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
728 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
729
730 result[n]->GetXaxis()->SetRangeUser(0, displayRange);
731 result[n]->SetStats(kFALSE);
732
733 // calculate ratio
734 TH1* ratio = (TH1*) result[n]->Clone("ratio");
735 ratio->Divide(mc[n], ratio, 1, 1, "B");
736 ratio->Add(ratioBase, -1);
737
738 ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
739 ratio->GetYaxis()->SetRangeUser(-1, 1);
740 ratio->SetLineColor(n+1);
741 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
742
743 // get average of ratio
744 Float_t sum = 0;
745 for (Int_t i=2; i<=ratioRange; ++i)
746 sum += TMath::Abs(ratio->GetBinContent(i));
747 sum /= ratioRange-1;
748
749 printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, ratioRange, sum);
750 }
751
752 TLine* line = new TLine(0, 0, displayRange, 0);
753 line->SetLineWidth(2);
754 line->Draw();
755
756 line = new TLine(0, 0.1, displayRange, 0.1);
757 line->SetLineWidth(2);
758 line->SetLineStyle(2);
759 line->Draw();
760 line = new TLine(0, -0.1, displayRange, -0.1);
761 line->SetLineWidth(2);
762 line->SetLineStyle(2);
763 line->Draw();
764
765 canvas->Modified();
766
767 canvas->SaveAs(epsName);
768 canvas->SaveAs(Form("%s.gif", epsName.Data()));
769
770 return canvas;
771}
772
773TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
774{
775 // draws the ratios of each mc to the corresponding result
776 // deducts from each ratio the ratio of mcBase / resultBase
777 // smoothens the ratios by a sliding window
778
779 // normalize
780 resultBase->Scale(1.0 / resultBase->Integral(2, 200));
781 mcBase->Scale(1.0 / mcBase->Integral(2, 200));
782
783 // calculate ratio
784 TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
785 ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
786
787 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
788 canvas->SetRightMargin(0.05);
789 canvas->SetTopMargin(0.05);
790
791 for (Int_t n=0; n<nResultSyst; ++n)
792 {
793 // normalize
794 result[n]->Scale(1.0 / result[n]->Integral(2, 200));
795 mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
796
797 result[n]->GetXaxis()->SetRangeUser(0, displayRange);
798 result[n]->SetStats(kFALSE);
799
800 // calculate ratio
801 TH1* ratio = (TH1*) result[n]->Clone("ratio");
802 ratio->Divide(mc[n], ratio, 1, 1, "B");
803 ratio->Add(ratioBase, -1);
804
805 //new TCanvas; ratio->DrawCopy();
806 // clear 0 bin
807 ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0);
808
809 Smooth(ratio);
810
811 //ratio->SetLineColor(1); ratio->DrawCopy("SAME");
812
813 canvas->cd();
814 ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)");
815 ratio->GetYaxis()->SetRangeUser(-0.3, 0.3);
816 ratio->SetLineColor((n / 2)+1);
817 ratio->SetLineStyle((n % 2)+1);
818 ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
819
820 // get average of ratio
821 Float_t sum = 0;
822 for (Int_t i=2; i<=150; ++i)
823 sum += TMath::Abs(ratio->GetBinContent(i));
824 sum /= 149;
825
826 printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
827 }
828
829 TLine* line = new TLine(0, 0, displayRange, 0);
830 line->SetLineWidth(2);
831 line->Draw();
832
833 line = new TLine(0, 0.1, displayRange, 0.1);
834 line->SetLineWidth(2);
835 line->SetLineStyle(2);
836 line->Draw();
837 line = new TLine(0, -0.1, displayRange, -0.1);
838 line->SetLineWidth(2);
839 line->SetLineStyle(2);
840 line->Draw();
841
842 canvas->Modified();
843
844 canvas->SaveAs(epsName);
845 canvas->SaveAs(Form("%s.gif", epsName.Data()));
846
847 return canvas;
848}
849
69b09e3b 850void DrawResiduals(const char* fileName, const char* epsName)
f3eb27f6 851{
69b09e3b 852 loadlibs();
853
854 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
855 TFile::Open(fileName);
856 mult->LoadHistograms("Multiplicity");
857
858 TH1* measured = mult->GetMultiplicityESD(etaRange)->ProjectionY("myesd", 1, 1);
859 TH1* unfoldedFolded = mult->CalculateMultiplicityESD(mult->GetMultiplicityESDCorrected(etaRange), etaRange)->ProjectionY("myfolded", 1, 1);
860
f3eb27f6 861 // normalize
69b09e3b 862 unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, displayRange+1));
863 unfoldedFolded->Scale(measured->Integral(2, displayRange+1));
f3eb27f6 864
865 TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
866 canvas->Range(0, 0, 1, 1);
867
69b09e3b 868 TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 1, 1);
f3eb27f6 869 pad1->Draw();
870 pad1->SetGridx();
871 pad1->SetGridy();
872
69b09e3b 873 TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 1, 0.5);
f3eb27f6 874 pad2->Draw();
875 pad2->SetGridx();
876 pad2->SetGridy();
877
878 TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
879 pad3->SetGridx();
880 pad3->SetGridy();
881 pad3->SetRightMargin(0.05);
882 pad3->SetTopMargin(0.05);
883 pad3->Draw();
884
885 pad1->SetRightMargin(0.05);
886 pad2->SetRightMargin(0.05);
887
888 // no border between them
889 pad1->SetBottomMargin(0);
890 pad2->SetTopMargin(0);
891
892 pad1->cd();
893
894 measured->GetXaxis()->SetLabelSize(0.06);
895 measured->GetYaxis()->SetLabelSize(0.06);
896 measured->GetXaxis()->SetTitleSize(0.06);
897 measured->GetYaxis()->SetTitleSize(0.06);
898 measured->GetYaxis()->SetTitleOffset(0.6);
899
69b09e3b 900 measured->GetXaxis()->SetRangeUser(0, displayRange);
f3eb27f6 901
69b09e3b 902 measured->SetTitle(Form(";%s;Entries", GetMultLabel(etaRange, kFALSE)));
f3eb27f6 903 measured->SetStats(kFALSE);
904
905 measured->DrawCopy("HIST");
906 gPad->SetLogy();
907
908 unfoldedFolded->SetMarkerStyle(5);
909 unfoldedFolded->SetMarkerColor(2);
69b09e3b 910 unfoldedFolded->SetLineColor(2);
911 unfoldedFolded->DrawCopy("SAME PHIST");
f3eb27f6 912
913 TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
69b09e3b 914 legend->AddEntry(measured, "Measured distribution", "L");
915 legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution", "P");
f3eb27f6 916 legend->SetFillColor(0);
69b09e3b 917 legend->SetTextSize(0.06);
f3eb27f6 918 legend->Draw();
919
920 pad2->cd();
921 pad2->SetBottomMargin(0.15);
922
923 // calculate ratio
924 measured->Sumw2();
925 TH1* residual = (TH1*) measured->Clone("residual");
926 unfoldedFolded->Sumw2();
927
928 residual->Add(unfoldedFolded, -1);
929
930 // projection
69b09e3b 931 TH1* residualHist = new TH1F("residualHist", ";", 11, -3, 3);
f3eb27f6 932
69b09e3b 933 for (Int_t i=1; i<=displayRange+1; ++i)
f3eb27f6 934 {
935 if (measured->GetBinError(i) > 0)
936 {
937 residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
938 residual->SetBinError(i, 1);
939
940 residualHist->Fill(residual->GetBinContent(i));
941 }
942 else
943 {
944 residual->SetBinContent(i, 0);
945 residual->SetBinError(i, 0);
946 }
947 }
948
69b09e3b 949 residual->GetYaxis()->SetTitle("Residuals: (1/e) (M - R #otimes U)");
f3eb27f6 950 residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
951 residual->DrawCopy();
952
69b09e3b 953 TLine* line = new TLine(-0.5, 0, displayRange + 0.5, 0);
f3eb27f6 954 line->SetLineWidth(2);
955 line->Draw();
956
957 pad3->cd();
958 residualHist->SetStats(kFALSE);
69b09e3b 959 residualHist->SetLabelSize(0.08, "xy");
f3eb27f6 960 residualHist->Fit("gaus");
961 residualHist->Draw();
962
963 canvas->Modified();
964 canvas->SaveAs(canvas->GetName());
965
966 //const char* epsName2 = "proj.eps";
967 //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
968 //canvas->SetGridx();
969 //canvas->SetGridy();
970
971 //canvas->SaveAs(canvas->GetName());
972}
973
69b09e3b 974void chi2FluctuationResult()
f3eb27f6 975{
69b09e3b 976 loadlibs();
f3eb27f6 977
69b09e3b 978 TFile::Open("chi2_noregularization.root");
f3eb27f6 979 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
980 mult->LoadHistograms("Multiplicity");
981
69b09e3b 982 TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
f3eb27f6 983
69b09e3b 984 mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
f3eb27f6 985
69b09e3b 986 TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_1");
987 ((TH1*) canvas->FindObject("proj"))->GetXaxis()->SetRangeUser(0, displayRange);
988 ((TH1*) canvas->FindObject("fCurrentESD"))->GetXaxis()->SetRangeUser(0, displayRange);
989 //((TH1*) canvas->FindObject("proj"))->GetXaxis()->SetTitle(GetMultTitle());
990 //((TH1*) canvas->FindObject("fCurrentESD"))->GetXaxis()->->SetTitle(GetMultTitle(etaRange, kFALSE));
991 canvas->SaveAs("chi2FluctuationResult.eps");
f3eb27f6 992}
993
69b09e3b 994void DrawUnfolded(const char* fileName, const char* eps)
f3eb27f6 995{
69b09e3b 996 loadlibs();
997
998 TFile::Open(fileName);
f3eb27f6 999
f3eb27f6 1000 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1001 mult->LoadHistograms("Multiplicity");
1002
69b09e3b 1003 TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, mult->GetMultiplicityVtx(etaRange)->GetNbinsX());
1004 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
f3eb27f6 1005
69b09e3b 1006 DrawResultRatio(mcHist, result, eps);
f3eb27f6 1007}
1008
69b09e3b 1009void minimizationInfluenceAlpha()
f3eb27f6 1010{
69b09e3b 1011 loadlibs();
f3eb27f6 1012
1013 TFile::Open(measuredFile);
1014 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1015 mult2->LoadHistograms("Multiplicity");
1016
69b09e3b 1017 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, mult2->GetMultiplicityVtx(etaRange)->GetNbinsX());
1018 mcHist->Scale(1.0 / mcHist->Integral());
1019 mcHist->SetStats(kFALSE);
1020 mcHist->SetTitle(";True multiplicity n in |#eta| < 1.0;P(n)");
f3eb27f6 1021
69b09e3b 1022 TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 350);
1023 canvas->Divide(3, 1, 0.005);
1024
1025 TFile::Open("chi2compare-influencealpha/EvaluateChi2Method.root");
1026
1027 TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_3.162278");
1028 TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_07_2_10000.000000");
1029 TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_13_2_10000000.000000");
1030
1031 /*mcHist->Rebin(2); mcHist->Scale(0.5);
1032 hist1->Rebin(2); hist1->Scale(0.5);
1033 hist2->Rebin(2); hist2->Scale(0.5);
1034 hist3->Rebin(2); hist3->Scale(0.5);*/
1035
1036 mcHist->GetXaxis()->SetRangeUser(0, displayRange);
1037 mcHist->SetLabelSize(0.06, "xy");
1038 mcHist->SetTitleSize(0.06, "xy");
1039 mcHist->GetYaxis()->SetTitleOffset(1.5);
1040
1041 canvas->cd(1);
1042
1043 gPad->SetLogy();
1044 gPad->SetRightMargin(0.03);
1045 gPad->SetLeftMargin(0.19);
1046 gPad->SetTopMargin(0.05);
1047 gPad->SetBottomMargin(0.13);
1048 gPad->SetGridx();
1049 gPad->SetGridy();
1050 mcHist->Draw();
1051 hist1->SetMarkerStyle(5);
1052 hist1->SetMarkerColor(2);
1053 hist1->Draw("SAME PE");
1054
1055 canvas->cd(2);
1056 gPad->SetRightMargin(0.03);
1057 gPad->SetLeftMargin(0.19);
1058 gPad->SetTopMargin(0.05);
1059 gPad->SetBottomMargin(0.13);
1060 gPad->SetLogy();
1061 gPad->SetGridx();
1062 gPad->SetGridy();
1063 mcHist->Draw();
1064 hist2->SetMarkerStyle(5);
1065 hist2->SetMarkerColor(2);
1066 hist2->Draw("SAME PE");
f3eb27f6 1067
69b09e3b 1068 canvas->cd(3);
1069 gPad->SetRightMargin(0.03);
1070 gPad->SetLeftMargin(0.19);
1071 gPad->SetTopMargin(0.05);
1072 gPad->SetBottomMargin(0.13);
1073 gPad->SetLogy();
1074 gPad->SetGridx();
1075 gPad->SetGridy();
1076 mcHist->Draw();
1077 hist3->SetMarkerStyle(5);
1078 hist3->SetMarkerColor(2);
1079 hist3->Draw("SAME PE");
f3eb27f6 1080
69b09e3b 1081 canvas->SaveAs("minimizationInfluenceAlpha.eps");
f3eb27f6 1082}
1083
69b09e3b 1084void NBDFit()
f3eb27f6 1085{
f3eb27f6 1086 gSystem->Load("libPWG0base");
1087
69b09e3b 1088 TFile::Open(correctionFile);
f3eb27f6 1089 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1090 mult->LoadHistograms("Multiplicity");
1091
1092 TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
1093 fCurrentESD->Sumw2();
1094 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1095
1096 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])");
1097 func->SetParNames("scaling", "averagen", "k");
1098 func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
1099 func->SetParLimits(1, 0.001, 1000);
1100 func->SetParLimits(2, 0.001, 1000);
1101 func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
1102
1103 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
1104 lognormal->SetParNames("scaling", "mean", "sigma");
1105 lognormal->SetParameters(1, 1, 1);
1106 lognormal->SetParLimits(0, 0, 10);
1107 lognormal->SetParLimits(1, 0, 100);
1108 lognormal->SetParLimits(2, 1e-3, 10);
1109
1110 TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
1111 fCurrentESD->SetStats(kFALSE);
1112 fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
1113 fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
1114 fCurrentESD->Draw("HIST");
1115 fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1116 fCurrentESD->Fit(func, "W0", "", 0, 50);
1117 func->SetRange(0, 100);
1118 func->Draw("SAME");
1119 printf("chi2 = %f\n", func->GetChisquare());
1120
1121 fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
1122 lognormal->SetLineColor(2);
1123 lognormal->SetLineStyle(2);
1124 lognormal->SetRange(0, 100);
1125 lognormal->Draw("SAME");
1126
1127 canvas->SaveAs("NBDFit.eps");
1128}
1129
f3eb27f6 1130void StartingConditions()
1131{
1132 // data generated by runMultiplicitySelector.C StartingConditions
1133
1134 const char* name = "StartingConditions";
1135
1136 TFile* file = TFile::Open(Form("%s.root", name));
1137
69b09e3b 1138 TCanvas* canvas = new TCanvas(name, name, 1200, 600);
f3eb27f6 1139 canvas->Divide(2, 1);
1140
1141 TH1* mc = (TH1*) file->Get("mc");
1142 mc->Sumw2();
69b09e3b 1143 mc->Scale(1.0 / mc->Integral(2, displayRange));
f3eb27f6 1144
1145 //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
1146
69b09e3b 1147 TLegend* legend = new TLegend(0.6, 0.7, 0.99, 0.99);
f3eb27f6 1148 legend->SetFillColor(0);
69b09e3b 1149 legend->SetTextSize(0.04);
1150
f3eb27f6 1151 const char* names[] = { "True", "Measured 1", "Measured 2", "Measured 3", "NBD", "Flat" };
1152
1153 for (Int_t i=0; i<6; ++i)
1154 {
1155 Int_t id = i;
1156 if (id > 2)
1157 id += 2;
1158
1159 TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
1160 TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
69b09e3b 1161
1162 chi2Result->Scale(1.0 / chi2Result->Integral(2, displayRange));
1163 bayesResult->Scale(1.0 / bayesResult->Integral(2, displayRange));
f3eb27f6 1164
69b09e3b 1165 chi2Result->Divide(mc, chi2Result, 1, 1, "");
1166 bayesResult->Divide(mc, bayesResult, 1, 1, "");
f3eb27f6 1167
69b09e3b 1168 chi2Result->SetTitle(Form("a) #chi^{2}-minimization;%s;MC / unfolded", GetMultLabel()));
1169 chi2Result->GetXaxis()->SetRangeUser(0, displayRange);
1170 chi2Result->GetYaxis()->SetRangeUser(0.7, 1.3);
1171 chi2Result->GetYaxis()->SetTitleOffset(1.7);
f3eb27f6 1172 //chi2Result->SetMarkerStyle(marker[i]);
1173 chi2Result->SetLineColor(i+1);
1174 chi2Result->SetMarkerColor(i+1);
1175 chi2Result->SetStats(kFALSE);
1176
69b09e3b 1177 bayesResult->SetTitle(Form("b) Bayesian unfolding;%s;MC / unfolded", GetMultLabel()));
1178 bayesResult->GetXaxis()->SetRangeUser(0, displayRange);
1179 bayesResult->GetYaxis()->SetRangeUser(0.7, 1.3);
1180 bayesResult->GetYaxis()->SetTitleOffset(1.7);
f3eb27f6 1181 bayesResult->SetStats(kFALSE);
1182 //bayesResult->SetLineColor(2);
1183 bayesResult->SetLineColor(i+1);
1184
1185 canvas->cd(1);
69b09e3b 1186 gPad->SetRightMargin(0.05);
f3eb27f6 1187 gPad->SetLeftMargin(0.12);
69b09e3b 1188 gPad->SetGridx();
1189 gPad->SetGridy();
f3eb27f6 1190 chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1191
1192 canvas->cd(2);
69b09e3b 1193 gPad->SetRightMargin(0.05);
f3eb27f6 1194 gPad->SetLeftMargin(0.12);
69b09e3b 1195 gPad->SetGridx();
1196 gPad->SetGridy();
f3eb27f6 1197 bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1198
1199 //TLine* line = new TLine(0, 1, 150, 1);
1200 //line->Draw();
1201
1202 legend->AddEntry(chi2Result, names[i]);
1203 }
1204
1205 canvas->cd(1);
1206 legend->Draw();
69b09e3b 1207 canvas->cd(2);
1208 legend->Draw();
f3eb27f6 1209 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1210}
1211
1212void StatisticsPlot()
1213{
1214 const char* name = "StatisticsPlot";
1215
1216 TFile* file = TFile::Open(Form("%s.root", name));
1217
1218 TCanvas* canvas = new TCanvas(name, name, 600, 400);
1219
1220 TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
1221 fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1222 fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1223 fitResultsChi2->Draw("AP");
1224
1225 TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1226 fitResultsChi2->Fit(f);
1227
1228 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1229
1230 TH1* mc[5];
1231 TH1* result[5];
1232
1233 const char* plotname = "chi2Result";
1234
1235 name = "StatisticsPlotRatios";
1236 canvas = new TCanvas(name, name, 600, 400);
1237
1238 for (Int_t i=0; i<5; ++i)
1239 {
1240 mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1241 result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1242
1243 result[i]->SetLineColor(i+1);
1244 result[i]->Draw(((i == 0) ? "" : "SAME"));
1245 }
1246}
1247
69b09e3b 1248void Draw2Unfolded(const char* file1, const char* file2, const char* output)
f3eb27f6 1249{
69b09e3b 1250 loadlibs();
1251
1252 TFile::Open(file1);
f3eb27f6 1253 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1254 mult->LoadHistograms("Multiplicity");
1255
69b09e3b 1256 // result with systematic effect
1257 TFile::Open(file2);
f3eb27f6 1258 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1259 mult2->LoadHistograms("Multiplicity");
1260
69b09e3b 1261 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1262 TH1* result1 = (TH1*) mult2->GetMultiplicityESDCorrected(etaRange)->Clone("result1"); // from file2 (with syst)
1263 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2"); // from file1 (without syst)
f3eb27f6 1264
69b09e3b 1265 DrawResultRatio(mcHist, result1, "tmp1.eps");
1266 DrawResultRatio(mcHist, result2, "tmp2.eps");
1267 Draw2ResultRatio(mcHist, result1, result2, output);
1268}
f3eb27f6 1269
69b09e3b 1270void PythiaPhojet()
1271{
1272 loadlibs();
1273
1274 displayRange = 55;
1275 Draw2Unfolded("self.root", "pythia.root", "test.eps");
1276
1277 canvas = (TCanvas*) gROOT->GetListOfCanvases()->Last();
1278 pad1 = (TCanvas*)canvas->GetListOfPrimitives()->First();
1279 pad2 = (TCanvas*)canvas->GetListOfPrimitives()->Last();
1280 legend = (TLegend*)pad1->GetListOfPrimitives()->Last();
1281
1282 ((TH1*)pad2->GetListOfPrimitives()->At(1))->GetYaxis()->SetTitle("Ratio (Pythia / Phojet)");
1283 ((TLegendEntry*)legend->GetListOfPrimitives()->At(1))->SetLabel("Unfolded distribution (Pythia)");
1284 ((TLegendEntry*)legend->GetListOfPrimitives()->At(2))->SetLabel("Unfolded distribution (Phojet)");
1285 canvas->SaveAs("PythiaPhojet.eps");
1286}
f3eb27f6 1287
69b09e3b 1288void Misalignment()
1289{
1290 loadlibs();
1291
1292 Draw2Unfolded("chi2_ideal.root", "chi2_misaligned.root", "test.eps");
1293
1294 canvas = (TCanvas*) gROOT->GetListOfCanvases()->Last();
1295 pad1 = (TCanvas*)canvas->GetListOfPrimitives()->First();
1296 pad2 = (TCanvas*)canvas->GetListOfPrimitives()->Last();
1297 legend = (TLegend*)pad1->GetListOfPrimitives()->Last();
1298
1299 ((TH1*)pad2->GetListOfPrimitives()->At(1))->GetYaxis()->SetTitle("Ratio (misaligned / realigned)");
1300 ((TLegendEntry*)legend->GetListOfPrimitives()->At(1))->SetLabel("Unfolded distribution (misaligned)");
1301 ((TLegendEntry*)legend->GetListOfPrimitives()->At(2))->SetLabel("Unfolded distribution (realigned)");
1302 canvas->SaveAs("SystematicMisalignment.eps");
1303}
f3eb27f6 1304
69b09e3b 1305void SystematicLowEfficiency()
1306{
1307 Draw2Unfolded("chi2.root", "chi2_loweff_5.root", "SystematicLowEfficiency.eps");
1308}
f3eb27f6 1309
69b09e3b 1310void SystematicLowEfficiency2()
1311{
1312 loadlibs();
1313
1314 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("chi2.root");
f3eb27f6 1315 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
69b09e3b 1316 TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1317 result2->Scale(1.0 / result2->Integral(2, displayRange));
1318 result2->Scale(mcHist->Integral(2, displayRange));
1319
1320 mult = AliMultiplicityCorrection::Open("chi2_loweff_5.root");
1321 TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
1322
1323 mult = AliMultiplicityCorrection::Open("chi2_loweff_2.root");
1324 TH1* ratio2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("ratio2");
1325 ratio2->Scale(1.0 / ratio2->Integral(2, displayRange));
1326 ratio2->Scale(mcHist->Integral(2, displayRange));
1327 ratio2->Divide(result2);
1328
1329 mult = AliMultiplicityCorrection::Open("chi2_loweff_1.root");
1330 TH1* ratio3 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("ratio3");
1331 ratio3->Scale(1.0 / ratio3->Integral(2, displayRange));
1332 ratio3->Scale(mcHist->Integral(2, displayRange));
1333 ratio3->Divide(result2);
1334
1335 Draw2ResultRatios(mcHist, result1, result2, ratio2, ratio3, "SystematicLowEfficiency2.eps");
f3eb27f6 1336}
1337
1338void SystematicMisalignment()
1339{
1340 gSystem->Load("libPWG0base");
1341
1342 TFile::Open(correctionFile);
1343 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1344 mult->LoadHistograms("Multiplicity");
1345
1346 TFile::Open("multiplicityMC_100k_fullmis.root");
1347 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1348 mult2->LoadHistograms("Multiplicity");
1349
1350 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1351 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1352 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1353
1354 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1355 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1356
1357 DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1358}
1359
1360void SystematicMisalignmentTPC()
1361{
1362 gSystem->Load("libPWG0base");
1363
1364 SetTPC();
1365
1366 TFile::Open(correctionFile);
1367 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1368 mult->LoadHistograms("Multiplicity");
1369
1370 TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
1371 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1372 mult2->LoadHistograms("Multiplicity");
1373
1374 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1375 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1376 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1377
1378 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1379 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1380
1381 DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
1382}
1383
69b09e3b 1384void LowMomentumEffectSPD()
1385{
1386 // this function increases/reduces the correction as function of pt between 0 and 0.2 gev by +-50% to 0% and checks the effect on the overall correction factor
1387 // only a normal acceptance region is considered to not get a bias by the edges
1388
1389 loadlibs();
1390 TFile::Open("multiplicity.root");
1391
1392
1393 AliCorrection* correction[8];
1394 Float_t values[3];
1395
1396 for (Int_t loop=0; loop<3; loop++)
1397 {
1398 Float_t sumGen = 0;
1399 Float_t sumMeas = 0;
1400
1401 Printf("loop %d", loop);
1402 for (Int_t i=0; i<4; ++i)
1403 {
1404 Printf("correction %d", i);
1405
1406 TString name; name.Form("correction_%d", i);
1407 correction[i] = new AliCorrection(name, name);
1408 correction[i]->LoadHistograms();
1409
1410 TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1411 TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1412
1413 Float_t vtxRange = 5.9;
1414 gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1415 meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1416
1417 Float_t etaRange = 0.99;
1418 gene->GetYaxis()->SetRangeUser(-etaRange, etaRange);
1419 meas->GetYaxis()->SetRangeUser(-etaRange, etaRange);
1420
1421 TH1* genePt = gene->Project3D(Form("z_%d", i));
1422 TH1* measPt = meas->Project3D(Form("z_%d", i));
1423
1424 if (loop > 0)
1425 {
1426 for (Int_t x=1; x<=genePt->GetNbinsX(); x++)
1427 {
1428 Float_t pt = genePt->GetXaxis()->GetBinCenter(x);
1429 //Printf("%f", pt);
1430 if (pt < 0.2)
1431 {
1432 Float_t factor = 1;
1433 if (loop == 1)
1434 factor = 1.5 - pt / 0.2 * 0.5;
1435 if (loop == 2)
1436 factor = 0.5 + pt / 0.2 * 0.5;
1437 //Printf("%f", factor);
1438 genePt->SetBinContent(x, genePt->GetBinContent(x) * factor);
1439 measPt->SetBinContent(x, measPt->GetBinContent(x) * factor);
1440 }
1441 }
1442 }
1443
1444 //new TCanvas; genePt->DrawCopy(); measPt->DrawCopy("SAME");
1445
1446 sumGen += genePt->Integral();
1447 sumMeas += measPt->Integral();
1448
1449 Float_t average = measPt->Integral() / genePt->Integral();
1450
1451 Printf("The average efficiency of this correction is %f", average);
1452 }
1453
1454 Float_t average = sumMeas / sumGen;
1455
1456 Printf("The average efficiency of all corrections is %f", average);
1457 values[loop] = average;
1458 }
1459
1460 Printf("relative is %f and %f", values[1] / values[0], values[2] / values[0]);
1461}
1462
1463
1464void EfficiencySpecies(Bool_t addDecayStopped = kFALSE)
f3eb27f6 1465{
5a6310fe 1466 loadlibs();
f3eb27f6 1467
69b09e3b 1468 Int_t marker[] = {24, 25, 26, 27};
1469 Int_t color[] = {1, 2, 4, 3};
f3eb27f6 1470
1471 // SPD TPC
5a6310fe 1472 //const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
69b09e3b 1473 //const char* fileName[] = { "spd/multiplicity.root", "tpc/multiplicity.root" };
1474 const char* fileName[] = { "multiplicity.root", "multiplicity.root" };
1475 Float_t etaRangeArr[] = {0.49, 0.9};
f3eb27f6 1476 const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1477
1478 TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1479 canvas->Divide(2, 1);
1480
69b09e3b 1481 TCanvas* canvas3 = new TCanvas("EfficiencySpecies_comb", "EfficiencySpecies_comb", 600, 600);
1482 gPad->SetGridx();
1483 gPad->SetGridy();
1484 gPad->SetRightMargin(0.05);
1485 gPad->SetTopMargin(0.05);
1486
1487 TLegend* legends[2];
1488
1489 for (Int_t loop=1; loop<2; ++loop)
f3eb27f6 1490 {
1491 Printf("%s", fileName[loop]);
1492
69b09e3b 1493 TCanvas* canvas2 = new TCanvas(Form("EfficiencySpecies_%d", loop), Form("EfficiencySpecies_%d", loop), 600, 600);
1494 gPad->SetGridx();
1495 gPad->SetGridy();
1496 gPad->SetRightMargin(0.05);
1497 gPad->SetTopMargin(0.05);
1498
1499 AliCorrection* correction[8];
f3eb27f6 1500
1501 canvas->cd(loop+1);
1502
1503 gPad->SetGridx();
1504 gPad->SetGridy();
1505 gPad->SetRightMargin(0.05);
f3eb27f6 1506
69b09e3b 1507 TLegend* legend = new TLegend(0.6, 0.4, 0.85, 0.6);
f3eb27f6 1508 legend->SetFillColor(0);
1509 legend->SetEntrySeparation(0.2);
69b09e3b 1510 legend->SetTextSize(gStyle->GetTextSize());
1511
1512 legends[loop] = new TLegend(0.4+loop*0.3, 0.2, 0.6+loop*0.3, 0.5);
1513 legends[loop]->SetFillColor(0);
1514 legends[loop]->SetEntrySeparation(0.2);
1515 legends[loop]->SetTextSize(gStyle->GetTextSize());
1516 legends[loop]->SetHeader((loop == 0) ? "SPD" : "TPC");
f3eb27f6 1517
1518 Float_t below = 0;
1519 Float_t total = 0;
1520
1521 TFile* file = TFile::Open(fileName[loop]);
1522 if (!file)
1523 {
1524 Printf("Could not open %s", fileName[loop]);
1525 return;
1526 }
1527
1528 Float_t sumGen = 0;
1529 Float_t sumMeas = 0;
1530
1531 for (Int_t i=0; i<3; ++i)
1532 {
1533 Printf("correction %d", i);
1534
1535 TString name; name.Form("correction_%d", i);
1536 correction[i] = new AliCorrection(name, name);
1537 correction[i]->LoadHistograms();
69b09e3b 1538
f3eb27f6 1539 TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1540 TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1541
1542 // limit vtx axis
69b09e3b 1543 Float_t vtxRange = 3.9;
1544 gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1545 meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
f3eb27f6 1546
1547 // 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)
1548 /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1549 for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1550 {
1551 gene->SetBinContent(x, 0, z, 0);
1552 gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1553 meas->SetBinContent(x, 0, z, 0);
1554 meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1555 }*/
1556
1557 // limit eta axis
69b09e3b 1558 Float_t etaBegin = -etaRangeArr[loop];
1559 Float_t etaEnd = etaRangeArr[loop];
1560 //etaBegin = 0.01;
1561 //etaEnd = -0.01;
1562 gene->GetYaxis()->SetRangeUser(etaBegin, etaEnd);
1563 meas->GetYaxis()->SetRangeUser(etaBegin, etaEnd);
f3eb27f6 1564
1565 TH1* genePt = gene->Project3D(Form("z_%d", i));
1566 TH1* measPt = meas->Project3D(Form("z_%d", i));
1567
1568 genePt->Sumw2();
1569 measPt->Sumw2();
69b09e3b 1570
1571 for (Int_t x=0; x<=genePt->GetNbinsX()+1; x++)
1572 {
1573 genePt->SetBinError(x, TMath::Sqrt(genePt->GetBinContent(x)));
1574 measPt->SetBinError(x, TMath::Sqrt(measPt->GetBinContent(x)));
1575 }
1576
f3eb27f6 1577 sumGen += genePt->Integral();
1578 sumMeas += measPt->Integral();
1579
1580 TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1581 effPt->Reset();
1582 effPt->Divide(measPt, genePt, 1, 1, "B");
1583
1584 Int_t bin = 0;
1585 for (bin=20; bin>=1; bin--)
1586 {
1587 if (effPt->GetBinContent(bin) < 0.5)
1588 break;
1589 }
1590
1591 Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1592
1593 Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1594 Printf("%.4f of the particles are below that momentum", fraction);
1595
1596 below += genePt->Integral(1, bin);
1597 total += genePt->Integral();
69b09e3b 1598
f3eb27f6 1599 effPt->SetLineColor(color[i]);
1600 effPt->SetMarkerColor(color[i]);
1601 effPt->SetMarkerStyle(marker[i]);
69b09e3b 1602 effPt->SetMarkerSize(2);
f3eb27f6 1603
69b09e3b 1604 effPt->GetXaxis()->SetRangeUser(0, 1);
1605 effPt->GetYaxis()->SetRangeUser(0.001, 1);
f3eb27f6 1606
5a6310fe 1607 effPt->GetXaxis()->SetTitleOffset(1.1);
f3eb27f6 1608 effPt->GetYaxis()->SetTitleOffset(1.2);
69b09e3b 1609
1610 effPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
f3eb27f6 1611
1612 effPt->SetStats(kFALSE);
1613 effPt->SetTitle(titles[loop]);
1614 effPt->GetYaxis()->SetTitle("Efficiency");
1615
69b09e3b 1616 canvas->cd(loop+1);
1617 effPt->DrawCopy((i == 0) ? "" : "SAME");
1618
1619 canvas2->cd();
1620 effPt->SetTitle("");
f3eb27f6 1621 effPt->DrawCopy((i == 0) ? "" : "SAME");
69b09e3b 1622
1623 canvas3->cd();
1624 effPtClone = (TH1*) effPt->Clone("effPtClone");
1625 effPtClone->SetMarkerStyle(marker[i]-4*loop);
1626 effPtClone->DrawCopy((i == 0 && loop == 0) ? "" : "SAME");
1627
1628 legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P");
1629 legends[loop]->AddEntry(effPtClone, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P");
1630 //legend2->AddEntry(effPt, Form("%s %s", (loop == 0) ? "SPD" : "TPC", ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}"))), "P");
f3eb27f6 1631
69b09e3b 1632 if (addDecayStopped)
1633 {
1634 name.Form("correction_%d", i+4);
1635 corr = new AliCorrection(name, name);
1636 corr->LoadHistograms();
1637
1638 TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
1639 TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
1640
1641 // limit axes
1642 gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1643 meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1644 gene->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]);
1645 meas->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]);
1646
1647 TH1* decayed = gene->Project3D(Form("z_%d", i+4));
1648 TH1* stopped = meas->Project3D(Form("z_%d", i+4));
1649
1650 Printf("%d: %d decayed, %d stopped, out of %d", i, (Int_t) decayed->Integral(), (Int_t) stopped->Integral(), (Int_t) genePt->Integral());
1651
1652 decayed->Divide(decayed, genePt, 1, 1, "B");
1653 stopped->Divide(stopped, genePt, 1, 1, "B");
1654
1655 decayed->SetMarkerStyle(20);
1656 stopped->SetMarkerStyle(21);
1657 stopped->SetMarkerColor(2);
1658
1659 new TCanvas(Form("all_%d_%d", loop, i), Form("all_%d_%d", loop, i), 600, 600);
1660 effPt->DrawCopy();
1661 decayed->DrawCopy("SAME");
1662 stopped->DrawCopy("SAME");
1663
1664 decayed->Add(stopped);
1665 decayed->Add(effPt);
1666 decayed->SetMarkerStyle(22);
1667 decayed->SetMarkerColor(4);
1668 decayed->DrawCopy("SAME");
1669 }
1670
f3eb27f6 1671 }
1672
1673 Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1674
1675 Printf("%f measured, %f generated, effiency: %f", sumGen, sumMeas, sumMeas / sumGen);
1676
69b09e3b 1677 canvas->cd(loop+1);
1678 legend->Draw();
1679
1680 canvas2->cd();
f3eb27f6 1681 legend->Draw();
69b09e3b 1682 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
f3eb27f6 1683 }
1684
69b09e3b 1685 canvas->cd();
f3eb27f6 1686 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
69b09e3b 1687
1688 canvas3->cd();
1689 legends[0]->Draw();
1690 legends[1]->Draw();
1691 canvas3->SaveAs(Form("%s.eps", canvas3->GetName()));
f3eb27f6 1692}
1693
69b09e3b 1694void DrawpTCutOff()
f3eb27f6 1695{
69b09e3b 1696/*
1697aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt_ref.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1698mv unfolded.root chi2_ptref.root
1699aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1700mv unfolded.root chi2_pt0.root
1701aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1702mv unfolded.root chi2_pt1.root
1703aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1704mv unfolded.root chi2_pt0_25.root
1705aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1706mv unfolded.root chi2_pt1_25.root
1707*/
f3eb27f6 1708
69b09e3b 1709 loadlibs();
1710
f3eb27f6 1711 TH1* results[10];
69b09e3b 1712 const char* files[] = { "chi2_ptref.root", "chi2_pt0.root", "chi2_pt1.root", "chi2_pt0_25.root", "chi2_pt1_25.root"};
1713
1714 Int_t nMax = 5;
f3eb27f6 1715 for (Int_t i = 0; i<nMax; ++i)
1716 {
69b09e3b 1717 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(files[i]);
1718 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1719 }
1720
1721 const char* legendStrings[] = { "Reduced 50%", "Enhanced 50%", "Reduced 25%", "Enhanced 25%" };
1722 DrawRatio(results[0], nMax-1, results+1, "LowMomentumSyst.eps", kFALSE, legendStrings);
1723}
f3eb27f6 1724
69b09e3b 1725void ParticleSpeciesComparison()
1726{
1727 loadlibs();
f3eb27f6 1728
69b09e3b 1729 TH1* results[10];
1730 TH1* mc = 0;
1731
1732 // loop over cases (normal, enhanced/reduced ratios)
1733 Int_t nMax = 9;
1734 for (Int_t i = 0; i<nMax; ++i)
1735 {
1736 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(Form("chi2_species_%d.root", i), Form("Multiplicity_%d", i));
1737 if (i == 0)
1738 mc = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymchist", 1, 1);
f3eb27f6 1739 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
f3eb27f6 1740 }
1741
69b09e3b 1742 DrawResultRatio(mc, results[0], "ParticleSpeciesComparison1_1.eps");
f3eb27f6 1743
1744 for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1745 {
1746 results[0]->SetBinError(i, 0);
1747 mc->SetBinError(i, 0);
1748 }
1749
69b09e3b 1750 const char* legendStrings[] = { "K #times 0.5", "K #times 1.5", "p #times 0.5", "p #times 1.5", "K #times 0.5, p #times 0.5", "K #times 1.5, p #times 1.5", "K #times 0.5, p #times 1.5", "K #times 1.5, p #times 0.5" };
f3eb27f6 1751
1752 DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
1753
1754 //not valid: draw chi2 uncertainty on top!
1755 /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1756 TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1757 errorHist->SetLineColor(1);
1758 errorHist->SetLineWidth(2);
1759 TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1760 for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1761 {
1762 errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1763 errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1764 }
1765
1766 errorHist->DrawCopy("SAME");
1767 errorHist2->DrawCopy("SAME");*/
1768
1769 //canvas->SaveAs(canvas->GetName());
1770
1771 DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
1772
1773 //errorHist->DrawCopy("SAME");
1774 //errorHist2->DrawCopy("SAME");
1775
1776 //canvas2->SaveAs(canvas2->GetName());
1777}
1778
1779/*void ParticleSpeciesComparison2()
1780{
1781 gSystem->Load("libPWG0base");
1782
1783 const char* fileNameMC = "multiplicityMC_400k_syst.root";
1784 const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1785 Bool_t chi2 = 0;
1786
1787 TFile::Open(fileNameMC);
1788 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1789 mult->LoadHistograms();
1790
1791 TH1* mc[10];
1792 TH1* results[10];
1793
1794 // loop over cases (normal, enhanced/reduced ratios)
1795 Int_t nMax = 7;
1796 for (Int_t i = 0; i<nMax; ++i)
1797 {
1798 TString folder;
1799 folder.Form("Multiplicity_%d", i);
1800
1801 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1802
1803 TFile::Open(fileNameESD);
1804 mult2->LoadHistograms();
1805
1806 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1807
1808 if (chi2)
1809 {
1810 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1811 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1812 //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1813 }
1814 else
1815 {
1816 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1817 //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1818 }
1819
1820 //Float_t averageRatio = 0;
1821 //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1822
1823 results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1824
1825 TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1826 mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1827
1828 //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1829 //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1830
1831 //Printf("Case %d. Average ratio is %f", i, averageRatio);
1832 }
1833
1834 DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1835}*/
1836
1837TH1* Invert(TH1* eff)
1838{
1839 // calculate corr = 1 / eff
1840
1841 TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1842 corr->Reset();
1843
1844 for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1845 {
1846 if (eff->GetBinContent(i) > 0)
1847 {
1848 corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1849 corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1850 }
1851 }
1852
1853 return corr;
1854}
1855
1856void TriggerVertexCorrection()
1857{
1858 //
1859 // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1860 //
1861
69b09e3b 1862 loadlibs();
f3eb27f6 1863
1864 TFile::Open(correctionFile);
1865 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1866 mult->LoadHistograms("Multiplicity");
1867
1868 TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
69b09e3b 1869 TH1* corrNSD = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kNSD));
f3eb27f6 1870 TH1* corrMB = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1871
69b09e3b 1872 TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 500);
1873 gPad->SetGridx();
1874 gPad->SetGridy();
1875 gPad->SetTopMargin(0.05);
1876 gPad->SetRightMargin(0.05);
f3eb27f6 1877
1878 corrINEL->SetStats(kFALSE);
69b09e3b 1879 corrINEL->GetXaxis()->SetRangeUser(0, 12);
1880 corrINEL->GetYaxis()->SetRangeUser(0, 8);
1881 corrINEL->SetTitle(Form(";%s;Correction factor", GetMultLabel()));
f3eb27f6 1882 corrINEL->SetMarkerStyle(22);
1883 corrINEL->Draw("PE");
1884
1885 corrMB->SetStats(kFALSE);
1886 corrMB->SetLineColor(2);
1887 corrMB->SetMarkerStyle(25);
1888 corrMB->SetMarkerColor(2);
1889 corrMB->Draw("SAME PE");
69b09e3b 1890
1891 corrNSD->SetLineColor(4);
1892 corrNSD->SetMarkerStyle(24);
1893 corrNSD->SetMarkerColor(4);
1894 corrNSD->Draw("SAME PE");
1895
1896 Printf(" MB INEL NSD");
1897 Printf("bin 0: %f %f %f", corrMB->GetBinContent(1), corrINEL->GetBinContent(1), corrNSD->GetBinContent(1));
1898 Printf("bin 1: %f %f %f", corrMB->GetBinContent(2), corrINEL->GetBinContent(2), corrNSD->GetBinContent(2));
f3eb27f6 1899
69b09e3b 1900 TLegend* legend = new TLegend(0.3, 0.6, 0.85, 0.85);
f3eb27f6 1901 legend->SetFillColor(0);
69b09e3b 1902 legend->AddEntry(corrINEL, "Correction to inelastic sample");
1903 legend->AddEntry(corrNSD, "Correction to NSD sample");
1904 legend->AddEntry(corrMB, "Correction to triggered sample");
1905 legend->SetTextSize(0.04);
f3eb27f6 1906
1907 legend->Draw();
1908
1909 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1910}
1911
1912void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
1913{
69b09e3b 1914 loadlibs();
f3eb27f6 1915
1916 TFile::Open(correctionFile);
1917 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1918 mult->LoadHistograms("Multiplicity");
1919
1920 TFile::Open(measuredFile);
1921 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1922 mult2->LoadHistograms("Multiplicity");
1923
1924 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1925
69b09e3b 1926 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
f3eb27f6 1927
69b09e3b 1928 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 1e5);
1929
1930 TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1931
1932 return;
1933
f3eb27f6 1934 TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1935
f3eb27f6 1936 TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
1937
1938 if (!mc)
1939 {
1940 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1941 DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
1942 }
1943
69b09e3b 1944 TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
1945 errorResponse->Write();
1946 errorMeasured->Write();
1947 errorBoth->Write();
1948 file->Close();
1949}
1950
1951void DrawStatisticalUncertainty()
1952{
1953 TFile::Open("StatisticalUncertainty.root");
1954
1955 errorResponse = (TH1*) gFile->Get("errorResponse");
1956 errorMeasured = (TH1*) gFile->Get("errorMeasured");
1957 errorBoth = (TH1*) gFile->Get("errorBoth");
1958
f3eb27f6 1959 TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
1960 canvas->SetGridx();
1961 canvas->SetGridy();
1962 canvas->SetRightMargin(0.05);
1963 canvas->SetTopMargin(0.05);
1964
1965 errorResponse->SetLineColor(1);
69b09e3b 1966 errorResponse->GetXaxis()->SetRangeUser(0, longDisplayRange);
f3eb27f6 1967 errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1968 errorResponse->SetStats(kFALSE);
1969 errorResponse->SetTitle(";true multiplicity;Uncertainty");
1970
1971 errorResponse->Draw();
1972
1973 errorMeasured->SetLineColor(2);
1974 errorMeasured->Draw("SAME");
1975
1976 errorBoth->SetLineColor(4);
1977 errorBoth->Draw("SAME");
1978
69b09e3b 1979 Printf("Average errorResponse: %f", errorResponse->Integral(2, displayRange) / (displayRange - 1));
1980 Printf("Average errorMeasured: %f", errorMeasured->Integral(2, displayRange) / (displayRange - 1));
1981 Printf("Average errorBoth: %f", errorBoth->Integral(2, displayRange) / (displayRange - 1));
f3eb27f6 1982
1983 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
f3eb27f6 1984}
1985
1986void StatisticalUncertaintyCompare(const char* det = "SPD")
1987{
1988 TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1989 TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1990 TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1991 TH1* errorBoth = (TH1*) file1->Get("errorBoth");
1992
1993 TString str;
1994 str.Form("StatisticalUncertaintyCompare%s", det);
1995
69b09e3b 1996 TCanvas* canvas = new TCanvas(str, str, 800, 500);
f3eb27f6 1997 canvas->SetGridx();
1998 canvas->SetGridy();
1999 canvas->SetRightMargin(0.05);
2000 canvas->SetTopMargin(0.05);
69b09e3b 2001
2002 errorResponse->Scale(1.0 / sqrt(2));
2003 errorMeasured->Scale(1.0 / sqrt(2));
2004 errorBoth->Scale(1.0 / sqrt(2));
f3eb27f6 2005
2006 errorResponse->SetLineColor(1);
69b09e3b 2007 errorResponse->GetXaxis()->SetRangeUser(0, displayRange);
2008 errorResponse->GetYaxis()->SetRangeUser(0, 0.18);
f3eb27f6 2009 errorResponse->SetStats(kFALSE);
2010 errorResponse->GetYaxis()->SetTitleOffset(1.2);
69b09e3b 2011 errorResponse->SetTitle(Form(";%s;#sqrt{2}^{-1} #sigma(unfolded - unfolded_{0}) / unfolded_{0}", GetMultLabel()));
f3eb27f6 2012
2013 errorResponse->Draw();
2014
2015 errorMeasured->SetLineColor(2);
2016 errorMeasured->Draw("SAME");
2017
2018 errorBoth->SetLineColor(4);
2019 errorBoth->Draw("SAME");
2020
2021 TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
69b09e3b 2022 TH1* errorResponse2 = (TH1*) file2->Get("errorResponse");
2023 TH1* errorMeasured2 = (TH1*) file2->Get("errorMeasured");
f3eb27f6 2024 TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
2025
69b09e3b 2026 errorResponse2->Scale(1.0 / sqrt(2));
2027 errorMeasured2->Scale(1.0 / sqrt(2));
2028 errorBoth2->Scale(1.0 / sqrt(2));
2029
2030 errorResponse2->SetLineStyle(2);
2031 errorResponse2->Draw("SAME");
2032
2033 errorMeasured2->SetLineColor(2);
2034 errorMeasured2->SetLineStyle(2);
2035 errorMeasured2->Draw("SAME");
2036
f3eb27f6 2037 errorBoth2->SetLineColor(4);
2038 errorBoth2->SetLineStyle(2);
2039 errorBoth2->Draw("SAME");
2040
69b09e3b 2041 TLegend* legend = new TLegend(0.2, 0.5, 0.8, 0.9);
f3eb27f6 2042 legend->SetFillColor(0);
69b09e3b 2043 legend->SetTextSize(0.04);
2044 legend->AddEntry(errorBoth, "Both (Bayesian unfolding)");
2045 legend->AddEntry(errorMeasured, "Measured (Bayesian unfolding)");
2046 legend->AddEntry(errorResponse, "Response matrix (Bayesian unfolding)");
2047 legend->AddEntry(errorBoth2, "Both (#chi^{2}-minimization)");
2048 legend->AddEntry(errorMeasured2, "Measured (#chi^{2}-minimization)");
2049 legend->AddEntry(errorResponse2, "Response matrix (#chi^{2}-minimization)");
f3eb27f6 2050 legend->Draw();
2051
2052 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2053}
2054
2055void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
2056{
69b09e3b 2057 const char* files[] = { "multiplicityMC_nd.root", "multiplicityMC_sd.root", "multiplicityMC_dd.root", "multiplicityMC_xsection.root" };
f3eb27f6 2058
69b09e3b 2059 loadlibs();
f3eb27f6 2060
2061 TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
2062 canvas->SetGridx();
2063 canvas->SetGridy();
2064 canvas->SetRightMargin(0.05);
2065 canvas->SetTopMargin(0.05);
2066
2067 AliMultiplicityCorrection* data[4];
2068 TH1* effArray[4];
69b09e3b 2069 TH1* effErrorArray[2];
f3eb27f6 2070
2071 Int_t markers[] = { 24, 25, 26, 5 };
69b09e3b 2072 //Int_t markers[] = { 2, 25, 24, 5 };
2073 Int_t colors[] = { 1, 2, 4, 6 };
2074 //Int_t colors[] = { 1, 1, 1, 1 };
f3eb27f6 2075
69b09e3b 2076 //TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
2077 TLegend* legend = new TLegend(0.3, 0.3, 0.9, 0.6);
2078 legend->SetTextSize(0.04);
f3eb27f6 2079 legend->SetFillColor(0);
69b09e3b 2080
f3eb27f6 2081 for (Int_t i=0; i<4; ++i)
2082 {
2083 TString name;
2084 name.Form("Multiplicity_%d", i);
2085
2086 TFile::Open(files[i]);
2087 data[i] = new AliMultiplicityCorrection(name, name);
2088
2089 if (i < 3)
2090 {
2091 data[i]->LoadHistograms("Multiplicity");
2092 }
2093 else
2094 data[i]->LoadHistograms("Multiplicity_0");
2095
69b09e3b 2096 TH1* eff = 0;
2097 if (eventType == -1)
2098 {
2099 eff = (TH1*) data[i]->GetTriggerEfficiency(etaRange)->Clone(Form("eff_%d", i));
2100 }
2101 else
2102 eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
f3eb27f6 2103 effArray[i] = eff;
2104
2105 eff->GetXaxis()->SetRangeUser(0, 15);
69b09e3b 2106 eff->GetYaxis()->SetRangeUser(0, 1.19);
f3eb27f6 2107 eff->SetStats(kFALSE);
69b09e3b 2108 eff->GetXaxis()->SetTitle(GetMultLabel());
2109 eff->GetYaxis()->SetTitle("Efficiency");
2110 eff->SetTitle("");
f3eb27f6 2111 eff->SetLineColor(colors[i]);
2112 eff->SetMarkerColor(colors[i]);
2113 eff->SetMarkerStyle(markers[i]);
2114
2115 if (i == 3)
2116 {
69b09e3b 2117 // once for INEL, once for NSD
2118 for (AliMultiplicityCorrection::EventType eventType2 = AliMultiplicityCorrection::kINEL; eventType2 <= AliMultiplicityCorrection::kNSD; eventType2++)
f3eb27f6 2119 {
69b09e3b 2120 effDiff = (TH1*) data[i]->GetEfficiency(etaRange, eventType2)->Clone(Form("effDiff_%d", i));
2121
2122 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2123 effDiff->SetBinError(bin, 0);
2124
2125 // loop over cross section combinations
2126 for (Int_t j=1; j<7; ++j)
f3eb27f6 2127 {
69b09e3b 2128 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
2129 mult->LoadHistograms(Form("Multiplicity_%d", j));
2130
2131 TH1* eff2 = mult->GetEfficiency(etaRange, eventType2);
2132
2133 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2134 {
2135 // TODO we could also do asymmetric errors here
2136 Float_t deviation = TMath::Abs(effDiff->GetBinContent(bin) - eff2->GetBinContent(bin));
2137
2138 effDiff->SetBinError(bin, TMath::Max(effDiff->GetBinError(bin), (Double_t) deviation));
2139 }
2140 }
2141
2142 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2143 {
2144 //if (eventType2 == AliMultiplicityCorrection::kINEL)
2145 //eff->SetBinError(bin, 0);
2146 //eff->SetBinError(bin, effDiff->GetBinError(bin));
2147 if (bin < 20 && effDiff->GetBinContent(bin) > 0)
2148 Printf("Bin %d: Error: %.2f", bin, 100.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin));
2149 }
2150
2151 if (uncertainty) {
2152 TH1* effError = (TH1*) effDiff->Clone(Form("effError_%s", (eventType2 == AliMultiplicityCorrection::kINEL) ? "INEL" : "NSD"));
2153 effErrorArray[eventType2 - AliMultiplicityCorrection::kINEL] = effError;
2154 effError->Reset();
2155
2156 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2157 if (effDiff->GetBinContent(bin) > 0)
2158 effError->SetBinContent(bin, 1.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin));
2159
2160 effError->SetLineColor(1);
2161 effError->SetMarkerStyle(1);
2162
2163 if (eventType2 == AliMultiplicityCorrection::kNSD)
2164 effError->SetLineStyle(2);
2165
2166 effError->DrawCopy("SAME HIST");
f3eb27f6 2167 }
f3eb27f6 2168 }
2169 }
2170
f3eb27f6 2171 canvas->cd();
2172 if (i == 0)
2173 {
2174 eff->DrawCopy("P");
2175 }
2176 else
2177 eff->DrawCopy("SAME P");
2178
69b09e3b 2179 legend->AddEntry(eff, (((i == 0) ? "Non-diffractive" : ((i == 1) ? "Single-diffractive" : ((i == 2) ? "Double-diffractive" : "Pythia combined")))));
f3eb27f6 2180 }
2181
2182 if (uncertainty)
69b09e3b 2183 {
2184 legend->AddEntry(effErrorArray[0], "Relative syst. uncertainty: inelastic");
2185 legend->AddEntry(effErrorArray[1], "Relative syst. uncertainty: NSD");
2186
2187 file = TFile::Open("uncertainty_xsection.root", "RECREATE");
2188 effErrorArray[0]->Write();
2189 effErrorArray[1]->Write();
2190 file->Close();
2191 }
f3eb27f6 2192
2193 legend->Draw();
2194
2195 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2196}
2197
2198void ModelDependencyPlot()
2199{
69b09e3b 2200 loadlibs();
f3eb27f6 2201
69b09e3b 2202 TFile::Open("multiplicityMC.root");
f3eb27f6 2203 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2204 mult->LoadHistograms("Multiplicity");
69b09e3b 2205
2206 hist = mult->GetCorrelation(etaRange);
2207
2208 for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
2209 {
2210 for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
2211 {
2212 hist->SetBinContent(0, y, z, 0);
2213 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2214 }
2215 }
f3eb27f6 2216
69b09e3b 2217 TH2* proj = (TH2*) hist->Project3D("zy");
f3eb27f6 2218
69b09e3b 2219 TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 1200, 600);
f3eb27f6 2220
2221 canvas->Divide(2, 1);
2222
2223 canvas->cd(2);
2224 gPad->SetLogy();
69b09e3b 2225 gPad->SetGridx();
2226 gPad->SetGridy();
2227 gPad->SetTopMargin(0.05);
2228 gPad->SetRightMargin(0.05);
f3eb27f6 2229
2230 Int_t selectedMult = 30;
69b09e3b 2231 Int_t yMax = 9e4;
f3eb27f6 2232
2233 TH1* full = proj->ProjectionX("full");
2234 TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
2235
2236 full->SetStats(kFALSE);
69b09e3b 2237 full->GetXaxis()->SetRangeUser(0, displayRange);
f3eb27f6 2238 full->GetYaxis()->SetRangeUser(5, yMax);
69b09e3b 2239 full->SetTitle(";Multiplicity;Entries");
f3eb27f6 2240
2241 selected->SetLineColor(0);
2242 selected->SetMarkerColor(2);
69b09e3b 2243 selected->SetMarkerStyle(5);
f3eb27f6 2244
2245 full->Draw();
2246 selected->Draw("SAME P");
2247
2248 TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
2249 legend->SetFillColor(0);
69b09e3b 2250 legend->AddEntry(full, "True");
2251 legend->AddEntry(selected, "Measured");
f3eb27f6 2252 legend->Draw();
2253
2254 TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
2255 line->SetLineWidth(2);
2256 line->Draw();
2257
2258 canvas->cd(1);
2259 gPad->SetLogy();
69b09e3b 2260 gPad->SetGridx();
2261 gPad->SetGridy();
2262 gPad->SetTopMargin(0.05);
2263 gPad->SetRightMargin(0.05);
f3eb27f6 2264
2265 full = proj->ProjectionY("full2");
2266 selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
2267
2268 full->SetStats(kFALSE);
69b09e3b 2269 full->GetXaxis()->SetRangeUser(0, displayRange);
f3eb27f6 2270 full->GetYaxis()->SetRangeUser(5, yMax);
69b09e3b 2271 full->SetTitle(";Multiplicity;Entries");
f3eb27f6 2272
2273 full->SetLineColor(0);
2274 full->SetMarkerColor(2);
69b09e3b 2275 full->SetMarkerStyle(5);
f3eb27f6 2276
2277 full->Draw("P");
2278 selected->Draw("SAME");
2279
2280 legend->Draw();
2281
2282 line = new TLine(selectedMult, 5, selectedMult, yMax);
2283 line->SetLineWidth(2);
2284 line->Draw();
2285
2286 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2287}
2288
2289void SystematicpTSpectrum()
2290{
2291 gSystem->Load("libPWG0base");
2292
2293 TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
2294 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2295 mult->LoadHistograms("Multiplicity");
2296
2297 TFile::Open("multiplicityMC_100k_syst.root");
2298 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
69b09e3b 2299 mult2->LoadHistograms("Multiplicity");
f3eb27f6 2300
69b09e3b 2301 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2302 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2303 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
f3eb27f6 2304
69b09e3b 2305 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2306 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
f3eb27f6 2307
69b09e3b 2308 DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
2309}
f3eb27f6 2310
69b09e3b 2311void FitPt(const char* fileName)
2312{
2313 // needs a MC file from the dNdEta analysis
f3eb27f6 2314
69b09e3b 2315 TFile::Open(fileName);
f3eb27f6 2316
69b09e3b 2317 TH1* genePt = (TH1*) gFile->Get("fdNdpT");
2318
2319 genePt->SetTitle(";p_{T} (GeV/c);1/p_{T} dN_{ch}/dp_{T} (GeV/c)^{-2}");
2320 // number of events
2321 genePt->Scale(1.0 / 287800);
2322 // bin width
2323 genePt->Scale(1.0 / genePt->GetXaxis()->GetBinWidth(1));
2324
2325 genePt->GetXaxis()->SetRangeUser(0, 0.4);
f3eb27f6 2326
69b09e3b 2327 TF1* func = new TF1("func", "[1]*x*exp(x*[0])");
2328 func->SetParameters(-1, 1);
f3eb27f6 2329
2330 genePt->SetMarkerStyle(25);
2331 genePt->SetTitle("");
2332 genePt->SetStats(kFALSE);
f3eb27f6 2333
2334 new TCanvas;
2335 genePt->DrawCopy("P");
f3eb27f6 2336 func->DrawCopy("SAME");
2337 gPad->SetLogy();
2338
69b09e3b 2339 genePt->Fit(func, "0", "", 0, 0.25);
2340 genePt->Fit(func, "0", "", 0, 0.25);
f3eb27f6 2341
69b09e3b 2342 TCanvas* canvas = new TCanvas("FitPt", "FitPt", 600, 600);
f3eb27f6 2343
2344 gPad->SetGridx();
2345 gPad->SetGridy();
2346 gPad->SetLeftMargin(0.13);
2347 gPad->SetRightMargin(0.05);
2348 gPad->SetTopMargin(0.05);
2349
69b09e3b 2350 //genePt->GetXaxis()->SetRangeUser(0, 1);
2351 genePt->GetYaxis()->SetRangeUser(2, 200);
f3eb27f6 2352 genePt->GetYaxis()->SetTitleOffset(1.4);
2353 genePt->GetXaxis()->SetTitleOffset(1.1);
2354 genePt->DrawCopy("P");
69b09e3b 2355 //func->SetRange(0, 0.3);
f3eb27f6 2356 func->DrawCopy("SAME");
2357 gPad->SetLogy();
2358
2359 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
69b09e3b 2360
f3eb27f6 2361 TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2362
2363 TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2364
2365 TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2366
69b09e3b 2367 for (Int_t param=0; param<2; param++)
f3eb27f6 2368 {
2369 for (Int_t sign=0; sign<2; sign++)
2370 {
2371 TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign)); // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2372 func2->SetParameters(func->GetParameters());
2373 //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2374
69b09e3b 2375 Float_t factor = ((sign == 0) ? 0.75 : 1.25);
f3eb27f6 2376 func2->SetParameter(param, func2->GetParameter(param) * factor);
2377 //func2->Print();
2378
69b09e3b 2379 canvas->cd();
2380 func2->SetLineWidth(2);
f3eb27f6 2381 func2->SetLineColor(2);
69b09e3b 2382 func2->SetLineStyle(2);
f3eb27f6 2383 func2->DrawCopy("SAME");
2384
2385 canvas2->cd();
2386 TH1* second = func2->GetHistogram();
2387 second->Divide(first);
2388 second->SetLineColor(param + 1);
69b09e3b 2389 // set to 1 above 0.2 GeV
2390 for (Int_t bin=second->GetXaxis()->FindBin(0.20001); bin<=second->GetNbinsX(); bin++)
2391 second->SetBinContent(bin, 1);
f3eb27f6 2392 second->GetYaxis()->SetRangeUser(0, 2);
2393 second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2394 second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2395 }
2396 }
2397
2398 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2399 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2400
2401 file->Close();
2402}
2403
2404void DrawSystematicpT()
2405{
2406 TFile* file = TFile::Open("SystematicpT.root");
2407
2408 TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2409 TH1* result2 = (TH1*) file->Get("result_unity");
2410
2411 TH1* mcHist[12];
2412 TH1* result[12];
2413
2414 Int_t nParams = 5;
2415
2416 for (Int_t id=0; id<nParams*2; ++id)
2417 {
2418 mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2419 result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2420 }
2421
2422 DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2423
2424 //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2425
2426 DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2427
2428 //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2429
2430 // does not make sense: mc is different
2431 //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2432}
2433
2434void SystematicpT(Bool_t chi2 = 1)
2435{
2436 gSystem->Load("libPWG0base");
2437
2438 TFile::Open("ptspectrum900.root");
2439 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2440 mult->LoadHistograms("Multiplicity");
2441
2442 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2443
2444 TH1* mcHist[12];
2445 TH1* result[12];
2446
2447 Int_t nParams = 5;
2448
2449 for (Int_t param=0; param<nParams; param++)
2450 {
2451 for (Int_t sign=0; sign<2; sign++)
2452 {
2453 // calculate result with systematic effect
2454 TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
2455 mult2->LoadHistograms("Multiplicity");
2456
2457 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2458
2459 if (chi2)
2460 {
2461 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2462 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2463 }
2464 else
2465 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2466
2467 Int_t id = param * 2 + sign;
2468
2469 mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2470 result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2471
2472 TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2473 DrawResultRatio(mcHist[id], result[id], tmp);
2474 }
2475 }
2476
2477 // calculate normal result
2478 TFile::Open("ptspectrum100_1.root");
2479 mult2->LoadHistograms("Multiplicity");
2480 TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
2481
2482 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2483
2484 if (chi2)
2485 {
2486 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2487 }
2488 else
2489 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2490
2491 TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
2492
2493 TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2494 mcHist2->Write();
2495 result2->Write();
2496 for (Int_t id=0; id<nParams*2; ++id)
2497 {
2498 mcHist[id]->Write();
2499 result[id]->Write();
2500 }
2501 file->Close();
2502
2503 DrawSystematicpT();
2504}
2505
2506void DrawSystematicpT2()
2507{
2508 //displayRange = 200;
2509
2510 // read from file
2511 TFile* file = TFile::Open("SystematicpT2.root");
2512 TH1* mcHist = (TH1*) file->Get("mymc");
2513 TH1* result[12];
2514 result[0] = (TH1*) file->Get("result_unity");
2515 Int_t nParams = 5;
2516 for (Int_t id=0; id<nParams*2; ++id)
2517 result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
2518
2519 DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
2520 DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
2521 DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
2522}
2523
2524void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
2525{
2526 gSystem->Load("libPWG0base");
2527
2528 if (tpc)
2529 {
2530 SetTPC();
2531 TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2532 }
2533 else
2534 TFile::Open("ptspectrum100_1.root");
2535
2536 AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2537 measured->LoadHistograms("Multiplicity");
2538 TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2539
2540 TH1* result[12];
2541
2542 Int_t nParams = 5;
2543
2544 // -1 = unity change, 0...4 parameters
2545 for (Int_t id=-1; id<nParams*2; id++)
2546 {
2547 Int_t param = id / 2;
2548 Int_t sign = id % 2;
2549
2550 TString idStr;
2551 if (id == -1)
2552 {
2553 idStr = "unity";
2554 }
2555 else
2556 idStr.Form("%d_%d", param, sign);
2557
2558 // calculate result with systematic effect
2559 if (tpc)
2560 {
2561 TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
2562 }
2563 else
2564 TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
2565
2566 AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2567 response->LoadHistograms("Multiplicity");
2568
2569 response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
2570
2571 if (chi2)
2572 {
2573 response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2574 response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2575 }
2576 else
2577 response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2578
2579 result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
2580
2581 TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
2582 DrawResultRatio(mcHist, result[id+1], tmp);
2583 }
2584
2585 TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
2586 mcHist->Write();
2587 for (Int_t id=0; id<nParams*2+1; ++id)
2588 result[id]->Write();
2589 file->Close();
2590
2591 DrawSystematicpT2();
2592}
2593
2594void SystematicpTCutOff(Bool_t chi2 = kTRUE)
2595{
2596 // only needed for TPC
2597 SetTPC();
2598
2599 gSystem->Load("libPWG0base");
2600
2601 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2602 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2603 mult->LoadHistograms("Multiplicity");
2604
2605 TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2606 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2607 mult2->LoadHistograms("Multiplicity");
2608
2609 // "normal" result
2610 mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2611
2612 if (chi2)
2613 {
2614 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2615 mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2616 }
2617 else
2618 mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2619
2620 TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2621 TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
2622
2623 TH1* syst[2];
2624
2625 // change of pt spectrum (down)
2626 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2627 AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2628 mult3->LoadHistograms("Multiplicity");
2629 mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2630 if (chi2)
2631 {
2632 mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2633 mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2634 }
2635 else
2636 mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2637 syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2638
2639 // change of pt spectrum (up)
2640 TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
2641 AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
2642 mult4->LoadHistograms("Multiplicity");
2643 mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2644 if (chi2)
2645 {
2646 mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2647 mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2648 }
2649 else
2650 mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2651 syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
2652
2653 DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
2654
2655 Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
2656 Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
2657}
2658
69b09e3b 2659TH1* SystematicsSummary(Bool_t tpc = 0, Bool_t nsd = kTRUE)
f3eb27f6 2660{
2661 Int_t nEffects = 7;
2662
2663 TH1* effects[10];
2664 const char** names = 0;
69b09e3b 2665 Int_t colors[] = { 1, 2, 4, 1, 2, 4 };
2666 Int_t styles[] = { 1, 2, 3, 1, 2, 3 };
2667 Int_t widths[] = { 1, 1, 1, 2, 2, 2 };
f3eb27f6 2668 Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
69b09e3b 2669
2670 TH1* dummy = new TH2F("dummy", Form(";%s;Uncertainty", GetMultLabel()), 202, -1.5, 200.5, 100, 0, 0.4);
2671 dummy->SetStats(0);
f3eb27f6 2672
2673 for (Int_t i=0; i<nEffects; ++i)
69b09e3b 2674 effects[i] = new TH1F("SystematicsSummary", Form(";%s;Uncertainty", GetMultLabel()), 201, -0.5, 200.5);
f3eb27f6 2675
2676 if (tpc)
2677 {
2678 SetTPC();
2679
2680 const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
2681 names = namesTPC;
2682
2683 // method
2684 TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
2685 TH1* hist = (TH1*) file->Get("errorBoth");
2686
2687 // smooth a bit, but skip 0 bin
2688 effects[0]->SetBinContent(2, hist->GetBinContent(2));
2689 for (Int_t i=3; i<=200; ++i)
2690 effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2691
2692 // relative x-section
2693 effects[1]->SetBinContent(2, 0.005);
2694 effects[1]->SetBinContent(3, 0.0025);
2695 effects[1]->SetBinContent(4, 0.0025);
2696
2697 // particle composition
2698 for (Int_t i=2; i<=101; ++i)
2699 {
2700 if (i < 41)
2701 {
2702 effects[2]->SetBinContent(i, 0.01);
2703 }
2704 else if (i < 76)
2705 {
2706 effects[2]->SetBinContent(i, 0.02);
2707 }
2708 else
2709 effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
2710 }
2711
2712 // pt cut off (only tpc)
2713 for (Int_t i=2; i<=101; ++i)
2714 {
2715 if (i < 11)
2716 {
2717 effects[3]->SetBinContent(i, 0.05);
2718 }
2719 else if (i < 51)
2720 {
2721 effects[3]->SetBinContent(i, 0.01);
2722 }
2723 else
2724 effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
2725 }
2726
2727 // track selection (only tpc)
2728 for (Int_t i=2; i<=101; ++i)
2729 effects[4]->SetBinContent(i, 0.03);
2730
2731 // secondaries
2732 for (Int_t i=2; i<=101; ++i)
2733 effects[5]->SetBinContent(i, 0.01);
2734
2735 // pt spectrum
2736 for (Int_t i=2; i<=101; ++i)
2737 {
2738 if (i < 21)
2739 {
2740 effects[6]->SetBinContent(i, 0.05);
2741 }
2742 else if (i < 51)
2743 {
2744 effects[6]->SetBinContent(i, 0.02);
2745 }
2746 else
2747 effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
2748 }
2749
2750 }
2751 else
2752 {
f3eb27f6 2753 nEffects = 5;
2754
69b09e3b 2755 //const char* namesSPD[] = { "Particle composition", "p_{t} cut-off", "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)"};
2756 const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)", "Particle composition", "p_{t} cut-off"};
f3eb27f6 2757 names = namesSPD;
2758
69b09e3b 2759 currentEffect = 0;
2760
f3eb27f6 2761 // method
2762 TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
2763 TH1* hist = (TH1*) file->Get("errorBoth");
69b09e3b 2764 hist->Scale(1.0 / sqrt(2));
f3eb27f6 2765
2766 // smooth a bit, but skip 0 bin
69b09e3b 2767 /*effects[currentEffect]->SetBinContent(1, hist->GetBinContent(1));
2768 for (Int_t i=2; i<=201; ++i)
2769 effects[currentEffect]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);*/
2770 effects[currentEffect] = hist;
f3eb27f6 2771
69b09e3b 2772 currentEffect++;
f3eb27f6 2773
69b09e3b 2774 // relative x-section
2775 file = TFile::Open("uncertainty_xsection.root");
2776 effects[currentEffect++] = (TH1*) file->Get("effError_INEL");
2777 effects[currentEffect] = (TH1*) file->Get("effError_NSD");
2778 effects[currentEffect]->SetLineStyle(1);
2779 //effects[2]->SetBinContent(1, 0.20);
2780 //effects[2]->SetBinContent(2, 0.01);
2781 //effects[2]->SetBinContent(3, 0.002);
2782
2783 currentEffect++;
2784
f3eb27f6 2785 // particle composition
69b09e3b 2786 effects[currentEffect]->SetBinContent(1, 0.16);
2787 for (Int_t i=2; i<=81; ++i)
f3eb27f6 2788 {
69b09e3b 2789 effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * i / 81);
f3eb27f6 2790 }
69b09e3b 2791
2792 currentEffect++;
f3eb27f6 2793
2794 // pt spectrum
69b09e3b 2795 effects[currentEffect]->SetBinContent(1, 0.06);
2796 effects[currentEffect]->SetBinContent(2, 0.03);
2797 for (Int_t i=3; i<=81; ++i)
f3eb27f6 2798 {
69b09e3b 2799 if (i <= 61)
f3eb27f6 2800 {
69b09e3b 2801 effects[currentEffect]->SetBinContent(i, 0.01);
f3eb27f6 2802 }
69b09e3b 2803 else if (i <= 81)
f3eb27f6 2804 {
69b09e3b 2805 effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * (i - 61) / 20);
f3eb27f6 2806 }
f3eb27f6 2807 }
69b09e3b 2808
2809// currentEffect++;
2810//
2811// // material budget
2812// for (Int_t i=1; i<=81; ++i)
2813// {
2814// if (i < 5)
2815// effects[currentEffect]->SetBinContent(i, 0.05 - 0.01 * i);
2816// if (i > 51)
2817// effects[currentEffect]->SetBinContent(i, 0.05 * (i - 50) / 30);
2818// }
2819//
2820 currentEffect++;
2821
f3eb27f6 2822 }
2823
69b09e3b 2824 TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 500);
2825 canvas->SetRightMargin(0.05);
f3eb27f6 2826 canvas->SetTopMargin(0.05);
69b09e3b 2827 //canvas->SetGridx();
2828 canvas->SetGridy();
2829 TLegend* legend = new TLegend(0.2, 0.4, 0.7, 0.4 + 0.5 * nEffects / 7);
f3eb27f6 2830 legend->SetFillColor(0);
69b09e3b 2831 legend->SetTextSize(0.04);
2832 dummy->Draw();
2833 dummy->GetXaxis()->SetRangeUser(0, displayRange);
f3eb27f6 2834
2835 for (Int_t i=0; i<nEffects; ++i)
2836 {
2837 TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
2838 /*current->Reset();
2839 for (Int_t j=0; j<nEffects-i; ++j)
2840 current->Add(effects[j]);*/
2841
2842 current->SetLineColor(colors[i]);
69b09e3b 2843 current->SetLineStyle(styles[i]);
2844 current->SetLineWidth(widths[i]);
f3eb27f6 2845 //current->SetFillColor(colors[i]);
2846 current->SetMarkerColor(colors[i]);
2847 //current->SetMarkerStyle(markers[i]);
2848
2849 current->SetStats(kFALSE);
2850 current->GetYaxis()->SetRangeUser(0, 0.4);
69b09e3b 2851 current->DrawCopy("SAME");
f3eb27f6 2852 legend->AddEntry(current, names[i]);
2853
69b09e3b 2854 //TLatex* text = new TLatex(displayRange+2, current->GetBinContent(displayRange+1), names[i]);
2855 TLatex* text = new TLatex(displayRange+2, 0.1 - i * 0.02, names[i]);
2856 text->SetTextSize(0.04);
f3eb27f6 2857 text->SetTextColor(colors[i]);
69b09e3b 2858 //text->Draw();
f3eb27f6 2859 }
2860
2861 // add total in square
69b09e3b 2862 TH1* totalINEL = (TH1*) effects[0]->Clone("totalINEL");
2863 totalINEL->Reset();
2864 TH1* totalNSD = (TH1*) totalINEL->Clone("totalNSD");
f3eb27f6 2865
2866 for (Int_t i=0; i<nEffects; ++i)
2867 {
2868 //Printf("%d %f", i, effects[i]->GetBinContent(20));
2869 effects[i]->Multiply(effects[i]);
69b09e3b 2870
2871 if (i != 2)
2872 totalINEL->Add(effects[i]);
2873 if (i != 1)
2874 totalNSD->Add(effects[i]);
2875 }
2876
2877 for (Int_t i=1; i<=totalINEL->GetNbinsX(); ++i)
2878 {
2879 totalINEL->SetBinError(i, 0);
2880 if (totalINEL->GetBinContent(i) > 0)
2881 totalINEL->SetBinContent(i, TMath::Min(sqrt(totalINEL->GetBinContent(i)), 1.0));
2882 totalNSD->SetBinError(i, 0);
2883 if (totalNSD->GetBinContent(i) > 0)
2884 totalNSD->SetBinContent(i, TMath::Min(sqrt(totalNSD->GetBinContent(i)), 1.0));
f3eb27f6 2885 }
f3eb27f6 2886
2887 //Printf("%f", total->GetBinContent(20));
2888
69b09e3b 2889 totalINEL->SetMarkerStyle(5);
2890 totalINEL->SetMarkerColor(1);
2891 legend->AddEntry(totalINEL, "Total (INEL)", "P");
2892
2893 totalNSD->SetMarkerStyle(24);
2894 totalNSD->SetMarkerColor(2);
2895 legend->AddEntry(totalNSD, "Total (NSD)", "P");
2896
2897 Printf("total in bin 0 is INEL: %f NSD: %f", totalINEL->GetBinContent(1), totalNSD->GetBinContent(1));
2898 totalINEL->DrawCopy("SAME P"); //->SetBinContent(1, 0);
2899 totalNSD->DrawCopy("SAME P"); //->SetBinContent(1, 0);
f3eb27f6 2900
2901 legend->Draw();
2902
2903 canvas->SaveAs(canvas->GetName());
2904
69b09e3b 2905 return (nsd) ? totalNSD : totalINEL;
f3eb27f6 2906}
2907
69b09e3b 2908void finalPlot(Bool_t tpc = 0, Bool_t small = kFALSE)
f3eb27f6 2909{
69b09e3b 2910 loadlibs();
f3eb27f6 2911
2912 if (tpc)
2913 SetTPC();
2914
69b09e3b 2915 //TH1* errorNSD = SystematicsSummary(tpc, 1);
f3eb27f6 2916
69b09e3b 2917 TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", (small) ? 600 : 800, 500);
f3eb27f6 2918 canvas->SetRightMargin(0.05);
2919 canvas->SetTopMargin(0.05);
69b09e3b 2920 canvas->SetGridx();
2921 canvas->SetGridy();
2922
2923 legend = new TLegend(0.5, 0.6, 0.9, 0.8);
2924 legend->SetFillColor(0);
2925 legend->SetTextSize(0.04);
2926
2927 for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++)
2928 {
2929 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root");
2930 TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc");
2931 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2932
2933 DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType));
f3eb27f6 2934
69b09e3b 2935 // normalize result
2936 //result->Scale(1.0 / result->Integral(2, displayRange));
2937
2938 result->GetXaxis()->SetRangeUser(0, displayRange);
2939 //result->SetBinContent(1, 0); result->SetBinError(1, 0);
2940 result->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (etaRange+1) * 0.5));
2941 result->SetMarkerStyle(0);
2942 result->SetLineColor(1);
2943 result->SetStats(kFALSE);
2944
2945 // systematic error
2946 TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD));
2947
2948 TH1* systError = (TH1*) result->Clone("systError");
2949 for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
2950 systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
2951
2952 // change error drawing style
2953 systError->SetFillColor(15);
2954
2955 if (eventType == AliMultiplicityCorrection::kNSD)
2956 {
2957 result->SetLineColor(2);
2958 result->SetMarkerColor(2);
2959 result->SetMarkerStyle(5);
2960 }
2961
2962 canvas->cd();
2963 systError->DrawCopy(Form("E2 ][ %s", (eventType == AliMultiplicityCorrection::kINEL) ? "" : "SAME"));
2964 result->DrawCopy("SAME E ][");
2965 canvas->SetLogy();
2966
2967 legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic cross-section" : "NSD cross-section", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P");
2968 }
2969
2970 legend->Draw();
2971 /*
f3eb27f6 2972 //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
2973 TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
2974 text->SetFillColor(0);
2975 text->SetTextAlign(12);
2976 text->AddText("Systematic errors summed quadratically");
2977 text->AddText("0.6 million minimum bias events");
2978 text->AddText("corrected to inelastic events");
2979 text->Draw("B");
2980
2981 TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
2982 text2->SetFillColor(0);
2983 text2->SetTextAlign(12);
2984 text2->AddText("#sqrt{s} = 14 TeV");
2985 if (tpc)
2986 {
2987 text2->AddText("|#eta| < 0.9");
2988 }
2989 else
2990 text2->AddText("|#eta| < 2.0");
2991 text2->AddText("simulated data (PYTHIA)");
2992 text2->Draw("B");
69b09e3b 2993
f3eb27f6 2994 if (tpc)
2995 {
2996 TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
2997 text3->SetNDC();
2998 text3->Draw();
2999 }
3000 else
3001 {
3002 TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
3003 text3->SetNDC();
3004 text3->Draw();
3005 }
3006
3007 // alice logo
3008 TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
3009 pad->Draw();
3010 pad->cd();
3011 TImage* img = TImage::Open("$HOME/alice.png");
3012 img->SetImageQuality(TAttImage::kImgBest);
3013 img->Draw();
3014
3015 canvas->Modified();
69b09e3b 3016 */
f3eb27f6 3017
3018/* TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
3019 text->SetTextSize(0.04);
3020 text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
3021 text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
3022 text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
3023 text->Draw();*/
3024
3025
3026 canvas->SaveAs(canvas->GetName());
3027}
3028
69b09e3b 3029TMatrixD* NonInvertable()
3030{
3031 const Int_t kSize = 5;
3032
3033 TMatrixD matrix(kSize, kSize);
3034 for (Int_t x=0; x<kSize; x++)
3035 {
3036 for (Int_t y=0; y<kSize; y++)
3037 {
3038 if (x == y)
3039 {
3040 if (x == 0 || x == kSize -1)
3041 {
3042 matrix(x, y) = 0.75;
3043 }
3044 else
3045 matrix(x, y) = 0.5;
3046 }
3047 else if (TMath::Abs(x - y) == 1)
3048 {
3049 matrix(x, y) = 0.25;
3050 }
3051 }
3052 }
3053
3054 matrix.Print();
3055
3056 //TMatrixD inverted(matrix);
3057 //inverted.Invert();
3058
3059 //inverted.Print();
3060
3061 return new TMatrixD(matrix);
3062}
3063
f3eb27f6 3064void BlobelUnfoldingExample()
3065{
3066 const Int_t kSize = 20;
3067
3068 TMatrixD matrix(kSize, kSize);
3069 for (Int_t x=0; x<kSize; x++)
3070 {
3071 for (Int_t y=0; y<kSize; y++)
3072 {
3073 if (x == y)
3074 {
3075 if (x == 0 || x == kSize -1)
3076 {
3077 matrix(x, y) = 0.75;
3078 }
3079 else
3080 matrix(x, y) = 0.5;
3081 }
3082 else if (TMath::Abs(x - y) == 1)
3083 {
3084 matrix(x, y) = 0.25;
3085 }
3086 }
3087 }
3088
3089 //matrix.Print();
3090
3091 TMatrixD inverted(matrix);
3092 inverted.Invert();
3093
3094 //inverted.Print();
3095
69b09e3b 3096 TH1F* inputDist = new TH1F("inputDist", ";t;Entries", kSize, -0.5, (Float_t) kSize - 0.5);
f3eb27f6 3097 TVectorD inputDistVector(kSize);
3098 TH1F* unfolded = inputDist->Clone("unfolded");
3099 TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
69b09e3b 3100 measuredIdealDist->SetTitle(";m;Entries");
f3eb27f6 3101 TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
3102
3103 TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
3104 // norm: 1/(sqrt(2pi)sigma)
3105 gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
3106 //gaus->Print();
3107
3108 for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
3109 {
3110 Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
3111 inputDist->SetBinContent(x, value);
3112 inputDistVector(x-1) = value;
3113 }
3114
3115 TVectorD measuredDistIdealVector = matrix * inputDistVector;
3116
3117 for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
3118 measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
3119
3120 measuredDist->FillRandom(measuredIdealDist, 10000);
3121
3122 // fill error matrix before scaling
3123 TMatrixD covarianceMatrixMeasured(kSize, kSize);
3124 for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
3125 covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
3126
3127 TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
3128 //covarianceMatrix.Print();
3129
3130 TVectorD measuredDistVector(kSize);
3131 for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
3132 measuredDistVector(x-1) = measuredDist->GetBinContent(x);
3133
3134 TVectorD unfoldedVector = inverted * measuredDistVector;
3135 for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
3136 unfolded->SetBinContent(x, unfoldedVector(x-1));
3137
69b09e3b 3138 TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1200, 600);
f3eb27f6 3139 canvas->SetTopMargin(0.05);
3140 canvas->Divide(2, 1);
3141
3142 canvas->cd(1);
3143 canvas->cd(1)->SetLeftMargin(0.15);
3144 canvas->cd(1)->SetRightMargin(0.05);
69b09e3b 3145 canvas->cd(1)->SetTopMargin(0.05);
3146 gPad->SetGridx();
3147 gPad->SetGridy();
3148 measuredDist->GetYaxis()->SetRangeUser(-600, 2799);
3149 measuredDist->GetYaxis()->SetTitleOffset(1.9);
f3eb27f6 3150 measuredDist->SetStats(0);
3151 measuredDist->DrawCopy();
3152 gaus->Draw("SAME");
3153
3154 canvas->cd(2);
3155 canvas->cd(2)->SetLeftMargin(0.15);
3156 canvas->cd(2)->SetRightMargin(0.05);
69b09e3b 3157 canvas->cd(2)->SetTopMargin(0.05);
3158 gPad->SetGridx();
3159 gPad->SetGridy();
3160 unfolded->GetYaxis()->SetRangeUser(-600, 2799);
3161 unfolded->GetYaxis()->SetTitleOffset(1.9);
f3eb27f6 3162 unfolded->SetStats(0);
3163 unfolded->DrawCopy();
3164 gaus->Draw("SAME");
3165
3166 canvas->SaveAs("BlobelUnfoldingExample.eps");
69b09e3b 3167
3168 return;
3169
3170 // now unfold this with Bayesian
3171 loadlibs();
3172
3173 // fill a multiplicity object
3174 mult = new AliMultiplicityCorrection("mult", "mult");
3175 for (Int_t x=0; x<kSize; x++)
3176 {
3177 mult->GetMultiplicityVtx(0)->SetBinContent(1, x+1, inputDistVector(x));
3178 mult->GetMultiplicityESD(0)->SetBinContent(1, x+1, measuredDistVector(x)*10000);
3179 for (Int_t y=0; y<kSize; y++)
3180 mult->GetCorrelation(0)->SetBinContent(1, x+1, y+1, matrix(x, y));
3181 }
3182
3183 //mult->DrawHistograms();
3184
3185 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 0);
3186 //mult->SetCreateBigBin(kFALSE);
3187 mult->ApplyMinuitFit(0, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); //hist2->ProjectionY("mymchist"));
3188
3189 //mult->ApplyBayesianMethod(0, kFALSE, AliMultiplicityCorrection::kTrVtx, 0, -1, 0, kFALSE);
3190
3191 mult->DrawComparison("BlobelExample", 0, kFALSE, kTRUE, mult->GetMultiplicityVtx(0)->ProjectionY("mcmchist", 1, mult->GetMultiplicityVtx(0)->GetNbinsX()));
3192
3193
f3eb27f6 3194}
3195
3196void E735Fit()
3197{
3198 TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
3199 fCurrentESD->Sumw2();
3200
3201 // Open the input stream
3202 ifstream in;
3203 in.open("e735data.txt");
3204
3205 while(in.good())
3206 {
3207 Float_t x, y, ye;
3208 in >> x >> y >> ye;
3209
3210 //Printf("%f %f %f", x, y, ye);
3211 fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
3212 fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
3213 }
3214
3215 in.close();
3216
3217 //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
3218
3219 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
3220
3221 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])");
3222 func->SetParNames("scaling", "averagen", "k");
3223 func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
3224 func->SetParLimits(1, 0.001, 1000);
3225 func->SetParLimits(2, 0.001, 1000);
3226 func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
3227
3228 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
3229 lognormal->SetParNames("scaling", "mean", "sigma");
3230 lognormal->SetParameters(1, 1, 1);
3231 lognormal->SetParLimits(0, 0, 10);
3232 lognormal->SetParLimits(1, 0, 100);
3233 lognormal->SetParLimits(2, 1e-3, 10);
3234
3235 TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
3236 fCurrentESD->SetStats(kFALSE);
69b09e3b 3237 fCurrentESD->SetMarkerStyle(0);
3238 fCurrentESD->SetLineColor(1);
f3eb27f6 3239 fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
3240 fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
3241 fCurrentESD->Draw("");
3242 fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
3243 fCurrentESD->Fit(func, "0", "", 0, 150);
3244 func->SetRange(0, 250);
3245 func->Draw("SAME");
3246 printf("chi2 = %f\n", func->GetChisquare());
3247
3248 fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
3249 lognormal->SetLineColor(2);
3250 lognormal->SetLineStyle(2);
3251 lognormal->SetRange(0, 250);
3252 lognormal->Draw("SAME");
3253
3254 gPad->SetLogy();
3255
3256 canvas->SaveAs("E735Fit.eps");
3257}
69b09e3b 3258
3259void DifferentSamples()
3260{
3261 loadlibs();
3262
3263 Int_t n = 2;
3264 const char* filesChi2[] = { "chi2_100k_1.root", "chi2_100k_2.root" };
3265 const char* filesBayesian[] = { "bayesian_100k_1.root", "bayesian_100k_2.root" };
3266
3267 TCanvas* canvas = new TCanvas("DifferentSamples", "DifferentSamples", 1200, 600);
3268 canvas->Divide(2, 1);
3269
3270 legend = new TLegend(0.15, 0.7, 0.65, 0.9);
3271 legend->SetFillColor(0);
3272 legend->SetTextSize(0.04);
3273
3274 for (Int_t i=0; i<n; i++)
3275 {
3276 AliMultiplicityCorrection* chi2 = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3277 TFile::Open(filesChi2[i]);
3278 chi2->LoadHistograms("Multiplicity");
3279
3280 AliMultiplicityCorrection* bayesian = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3281 TFile::Open(filesBayesian[i]);
3282 bayesian->LoadHistograms("Multiplicity");
3283
3284 chi2Hist = chi2->GetMultiplicityESDCorrected(etaRange);
3285 bayesianHist = bayesian->GetMultiplicityESDCorrected(etaRange);
3286
3287 mc = chi2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
3288
3289 // normalize and divide
3290 chi2Hist->Scale(1.0 / chi2Hist->Integral(2, displayRange+1) * mc->Integral(2, displayRange));
3291 bayesianHist->Scale(1.0 / bayesianHist->Integral(2, displayRange+1) * mc->Integral(2, displayRange));
3292
3293 chi2Hist->Divide(mc, chi2Hist);
3294 bayesianHist->Divide(mc, bayesianHist);
3295
3296 canvas->cd(i+1);
3297 gPad->SetTopMargin(0.05);
3298 gPad->SetRightMargin(0.05);
3299 //gPad->SetLeftMargin(0.12);
3300 gPad->SetGridx();
3301 gPad->SetGridy();
3302
3303 chi2Hist->GetXaxis()->SetRangeUser(0, displayRange);
3304 chi2Hist->GetYaxis()->SetTitleOffset(1.3);
3305 chi2Hist->SetStats(0);
3306 chi2Hist->SetTitle(Form(";%s;MC / unfolded", GetMultLabel()));
3307 chi2Hist->GetYaxis()->SetRangeUser(0.2, 1.8);
3308 chi2Hist->Draw("HIST");
3309
3310 for (Int_t x=1; x<=bayesianHist->GetNbinsX(); x++)
3311 bayesianHist->SetBinError(x, 1e-6);
3312
3313 bayesianHist->SetLineColor(2);
3314 bayesianHist->SetMarkerColor(2);
3315 bayesianHist->SetMarkerStyle(5);
3316 bayesianHist->Draw("HIST E SAME");
3317
3318 if (i == 0)
3319 {
3320 legend->AddEntry(chi2Hist, "#chi^{2}-minimization", "L");
3321 legend->AddEntry(bayesianHist, "Bayesian unfolding", "LP");
3322 }
3323 legend->Draw();
3324 }
3325
3326 canvas->SaveAs("DifferentSamples.eps");
3327}
3328
3329void PileUp()
3330{
3331 loadlibs();
3332
3333 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root");
3334 hist2d = mult->GetMultiplicityMC(etaRange, AliMultiplicityCorrection::kINEL);
3335 mult1 = hist2d->ProjectionY("mult1", 1, hist2d->GetNbinsX());
3336
3337 conv = (TH1*) mult1->Clone("conv");
3338 conv->Reset();
3339
3340 mult1->Scale(1.0 / mult1->Integral());
3341
3342 for (Int_t i=1; i<=mult1->GetNbinsX(); i++)
3343 for (Int_t j=1; j<=mult1->GetNbinsX(); j++)
3344 conv->Fill(mult1->GetBinCenter(i)+mult1->GetBinCenter(j), mult1->GetBinContent(i) * mult1->GetBinContent(j));
3345
3346 conv->Scale(1.0 / conv->Integral());
3347
3348 c = new TCanvas("c", "c", 800, 500);
3349 gPad->SetLogy();
3350 gPad->SetTopMargin(0.05);
3351 gPad->SetRightMargin(0.05);
3352 mult1->SetTitle(Form(";%s;Probability", GetMultLabel()));
3353 mult1->SetStats(0);
3354 gPad->SetGridx();
3355 gPad->SetGridy();
3356 mult1->Draw();
3357 mult1->GetYaxis()->SetRangeUser(1e-7, 2 * mult1->GetMaximum());
3358 mult1->GetXaxis()->SetRangeUser(0, displayRange);
3359 mult1->GetXaxis()->SetTitleOffset(1.15);
3360 conv->SetLineColor(2);
3361 conv->SetMarkerColor(2);
3362 conv->SetMarkerStyle(5);
3363 conv->DrawCopy("SAME P");
3364
3365 conv->Scale(0.00058);
3366 conv->DrawCopy("SAME P");
3367
3368 legend = new TLegend(0.73, 0.73, 0.93, 0.93);
3369 legend->SetFillColor(0);
3370 legend->SetTextSize(0.04);
3371 legend->AddEntry(mult1, "1 collision");
3372 legend->AddEntry(conv, "2 collisions", "P");
3373 legend->Draw();
3374
3375 c->SaveAs("pileup.eps");
3376
3377 new TCanvas;
3378 conv->Divide(mult1);
3379 conv->Draw();
3380}
3381
3382void TestErrorDetermination(Int_t nRandomizations)
3383{
3384 TF1* func = new TF1("nbd", "exp(log([0]) + TMath::LnGamma([2]+x) - TMath::LnGamma([2]) - TMath::LnGamma(x+1) + log([1] / ([1]+[2])) * x + log(1.0 + [1]/[2]) * -[2])", 0, 100);
3385 func->SetParNames("scaling", "averagen", "k");
3386 func->SetParameters(1, 15, 2);
3387
3388 TF1* func2 = new TF1("nbd2", "exp(log([0]) + TMath::LnGamma([2]+x) - TMath::LnGamma([2]) - TMath::LnGamma(x+1) + log([1] / ([1]+[2])) * x + log(1.0 + [1]/[2]) * -[2])", 0, 100);
3389 func2->SetParNames("scaling", "averagen", "k");
3390 func2->SetParLimits(0, 0.5, 2);
3391 func2->SetParLimits(1, 1, 50);
3392 func2->SetParLimits(2, 1, 10);
3393 func2->SetParameters(1, 15, 2);
3394 //func2->FixParameter(0, 1);
3395
3396 //new TCanvas; func->Draw("L");
3397
3398 hist1 = new TH1F("hist1", "", 100, 0.5, 100.5);
3399 hist2 = new TH1F("hist2", "", 100, 0.5, 100.5);
3400 hist1->Sumw2();
3401
3402 TH1* params[3];
3403 params[0] = new TH1F("param_0", Form("param_%d", 0), 100, 0.95, 1.05);
3404 params[1] = new TH1F("param_1", Form("param_%d", 1), 100, 14, 16);
3405 params[2] = new TH1F("param_2", Form("param_%d", 2), 100, 1.8, 2.2);
3406
3407 const Int_t nTries = 1000;
3408 for (Int_t i=0; i<nTries; i++)
3409 {
3410 hist1->Reset();
3411
3412 if (nRandomizations == 1)
3413 {
3414 hist1->FillRandom("nbd", 10000);
3415 }
3416 else if (nRandomizations == 2)
3417 {
3418 hist2->Reset();
3419 hist2->FillRandom("nbd", 10000);
3420 hist1->FillRandom(hist2, 10000);
3421 }
3422 else if (nRandomizations == 3)
3423 {
3424 hist2->Reset();
3425 hist1->FillRandom("nbd", 10000);
3426 hist2->FillRandom(hist1, 10000);
3427 hist1->Reset();
3428 hist1->FillRandom(hist2, 10000);
3429 }
3430 else
3431 return;
3432
3433 //new TCanvas; hist1->Draw();
3434
3435 hist1->Scale(1.0 / hist1->Integral());
3436 hist1->Fit(func2, "NQ");
3437 hist1->Fit(func2, "NQ");
3438 for (Int_t j=0; j<3; j++)
3439 params[j]->Fill(func2->GetParameter(j));
3440 }
3441
3442 for (Int_t j=0; j<3; j++)
3443 {
3444 new TCanvas; params[j]->Draw();
3445 params[j]->Fit("gaus");
3446 Printf("sigma of param %d if %f", j, ((TF1*) params[j]->FindObject("gaus"))->GetParameter(2));
3447 }
3448}
3449
3450void DrawRawDistributions(const char* fileName = "multiplicityESD.root")
3451{
3452 loadlibs();
3453
3454 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileName);
3455
3456 c = new TCanvas("c", "c", 600, 600);
3457
3458 dummy = new TH2F("dummy", ";measured multiplicity", 100, -0.5, 149.5, 100, 0.5, 4e4);
3459 dummy->SetStats(0);
3460 dummy->Draw();
3461 gPad->SetGridx();
3462 gPad->SetGridy();
3463 gPad->SetLogy();
3464
3465 Int_t colors[] = { 1, 2, 4 };
3466
3467 for (Int_t i=2; i>=0; i--)
3468 {
3469 hist = mult->GetMultiplicityESD(i)->ProjectionY();
3470
3471 hist->SetLineColor(colors[i]);
3472 hist->DrawCopy("SAME");
3473 }
3474
3475
3476}