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