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