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