skeleton for documentation
[u/mrichter/AliRoot.git] / PWG0 / multiplicity / correct.C
CommitLineData
0ab29cfa 1/* $Id$ */
2
3//
dca331bb 4// script to correct the multiplicity spectrum + helpers
0ab29cfa 5//
6
0b4bfd98 7void SetTPC()
0a173978 8{
9 gSystem->Load("libPWG0base");
0b4bfd98 10 AliMultiplicityCorrection::SetQualityRegions(kFALSE);
11}
0a173978 12
2440928d 13void draw(const char* fileName = "multiplicity.root", const char* folder = "Multiplicity")
0b4bfd98 14{
2440928d 15 loadlibs();
0a173978 16
0b4bfd98 17 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
0a173978 18
0b4bfd98 19 TFile::Open(fileName);
20 mult->LoadHistograms();
0a173978 21 mult->DrawHistograms();
447c325d 22
f3eb27f6 23 TH2* hist = (TH2*) gROOT->FindObject("fCorrelation2_zy");
447c325d 24 canvas = new TCanvas("c1", "c1", 600, 500);
f3eb27f6 25 canvas->SetTopMargin(0.05);
447c325d 26 hist->SetStats(kFALSE);
27 hist->Draw("COLZ");
f3eb27f6 28 hist->SetTitle(";true multiplicity in |#eta| < 1.5;measured multiplicity in |#eta| < 1.5");
447c325d 29 hist->GetYaxis()->SetTitleOffset(1.1);
f3eb27f6 30 gPad->SetRightMargin(0.12);
447c325d 31 gPad->SetLogz();
32
f3eb27f6 33 canvas->SaveAs("responsematrix.eps");
0ab29cfa 34}
35
0f67a57c 36void loadlibs()
9ca6be09 37{
0f67a57c 38 gSystem->Load("libTree");
39 gSystem->Load("libVMC");
40
41 gSystem->Load("libSTEERBase");
42 gSystem->Load("libANALYSIS");
2440928d 43 gSystem->Load("libANALYSISalice");
9ca6be09 44 gSystem->Load("libPWG0base");
0f67a57c 45}
9ca6be09 46
0a0f4adb 47void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityESD.root", Bool_t chi2 = kTRUE, Int_t histID = 2, Bool_t fullPhaseSpace = kFALSE, Float_t beta = 1e3, Int_t eventType = 0 /* AliMultiplicityCorrection::kTrVtx */)
0f67a57c 48{
49 loadlibs();
0a0f4adb 50
0b4bfd98 51 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
9ca6be09 52 TFile::Open(fileNameMC);
0b4bfd98 53 mult->LoadHistograms();
9ca6be09 54
f3eb27f6 55 AliMultiplicityCorrection* esd = new AliMultiplicityCorrection(folder, folder);
9ca6be09 56 TFile::Open(fileNameESD);
f3eb27f6 57 esd->LoadHistograms();
58
59 TH2F* hist = esd->GetMultiplicityESD(histID);
60 TH2F* hist2 = esd->GetMultiplicityMC(histID, eventType);
9ca6be09 61
cfc19dd5 62 mult->SetMultiplicityESD(histID, hist);
9ca6be09 63
447c325d 64 // small hack to get around charge conservation for full phase space ;-)
65 if (fullPhaseSpace)
66 {
67 TH1* corr = mult->GetCorrelation(histID + 4);
dd701109 68
447c325d 69 for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
70 for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
71 {
72 corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
73 corr->SetBinError(i, j, corr->GetBinError(i-1, j));
74 }
75 }
dd701109 76
447c325d 77 if (chi2)
78 {
0f67a57c 79 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, beta);
80 //mult->SetCreateBigBin(kFALSE);
81 //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); //mult->SetCreateBigBin(kFALSE);
82 //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE);
83 //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e5);
0b4bfd98 84 //mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5);
44466df2 85 //mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kTRUE, hist2->ProjectionY("mymchist"));
0a0f4adb 86 mult->ApplyMinuitFit(histID, fullPhaseSpace, eventType, kFALSE); //hist2->ProjectionY("mymchist"));
447c325d 87 }
88 else
dd701109 89 {
0a0f4adb 90 mult->ApplyBayesianMethod(histID, fullPhaseSpace, eventType, 0.2, 100);
dd701109 91 }
92
dca331bb 93 TFile* file = TFile::Open("unfolded.root", "RECREATE");
144ff489 94 mult->SaveHistograms();
95 file->Write();
96 file->Close();
0a0f4adb 97
98 mult->DrawComparison((chi2) ? "MinuitChi2" : "Bayesian", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
99}
100
f3eb27f6 101void CompareChi2Bayesian(Int_t histID = 2, const char* chi2File = "chi2.root", const char* bayesianFile = "bayesian.root", const char* label1 = "Chi2", const char* label2 = "Bayesian", const char* mcFile = 0, Float_t simpleCorrect = 0)
0a0f4adb 102{
103 const char* folder = "Multiplicity";
104
105 loadlibs();
106
107 AliMultiplicityCorrection* chi2 = new AliMultiplicityCorrection(folder, folder);
108 TFile::Open(chi2File);
109 chi2->LoadHistograms();
110
111 AliMultiplicityCorrection* bayesian = new AliMultiplicityCorrection(folder, folder);
112 TFile::Open(bayesianFile);
113 bayesian->LoadHistograms();
114
f3eb27f6 115 histRAW = chi2->GetMultiplicityESD(histID)->ProjectionY("raw", 1, chi2->GetMultiplicityESD(histID)->GetNbinsX());
116 histRAW->Scale(1.0 / histRAW->Integral());
117
0a0f4adb 118 histC = chi2->GetMultiplicityESDCorrected(histID);
119 histB = bayesian->GetMultiplicityESDCorrected(histID);
120
121 c = new TCanvas("CompareChi2Bayesian", "CompareChi2Bayesian", 800, 600);
122 c->SetRightMargin(0.05);
123 c->SetTopMargin(0.05);
124 c->SetLogy();
f3eb27f6 125 c->SetGridx();
126 c->SetGridy();
0a0f4adb 127
128 histC->SetTitle(";N;P(N)");
129 histC->SetStats(kFALSE);
130 histC->GetXaxis()->SetRangeUser(0, 100);
131
132 histC->SetLineColor(1);
133 histB->SetLineColor(2);
f3eb27f6 134 histRAW->SetLineColor(3);
0a0f4adb 135
51f6de65 136 histC->DrawCopy("HISTE");
137 histB->DrawCopy("HISTE SAME");
f3eb27f6 138 histRAW->DrawCopy("SAME");
0a0f4adb 139
140 legend = new TLegend(0.2, 0.2, 0.4, 0.4);
141 legend->SetFillColor(0);
142
f3eb27f6 143 legend->AddEntry(histC, label1);
144 legend->AddEntry(histB, label2);
145 legend->AddEntry(histRAW, "raw ESD");
0a0f4adb 146
51f6de65 147 histGraph = 0;
f3eb27f6 148 if (simpleCorrect > 0)
149 {
150 graph = new TGraph;
151 graph->SetMarkerStyle(25);
152 graph->SetFillColor(0);
153 for (Int_t bin=1; bin<=histRAW->GetNbinsX(); bin++)
154 graph->SetPoint(graph->GetN(), histRAW->GetXaxis()->GetBinCenter(bin) * simpleCorrect, histRAW->GetBinContent(bin));
155
156 graph->Draw("PSAME");
157 legend->AddEntry(graph, "weighting");
158
159 // now create histogram from graph and normalize
160 histGraph = (TH1*) histRAW->Clone();
161 histGraph->Reset();
162 for (Int_t bin=1; bin<=histGraph->GetNbinsX(); bin++)
163 {
164 Int_t j=1;
165 for (j=1; j<graph->GetN(); j++)
166 if (graph->GetX()[j] > histGraph->GetXaxis()->GetBinCenter(bin))
167 break;
168 if (j == graph->GetN())
169 continue;
170 if (histGraph->GetXaxis()->GetBinCenter(bin) - graph->GetX()[j] < graph->GetX()[j-1] - histGraph->GetXaxis()->GetBinCenter(bin))
171 j--;
172 histGraph->SetBinContent(bin, graph->GetY()[j]);
173 }
174
175 Printf("Integral = %f", histGraph->Integral());
176 histGraph->Scale(1.0 / histGraph->Integral());
177
178 histGraph->SetLineColor(6);
179 histGraph->DrawCopy("SAME");
180 legend->AddEntry(histGraph, "weighting normalized");
181 }
182
183 if (mcFile)
0a0f4adb 184 {
185 AliMultiplicityCorrection* mc = new AliMultiplicityCorrection(folder, folder);
186 TFile::Open(mcFile);
187 mc->LoadHistograms();
188
189 histMC = mc->GetMultiplicityVtx(histID)->ProjectionY("mc", 1, mc->GetMultiplicityVtx(histID)->GetNbinsX());
f3eb27f6 190 histMC->Sumw2();
0a0f4adb 191 histMC->Scale(1.0 / histMC->Integral());
192
51f6de65 193 histMC->Draw("HISTE SAME");
0a0f4adb 194 histMC->SetLineColor(4);
195 legend->AddEntry(histMC, "MC");
196 }
197
198 legend->Draw();
199
f3eb27f6 200 c->SaveAs(Form("%s.png", c->GetName()));
201 c->SaveAs(Form("%s.eps", c->GetName()));
202
203 if (!mcFile)
204 return;
205
206 // build ratios
207
208 c = new TCanvas("CompareChi2BayesianRatio", "CompareChi2BayesianRatio", 800, 600);
209 c->SetRightMargin(0.05);
210 c->SetTopMargin(0.05);
211 c->SetGridx();
212 c->SetGridy();
213
214 for (Int_t bin=1; bin<=histC->GetNbinsX(); bin++)
215 {
216 if (histMC->GetBinContent(bin) > 0)
217 {
218 histC->SetBinContent(bin, histC->GetBinContent(bin) / histMC->GetBinContent(bin));
219 histB->SetBinContent(bin, histB->GetBinContent(bin) / histMC->GetBinContent(bin));
51f6de65 220
221 /*
222 if (histGraph)
223 {
224 histGraph->SetBinContent(bin, histGraph->GetBinContent(bin) / histMC->GetBinContent(bin));
225 histGraph->SetBinError(bin, 0);
226 }
227 */
f3eb27f6 228
229 // TODO errors?
230 histC->SetBinError(bin, 0);
231 histB->SetBinError(bin, 0);
f3eb27f6 232 }
233 }
234
235 histC->GetYaxis()->SetRangeUser(0.5, 2);
236 histC->GetYaxis()->SetTitle("Unfolded / MC");
237
238 histC->Draw("HIST");
239 histB->Draw("HIST SAME");
51f6de65 240 /*
241 if (histGraph)
242 histGraph->Draw("HIST SAME");
243 */
f3eb27f6 244
245 /*
246 if (simpleCorrect > 0)
247 {
248 graph2 = new TGraph;
249 graph2->SetMarkerStyle(25);
250 graph2->SetFillColor(0);
251 for (Int_t i=0; i<graph->GetN(); i++)
252 {
253 Float_t mcValue = histMC->GetBinContent(histMC->FindBin(graph->GetX()[i]));
254 Float_t mcError = histMC->GetBinError(histMC->FindBin(graph->GetX()[i]));
255 if (mcValue > 0)
256 graph2->SetPoint(graph2->GetN(), graph->GetX()[i], graph->GetY()[i] / mcValue);
257 }
258 graph2->Draw("PSAME");
259 }
260 */
261
262 c->SaveAs(Form("%s.png", c->GetName()));
263 c->SaveAs(Form("%s.eps", c->GetName()));
264}
265
266
267void CompareMC(Int_t histID, Int_t eventType, const char* file1, const char* file2, const char* label1, const char* label2)
268{
269 const char* folder = "Multiplicity";
270
271 loadlibs();
272
273 AliMultiplicityCorrection* mc1 = new AliMultiplicityCorrection(folder, folder);
274 TFile::Open(file1);
275 mc1->LoadHistograms();
276 histMC1 = mc1->GetMultiplicityMC(histID, eventType)->ProjectionY("mc1", 1, mc1->GetMultiplicityMC(histID, eventType)->GetNbinsX());
277 histMC1->Scale(1.0 / histMC1->Integral());
278
279 AliMultiplicityCorrection* mc2 = new AliMultiplicityCorrection(folder, folder);
280 TFile::Open(file2);
281 mc2->LoadHistograms();
282 histMC2 = mc2->GetMultiplicityMC(histID, eventType)->ProjectionY("mc2", 1, mc2->GetMultiplicityMC(histID, eventType)->GetNbinsX());
283 histMC2->Scale(1.0 / histMC2->Integral());
284
285 c = new TCanvas("CompareMC", "CompareMC", 800, 600);
286 c->SetRightMargin(0.05);
287 c->SetTopMargin(0.05);
288 c->SetLogy();
289
290 histMC1->SetTitle(";N;P(N)");
291 histMC1->SetStats(kFALSE);
292 histMC1->GetXaxis()->SetRangeUser(0, 100);
293
294 histMC1->SetLineColor(1);
295 histMC2->SetLineColor(2);
296
297 histMC1->Draw();
298 histMC2->Draw("SAME");
299
300 legend = new TLegend(0.2, 0.2, 0.4, 0.4);
301 legend->SetFillColor(0);
302
303 legend->AddEntry(histMC1, label1);
304 legend->AddEntry(histMC2, label2);
305
306 legend->Draw();
307
0a0f4adb 308 c->SaveAs(Form("%s.gif", c->GetName()));
309 c->SaveAs(Form("%s.eps", c->GetName()));
447c325d 310}
311
312void* fit2Step(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE)
313{
314 gSystem->Load("libPWG0base");
315
316 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
317
318 TFile::Open(fileNameMC);
319 mult->LoadHistograms("Multiplicity");
320
321 TFile::Open(fileNameESD);
322 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
323 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID)));
324 //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID));
325
326 mult->SetMultiplicityESD(histID, hist);
327
328 // small hack to get around charge conservation for full phase space ;-)
329 if (fullPhaseSpace)
330 {
331 TH1* corr = mult->GetCorrelation(histID + 4);
332
333 for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
334 for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
335 {
336 corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
337 corr->SetBinError(i, j, corr->GetBinError(i-1, j));
338 }
339 }
340
341 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
342 mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
343 mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
344
345 TH1* result = (TH1*) mult->GetMultiplicityESDCorrected((fullPhaseSpace) ? 4 : histID))->Clone("firstresult");
346
347 mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 100000);
348 mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE, result);
349 mult->DrawComparison("MinuitChi2_Step2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
0a173978 350
351 return mult;
352}
cfc19dd5 353
354const char* GetRegName(Int_t type)
355{
356 switch (type)
357 {
358 case AliMultiplicityCorrection::kNone: return "None"; break;
359 case AliMultiplicityCorrection::kPol0: return "Pol0"; break;
360 case AliMultiplicityCorrection::kPol1: return "Pol1"; break;
361 case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break;
362 case AliMultiplicityCorrection::kEntropy: return "Reduced cross-entropy"; break;
0b4bfd98 363 case AliMultiplicityCorrection::kLog : return "Log"; break;
dd701109 364 }
365 return 0;
366}
367
368const char* GetEventTypeName(Int_t type)
369{
370 switch (type)
371 {
372 case AliMultiplicityCorrection::kTrVtx: return "trigger, vertex"; break;
373 case AliMultiplicityCorrection::kMB: return "minimum bias"; break;
374 case AliMultiplicityCorrection::kINEL: return "inelastic"; break;
cfc19dd5 375 }
376 return 0;
377}
378
447c325d 379void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
380{
381 gSystem->mkdir(targetDir);
382
383 gSystem->Load("libPWG0base");
384
385 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
386 TFile::Open(fileNameMC);
387 mult->LoadHistograms("Multiplicity");
388
389 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
390 TFile::Open(fileNameESD);
391 multESD->LoadHistograms("Multiplicity");
392 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
393
394 Int_t count = 0; // just to order the saved images...
395
396 TFile* graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir), "RECREATE");
397
0b4bfd98 398 Int_t colors[3] = {1, 2, 4};
0f67a57c 399 Int_t markers[20] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3, 4, 5, 6};
0b4bfd98 400
401 for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
402 //for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
447c325d 403 {
404 TString tmp;
405 tmp.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
406
407 TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600);
408
0f67a57c 409 for (Int_t i = 1; i <= 5; i++)
447c325d 410 {
0f67a57c 411 Int_t iterArray[5] = {5, 20, 50, 100, -1};
0b4bfd98 412 //Int_t iter = i * 40 - 20;
413 Int_t iter = iterArray[i-1];
414
415 TGraph* fitResultsMC[3];
416 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
417 {
418 fitResultsMC[region] = new TGraph;
419 fitResultsMC[region]->SetTitle(Form("%d iter. - reg. %d", iter, region+1));
420 fitResultsMC[region]->GetXaxis()->SetTitle("smoothing parameter #alpha");
421 fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region));
422 fitResultsMC[region]->SetName(Form("%s_MC_%d", tmp.Data(), i * AliMultiplicityCorrection::kQualityRegions + region - 2));
423 fitResultsMC[region]->SetFillColor(0);
0f67a57c 424 //fitResultsMC[region]->SetMarkerStyle(markers[(i-1) * AliMultiplicityCorrection::kQualityRegions + region]);
425 fitResultsMC[region]->SetMarkerStyle(markers[(i-1)]);
0b4bfd98 426 fitResultsMC[region]->SetLineColor(colors[region]);
427 }
428
447c325d 429 TGraph* fitResultsRes = new TGraph;
430 fitResultsRes->SetTitle(Form("%d iterations", iter));
431 fitResultsRes->SetName(Form("%s_Res_%d", tmp.Data(), i));
432 fitResultsRes->GetXaxis()->SetTitle("smoothing parameter");
433 fitResultsRes->GetYaxis()->SetTitle("P_{2}");
434
447c325d 435 fitResultsRes->SetFillColor(0);
447c325d 436 fitResultsRes->SetMarkerStyle(19+i);
0b4bfd98 437 fitResultsRes->SetMarkerColor(1);
438 fitResultsRes->SetLineColor(1);
447c325d 439
440 for (Float_t weight = 0.0; weight < 1.01; weight += 0.2)
441 {
442 Float_t chi2MC = 0;
443 Float_t residuals = 0;
444
0b4bfd98 445 mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter, 0, kFALSE);
447c325d 446 mult->DrawComparison(Form("%s/BayesianIterSmooth_%03d_%d_%d_%f", targetDir, count++, type, iter, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
447 mult->GetComparisonResults(&chi2MC, 0, &residuals);
448
0b4bfd98 449 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
450 fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region));
451
447c325d 452 fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
453 }
454
0b4bfd98 455 graphFile->cd();
456 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
457 fitResultsMC[region]->Write();
458
447c325d 459 fitResultsRes->Write();
460 }
461 }
0b4bfd98 462
463 graphFile->Close();
447c325d 464}
465
466void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = kTRUE)
467{
468 gSystem->Load("libPWG0base");
469
470 TString name;
471 TFile* graphFile = 0;
472 if (type == -1)
473 {
474 name = "EvaluateChi2Method";
475 graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir));
476 }
477 else
478 {
479 name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
480 graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir));
481 }
482
483 TCanvas* canvas = new TCanvas(name, name, 800, 500);
484 if (type == -1)
0b4bfd98 485 {
447c325d 486 canvas->SetLogx();
0b4bfd98 487 canvas->SetLogy();
488 }
489 canvas->SetTopMargin(0.05);
490 canvas->SetGridx();
491 canvas->SetGridy();
447c325d 492
0b4bfd98 493 TLegend* legend = new TLegend(0.8, 0.15, 0.98, 0.98);
447c325d 494 legend->SetFillColor(0);
495
496 Int_t count = 1;
497
498 Float_t xMin = 1e20;
499 Float_t xMax = 0;
500
501 Float_t yMin = 1e20;
502 Float_t yMax = 0;
503
0b4bfd98 504 Float_t yMinRegion[3];
505 for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
506 yMinRegion[i] = 1e20;
507
447c325d 508 TString xaxis, yaxis;
509
510 while (1)
511 {
512 TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
513 TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
514
0b4bfd98 515 if (!mc)
447c325d 516 break;
517
518 xaxis = mc->GetXaxis()->GetTitle();
519 yaxis = mc->GetYaxis()->GetTitle();
520
521 mc->Print();
447c325d 522
0b4bfd98 523 if (res)
524 res->Print();
447c325d 525
0b4bfd98 526 xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
527 yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin());
447c325d 528
0b4bfd98 529 xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
530 yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
531
532 if (plotRes && res)
447c325d 533 {
0b4bfd98 534 xMin = TMath::Min(xMin, res->GetXaxis()->GetXmin());
535 yMin = TMath::Min(yMin, res->GetYaxis()->GetXmin());
447c325d 536
0b4bfd98 537 xMax = TMath::Max(xMax, res->GetXaxis()->GetXmax());
538 yMax = TMath::Max(yMax, res->GetYaxis()->GetXmax());
447c325d 539 }
0b4bfd98 540
541 for (Int_t i=0; i<mc->GetN(); ++i)
542 yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]);
543
544 count++;
447c325d 545 }
546
0b4bfd98 547 for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
548 Printf("Minimum for region %d is %f", i, yMinRegion[i]);
549
447c325d 550 if (type >= 0)
551 {
552 xaxis = "smoothing parameter";
553 }
554 else if (type == -1)
555 {
556 xaxis = "weight parameter";
0b4bfd98 557 xMax *= 5;
447c325d 558 }
0b4bfd98 559 //yaxis = "P_{1} (2 <= t <= 150)";
447c325d 560
561 printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
562
563 TGraph* dummy = new TGraph;
564 dummy->SetPoint(0, xMin, yMin);
565 dummy->SetPoint(1, xMax, yMax);
566 dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data()));
567
568 dummy->SetMarkerColor(0);
569 dummy->Draw("AP");
0b4bfd98 570 dummy->GetYaxis()->SetMoreLogLabels(1);
447c325d 571
572 count = 1;
573
574 while (1)
575 {
576 TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
577 TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
578
0b4bfd98 579 //printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res);
447c325d 580
0b4bfd98 581 if (!mc)
447c325d 582 break;
583
584 printf("Loaded %d sucessful.\n", count);
585
44466df2 586 if (type == -1)
447c325d 587 {
44466df2 588 legend->AddEntry(mc, Form("Eq. (%d) - reg. %d", 10 + (count-1) / 3, 1+ (count-1) % 3));
447c325d 589 }
44466df2 590 else
447c325d 591 legend->AddEntry(mc);
592
593 mc->Draw("SAME PC");
594
0b4bfd98 595 if (plotRes && res)
447c325d 596 {
597 legend->AddEntry(res);
598 res->Draw("SAME PC");
599 }
600
601 count++;
602 }
603
604 legend->Draw();
605
606 canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
607 canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
608}
609
0f67a57c 610void EvaluateDrawResultRegions(const char* targetDir, Int_t type = 0)
611{
612 gSystem->Load("libPWG0base");
613
614 TString name;
615 TFile* graphFile = 0;
616 if (type == -1)
617 {
618 name = "EvaluateChi2Method";
619 graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir));
620 }
621 else
622 {
623 name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
624 graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir));
625 }
626
627 TCanvas* canvas = new TCanvas(name, name, 800, 1200);
628 //canvas->Divide(1, AliMultiplicityCorrection::kQualityRegions, 0, 0);
629 canvas->Range(0, 0, 1, 1);
630
631 TPad* pad[3];
632 pad[0] = new TPad(Form("%s_pad1", name.Data()), "", 0.02, 0.05, 0.98, 0.35);
633 pad[1] = new TPad(Form("%s_pad2", name.Data()), "", 0.02, 0.35, 0.98, 0.65);
634 pad[2] = new TPad(Form("%s_pad3", name.Data()), "", 0.02, 0.65, 0.98, 0.95);
635
636 Float_t yMinRegion[3];
637 for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
638 yMinRegion[i] = 1e20;
639
640 for (Int_t region = 1; region <= AliMultiplicityCorrection::kQualityRegions; region++)
641 {
642 canvas->cd();
643 pad[region-1]->Draw();
644 pad[region-1]->cd();
645 pad[region-1]->SetRightMargin(0.05);
646
647 if (region != 1)
648 pad[region-1]->SetBottomMargin(0);
649 if (region != AliMultiplicityCorrection::kQualityRegions)
650 pad[region-1]->SetTopMargin(0);
651
652 if (type == -1)
653 {
654 pad[region-1]->SetLogx();
655 pad[region-1]->SetLogy();
656 }
657 pad[region-1]->SetTopMargin(0.05);
658 pad[region-1]->SetGridx();
659 pad[region-1]->SetGridy();
660
661 TLegend* legend = new TLegend(0.8, 0.15, 0.98, 0.98);
662 legend->SetFillColor(0);
663
664 Int_t count = 0;
665
666 Float_t xMin = 1e20;
667 Float_t xMax = 0;
668
669 Float_t yMin = 1e20;
670 Float_t yMax = 0;
671
672 TString xaxis, yaxis;
673
674 while (1)
675 {
676 count++;
677
678 TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
679 if (!mc)
680 break;
681
682 if (TString(mc->GetTitle()).Contains(Form("reg. %d", region)) == kFALSE)
683 continue;
684
685 xaxis = mc->GetXaxis()->GetTitle();
686 yaxis = mc->GetYaxis()->GetTitle();
687
688 mc->Print();
689
690 xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
691 yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin());
692
693 xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
694 yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
695
696 for (Int_t i=0; i<mc->GetN(); ++i)
697 yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]);
698 }
699
700 if (type >= 0)
701 {
702 xaxis = "smoothing parameter";
703 }
704 else if (type == -1)
705 {
706 xaxis = "weight parameter";
707 xMax *= 5;
708 }
709 yaxis = "P_{1}"; // (2 <= t <= 150)";
710
711 printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
712
713 TGraph* dummy = new TGraph;
714 dummy->SetPoint(0, xMin, yMin);
715 dummy->SetPoint(1, xMax, yMax);
716 dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data()));
717
718 dummy->SetMarkerColor(0);
719 dummy->Draw("AP");
720 dummy->GetYaxis()->SetMoreLogLabels(1);
721
722 count = 0;
723
724 while (1)
725 {
726 count++;
727
728 TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
729 if (!mc)
730 break;
731 if (TString(mc->GetTitle()).Contains(Form("reg. %d", region)) == kFALSE)
732 continue;
733
734 printf("Loaded %d sucessful.\n", count);
735
736 if (type == -1)
737 {
738 legend->AddEntry(mc, Form("Eq. (%d) - reg. %d", 10 + (count-1) / 3, 1+ (count-1) % 3));
739 }
740 else
741 legend->AddEntry(mc);
742
743 mc->Draw("SAME PC");
744
745 }
746
747 legend->Draw();
748 }
749
750 for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
751 Printf("Minimum for region %d is %f", i, yMinRegion[i]);
752
753 canvas->Modified();
754 canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
755 canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
756}
757
0b4bfd98 758void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", const char* targetDir, Int_t histID = 3)
dd701109 759{
760 gSystem->mkdir(targetDir);
761
cfc19dd5 762 gSystem->Load("libPWG0base");
763
764 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
765
766 TFile::Open(fileNameMC);
767 mult->LoadHistograms("Multiplicity");
768
769 TFile::Open(fileNameESD);
770 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
771 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
772
773 mult->SetMultiplicityESD(histID, hist);
774
dd701109 775 Int_t count = 0; // just to order the saved images...
0b4bfd98 776 Int_t colors[3] = {1, 2, 4};
777 Int_t markers[12] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3};
dd701109 778
447c325d 779 TGraph* fitResultsRes = 0;
780
781 TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE");
cfc19dd5 782
0f67a57c 783 for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type)
0b4bfd98 784// for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kEntropy; type <= AliMultiplicityCorrection::kEntropy; ++type)
0f67a57c 785 //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
cfc19dd5 786 {
0b4bfd98 787 TGraph* fitResultsMC[3];
788 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
789 {
790 fitResultsMC[region] = new TGraph;
44466df2 791 fitResultsMC[region]->SetTitle(Form("Eq. (%d) - reg. %d", type+9, region+1));
0b4bfd98 792 fitResultsMC[region]->GetXaxis()->SetTitle("weight parameter #alpha");
793 fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region));
794 fitResultsMC[region]->SetName(Form("EvaluateChi2Method_MC_%d", type * AliMultiplicityCorrection::kQualityRegions + region - 2));
795 fitResultsMC[region]->SetFillColor(0);
796 fitResultsMC[region]->SetMarkerStyle(markers[(type-1) * AliMultiplicityCorrection::kQualityRegions + region]);
797 fitResultsMC[region]->SetLineColor(colors[region]);
798 }
799
447c325d 800 fitResultsRes = new TGraph;
801 fitResultsRes->SetTitle(Form("%s residual chi2", GetRegName(type)));
802 fitResultsRes->SetName(Form("EvaluateChi2Method_Res_%d", type));
803 fitResultsRes->GetXaxis()->SetTitle("Weight Parameter");
cfc19dd5 804
cfc19dd5 805 fitResultsRes->SetFillColor(0);
dd701109 806 fitResultsRes->SetMarkerStyle(23+type);
cfc19dd5 807 fitResultsRes->SetMarkerColor(kRed);
808 fitResultsRes->SetLineColor(kRed);
809
0b4bfd98 810 for (Int_t i=0; i<7; ++i)
447c325d 811 {
812 Float_t weight = TMath::Power(TMath::Sqrt(10), i+6);
813 //Float_t weight = TMath::Power(10, i+2);
cfc19dd5 814
447c325d 815 //if (type == AliMultiplicityCorrection::kEntropy) weight = 1e4 * (i+1) * 1.5;
cfc19dd5 816
cfc19dd5 817 Float_t chi2MC = 0;
818 Float_t residuals = 0;
447c325d 819 Float_t chi2Limit = 0;
820
821 TString runName;
822 runName.Form("MinuitChi2_%02d_%d_%f", count++, type, weight);
cfc19dd5 823
824 mult->SetRegularizationParameters(type, weight);
447c325d 825 Int_t status = mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
826 mult->DrawComparison(Form("%s/%s", targetDir, runName.Data()), histID, kFALSE, kTRUE, hist2->ProjectionY());
827 if (status != 0)
828 {
829 printf("MINUIT did not succeed. Skipping...\n");
830 continue;
831 }
832
dd701109 833 mult->GetComparisonResults(&chi2MC, 0, &residuals);
447c325d 834 TH1* result = mult->GetMultiplicityESDCorrected(histID);
835 result->SetName(runName);
836 result->Write();
cfc19dd5 837
0b4bfd98 838 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
839 fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region));
840
cfc19dd5 841 fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
cfc19dd5 842 }
843
447c325d 844 graphFile->cd();
0b4bfd98 845 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
846 fitResultsMC[region]->Write();
447c325d 847 fitResultsRes->Write();
cfc19dd5 848 }
849
447c325d 850 graphFile->Close();
dd701109 851}
852
853void EvaluateChi2MethodAll()
854{
855 EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M");
856 EvaluateChi2Method("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M");
857 EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD");
858 EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M");
859 EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
860}
861
862void EvaluateBayesianMethodAll()
863{
864 EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M");
865 EvaluateBayesianMethod("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M");
866 EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD");
867 EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M");
868 EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
869}
870
447c325d 871void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
dd701109 872{
873 gSystem->mkdir(targetDir);
874
875 gSystem->Load("libPWG0base");
876
877 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
878
879 TFile::Open(fileNameMC);
880 mult->LoadHistograms("Multiplicity");
881
882 TFile::Open(fileNameESD);
883 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
884 multESD->LoadHistograms("Multiplicity");
885
886 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
887
447c325d 888 TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 1200);
889 canvas->Divide(3, 3);
dd701109 890
891 Int_t count = 0;
892
447c325d 893 for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
dd701109 894 {
447c325d 895 TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY("mymc");
896 mc->Sumw2();
dd701109 897
447c325d 898 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
dd701109 899 mult->ApplyMinuitFit(histID, kFALSE, type);
900 mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc);
901 TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
902
903 mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1);
904 mult->DrawComparison(Form("%s/CompareMethods_%d_Bayesian", targetDir, type), histID, kFALSE, kTRUE, mc);
905 TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
906
907 mc->GetXaxis()->SetRangeUser(0, 150);
908 chi2Result->GetXaxis()->SetRangeUser(0, 150);
909
447c325d 910/* // skip errors for now
dd701109 911 for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
912 {
913 chi2Result->SetBinError(i, 0);
914 bayesResult->SetBinError(i, 0);
447c325d 915 }*/
dd701109 916
917 canvas->cd(++count);
918 mc->SetFillColor(kYellow);
919 mc->DrawCopy();
920 chi2Result->SetLineColor(kRed);
921 chi2Result->DrawCopy("SAME");
922 bayesResult->SetLineColor(kBlue);
923 bayesResult->DrawCopy("SAME");
924 gPad->SetLogy();
925
926 canvas->cd(count + 3);
447c325d 927 chi2ResultRatio = (TH1*) chi2Result->Clone("chi2ResultRatio");
928 bayesResultRatio = (TH1*) bayesResult->Clone("bayesResultRatio");
929 chi2ResultRatio->Divide(chi2Result, mc, 1, 1, "");
930 bayesResultRatio->Divide(bayesResult, mc, 1, 1, "");
dd701109 931
447c325d 932 chi2ResultRatio->GetYaxis()->SetRangeUser(0.5, 1.5);
dd701109 933
447c325d 934 chi2ResultRatio->DrawCopy("HIST");
935 bayesResultRatio->DrawCopy("SAME HIST");
dd701109 936
447c325d 937 canvas->cd(count + 6);
938 chi2Result->Divide(chi2Result, bayesResult, 1, 1, "");
939 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
940 chi2Result->DrawCopy("HIST");
dd701109 941 }
942
943 canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
944 canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
945}
946
447c325d 947void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
dd701109 948{
949 gSystem->Load("libPWG0base");
950
951 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
952
953 TFile::Open(fileNameMC);
954 mult->LoadHistograms("Multiplicity");
955
447c325d 956 const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_200k.root", "multiplicityMC_400k.root", "multiplicityMC_600k.root", "multiplicityMC_800k.root" };
dd701109 957
0b4bfd98 958 TGraph* fitResultsChi2[3];
959 TGraph* fitResultsBayes[3];
960
961 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
962 {
963 fitResultsChi2[region] = new TGraph;
964 fitResultsChi2[region]->SetTitle(";Nevents;Chi2");
965 fitResultsChi2[region]->SetName(Form("fitResultsChi2_%d", region));
966 fitResultsChi2[region]->SetMarkerStyle(20+region);
967
968 fitResultsBayes[region] = new TGraph;
969 fitResultsBayes[region]->SetTitle(";Nevents;Chi2");
970 fitResultsBayes[region]->SetName(Form("fitResultsBayes_%d", region));
971 fitResultsBayes[region]->SetMarkerStyle(20+region);
972 fitResultsBayes[region]->SetMarkerColor(2);
973 }
974
447c325d 975 TGraph* fitResultsChi2Limit = new TGraph; fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
976 TGraph* fitResultsBayesLimit = new TGraph; fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
0b4bfd98 977 TGraph* fitResultsChi2Res = new TGraph; fitResultsChi2Res->SetTitle(";Nevents;Chi2");
978 TGraph* fitResultsBayesRes = new TGraph; fitResultsBayesRes->SetTitle(";Nevents;Chi2");
447c325d 979
447c325d 980 fitResultsChi2Limit->SetName("fitResultsChi2Limit");
981 fitResultsBayesLimit->SetName("fitResultsBayesLimit");
0b4bfd98 982 fitResultsChi2Res->SetName("fitResultsChi2Res");
983 fitResultsBayesRes->SetName("fitResultsBayesRes");
dd701109 984
985 TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600);
986 canvas->Divide(5, 2);
987
988 Float_t min = 1e10;
989 Float_t max = 0;
990
447c325d 991 TFile* file = TFile::Open("StatisticsPlot.root", "RECREATE");
992 file->Close();
993
dd701109 994 for (Int_t i=0; i<5; ++i)
995 {
996 TFile::Open(files[i]);
997 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
998 multESD->LoadHistograms("Multiplicity");
999
1000 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
1001 Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries();
447c325d 1002 TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i));
dd701109 1003
447c325d 1004 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
dd701109 1005 mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
1006 mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
1007
dd701109 1008 Int_t chi2MCLimit = 0;
0b4bfd98 1009 Float_t chi2Residuals = 0;
1010 mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals);
1011 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1012 {
1013 fitResultsChi2[region]->SetPoint(fitResultsChi2[region]->GetN(), nEntries, mult->GetQuality(region));
1014 min = TMath::Min(min, mult->GetQuality(region));
1015 max = TMath::Max(max, mult->GetQuality(region));
1016 }
dd701109 1017 fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit);
0b4bfd98 1018 fitResultsChi2Res->SetPoint(fitResultsChi2Res->GetN(), nEntries, chi2Residuals);
dd701109 1019
447c325d 1020 TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
dd701109 1021
0b4bfd98 1022 mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0, kFALSE);
dd701109 1023 mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
447c325d 1024 TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
0b4bfd98 1025 mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals);
1026 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1027 {
1028 fitResultsBayes[region]->SetPoint(fitResultsBayes[region]->GetN(), nEntries, mult->GetQuality(region));
1029 min = TMath::Min(min, mult->GetQuality(region));
1030 max = TMath::Max(max, mult->GetQuality(region));
1031 }
dd701109 1032 fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit);
0b4bfd98 1033 fitResultsBayesRes->SetPoint(fitResultsBayesRes->GetN(), nEntries, chi2Residuals);
dd701109 1034
dd701109 1035 mc->GetXaxis()->SetRangeUser(0, 150);
1036 chi2Result->GetXaxis()->SetRangeUser(0, 150);
1037
1038 // skip errors for now
1039 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
1040 {
1041 chi2Result->SetBinError(j, 0);
1042 bayesResult->SetBinError(j, 0);
1043 }
1044
1045 canvas->cd(i+1);
1046 mc->SetFillColor(kYellow);
1047 mc->DrawCopy();
1048 chi2Result->SetLineColor(kRed);
1049 chi2Result->DrawCopy("SAME");
1050 bayesResult->SetLineColor(kBlue);
1051 bayesResult->DrawCopy("SAME");
1052 gPad->SetLogy();
1053
1054 canvas->cd(i+6);
1055 chi2Result->Divide(chi2Result, mc, 1, 1, "B");
1056 bayesResult->Divide(bayesResult, mc, 1, 1, "B");
1057
1058 // skip errors for now
1059 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
1060 {
1061 chi2Result->SetBinError(j, 0);
1062 bayesResult->SetBinError(j, 0);
1063 }
1064
1065 chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
1066 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
1067
1068 chi2Result->DrawCopy("");
1069 bayesResult->DrawCopy("SAME");
447c325d 1070
1071 TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
1072 mc->Write();
1073 chi2Result->Write();
1074 bayesResult->Write();
1075 file->Close();
dd701109 1076 }
1077
1078 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1079 canvas->SaveAs(Form("%s.C", canvas->GetName()));
1080
1081 TCanvas* canvas2 = new TCanvas("StatisticsPlot2", "StatisticsPlot2", 800, 400);
1082 canvas2->Divide(2, 1);
1083
1084 canvas2->cd(1);
dd701109 1085
0b4bfd98 1086 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1087 {
1088 fitResultsChi2[region]->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
1089 fitResultsChi2[region]->Draw(((region == 0) ? "AP" : "P SAME"));
1090
1091 fitResultsBayes[region]->Draw("P SAME");
1092 }
dd701109 1093
1094 gPad->SetLogy();
1095
1096 canvas2->cd(2);
1097 fitResultsChi2Limit->SetMarkerStyle(20);
1098 fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
1099 fitResultsChi2Limit->Draw("AP");
1100
1101 fitResultsBayesLimit->SetMarkerStyle(3);
1102 fitResultsBayesLimit->SetMarkerColor(2);
1103 fitResultsBayesLimit->Draw("P SAME");
1104
1105 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
1106 canvas2->SaveAs(Form("%s.C", canvas2->GetName()));
447c325d 1107
1108 TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
0b4bfd98 1109
1110 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1111 {
1112 fitResultsChi2[region]->Write();
1113 fitResultsBayes[region]->Write();
1114 }
447c325d 1115 fitResultsChi2Limit->Write();
1116 fitResultsBayesLimit->Write();
0b4bfd98 1117 fitResultsChi2Res->Write();
1118 fitResultsBayesRes->Write();
447c325d 1119 file->Close();
dd701109 1120}
1121
447c325d 1122void StartingConditions(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
dd701109 1123{
1124 gSystem->Load("libPWG0base");
1125
1126 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1127
1128 TFile::Open(fileNameMC);
1129 mult->LoadHistograms("Multiplicity");
1130
447c325d 1131 const char* files[] = { "multiplicityMC_1M_3.root", "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root" }
dd701109 1132
1133 // this one we try to unfold
1134 TFile::Open(files[0]);
1135 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
1136 multESD->LoadHistograms("Multiplicity");
1137 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
447c325d 1138 TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY("mc");
dd701109 1139
1140 TGraph* fitResultsChi2 = new TGraph;
1141 fitResultsChi2->SetTitle(";Input Dist ID;Chi2");
1142 TGraph* fitResultsBayes = new TGraph;
1143 fitResultsBayes->SetTitle(";Input Dist ID;Chi2");
1144 TGraph* fitResultsChi2Limit = new TGraph;
1145 fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach");
1146 TGraph* fitResultsBayesLimit = new TGraph;
1147 fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach");
1148
1149 TCanvas* canvas = new TCanvas("StartingConditions", "StartingConditions", 1200, 600);
1150 canvas->Divide(8, 2);
1151
1152 TCanvas* canvas3 = new TCanvas("StartingConditions3", "StartingConditions3", 1000, 400);
1153 canvas3->Divide(2, 1);
1154
1155 Float_t min = 1e10;
1156 Float_t max = 0;
1157
1158 TH1* firstChi = 0;
1159 TH1* firstBayesian = 0;
1160 TH1* startCond = multESD->GetMultiplicityESD(histID)->ProjectionY("startCond");
1161
1162 TLegend* legend = new TLegend(0.7, 0.7, 1, 1);
1163
447c325d 1164 TFile* file = TFile::Open("StartingConditions.root", "RECREATE");
1165 mc->Write();
1166 file->Close();
1167
dd701109 1168 for (Int_t i=0; i<8; ++i)
1169 {
1170 if (i == 0)
1171 {
1172 startCond = (TH1*) mc->Clone("startCond2");
1173 }
1174 else if (i < 6)
1175 {
1176 TFile::Open(files[i-1]);
1177 AliMultiplicityCorrection* multESD2 = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2");
1178 multESD2->LoadHistograms("Multiplicity");
1179 startCond = multESD2->GetMultiplicityESD(histID)->ProjectionY("startCond");
1180 }
1181 else if (i == 6)
1182 {
1183 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])", 0, 50);
1184 func->SetParNames("scaling", "averagen", "k");
1185 func->SetParLimits(0, 1e-3, 1e10);
1186 func->SetParLimits(1, 0.001, 1000);
1187 func->SetParLimits(2, 0.001, 1000);
1188
1189 func->SetParameters(1, 10, 2);
1190 for (Int_t j=2; j<=startCond->GetNbinsX(); j++)
1191 startCond->SetBinContent(j, func->Eval(j-1));
1192 }
1193 else if (i == 7)
1194 {
1195 for (Int_t j=1; j<=startCond->GetNbinsX(); j++)
1196 startCond->SetBinContent(j, 1);
1197 }
1198
447c325d 1199 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
dd701109 1200 mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond);
1201 mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
1202
1203 Float_t chi2MC = 0;
1204 Int_t chi2MCLimit = 0;
1205 mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
1206 fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC);
1207 fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit);
1208 min = TMath::Min(min, chi2MC);
1209 max = TMath::Max(max, chi2MC);
1210
447c325d 1211 TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
dd701109 1212 if (!firstChi)
1213 firstChi = (TH1*) chi2Result->Clone("firstChi");
1214
447c325d 1215 mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, startCond);
dd701109 1216 mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
447c325d 1217 TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
dd701109 1218 if (!firstBayesian)
1219 firstBayesian = (TH1*) bayesResult->Clone("firstBayesian");
1220
1221 mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
1222 fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC);
1223 fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit);
1224
447c325d 1225 TFile* file = TFile::Open("StartingConditions.root", "UPDATE");
1226 chi2Result->Write();
1227 bayesResult->Write();
1228 file->Close();
1229
dd701109 1230 min = TMath::Min(min, chi2MC);
1231 max = TMath::Max(max, chi2MC);
1232 mc->GetXaxis()->SetRangeUser(0, 150);
1233 chi2Result->GetXaxis()->SetRangeUser(0, 150);
1234
1235 // skip errors for now
1236 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
1237 {
1238 chi2Result->SetBinError(j, 0);
1239 bayesResult->SetBinError(j, 0);
1240 }
1241
1242 canvas3->cd(1);
1243 TH1* tmp = (TH1*) chi2Result->Clone("tmp");
1244 tmp->SetTitle("Difference to best initial conditions;Npart;Ratio");
1245 tmp->Divide(firstChi);
1246 tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
1247 tmp->GetXaxis()->SetRangeUser(0, 200);
1248 tmp->SetLineColor(i+1);
1249 legend->AddEntry(tmp, Form("%d", i));
1250 tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
1251
1252 canvas3->cd(2);
1253 tmp = (TH1*) bayesResult->Clone("tmp");
1254 tmp->SetTitle("Difference to best initial conditions;Npart;Ratio");
1255 tmp->Divide(firstBayesian);
1256 tmp->SetLineColor(i+1);
1257 tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
1258 tmp->GetXaxis()->SetRangeUser(0, 200);
1259 tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
1260
1261 canvas->cd(i+1);
1262 mc->SetFillColor(kYellow);
1263 mc->DrawCopy();
1264 chi2Result->SetLineColor(kRed);
1265 chi2Result->DrawCopy("SAME");
1266 bayesResult->SetLineColor(kBlue);
1267 bayesResult->DrawCopy("SAME");
1268 gPad->SetLogy();
1269
1270 canvas->cd(i+9);
1271 chi2Result->Divide(chi2Result, mc, 1, 1, "B");
1272 bayesResult->Divide(bayesResult, mc, 1, 1, "B");
1273
1274 // skip errors for now
1275 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
1276 {
1277 chi2Result->SetBinError(j, 0);
1278 bayesResult->SetBinError(j, 0);
1279 }
1280
1281 chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
1282 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
1283
1284 chi2Result->DrawCopy("");
1285 bayesResult->DrawCopy("SAME");
1286 }
1287
1288 canvas3->cd(1);
1289 legend->Draw();
1290
cfc19dd5 1291 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
dd701109 1292
1293 TCanvas* canvas2 = new TCanvas("StartingConditions2", "StartingConditions2", 800, 400);
1294 canvas2->Divide(2, 1);
1295
1296 canvas2->cd(1);
1297 fitResultsChi2->SetMarkerStyle(20);
1298 fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
1299 fitResultsChi2->Draw("AP");
1300
1301 fitResultsBayes->SetMarkerStyle(3);
1302 fitResultsBayes->SetMarkerColor(2);
1303 fitResultsBayes->Draw("P SAME");
1304
1305 gPad->SetLogy();
1306
1307 canvas2->cd(2);
1308 fitResultsChi2Limit->SetMarkerStyle(20);
1309 fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
1310 fitResultsChi2Limit->Draw("AP");
1311
1312 fitResultsBayesLimit->SetMarkerStyle(3);
1313 fitResultsBayesLimit->SetMarkerColor(2);
1314 fitResultsBayesLimit->Draw("P SAME");
1315
1316 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
1317 canvas3->SaveAs(Form("%s.gif", canvas3->GetName()));
cfc19dd5 1318}
1319
447c325d 1320void DifferentSamples(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
1321{
1322 gSystem->Load("libPWG0base");
1323
1324 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1325
1326 TFile::Open(fileNameMC);
1327 mult->LoadHistograms("Multiplicity");
1328
1329 const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" };
1330
1331 TGraph* fitResultsChi2 = new TGraph;
1332 fitResultsChi2->SetTitle(";Input Dist ID;Chi2");
1333 TGraph* fitResultsBayes = new TGraph;
1334 fitResultsBayes->SetTitle(";Input Dist ID;Chi2");
1335 TGraph* fitResultsChi2Limit = new TGraph;
1336 fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach");
1337 TGraph* fitResultsBayesLimit = new TGraph;
1338 fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach");
1339
1340 TCanvas* canvasA = new TCanvas("DifferentSamplesA", "DifferentSamplesA", 1200, 600);
1341 canvasA->Divide(4, 2);
1342
1343 TCanvas* canvasB = new TCanvas("DifferentSamplesB", "DifferentSamplesB", 1200, 600);
1344 canvasB->Divide(4, 2);
1345
1346 TCanvas* canvas4 = new TCanvas("DifferentSamples4", "DifferentSamples4", 1000, 400);
1347 canvas4->Divide(2, 1);
1348
1349 TCanvas* canvas3 = new TCanvas("DifferentSamples3", "DifferentSamples3", 1000, 400);
1350 canvas3->Divide(2, 1);
1351
1352 Float_t min = 1e10;
1353 Float_t max = 0;
1354
1355 TH1* firstChi = 0;
1356 TH1* firstBayesian = 0;
1357
1358 TLegend* legend = new TLegend(0.7, 0.7, 1, 1);
1359
1360 TFile* file = TFile::Open("DifferentSamples.root", "RECREATE");
1361 file->Close();
1362
1363 for (Int_t i=0; i<8; ++i)
1364 {
1365 TFile::Open(files[i]);
1366 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2");
1367 multESD->LoadHistograms("Multiplicity");
1368 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
1369 TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i));
1370 mc->Sumw2();
1371
1372 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1373 mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1374 mult->DrawComparison(Form("DifferentSamples_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
1375
1376 Float_t chi2MC = 0;
1377 Int_t chi2MCLimit = 0;
1378 mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
1379 fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC);
1380 fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit);
1381 min = TMath::Min(min, chi2MC);
1382 max = TMath::Max(max, chi2MC);
1383
1384 TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
1385 if (!firstChi)
1386 firstChi = (TH1*) chi2Result->Clone("firstChi");
1387
1388 mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1389 mult->DrawComparison(Form("DifferentSamples_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
1390 TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
1391 if (!firstBayesian)
1392 firstBayesian = (TH1*) bayesResult->Clone("firstBayesian");
1393
1394 TFile* file = TFile::Open("DifferentSamples.root", "UPDATE");
1395 mc->Write();
1396 chi2Result->Write();
1397 bayesResult->Write();
1398 file->Close();
1399
1400 mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
1401 fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC);
1402 fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit);
1403
1404 min = TMath::Min(min, chi2MC);
1405 max = TMath::Max(max, chi2MC);
1406 mc->GetXaxis()->SetRangeUser(0, 150);
1407 chi2Result->GetXaxis()->SetRangeUser(0, 150);
1408
1409 // skip errors for now
1410 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
1411 {
1412 chi2Result->SetBinError(j, 0);
1413 bayesResult->SetBinError(j, 0);
1414 }
1415
1416 canvas4->cd(1);
1417 TH1* tmp = (TH1*) chi2Result->Clone("tmp");
1418 tmp->SetTitle("Unfolded/MC;Npart;Ratio");
1419 tmp->Divide(mc);
1420 tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
1421 tmp->GetXaxis()->SetRangeUser(0, 200);
1422 tmp->SetLineColor(i+1);
1423 tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
1424
1425 canvas4->cd(2);
1426 tmp = (TH1*) bayesResult->Clone("tmp");
1427 tmp->SetTitle("Unfolded/MC;Npart;Ratio");
1428 tmp->Divide(mc);
1429 tmp->SetLineColor(i+1);
1430 tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
1431 tmp->GetXaxis()->SetRangeUser(0, 200);
1432 tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
1433
1434 canvas3->cd(1);
1435 TH1* tmp = (TH1*) chi2Result->Clone("tmp");
1436 tmp->SetTitle("Ratio to first result;Npart;Ratio");
1437 tmp->Divide(firstChi);
1438 tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
1439 tmp->GetXaxis()->SetRangeUser(0, 200);
1440 tmp->SetLineColor(i+1);
1441 legend->AddEntry(tmp, Form("%d", i));
1442 tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
1443
1444 canvas3->cd(2);
1445 tmp = (TH1*) bayesResult->Clone("tmp");
1446 tmp->SetTitle("Ratio to first result;Npart;Ratio");
1447 tmp->Divide(firstBayesian);
1448 tmp->SetLineColor(i+1);
1449 tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
1450 tmp->GetXaxis()->SetRangeUser(0, 200);
1451 tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
1452
1453 if (i < 4)
1454 {
1455 canvasA->cd(i+1);
1456 }
1457 else
1458 canvasB->cd(i+1-4);
1459
1460 mc->SetFillColor(kYellow);
1461 mc->DrawCopy();
1462 chi2Result->SetLineColor(kRed);
1463 chi2Result->DrawCopy("SAME");
1464 bayesResult->SetLineColor(kBlue);
1465 bayesResult->DrawCopy("SAME");
1466 gPad->SetLogy();
1467
1468 if (i < 4)
1469 {
1470 canvasA->cd(i+5);
1471 }
1472 else
1473 canvasB->cd(i+5-4);
1474
1475 chi2Result->Divide(chi2Result, mc, 1, 1, "B");
1476 bayesResult->Divide(bayesResult, mc, 1, 1, "B");
1477
1478 // skip errors for now
1479 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
1480 {
1481 chi2Result->SetBinError(j, 0);
1482 bayesResult->SetBinError(j, 0);
1483 }
1484
1485 chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
1486 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
1487
1488 chi2Result->DrawCopy("");
1489 bayesResult->DrawCopy("SAME");
1490 }
1491
1492 canvas3->cd(1);
1493 legend->Draw();
1494
1495 canvasA->SaveAs(Form("%s.gif", canvasA->GetName()));
1496 canvasB->SaveAs(Form("%s.gif", canvasB->GetName()));
1497
1498 TCanvas* canvas2 = new TCanvas("DifferentSamples2", "DifferentSamples2", 800, 400);
1499 canvas2->Divide(2, 1);
1500
1501 canvas2->cd(1);
1502 fitResultsChi2->SetMarkerStyle(20);
1503 fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
1504 fitResultsChi2->Draw("AP");
1505
1506 fitResultsBayes->SetMarkerStyle(3);
1507 fitResultsBayes->SetMarkerColor(2);
1508 fitResultsBayes->Draw("P SAME");
1509
1510 gPad->SetLogy();
1511
1512 canvas2->cd(2);
1513 fitResultsChi2Limit->SetMarkerStyle(20);
1514 fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
1515 fitResultsChi2Limit->Draw("AP");
1516
1517 fitResultsBayesLimit->SetMarkerStyle(3);
1518 fitResultsBayesLimit->SetMarkerColor(2);
1519 fitResultsBayesLimit->Draw("P SAME");
1520
1521 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
1522 canvas3->SaveAs(Form("%s.gif", canvas3->GetName()));
1523 canvas4->SaveAs(Form("%s.gif", canvas4->GetName()));
1524}
1525
cfc19dd5 1526void Merge(Int_t n, const char** files, const char* output)
1527{
447c325d 1528 // const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" };
1529
1530
cfc19dd5 1531 gSystem->Load("libPWG0base");
1532
1533 AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n];
1534 TList list;
1535 for (Int_t i=0; i<n; ++i)
1536 {
1537 TString name("Multiplicity");
1538 if (i > 0)
1539 name.Form("Multiplicity%d", i);
1540
1541 TFile::Open(files[i]);
1542 data[i] = new AliMultiplicityCorrection(name, name);
1543 data[i]->LoadHistograms("Multiplicity");
1544 if (i > 0)
1545 list.Add(data[i]);
1546 }
1547
1548 data[0]->Merge(&list);
1549
447c325d 1550 //data[0]->DrawHistograms();
cfc19dd5 1551
1552 TFile::Open(output, "RECREATE");
1553 data[0]->SaveHistograms();
1554 gFile->Close();
1555}
1556
1557void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
1558{
1559 gSystem->Load("libPWG0base");
1560
1561 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1562
1563 TFile::Open(fileName);
1564 mult->LoadHistograms("Multiplicity");
1565
1566 TF1* func = 0;
1567
1568 if (caseNo >= 4)
1569 {
0b4bfd98 1570 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])", 0, 500);
cfc19dd5 1571 func->SetParNames("scaling", "averagen", "k");
1572 }
1573
1574 switch (caseNo)
1575 {
0b4bfd98 1576 case 0: func = new TF1("flat", "1000"); break;
cfc19dd5 1577 case 1: func = new TF1("flat", "501-x"); break;
1578 case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
1579 case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
dd701109 1580 case 4: func->SetParameters(1e7, 10, 2); break;
447c325d 1581 case 5: func->SetParameters(1, 13, 7); break;
dd701109 1582 case 6: func->SetParameters(1e7, 30, 4); break;
0b4bfd98 1583 case 7: func->SetParameters(1e7, 30, 2); break; // ***
dd701109 1584 case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break;
cfc19dd5 1585
1586 default: return;
1587 }
1588
447c325d 1589 new TCanvas;
1590 func->Draw();
1591
1592 mult->SetGenMeasFromFunc(func, 3);
cfc19dd5 1593
dd701109 1594 TFile::Open("out.root", "RECREATE");
1595 mult->SaveHistograms();
1596
0b4bfd98 1597 new TCanvas; mult->GetMultiplicityESD(3)->ProjectionY()->DrawCopy();
1598 new TCanvas; mult->GetMultiplicityVtx(3)->ProjectionY()->DrawCopy();
1599
dd701109 1600 //mult->ApplyBayesianMethod(2, kFALSE);
1601 //mult->ApplyMinuitFit(2, kFALSE);
cfc19dd5 1602 //mult->ApplyGaussianMethod(2, kFALSE);
447c325d 1603 //mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx);
dd701109 1604}
1605
1606void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t corrMatrix = 2)
1607{
1608 gSystem->Load("libPWG0base");
1609
1610 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1611
1612 TFile::Open(fileName);
1613 mult->LoadHistograms("Multiplicity");
1614
1615 // empty under/overflow bins in x, otherwise Project3D takes them into account
1616 TH3* corr = mult->GetCorrelation(corrMatrix);
1617 for (Int_t j=1; j<=corr->GetYaxis()->GetNbins(); ++j)
1618 {
1619 for (Int_t k=1; k<=corr->GetZaxis()->GetNbins(); ++k)
1620 {
1621 corr->SetBinContent(0, j, k, 0);
1622 corr->SetBinContent(corr->GetXaxis()->GetNbins()+1, j, k, 0);
1623 }
1624 }
1625
1626 TH2* proj = (TH2*) corr->Project3D("zy");
1627
1628 // normalize correction for given nPart
1629 for (Int_t i=1; i<=proj->GetNbinsX(); ++i)
1630 {
1631 Double_t sum = proj->Integral(i, i, 1, proj->GetNbinsY());
1632 if (sum <= 0)
1633 continue;
1634
1635 for (Int_t j=1; j<=proj->GetNbinsY(); ++j)
1636 {
1637 // npart sum to 1
1638 proj->SetBinContent(i, j, proj->GetBinContent(i, j) / sum);
1639 proj->SetBinError(i, j, proj->GetBinError(i, j) / sum);
1640 }
1641 }
1642
1643 new TCanvas;
447c325d 1644 proj->DrawCopy("COLZ");
dd701109 1645
1646 TH1* scaling = proj->ProjectionY("scaling", 1, 1);
1647 scaling->Reset();
1648 scaling->SetMarkerStyle(3);
1649 //scaling->GetXaxis()->SetRangeUser(0, 50);
1650 TH1* mean = (TH1F*) scaling->Clone("mean");
1651 TH1* width = (TH1F*) scaling->Clone("width");
1652
1653 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
1654 lognormal->SetParNames("scaling", "mean", "sigma");
1655 lognormal->SetParameters(1, 1, 1);
1656 lognormal->SetParLimits(0, 1, 1);
1657 lognormal->SetParLimits(1, 0, 100);
1658 lognormal->SetParLimits(2, 1e-3, 1);
1659
1660 TF1* nbd = 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])", 0, 50);
1661 nbd->SetParNames("scaling", "averagen", "k");
1662 nbd->SetParameters(1, 13, 5);
1663 nbd->SetParLimits(0, 1, 1);
1664 nbd->SetParLimits(1, 1, 100);
1665 nbd->SetParLimits(2, 1, 1e8);
1666
1667 TF1* poisson = new TF1("poisson", "[0] * exp(-(x+[2])) * (x+[2])**[1] / TMath::Factorial([1])", 0.01, 50);
1668 poisson->SetParNames("scaling", "k", "deltax");
1669 poisson->SetParameters(1, 1, 0);
1670 poisson->SetParLimits(0, 0, 10);
1671 poisson->SetParLimits(1, 0.01, 100);
1672 poisson->SetParLimits(2, 0, 10);
1673
1674 TF1* mygaus = new TF1("mygaus", "[0] * exp(-(x-[1])**2 / 2 / [2] - [3] * log(x + [4]) / [5])", 0.01, 50);
1675 mygaus->SetParNames("scaling", "mean", "width", "scale2log", "logmean", "logwidth");
1676 mygaus->SetParameters(1, 0, 1, 1, 0, 1);
1677 mygaus->SetParLimits(2, 1e-5, 10);
1678 mygaus->SetParLimits(4, 1, 1);
1679 mygaus->SetParLimits(5, 1e-5, 10);
1680
1681 //TF1* sqrt = new TF1("sqrt", "[0] + [1] * sqrt((x + [3]) * [2])", 0, 50);
1682 TF1* sqrt = new TF1("sqrt", "[0] + (x + [1])**[2]", 0, 50);
1683 sqrt->SetParNames("ydelta", "exp", "xdelta");
1684 sqrt->SetParameters(0, 0, 1);
1685 sqrt->SetParLimits(1, 0, 10);
1686
1687 const char* fitWith = "gaus";
1688
1689 for (Int_t i=1; i<=150; ++i)
1690 {
1691 printf("Fitting %d...\n", i);
1692
1693 TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e");
447c325d 1694
dd701109 1695 //hist->GetXaxis()->SetRangeUser(0, 50);
1696 //lognormal->SetParameter(0, hist->GetMaximum());
1697 hist->Fit(fitWith, "0 M", "");
1698
1699 TF1* func = hist->GetListOfFunctions()->FindObject(fitWith);
1700
447c325d 1701 if (0 && (i % 5 == 0))
dd701109 1702 {
447c325d 1703 pad = new TCanvas;
dd701109 1704 hist->Draw();
1705 func->Clone()->Draw("SAME");
447c325d 1706 pad->SetLogy();
dd701109 1707 }
1708
1709 scaling->Fill(i, func->GetParameter(0));
1710 mean->Fill(i, func->GetParameter(1));
1711 width->Fill(i, func->GetParameter(2));
1712 }
1713
1714 TF1* log = new TF1("log", "[0] + [1] * log([2] * x)", 0.01, 500);
1715 log->SetParameters(0, 1, 1);
1716 log->SetParLimits(1, 0, 100);
1717 log->SetParLimits(2, 1e-3, 10);
1718
1719 TF1* over = new TF1("over", "[0] + [1] / (x+[2])", 0.01, 500);
1720 over->SetParameters(0, 1, 0);
1721 //over->SetParLimits(0, 0, 100);
1722 over->SetParLimits(1, 1e-3, 10);
1723 over->SetParLimits(2, 0, 100);
1724
1725 c1 = new TCanvas("fitparams", "fitparams", 1200, 400);
1726 c1->Divide(3, 1);
1727
1728 c1->cd(1);
1729 scaling->Draw("P");
1730
1731 //TF1* scalingFit = new TF1("mypol0", "[0]");
1732 TF1* scalingFit = over;
447c325d 1733 scaling->Fit(scalingFit, "", "", 3, 140);
1734 scalingFit->SetRange(0, 200);
1735 scalingFit->Draw("SAME");
dd701109 1736
1737 c1->cd(2);
1738 mean->Draw("P");
1739
1740 //TF1* meanFit = log;
1741 TF1* meanFit = new TF1("mypol1", "[0]+[1]*x");
447c325d 1742 mean->Fit(meanFit, "", "", 3, 140);
1743 meanFit->SetRange(0, 200);
1744 meanFit->Draw("SAME");
dd701109 1745
1746 c1->cd(3);
1747 width->Draw("P");
1748
1749 //TF1* widthFit = over;
447c325d 1750 TF1* widthFit = new TF1("mypol", "[0]+[1]*TMath::Sqrt([2]*x)");
1751 widthFit->SetParLimits(2, 1e-5, 1e5);
1752 width->Fit(widthFit, "", "", 5, 140);
1753 widthFit->SetRange(0, 200);
1754 widthFit->Draw("SAME");
dd701109 1755
1756 // build new correction matrix
1757 TH2* new = (TH2*) proj->Clone("new");
1758 new->Reset();
1759 Float_t x, y;
1760 for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1)
1761 {
1762 TF1* func = (TF1*) gROOT->FindObject(fitWith);
1763 x = new->GetXaxis()->GetBinCenter(i);
1764 //if (i == 1)
1765 // x = 0.1;
1766 x++;
1767 func->SetParameters(scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x));
1768 printf("%f %f %f %f\n", x, scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x));
1769
1770 for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
1771 {
447c325d 1772 if (i < 11)
dd701109 1773 {
1774 // leave bins 1..20 untouched
1775 new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j));
1776 }
1777 else
1778 {
1779 y = new->GetYaxis()->GetBinCenter(j);
1780 if (j == 1)
1781 y = 0.1;
1782 if (func->Eval(y) > 1e-4)
1783 new->SetBinContent(i, j, func->Eval(y));
1784 }
1785 }
1786 }
1787
1788 // fill 0 multiplicity bins, this cannot be done with the function because it does not accept 0
1789 // we take the values from the old response matrix
1790 //for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1)
1791 // new->SetBinContent(i, 1, proj->GetBinContent(i, 1));
1792
1793 //for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
1794 // new->SetBinContent(1, j, proj->GetBinContent(1, j));
1795
1796 // normalize correction for given nPart
1797 for (Int_t i=1; i<=new->GetNbinsX(); ++i)
1798 {
1799 Double_t sum = new->Integral(i, i, 1, proj->GetNbinsY());
1800 if (sum <= 0)
1801 continue;
1802
1803 for (Int_t j=1; j<=new->GetNbinsY(); ++j)
1804 {
1805 // npart sum to 1
1806 new->SetBinContent(i, j, new->GetBinContent(i, j) / sum);
1807 new->SetBinError(i, j, new->GetBinError(i, j) / sum);
1808 }
1809 }
1810
1811 new TCanvas;
1812 new->Draw("COLZ");
1813
1814 TH2* diff = (TH2*) new->Clone("diff");
1815 diff->Add(proj, -1);
1816
1817 new TCanvas;
1818 diff->Draw("COLZ");
1819 diff->SetMinimum(-0.05);
1820 diff->SetMaximum(0.05);
1821
1822 corr->Reset();
1823
1824 for (Int_t i=1; i<=new->GetNbinsX(); ++i)
1825 for (Int_t j=1; j<=new->GetNbinsY(); ++j)
1826 corr->SetBinContent(corr->GetXaxis()->GetNbins() / 2, i, j, new->GetBinContent(i, j));
1827
1828 new TCanvas;
1829 corr->Project3D("zy")->Draw("COLZ");
1830
1831 TFile::Open("out.root", "RECREATE");
1832 mult->SaveHistograms();
1833
1834 TH1* proj1 = proj->ProjectionY("proj1", 36, 36);
1835 TH1* proj2 = new->ProjectionY("proj2", 36, 36);
1836 proj2->SetLineColor(2);
1837
1838 new TCanvas;
1839 proj1->Draw();
1840 proj2->Draw("SAME");
cfc19dd5 1841}
447c325d 1842
0b4bfd98 1843void buildCorrelationMap(const char* fileName = "multiplicityMC_2M.root", Int_t corrMatrix = 3)
1844{
1845 gSystem->Load("libPWG0base");
1846
1847 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1848
1849 TFile::Open(fileName);
1850 mult->LoadHistograms("Multiplicity");
1851
1852 TH3F* new = mult->GetCorrelation(corrMatrix);
1853 new->Reset();
1854
1855 TF1* func = new TF1("func", "gaus(0)");
1856
1857 Int_t vtxBin = new->GetNbinsX() / 2;
1858 if (vtxBin == 0)
1859 vtxBin = 1;
1860
1861 Float_t sigma = 2;
1862 for (Int_t i=1; i<=new->GetYaxis()->GetNbins(); i+=1)
1863 {
1864 Float_t x = new->GetYaxis()->GetBinCenter(i);
1865 func->SetParameters(1, x * 0.8, sigma);
1866 //func->SetParameters(1, x, sigma);
1867
1868 for (Int_t j=1; j<=new->GetZaxis()->GetNbins(); j+=1)
1869 {
1870 Float_t y = new->GetYaxis()->GetBinCenter(j);
1871
1872 // cut at 1 sigma
1873 if (TMath::Abs(y-x*0.8) < sigma)
1874 new->SetBinContent(vtxBin, i, j, func->Eval(y));
1875
1876 // test only bin 40 has smearing
1877 //if (x != 40)
1878 // new->SetBinContent(vtxBin, i, j, (i == j));
1879 }
1880 }
1881
1882 new TCanvas;
1883 new->Project3D("zy")->DrawCopy("COLZ");
1884
1885 TFile* file = TFile::Open("out.root", "RECREATE");
1886 mult->SetCorrelation(corrMatrix, new);
1887 mult->SaveHistograms();
1888 file->Close();
1889}
1890
447c325d 1891void GetCrossSections(const char* fileName)
1892{
1893 gSystem->Load("libPWG0base");
1894
1895 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1896
1897 TFile::Open(fileName);
1898 mult->LoadHistograms("Multiplicity");
1899
1900 TH1* xSection2 = mult->GetCorrelation(3)->Project3D("y")->Clone("xSection2");
1901 xSection2->Sumw2();
1902 xSection2->Scale(1.0 / xSection2->Integral());
1903
1904 TH1* xSection15 = mult->GetCorrelation(2)->Project3D("y")->Clone("xSection15");
1905 xSection15->Sumw2();
1906 xSection15->Scale(1.0 / xSection15->Integral());
1907
1908 TFile::Open("crosssection.root", "RECREATE");
1909 xSection2->Write();
1910 xSection15->Write();
1911 gFile->Close();
1912}
0b4bfd98 1913
44466df2 1914void AnalyzeSpeciesTree(const char* fileName)
1915{
1916 //
1917 // prints statistics about fParticleSpecies
1918 //
1919
1920 gSystem->Load("libPWG0base");
1921
1922 TFile::Open(fileName);
1923 TNtuple* fParticleSpecies = (TNtuple*) gFile->Get("fParticleSpecies");
1924
1925 const Int_t nFields = 8;
1926 Long_t totals[8];
1927 for (Int_t i=0; i<nFields; i++)
1928 totals[i] = 0;
1929
1930 for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++)
1931 {
1932 fParticleSpecies->GetEvent(i);
1933
1934 Float_t* f = fParticleSpecies->GetArgs();
1935
1936 for (Int_t j=0; j<nFields; j++)
1937 totals[j] += f[j+1];
1938 }
1939
1940 for (Int_t i=0; i<nFields; i++)
1941 Printf("%d --> %ld", i, totals[i]);
1942}
1943
0b4bfd98 1944void BuildResponseFromTree(const char* fileName, const char* target)
1945{
1946 //
1947 // builds several response matrices with different particle ratios (systematic study)
1948 //
1949
1950 gSystem->Load("libPWG0base");
1951
1952 TFile::Open(fileName);
1953 TNtuple* fParticleSpecies = (TNtuple*) gFile->Get("fParticleSpecies");
1954
1955 TFile* file = TFile::Open(target, "RECREATE");
1956 file->Close();
1957
1958 Int_t tracks = 0; // control variables
1959 Int_t noLabel = 0;
6d81c2de 1960 Int_t secondaries = 0;
0b4bfd98 1961 Int_t doubleCount = 0;
1962
44466df2 1963 for (Int_t num = 0; num < 7; num++)
0b4bfd98 1964 {
1965 AliMultiplicityCorrection* fMultiplicity = new AliMultiplicityCorrection(Form("Multiplicity_%d", num), Form("Multiplicity_%d", num));
1966
1967 Float_t ratio[4]; // pi, K, p, other
1968 for (Int_t i = 0; i < 4; i++)
1969 ratio[i] = 1;
1970
1971 switch (num)
1972 {
1973 case 1 : ratio[1] = 0.5; break;
1974 case 2 : ratio[2] = 0.5; break;
1975 case 3 : ratio[1] = 1.5; break;
1976 case 4 : ratio[2] = 1.5; break;
1977 case 5 : ratio[1] = 0.5; ratio[2] = 0.5; break;
1978 case 6 : ratio[1] = 1.5; ratio[2] = 1.5; break;
1979 }
1980
1981 for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++)
1982 {
1983 fParticleSpecies->GetEvent(i);
1984
1985 Float_t* f = fParticleSpecies->GetArgs();
1986
1987 Float_t gene = 0;
1988 Float_t meas = 0;
1989
1990 for (Int_t j = 0; j < 4; j++)
1991 {
1992 gene += ratio[j] * f[j+1];
1993 meas += ratio[j] * f[j+1+4];
1994 tracks += f[j+1+4];
1995 }
1996
1997 // add the ones w/o label
1998 tracks += f[9];
1999 noLabel += f[9];
2000
6d81c2de 2001 // secondaries are already part of meas!
2002 secondaries += f[10];
2003
2004 // double counted are already part of meas!
2005 doubleCount += f[11];
0b4bfd98 2006
6d81c2de 2007 // ones w/o label are added without weighting to allow comparison to default analysis. however this is only valid when their fraction is low enough!
0b4bfd98 2008 meas += f[9];
0b4bfd98 2009
2010 //Printf("%.f %.f %.f %.f %.f", f[5], f[6], f[7], f[8], f[9]);
2011
6d81c2de 2012 fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, 0, meas, meas, meas, meas);
2013 fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, gene, gene, gene, gene, 0);
2014 fMultiplicity->FillMeasured(f[0], meas, meas, meas, meas);
0b4bfd98 2015 }
2016
2017 //fMultiplicity->DrawHistograms();
2018
2019 TFile* file = TFile::Open(target, "UPDATE");
2020 fMultiplicity->SaveHistograms();
2021 file->Close();
2022
2023 if (num == 0)
6d81c2de 2024 {
2025 Printf("%d total tracks, %d w/o label = %.2f %%, %d double counted = %.2f %%, secondaries = %.2f %%", tracks, noLabel, 100.0 * noLabel / tracks, doubleCount, 100.0 * doubleCount / tracks, 100.0 * secondaries / tracks);
2026 if ((Float_t) noLabel / tracks > 0.02)
2027 Printf("WARNING: More than 2%% of tracks without label, this might bias the study!");
2028 }
0b4bfd98 2029 }
2030}
2031
2032void MergeModifyCrossSection(const char* output)
2033{
2034 const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root" };
2035
2036 gSystem->Load("libPWG0base");
2037
2038 TFile::Open(output, "RECREATE");
2039 gFile->Close();
2040
2041 for (Int_t num=0; num<7; ++num)
2042 {
2043 AliMultiplicityCorrection* data[3];
2044 TList list;
2045
2046 Float_t ratio[3];
2047 switch (num)
2048 {
2049 case 0: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.0; break;
2050 case 1: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.0; break;
2051 case 2: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 1.0; break;
2052 case 3: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.5; break;
2053 case 4: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 0.5; break;
2054 case 5: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.5; break;
2055 case 6: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 0.5; break;
2056 default: return;
2057 }
2058
2059 for (Int_t i=0; i<3; ++i)
2060 {
2061 TString name;
2062 name.Form("Multiplicity_%d", num);
2063 if (i > 0)
2064 name.Form("Multiplicity_%d_%d", num, i);
2065
2066 TFile::Open(files[i]);
2067 data[i] = new AliMultiplicityCorrection(name, name);
2068 data[i]->LoadHistograms("Multiplicity");
2069
2070 // modify x-section
2071 for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++)
2072 {
2073 data[i]->GetMultiplicityVtx(j)->Scale(ratio[i]);
2074 data[i]->GetMultiplicityMB(j)->Scale(ratio[i]);
2075 data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]);
2076 }
2077
2078 for (Int_t j=0; j<AliMultiplicityCorrection::kESDHists; j++)
2079 data[i]->GetMultiplicityESD(j)->Scale(ratio[i]);
2080
2081 for (Int_t j=0; j<AliMultiplicityCorrection::kCorrHists; j++)
2082 data[i]->GetCorrelation(j)->Scale(ratio[i]);
2083
2084 if (i > 0)
2085 list.Add(data[i]);
2086 }
2087
2088 printf("Case %d, %s: Entries in response matrix 3: ND: %.2f SD: %.2f DD: %.2f", num, data[0]->GetName(), data[0]->GetCorrelation(3)->Integral(), data[1]->GetCorrelation(3)->Integral(), data[2]->GetCorrelation(3)->Integral());
2089
2090 data[0]->Merge(&list);
2091
2092 Printf(" Total: %.2f", data[0]->GetCorrelation(3)->Integral());
2093
2094 TFile::Open(output, "UPDATE");
2095 data[0]->SaveHistograms();
2096 gFile->Close();
2097
2098 list.Clear();
2099
2100 for (Int_t i=0; i<3; ++i)
2101 delete data[i];
2102 }
2103}
44466df2 2104
2105void Rebin(const char* fileName = "multiplicityMC_3M.root", Int_t corrMatrix = 3)
2106{
2107 // rebins MC axis of correlation map, MC and histogram for corrected (for evaluation of effect of regularization)
2108 // rebin does not exist for 3D hists, so we convert to 2D and then back to 3D (loosing the vertex information)
2109
2110 Printf("WARNING: Vertex information is lost in this process. Use result only for evaluation of errors.");
2111
2112 gSystem->Load("libPWG0base");
2113
2114 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2115
2116 TFile::Open(fileName);
2117 mult->LoadHistograms("Multiplicity");
2118
2119 // rebin correlation
2120 TH3* old = mult->GetCorrelation(corrMatrix);
2121
2122 // empty under/overflow bins in x, otherwise Project3D takes them into account
2123 for (Int_t y=1; y<=old->GetYaxis()->GetNbins(); ++y)
2124 {
2125 for (Int_t z=1; z<=old->GetZaxis()->GetNbins(); ++z)
2126 {
2127 old->SetBinContent(0, y, z, 0);
2128 old->SetBinContent(old->GetXaxis()->GetNbins()+1, y, z, 0);
2129 }
2130 }
2131
2132 TH2* response = (TH2*) old->Project3D("zy");
2133 response->RebinX(2);
2134
2135 TH3F* new = new TH3F(old->GetName(), old->GetTitle(),
2136 old->GetXaxis()->GetNbins(), old->GetXaxis()->GetBinLowEdge(1), old->GetXaxis()->GetBinUpEdge(old->GetXaxis()->GetNbins()),
2137 old->GetYaxis()->GetNbins() / 2, old->GetYaxis()->GetBinLowEdge(1), old->GetYaxis()->GetBinUpEdge(old->GetYaxis()->GetNbins()),
2138 old->GetZaxis()->GetNbins(), old->GetZaxis()->GetBinLowEdge(1), old->GetZaxis()->GetBinUpEdge(old->GetZaxis()->GetNbins()));
2139 new->Reset();
2140
2141 Int_t vtxBin = new->GetNbinsX() / 2;
2142 if (vtxBin == 0)
2143 vtxBin = 1;
2144
2145 for (Int_t i=1; i<=new->GetYaxis()->GetNbins(); i+=1)
2146 for (Int_t j=1; j<=new->GetZaxis()->GetNbins(); j+=1)
2147 new->SetBinContent(vtxBin, i, j, response->GetBinContent(i, j));
2148
2149 // rebin MC + hist for corrected
2150 for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kTrVtx; eventType <= AliMultiplicityCorrection::kINEL; eventType++)
2151 mult->GetMultiplicityMC(corrMatrix, eventType)->RebinY(2);
2152
2153 mult->GetMultiplicityESDCorrected(corrMatrix)->Rebin(2);
2154
2155 // recreate measured from correlation matrix to get rid of vertex shift effect
2156 TH2* newMeasured = (TH2*) old->Project3D("zx");
2157 TH2* esd = mult->GetMultiplicityESD(corrMatrix);
2158 esd->Reset();
2159
2160 // transfer from TH2D to TH2F
2161 for (Int_t i=0; i<=new->GetXaxis()->GetNbins()+1; i+=1)
2162 for (Int_t j=0; j<=new->GetYaxis()->GetNbins()+1; j+=1)
2163 esd->SetBinContent(i, j, newMeasured->GetBinContent(i, j));
2164
2165 new TCanvas;
2166 new->Project3D("zy")->DrawCopy("COLZ");
2167
2168 TFile* file = TFile::Open("out.root", "RECREATE");
2169 mult->SetCorrelation(corrMatrix, new);
2170 mult->SaveHistograms();
2171 file->Close();
2172}
2173
2174void EvaluateRegularizationEffect(Int_t step, const char* fileNameRebinned = "multiplicityMC_3M_rebinned.root", const char* fileNameNormal = "multiplicityMC_3M.root", Int_t histID = 3)
2175{
2176 // due to some static members in AliMultiplicityCorrection, the session has to be restarted after changing the number of parameters, to be fixed
2177 // that is why this is done in 2 steps
2178
2179 gSystem->Load("libPWG0base");
2180
2181 Bool_t fullPhaseSpace = kFALSE;
2182
2183 if (step == 1)
2184 {
2185 // first step: unfold without regularization and rebinned histogram ("N=M")
2186 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2187 TFile::Open(fileNameRebinned);
2188 mult->LoadHistograms();
2189
2190 mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125);
2191 mult->SetCreateBigBin(kFALSE);
2192
2193 mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
2194 mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, mult->GetMultiplicityVtx(histID)->ProjectionY("mymchist"));
2195
2196 TFile* file = TFile::Open("EvaluateRegularizationEffect1.root", "RECREATE");
2197 mult->SaveHistograms();
2198 file->Close();
2199 }
2200 else if (step == 2)
2201 {
2202 // second step: unfold with regularization and normal histogram
2203 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2204 TFile::Open(fileNameNormal);
2205 mult2->LoadHistograms();
2206
2207 mult2->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
2208 mult2->SetCreateBigBin(kTRUE);
2209 mult2->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
2210 mult2->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, mult2->GetMultiplicityVtx(histID)->ProjectionY("mymchist"));
2211
2212 TH1* result2 = mult2->GetMultiplicityESDCorrected(histID);
2213
2214 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2215 TFile* file = TFile::Open("EvaluateRegularizationEffect1.root");
2216 mult->LoadHistograms();
2217
2218 TH1* result1 = mult->GetMultiplicityESDCorrected(histID);
2219
2220 // compare results
2221 TCanvas* canvas = new TCanvas("EvaluateRegularizationEffect", "EvaluateRegularizationEffect", 1000, 800);
2222 canvas->Divide(2, 2);
2223
2224 canvas->cd(1);
2225 result1->SetLineColor(1);
2226 result1->DrawCopy();
2227 result2->SetLineColor(2);
2228 result2->DrawCopy("SAME");
2229 gPad->SetLogy();
2230
2231 result2->Rebin(2);
2232 result1->Scale(1.0 / result1->Integral());
2233 result2->Scale(1.0 / result2->Integral());
2234
2235 canvas->cd(2);
2236 result1->DrawCopy();
2237 result2->DrawCopy("SAME");
2238 gPad->SetLogy();
2239
2240 TH1* diff = (TH1*) result1->Clone("diff");
2241 diff->Add(result2, -1);
2242
2243 canvas->cd(3);
2244 diff->DrawCopy("HIST");
2245
2246 canvas->cd(4);
2247 diff->Divide(result1);
2248 diff->GetYaxis()->SetRangeUser(-0.3, 0.3);
2249 diff->DrawCopy("HIST");
2250
2251 Double_t chi2 = 0;
2252 for (Int_t i=1; i<=diff->GetNbinsX(); i++)
2253 chi2 += diff->GetBinContent(i) * diff->GetBinContent(i);
2254
2255 Printf("Chi2 is %e", chi2);
2256
2257 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2258 }
2259}