PWG0/dNdEta/runMultiplicitySelector.C

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