PWG0/multiplicity/correct.C

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