PWG0/multiplicity/correct.C

   1 /* $Id$ */
   2
   3 //
   4 // script to correct the multiplicity spectrum + helpers
   5 //
   6
   7 Bool_t is900GeV = 1;
   8 Bool_t is2360TeV = 0;
   9
  10 // 900 GeV, MC
  11 const Int_t kBinLimits[]   = { 42, 57, 60 };
  12 const Int_t kTrustLimits[] = { 26, 42, 54 };
  13
  14 // 2.36 TeV
  15 //const Int_t kBinLimits[]   = { 43, 67, 83 };
  16 //const Int_t kTrustLimits[] = { 33, 57, 73 };
  17
  18 // 7 TeV
  19 //const Int_t kBinLimits[]   = { 60, 120, 60 };
  20 //const Int_t kTrustLimits[] = { 26, 70, 54 };
  21
  22 void SetTPC()
  23 {
  24   gSystem->Load("libPWG0base");
  25   AliMultiplicityCorrection::SetQualityRegions(kFALSE);
  26 }
  27
  28 const char* GetMultLabel(Int_t etaR = -1, Bool_t trueM = kTRUE)
  29 {
  30         if (etaR == -1)
  31                 etaR = etaRange;
  32
  33         TString tmpStr((trueM) ? "True " : "Measured ");
  34
  35         if (etaR == 4)
  36         {
  37                 tmpStr += "multiplicity (full phase space)";
  38         }
  39         else
  40                 tmpStr += Form("multiplicity in |#eta| < %.1f", (etaR+1)* 0.5);
  41         return Form("%s", tmpStr.Data());
  42 }
  43
  44 void draw(const char* fileName = "multiplicity.root", const char* folder = "Multiplicity")
  45 {
  46   loadlibs();
  47
  48   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
  49
  50   TFile::Open(fileName);
  51   mult->LoadHistograms();
  52   mult->DrawHistograms();
  53
  54   TH2* hist = (TH2*) gROOT->FindObject("fCorrelation2_zy");
  55   canvas = new TCanvas("c1", "c1", 600, 500);
  56   canvas->SetTopMargin(0.05);
  57   hist->SetStats(kFALSE);
  58   hist->Draw("COLZ");
  59   hist->SetTitle(";true multiplicity in |#eta| < 1.5;measured multiplicity in |#eta| < 1.5");
  60   hist->GetYaxis()->SetTitleOffset(1.1);
  61   gPad->SetRightMargin(0.12);
  62   gPad->SetLogz();
  63
  64   canvas->SaveAs("responsematrix.eps");
  65 }
  66
  67 void loadlibs()
  68 {
  69   gSystem->Load("libTree");
  70   gSystem->Load("libVMC");
  71
  72   gSystem->Load("libSTEERBase");
  73   gSystem->Load("libESD.so");
  74   gSystem->Load("libAOD.so");
  75   gSystem->Load("libANALYSIS");
  76   gSystem->Load("libANALYSISalice");
  77   gSystem->Load("libPWG0base");
  78 }
  79
  80 void LoadAndInitialize(void* multVoid, void* esdVoid, void* multTriggerVoid, Int_t histID, Bool_t fullPhaseSpace, Int_t* geneLimits)
  81 {
  82   AliMultiplicityCorrection* mult = (AliMultiplicityCorrection*) multVoid;
  83   AliMultiplicityCorrection* esd = (AliMultiplicityCorrection*) esdVoid;
  84   AliMultiplicityCorrection* multTrigger = (AliMultiplicityCorrection*) multTriggerVoid;
  85
  86   for (Int_t i=0; i<3; i++)
  87     geneLimits[i] = kBinLimits[i];
  88
  89   // REBINNING
  90   if (1)
  91   {
  92     // 900 GeV
  93     if (is900GeV)
  94     {
  95       if (1)
  96       {
  97         Int_t bins05 = 31;
  98         Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
  99                   10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
 100                   20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 34.5, 38.5, 42.5, 50.5,
 101                   100.5 };
 102       }
 103
 104       if (0)
 105       {
 106         Int_t bins05 = 29;
 107         Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
 108                 10.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
 109                 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 34.5, 38.5, 42.5, 50.5,
 110                 100.5 };
 111       }
 112
 113       if (0)
 114       {
 115         Int_t bins05 = 25;
 116         Double_t* newBinsEta05 = new Double_t[bins05+1];
 117
 118         //newBinsEta05[0] = -0.5;
 119         //newBinsEta05[1] = 0.5;
 120         //newBinsEta05[2] = 1.5;
 121         //newBinsEta05[3] = 2.5;
 122
 123         for (Int_t i=0; i<bins05+1; i++)
 124           newBinsEta05[i] = -0.5 + i*2;
 125           //newBinsEta05[i] = 4.5 + (i-4)*2;
 126       }
 127
 128       Int_t bins10 = 54;
 129       Double_t newBinsEta10[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
 130                                   10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
 131                                   20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
 132                                   30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
 133                                   40.5, 42.5, 44.5, 46.5, 48.5, 50.5, 52.5, 54.5, 56.5, 58.5,
 134                                   60.5, 65.5, 70.5, 100.5 };
 135
 136       geneLimits[0] = bins05;
 137       geneLimits[1] = bins10;
 138       geneLimits[2] = bins10;
 139     }
 140
 141     // 2.36 TeV
 142     if (is2360TeV)
 143     {
 144       Int_t bins05 = 36;
 145       Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
 146                                   10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
 147                                   20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 32.5, 34.5, 36.5, 38.5,
 148                                   40.5, 45.5, 50.5, 55.5, 60.5, 100.5 };
 149       Int_t bins10 = 64;
 150       Double_t newBinsEta10[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
 151                                   10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
 152                                   20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
 153                                   30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
 154                                   40.5, 42.5, 44.5, 46.5, 48.5, 50.5, 52.5, 54.5, 56.5, 58.5,
 155                                   60.5, 62.5, 64.5, 66.5, 68.5, 70.5, 72.5, 74.5, 76.5, 78.5,
 156                                   80.5, 85.5, 90.5, 100.5 };
 157
 158       geneLimits[0] = bins05;
 159       geneLimits[1] = bins10;
 160       geneLimits[2] = bins10;
 161     }
 162
 163     // 7 TeV
 164     if (!is900GeV && !is2360TeV)
 165     {
 166       Int_t bins05 = 36;
 167       Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
 168                                   10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
 169                                   20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 32.5, 34.5, 36.5, 38.5,
 170                                   40.5, 45.5, 50.5, 55.5, 60.5, 100.5 };
 171
 172       Int_t bins10 = 85;
 173       Double_t newBinsEta10[] = {
 174           -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
 175                                   10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
 176                                   20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
 177                                   30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
 178                                   40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
 179                                   50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5,
 180                                   60.5, 62.5, 64.5, 66.5, 68.5, 70.5, 72.5, 74.5, 76.5, 78.5,
 181                                   80.5, 82.5, 84.5, 86.5, 88.5, 90.5, 92.5, 94.5, 96.5, 98.5,
 182                                   100.5, 105.5, 110.5, 115.5, 120.5 };
 183
 184       geneLimits[0] = bins05;
 185       geneLimits[1] = bins10;
 186       geneLimits[2] = bins10;
 187     }
 188
 189     esd->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
 190     mult->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
 191     multTrigger->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
 192   }
 193
 194   for (Int_t hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
 195   {
 196     TH2F* hist = esd->GetMultiplicityESD(hID);
 197
 198     mult->SetMultiplicityESD(hID, hist);
 199     mult->SetTriggeredEvents(hID, esd->GetTriggeredEvents(hID));
 200
 201     // insert trigger efficiency in flat response matrix
 202     for (Int_t evType = AliMultiplicityCorrection::kTrVtx; evType <= AliMultiplicityCorrection::kNSD; evType++)
 203       mult->SetMultiplicityMC(hID, evType, multTrigger->GetMultiplicityMC(hID, evType));
 204
 205     mult->SetNoVertexEvents(hID, multTrigger->GetNoVertexEvents(hID));
 206     mult->FixTriggerEfficiencies(10);
 207
 208     // small hack to get around charge conservation for full phase space ;-)
 209     if (fullPhaseSpace)
 210     {
 211       TH1* corr = mult->GetCorrelation(histID + 4);
 212
 213       for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
 214         for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
 215         {
 216           corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
 217           corr->SetBinError(i, j, corr->GetBinError(i-1, j));
 218         }
 219     }
 220   }
 221 }
 222
 223 void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityESD.root", Bool_t chi2 = 1, Int_t histID = 1, Bool_t fullPhaseSpace = kFALSE, Float_t beta = -25, Int_t eventType = 2 /* AliMultiplicityCorrection::kTrVtx */, const char* targetFile = "unfolded.root", const char* folderESD = "Multiplicity", Bool_t calcBias = kFALSE)
 224 {
 225   loadlibs();
 226
 227   Int_t geneLimits[] = { 0, 0, 0 };
 228
 229   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
 230   AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
 231   AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
 232
 233   LoadAndInitialize(mult, esd, multTrigger, histID, fullPhaseSpace, geneLimits);
 234
 235   //AliUnfolding::SetDebug(1);
 236
 237   for (Int_t hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
 238   {
 239     AliUnfolding::SetNbins(kBinLimits[hID], geneLimits[hID]);
 240
 241     TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
 242
 243     if (chi2)
 244     {
 245       if (is900GeV)
 246       {
 247         //Float_t regParamPol0[] = { 2, 2, 10 };   // TPCITS
 248         Float_t regParamPol0[] = { 5, 15, 25 }; // SPD
 249         Float_t regParamPol1[] =  { 0.15, 0.25, 0.5 };
 250       }
 251       else if (is2360TeV)
 252       {
 253         Float_t regParamPol0[] = { 10, 12, 40 };
 254         Float_t regParamPol1[] =  { 0.25, 0.25, 2 };
 255       }
 256       else
 257       {
 258         Float_t regParamPol0[] = { 1, 25, 10 };
 259         Float_t regParamPol1[] =  { 0.15, 0.5, 0.5 };
 260         AliUnfolding::SetSkipBinsBegin(3);
 261       }
 262
 263       Int_t reg = AliUnfolding::kPol0;
 264       if (beta > 0)
 265         reg = AliUnfolding::kPol1;
 266
 267       Float_t regParam = TMath::Abs(beta);
 268       if (histID == -1)
 269       {
 270         if (beta > 0)
 271           regParam = regParamPol1[hID];
 272         else
 273           regParam = regParamPol0[hID];
 274       }
 275       AliUnfolding::SetChi2Regularization(reg, regParam);
 276
 277       //AliUnfolding::SetChi2Regularization(AliUnfolding::kLog, 1000000);
 278       //AliUnfolding::SetChi2Regularization(AliUnfolding::kRatio, 10);
 279       //TVirtualFitter::SetDefaultFitter("Minuit2");
 280
 281       //AliUnfolding::SetCreateOverflowBin(kFALSE);
 282       //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); //mult->SetCreateBigBin(kFALSE);
 283       //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE);
 284      // AliUnfolding::SetChi2Regularization(AliUnfolding::kEntropy, beta);
 285       //mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5);
 286
 287       if (0)
 288       {
 289         // part for checking
 290         TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
 291         mcCompare->Sumw2();
 292         mult->ApplyMinuitFit(hID, fullPhaseSpace, AliMultiplicityCorrection::kMB, 0, kTRUE, mcCompare);
 293         mult->SetMultiplicityESDCorrected(hID, (TH1F*) mcCompare);
 294       }
 295       else
 296       {
 297         // Zero Bin
 298         Int_t zeroBin = 0;
 299         if (is900GeV) // from data
 300         {
 301           // background subtraction
 302           Int_t background = 0;
 303
 304           //background = 1091 + 4398; // MB1 for run 104824 - 52 (SPD)
 305           //background = 1087 + 4398; // MB1 for run 104824 - 52 (TPCITS)
 306
 307           background = 417 + 1758; // MB1 for run 104867 - 92 (SPD)
 308           //background = 1162+422;     // MB1 for run 104892 (TPCITS)
 309           //background = 5830 + 1384;    // MB1 for run 104824,25,45,52,67,92 (TPC runs) (TPCITS)
 310
 311           //background = 20 + 0; // V0AND for run 104824 - 52
 312           //background = 10 + 0; // V0AND for run 104867 - 92
 313
 314           Printf("NOTE: Subtracting %d background events", background);
 315           gSystem->Sleep(1000);
 316
 317           zeroBin = mult->GetTriggeredEvents(hID)->GetBinContent(1) - background - mult->GetMultiplicityESD(hID)->Integral(0, 2, 1, 1);
 318         }
 319         else if (is2360TeV)
 320         {
 321           // from MC
 322           Float_t fractionZeroBin = (multTrigger->GetTriggeredEvents(hID)->GetBinContent(1) - multTrigger->GetMultiplicityESD(hID)->Integral(0, 2, 1, 1)) / multTrigger->GetMultiplicityESD(hID)->Integral(1, 1, 1, mult->GetMultiplicityESD(hID)->GetNbinsY());
 323           zeroBin = fractionZeroBin * mult->GetMultiplicityESD(hID)->Integral(1, 1, 1, mult->GetMultiplicityESD(hID)->GetNbinsY());
 324
 325           Printf("Zero bin from MC: Estimating %d events with trigger but without vertex", zeroBin);
 326           gSystem->Sleep(1000);
 327         }
 328         else
 329         {
 330           AliUnfolding::SetSkip0BinInChi2(kTRUE);
 331         }
 332
 333         //mult->SetVertexRange(3, 4);
 334         mult->ApplyMinuitFit(hID, fullPhaseSpace, eventType, zeroBin, kFALSE, 0, calcBias);
 335       }
 336     }
 337     else
 338     {
 339       // HACK
 340       //mult->GetMultiplicityESD(hID)->SetBinContent(1, 1, 0);
 341       //for (Int_t bin=1; bin<=mult->GetCorrelation(hID)->GetNbinsY(); bin++)
 342       //  mult->GetCorrelation(hID)->SetBinContent(1, bin, 1, 0);
 343       AliUnfolding::SetChi2Regularization(AliUnfolding::kNone, 0);
 344       mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, beta, 0, kTRUE);
 345       //mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, 5, 0, kFALSE);
 346     }
 347   }
 348
 349   Printf("Writing result to %s", targetFile);
 350   TFile* file = TFile::Open(targetFile, "RECREATE");
 351   mult->SaveHistograms("Multiplicity");
 352   file->Write();
 353   file->Close();
 354
 355   if (histID == -1)
 356     return;
 357
 358   for (hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
 359   {
 360     TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
 361     TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
 362     mult->DrawComparison(Form("%s_%d_%f", (chi2) ? "MinuitChi2" : "Bayesian", hID, beta), hID, fullPhaseSpace, kTRUE, mcCompare, kFALSE, eventType);
 363
 364     Printf("<n> = %f", mult->GetMultiplicityESDCorrected(hID)->GetMean());
 365   }
 366 }
 367
 368 void correctAll(Int_t eventType)
 369 {
 370   const char* suffix = "";
 371   switch (eventType)
 372   {
 373     case 0: suffix = "trvtx"; break;
 374     case 1: suffix = "mb"; break;
 375     case 2: suffix = "inel"; break;
 376     case 3: suffix = "nsd"; break;
 377   }
 378
 379   correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, -1, 0, -1, eventType, TString(Form("chi2a_%s.root", suffix)));
 380   correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, -1, 0,  1, eventType, TString(Form("chi2b_%s.root", suffix)));
 381   correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 0, -1, 0, 40, eventType, TString(Form("bayes_%s.root", suffix)));
 382
 383   if (eventType == 3)
 384     drawAll(1);
 385   else if (eventType == 2)
 386     drawAllINEL();
 387 }
 388
 389 void drawAll(Bool_t showUA5 = kFALSE)
 390 {
 391   const char* files[] = { "chi2a_nsd.root", "chi2b_nsd.root", "bayes_nsd.root" };
 392   drawAll(files, showUA5);
 393 }
 394
 395 void drawAllINEL()
 396 {
 397   const char* files[] = { "chi2a_inel.root", "chi2b_inel.root", "bayes_inel.root" };
 398   drawAll(files);
 399 }
 400
 401 void drawAllMB()
 402 {
 403   const char* files[] = { "chi2a_mb.root", "chi2b_mb.root", "bayes_mb.root" };
 404   drawAll(files);
 405 }
 406
 407 void drawAll(const char** files, Bool_t showUA5 = kFALSE, Bool_t normalize = kFALSE)
 408 {
 409   loadlibs();
 410
 411   Int_t colors[] = { 1, 2, 4 };
 412
 413   c = new TCanvas(Form("c%d", gRandom->Uniform(100)), "c", 1800, 600);
 414   c->Divide(3, 1);
 415
 416   l = new TLegend(0.6, 0.6, 0.99, 0.99);
 417   l->SetFillColor(0);
 418
 419   TH1* hist0 = 0;
 420
 421   for (Int_t hID=0; hID<3; hID++)
 422   {
 423     c->cd(hID+1)->SetLogy();
 424     gPad->SetGridx();
 425     gPad->SetGridy();
 426     for (Int_t i=0; i<3; i++)
 427     {
 428       TFile::Open(files[i]);
 429
 430       hist = (TH1*) gFile->Get(Form("Multiplicity/fMultiplicityESDCorrected%d", hID));
 431
 432       Float_t average0 = hist->GetMean();
 433       hist2 = (TH1*) hist->Clone("temp");
 434       hist2->SetBinContent(1, 0);
 435       Float_t average1 = hist2->GetMean();
 436       Printf("%d: <N> = %.2f <N>/(eta) = %.2f | without 0: <N> = %.2f <N>/(eta) = %.2f", hID, average0, average0 / ((hID+1) - 0.4 * (hID / 2)), average1, average1 / ((hID+1) - 0.4 * (hID / 2)));
 437
 438       hist->SetLineColor(colors[i]);
 439
 440       hist->SetStats(0);
 441       hist->GetXaxis()->SetRangeUser(0, kBinLimits[hID]);
 442
 443       Float_t min = 0.1;
 444       Float_t max = hist->GetMaximum() * 1.5;
 445
 446       if (normalize)
 447         min = 1e-6;
 448
 449       hist->GetYaxis()->SetRangeUser(min, max);
 450       hist->SetTitle(";unfolded multiplicity;events");
 451
 452       if (hID == 0)
 453       {
 454         l->AddEntry(hist, files[i]);
 455       }
 456
 457       if (hist->Integral() <= 0)
 458         continue;
 459
 460       if (normalize)
 461         hist->Scale(1.0 / hist->Integral());
 462
 463       AliPWG0Helper::NormalizeToBinWidth(hist);
 464
 465       hist->DrawCopy((i == 0) ? "" : "SAME");
 466
 467       if (!hist0)
 468         hist0 = (TH1*) hist->Clone("hist0");
 469     }
 470
 471     if (hist0)
 472     {
 473       line = new TLine(kTrustLimits[hID], hist0->GetMinimum(), kTrustLimits[hID], hist0->GetMaximum());
 474       line->SetLineWidth(3);
 475       line->Draw();
 476     }
 477   }
 478
 479   c->cd(1);
 480   l->Draw();
 481
 482   if (showUA5)
 483   {
 484     TGraphErrors *gre = new TGraphErrors(23);
 485     gre->SetName("Graph");
 486     gre->SetTitle("UA5");
 487     gre->SetFillColor(1);
 488     gre->SetMarkerStyle(24);
 489     gre->SetPoint(0,0,0.15);
 490     gre->SetPointError(0,0.5,0.01);
 491     gre->SetPoint(1,1,0.171);
 492     gre->SetPointError(1,0.5,0.007);
 493     gre->SetPoint(2,2,0.153);
 494     gre->SetPointError(2,0.5,0.007);
 495     gre->SetPoint(3,3,0.124);
 496     gre->SetPointError(3,0.5,0.006);
 497     gre->SetPoint(4,4,0.099);
 498     gre->SetPointError(4,0.5,0.005);
 499     gre->SetPoint(5,5,0.076);
 500     gre->SetPointError(5,0.5,0.005);
 501     gre->SetPoint(6,6,0.057);
 502     gre->SetPointError(6,0.5,0.004);
 503     gre->SetPoint(7,7,0.043);
 504     gre->SetPointError(7,0.5,0.004);
 505     gre->SetPoint(8,8,0.032);
 506     gre->SetPointError(8,0.5,0.003);
 507     gre->SetPoint(9,9,0.024);
 508     gre->SetPointError(9,0.5,0.003);
 509     gre->SetPoint(10,10,0.018);
 510     gre->SetPointError(10,0.5,0.002);
 511     gre->SetPoint(11,11,0.013);
 512     gre->SetPointError(11,0.5,0.002);
 513     gre->SetPoint(12,12,0.01);
 514     gre->SetPointError(12,0.5,0.002);
 515     gre->SetPoint(13,13,0.007);
 516     gre->SetPointError(13,0.5,0.001);
 517     gre->SetPoint(14,14,0.005);
 518     gre->SetPointError(14,0.5,0.001);
 519     gre->SetPoint(15,15,0.004);
 520     gre->SetPointError(15,0.5,0.001);
 521     gre->SetPoint(16,16,0.0027);
 522     gre->SetPointError(16,0.5,0.0009);
 523     gre->SetPoint(17,17,0.0021);
 524     gre->SetPointError(17,0.5,0.0008);
 525     gre->SetPoint(18,18,0.0015);
 526     gre->SetPointError(18,0.5,0.0007);
 527     gre->SetPoint(19,19,0.0011);
 528     gre->SetPointError(19,0.5,0.0006);
 529     gre->SetPoint(20,20,0.0008);
 530     gre->SetPointError(20,0.5,0.0005);
 531     gre->SetPoint(21,22,0.00065);
 532     gre->SetPointError(21,1,0.0003);
 533     gre->SetPoint(22,27.5,0.00013);
 534     gre->SetPointError(22,3.5,7.1e-05);
 535
 536     for (Int_t i=0; i<gre->GetN(); i++)
 537     {
 538       gre->GetY()[i] *= hist0->Integral();
 539       gre->GetEY()[i] *= hist0->Integral();
 540     }
 541
 542     gre->Draw("P");
 543
 544     ratio = (TGraphErrors*) gre->Clone("ratio");
 545
 546     for (Int_t i=0; i<gre->GetN(); i++)
 547     {
 548       //Printf("%f %d", hist0->GetBinContent(hist0->FindBin(ratio->GetX()[i])), hist0->FindBin(ratio->GetX()[i]));
 549       Int_t bin = hist0->FindBin(gre->GetX()[i]);
 550       if (hist0->GetBinContent(bin) > 0)
 551       {
 552         ratio->GetEY()[i] = TMath::Sqrt((ratio->GetEY()[i] / ratio->GetY()[i])**2 + (hist0->GetBinError(bin) / hist0->GetBinContent(bin))**2);
 553         ratio->GetY()[i] /= hist0->GetBinContent(bin);
 554         ratio->GetEY()[i] *= ratio->GetY()[i];
 555       }
 556     }
 557     new TCanvas;
 558     gPad->SetGridx();
 559     gPad->SetGridy();
 560     //ratio->Print();
 561     ratio->Draw("AP");
 562     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
 563   }
 564
 565   c->SaveAs("draw_all.png");
 566 }
 567
 568 TH1* GetChi2Bias(Float_t alpha = -5)
 569 {
 570   loadlibs();
 571
 572   AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kMB;
 573   Int_t hID = 1;
 574
 575   correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, hID, 0, alpha, eventType, "nobias.root", "Multiplicity", kFALSE);
 576
 577   correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, hID, 0, alpha, eventType, "withbias.root", "Multiplicity", kTRUE);
 578
 579   // without bias calculation
 580   mult = AliMultiplicityCorrection::Open("nobias.root");
 581   baseold = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("baseold");
 582
 583   // with bias calculation
 584   mult = AliMultiplicityCorrection::Open("withbias.root");
 585   base = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("base");
 586
 587   // relative error plots
 588   ratio1 = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("hist1");
 589   ratio1->SetStats(0);
 590   ratio1->SetTitle(";unfolded multiplicity;relative error");
 591   ratio1->Reset();
 592   ratio2 = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("hist1");
 593   ratio2->Reset();
 594
 595   for (Int_t t = 0; t<baseold->GetNbinsX(); t++)
 596   {
 597     Printf("Bin %d; Content: %f; Chi2 Error: %f; Bias: %f; In percent: %.2f %.2f", t+1, base->GetBinContent(t+1), baseold->GetBinError(t+1), base->GetBinError(t+1), (base->GetBinContent(t+1) > 0) ? 100.0 * baseold->GetBinError(t+1) / base->GetBinContent(t+1) : -1, (base->GetBinContent(t+1) > 0) ? 100.0 * base->GetBinError(t+1) / base->GetBinContent(t+1) : -1);
 598     if (base->GetBinContent(t+1) <= 0)
 599       continue;
 600     ratio1->SetBinContent(t+1, baseold->GetBinError(t+1) / base->GetBinContent(t+1));
 601     ratio2->SetBinContent(t+1, base->GetBinError(t+1) / base->GetBinContent(t+1));
 602   }
 603
 604   //new TCanvas; base->Draw(); gPad->SetLogy();
 605
 606   c = new TCanvas;
 607   ratio1->GetYaxis()->SetRangeUser(0, 1);
 608   ratio1->GetXaxis()->SetRangeUser(0, 50);
 609   ratio1->Draw();
 610   ratio2->SetLineColor(2);
 611   ratio2->Draw("SAME");
 612
 613   return base;
 614 }
 615
 616 void CheckBayesianBias()
 617 {
 618   loadlibs();
 619
 620   const char* fileNameMC = "multiplicityMC.root";
 621   const char* folder = "Multiplicity";
 622   const char* fileNameESD = "multiplicityESD.root";
 623   const char* folderESD = "Multiplicity";
 624
 625   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
 626   AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
 627
 628   const Int_t kMaxBins = 35;
 629   AliUnfolding::SetNbins(kMaxBins, kMaxBins);
 630   //AliUnfolding::SetDebug(1);
 631
 632   Int_t hID = 0;
 633
 634   TH2F* hist = esd->GetMultiplicityESD(hID);
 635   mult->SetMultiplicityESD(hID, hist);
 636
 637   // without bias calculation
 638   mult->ApplyBayesianMethod(hID, kFALSE, 0, 1, 10, 0, 1);
 639   baseold = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("baseold");
 640
 641   // with bias calculation
 642   mult->ApplyBayesianMethod(hID, kFALSE, 0, 1, 10, 0, 2);
 643   base = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("base");
 644
 645   TH1* hist1 = new TH1F("hist1", "", 100, 0, 20);
 646   TH1* hist2 = new TH1F("hist2", "", 100, 0, 20);
 647
 648   for (Int_t t = 0; t<kMaxBins; t++)
 649   {
 650     Printf("Bin %d; Content: %f; Randomization Error: %f; Bias: %f; In percent: %.2f %.2f", t+1, base->GetBinContent(t+1), baseold->GetBinError(t+1), base->GetBinError(t+1), (base->GetBinContent(t+1) > 0) ? 100.0 * baseold->GetBinError(t+1) / base->GetBinContent(t+1) : -1, (base->GetBinContent(t+1) > 0) ? 100.0 * base->GetBinError(t+1) / base->GetBinContent(t+1) : -1);
 651
 652     hist1->Fill((base->GetBinContent(t+1) > 0) ? 100.0 * baseold->GetBinError(t+1) / base->GetBinContent(t+1) : -1);
 653     hist2->Fill((base->GetBinContent(t+1) > 0) ? 100.0 * base->GetBinError(t+1) / base->GetBinContent(t+1) : -1);
 654   }
 655
 656   new TCanvas;
 657   hist1->Draw();
 658   hist2->SetLineColor(2);
 659   hist2->Draw("SAME");
 660 }
 661
 662 void DrawUnfoldingLimit(Int_t hID, Float_t min, Float_t max)
 663 {
 664   line = new TLine(kTrustLimits[hID], min, kTrustLimits[hID], max);
 665   line->SetLineWidth(2);
 666   line->Draw();
 667 }
 668
 669 void DataScan(Int_t hID, Bool_t redo = kTRUE)
 670 {
 671   // makes a set of unfolded spectra and compares
 672   // don't forget FindUnfoldedLimit in plots.C
 673
 674   loadlibs();
 675
 676   // files...
 677   Bool_t mc = 1;
 678   const char* fileNameMC = "multiplicityMC.root";
 679   const char* folder = "Multiplicity";
 680   const char* fileNameESD = "multiplicityESD.root";
 681   const char* folderESD = "Multiplicity";
 682   const Int_t kMaxBins = kBinLimits[hID];
 683
 684   AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kTrVtx;
 685   Bool_t evaluteBias = kFALSE;
 686
 687   Int_t referenceCase = 2; // all results are shown as ratio to this case
 688
 689   // chi2 range
 690   AliUnfolding::RegularizationType regBegin = AliUnfolding::kPol0;
 691   AliUnfolding::RegularizationType regEnd = AliUnfolding::kPol1;
 692   Float_t regParamBegin[] = { 0, 1, 0.2, 3 };
 693   Float_t regParamEnd[] = { 0, 11, 0.5, 31 };
 694   Float_t regParamStep[] = { 0, 2, 2, TMath::Sqrt(10) };
 695
 696   // bayesian range
 697   Int_t iterBegin = 40;
 698   Int_t iterEnd = 41;
 699   Int_t iterStep = 20;
 700
 701   TList labels;
 702
 703   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
 704   AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
 705
 706   mult->SetMultiplicityESD(hID, esd->GetMultiplicityESD(hID));
 707
 708   // insert trigger efficiency in flat response matrix
 709   AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
 710   for (Int_t evType = AliMultiplicityCorrection::kTrVtx; evType <= AliMultiplicityCorrection::kNSD; evType++)
 711     mult->SetMultiplicityMC(hID, evType, multTrigger->GetMultiplicityMC(hID, evType));
 712
 713   mult->SetNoVertexEvents(hID, multTrigger->GetNoVertexEvents(hID));
 714   mult->FixTriggerEfficiencies(10);
 715
 716   Float_t fraction0Generated = multTrigger->GetFraction0Generated(hID);
 717
 718   if (mc)
 719     mult->SetMultiplicityMC(hID, eventType, esd->GetMultiplicityMC(hID, eventType));
 720
 721   AliUnfolding::SetNbins(kMaxBins, kMaxBins);
 722   //AliUnfolding::SetDebug(1);
 723
 724   Int_t count = -1;
 725
 726   for (AliUnfolding::RegularizationType reg = regBegin; reg <= regEnd; reg++)
 727   {
 728     for (Float_t regParam = regParamBegin[reg]; regParam < regParamEnd[reg]; regParam *= regParamStep[reg])
 729     {
 730       count++;
 731       labels.Add(new TObjString(Form("#chi^{2} Reg %d #beta = %.2f", (Int_t) reg, regParam)));
 732
 733       if (!redo)
 734         continue;
 735
 736       AliUnfolding::SetChi2Regularization(reg, regParam);
 737
 738       mult->ApplyMinuitFit(hID, kFALSE, eventType, fraction0Generated, kFALSE, 0, evaluteBias);
 739
 740       file = TFile::Open(Form("datascan_%d.root", count), "RECREATE");
 741       mult->SaveHistograms();
 742       file->Close();
 743     }
 744   }
 745
 746   for (Int_t iter = iterBegin; iter <= iterEnd; iter += iterStep)
 747   {
 748     count++;
 749     labels.Add(new TObjString(Form("Bayesian Iter = %d", iter)));
 750
 751     if (!redo)
 752       continue;
 753
 754     mult->ApplyBayesianMethod(hID, kFALSE, eventType, 1, iter, 0, kTRUE);
 755
 756     file = TFile::Open(Form("datascan_%d.root", count), "RECREATE");
 757     mult->SaveHistograms();
 758     file->Close();
 759   }
 760
 761   // 1) all distributions
 762   // 2) ratio to MC
 763   // 3) ratio to ref point
 764   // 4) relative error
 765   // 5) residuals
 766   c = new TCanvas("c", "c", 1200, 800);
 767   c->Divide(3, 2);
 768
 769   c->cd(1)->SetLogy();
 770
 771   // get reference point
 772   mult = AliMultiplicityCorrection::Open(Form("datascan_%d.root", referenceCase));
 773   if (!mult)
 774     return;
 775   refPoint = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("refPoint");
 776
 777   legend = new TLegend(0.5, 0.5, 0.99, 0.99);
 778   legend->SetFillColor(0);
 779
 780   TH1* residualSum = new TH1F("residualSum", ";;residual squared sum", count + 1, 0.5, count + 1.5);
 781
 782   Int_t allowedCount = count;
 783   count = 0;
 784   while (1)
 785   {
 786     mult = AliMultiplicityCorrection::Open(Form("datascan_%d.root", count++));
 787     if (!mult)
 788       break;
 789     if (count > allowedCount+1)
 790       break;
 791
 792     hist = (TH1*) mult->GetMultiplicityESDCorrected(hID);
 793     hist->SetLineColor((count-1) % 4 + 1);
 794     hist->SetLineStyle((count-1) / 4 + 1);
 795     hist->GetXaxis()->SetRangeUser(0, kMaxBins);
 796     hist->GetYaxis()->SetRangeUser(0.1, hist->GetMaximum() * 1.5);
 797     hist->SetStats(kFALSE);
 798     hist->SetTitle("");
 799
 800     legend->AddEntry(hist->Clone(), labels.At(count-1)->GetName());
 801
 802     c->cd(1);
 803     hist->DrawCopy((count == 1) ? "HIST" : "HISTSAME");
 804
 805     DrawUnfoldingLimit(hID, 0.1, hist->GetMaximum());
 806
 807     if (mc)
 808     {
 809       TH2* mcHist2d = mult->GetMultiplicityMC(hID, eventType);
 810       mcHist = mcHist2d->ProjectionY("mcmchist", 1, mcHist2d->GetNbinsX());
 811
 812       c->cd(1);
 813       mcHist->SetMarkerStyle(24);
 814       mcHist->Draw("P SAME");
 815
 816       c->cd(2);
 817       // calculate ratio using only the error on the mc bin
 818       ratio = (TH1*) hist->Clone("ratio");
 819       for (Int_t bin=1; bin<=ratio->GetNbinsX(); bin++)
 820       {
 821         if (mcHist->GetBinContent(bin) <= 0)
 822           continue;
 823         ratio->SetBinContent(bin, hist->GetBinContent(bin) / mcHist->GetBinContent(bin));
 824         ratio->SetBinError(bin, mcHist->GetBinError(bin) / mcHist->GetBinContent(bin) * ratio->GetBinContent(bin));
 825       }
 826
 827       ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
 828       ratio->GetYaxis()->SetTitle("Ratio unfolded / MC");
 829       ratio->DrawCopy((count == 1) ? "HIST" : "HISTSAME");
 830
 831       DrawUnfoldingLimit(hID, 0.5, 1.5);
 832     }
 833
 834     c->cd(3);
 835     // calculate ratio using no errors for now
 836     ratio = (TH1*) hist->Clone("ratio");
 837     for (Int_t bin=1; bin<=ratio->GetNbinsX(); bin++)
 838     {
 839       if (refPoint->GetBinContent(bin) <= 0)
 840         continue;
 841       ratio->SetBinContent(bin, hist->GetBinContent(bin) / refPoint->GetBinContent(bin));
 842       ratio->SetBinError(bin, 0);
 843     }
 844
 845     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
 846     ratio->GetYaxis()->SetTitle("Ratio unfolded / unfolded reference case");
 847     ratio->DrawCopy((count == 1) ? "" : "SAME");
 848
 849     DrawUnfoldingLimit(hID, 0.5, 1.5);
 850
 851     c->cd(4);
 852     ratio = (TH1*) hist->Clone("ratio");
 853     for (Int_t bin=1; bin<=ratio->GetNbinsX(); bin++)
 854     {
 855       if (ratio->GetBinContent(bin) <= 0)
 856         continue;
 857       ratio->SetBinContent(bin, ratio->GetBinError(bin) / ratio->GetBinContent(bin));
 858       ratio->SetBinError(bin, 0);
 859     }
 860     ratio->GetYaxis()->SetRangeUser(0, 1);
 861     ratio->GetYaxis()->SetTitle("Relative error");
 862     ratio->DrawCopy((count == 1) ? "" : "SAME");
 863
 864     DrawUnfoldingLimit(hID, 0, 1);
 865
 866     c->cd(5);
 867     Float_t sum;
 868     residuals = mult->GetResiduals(hID, AliMultiplicityCorrection::kTrVtx, sum);
 869     residuals->SetStats(0);
 870     residuals->GetXaxis()->SetRangeUser(0, kMaxBins);
 871     residuals->SetStats(kFALSE);
 872     residuals->SetLineColor((count-1) % 4 + 1);
 873     residuals->SetLineStyle((count-1) / 4 + 1);
 874     residuals->GetYaxis()->SetRangeUser(-5, 5);
 875     residuals->DrawCopy((count == 1) ? "" : "SAME");
 876
 877     residualSum->Fill(count, sum);
 878     residualSum->GetXaxis()->SetBinLabel(count, labels.At(count-1)->GetName());
 879   }
 880
 881   c->cd(6);
 882   residualSum->Draw();
 883   legend->Draw();
 884 }
 885
 886 void DrawBayesianIterations(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", Int_t histID = 1, Int_t eventType = 0 /* AliMultiplicityCorrection::kTrVtx */)
 887 {
 888   loadlibs();
 889
 890   const char* folder = "Multiplicity";
 891
 892   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
 893   TFile::Open(fileNameMC);
 894   mult->LoadHistograms();
 895
 896   AliMultiplicityCorrection* esd = new AliMultiplicityCorrection(folder, folder);
 897   TFile::Open(fileNameESD);
 898   esd->LoadHistograms();
 899
 900   TH2F* hist = esd->GetMultiplicityESD(histID);
 901   TH2F* hist2 = esd->GetMultiplicityMC(histID, eventType);
 902
 903   mult->SetMultiplicityESD(histID, hist);
 904   mult->SetMultiplicityMC(histID, eventType, hist2);
 905
 906   TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
 907   //mcCompare->Scale(1.0 / mcCompare->Integral());
 908
 909   TH1* esdHist = (TH1*) hist->ProjectionY("myesd", 1, 1)->Clone("myesd");
 910   //esdHist->Scale(1.0 / esdHist->Integral());
 911
 912   c = new TCanvas("c", "c", 1200, 600);
 913   c->Divide(2, 1);
 914
 915   c->cd(1);
 916   gPad->SetLeftMargin(0.12);
 917   gPad->SetTopMargin(0.05);
 918   gPad->SetRightMargin(0.05);
 919   gPad->SetLogy();
 920   gPad->SetGridx();
 921   gPad->SetGridy();
 922
 923   mcCompare->GetXaxis()->SetRangeUser(0, 80);
 924   mcCompare->SetStats(0);
 925   mcCompare->SetFillColor(5);
 926   mcCompare->SetLineColor(5);
 927   mcCompare->SetTitle(Form(";%s;Entries", GetMultLabel(1)));
 928   mcCompare->GetYaxis()->SetRangeUser(2, esdHist->GetMaximum() * 2);
 929   mcCompare->GetYaxis()->SetTitleOffset(1.3);
 930   mcCompare->Draw("HIST");
 931   esdHist->SetMarkerStyle(5);
 932   esdHist->Draw("P HIST SAME");
 933
 934   c->cd(2);
 935   gPad->SetTopMargin(0.05);
 936   gPad->SetRightMargin(0.05);
 937   gPad->SetGridx();
 938   gPad->SetGridy();
 939
 940   trueMeasuredRatio = (TH1*) mcCompare->Clone("trueMeasuredRatio");
 941   trueMeasuredRatio->Divide(esdHist);
 942   trueMeasuredRatio->SetStats(0);
 943   trueMeasuredRatio->SetTitle(Form(";%s;Ratio", GetMultLabel(1)));
 944   trueMeasuredRatio->GetYaxis()->SetTitleOffset(1.3);
 945   trueMeasuredRatio->GetYaxis()->SetRangeUser(0, 2);
 946   // ROOT displays all values > 2 at 2 which looks weird
 947   for (Int_t i=1; i<=trueMeasuredRatio->GetNbinsX(); i++)
 948     if (trueMeasuredRatio->GetBinContent(i) > 2)
 949       trueMeasuredRatio->SetBinContent(i, 0);
 950   trueMeasuredRatio->SetMarkerStyle(5);
 951   trueMeasuredRatio->Draw("P");
 952
 953   Int_t colors[] = {1, 2, 4, 6};
 954
 955   legend = new TLegend(0.15, 0.13, 0.5, 0.5);
 956   legend->SetFillColor(0);
 957   legend->SetTextSize(0.04);
 958
 959   legend->AddEntry(mcCompare, "True", "F");
 960   legend->AddEntry(esdHist, "Measured", "P");
 961
 962   legend2 = new TLegend(0.15, 0.13, 0.5, 0.4);
 963   legend2->SetFillColor(0);
 964   legend2->SetTextSize(0.04);
 965
 966   legend2->AddEntry(trueMeasuredRatio, "Measured", "P");
 967
 968   Int_t iters[] = {1, 3, 10, -1};
 969   for (Int_t i = 0; i<4; i++)
 970   {
 971     Int_t iter = iters[i];
 972     mult->ApplyBayesianMethod(histID, kFALSE, eventType, 1, iter, 0, 0);
 973     corr = mult->GetMultiplicityESDCorrected(histID);
 974     corr->Scale(1.0 / corr->Integral());
 975     corr->Scale(mcCompare->Integral());
 976     corr->GetXaxis()->SetRangeUser(0, 80);
 977     //corr->SetMarkerStyle(20+iter);
 978     corr->SetLineColor(colors[i]);
 979     corr->SetLineStyle(i+1);
 980     corr->SetLineWidth(2);
 981
 982     c->cd(1);
 983     corr->DrawCopy("SAME HIST");
 984
 985     c->cd(2);
 986     clone = (TH1*) corr->Clone("clone");
 987     clone->Divide(mcCompare, corr);
 988     clone->GetYaxis()->SetRangeUser(0, 2);
 989     clone->DrawCopy("SAME HIST");
 990
 991     legend->AddEntry((TH1*) corr->Clone(), (iter > -1) ? Form("%d iteration%s", iter, (iter == 1) ? "" : "s") : "Convergence", "L");
 992     legend2->AddEntry((TH1*) corr->Clone(), (iter > -1) ? Form("%d iteration%s", iter, (iter == 1) ? "" : "s") : "Convergence", "L");
 993   }
 994
 995   c->cd(1);
 996   legend->Draw();
 997
 998   c->cd(2);
 999   legend2->Draw();
1000
1001   c->SaveAs("bayesian_iterations.eps");
1002 }
1003
1004 void CompareChi2Bayesian(Int_t histID = 2, const char* chi2File = "chi2.root", const char* bayesianFile = "bayesian.root", const char* label1 = "Chi2", const char* label2 = "Bayesian", const char* mcFile = 0, Float_t simpleCorrect = 0)
1005 {
1006   const char* folder = "Multiplicity";
1007
1008   loadlibs();
1009
1010   AliMultiplicityCorrection* chi2 = new AliMultiplicityCorrection(folder, folder);
1011   TFile::Open(chi2File);
1012   chi2->LoadHistograms();
1013
1014   AliMultiplicityCorrection* bayesian = new AliMultiplicityCorrection(folder, folder);
1015   TFile::Open(bayesianFile);
1016   bayesian->LoadHistograms();
1017
1018   histRAW = chi2->GetMultiplicityESD(histID)->ProjectionY("raw", 1, chi2->GetMultiplicityESD(histID)->GetNbinsX());
1019   histRAW->Scale(1.0 / histRAW->Integral());
1020
1021   histC = chi2->GetMultiplicityESDCorrected(histID);
1022   histB = bayesian->GetMultiplicityESDCorrected(histID);
1023
1024   c = new TCanvas("CompareChi2Bayesian", "CompareChi2Bayesian", 800, 600);
1025   c->SetRightMargin(0.05);
1026   c->SetTopMargin(0.05);
1027   c->SetLogy();
1028   c->SetGridx();
1029   c->SetGridy();
1030
1031   histC->SetTitle(";N;P(N)");
1032   histC->SetStats(kFALSE);
1033   histC->GetXaxis()->SetRangeUser(0, 100);
1034
1035   histC->SetLineColor(1);
1036   histB->SetLineColor(2);
1037   histRAW->SetLineColor(3);
1038
1039   histC->DrawCopy("HISTE");
1040   histB->DrawCopy("HISTE SAME");
1041   histRAW->DrawCopy("SAME");
1042
1043   legend = new TLegend(0.2, 0.2, 0.4, 0.4);
1044   legend->SetFillColor(0);
1045
1046   legend->AddEntry(histC, label1);
1047   legend->AddEntry(histB, label2);
1048   legend->AddEntry(histRAW, "raw ESD");
1049
1050   histGraph = 0;
1051   if (simpleCorrect > 0)
1052   {
1053     graph = new TGraph;
1054     graph->SetMarkerStyle(25);
1055     graph->SetFillColor(0);
1056     for (Int_t bin=1; bin<=histRAW->GetNbinsX(); bin++)
1057       graph->SetPoint(graph->GetN(), histRAW->GetXaxis()->GetBinCenter(bin) * simpleCorrect, histRAW->GetBinContent(bin));
1058
1059     graph->Draw("PSAME");
1060     legend->AddEntry(graph, "weighting");
1061
1062     // now create histogram from graph and normalize
1063     histGraph = (TH1*) histRAW->Clone();
1064     histGraph->Reset();
1065     for (Int_t bin=1; bin<=histGraph->GetNbinsX(); bin++)
1066     {
1067       Int_t j=1;
1068       for (j=1; j<graph->GetN(); j++)
1069         if (graph->GetX()[j] > histGraph->GetXaxis()->GetBinCenter(bin))
1070           break;
1071       if (j == graph->GetN())
1072         continue;
1073       if (histGraph->GetXaxis()->GetBinCenter(bin) - graph->GetX()[j] < graph->GetX()[j-1] - histGraph->GetXaxis()->GetBinCenter(bin))
1074         j--;
1075       histGraph->SetBinContent(bin, graph->GetY()[j]);
1076     }
1077
1078     Printf("Integral = %f", histGraph->Integral());
1079     histGraph->Scale(1.0 / histGraph->Integral());
1080
1081     histGraph->SetLineColor(6);
1082     histGraph->DrawCopy("SAME");
1083     legend->AddEntry(histGraph, "weighting normalized");
1084   }
1085
1086   if (mcFile)
1087   {
1088     AliMultiplicityCorrection* mc = new AliMultiplicityCorrection(folder, folder);
1089     TFile::Open(mcFile);
1090     mc->LoadHistograms();
1091
1092     histMC = mc->GetMultiplicityVtx(histID)->ProjectionY("mc", 1, mc->GetMultiplicityVtx(histID)->GetNbinsX());
1093     histMC->Sumw2();
1094     histMC->Scale(1.0 / histMC->Integral());
1095
1096     histMC->Draw("HISTE SAME");
1097     histMC->SetLineColor(4);
1098     legend->AddEntry(histMC, "MC");
1099   }
1100
1101   legend->Draw();
1102
1103   c->SaveAs(Form("%s.png", c->GetName()));
1104   c->SaveAs(Form("%s.eps", c->GetName()));
1105
1106   if (!mcFile)
1107     return;
1108
1109   // build ratios
1110
1111   c = new TCanvas("CompareChi2BayesianRatio", "CompareChi2BayesianRatio", 800, 600);
1112   c->SetRightMargin(0.05);
1113   c->SetTopMargin(0.05);
1114   c->SetGridx();
1115   c->SetGridy();
1116
1117   for (Int_t bin=1; bin<=histC->GetNbinsX(); bin++)
1118   {
1119     if (histMC->GetBinContent(bin) > 0)
1120     {
1121       histC->SetBinContent(bin, histC->GetBinContent(bin) / histMC->GetBinContent(bin));
1122       histB->SetBinContent(bin, histB->GetBinContent(bin) / histMC->GetBinContent(bin));
1123
1124       /*
1125       if (histGraph)
1126       {
1127         histGraph->SetBinContent(bin, histGraph->GetBinContent(bin) / histMC->GetBinContent(bin));
1128         histGraph->SetBinError(bin, 0);
1129       }
1130       */
1131
1132       // TODO errors?
1133       histC->SetBinError(bin, 0);
1134       histB->SetBinError(bin, 0);
1135     }
1136   }
1137
1138   histC->GetYaxis()->SetRangeUser(0.5, 2);
1139   histC->GetYaxis()->SetTitle("Unfolded / MC");
1140
1141   histC->Draw("HIST");
1142   histB->Draw("HIST SAME");
1143   /*
1144   if (histGraph)
1145     histGraph->Draw("HIST SAME");
1146   */
1147
1148   /*
1149   if (simpleCorrect > 0)
1150   {
1151     graph2 = new TGraph;
1152     graph2->SetMarkerStyle(25);
1153     graph2->SetFillColor(0);
1154     for (Int_t i=0; i<graph->GetN(); i++)
1155     {
1156       Float_t mcValue = histMC->GetBinContent(histMC->FindBin(graph->GetX()[i]));
1157       Float_t mcError = histMC->GetBinError(histMC->FindBin(graph->GetX()[i]));
1158       if (mcValue > 0)
1159         graph2->SetPoint(graph2->GetN(), graph->GetX()[i], graph->GetY()[i] / mcValue);
1160     }
1161     graph2->Draw("PSAME");
1162   }
1163   */
1164
1165   c->SaveAs(Form("%s.png", c->GetName()));
1166   c->SaveAs(Form("%s.eps", c->GetName()));
1167 }
1168
1169
1170 void CompareMC(Int_t histID, Int_t eventType, const char* file1, const char* file2, const char* label1, const char* label2)
1171 {
1172   const char* folder = "Multiplicity";
1173
1174   loadlibs();
1175
1176   AliMultiplicityCorrection* mc1 = new AliMultiplicityCorrection(folder, folder);
1177   TFile::Open(file1);
1178   mc1->LoadHistograms();
1179   histMC1 = mc1->GetMultiplicityMC(histID, eventType)->ProjectionY("mc1", 1, mc1->GetMultiplicityMC(histID, eventType)->GetNbinsX());
1180   histMC1->Scale(1.0 / histMC1->Integral());
1181
1182   AliMultiplicityCorrection* mc2 = new AliMultiplicityCorrection(folder, folder);
1183   TFile::Open(file2);
1184   mc2->LoadHistograms();
1185   histMC2 = mc2->GetMultiplicityMC(histID, eventType)->ProjectionY("mc2", 1, mc2->GetMultiplicityMC(histID, eventType)->GetNbinsX());
1186   histMC2->Scale(1.0 / histMC2->Integral());
1187
1188   c = new TCanvas("CompareMC", "CompareMC", 800, 600);
1189   c->SetRightMargin(0.05);
1190   c->SetTopMargin(0.05);
1191   c->SetLogy();
1192
1193   histMC1->SetTitle(";N;P(N)");
1194   histMC1->SetStats(kFALSE);
1195   histMC1->GetXaxis()->SetRangeUser(0, 100);
1196
1197   histMC1->SetLineColor(1);
1198   histMC2->SetLineColor(2);
1199
1200   histMC1->Draw();
1201   histMC2->Draw("SAME");
1202
1203   legend = new TLegend(0.2, 0.2, 0.4, 0.4);
1204   legend->SetFillColor(0);
1205
1206   legend->AddEntry(histMC1, label1);
1207   legend->AddEntry(histMC2, label2);
1208
1209   legend->Draw();
1210
1211   c->SaveAs(Form("%s.gif", c->GetName()));
1212   c->SaveAs(Form("%s.eps", c->GetName()));
1213 }
1214
1215 void* fit2Step(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE)
1216 {
1217   gSystem->Load("libPWG0base");
1218
1219   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1220
1221   TFile::Open(fileNameMC);
1222   mult->LoadHistograms("Multiplicity");
1223
1224   TFile::Open(fileNameESD);
1225   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
1226   TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID)));
1227   //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID));
1228
1229   mult->SetMultiplicityESD(histID, hist);
1230
1231   // small hack to get around charge conservation for full phase space ;-)
1232   if (fullPhaseSpace)
1233   {
1234     TH1* corr = mult->GetCorrelation(histID + 4);
1235
1236     for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
1237       for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
1238       {
1239         corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
1240         corr->SetBinError(i, j, corr->GetBinError(i-1, j));
1241       }
1242   }
1243
1244   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1245   mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
1246   mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
1247
1248   TH1* result = (TH1*) mult->GetMultiplicityESDCorrected((fullPhaseSpace) ? 4 : histID))->Clone("firstresult");
1249
1250   mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 100000);
1251   mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE, result);
1252   mult->DrawComparison("MinuitChi2_Step2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
1253
1254   return mult;
1255 }
1256
1257 const char* GetRegName(Int_t type)
1258 {
1259   switch (type)
1260   {
1261     case AliMultiplicityCorrection::kNone:      return "None"; break;
1262     case AliMultiplicityCorrection::kPol0:      return "Pol0"; break;
1263     case AliMultiplicityCorrection::kPol1:      return "Pol1"; break;
1264     case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break;
1265     case AliMultiplicityCorrection::kEntropy:   return "Reduced cross-entropy"; break;
1266     case AliMultiplicityCorrection::kLog   :    return "Log"; break;
1267   }
1268   return 0;
1269 }
1270
1271 const char* GetEventTypeName(Int_t type)
1272 {
1273   switch (type)
1274   {
1275     case AliMultiplicityCorrection::kTrVtx:   return "trigger, vertex"; break;
1276     case AliMultiplicityCorrection::kMB:      return "minimum bias"; break;
1277     case AliMultiplicityCorrection::kINEL:    return "inelastic"; break;
1278   }
1279   return 0;
1280 }
1281
1282 void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", const char* targetDir = "bayesian", Int_t histID = 1)
1283 {
1284   gSystem->mkdir(targetDir);
1285
1286   gSystem->Load("libPWG0base");
1287
1288   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1289   TFile::Open(fileNameMC);
1290   mult->LoadHistograms("Multiplicity");
1291
1292   AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
1293   TFile::Open(fileNameESD);
1294   multESD->LoadHistograms("Multiplicity");
1295   mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
1296
1297   Int_t count = 0; // just to order the saved images...
1298
1299   TFile* graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir), "RECREATE");
1300
1301   Int_t colors[] =  {1, 2, 3, 4, 6};
1302   Int_t markers[] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3, 4, 5, 6};
1303
1304   for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
1305   //for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
1306   {
1307     TString tmp;
1308     tmp.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
1309
1310     TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600);
1311
1312     for (Int_t i = 1; i <= 5; i++)
1313     {
1314       Int_t iterArray[5] = {2, 5, 10, 20, -1};
1315       Int_t iter = iterArray[i-1];
1316
1317       TGraph* fitResultsMC[3];
1318       for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1319       {
1320         fitResultsMC[region] = new TGraph;
1321         fitResultsMC[region]->SetTitle(Form("%d iterations - reg. %d", iter, region+1));
1322         fitResultsMC[region]->GetXaxis()->SetTitle("smoothing parameter #alpha");
1323         fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region));
1324         fitResultsMC[region]->SetName(Form("%s_MC_%d", tmp.Data(), i * AliMultiplicityCorrection::kQualityRegions + region - 2));
1325         fitResultsMC[region]->SetFillColor(0);
1326         //fitResultsMC[region]->SetMarkerStyle(markers[(i-1) * AliMultiplicityCorrection::kQualityRegions + region]);
1327         fitResultsMC[region]->SetMarkerStyle(markers[i-1]);
1328         fitResultsMC[region]->SetLineColor(colors[i-1]);
1329         fitResultsMC[region]->SetMarkerColor(colors[i-1]);
1330       }
1331
1332       TGraph* fitResultsRes = new TGraph;
1333       fitResultsRes->SetTitle(Form("%d iterations", iter));
1334       fitResultsRes->SetName(Form("%s_Res_%d", tmp.Data(), i));
1335       fitResultsRes->GetXaxis()->SetTitle("smoothing parameter");
1336       fitResultsRes->GetYaxis()->SetTitle("P_{2}");
1337
1338       fitResultsRes->SetFillColor(0);
1339       fitResultsRes->SetMarkerStyle(markers[i-1]);
1340       fitResultsRes->SetMarkerColor(colors[i-1]);
1341       fitResultsRes->SetLineColor(colors[i-1]);
1342
1343       for (Float_t weight = 0.0; weight < 1.01; weight += 0.2)
1344       {
1345         Float_t chi2MC = 0;
1346         Float_t residuals = 0;
1347
1348         mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter, 0, kFALSE);
1349         mult->DrawComparison(Form("%s/BayesianIterSmooth_%03d_%d_%d_%f", targetDir, count++, type, iter, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
1350         mult->GetComparisonResults(&chi2MC, 0, &residuals);
1351
1352         for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1353           fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region));
1354
1355         fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
1356       }
1357
1358       graphFile->cd();
1359       for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1360         fitResultsMC[region]->Write();
1361
1362       fitResultsRes->Write();
1363     }
1364   }
1365
1366   graphFile->Close();
1367 }
1368
1369 void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = kTRUE)
1370 {
1371   gSystem->Load("libPWG0base");
1372
1373   TString name;
1374   TFile* graphFile = 0;
1375   if (type == -1)
1376   {
1377     name = "EvaluateChi2Method";
1378     graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir));
1379   }
1380   else
1381   {
1382     name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
1383     graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir));
1384   }
1385
1386   TCanvas* canvas = new TCanvas(name, name, 800, 500);
1387   if (type == -1)
1388   {
1389     canvas->SetLogx();
1390     canvas->SetLogy();
1391   }
1392   canvas->SetTopMargin(0.05);
1393   canvas->SetGridx();
1394   canvas->SetGridy();
1395
1396   TLegend* legend = new TLegend(0.8, 0.15, 0.98, 0.98);
1397   legend->SetFillColor(0);
1398
1399   Int_t count = 1;
1400
1401   Float_t xMin = 1e20;
1402   Float_t xMax = 0;
1403
1404   Float_t yMin = 1e20;
1405   Float_t yMax = 0;
1406
1407   Float_t yMinRegion[3];
1408   for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
1409     yMinRegion[i] = 1e20;
1410
1411   TString xaxis, yaxis;
1412
1413   while (1)
1414   {
1415     TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
1416     TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
1417
1418     if (!mc)
1419       break;
1420
1421     xaxis = mc->GetXaxis()->GetTitle();
1422     yaxis = mc->GetYaxis()->GetTitle();
1423
1424     mc->Print();
1425
1426     if (res)
1427       res->Print();
1428
1429     xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
1430     yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin());
1431
1432     xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
1433     yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
1434
1435     if (plotRes && res)
1436     {
1437       xMin = TMath::Min(xMin, res->GetXaxis()->GetXmin());
1438       yMin = TMath::Min(yMin, res->GetYaxis()->GetXmin());
1439
1440       xMax = TMath::Max(xMax, res->GetXaxis()->GetXmax());
1441       yMax = TMath::Max(yMax, res->GetYaxis()->GetXmax());
1442     }
1443
1444     for (Int_t i=0; i<mc->GetN(); ++i)
1445       yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]);
1446
1447     count++;
1448   }
1449
1450   for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
1451     Printf("Minimum for region %d is %f", i, yMinRegion[i]);
1452
1453   if (type >= 0)
1454   {
1455     xaxis = "smoothing parameter";
1456   }
1457   else if (type == -1)
1458   {
1459     xaxis = "weight parameter";
1460     xMax *= 5;
1461   }
1462   //yaxis = "P_{1} (2 <= t <= 150)";
1463
1464   printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
1465
1466   TGraph* dummy = new TGraph;
1467   dummy->SetPoint(0, xMin, yMin);
1468   dummy->SetPoint(1, xMax, yMax);
1469   dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data()));
1470
1471   dummy->SetMarkerColor(0);
1472   dummy->Draw("AP");
1473   dummy->GetYaxis()->SetMoreLogLabels(1);
1474
1475   count = 1;
1476
1477   while (1)
1478   {
1479     TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
1480     TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
1481
1482     //printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res);
1483
1484     if (!mc)
1485       break;
1486
1487     printf("Loaded %d sucessful.\n", count);
1488
1489     if (type == -1)
1490     {
1491       legend->AddEntry(mc, Form("Eq. (%d) - reg. %d", 10 + (count-1) / 3, 1+ (count-1) % 3));
1492     }
1493     else
1494       legend->AddEntry(mc);
1495
1496     mc->Draw("SAME PC");
1497
1498     if (plotRes && res)
1499     {
1500       legend->AddEntry(res);
1501       res->Draw("SAME PC");
1502     }
1503
1504     count++;
1505   }
1506
1507   legend->Draw();
1508
1509   canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
1510   canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
1511 }
1512
1513 void EvaluateDrawResultRegions(const char* targetDir, Int_t type = 0)
1514 {
1515   loadlibs();
1516
1517   TString name;
1518   TFile* graphFile = 0;
1519   if (type == -1)
1520   {
1521     name = "EvaluateChi2Method";
1522     graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir));
1523   }
1524   else
1525   {
1526     name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
1527     graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir));
1528   }
1529
1530   TCanvas* canvas = new TCanvas(name, name, 1200, 800);
1531   //canvas->Divide(1, AliMultiplicityCorrection::kQualityRegions, 0, 0);
1532   canvas->SetTopMargin(0);
1533   canvas->SetBottomMargin(0);
1534   canvas->SetRightMargin(0.05);
1535   canvas->Range(0, 0, 1, 1);
1536
1537   TPad* pad[4];
1538   pad[3] = new TPad(Form("%s_pad1", name.Data()), "", 0, 0, 0.5, 0.5); pad[3]->SetTopMargin(0); pad[3]->SetBottomMargin(0.15);
1539   pad[1] = new TPad(Form("%s_pad2", name.Data()), "", 0, 0.5, 0.5, 1); pad[1]->SetBottomMargin(0);
1540   pad[0] = new TPad(Form("%s_pad3", name.Data()), "", 0.5, 0, 1, 0.5); pad[0]->SetTopMargin(0); pad[0]->SetBottomMargin(0.15);
1541   pad[2] = new TPad(Form("%s_pad4", name.Data()), "", 0.5, 0.5, 1, 1); pad[2]->SetBottomMargin(0);
1542
1543   Float_t yMinRegion[3];
1544   for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
1545     yMinRegion[i] = 1e20;
1546
1547   for (Int_t region = 0; region <= AliMultiplicityCorrection::kQualityRegions; region++)
1548   {
1549     canvas->cd();
1550     pad[region]->Draw();
1551     pad[region]->cd();
1552     pad[region]->SetRightMargin(0.05);
1553
1554     if (type == -1)
1555     {
1556       pad[region]->SetLogx();
1557     }
1558     pad[region]->SetLogy();
1559
1560     pad[region]->SetGridx();
1561     pad[region]->SetGridy();
1562
1563     TLegend* legend = new TLegend(0.5, 0.4, 0.85, 0.85);
1564     legend->SetFillColor(0);
1565     //legend->SetNColumns(3);
1566     legend->SetTextSize(0.06);
1567     Int_t count = 0;
1568
1569     Float_t xMin = 1e20;
1570     Float_t xMax = 0;
1571
1572     Float_t yMin = 1e20;
1573     Float_t yMax = 0;
1574
1575     TString xaxis, yaxis;
1576
1577     while (1)
1578     {
1579       count++;
1580
1581       if (region > 0)
1582       {
1583         TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
1584         if (!mc)
1585           break;
1586
1587         if (TString(mc->GetTitle()).Contains(Form("reg. %d", region)) == kFALSE)
1588           continue;
1589
1590         for (Int_t i=0; i<mc->GetN(); ++i)
1591           yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]);
1592       }
1593       else
1594       {
1595         TGraph* mc = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
1596         if (!mc)
1597           break;
1598       }
1599
1600       xaxis = mc->GetXaxis()->GetTitle();
1601       yaxis = mc->GetYaxis()->GetTitle();
1602
1603       //mc->Print();
1604
1605       xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
1606       yMin = TMath::Min(yMin, TMath::MinElement(mc->GetN(), mc->GetY()));
1607
1608       xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
1609       yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
1610
1611       //Printf("%f %f %f %f", xMin, xMax, yMin, yMax);
1612     }
1613
1614     if (type >= 0)
1615     {
1616       xaxis = "Smoothing parameter #alpha";
1617     }
1618     else if (type == -1)
1619     {
1620       xaxis = "Weight parameter #beta";
1621       //xMax *= 5;
1622     }
1623     if (region > 0)
1624     {
1625       yaxis = Form("Q_{1} in region %d ", region); // (2 <= t <= 150)";
1626     }
1627     else
1628       yaxis = "Q_{2} ";
1629
1630     printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
1631
1632     if (region > 0)
1633     {
1634       if (type == -1)
1635       {
1636         yMin = 0.534622;
1637         yMax = 19.228941;
1638       }
1639       else
1640       {
1641         yMin = 0.759109;
1642         yMax = 10;
1643       }
1644     }
1645
1646     if (type != -1)
1647       xMax = 1.0;
1648
1649     TH1* dummy = new TH2F("dummy", "dummy", 100, xMin, xMax, 100, yMin * 0.9, yMax * 1.1);
1650     dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data()));
1651
1652     if (region == 0 && type != -1)
1653       dummy->GetYaxis()->SetMoreLogLabels(1);
1654     dummy->SetLabelSize(0.06, "xy");
1655     dummy->SetTitleSize(0.06, "xy");
1656     dummy->GetXaxis()->SetTitleOffset(1.2);
1657     dummy->GetYaxis()->SetTitleOffset(0.8);
1658     dummy->SetStats(0);
1659
1660     dummy->DrawCopy();
1661
1662
1663     count = 0;
1664
1665     while (1)
1666     {
1667       count++;
1668
1669       if (region > 0)
1670       {
1671         TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
1672         if (mc && TString(mc->GetTitle()).Contains(Form("reg. %d", region)) == kFALSE)
1673           continue;
1674       }
1675       else
1676         TGraph* mc = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
1677
1678       if (!mc)
1679         break;
1680       //printf("Loaded %d sucessful.\n", count);
1681
1682       if (type == -1)
1683       {
1684         //legend->AddEntry(mc, Form("Eq. (%d) - reg. %d", 10 + (count-1) / 3, 1+ (count-1) % 3));
1685         //legend->AddEntry(mc, Form("Eq. (7.%d)", 11 + (count-1) / 3));
1686         const char* names[] = { "Const", "Linear", "Log" };
1687         legend->AddEntry(mc, names[(count-1) / 3], "PL");
1688       }
1689       else
1690       {
1691         TString label(mc->GetTitle());
1692         TString newLabel(label(0, label.Index(" ")));
1693         //Printf("%s", newLabel.Data());
1694         if (newLabel.Atoi() == -1)
1695         {
1696           newLabel = "Convergence";
1697         }
1698         else
1699           newLabel = label(0, label.Index("iterations") + 10);
1700
1701         legend->AddEntry(mc, newLabel, "PL");
1702       }
1703
1704       mc->Draw("SAME PC");
1705
1706     }
1707
1708     if (region == 2)
1709       legend->Draw();
1710   }
1711
1712   for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
1713     Printf("Minimum for region %d is %f", i, yMinRegion[i]);
1714
1715   canvas->Modified();
1716   canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
1717   canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
1718 }
1719
1720 void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", const char* targetDir = "chi2compare", Int_t histID = 1)
1721 {
1722   gSystem->mkdir(targetDir);
1723
1724   gSystem->Load("libPWG0base");
1725
1726   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1727
1728   TFile::Open(fileNameMC);
1729   mult->LoadHistograms("Multiplicity");
1730
1731   TFile::Open(fileNameESD);
1732   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
1733   TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
1734
1735   mult->SetMultiplicityESD(histID, hist);
1736
1737   Int_t count = 0; // just to order the saved images...
1738   Int_t colors[] = {1, 2, 4, 3};
1739   Int_t markers[] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3};
1740
1741   TGraph* fitResultsRes = 0;
1742
1743   TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE");
1744
1745   //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kLog; ++type)
1746   //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kEntropy; type <= AliMultiplicityCorrection::kEntropy; ++type)
1747   for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
1748   {
1749     TGraph* fitResultsMC[3];
1750     for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1751     {
1752       fitResultsMC[region] = new TGraph;
1753       fitResultsMC[region]->SetTitle(Form("Eq. (%d) - reg. %d", type+9, region+1));
1754       fitResultsMC[region]->GetXaxis()->SetTitle("weight parameter #alpha");
1755       fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region));
1756       fitResultsMC[region]->SetName(Form("EvaluateChi2Method_MC_%d", type * AliMultiplicityCorrection::kQualityRegions + region - 2));
1757       fitResultsMC[region]->SetFillColor(0);
1758       //fitResultsMC[region]->SetMarkerStyle(markers[(type-1) * AliMultiplicityCorrection::kQualityRegions + region]);
1759       //fitResultsMC[region]->SetLineColor(colors[region]);
1760       fitResultsMC[region]->SetMarkerStyle(markers[type-1]);
1761       fitResultsMC[region]->SetMarkerColor(colors[type-1]);
1762       fitResultsMC[region]->SetLineColor(colors[type-1]);
1763     }
1764
1765     fitResultsRes = new TGraph;
1766     fitResultsRes->SetTitle(Form("%s residual chi2", GetRegName(type)));
1767     fitResultsRes->SetName(Form("EvaluateChi2Method_Res_%d", type));
1768     fitResultsRes->GetXaxis()->SetTitle("Weight Parameter");
1769
1770     fitResultsRes->SetFillColor(0);
1771     fitResultsRes->SetMarkerStyle(markers[type-1]);
1772     fitResultsRes->SetMarkerColor(colors[type-1]);
1773     fitResultsRes->SetLineColor(colors[type-1]);
1774
1775     for (Int_t i=0; i<15; ++i)
1776     {
1777       Float_t weight = TMath::Power(TMath::Sqrt(10), i+1);
1778       //Float_t weight = TMath::Power(10, i+2);
1779
1780       //if (type == AliMultiplicityCorrection::kEntropy) weight = 1e4 * (i+1) * 1.5;
1781
1782       Float_t chi2MC = 0;
1783       Float_t residuals = 0;
1784       Float_t chi2Limit = 0;
1785
1786       TString runName;
1787       runName.Form("MinuitChi2_%02d_%d_%f", count++, type, weight);
1788
1789       mult->SetRegularizationParameters(type, weight);
1790       Int_t status = mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
1791       mult->DrawComparison(Form("%s/%s", targetDir, runName.Data()), histID, kFALSE, kTRUE, hist2->ProjectionY("mymc", 1, hist2->GetNbinsX()));
1792       if (status != 0)
1793       {
1794         printf("MINUIT did not succeed. Skipping...\n");
1795         continue;
1796       }
1797
1798       mult->GetComparisonResults(&chi2MC, 0, &residuals);
1799       TH1* result = mult->GetMultiplicityESDCorrected(histID);
1800       result->SetName(runName);
1801       result->Write();
1802
1803       for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1804         fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region));
1805
1806       fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
1807     }
1808
1809     graphFile->cd();
1810     for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1811       fitResultsMC[region]->Write();
1812     fitResultsRes->Write();
1813   }
1814
1815   graphFile->Close();
1816 }
1817
1818 void EvaluateChi2MethodAll()
1819 {
1820   EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M");
1821   EvaluateChi2Method("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M");
1822   EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD");
1823   EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M");
1824   EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
1825 }
1826
1827 void EvaluateBayesianMethodAll()
1828 {
1829   EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M");
1830   EvaluateBayesianMethod("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M");
1831   EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD");
1832   EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M");
1833   EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
1834 }
1835
1836 void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
1837 {
1838   gSystem->mkdir(targetDir);
1839
1840   gSystem->Load("libPWG0base");
1841
1842   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1843
1844   TFile::Open(fileNameMC);
1845   mult->LoadHistograms("Multiplicity");
1846
1847   TFile::Open(fileNameESD);
1848   AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
1849   multESD->LoadHistograms("Multiplicity");
1850
1851   mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
1852
1853   TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 1200);
1854   canvas->Divide(3, 3);
1855
1856   Int_t count = 0;
1857
1858   for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
1859   {
1860     TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY("mymc");
1861     mc->Sumw2();
1862
1863     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1864     mult->ApplyMinuitFit(histID, kFALSE, type);
1865     mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc);
1866     TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
1867
1868     mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1);
1869     mult->DrawComparison(Form("%s/CompareMethods_%d_Bayesian", targetDir, type), histID, kFALSE, kTRUE, mc);
1870     TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
1871
1872     mc->GetXaxis()->SetRangeUser(0, 150);
1873     chi2Result->GetXaxis()->SetRangeUser(0, 150);
1874
1875 /*    // skip errors for now
1876     for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
1877     {
1878       chi2Result->SetBinError(i, 0);
1879       bayesResult->SetBinError(i, 0);
1880     }*/
1881
1882     canvas->cd(++count);
1883     mc->SetFillColor(kYellow);
1884     mc->DrawCopy();
1885     chi2Result->SetLineColor(kRed);
1886     chi2Result->DrawCopy("SAME");
1887     bayesResult->SetLineColor(kBlue);
1888     bayesResult->DrawCopy("SAME");
1889     gPad->SetLogy();
1890
1891     canvas->cd(count + 3);
1892     chi2ResultRatio = (TH1*) chi2Result->Clone("chi2ResultRatio");
1893     bayesResultRatio = (TH1*) bayesResult->Clone("bayesResultRatio");
1894     chi2ResultRatio->Divide(chi2Result, mc, 1, 1, "");
1895     bayesResultRatio->Divide(bayesResult, mc, 1, 1, "");
1896
1897     chi2ResultRatio->GetYaxis()->SetRangeUser(0.5, 1.5);
1898
1899     chi2ResultRatio->DrawCopy("HIST");
1900     bayesResultRatio->DrawCopy("SAME HIST");
1901
1902     canvas->cd(count + 6);
1903     chi2Result->Divide(chi2Result, bayesResult, 1, 1, "");
1904     chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
1905     chi2Result->DrawCopy("HIST");
1906   }
1907
1908   canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
1909   canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
1910 }
1911
1912 void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
1913 {
1914   gSystem->Load("libPWG0base");
1915
1916   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1917
1918   TFile::Open(fileNameMC);
1919   mult->LoadHistograms("Multiplicity");
1920
1921   const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_200k.root", "multiplicityMC_400k.root", "multiplicityMC_600k.root", "multiplicityMC_800k.root" };
1922
1923   TGraph* fitResultsChi2[3];
1924   TGraph* fitResultsBayes[3];
1925
1926   for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1927   {
1928     fitResultsChi2[region] = new TGraph;
1929     fitResultsChi2[region]->SetTitle(";Nevents;Chi2");
1930     fitResultsChi2[region]->SetName(Form("fitResultsChi2_%d", region));
1931     fitResultsChi2[region]->SetMarkerStyle(20+region);
1932
1933     fitResultsBayes[region] = new TGraph;
1934     fitResultsBayes[region]->SetTitle(";Nevents;Chi2");
1935     fitResultsBayes[region]->SetName(Form("fitResultsBayes_%d", region));
1936     fitResultsBayes[region]->SetMarkerStyle(20+region);
1937     fitResultsBayes[region]->SetMarkerColor(2);
1938   }
1939
1940   TGraph* fitResultsChi2Limit = new TGraph;  fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
1941   TGraph* fitResultsBayesLimit = new TGraph; fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
1942   TGraph* fitResultsChi2Res = new TGraph;       fitResultsChi2Res->SetTitle(";Nevents;Chi2");
1943   TGraph* fitResultsBayesRes = new TGraph;      fitResultsBayesRes->SetTitle(";Nevents;Chi2");
1944
1945   fitResultsChi2Limit->SetName("fitResultsChi2Limit");
1946   fitResultsBayesLimit->SetName("fitResultsBayesLimit");
1947   fitResultsChi2Res->SetName("fitResultsChi2Res");
1948   fitResultsBayesRes->SetName("fitResultsBayesRes");
1949
1950   TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600);
1951   canvas->Divide(5, 2);
1952
1953   Float_t min = 1e10;
1954   Float_t max = 0;
1955
1956   TFile* file = TFile::Open("StatisticsPlot.root", "RECREATE");
1957   file->Close();
1958
1959   for (Int_t i=0; i<5; ++i)
1960   {
1961     TFile::Open(files[i]);
1962     AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
1963     multESD->LoadHistograms("Multiplicity");
1964
1965     mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
1966     Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries();
1967     TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i));
1968
1969     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1970     mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
1971     mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
1972
1973     Int_t chi2MCLimit = 0;
1974     Float_t chi2Residuals = 0;
1975     mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals);
1976     for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1977     {
1978       fitResultsChi2[region]->SetPoint(fitResultsChi2[region]->GetN(), nEntries, mult->GetQuality(region));
1979       min = TMath::Min(min, mult->GetQuality(region));
1980       max = TMath::Max(max, mult->GetQuality(region));
1981     }
1982     fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit);
1983     fitResultsChi2Res->SetPoint(fitResultsChi2Res->GetN(), nEntries, chi2Residuals);
1984
1985     TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
1986
1987     mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0, kFALSE);
1988     mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
1989     TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
1990     mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals);
1991     for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1992     {
1993       fitResultsBayes[region]->SetPoint(fitResultsBayes[region]->GetN(), nEntries, mult->GetQuality(region));
1994       min = TMath::Min(min, mult->GetQuality(region));
1995       max = TMath::Max(max, mult->GetQuality(region));
1996     }
1997     fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit);
1998     fitResultsBayesRes->SetPoint(fitResultsBayesRes->GetN(), nEntries, chi2Residuals);
1999
2000     mc->GetXaxis()->SetRangeUser(0, 150);
2001     chi2Result->GetXaxis()->SetRangeUser(0, 150);
2002
2003     // skip errors for now
2004     for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
2005     {
2006       chi2Result->SetBinError(j, 0);
2007       bayesResult->SetBinError(j, 0);
2008     }
2009
2010     canvas->cd(i+1);
2011     mc->SetFillColor(kYellow);
2012     mc->DrawCopy();
2013     chi2Result->SetLineColor(kRed);
2014     chi2Result->DrawCopy("SAME");
2015     bayesResult->SetLineColor(kBlue);
2016     bayesResult->DrawCopy("SAME");
2017     gPad->SetLogy();
2018
2019     canvas->cd(i+6);
2020     chi2Result->Divide(chi2Result, mc, 1, 1, "B");
2021     bayesResult->Divide(bayesResult, mc, 1, 1, "B");
2022
2023     // skip errors for now
2024     for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
2025     {
2026       chi2Result->SetBinError(j, 0);
2027       bayesResult->SetBinError(j, 0);
2028     }
2029
2030     chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
2031     chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
2032
2033     chi2Result->DrawCopy("");
2034     bayesResult->DrawCopy("SAME");
2035
2036     TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
2037     mc->Write();
2038     chi2Result->Write();
2039     bayesResult->Write();
2040     file->Close();
2041   }
2042
2043   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
2044   canvas->SaveAs(Form("%s.C", canvas->GetName()));
2045
2046   TCanvas* canvas2 = new TCanvas("StatisticsPlot2", "StatisticsPlot2", 800, 400);
2047   canvas2->Divide(2, 1);
2048
2049   canvas2->cd(1);
2050
2051   for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
2052   {
2053     fitResultsChi2[region]->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
2054     fitResultsChi2[region]->Draw(((region == 0) ? "AP" : "P SAME"));
2055
2056     fitResultsBayes[region]->Draw("P SAME");
2057   }
2058
2059   gPad->SetLogy();
2060
2061   canvas2->cd(2);
2062   fitResultsChi2Limit->SetMarkerStyle(20);
2063   fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
2064   fitResultsChi2Limit->Draw("AP");
2065
2066   fitResultsBayesLimit->SetMarkerStyle(3);
2067   fitResultsBayesLimit->SetMarkerColor(2);
2068   fitResultsBayesLimit->Draw("P SAME");
2069
2070   canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
2071   canvas2->SaveAs(Form("%s.C", canvas2->GetName()));
2072
2073   TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
2074
2075   for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
2076   {
2077     fitResultsChi2[region]->Write();
2078     fitResultsBayes[region]->Write();
2079   }
2080   fitResultsChi2Limit->Write();
2081   fitResultsBayesLimit->Write();
2082   fitResultsChi2Res->Write();
2083   fitResultsBayesRes->Write();
2084   file->Close();
2085 }
2086
2087 void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 1)
2088 {
2089   loadlibs();
2090
2091   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2092
2093   TFile::Open(fileNameMC);
2094   mult->LoadHistograms("Multiplicity");
2095
2096   const char* files[] = { "multiplicityESD.root", "multiplicityESD_100k_1.root", "multiplicityESD_100k_2.root", "multiplicityESD_100k_1.root", "multiplicityESD_100k_2.root" }
2097
2098   // this one we try to unfold
2099   TFile::Open(files[0]);
2100   AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
2101   multESD->LoadHistograms("Multiplicity");
2102   mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
2103   TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY("mc", 1, 1);
2104
2105   TGraph* fitResultsChi2 = new TGraph;
2106   fitResultsChi2->SetTitle(";Input Dist ID;Chi2");
2107   TGraph* fitResultsBayes = new TGraph;
2108   fitResultsBayes->SetTitle(";Input Dist ID;Chi2");
2109   TGraph* fitResultsChi2Limit = new TGraph;
2110   fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach");
2111   TGraph* fitResultsBayesLimit = new TGraph;
2112   fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach");
2113
2114   TCanvas* canvas = new TCanvas("StartingConditions", "StartingConditions", 1200, 600);
2115   canvas->Divide(8, 2);
2116
2117   TCanvas* canvas3 = new TCanvas("StartingConditions3", "StartingConditions3", 1000, 400);
2118   canvas3->Divide(2, 1);
2119
2120   Float_t min = 1e10;
2121   Float_t max = 0;
2122
2123   TH1* firstChi = 0;
2124   TH1* firstBayesian = 0;
2125   TH1* startCond = multESD->GetMultiplicityESD(histID)->ProjectionY("startCond");
2126
2127   TLegend* legend = new TLegend(0.7, 0.7, 1, 1);
2128
2129   TFile* file = TFile::Open("StartingConditions.root", "RECREATE");
2130   mc->Write();
2131   file->Close();
2132
2133   for (Int_t i=0; i<8; ++i)
2134   {
2135     if (i == 0)
2136     {
2137       startCond = (TH1*) mc->Clone("startCond2");
2138     }
2139     else if (i < 6)
2140     {
2141       TFile::Open(files[i-1]);
2142       AliMultiplicityCorrection* multESD2 = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2");
2143       multESD2->LoadHistograms("Multiplicity");
2144       startCond = multESD2->GetMultiplicityESD(histID)->ProjectionY("startCond");
2145     }
2146     else if (i == 6)
2147     {
2148       func = new TF1("nbd", "exp(log([0]) + TMath::LnGamma([2]+x) - TMath::LnGamma([2]) - TMath::LnGamma(x+1) + log([1] / ([1]+[2])) * x + log(1.0 + [1]/[2]) * -[2])", 0, 50);
2149
2150       func->SetParNames("scaling", "averagen", "k");
2151       func->SetParLimits(0, 1e-3, 1e10);
2152       func->SetParLimits(1, 0.001, 1000);
2153       func->SetParLimits(2, 0.001, 1000);
2154
2155       func->SetParameters(1, 10, 2);
2156       for (Int_t j=2; j<=startCond->GetNbinsX(); j++)
2157         startCond->SetBinContent(j, func->Eval(j-1));
2158     }
2159     else if (i == 7)
2160     {
2161       for (Int_t j=1; j<=startCond->GetNbinsX(); j++)
2162         startCond->SetBinContent(j, 1);
2163     }
2164
2165     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 1e5);
2166     mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond);
2167     //mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
2168
2169     Float_t chi2MC = 0;
2170     Int_t chi2MCLimit = 0;
2171     mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
2172     fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC);
2173     fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit);
2174     min = TMath::Min(min, chi2MC);
2175     max = TMath::Max(max, chi2MC);
2176
2177     TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
2178     if (!firstChi)
2179       firstChi = (TH1*) chi2Result->Clone("firstChi");
2180
2181     mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 10, startCond);
2182     //mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
2183     TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
2184     if (!firstBayesian)
2185       firstBayesian = (TH1*) bayesResult->Clone("firstBayesian");
2186
2187     mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
2188     fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC);
2189     fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit);
2190
2191     TFile* file = TFile::Open("StartingConditions.root", "UPDATE");
2192     chi2Result->Write();
2193     bayesResult->Write();
2194     file->Close();
2195
2196     min = TMath::Min(min, chi2MC);
2197     max = TMath::Max(max, chi2MC);
2198     mc->GetXaxis()->SetRangeUser(0, 150);
2199     chi2Result->GetXaxis()->SetRangeUser(0, 150);
2200
2201     // skip errors for now
2202     for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
2203     {
2204       chi2Result->SetBinError(j, 0);
2205       bayesResult->SetBinError(j, 0);
2206     }
2207
2208     canvas3->cd(1);
2209     TH1* tmp = (TH1*) chi2Result->Clone("tmp");
2210     tmp->SetTitle("Difference to best initial conditions;Npart;Ratio");
2211     tmp->Divide(firstChi);
2212     tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
2213     tmp->GetXaxis()->SetRangeUser(0, 200);
2214     tmp->SetLineColor(i+1);
2215     legend->AddEntry(tmp, Form("%d", i));
2216     tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
2217
2218     canvas3->cd(2);
2219     tmp = (TH1*) bayesResult->Clone("tmp");
2220     tmp->SetTitle("Difference to best initial conditions;Npart;Ratio");
2221     tmp->Divide(firstBayesian);
2222     tmp->SetLineColor(i+1);
2223     tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
2224     tmp->GetXaxis()->SetRangeUser(0, 200);
2225     tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
2226
2227     canvas->cd(i+1);
2228     mc->SetFillColor(kYellow);
2229     mc->DrawCopy();
2230     chi2Result->SetLineColor(kRed);
2231     chi2Result->DrawCopy("SAME");
2232     bayesResult->SetLineColor(kBlue);
2233     bayesResult->DrawCopy("SAME");
2234     gPad->SetLogy();
2235
2236     canvas->cd(i+9);
2237     chi2Result->Divide(chi2Result, mc, 1, 1, "B");
2238     bayesResult->Divide(bayesResult, mc, 1, 1, "B");
2239
2240     // skip errors for now
2241     for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
2242     {
2243       chi2Result->SetBinError(j, 0);
2244       bayesResult->SetBinError(j, 0);
2245     }
2246
2247     chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
2248     chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
2249
2250     chi2Result->DrawCopy("");
2251     bayesResult->DrawCopy("SAME");
2252   }
2253
2254   canvas3->cd(1);
2255   legend->Draw();
2256
2257   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
2258
2259   TCanvas* canvas2 = new TCanvas("StartingConditions2", "StartingConditions2", 800, 400);
2260   canvas2->Divide(2, 1);
2261
2262   canvas2->cd(1);
2263   fitResultsChi2->SetMarkerStyle(20);
2264   fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
2265   fitResultsChi2->Draw("AP");
2266
2267   fitResultsBayes->SetMarkerStyle(3);
2268   fitResultsBayes->SetMarkerColor(2);
2269   fitResultsBayes->Draw("P SAME");
2270
2271   gPad->SetLogy();
2272
2273   canvas2->cd(2);
2274   fitResultsChi2Limit->SetMarkerStyle(20);
2275   fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
2276   fitResultsChi2Limit->Draw("AP");
2277
2278   fitResultsBayesLimit->SetMarkerStyle(3);
2279   fitResultsBayesLimit->SetMarkerColor(2);
2280   fitResultsBayesLimit->Draw("P SAME");
2281
2282   canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
2283   canvas3->SaveAs(Form("%s.gif", canvas3->GetName()));
2284 }
2285
2286 void Merge(Int_t n, const char** files, const char* output)
2287 {
2288   // 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" };
2289
2290
2291   gSystem->Load("libPWG0base");
2292
2293   AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n];
2294   TList list;
2295   for (Int_t i=0; i<n; ++i)
2296   {
2297     TString name("Multiplicity");
2298     if (i > 0)
2299       name.Form("Multiplicity%d", i);
2300
2301     TFile::Open(files[i]);
2302     Printf("Loading from file %s", files[i]);
2303     data[i] = new AliMultiplicityCorrection(name, name);
2304     data[i]->LoadHistograms("Multiplicity");
2305     if (i > 0)
2306       list.Add(data[i]);
2307   }
2308
2309   data[0]->Merge(&list);
2310
2311   //data[0]->DrawHistograms();
2312
2313   TFile::Open(output, "RECREATE");
2314   data[0]->SaveHistograms();
2315   gFile->Close();
2316 }
2317
2318 void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
2319 {
2320   loadlibs();
2321
2322   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2323
2324   TFile::Open(fileName);
2325   mult->LoadHistograms("Multiplicity");
2326
2327   TF1* func = 0;
2328
2329   if (caseNo >= 4)
2330   {
2331     TF1* func = new TF1("nbd", "exp(log([0]) + TMath::LnGamma([2]+x) - TMath::LnGamma([2]) - TMath::LnGamma(x+1) + log([1] / ([1]+[2])) * x + log(1.0 + [1]/[2]) * -[2])", 0, 200);
2332     //func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 200);
2333     func->SetParNames("scaling", "averagen", "k");
2334   }
2335
2336   switch (caseNo)
2337   {
2338     case 0: func = new TF1("flat", "1000"); break;
2339     case 1: func = new TF1("flat", "501-x"); break;
2340     case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
2341     case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
2342     case 4: func->SetParameters(2e5, 15, 2); break;
2343     case 5: func->SetParameters(1, 13, 7); break;
2344     case 6: func->SetParameters(1e7, 30, 4); break;
2345     case 7: func->SetParameters(1e7, 30, 2); break; // ***
2346     case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break;
2347
2348     default: return;
2349   }
2350
2351   new TCanvas;
2352   func->Draw();
2353
2354   mult->SetGenMeasFromFunc(func, 1);
2355
2356   TFile::Open("out.root", "RECREATE");
2357   mult->SaveHistograms();
2358
2359   new TCanvas; mult->GetMultiplicityESD(1)->ProjectionY()->DrawCopy();
2360   new TCanvas; mult->GetMultiplicityVtx(1)->ProjectionY("mc", 1, mult->GetMultiplicityVtx(1)->GetNbinsX())->DrawCopy();
2361
2362   //mult->ApplyBayesianMethod(2, kFALSE);
2363   //mult->ApplyMinuitFit(2, kFALSE);
2364   //mult->ApplyGaussianMethod(2, kFALSE);
2365   //mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx);
2366 }
2367
2368 void smoothCorrelationMap(const char* fileName = "multiplicity.root", Int_t corrMatrix = 1)
2369 {
2370   gSystem->Load("libPWG0base");
2371
2372   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2373
2374   TFile::Open(fileName);
2375   mult->LoadHistograms("Multiplicity");
2376
2377   // empty under/overflow bins in x, otherwise Project3D takes them into account
2378   TH3* corr = mult->GetCorrelation(corrMatrix);
2379   for (Int_t j=1; j<=corr->GetYaxis()->GetNbins(); ++j)
2380   {
2381     for (Int_t k=1; k<=corr->GetZaxis()->GetNbins(); ++k)
2382     {
2383       corr->SetBinContent(0, j, k, 0);
2384       corr->SetBinContent(corr->GetXaxis()->GetNbins()+1, j, k, 0);
2385     }
2386   }
2387
2388   TH2* proj = (TH2*) corr->Project3D("zy");
2389
2390   // normalize correction for given nPart
2391   for (Int_t i=1; i<=proj->GetNbinsX(); ++i)
2392   {
2393     Double_t sum = proj->Integral(i, i, 1, proj->GetNbinsY());
2394     if (sum <= 0)
2395       continue;
2396
2397     for (Int_t j=1; j<=proj->GetNbinsY(); ++j)
2398     {
2399       // npart sum to 1
2400       proj->SetBinContent(i, j, proj->GetBinContent(i, j) / sum);
2401       proj->SetBinError(i, j, proj->GetBinError(i, j) / sum);
2402     }
2403   }
2404
2405   new TCanvas;
2406   proj->DrawCopy("COLZ");
2407
2408   TH1* scaling = proj->ProjectionY("scaling", 1, 1);
2409   scaling->Reset();
2410   scaling->SetMarkerStyle(3);
2411   //scaling->GetXaxis()->SetRangeUser(0, 50);
2412   TH1* mean = (TH1F*) scaling->Clone("mean");
2413   TH1* width = (TH1F*) scaling->Clone("width");
2414
2415   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
2416   lognormal->SetParNames("scaling", "mean", "sigma");
2417   lognormal->SetParameters(1, 1, 1);
2418   lognormal->SetParLimits(0, 1, 1);
2419   lognormal->SetParLimits(1, 0, 100);
2420   lognormal->SetParLimits(2, 1e-3, 1);
2421
2422   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);
2423   nbd->SetParNames("scaling", "averagen", "k");
2424   nbd->SetParameters(1, 13, 5);
2425   nbd->SetParLimits(0, 1, 1);
2426   nbd->SetParLimits(1, 1, 100);
2427   nbd->SetParLimits(2, 1, 1e8);
2428
2429   TF1* poisson = new TF1("poisson", "[0] * exp(-(x+[2])) * (x+[2])**[1] / TMath::Factorial([1])", 0.01, 50);
2430   poisson->SetParNames("scaling", "k", "deltax");
2431   poisson->SetParameters(1, 1, 0);
2432   poisson->SetParLimits(0, 0, 10);
2433   poisson->SetParLimits(1, 0.01, 100);
2434   poisson->SetParLimits(2, 0, 10);
2435
2436   TF1* mygaus = new TF1("mygaus", "[0] * exp(-(x-[1])**2 / 2 / [2] - [3] * log(x + [4]) / [5])", 0.01, 50);
2437   mygaus->SetParNames("scaling", "mean", "width", "scale2log", "logmean", "logwidth");
2438   mygaus->SetParameters(1, 0, 1, 1, 0, 1);
2439   mygaus->SetParLimits(2, 1e-5, 10);
2440   mygaus->SetParLimits(4, 1, 1);
2441   mygaus->SetParLimits(5, 1e-5, 10);
2442
2443   //TF1* sqrt = new TF1("sqrt", "[0] + [1] * sqrt((x + [3]) * [2])", 0, 50);
2444   TF1* sqrt = new TF1("sqrt", "[0] + (x + [1])**[2]", 0, 50);
2445   sqrt->SetParNames("ydelta", "exp", "xdelta");
2446   sqrt->SetParameters(0, 0, 1);
2447   sqrt->SetParLimits(1, 0, 10);
2448
2449   const char* fitWith = "gaus";
2450
2451   for (Int_t i=1; i<=80; ++i)
2452   {
2453     printf("Fitting %d...\n", i);
2454
2455     TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e");
2456
2457     //hist->GetXaxis()->SetRangeUser(0, 50);
2458     //lognormal->SetParameter(0, hist->GetMaximum());
2459     hist->Fit(fitWith, "0 M", "");
2460
2461     TF1* func = hist->GetListOfFunctions()->FindObject(fitWith);
2462
2463     if (0 && (i % 5 == 0))
2464     {
2465       pad = new TCanvas;
2466       hist->Draw();
2467       func->Clone()->Draw("SAME");
2468       pad->SetLogy();
2469     }
2470
2471     scaling->SetBinContent(i, func->GetParameter(0));
2472     mean->SetBinContent(i, func->GetParameter(1));
2473     width->SetBinContent(i, func->GetParameter(2));
2474
2475     scaling->SetBinError(i, func->GetParError(0));
2476     mean->SetBinError(i, func->GetParError(1));
2477     width->SetBinError(i, func->GetParError(2));
2478   }
2479
2480   TF1* log = new TF1("log", "[0] + [1] * log([2] * x)", 0.01, 500);
2481   log->SetParameters(0, 1, 1);
2482   log->SetParLimits(1, 0, 100);
2483   log->SetParLimits(2, 1e-3, 10);
2484
2485   TF1* over = new TF1("over", "[0] + [1] / (x+[2])", 0.01, 500);
2486   over->SetParameters(0, 1, 0);
2487   //over->SetParLimits(0, 0, 100);
2488   over->SetParLimits(1, 1e-3, 10);
2489   over->SetParLimits(2, 0, 100);
2490
2491   c1 = new TCanvas("fitparams", "fitparams", 1200, 400);
2492   c1->Divide(3, 1);
2493
2494   c1->cd(1);
2495   scaling->Draw("P");
2496
2497   //TF1* scalingFit = new TF1("mypol0", "[0]");
2498   TF1* scalingFit = over;
2499   scaling->Fit(scalingFit, "", "", 3, 140);
2500   scalingFit->SetRange(0, 200);
2501   scalingFit->Draw("SAME");
2502
2503   c1->cd(2);
2504   mean->Draw("P");
2505
2506   //TF1* meanFit = log;
2507   TF1* meanFit = new TF1("mypol1", "[0]+[1]*x");
2508   mean->Fit(meanFit, "", "", 3, 140);
2509   meanFit->SetRange(0, 200);
2510   meanFit->Draw("SAME");
2511
2512   c1->cd(3);
2513   width->Draw("P");
2514
2515   //TF1* widthFit = over;
2516   TF1* widthFit = new TF1("mypol", "[0]+[1]*TMath::Sqrt([2]*x)");
2517   widthFit->SetParLimits(2, 1e-5, 1e5);
2518   width->Fit(widthFit, "", "", 5, 140);
2519   widthFit->SetRange(0, 200);
2520   widthFit->Draw("SAME");
2521
2522   // build new correction matrix
2523   TH2* new = (TH2*) proj->Clone("new");
2524   new->Reset();
2525   Float_t x, y;
2526   for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1)
2527   {
2528     TF1* func = (TF1*) gROOT->FindObject(fitWith);
2529     x = new->GetXaxis()->GetBinCenter(i);
2530     //if (i == 1)
2531     //  x = 0.1;
2532     x++;
2533     func->SetParameters(scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x));
2534     printf("%f %f %f %f\n", x, scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x));
2535
2536     for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
2537     {
2538       if (i < 11)
2539       {
2540         // leave bins 1..20 untouched
2541         new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j));
2542       }
2543       else
2544       {
2545         y = new->GetYaxis()->GetBinCenter(j);
2546         if (j == 1)
2547           y = 0.1;
2548         if (func->Eval(y) > 1e-4)
2549           new->SetBinContent(i, j, func->Eval(y));
2550       }
2551     }
2552   }
2553
2554   // fill 0 multiplicity bins, this cannot be done with the function because it does not accept 0
2555   // we take the values from the old response matrix
2556   //for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1)
2557   //  new->SetBinContent(i, 1, proj->GetBinContent(i, 1));
2558
2559   //for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
2560   //  new->SetBinContent(1, j, proj->GetBinContent(1, j));
2561
2562   // normalize correction for given nPart
2563   for (Int_t i=1; i<=new->GetNbinsX(); ++i)
2564   {
2565     Double_t sum = new->Integral(i, i, 1, proj->GetNbinsY());
2566     if (sum <= 0)
2567       continue;
2568
2569     for (Int_t j=1; j<=new->GetNbinsY(); ++j)
2570     {
2571       // npart sum to 1
2572       new->SetBinContent(i, j, new->GetBinContent(i, j) / sum);
2573       new->SetBinError(i, j, new->GetBinError(i, j) / sum);
2574     }
2575   }
2576
2577   new TCanvas;
2578   new->Draw("COLZ");
2579
2580   TH2* diff = (TH2*) new->Clone("diff");
2581   diff->Add(proj, -1);
2582
2583   new TCanvas;
2584   diff->Draw("COLZ");
2585   diff->SetMinimum(-0.05);
2586   diff->SetMaximum(0.05);
2587
2588   corr->Reset();
2589
2590   for (Int_t i=1; i<=new->GetNbinsX(); ++i)
2591     for (Int_t j=1; j<=new->GetNbinsY(); ++j)
2592       corr->SetBinContent(corr->GetXaxis()->GetNbins() / 2 + 1, i, j, new->GetBinContent(i, j));
2593
2594   new TCanvas;
2595   corr->Project3D("zy")->Draw("COLZ");
2596
2597   TFile::Open("out.root", "RECREATE");
2598   mult->SaveHistograms();
2599
2600   TH1* proj1 = proj->ProjectionY("proj1", 36, 36);
2601   TH1* proj2 = new->ProjectionY("proj2", 36, 36);
2602   proj2->SetLineColor(2);
2603
2604   new TCanvas;
2605   proj1->Draw();
2606   proj2->Draw("SAME");
2607 }
2608
2609 void buildCorrelationMap(const char* fileName = "multiplicityMC_2M.root", Int_t corrMatrix = 3)
2610 {
2611   gSystem->Load("libPWG0base");
2612
2613   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2614
2615   TFile::Open(fileName);
2616   mult->LoadHistograms("Multiplicity");
2617
2618   TH3F* new = mult->GetCorrelation(corrMatrix);
2619   new->Reset();
2620
2621   TF1* func = new TF1("func", "gaus(0)");
2622
2623   Int_t vtxBin = new->GetNbinsX() / 2;
2624   if (vtxBin == 0)
2625     vtxBin = 1;
2626
2627   Float_t sigma = 2;
2628   for (Int_t i=1; i<=new->GetYaxis()->GetNbins(); i+=1)
2629   {
2630     Float_t x = new->GetYaxis()->GetBinCenter(i);
2631     func->SetParameters(1, x * 0.8, sigma);
2632     //func->SetParameters(1, x, sigma);
2633
2634     for (Int_t j=1; j<=new->GetZaxis()->GetNbins(); j+=1)
2635     {
2636       Float_t y = new->GetYaxis()->GetBinCenter(j);
2637
2638       // cut at 1 sigma
2639       if (TMath::Abs(y-x*0.8) < sigma)
2640         new->SetBinContent(vtxBin, i, j, func->Eval(y));
2641
2642       // test only bin 40 has smearing
2643       //if (x != 40)
2644       //  new->SetBinContent(vtxBin, i, j, (i == j));
2645     }
2646   }
2647
2648   new TCanvas;
2649   new->Project3D("zy")->DrawCopy("COLZ");
2650
2651   TFile* file = TFile::Open("out.root", "RECREATE");
2652   mult->SetCorrelation(corrMatrix, new);
2653   mult->SaveHistograms();
2654   file->Close();
2655 }
2656
2657 void GetCrossSections(const char* fileName)
2658 {
2659   gSystem->Load("libPWG0base");
2660
2661   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2662
2663   TFile::Open(fileName);
2664   mult->LoadHistograms("Multiplicity");
2665
2666   TH1* xSection2 = mult->GetCorrelation(3)->Project3D("y")->Clone("xSection2");
2667   xSection2->Sumw2();
2668   xSection2->Scale(1.0 / xSection2->Integral());
2669
2670   TH1* xSection15 = mult->GetCorrelation(2)->Project3D("y")->Clone("xSection15");
2671   xSection15->Sumw2();
2672   xSection15->Scale(1.0 / xSection15->Integral());
2673
2674   TFile::Open("crosssection.root", "RECREATE");
2675   xSection2->Write();
2676   xSection15->Write();
2677   gFile->Close();
2678 }
2679
2680 void AnalyzeSpeciesTree(const char* fileName)
2681 {
2682   //
2683   // prints statistics about fParticleSpecies
2684   //
2685
2686   gSystem->Load("libPWG0base");
2687
2688   TFile::Open(fileName);
2689   TNtuple* fParticleSpecies = (TNtuple*) gFile->Get("fParticleSpecies");
2690
2691   const Int_t nFields = 8;
2692   Long_t totals[8];
2693   for (Int_t i=0; i<nFields; i++)
2694     totals[i] = 0;
2695
2696   for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++)
2697   {
2698     fParticleSpecies->GetEvent(i);
2699
2700     Float_t* f = fParticleSpecies->GetArgs();
2701
2702     for (Int_t j=0; j<nFields; j++)
2703       totals[j] += f[j+1];
2704   }
2705
2706   for (Int_t i=0; i<nFields; i++)
2707     Printf("%d --> %ld", i, totals[i]);
2708 }
2709
2710 void BuildResponseFromTree(const char* fileName, const char* target)
2711 {
2712   //
2713   // builds several response matrices with different particle ratios (systematic study)
2714   // particle-species study
2715   //
2716   // WARNING doesn't work uncompiled, see test.C
2717
2718   loadlibs();
2719
2720   TFile::Open(fileName);
2721   TNtuple* fParticleSpecies = (TNtuple*) gFile->Get("fParticleSpecies");
2722
2723   TFile* file = TFile::Open(target, "RECREATE");
2724   file->Close();
2725
2726   Int_t tracks = 0; // control variables
2727   Int_t noLabel = 0;
2728   Int_t secondaries = 0;
2729   Int_t doubleCount = 0;
2730
2731   for (Int_t num = 0; num < 7; num++)
2732   {
2733     Printf("%d", num);
2734
2735     TString name;
2736     name.Form("Multiplicity_%d", num);
2737     AliMultiplicityCorrection* fMultiplicity = new AliMultiplicityCorrection(name, name);
2738
2739     Float_t ratio[4]; // pi, K, p, other
2740     for (Int_t i = 0; i < 4; i++)
2741       ratio[i] = 1;
2742
2743     if (num == 1)
2744       ratio[1] = 0.7;
2745     if (num == 2)
2746       ratio[1] = 1.3;
2747
2748     if (num == 3)
2749       ratio[2] = 0.7;
2750     if (num == 4)
2751       ratio[2] = 1.3;
2752
2753     if (num == 5)
2754       ratio[3] = 0.7;
2755     if (num == 6)
2756       ratio[3] = 1.3;
2757
2758     for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++)
2759     {
2760       fParticleSpecies->GetEvent(i);
2761
2762
2763       Float_t* f = fParticleSpecies->GetArgs();
2764
2765       Float_t gene = 0;
2766       Float_t meas = 0;
2767
2768       for (Int_t j = 0; j < 4; j++)
2769       {
2770         if (ratio[j] == 1)
2771         {
2772           gene += f[j+1];
2773         }
2774         else
2775         {
2776           for (Int_t k=0; k<f[j+1]; k++)
2777           {
2778             if (ratio[j] < 1 && gRandom->Uniform() < ratio[j])
2779             {
2780               gene += 1;
2781             }
2782             else if (ratio[j] > 1)
2783             {
2784               gene += 1;
2785               if (gRandom->Uniform() < ratio[j] - 1)
2786                 gene += 1;
2787             }
2788           }
2789         }
2790
2791         if (ratio[j] == 1)
2792         {
2793           meas += f[j+1+4];
2794         }
2795         else
2796         {
2797           for (Int_t k=0; k<f[j+1+4]; k++)
2798           {
2799             if (ratio[j] < 1 && gRandom->Uniform() < ratio[j])
2800             {
2801               meas += 1;
2802             }
2803             else if (ratio[j] > 1)
2804             {
2805               meas += 1;
2806               if (gRandom->Uniform() < ratio[j] - 1)
2807                 meas += 1;
2808             }
2809           }
2810         }
2811
2812         tracks += f[j+1+4];
2813       }
2814
2815       // add the ones w/o label
2816       tracks += f[9];
2817       noLabel += f[9];
2818
2819       // secondaries are already part of meas!
2820       secondaries += f[10];
2821
2822       // double counted are already part of meas!
2823       doubleCount += f[11];
2824
2825       // 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!
2826       meas += f[9];
2827
2828       //Printf("%.f %.f %.f %.f %.f", f[5], f[6], f[7], f[8], f[9]);
2829
2830       fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, meas, meas, meas);
2831       // HACK all as kND = 1
2832       fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, 1, gene, gene, gene, gene);
2833       fMultiplicity->FillMeasured(f[0], meas, meas, meas);
2834     }
2835
2836     //fMultiplicity->DrawHistograms();
2837
2838     TFile* file = TFile::Open(target, "UPDATE");
2839     fMultiplicity->SaveHistograms();
2840     file->Close();
2841
2842     if (num == 0)
2843     {
2844       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);
2845       if ((Float_t) noLabel / tracks > 0.02)
2846         Printf("WARNING: More than 2%% of tracks without label, this might bias the study!");
2847     }
2848   }
2849 }
2850
2851 void ParticleRatioStudy()
2852 {
2853   loadlibs();
2854
2855   for (Int_t num = 0; num < 7; num++)
2856   {
2857     TString target;
2858     target.Form("chi2a_inel_species_%d.root", num);
2859
2860     correct("compositions.root", Form("Multiplicity_%d", num), "multiplicityESD.root", 1, -1, 0, -1, 2, target);
2861   }
2862 }
2863
2864 void MergeCrossSection(Int_t xsectionID, const char* output = "multiplicityMC_xsection.root")
2865 {
2866   const char* files[] = { "multiplicitySD.root", "multiplicityDD.root", "multiplicityND.root" };
2867
2868   loadlibs();
2869
2870   TFile::Open(output, "RECREATE");
2871   gFile->Close();
2872
2873   AliMultiplicityCorrection* data[3];
2874   TList list;
2875
2876   Float_t ref_SD = -1;
2877   Float_t ref_DD = -1;
2878   Float_t ref_ND = -1;
2879
2880   Float_t error_SD = -1;
2881   Float_t error_DD = -1;
2882   Float_t error_ND = -1;
2883
2884   gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawSystematics.C"));
2885   GetRelativeFractions(xsectionID, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
2886
2887   for (Int_t i=0; i<3; ++i)
2888   {
2889     TString name;
2890     name.Form("Multiplicity");
2891     if (i > 0)
2892       name.Form("Multiplicity_%d", i);
2893
2894     TFile::Open(files[i]);
2895     data[i] = new AliMultiplicityCorrection(name, name);
2896     data[i]->LoadHistograms("Multiplicity");
2897   }
2898
2899   // TODO is the under/overflow bin scaled as well? --> seems like, verify anyway!
2900
2901   // calculating relative
2902   Float_t sd = data[0]->GetMultiplicityINEL(0)->Integral(0, data[0]->GetMultiplicityINEL(0)->GetNbinsX()+1);
2903   Float_t dd = data[1]->GetMultiplicityINEL(0)->Integral(0, data[1]->GetMultiplicityINEL(0)->GetNbinsX()+1);
2904   Float_t nd = data[2]->GetMultiplicityINEL(0)->Integral(0, data[2]->GetMultiplicityINEL(0)->GetNbinsX()+1);
2905   Float_t total = nd + dd + sd;
2906
2907   nd /= total;
2908   sd /= total;
2909   dd /= total;
2910
2911   Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
2912
2913   Float_t ratio[3];
2914   ratio[0] = ref_SD / sd;
2915   ratio[1] = ref_DD / dd;
2916   ratio[2] = ref_ND / nd;
2917
2918   Printf("SD=%.2f, DD=%.2f, ND=%.2f",ratio[0], ratio[1], ratio[2]);
2919
2920   for (Int_t i=0; i<3; ++i)
2921   {
2922     // modify x-section
2923     for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++)
2924     {
2925       data[i]->GetMultiplicityVtx(j)->Scale(ratio[i]);
2926       data[i]->GetMultiplicityMB(j)->Scale(ratio[i]);
2927       data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]);
2928       data[i]->GetMultiplicityNSD(j)->Scale(ratio[i]);
2929     }
2930
2931     for (Int_t j=0; j<AliMultiplicityCorrection::kESDHists; j++)
2932     {
2933       data[i]->GetMultiplicityESD(j)->Scale(ratio[i]);
2934       data[i]->GetTriggeredEvents(j)->Scale(ratio[i]);
2935       data[i]->GetNoVertexEvents(j)->Scale(ratio[i]);
2936     }
2937
2938     for (Int_t j=0; j<AliMultiplicityCorrection::kCorrHists; j++)
2939       data[i]->GetCorrelation(j)->Scale(ratio[i]);
2940
2941     if (i > 0)
2942       list.Add(data[i]);
2943   }
2944
2945   printf("Case %d, %s: Entries in response matrix: SD: %.2f DD: %.2f ND: %.2f", xsectionID, data[0]->GetName(), data[0]->GetCorrelation(0)->Integral(), data[1]->GetCorrelation(0)->Integral(), data[2]->GetCorrelation(0)->Integral());
2946
2947   data[0]->Merge(&list);
2948
2949   Printf(" Total: %.2f", data[0]->GetCorrelation(0)->Integral());
2950
2951   TFile::Open(output, "RECREATE");
2952   data[0]->SaveHistograms();
2953   gFile->Close();
2954
2955   list.Clear();
2956
2957   for (Int_t i=0; i<3; ++i)
2958     delete data[i];
2959 }
2960
2961 void Rebin(const char* fileName = "multiplicityMC_3M.root", Int_t corrMatrix = 3)
2962 {
2963   // rebins MC axis of correlation map, MC and histogram for corrected (for evaluation of effect of regularization)
2964   // rebin does not exist for 3D hists, so we convert to 2D and then back to 3D (loosing the vertex information)
2965
2966   Printf("WARNING: Vertex information is lost in this process. Use result only for evaluation of errors.");
2967
2968   gSystem->Load("libPWG0base");
2969
2970   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2971
2972   TFile::Open(fileName);
2973   mult->LoadHistograms("Multiplicity");
2974
2975   // rebin correlation
2976   TH3* old = mult->GetCorrelation(corrMatrix);
2977
2978   // empty under/overflow bins in x, otherwise Project3D takes them into account
2979   for (Int_t y=1; y<=old->GetYaxis()->GetNbins(); ++y)
2980   {
2981     for (Int_t z=1; z<=old->GetZaxis()->GetNbins(); ++z)
2982     {
2983       old->SetBinContent(0, y, z, 0);
2984       old->SetBinContent(old->GetXaxis()->GetNbins()+1, y, z, 0);
2985     }
2986   }
2987
2988   TH2* response = (TH2*) old->Project3D("zy");
2989   response->RebinX(2);
2990
2991   TH3F* new = new TH3F(old->GetName(), old->GetTitle(),
2992     old->GetXaxis()->GetNbins(), old->GetXaxis()->GetBinLowEdge(1), old->GetXaxis()->GetBinUpEdge(old->GetXaxis()->GetNbins()),
2993     old->GetYaxis()->GetNbins() / 2, old->GetYaxis()->GetBinLowEdge(1), old->GetYaxis()->GetBinUpEdge(old->GetYaxis()->GetNbins()),
2994     old->GetZaxis()->GetNbins(), old->GetZaxis()->GetBinLowEdge(1), old->GetZaxis()->GetBinUpEdge(old->GetZaxis()->GetNbins()));
2995   new->Reset();
2996
2997   Int_t vtxBin = new->GetNbinsX() / 2;
2998   if (vtxBin == 0)
2999     vtxBin = 1;
3000
3001   for (Int_t i=1; i<=new->GetYaxis()->GetNbins(); i+=1)
3002     for (Int_t j=1; j<=new->GetZaxis()->GetNbins(); j+=1)
3003       new->SetBinContent(vtxBin, i, j, response->GetBinContent(i, j));
3004
3005   // rebin MC + hist for corrected
3006   for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kTrVtx; eventType <= AliMultiplicityCorrection::kINEL; eventType++)
3007     mult->GetMultiplicityMC(corrMatrix, eventType)->RebinY(2);
3008
3009   mult->GetMultiplicityESDCorrected(corrMatrix)->Rebin(2);
3010
3011   // recreate measured from correlation matrix to get rid of vertex shift effect
3012   TH2* newMeasured = (TH2*) old->Project3D("zx");
3013   TH2* esd = mult->GetMultiplicityESD(corrMatrix);
3014   esd->Reset();
3015
3016   // transfer from TH2D to TH2F
3017   for (Int_t i=0; i<=new->GetXaxis()->GetNbins()+1; i+=1)
3018     for (Int_t j=0; j<=new->GetYaxis()->GetNbins()+1; j+=1)
3019       esd->SetBinContent(i, j, newMeasured->GetBinContent(i, j));
3020
3021   new TCanvas;
3022   new->Project3D("zy")->DrawCopy("COLZ");
3023
3024   TFile* file = TFile::Open("out.root", "RECREATE");
3025   mult->SetCorrelation(corrMatrix, new);
3026   mult->SaveHistograms();
3027   file->Close();
3028 }
3029
3030 void EvaluateRegularizationEffect(Int_t step, const char* fileNameRebinned = "multiplicityMC_3M_rebinned.root", const char* fileNameNormal = "multiplicityMC_3M.root", Int_t histID = 3)
3031 {
3032   // due to some static members in AliMultiplicityCorrection, the session has to be restarted after changing the number of parameters, to be fixed
3033   // that is why this is done in 2 steps
3034
3035   gSystem->Load("libPWG0base");
3036
3037   Bool_t fullPhaseSpace = kFALSE;
3038
3039   if (step == 1)
3040   {
3041     // first step: unfold without regularization and rebinned histogram ("N=M")
3042     AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3043     TFile::Open(fileNameRebinned);
3044     mult->LoadHistograms();
3045
3046     mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125);
3047     mult->SetCreateBigBin(kFALSE);
3048
3049     mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
3050     mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, mult->GetMultiplicityVtx(histID)->ProjectionY("mymchist"));
3051
3052     TFile* file = TFile::Open("EvaluateRegularizationEffect1.root", "RECREATE");
3053     mult->SaveHistograms();
3054     file->Close();
3055   }
3056   else if (step == 2)
3057   {
3058     // second step: unfold with regularization and normal histogram
3059     AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3060     TFile::Open(fileNameNormal);
3061     mult2->LoadHistograms();
3062
3063     mult2->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
3064     mult2->SetCreateBigBin(kTRUE);
3065     mult2->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
3066     mult2->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, mult2->GetMultiplicityVtx(histID)->ProjectionY("mymchist"));
3067
3068     TH1* result2 = mult2->GetMultiplicityESDCorrected(histID);
3069
3070     AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3071     TFile* file = TFile::Open("EvaluateRegularizationEffect1.root");
3072     mult->LoadHistograms();
3073
3074     TH1* result1 = mult->GetMultiplicityESDCorrected(histID);
3075
3076     // compare results
3077     TCanvas* canvas = new TCanvas("EvaluateRegularizationEffect", "EvaluateRegularizationEffect", 1000, 800);
3078     canvas->Divide(2, 2);
3079
3080     canvas->cd(1);
3081     result1->SetLineColor(1);
3082     result1->DrawCopy();
3083     result2->SetLineColor(2);
3084     result2->DrawCopy("SAME");
3085     gPad->SetLogy();
3086
3087     result2->Rebin(2);
3088     result1->Scale(1.0 / result1->Integral());
3089     result2->Scale(1.0 / result2->Integral());
3090
3091     canvas->cd(2);
3092     result1->DrawCopy();
3093     result2->DrawCopy("SAME");
3094     gPad->SetLogy();
3095
3096     TH1* diff = (TH1*) result1->Clone("diff");
3097     diff->Add(result2, -1);
3098
3099     canvas->cd(3);
3100     diff->DrawCopy("HIST");
3101
3102     canvas->cd(4);
3103     diff->Divide(result1);
3104     diff->GetYaxis()->SetRangeUser(-0.3, 0.3);
3105     diff->DrawCopy("HIST");
3106
3107     Double_t chi2 = 0;
3108     for (Int_t i=1; i<=diff->GetNbinsX(); i++)
3109       chi2 += diff->GetBinContent(i) * diff->GetBinContent(i);
3110
3111     Printf("Chi2 is %e", chi2);
3112
3113     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
3114   }
3115 }
3116
3117 void MergeDistributions()
3118 {
3119   loadlibs();
3120
3121   const char* dir1 = "run000104824-52_pass4";
3122   const char* dir2 = "run000104867_90_92_pass4";
3123
3124   for (Int_t evType = 0; evType < 2; evType++)
3125   {
3126     Printf("%d", evType);
3127
3128     const char* evTypeStr = ((evType == 0) ? "inel" : "nsd");
3129
3130     const char* id = "chi2a";
3131     //const char* id = "bayes";
3132
3133     TString suffix;
3134     suffix.Form("/all/spd/%s_", id);
3135     if (evType == 1)
3136       suffix.Form("/v0and/spd/%s_", id);
3137
3138     TString file1, file2;
3139     file1.Form("%s%s%%s.root", dir1, suffix.Data());
3140     file2.Form("%s%s%%s.root", dir2, suffix.Data());
3141
3142     if (1)
3143     {
3144       const char* files[] = { Form(file1.Data(), evTypeStr), Form(file2.Data(), evTypeStr) };
3145       Merge(2, files, Form("merged/%s_%s.root", id, evTypeStr));
3146     }
3147     else
3148     {
3149       AliMultiplicityCorrection* mult1 = AliMultiplicityCorrection::Open(Form(file1.Data(), evTypeStr));
3150       AliMultiplicityCorrection* mult2 = AliMultiplicityCorrection::Open(Form(file2.Data(), evTypeStr));
3151
3152       AliMultiplicityCorrection* target = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3153
3154       for (Int_t etaRange = 0; etaRange < 3; etaRange++)
3155       {
3156         hist1 = mult1->GetMultiplicityESDCorrected(etaRange);
3157         hist2 = mult2->GetMultiplicityESDCorrected(etaRange);
3158         targetHist = target->GetMultiplicityESDCorrected(etaRange);
3159
3160         //hist1->Scale(1.0 / hist1->Integral());
3161         //hist2->Scale(1.0 / hist2->Integral());
3162
3163         //targetHist->SetBinContent(1, hist1->GetBinContent(1) * 0.5 + hist2->GetBinContent(1) * 0.5);
3164         targetHist->SetBinContent(1, hist1->GetBinContent(1) + hist2->GetBinContent(1));
3165         for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
3166         {
3167           if (hist1->GetBinError(i) > 0 && hist2->GetBinError(i) > 0)
3168           {
3169             Float_t w1 = 1.0 / hist1->GetBinError(i) / hist1->GetBinError(i);
3170             Float_t w2 = 1.0 / hist2->GetBinError(i) / hist2->GetBinError(i);
3171
3172             Float_t average = (hist1->GetBinContent(i) * w1 + hist2->GetBinContent(i) * w2) / (w1 + w2);
3173
3174             //targetHist->SetBinContent(i, average);
3175             //targetHist->SetBinError(i, TMath::Max(mult1->GetBinError(i), mult2->GetBinError(i)));
3176
3177             targetHist->SetBinContent(i, hist1->GetBinContent(i) + hist2->GetBinContent(i));
3178             targetHist->SetBinError(i, hist1->GetBinError(i) + hist2->GetBinError(i));
3179           }
3180         }
3181       }
3182
3183       file = TFile::Open(Form("merged/%s_%s.root", id, evTypeStr), "RECREATE");
3184       target->SaveHistograms();
3185       file->Close();
3186     }
3187
3188     const char* files2[] = { Form(file1.Data(), evTypeStr), Form(file2.Data(), evTypeStr), Form("merged/%s_%s.root", id, evTypeStr) };
3189     drawAll(files2, (evType == 1), kTRUE);
3190   }
3191 }
3192
3193 void CompareDistributions(Int_t evType)
3194 {
3195   loadlibs();
3196   gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawPlots.C"));
3197
3198   const char* dir1 = "run000104824-52_pass4";
3199   const char* dir2 = "run000104867_90_92_pass4";
3200
3201   const char* evTypeStr = (evType == 0) ? "inel" : "nsd";
3202
3203   const char* suffix = "/all/spd/chi2a_";
3204   if (evType == 1)
3205     suffix = "/v0and/spd/chi2a_";
3206
3207   TString file1, file2;
3208   file1.Form("%s%s%s.root", dir1, suffix, evTypeStr);
3209   file2.Form("%s%s%s.root", dir2, suffix, evTypeStr);
3210
3211   for (Int_t hID = 0; hID < 3; hID++)
3212     CompareQualityHists(file1, file2, Form("Multiplicity/fMultiplicityESDCorrected%d", hID), 1, 1);
3213 }
3214
3215 void StatisticalUncertainty(Int_t methodType, Int_t etaRange, Bool_t mc = kFALSE)
3216 {
3217   loadlibs();
3218
3219   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root");
3220   AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open("multiplicityESD.root");
3221   AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
3222
3223   TH1* mcHist = esd->GetMultiplicityMB(etaRange)->ProjectionY("mymc", 1, 1);
3224
3225   Int_t geneLimits[] = { 0, 0, 0 };
3226
3227   LoadAndInitialize(mult, esd, multTrigger, etaRange, kFALSE, geneLimits);
3228
3229   AliUnfolding::SetNbins(kBinLimits[etaRange], geneLimits[etaRange]);
3230   AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, 5);
3231   AliUnfolding::SetBayesianParameters(1, 10);
3232
3233   // background subtraction
3234   Int_t background = 0;
3235
3236   //background = 1091 + 4398; // MB1 for run 104824 - 52
3237   background = 417 + 1758; // MB1 for run 104867 - 92
3238
3239   //background = 20 + 0; // V0AND for run 104824 - 52
3240   //background = 10 + 0; // V0AND for run 104867 - 92
3241
3242   Printf("NOTE: Subtracting %d background events", background);
3243
3244   Int_t zeroBin = mult->GetTriggeredEvents(etaRange)->GetBinContent(1) - background - mult->GetMultiplicityESD(etaRange)->Integral(0, 2, 1, 1);
3245
3246   TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
3247
3248   new TCanvas; errorMeasured->Draw();
3249   new TCanvas;
3250
3251   AliPWG0Helper::NormalizeToBinWidth(mult->GetMultiplicityESDCorrected(etaRange));
3252   mult->GetMultiplicityESDCorrected(etaRange)->Draw();
3253
3254   if (mc)
3255   {
3256     AliPWG0Helper::NormalizeToBinWidth(mcHist);
3257     mcHist->Scale(mult->GetMultiplicityESDCorrected(etaRange)->Integral() / mcHist->Integral());
3258     mcHist->SetLineColor(2);
3259     mcHist->Draw("SAME");
3260   }
3261   gPad->SetLogy();
3262
3263   TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
3264
3265   TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
3266
3267   if (mc)
3268   {
3269     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
3270     DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
3271   }
3272
3273   TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
3274   errorResponse->Write();
3275   errorMeasured->Write();
3276   errorBoth->Write();
3277   file->Close();
3278 }
3279
3280 void ChangeResponseMatrixEfficiency(Bool_t reduceEff = kTRUE, Float_t factor = 0.01625)
3281 {
3282   loadlibs();
3283
3284   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicity.root");
3285
3286   for (Int_t etaR = 0; etaR < 3; etaR++)
3287   {
3288     TH3* corr = mult->GetCorrelation(etaR);
3289     TH3* corrOld = (TH3*) corr->Clone("corrOld");
3290
3291     corr->Reset();
3292
3293     Int_t total = 0;
3294     Int_t change = 0;
3295     for (Int_t x=0; x<=corr->GetNbinsX()+1; x++)
3296     {
3297       for (Int_t y=0; y<=corr->GetNbinsY()+1; y++)
3298       {
3299         Printf("%d", y);
3300         for (Int_t z=0; z<=corr->GetNbinsZ()+1; z++)
3301         {
3302           Float_t tracklets = corr->GetZaxis()->GetBinCenter(z);
3303
3304           for (Int_t n=0; n<corrOld->GetBinContent(x, y, z); n++)
3305           {
3306             total += tracklets;
3307             Float_t trackletsNew = tracklets;
3308
3309             for (Int_t i=0; i<tracklets; i++)
3310             {
3311               if (gRandom->Uniform() < factor)
3312               {
3313                 trackletsNew += (reduceEff) ? -1 : 1;
3314                 change += (reduceEff) ? -1 : 1;
3315               }
3316             }
3317
3318             corr->Fill(corr->GetXaxis()->GetBinCenter(x), corr->GetYaxis()->GetBinCenter(y), trackletsNew);
3319           }
3320         }
3321       }
3322     }
3323
3324     Printf("%d: Changed %d out of %d total tracklets", etaR, change, total);
3325   }
3326
3327   file = TFile::Open("out.root", "RECREATE");
3328   mult->SaveHistograms();
3329   file->Close();
3330 }
3331