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