PWG0/dNdEta/runMultiplicitySelector.C

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