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