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