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