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