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