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