PWG0/dNdEta/AliMultiplicityCorrection.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // This class is used to store correction maps, raw input and results of the multiplicity
  19 // measurement with the ITS or TPC
  20 // It also contains functions to correct the spectrum using different methods.
  21 //
  22 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
  23
  24 #include "AliMultiplicityCorrection.h"
  25
  26 #include <TFile.h>
  27 #include <TH1F.h>
  28 #include <TH2F.h>
  29 #include <TH3F.h>
  30 #include <TDirectory.h>
  31 #include <TVirtualFitter.h>
  32 #include <TCanvas.h>
  33 #include <TString.h>
  34 #include <TF1.h>
  35 #include <TMath.h>
  36 #include <TCollection.h>
  37 #include <TLegend.h>
  38 #include <TLine.h>
  39
  40 ClassImp(AliMultiplicityCorrection)
  41
  42 const Int_t AliMultiplicityCorrection::fgMaxInput  = 250;  // bins in measured histogram
  43 const Int_t AliMultiplicityCorrection::fgMaxParams = 250;  // bins in unfolded histogram = number of fit params
  44
  45 TH1* AliMultiplicityCorrection::fCurrentESD = 0;
  46 TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
  47 TH1* AliMultiplicityCorrection::fCurrentEfficiency = 0;
  48 TMatrixD* AliMultiplicityCorrection::fCorrelationMatrix = 0;
  49 TMatrixD* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0;
  50 TVectorD* AliMultiplicityCorrection::fCurrentESDVector = 0;
  51 TVectorD* AliMultiplicityCorrection::fEntropyAPriori = 0;
  52 AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegularizationType = AliMultiplicityCorrection::kPol1;
  53 Float_t AliMultiplicityCorrection::fRegularizationWeight = 5000;
  54 TF1* AliMultiplicityCorrection::fNBD = 0;
  55
  56 //____________________________________________________________________
  57 AliMultiplicityCorrection::AliMultiplicityCorrection() :
  58   TNamed(), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0)
  59 {
  60   //
  61   // default constructor
  62   //
  63
  64   for (Int_t i = 0; i < kESDHists; ++i)
  65     fMultiplicityESD[i] = 0;
  66
  67   for (Int_t i = 0; i < kMCHists; ++i)
  68   {
  69     fMultiplicityVtx[i] = 0;
  70     fMultiplicityMB[i] = 0;
  71     fMultiplicityINEL[i] = 0;
  72   }
  73
  74   for (Int_t i = 0; i < kCorrHists; ++i)
  75   {
  76     fCorrelation[i] = 0;
  77     fMultiplicityESDCorrected[i] = 0;
  78   }
  79 }
  80
  81 //____________________________________________________________________
  82 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
  83   TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0)
  84 {
  85   //
  86   // named constructor
  87   //
  88
  89   // do not add this hists to the directory
  90   Bool_t oldStatus = TH1::AddDirectoryStatus();
  91   TH1::AddDirectory(kFALSE);
  92
  93   /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
  94   Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
  95                           10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
  96                           20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
  97                           30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
  98                           40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
  99                           50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
 100                           100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
 101                           150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
 102                           250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
 103                           500.5 };
 104                                 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
 105                           //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
 106                           //1000.5 };
 107
 108   #define VTXBINNING 10, binLimitsVtx
 109   #define NBINNING fgMaxParams, binLimitsN*/
 110
 111   #define NBINNING 501, -0.5, 500.5
 112   #define VTXBINNING 1, -10, 10
 113
 114   for (Int_t i = 0; i < kESDHists; ++i)
 115     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
 116
 117   for (Int_t i = 0; i < kMCHists; ++i)
 118   {
 119     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
 120     fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
 121
 122     fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
 123     fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
 124
 125     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
 126     fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
 127   }
 128
 129   for (Int_t i = 0; i < kCorrHists; ++i)
 130   {
 131     fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
 132     fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
 133   }
 134
 135   TH1::AddDirectory(oldStatus);
 136 }
 137
 138 //____________________________________________________________________
 139 AliMultiplicityCorrection::~AliMultiplicityCorrection()
 140 {
 141   //
 142   // Destructor
 143   //
 144 }
 145
 146 //____________________________________________________________________
 147 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
 148 {
 149   // Merge a list of AliMultiplicityCorrection objects with this (needed for
 150   // PROOF).
 151   // Returns the number of merged objects (including this).
 152
 153   if (!list)
 154     return 0;
 155
 156   if (list->IsEmpty())
 157     return 1;
 158
 159   TIterator* iter = list->MakeIterator();
 160   TObject* obj;
 161
 162   // collections of all histograms
 163   TList collections[kESDHists+kMCHists*3+kCorrHists*2];
 164
 165   Int_t count = 0;
 166   while ((obj = iter->Next())) {
 167
 168     AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
 169     if (entry == 0)
 170       continue;
 171
 172     for (Int_t i = 0; i < kESDHists; ++i)
 173       collections[i].Add(entry->fMultiplicityESD[i]);
 174
 175     for (Int_t i = 0; i < kMCHists; ++i)
 176     {
 177       collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
 178       collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
 179       collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
 180     }
 181
 182     for (Int_t i = 0; i < kCorrHists; ++i)
 183       collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]);
 184
 185     for (Int_t i = 0; i < kCorrHists; ++i)
 186       collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
 187
 188     count++;
 189   }
 190
 191   for (Int_t i = 0; i < kESDHists; ++i)
 192     fMultiplicityESD[i]->Merge(&collections[i]);
 193
 194   for (Int_t i = 0; i < kMCHists; ++i)
 195   {
 196     fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
 197     fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
 198     fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
 199   }
 200
 201   for (Int_t i = 0; i < kCorrHists; ++i)
 202     fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
 203
 204   for (Int_t i = 0; i < kCorrHists; ++i)
 205     fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]);
 206
 207   delete iter;
 208
 209   return count+1;
 210 }
 211
 212 //____________________________________________________________________
 213 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
 214 {
 215   //
 216   // loads the histograms from a file
 217   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
 218   //
 219
 220   if (!dir)
 221     dir = GetName();
 222
 223   if (!gDirectory->cd(dir))
 224     return kFALSE;
 225
 226   // TODO memory leak. old histograms needs to be deleted.
 227
 228   Bool_t success = kTRUE;
 229
 230   for (Int_t i = 0; i < kESDHists; ++i)
 231   {
 232     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
 233     if (!fMultiplicityESD[i])
 234       success = kFALSE;
 235   }
 236
 237   for (Int_t i = 0; i < kMCHists; ++i)
 238   {
 239     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
 240     fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
 241     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
 242     if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
 243       success = kFALSE;
 244   }
 245
 246   for (Int_t i = 0; i < kCorrHists; ++i)
 247   {
 248     fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
 249     if (!fCorrelation[i])
 250       success = kFALSE;
 251     fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
 252     if (!fMultiplicityESDCorrected[i])
 253       success = kFALSE;
 254   }
 255
 256   gDirectory->cd("..");
 257
 258   return success;
 259 }
 260
 261 //____________________________________________________________________
 262 void AliMultiplicityCorrection::SaveHistograms()
 263 {
 264   //
 265   // saves the histograms
 266   //
 267
 268   gDirectory->mkdir(GetName());
 269   gDirectory->cd(GetName());
 270
 271   for (Int_t i = 0; i < kESDHists; ++i)
 272     if (fMultiplicityESD[i])
 273       fMultiplicityESD[i]->Write();
 274
 275   for (Int_t i = 0; i < kMCHists; ++i)
 276   {
 277     if (fMultiplicityVtx[i])
 278       fMultiplicityVtx[i]->Write();
 279     if (fMultiplicityMB[i])
 280       fMultiplicityMB[i]->Write();
 281     if (fMultiplicityINEL[i])
 282       fMultiplicityINEL[i]->Write();
 283   }
 284
 285   for (Int_t i = 0; i < kCorrHists; ++i)
 286   {
 287     if (fCorrelation[i])
 288       fCorrelation[i]->Write();
 289     if (fMultiplicityESDCorrected[i])
 290       fMultiplicityESDCorrected[i]->Write();
 291   }
 292
 293   gDirectory->cd("..");
 294 }
 295
 296 //____________________________________________________________________
 297 void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
 298 {
 299   //
 300   // Fills an event from MC
 301   //
 302
 303   if (triggered)
 304   {
 305     fMultiplicityMB[0]->Fill(vtx, generated05);
 306     fMultiplicityMB[1]->Fill(vtx, generated10);
 307     fMultiplicityMB[2]->Fill(vtx, generated15);
 308     fMultiplicityMB[3]->Fill(vtx, generated20);
 309     fMultiplicityMB[4]->Fill(vtx, generatedAll);
 310
 311     if (vertex)
 312     {
 313       fMultiplicityVtx[0]->Fill(vtx, generated05);
 314       fMultiplicityVtx[1]->Fill(vtx, generated10);
 315       fMultiplicityVtx[2]->Fill(vtx, generated15);
 316       fMultiplicityVtx[3]->Fill(vtx, generated20);
 317       fMultiplicityVtx[4]->Fill(vtx, generatedAll);
 318     }
 319   }
 320
 321   fMultiplicityINEL[0]->Fill(vtx, generated05);
 322   fMultiplicityINEL[1]->Fill(vtx, generated10);
 323   fMultiplicityINEL[2]->Fill(vtx, generated15);
 324   fMultiplicityINEL[3]->Fill(vtx, generated20);
 325   fMultiplicityINEL[4]->Fill(vtx, generatedAll);
 326 }
 327
 328 //____________________________________________________________________
 329 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
 330 {
 331   //
 332   // Fills an event from ESD
 333   //
 334
 335   fMultiplicityESD[0]->Fill(vtx, measured05);
 336   fMultiplicityESD[1]->Fill(vtx, measured10);
 337   fMultiplicityESD[2]->Fill(vtx, measured15);
 338   fMultiplicityESD[3]->Fill(vtx, measured20);
 339 }
 340
 341 //____________________________________________________________________
 342 void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
 343 {
 344   //
 345   // Fills an event into the correlation map with the information from MC and ESD
 346   //
 347
 348   fCorrelation[0]->Fill(vtx, generated05, measured05);
 349   fCorrelation[1]->Fill(vtx, generated10, measured10);
 350   fCorrelation[2]->Fill(vtx, generated15, measured15);
 351   fCorrelation[3]->Fill(vtx, generated20, measured20);
 352
 353   fCorrelation[4]->Fill(vtx, generatedAll, measured05);
 354   fCorrelation[5]->Fill(vtx, generatedAll, measured10);
 355   fCorrelation[6]->Fill(vtx, generatedAll, measured15);
 356   fCorrelation[7]->Fill(vtx, generatedAll, measured20);
 357 }
 358
 359 //____________________________________________________________________
 360 Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
 361 {
 362   // homogenity term for minuit fitting
 363   // pure function of the parameters
 364   // prefers constant function (pol0)
 365
 366   Double_t chi2 = 0;
 367
 368   // ignore the first bin here. on purpose...
 369   for (Int_t i=2; i<fgMaxParams; ++i)
 370   {
 371     Double_t right  = params[i];
 372     Double_t left   = params[i-1];
 373
 374     if (right != 0)
 375     {
 376       Double_t diff = 1 - left / right;
 377       chi2 += diff * diff;
 378     }
 379   }
 380
 381   return chi2 / 100.0;
 382 }
 383
 384 //____________________________________________________________________
 385 Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
 386 {
 387   // homogenity term for minuit fitting
 388   // pure function of the parameters
 389   // prefers linear function (pol1)
 390
 391   Double_t chi2 = 0;
 392
 393   // ignore the first bin here. on purpose...
 394   for (Int_t i=2+1; i<fgMaxParams; ++i)
 395   {
 396     if (params[i-1] == 0)
 397       continue;
 398
 399     Double_t right  = params[i];
 400     Double_t middle = params[i-1];
 401     Double_t left   = params[i-2];
 402
 403     Double_t der1 = (right - middle);
 404     Double_t der2 = (middle - left);
 405
 406     Double_t diff = (der1 - der2) / middle;
 407
 408     chi2 += diff * diff;
 409   }
 410
 411   return chi2;
 412 }
 413
 414 //____________________________________________________________________
 415 Double_t AliMultiplicityCorrection::RegularizationTest(TVectorD& params)
 416 {
 417   // homogenity term for minuit fitting
 418   // pure function of the parameters
 419   // prefers linear function (pol1)
 420
 421   Double_t chi2 = 0;
 422
 423   Float_t der[fgMaxParams];
 424
 425   for (Int_t i=0; i<fgMaxParams-1; ++i)
 426     der[i] = params[i+1] - params[i];
 427
 428   for (Int_t i=0; i<fgMaxParams-6; ++i)
 429   {
 430     for (Int_t j=1; j<=5; ++j)
 431     {
 432       Double_t diff = der[i+j] - der[i];
 433       chi2 += diff * diff;
 434     }
 435   }
 436
 437   return chi2 * 100; // 10000
 438 }
 439
 440 //____________________________________________________________________
 441 Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
 442 {
 443   // homogenity term for minuit fitting
 444   // pure function of the parameters
 445   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
 446   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
 447
 448   Double_t chi2 = 0;
 449
 450   // ignore the first bin here. on purpose...
 451   for (Int_t i=3; i<fgMaxParams; ++i)
 452   {
 453     Double_t right  = params[i];
 454     Double_t middle = params[i-1];
 455     Double_t left   = params[i-2];
 456
 457     Double_t der1 = (right - middle);
 458     Double_t der2 = (middle - left);
 459
 460     Double_t diff = (der1 - der2);
 461
 462     chi2 += diff * diff;
 463   }
 464
 465   return chi2 * 1e4;
 466 }
 467
 468 //____________________________________________________________________
 469 Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
 470 {
 471   // homogenity term for minuit fitting
 472   // pure function of the parameters
 473   // calculates entropy, from
 474   // The method of reduced cross-entropy (M. Schmelling 1993)
 475
 476   Double_t paramSum = 0;
 477   // ignore the first bin here. on purpose...
 478   for (Int_t i=1; i<fgMaxParams; ++i)
 479     paramSum += params[i];
 480
 481   Double_t chi2 = 0;
 482   for (Int_t i=1; i<fgMaxParams; ++i)
 483   {
 484     Double_t tmp = params[i] / paramSum;
 485     if (tmp > 0 && (*fEntropyAPriori)[i] > 0)
 486       chi2 += tmp * log(tmp / (*fEntropyAPriori)[i]);
 487   }
 488
 489   return 10.0 + chi2;
 490 }
 491
 492 //____________________________________________________________________
 493 void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
 494 {
 495   //
 496   // fit function for minuit
 497   // does: nbd
 498   // func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
 499   // func->SetParNames("scaling", "averagen", "k");
 500   // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
 501   // func->SetParLimits(1, 0.001, 1000);
 502   // func->SetParLimits(2, 0.001, 1000);
 503   //
 504
 505   fNBD->SetParameters(params[0], params[1], params[2]);
 506
 507   Double_t params2[fgMaxParams];
 508
 509   for (Int_t i=0; i<fgMaxParams; ++i)
 510     params2[i] = fNBD->Eval(i);
 511
 512   MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
 513
 514   printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2);
 515 }
 516
 517 //____________________________________________________________________
 518 void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
 519 {
 520   //
 521   // fit function for minuit
 522   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
 523   //
 524
 525   // d
 526   TVectorD paramsVector(fgMaxParams);
 527   for (Int_t i=0; i<fgMaxParams; ++i)
 528     paramsVector[i] = params[i] * params[i];
 529
 530   // calculate penalty factor
 531   Double_t penaltyVal = 0;
 532   switch (fRegularizationType)
 533   {
 534     case kNone:       break;
 535     case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
 536     case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
 537     case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
 538     case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
 539     case kTest:       penaltyVal = RegularizationTest(paramsVector); break;
 540   }
 541
 542   //if (penaltyVal > 1e10)
 543   //  paramsVector2.Print();
 544
 545   penaltyVal *= fRegularizationWeight;
 546
 547   // Ad
 548   TVectorD measGuessVector(fgMaxInput);
 549   measGuessVector = (*fCorrelationMatrix) * paramsVector;
 550
 551   // Ad - m
 552   measGuessVector -= (*fCurrentESDVector);
 553
 554   TVectorD copy(measGuessVector);
 555
 556   // (Ad - m) W
 557   measGuessVector *= (*fCorrelationCovarianceMatrix);
 558
 559   //measGuessVector.Print();
 560
 561   // (Ad - m) W (Ad - m)
 562   Double_t chi2FromFit = measGuessVector * copy * 1e6;
 563
 564   chi2 = chi2FromFit + penaltyVal;
 565
 566   static Int_t callCount = 0;
 567   if ((callCount++ % 10000) == 0)
 568     printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2);
 569 }
 570
 571 //____________________________________________________________________
 572 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
 573 {
 574   //
 575   // fills fCurrentESD, fCurrentCorrelation
 576   // resets fMultiplicityESDCorrected
 577   //
 578
 579   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 580   fMultiplicityESDCorrected[correlationID]->Reset();
 581
 582   fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY();
 583   fCurrentESD->Sumw2();
 584
 585   // empty under/overflow bins in x, otherwise Project3D takes them into account
 586   TH3* hist = fCorrelation[correlationID];
 587   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
 588   {
 589     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
 590     {
 591       hist->SetBinContent(0, y, z, 0);
 592       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
 593     }
 594   }
 595
 596   fCurrentCorrelation = hist->Project3D("zy");
 597   //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1);
 598   //fMultiplicityESDCorrected[correlationID]->Rebin(2);
 599   fCurrentCorrelation->Sumw2();
 600
 601   if (createBigBin)
 602   {
 603     Int_t maxBin = 0;
 604     for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
 605     {
 606       if (fCurrentESD->GetBinContent(i) <= 5)
 607       {
 608         maxBin = i;
 609         break;
 610       }
 611     }
 612
 613     if (maxBin > 0)
 614     {
 615       TCanvas* canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
 616       canvas->Divide(2, 2);
 617
 618       canvas->cd(1);
 619       fCurrentESD->DrawCopy();
 620       gPad->SetLogy();
 621
 622       canvas->cd(2);
 623       fCurrentCorrelation->DrawCopy("COLZ");
 624
 625       printf("Bin limit in measured spectrum is %d.\n", maxBin);
 626       fCurrentESD->SetBinContent(maxBin, fCurrentESD->Integral(maxBin, fCurrentESD->GetNbinsX()));
 627       for (Int_t i=maxBin+1; i<=fCurrentESD->GetNbinsX(); ++i)
 628       {
 629         fCurrentESD->SetBinContent(i, 0);
 630         fCurrentESD->SetBinError(i, 0);
 631       }
 632       // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
 633       fCurrentESD->SetBinError(maxBin, TMath::Sqrt(fCurrentESD->GetBinContent(maxBin)));
 634
 635       printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(maxBin), fCurrentESD->GetBinError(maxBin));
 636
 637       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
 638       {
 639         fCurrentCorrelation->SetBinContent(i, maxBin, fCurrentCorrelation->Integral(i, i, maxBin, fCurrentCorrelation->GetNbinsY()));
 640         // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
 641         fCurrentCorrelation->SetBinError(i, maxBin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, maxBin)));
 642
 643         for (Int_t j=maxBin+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
 644         {
 645           fCurrentCorrelation->SetBinContent(i, j, 0);
 646           fCurrentCorrelation->SetBinError(i, j, 0);
 647         }
 648       }
 649
 650       printf("Adjusted correlation matrix!\n");
 651
 652       canvas->cd(3);
 653       fCurrentESD->DrawCopy();
 654       gPad->SetLogy();
 655
 656       canvas->cd(4);
 657       fCurrentCorrelation->DrawCopy("COLZ");
 658     }
 659   }
 660
 661   //normalize ESD
 662   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 663
 664   fCurrentEfficiency = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency");
 665   TH2* divisor = 0;
 666   switch (eventType)
 667   {
 668     case kTrVtx : divisor = fMultiplicityVtx[inputRange]; break;
 669     case kMB: divisor = fMultiplicityMB[inputRange]; break;
 670     case kINEL: divisor = fMultiplicityINEL[inputRange]; break;
 671   }
 672   fCurrentEfficiency->Divide(divisor->ProjectionY());
 673   //fCurrentEfficiency->Rebin(2);
 674   //fCurrentEfficiency->Scale(0.5);
 675 }
 676
 677 //____________________________________________________________________
 678 void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight)
 679 {
 680   fRegularizationType = type;
 681   fRegularizationWeight = weight;
 682
 683   printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight);
 684 }
 685
 686 //____________________________________________________________________
 687 Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* inputDist)
 688 {
 689   //
 690   // correct spectrum using minuit chi2 method
 691   //
 692   // if check is kTRUE the input MC solution (by definition the right solution) is used
 693   // no fit is made and just the chi2 is calculated
 694   //
 695
 696   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 697   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
 698
 699   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); // TODO FAKE kTRUE
 700
 701   fCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
 702   fCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
 703   fCurrentESDVector = new TVectorD(fgMaxInput);
 704   fEntropyAPriori = new TVectorD(fgMaxParams);
 705
 706   /*new TCanvas; fCurrentESD->DrawCopy();
 707   fCurrentESD = ((TH2*) fCurrentCorrelation)->ProjectionY("check-proj2");
 708   fCurrentESD->Sumw2();
 709   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 710   fCurrentESD->SetLineColor(2);
 711   fCurrentESD->DrawCopy("SAME");*/
 712
 713   // normalize correction for given nPart
 714   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
 715   {
 716     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
 717     if (sum <= 0)
 718       continue;
 719     Float_t maxValue = 0;
 720     Int_t maxBin = -1;
 721     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
 722     {
 723       // find most probably value
 724       if (maxValue < fCurrentCorrelation->GetBinContent(i, j))
 725       {
 726         maxValue = fCurrentCorrelation->GetBinContent(i, j);
 727         maxBin = j;
 728       }
 729
 730       // npart sum to 1
 731       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
 732       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
 733
 734       if (i <= fgMaxParams && j <= fgMaxInput)
 735         (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
 736     }
 737
 738     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
 739   }
 740
 741   // Initialize TMinuit via generic fitter interface
 742   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
 743   Double_t arglist[100];
 744   arglist[0] = 0;
 745   minuit->ExecuteCommand("SET PRINT", arglist, 1);
 746
 747   minuit->SetFCN(MinuitFitFunction);
 748
 749   for (Int_t i=0; i<fgMaxParams; i++)
 750     (*fEntropyAPriori)[i] = 1;
 751
 752   if (inputDist)
 753   {
 754     printf("Using different starting conditions...\n");
 755     new TCanvas;
 756     inputDist->DrawCopy();
 757
 758     inputDist->Scale(1.0 / inputDist->Integral());
 759     for (Int_t i=0; i<fgMaxParams; i++)
 760       if (inputDist->GetBinContent(i+1) > 0)
 761         (*fEntropyAPriori)[i] = inputDist->GetBinContent(i+1);
 762   }
 763   else
 764     inputDist = fCurrentESD;
 765
 766
 767   //Float_t minStartValue = 0; //1e-3;
 768
 769   //new TCanvas; fMultiplicityVtx[mcTarget]->Draw("COLZ");
 770   TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY("check-proj");
 771   //proj->Rebin(2);
 772   proj->Scale(1.0 / proj->Integral());
 773
 774   Double_t results[fgMaxParams];
 775   for (Int_t i=0; i<fgMaxParams; ++i)
 776   {
 777     results[i] = inputDist->GetBinContent(i+1);
 778
 779     if (check)
 780       results[i] = proj->GetBinContent(i+1);
 781
 782     // minimum value
 783     results[i] = TMath::Max(results[i], 1e-3);
 784
 785     Float_t stepSize = 0.1;
 786
 787     // minuit sees squared values to prevent it from going negative...
 788     results[i] = TMath::Sqrt(results[i]);
 789     //stepSize /= results[i]; // keep relative step size
 790
 791     minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
 792   }
 793   // bin 0 is filled with value from bin 1 (otherwise it's 0)
 794   //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0);
 795   //results[0] = 0;
 796   //minuit->SetParameter(0, "param0", 0, 0, 0, 0);
 797
 798   for (Int_t i=0; i<fgMaxInput; ++i)
 799   {
 800     //printf("%f %f %f\n", inputDist->GetBinContent(i+1), TMath::Sqrt(inputDist->GetBinContent(i+1)), inputDist->GetBinError(i+1));
 801
 802     (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
 803     if (fCurrentESD->GetBinError(i+1) > 0)
 804       (*fCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / fCurrentESD->GetBinError(i+1) / fCurrentESD->GetBinError(i+1);
 805
 806     if ((*fCorrelationCovarianceMatrix)(i, i) > 1e7)
 807       (*fCorrelationCovarianceMatrix)(i, i) = 0;
 808
 809     //printf("%d --> %e\n", i, (*fCorrelationCovarianceMatrix)(i, i));
 810   }
 811
 812   Int_t dummy;
 813   Double_t chi2;
 814   MinuitFitFunction(dummy, 0, chi2, results, 0);
 815   printf("Chi2 of initial parameters is = %f\n", chi2);
 816
 817   if (check)
 818     return -1;
 819
 820   // first param is number of iterations, second is precision....
 821   arglist[0] = 1e6;
 822   //arglist[1] = 1e-5;
 823   //minuit->ExecuteCommand("SCAN", arglist, 0);
 824   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
 825   printf("MINUIT status is %d\n", status);
 826   //minuit->ExecuteCommand("MIGRAD", arglist, 1);
 827   //minuit->ExecuteCommand("MIGRAD", arglist, 1);
 828   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
 829   //minuit->ExecuteCommand("MINOS", arglist, 0);
 830
 831   for (Int_t i=0; i<fgMaxParams; ++i)
 832   {
 833     if (fCurrentEfficiency->GetBinContent(i+1) > 0)
 834     {
 835       fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) * minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1));
 836       // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
 837       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, minuit->GetParError(i) * minuit->GetParameter(i) /  fCurrentEfficiency->GetBinContent(i+1));
 838     }
 839     else
 840     {
 841       fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, 0);
 842       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
 843     }
 844   }
 845
 846   return status;
 847 }
 848
 849 //____________________________________________________________________
 850 void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
 851 {
 852   //
 853   // correct spectrum using minuit chi2 method applying a NBD fit
 854   //
 855
 856   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 857
 858   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
 859
 860   fCorrelationMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
 861   fCorrelationCovarianceMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
 862   fCurrentESDVector = new TVectorD(fgMaxParams);
 863
 864   // normalize correction for given nPart
 865   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
 866   {
 867     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
 868     if (sum <= 0)
 869       continue;
 870     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
 871     {
 872       // npart sum to 1
 873       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
 874       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
 875
 876       if (i <= fgMaxParams && j <= fgMaxParams)
 877         (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
 878     }
 879   }
 880
 881   for (Int_t i=0; i<fgMaxParams; ++i)
 882   {
 883     (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
 884     if (fCurrentESD->GetBinError(i+1) > 0)
 885       (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
 886   }
 887
 888   // Create NBD function
 889   if (!fNBD)
 890   {
 891     fNBD = 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, 250);
 892     fNBD->SetParNames("scaling", "averagen", "k");
 893   }
 894
 895   // Initialize TMinuit via generic fitter interface
 896   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3);
 897
 898   minuit->SetFCN(MinuitNBD);
 899   minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000);
 900   minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000);
 901   minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000);
 902
 903   Double_t arglist[100];
 904   arglist[0] = 0;
 905   minuit->ExecuteCommand("SET PRINT", arglist, 1);
 906   minuit->ExecuteCommand("MIGRAD", arglist, 0);
 907   //minuit->ExecuteCommand("MINOS", arglist, 0);
 908
 909   fNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
 910
 911   for (Int_t i=0; i<fgMaxParams; ++i)
 912   {
 913     printf("%d : %f\n", i, fNBD->Eval(i));
 914     if (fNBD->Eval(i) > 0)
 915     {
 916       fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fNBD->Eval(i));
 917       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
 918     }
 919   }
 920 }
 921
 922 //____________________________________________________________________
 923 void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
 924 {
 925   //
 926   // normalizes a 1-d histogram to its bin width
 927   //
 928
 929   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
 930   {
 931     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
 932     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
 933   }
 934 }
 935
 936 //____________________________________________________________________
 937 void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
 938 {
 939   //
 940   // normalizes a 2-d histogram to its bin width (x width * y width)
 941   //
 942
 943   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
 944     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
 945     {
 946       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
 947       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
 948       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
 949     }
 950 }
 951
 952 //____________________________________________________________________
 953 void AliMultiplicityCorrection::DrawHistograms()
 954 {
 955   printf("ESD:\n");
 956
 957   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
 958   canvas1->Divide(3, 2);
 959   for (Int_t i = 0; i < kESDHists; ++i)
 960   {
 961     canvas1->cd(i+1);
 962     fMultiplicityESD[i]->DrawCopy("COLZ");
 963     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
 964   }
 965
 966   printf("Vtx:\n");
 967
 968   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
 969   canvas2->Divide(3, 2);
 970   for (Int_t i = 0; i < kMCHists; ++i)
 971   {
 972     canvas2->cd(i+1);
 973     fMultiplicityVtx[i]->DrawCopy("COLZ");
 974     printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
 975   }
 976
 977   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
 978   canvas3->Divide(3, 3);
 979   for (Int_t i = 0; i < kCorrHists; ++i)
 980   {
 981     canvas3->cd(i+1);
 982     TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
 983     for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
 984     {
 985       for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
 986       {
 987         hist->SetBinContent(0, y, z, 0);
 988         hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
 989       }
 990     }
 991     TH1* proj = hist->Project3D("zy");
 992     proj->DrawCopy("COLZ");
 993   }
 994 }
 995
 996 //____________________________________________________________________
 997 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple)
 998 {
 999   //mcHist->Rebin(2);
1000
1001   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1002
1003   TString tmpStr;
1004   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
1005
1006   if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
1007   {
1008     printf("ERROR. Unfolded histogram is empty\n");
1009     return;
1010   }
1011
1012   //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
1013   fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY();
1014   fCurrentESD->Sumw2();
1015   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1016
1017   // normalize unfolded result to 1
1018   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
1019
1020   //fCurrentESD->Scale(mcHist->Integral(2, 200));
1021
1022   //new TCanvas;
1023   /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1024   ratio->Divide(mcHist);
1025   ratio->Draw("HIST");
1026   ratio->Fit("pol0", "W0", "", 20, 120);
1027   Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
1028   delete ratio;
1029   mcHist->Scale(scalingFactor);*/
1030
1031   mcHist->Scale(1.0 / mcHist->Integral());
1032
1033   // calculate residual
1034
1035   // for that we convolute the response matrix with the gathered result
1036   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
1037   TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
1038   tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
1039   TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
1040   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
1041   if (convolutedProj->Integral() > 0)
1042   {
1043     convolutedProj->Scale(1.0 / convolutedProj->Integral());
1044   }
1045   else
1046     printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
1047
1048   //NormalizeToBinWidth(proj2);
1049
1050   TH1* residual = (TH1*) convolutedProj->Clone("residual");
1051   residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
1052
1053   residual->Add(fCurrentESD, -1);
1054   //residual->Divide(residual, fCurrentESD, 1, 1, "B");
1055
1056   TH1* residualHist = new TH1F("residualHist", "residualHist", 50, -5, 5);
1057
1058   // TODO fix errors
1059   Float_t chi2 = 0;
1060   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
1061   {
1062     if (fCurrentESD->GetBinError(i) > 0)
1063     {
1064       Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
1065       if (i > 1)
1066         chi2 += value * value;
1067       residual->SetBinContent(i, value);
1068       residualHist->Fill(value);
1069     }
1070     else
1071     {
1072       //printf("Residual bin %d set to 0\n", i);
1073       residual->SetBinContent(i, 0);
1074     }
1075     convolutedProj->SetBinError(i, 0);
1076     residual->SetBinError(i, 0);
1077
1078     if (i == 200)
1079       fLastChi2Residuals = chi2;
1080   }
1081
1082   residualHist->Fit("gaus", "N");
1083   delete residualHist;
1084
1085   printf("Difference (Residuals) is %f for bin 2-200\n", chi2);
1086
1087   TCanvas* canvas1 = 0;
1088   if (simple)
1089   {
1090     canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
1091     canvas1->Divide(2, 1);
1092   }
1093   else
1094   {
1095     canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
1096     canvas1->Divide(2, 3);
1097   }
1098
1099   canvas1->cd(1);
1100   TH1* proj = (TH1*) mcHist->Clone("proj");
1101   NormalizeToBinWidth(proj);
1102
1103   if (normalizeESD)
1104     NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
1105
1106   proj->GetXaxis()->SetRangeUser(0, 200);
1107   proj->SetTitle(";true multiplicity;Entries");
1108   proj->SetStats(kFALSE);
1109   proj->DrawCopy("HIST");
1110   gPad->SetLogy();
1111
1112   //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
1113   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
1114   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
1115
1116   TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1117   legend->AddEntry(proj, "true distribution");
1118   legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
1119   legend->SetFillColor(0);
1120   legend->Draw();
1121   // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
1122
1123   canvas1->cd(2);
1124
1125   gPad->SetLogy();
1126   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1127   //fCurrentESD->SetLineColor(2);
1128   fCurrentESD->SetTitle(";measured multiplicity;Entries");
1129   fCurrentESD->SetStats(kFALSE);
1130   fCurrentESD->DrawCopy("HIST E");
1131
1132   convolutedProj->SetLineColor(2);
1133   //proj2->SetMarkerColor(2);
1134   //proj2->SetMarkerStyle(5);
1135   convolutedProj->DrawCopy("HIST SAME");
1136
1137   legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1138   legend->AddEntry(fCurrentESD, "measured distribution");
1139   legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
1140   legend->SetFillColor(0);
1141   legend->Draw();
1142
1143   TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
1144   diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1145
1146   /*Float_t chi2 = 0;
1147   Float_t chi = 0;
1148   fLastChi2MCLimit = 0;
1149   Int_t limit = 0;
1150   for (Int_t i=2; i<=200; ++i)
1151   {
1152     if (proj->GetBinContent(i) != 0)
1153     {
1154       Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
1155       chi2 += value * value;
1156       chi += TMath::Abs(value);
1157
1158       //printf("%d %f\n", i, chi);
1159
1160       if (chi2 < 0.2)
1161         fLastChi2MCLimit = i;
1162
1163       if (chi < 3)
1164         limit = i;
1165
1166     }
1167   }*/
1168
1169   chi2 = 0;
1170   Float_t chi = 0;
1171   Int_t limit = 0;
1172   for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
1173   {
1174     if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
1175     {
1176       Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
1177       if (value > 1e8)
1178         value = 1e8; //prevent arithmetic exception
1179       else if (value < -1e8)
1180         value = -1e8;
1181       if (i > 1)
1182       {
1183         chi2 += value * value;
1184         chi += TMath::Abs(value);
1185       }
1186       diffMCUnfolded->SetBinContent(i, value);
1187     }
1188     else
1189     {
1190       //printf("diffMCUnfolded bin %d set to 0\n", i);
1191       diffMCUnfolded->SetBinContent(i, 0);
1192     }
1193     if (chi2 < 1000)
1194       fLastChi2MCLimit = i;
1195     if (chi < 1000)
1196       limit = i;
1197     if (i == 150)
1198       fLastChi2MC = chi2;
1199   }
1200
1201   printf("limits %d %d\n", fLastChi2MCLimit, limit);
1202   fLastChi2MCLimit = limit;
1203
1204   printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);
1205
1206   if (!simple)
1207   {
1208     canvas1->cd(3);
1209
1210     diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e");
1211     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1212     diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
1213     diffMCUnfolded->DrawCopy("HIST");
1214
1215     TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
1216     for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
1217       fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
1218
1219     new TCanvas;
1220     fluctuation->Draw();
1221
1222     /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1223     legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1224     legend->AddEntry(fMultiplicityMC, "MC");
1225     legend->AddEntry(fMultiplicityESD, "ESD");
1226     legend->Draw();*/
1227
1228     canvas1->cd(4);
1229     //residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
1230     residual->GetXaxis()->SetRangeUser(0, 200);
1231     residual->DrawCopy();
1232
1233     canvas1->cd(5);
1234
1235     TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1236     ratio->Divide(mcHist);
1237     ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1238     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1239     ratio->GetXaxis()->SetRangeUser(0, 200);
1240     ratio->SetStats(kFALSE);
1241     ratio->Draw("HIST");
1242
1243     Double_t ratioChi2 = 0;
1244     fLastChi2MCLimit = 0;
1245     for (Int_t i=2; i<=150; ++i)
1246     {
1247       Float_t value = ratio->GetBinContent(i) - 1;
1248       if (value > 1e8)
1249         value = 1e8; //prevent arithmetic exception
1250       else if (value < -1e8)
1251         value = -1e8;
1252
1253       ratioChi2 += value * value;
1254
1255       if (ratioChi2 < 0.1)
1256         fLastChi2MCLimit = i;
1257     }
1258
1259     printf("Sum over (ratio-1)^2 (2..150) is %f\n", ratioChi2);
1260     // TODO FAKE
1261     fLastChi2MC = ratioChi2;
1262   }
1263
1264   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
1265 }
1266
1267 //____________________________________________________________________
1268 void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals)
1269 {
1270   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1271   // These values are computed during DrawComparison, thus this function picks up the
1272   // last calculation
1273
1274   if (mc)
1275     *mc = fLastChi2MC;
1276   if (mcLimit)
1277     *mcLimit = fLastChi2MCLimit;
1278   if (residuals)
1279     *residuals = fLastChi2Residuals;
1280 }
1281
1282
1283 //____________________________________________________________________
1284 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1285 {
1286   //
1287   // returns the corresponding MC spectrum
1288   //
1289
1290   switch (eventType)
1291   {
1292     case kTrVtx : return fMultiplicityVtx[i]; break;
1293     case kMB : return fMultiplicityMB[i]; break;
1294     case kINEL : return fMultiplicityINEL[i]; break;
1295   }
1296
1297   return 0;
1298 }
1299
1300 //____________________________________________________________________
1301 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* inputDist)
1302 {
1303   //
1304   // correct spectrum using bayesian method
1305   //
1306
1307   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1308
1309   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
1310
1311   // smooth efficiency
1312   //TH1* tmp = (TH1*) fCurrentEfficiency->Clone("eff_clone");
1313   //for (Int_t i=2; i<fCurrentEfficiency->GetNbinsX(); ++i)
1314   //  fCurrentEfficiency->SetBinContent(i, (tmp->GetBinContent(i-1) + tmp->GetBinContent(i) + tmp->GetBinContent(i+1)) / 3);
1315
1316   // set efficiency to 1 above 150
1317   // FAKE TEST
1318   //for (Int_t i=150; i<=fCurrentEfficiency->GetNbinsX(); ++i)
1319   //  fCurrentEfficiency->SetBinContent(i, 1);
1320
1321   // normalize correction for given nPart
1322   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1323   {
1324     // with this it is normalized to 1
1325     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1326
1327     // with this normalized to the given efficiency
1328     if (fCurrentEfficiency->GetBinContent(i) > 0)
1329       sum /= fCurrentEfficiency->GetBinContent(i);
1330     else
1331       sum = 0;
1332
1333     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1334     {
1335       if (sum > 0)
1336       {
1337         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1338         fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1339       }
1340       else
1341       {
1342         fCurrentCorrelation->SetBinContent(i, j, 0);
1343         fCurrentCorrelation->SetBinError(i, j, 0);
1344       }
1345     }
1346   }
1347
1348   // normalize nTrack
1349   /*for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1350   {
1351     // with this it is normalized to 1
1352     Double_t sum = fCurrentCorrelation->Integral(1, fCurrentCorrelation->GetNbinsX(), j, j);
1353
1354     for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1355     {
1356       if (sum > 0)
1357       {
1358         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1359         fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1360       }
1361       else
1362       {
1363         fCurrentCorrelation->SetBinContent(i, j, 0);
1364         fCurrentCorrelation->SetBinError(i, j, 0);
1365       }
1366     }
1367   }*/
1368
1369   // smooth input spectrum
1370   /*
1371   TH1* esdClone = (TH1*) fCurrentESD->Clone("esdClone");
1372
1373   for (Int_t i=2; i<fCurrentESD->GetNbinsX(); ++i)
1374     if (esdClone->GetBinContent(i) < 1e-3)
1375       fCurrentESD->SetBinContent(i, (esdClone->GetBinContent(i-1) + esdClone->GetBinContent(i) + esdClone->GetBinContent(i+1)) / 3);
1376
1377   delete esdClone;
1378
1379   // rescale to 1
1380   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1381   */
1382
1383   /*new TCanvas;
1384   fCurrentEfficiency->Draw();
1385
1386   new TCanvas;
1387   fCurrentCorrelation->DrawCopy("COLZ");
1388
1389   new TCanvas;
1390   ((TH2*) fCurrentCorrelation)->ProjectionX()->DrawCopy();
1391
1392   new TCanvas;
1393   ((TH2*) fCurrentCorrelation)->ProjectionY()->DrawCopy();*/
1394
1395   // pick prior distribution
1396   TH1* hPrior = 0;
1397   if (inputDist)
1398   {
1399     printf("Using different starting conditions...\n");
1400     hPrior = (TH1*)inputDist->Clone("prior");
1401   }
1402   else
1403     hPrior = (TH1*)fCurrentESD->Clone("prior");
1404   Float_t norm = hPrior->Integral();
1405   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1406     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1407
1408   // define temp hist
1409   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1410   hTemp->Reset();
1411
1412   // just a shortcut
1413   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
1414
1415   // histogram to store the inverse response calculated from bayes theorem (from prior and response) IR
1416   TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji");
1417   hInverseResponseBayes->Reset();
1418
1419   TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
1420
1421   const Int_t kStartBin = 1;
1422
1423   // unfold...
1424   for (Int_t i=0; i<nIterations; i++)
1425   {
1426     //printf(" iteration %i \n", i);
1427
1428     // calculate IR from Beyes theorem
1429     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
1430
1431     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1432     {
1433       Float_t norm = 0;
1434       for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
1435         norm += hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t);
1436       for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
1437       {
1438         if (norm==0)
1439           hInverseResponseBayes->SetBinContent(t,m,0);
1440         else
1441           hInverseResponseBayes->SetBinContent(t,m, hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t)/norm);
1442       }
1443     }
1444
1445     for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
1446     {
1447       Float_t value = 0;
1448       for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1449         value += fCurrentESD->GetBinContent(m) * hInverseResponseBayes->GetBinContent(t,m);
1450
1451       if (fCurrentEfficiency->GetBinContent(t) > 0)
1452         hTemp->SetBinContent(t, value / fCurrentEfficiency->GetBinContent(t));
1453       else
1454         hTemp->SetBinContent(t, 0);
1455     }
1456
1457     // this is the last guess, fill before (!) smoothing
1458     for (Int_t t=kStartBin; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1459     {
1460       //as bin error put the difference to the last iteration
1461       //fMultiplicityESDCorrected[correlationID]->SetBinError(t, hTemp->GetBinContent(t) - fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
1462       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1463       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t));
1464
1465       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
1466     }
1467
1468     /*new TCanvas;
1469     fMultiplicityESDCorrected[correlationID]->DrawCopy();
1470     gPad->SetLogy();*/
1471
1472     // regularization (simple smoothing)
1473     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
1474
1475     for (Int_t t=kStartBin+2; t<hTrueSmoothed->GetNbinsX(); t++)
1476     {
1477       Float_t average = (hTemp->GetBinContent(t-1)
1478                          + hTemp->GetBinContent(t)
1479                          + hTemp->GetBinContent(t+1)
1480                          ) / 3.;
1481
1482       // weight the average with the regularization parameter
1483       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
1484     }
1485
1486     // calculate chi2 (change from last iteration)
1487     Double_t chi2 = 0;
1488
1489     // use smoothed true (normalized) as new prior
1490     Float_t norm = 1;
1491     //for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1492     //  norm = norm + hTrueSmoothed->GetBinContent(t);
1493
1494     for (Int_t t=kStartBin; t<=hTrueSmoothed->GetNbinsX(); t++)
1495     {
1496       Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm;
1497       Float_t diff = hPrior->GetBinContent(t) - newValue;
1498       chi2 += (Double_t) diff * diff;
1499
1500       hPrior->SetBinContent(t, newValue);
1501     }
1502
1503     //printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1504
1505     convergence->Fill(i+1, chi2);
1506
1507     //if (i % 5 == 0)
1508     //  DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
1509
1510     delete hTrueSmoothed;
1511   } // end of iterations
1512
1513   //new TCanvas;
1514   //convergence->DrawCopy();
1515   //gPad->SetLogy();
1516
1517   delete convergence;
1518
1519   return;
1520
1521   // ********
1522   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
1523
1524   /*printf("Calculating covariance matrix. This may take some time...\n");
1525
1526   // TODO check if this is the right one...
1527   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
1528
1529   Int_t xBins = hInverseResponseBayes->GetNbinsX();
1530   Int_t yBins = hInverseResponseBayes->GetNbinsY();
1531
1532   // calculate "unfolding matrix" Mij
1533   Float_t matrixM[251][251];
1534   for (Int_t i=1; i<=xBins; i++)
1535   {
1536     for (Int_t j=1; j<=yBins; j++)
1537     {
1538       if (fCurrentEfficiency->GetBinContent(i) > 0)
1539         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
1540       else
1541         matrixM[i-1][j-1] = 0;
1542     }
1543   }
1544
1545   Float_t* vectorn = new Float_t[yBins];
1546   for (Int_t j=1; j<=yBins; j++)
1547     vectorn[j-1] = fCurrentESD->GetBinContent(j);
1548
1549   // first part of covariance matrix, depends on input distribution n(E)
1550   Float_t cov1[251][251];
1551
1552   Float_t nEvents = fCurrentESD->Integral(); // N
1553
1554   xBins = 20;
1555   yBins = 20;
1556
1557   for (Int_t k=0; k<xBins; k++)
1558   {
1559     printf("In Cov1: %d\n", k);
1560     for (Int_t l=0; l<yBins; l++)
1561     {
1562       cov1[k][l] = 0;
1563
1564       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
1565       for (Int_t j=0; j<yBins; j++)
1566         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
1567           * (1.0 - vectorn[j] / nEvents);
1568
1569       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
1570       for (Int_t i=0; i<yBins; i++)
1571         for (Int_t j=0; j<yBins; j++)
1572         {
1573           if (i == j)
1574             continue;
1575           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
1576             * vectorn[j] / nEvents;
1577          }
1578     }
1579   }
1580
1581   printf("Cov1 finished\n");
1582
1583   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
1584   cov->Reset();
1585
1586   for (Int_t i=1; i<=xBins; i++)
1587     for (Int_t j=1; j<=yBins; j++)
1588       cov->SetBinContent(i, j, cov1[i-1][j-1]);
1589
1590   new TCanvas;
1591   cov->Draw("COLZ");
1592
1593   // second part of covariance matrix, depends on response matrix
1594   Float_t cov2[251][251];
1595
1596   // Cov[P(Er|Cu), P(Es|Cu)] term
1597   Float_t covTerm[100][100][100];
1598   for (Int_t r=0; r<yBins; r++)
1599     for (Int_t u=0; u<xBins; u++)
1600       for (Int_t s=0; s<yBins; s++)
1601       {
1602         if (r == s)
1603           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
1604             * (1.0 - hResponse->GetBinContent(u+1, r+1));
1605         else
1606           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
1607             * hResponse->GetBinContent(u+1, s+1);
1608       }
1609
1610   for (Int_t k=0; k<xBins; k++)
1611     for (Int_t l=0; l<yBins; l++)
1612     {
1613       cov2[k][l] = 0;
1614       printf("In Cov2: %d %d\n", k, l);
1615       for (Int_t i=0; i<yBins; i++)
1616         for (Int_t j=0; j<yBins; j++)
1617         {
1618           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
1619           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
1620           Float_t tmpCov = 0;
1621           for (Int_t r=0; r<yBins; r++)
1622             for (Int_t u=0; u<xBins; u++)
1623               for (Int_t s=0; s<yBins; s++)
1624               {
1625                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
1626                   || hResponse->GetBinContent(u+1, i+1) == 0)
1627                   continue;
1628
1629                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
1630                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
1631                   * covTerm[r][u][s];
1632               }
1633
1634           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
1635         }
1636     }
1637
1638   printf("Cov2 finished\n");
1639
1640   for (Int_t i=1; i<=xBins; i++)
1641     for (Int_t j=1; j<=yBins; j++)
1642       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
1643
1644   new TCanvas;
1645   cov->Draw("COLZ");*/
1646 }
1647
1648 //____________________________________________________________________
1649 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse,
1650   TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u)
1651 {
1652   //
1653   // helper function for the covariance matrix of the bayesian method
1654   //
1655
1656   Float_t result = 0;
1657
1658   if (k == u && r == i)
1659     result += 1.0 / hResponse->GetBinContent(u+1, r+1);
1660
1661   if (k == u)
1662     result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
1663
1664   if (r == i)
1665     result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
1666
1667   result *= matrixM[k][i];
1668
1669   return result;
1670 }
1671
1672 //____________________________________________________________________
1673 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
1674 {
1675   //
1676   // correct spectrum using bayesian method
1677   //
1678
1679   Float_t regPar = 0;
1680
1681   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1682   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1683
1684   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
1685
1686   // TODO should be taken from correlation map
1687   //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
1688
1689   // normalize correction for given nPart
1690   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1691   {
1692     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1693     //Double_t sum = sumHist->GetBinContent(i);
1694     if (sum <= 0)
1695       continue;
1696     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1697     {
1698       // npart sum to 1
1699       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1700       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1701     }
1702   }
1703
1704   new TCanvas;
1705   fCurrentCorrelation->Draw("COLZ");
1706
1707   // FAKE
1708   fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
1709
1710   // pick prior distribution
1711   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
1712   Float_t norm = 1; //hPrior->Integral();
1713   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1714     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1715
1716   // zero distribution
1717   TH1F* zero =  (TH1F*)hPrior->Clone("zero");
1718
1719   // define temp hist
1720   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1721   hTemp->Reset();
1722
1723   // just a shortcut
1724   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
1725
1726   // unfold...
1727   Int_t iterations = 25;
1728   for (Int_t i=0; i<iterations; i++)
1729   {
1730     //printf(" iteration %i \n", i);
1731
1732     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1733     {
1734       Float_t value = 0;
1735       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1736         value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
1737       hTemp->SetBinContent(m, value);
1738       //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
1739     }
1740
1741     // regularization (simple smoothing)
1742     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
1743
1744     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
1745     {
1746       Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
1747                          + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
1748                          + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
1749       average *= hTrueSmoothed->GetBinWidth(t);
1750
1751       // weight the average with the regularization parameter
1752       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
1753     }
1754
1755     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1756       hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
1757
1758     // fill guess
1759     for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1760     {
1761       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1762       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
1763
1764       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
1765     }
1766
1767
1768     // calculate chi2 (change from last iteration)
1769     Double_t chi2 = 0;
1770
1771     // use smoothed true (normalized) as new prior
1772     Float_t norm = 1; //hTrueSmoothed->Integral();
1773
1774     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
1775     {
1776       Float_t newValue = hTemp->GetBinContent(t)/norm;
1777       Float_t diff = hPrior->GetBinContent(t) - newValue;
1778       chi2 += (Double_t) diff * diff;
1779
1780       hPrior->SetBinContent(t, newValue);
1781     }
1782
1783     printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1784
1785     //if (i % 5 == 0)
1786       DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1787
1788     delete hTrueSmoothed;
1789   } // end of iterations
1790
1791   DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1792 }
1793
1794 //____________________________________________________________________
1795 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
1796 {
1797   //
1798   // correct spectrum using a simple Gaussian approach, that is model-dependent
1799   //
1800
1801   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1802   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1803
1804   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
1805
1806   NormalizeToBinWidth((TH2*) fCurrentCorrelation);
1807
1808   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
1809   correction->SetTitle("GaussianMean");
1810
1811   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
1812   correction->SetTitle("GaussianWidth");
1813
1814   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1815   {
1816     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
1817     proj->Fit("gaus", "0Q");
1818     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1819     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1820
1821     continue;
1822
1823     // draw for debugging
1824     new TCanvas;
1825     proj->DrawCopy();
1826     proj->GetFunction("gaus")->DrawCopy("SAME");
1827   }
1828
1829   TH1* target = fMultiplicityESDCorrected[correlationID];
1830
1831   Int_t nBins = target->GetNbinsX()*10+1;
1832   Float_t* binning = new Float_t[nBins];
1833   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1834     for (Int_t j=0; j<10; ++j)
1835       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
1836
1837   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
1838
1839   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
1840
1841   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1842   {
1843     Float_t mean = correction->GetBinContent(i);
1844     Float_t width = correctionWidth->GetBinContent(i);
1845
1846     Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
1847     Int_t fillEnd   = fineBinned->FindBin(mean + width * 5);
1848     //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1849
1850     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1851     {
1852       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
1853     }
1854   }
1855
1856   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1857   {
1858     Float_t sum = 0;
1859     for (Int_t j=1; j<=10; ++j)
1860       sum += fineBinned->GetBinContent((i-1)*10 + j);
1861     target->SetBinContent(i, sum / 10);
1862   }
1863
1864   delete[] binning;
1865
1866   DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
1867 }
1868
1869 //____________________________________________________________________
1870 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
1871 {
1872   // runs the distribution given in inputMC through the response matrix identified by
1873   // correlationMap and produces a measured distribution
1874   // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
1875   // if normalized is set, inputMC is expected to be normalized to the bin width
1876
1877   if (!inputMC)
1878     return 0;
1879
1880   if (correlationMap < 0 || correlationMap >= kCorrHists)
1881     return 0;
1882
1883   // empty under/overflow bins in x, otherwise Project3D takes them into account
1884   TH3* hist = fCorrelation[correlationMap];
1885   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1886   {
1887     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1888     {
1889       hist->SetBinContent(0, y, z, 0);
1890       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1891     }
1892   }
1893
1894   TH2* corr = (TH2*) hist->Project3D("zy");
1895   //corr->Rebin2D(2, 1);
1896   corr->Sumw2();
1897
1898   // normalize correction for given nPart
1899   for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
1900   {
1901     Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
1902     if (sum <= 0)
1903       continue;
1904
1905     for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
1906     {
1907       // npart sum to 1
1908       corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
1909       corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
1910     }
1911   }
1912
1913   TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
1914   target->Reset();
1915
1916   for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
1917   {
1918     Float_t measured = 0;
1919     Float_t error = 0;
1920
1921     for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
1922     {
1923       Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
1924
1925       measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
1926       error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
1927     }
1928
1929     //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
1930
1931     target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
1932     target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
1933   }
1934
1935   return target;
1936 }
1937
1938 //____________________________________________________________________
1939 void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
1940 {
1941   // uses the given function to fill the input MC histogram and generates from that
1942   // the measured histogram by applying the response matrix
1943   // this can be used to evaluate if the methods work indepedently of the input
1944   // distribution
1945   // WARNING does not respect the vertex distribution, just fills central vertex bin
1946
1947   if (!inputMC)
1948     return;
1949
1950   if (id < 0 || id >= kESDHists)
1951     return;
1952
1953   TH2F* mc = fMultiplicityVtx[id];
1954
1955   mc->Reset();
1956
1957   for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
1958   {
1959     mc->SetBinContent(mc->GetNbinsX() / 2, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
1960     mc->SetBinError(mc->GetNbinsX() / 2, i, 0);
1961   }
1962
1963   //new TCanvas;
1964   //mc->Draw("COLZ");
1965
1966   TH1* proj = mc->ProjectionY();
1967
1968   TString tmp(fMultiplicityESD[id]->GetName());
1969   delete fMultiplicityESD[id];
1970   fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
1971   fMultiplicityESD[id]->SetName(tmp);
1972 }