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