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 // e.g. chi2 minimization and bayesian unfolding
  22 //
  23 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
  24
  25 #include "AliMultiplicityCorrection.h"
  26
  27 #include <TFile.h>
  28 #include <TH1F.h>
  29 #include <TH2F.h>
  30 #include <TH3F.h>
  31 #include <TDirectory.h>
  32 #include <TVirtualFitter.h>
  33 #include <TCanvas.h>
  34 #include <TString.h>
  35 #include <TF1.h>
  36 #include <TMath.h>
  37 #include <TCollection.h>
  38 #include <TLegend.h>
  39 #include <TLine.h>
  40 #include <TRandom.h>
  41 #include <TProfile.h>
  42 #include <TProfile2D.h>
  43
  44 ClassImp(AliMultiplicityCorrection)
  45
  46 const Int_t AliMultiplicityCorrection::fgkMaxInput  = 250;  // bins in measured histogram
  47 const Int_t AliMultiplicityCorrection::fgkMaxParams = 250;  // bins in unfolded histogram = number of fit params
  48
  49 TMatrixD* AliMultiplicityCorrection::fgCorrelationMatrix = 0;
  50 TMatrixD* AliMultiplicityCorrection::fgCorrelationCovarianceMatrix = 0;
  51 TVectorD* AliMultiplicityCorrection::fgCurrentESDVector = 0;
  52 TVectorD* AliMultiplicityCorrection::fgEntropyAPriori = 0;
  53
  54 AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fgRegularizationType = AliMultiplicityCorrection::kPol1;
  55 Float_t AliMultiplicityCorrection::fgRegularizationWeight = 10000;
  56
  57 TF1* AliMultiplicityCorrection::fgNBD = 0;
  58
  59 Float_t AliMultiplicityCorrection::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
  60 Int_t   AliMultiplicityCorrection::fgBayesianIterations = 100;         // number of iterations in Bayesian method
  61
  62 // These are the areas where the quality of the unfolding results are evaluated
  63 // Default defined here, call SetQualityRegions to change them
  64 // unit is in multiplicity (not in bin!)
  65
  66 // SPD:   peak area - flat area - low stat area
  67 Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {4,  60,  180};
  68 Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {20, 140, 210};
  69
  70 //____________________________________________________________________
  71 void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
  72 {
  73   //
  74   // sets the quality region definition to TPC or SPD
  75   //
  76
  77   if (SPDStudy)
  78   {
  79     // SPD:   peak area - flat area - low stat area
  80     fgQualityRegionsB[0] = 4;
  81     fgQualityRegionsE[0] = 20;
  82
  83     fgQualityRegionsB[1] = 60;
  84     fgQualityRegionsE[1] = 140;
  85
  86     fgQualityRegionsB[2] = 180;
  87     fgQualityRegionsE[2] = 210;
  88
  89     Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD");
  90   }
  91   else
  92   {
  93     // TPC:   peak area - flat area - low stat area
  94     fgQualityRegionsB[0] = 4;
  95     fgQualityRegionsE[0] = 12;
  96
  97     fgQualityRegionsB[1] = 25;
  98     fgQualityRegionsE[1] = 55;
  99
 100     fgQualityRegionsB[2] = 88;
 101     fgQualityRegionsE[2] = 108;
 102
 103     Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for TPC");
 104   }
 105 }
 106
 107 //____________________________________________________________________
 108 AliMultiplicityCorrection::AliMultiplicityCorrection() :
 109   TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
 110 {
 111   //
 112   // default constructor
 113   //
 114
 115   for (Int_t i = 0; i < kESDHists; ++i)
 116     fMultiplicityESD[i] = 0;
 117
 118   for (Int_t i = 0; i < kMCHists; ++i)
 119   {
 120     fMultiplicityVtx[i] = 0;
 121     fMultiplicityMB[i] = 0;
 122     fMultiplicityINEL[i] = 0;
 123   }
 124
 125   for (Int_t i = 0; i < kCorrHists; ++i)
 126   {
 127     fCorrelation[i] = 0;
 128     fMultiplicityESDCorrected[i] = 0;
 129   }
 130
 131   for (Int_t i = 0; i < kQualityRegions; ++i)
 132     fQuality[i] = 0;
 133 }
 134
 135 //____________________________________________________________________
 136 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
 137   TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
 138 {
 139   //
 140   // named constructor
 141   //
 142
 143   // do not add this hists to the directory
 144   Bool_t oldStatus = TH1::AddDirectoryStatus();
 145   TH1::AddDirectory(kFALSE);
 146
 147   /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
 148   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,
 149                           10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
 150                           20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
 151                           30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
 152                           40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
 153                           50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
 154                           100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
 155                           150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
 156                           250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
 157                           500.5 };
 158                                 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
 159                           //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
 160                           //1000.5 };
 161
 162   #define VTXBINNING 10, binLimitsVtx
 163   #define NBINNING fgkMaxParams, binLimitsN*/
 164
 165   #define NBINNING 501, -0.5, 500.5
 166   #define VTXBINNING 1, -10, 10
 167
 168   for (Int_t i = 0; i < kESDHists; ++i)
 169     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
 170
 171   for (Int_t i = 0; i < kMCHists; ++i)
 172   {
 173     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
 174     fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
 175
 176     fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
 177     fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
 178
 179     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
 180     fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
 181   }
 182
 183   for (Int_t i = 0; i < kCorrHists; ++i)
 184   {
 185     fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
 186     fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
 187   }
 188
 189   TH1::AddDirectory(oldStatus);
 190 }
 191
 192 //____________________________________________________________________
 193 AliMultiplicityCorrection::~AliMultiplicityCorrection()
 194 {
 195   //
 196   // Destructor
 197   //
 198
 199   for (Int_t i = 0; i < kESDHists; ++i)
 200   {
 201     if (fMultiplicityESD[i])
 202       delete fMultiplicityESD[i];
 203     fMultiplicityESD[i] = 0;
 204   }
 205
 206   for (Int_t i = 0; i < kMCHists; ++i)
 207   {
 208     if (fMultiplicityVtx[i])
 209       delete fMultiplicityVtx[i];
 210     fMultiplicityVtx[i] = 0;
 211
 212     if (fMultiplicityMB[i])
 213       delete fMultiplicityMB[i];
 214     fMultiplicityMB[i] = 0;
 215
 216     if (fMultiplicityINEL[i])
 217       delete fMultiplicityINEL[i];
 218     fMultiplicityINEL[i] = 0;
 219   }
 220
 221   for (Int_t i = 0; i < kCorrHists; ++i)
 222   {
 223     if (fCorrelation[i])
 224       delete fCorrelation[i];
 225     fCorrelation[i] = 0;
 226
 227     if (fMultiplicityESDCorrected[i])
 228       delete fMultiplicityESDCorrected[i];
 229     fMultiplicityESDCorrected[i] = 0;
 230   }
 231 }
 232
 233 //____________________________________________________________________
 234 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
 235 {
 236   // Merge a list of AliMultiplicityCorrection objects with this (needed for
 237   // PROOF).
 238   // Returns the number of merged objects (including this).
 239
 240   if (!list)
 241     return 0;
 242
 243   if (list->IsEmpty())
 244     return 1;
 245
 246   TIterator* iter = list->MakeIterator();
 247   TObject* obj;
 248
 249   // collections of all histograms
 250   TList collections[kESDHists+kMCHists*3+kCorrHists*2];
 251
 252   Int_t count = 0;
 253   while ((obj = iter->Next())) {
 254
 255     AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
 256     if (entry == 0)
 257       continue;
 258
 259     for (Int_t i = 0; i < kESDHists; ++i)
 260       collections[i].Add(entry->fMultiplicityESD[i]);
 261
 262     for (Int_t i = 0; i < kMCHists; ++i)
 263     {
 264       collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
 265       collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
 266       collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
 267     }
 268
 269     for (Int_t i = 0; i < kCorrHists; ++i)
 270       collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]);
 271
 272     for (Int_t i = 0; i < kCorrHists; ++i)
 273       collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
 274
 275     count++;
 276   }
 277
 278   for (Int_t i = 0; i < kESDHists; ++i)
 279     fMultiplicityESD[i]->Merge(&collections[i]);
 280
 281   for (Int_t i = 0; i < kMCHists; ++i)
 282   {
 283     fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
 284     fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
 285     fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
 286   }
 287
 288   for (Int_t i = 0; i < kCorrHists; ++i)
 289     fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
 290
 291   for (Int_t i = 0; i < kCorrHists; ++i)
 292     fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]);
 293
 294   delete iter;
 295
 296   return count+1;
 297 }
 298
 299 //____________________________________________________________________
 300 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
 301 {
 302   //
 303   // loads the histograms from a file
 304   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
 305   //
 306
 307   if (!dir)
 308     dir = GetName();
 309
 310   if (!gDirectory->cd(dir))
 311     return kFALSE;
 312
 313   // TODO memory leak. old histograms needs to be deleted.
 314
 315   Bool_t success = kTRUE;
 316
 317   for (Int_t i = 0; i < kESDHists; ++i)
 318   {
 319     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
 320     if (!fMultiplicityESD[i])
 321       success = kFALSE;
 322   }
 323
 324   for (Int_t i = 0; i < kMCHists; ++i)
 325   {
 326     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
 327     fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
 328     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
 329     if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
 330       success = kFALSE;
 331   }
 332
 333   for (Int_t i = 0; i < kCorrHists; ++i)
 334   {
 335     fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
 336     if (!fCorrelation[i])
 337       success = kFALSE;
 338     fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
 339     if (!fMultiplicityESDCorrected[i])
 340       success = kFALSE;
 341   }
 342
 343   gDirectory->cd("..");
 344
 345   return success;
 346 }
 347
 348 //____________________________________________________________________
 349 void AliMultiplicityCorrection::SaveHistograms()
 350 {
 351   //
 352   // saves the histograms
 353   //
 354
 355   gDirectory->mkdir(GetName());
 356   gDirectory->cd(GetName());
 357
 358   for (Int_t i = 0; i < kESDHists; ++i)
 359     if (fMultiplicityESD[i])
 360       fMultiplicityESD[i]->Write();
 361
 362   for (Int_t i = 0; i < kMCHists; ++i)
 363   {
 364     if (fMultiplicityVtx[i])
 365       fMultiplicityVtx[i]->Write();
 366     if (fMultiplicityMB[i])
 367       fMultiplicityMB[i]->Write();
 368     if (fMultiplicityINEL[i])
 369       fMultiplicityINEL[i]->Write();
 370   }
 371
 372   for (Int_t i = 0; i < kCorrHists; ++i)
 373   {
 374     if (fCorrelation[i])
 375       fCorrelation[i]->Write();
 376     if (fMultiplicityESDCorrected[i])
 377       fMultiplicityESDCorrected[i]->Write();
 378   }
 379
 380   gDirectory->cd("..");
 381 }
 382
 383 //____________________________________________________________________
 384 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)
 385 {
 386   //
 387   // Fills an event from MC
 388   //
 389
 390   if (triggered)
 391   {
 392     fMultiplicityMB[0]->Fill(vtx, generated05);
 393     fMultiplicityMB[1]->Fill(vtx, generated10);
 394     fMultiplicityMB[2]->Fill(vtx, generated15);
 395     fMultiplicityMB[3]->Fill(vtx, generated20);
 396     fMultiplicityMB[4]->Fill(vtx, generatedAll);
 397
 398     if (vertex)
 399     {
 400       fMultiplicityVtx[0]->Fill(vtx, generated05);
 401       fMultiplicityVtx[1]->Fill(vtx, generated10);
 402       fMultiplicityVtx[2]->Fill(vtx, generated15);
 403       fMultiplicityVtx[3]->Fill(vtx, generated20);
 404       fMultiplicityVtx[4]->Fill(vtx, generatedAll);
 405     }
 406   }
 407
 408   fMultiplicityINEL[0]->Fill(vtx, generated05);
 409   fMultiplicityINEL[1]->Fill(vtx, generated10);
 410   fMultiplicityINEL[2]->Fill(vtx, generated15);
 411   fMultiplicityINEL[3]->Fill(vtx, generated20);
 412   fMultiplicityINEL[4]->Fill(vtx, generatedAll);
 413 }
 414
 415 //____________________________________________________________________
 416 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
 417 {
 418   //
 419   // Fills an event from ESD
 420   //
 421
 422   fMultiplicityESD[0]->Fill(vtx, measured05);
 423   fMultiplicityESD[1]->Fill(vtx, measured10);
 424   fMultiplicityESD[2]->Fill(vtx, measured15);
 425   fMultiplicityESD[3]->Fill(vtx, measured20);
 426 }
 427
 428 //____________________________________________________________________
 429 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)
 430 {
 431   //
 432   // Fills an event into the correlation map with the information from MC and ESD
 433   //
 434
 435   fCorrelation[0]->Fill(vtx, generated05, measured05);
 436   fCorrelation[1]->Fill(vtx, generated10, measured10);
 437   fCorrelation[2]->Fill(vtx, generated15, measured15);
 438   fCorrelation[3]->Fill(vtx, generated20, measured20);
 439
 440   fCorrelation[4]->Fill(vtx, generatedAll, measured05);
 441   fCorrelation[5]->Fill(vtx, generatedAll, measured10);
 442   fCorrelation[6]->Fill(vtx, generatedAll, measured15);
 443   fCorrelation[7]->Fill(vtx, generatedAll, measured20);
 444 }
 445
 446 //____________________________________________________________________
 447 Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
 448 {
 449   // homogenity term for minuit fitting
 450   // pure function of the parameters
 451   // prefers constant function (pol0)
 452
 453   Double_t chi2 = 0;
 454
 455   // ignore the first bin here. on purpose...
 456   for (Int_t i=2; i<fgkMaxParams; ++i)
 457   {
 458     Double_t right  = params[i];
 459     Double_t left   = params[i-1];
 460
 461     if (right != 0)
 462     {
 463       Double_t diff = 1 - left / right;
 464       chi2 += diff * diff;
 465     }
 466   }
 467
 468   return chi2 / 100.0;
 469 }
 470
 471 //____________________________________________________________________
 472 Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
 473 {
 474   // homogenity term for minuit fitting
 475   // pure function of the parameters
 476   // prefers linear function (pol1)
 477
 478   Double_t chi2 = 0;
 479
 480   // ignore the first bin here. on purpose...
 481   for (Int_t i=2+1; i<fgkMaxParams; ++i)
 482   {
 483     if (params[i-1] == 0)
 484       continue;
 485
 486     Double_t right  = params[i];
 487     Double_t middle = params[i-1];
 488     Double_t left   = params[i-2];
 489
 490     Double_t der1 = (right - middle);
 491     Double_t der2 = (middle - left);
 492
 493     Double_t diff = (der1 - der2) / middle;
 494
 495     chi2 += diff * diff;
 496   }
 497
 498   return chi2;
 499 }
 500
 501 //____________________________________________________________________
 502 Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
 503 {
 504   // homogenity term for minuit fitting
 505   // pure function of the parameters
 506   // prefers linear function (pol1)
 507
 508   Double_t chi2 = 0;
 509
 510   /*Float_t der[fgkMaxParams];
 511
 512   for (Int_t i=0; i<fgkMaxParams-1; ++i)
 513     der[i] = params[i+1] - params[i];
 514
 515   for (Int_t i=0; i<fgkMaxParams-6; ++i)
 516   {
 517     for (Int_t j=1; j<=5; ++j)
 518     {
 519       Double_t diff = der[i+j] - der[i];
 520       chi2 += diff * diff;
 521     }
 522   }*/
 523
 524   // ignore the first bin here. on purpose...
 525   for (Int_t i=2+1; i<fgkMaxParams; ++i)
 526   {
 527     if (params[i-1] == 0)
 528       continue;
 529
 530     Double_t right  = log(params[i]);
 531     Double_t middle = log(params[i-1]);
 532     Double_t left   = log(params[i-2]);
 533
 534     Double_t der1 = (right - middle);
 535     Double_t der2 = (middle - left);
 536
 537     Double_t diff = (der1 - der2) / middle;
 538
 539     chi2 += diff * diff;
 540   }
 541
 542   return chi2 * 100;
 543 }
 544
 545 //____________________________________________________________________
 546 Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
 547 {
 548   // homogenity term for minuit fitting
 549   // pure function of the parameters
 550   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
 551   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
 552
 553   Double_t chi2 = 0;
 554
 555   // ignore the first bin here. on purpose...
 556   for (Int_t i=3; i<fgkMaxParams; ++i)
 557   {
 558     Double_t right  = params[i];
 559     Double_t middle = params[i-1];
 560     Double_t left   = params[i-2];
 561
 562     Double_t der1 = (right - middle);
 563     Double_t der2 = (middle - left);
 564
 565     Double_t diff = (der1 - der2);
 566
 567     chi2 += diff * diff;
 568   }
 569
 570   return chi2 * 1e4;
 571 }
 572
 573 //____________________________________________________________________
 574 Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
 575 {
 576   // homogenity term for minuit fitting
 577   // pure function of the parameters
 578   // calculates entropy, from
 579   // The method of reduced cross-entropy (M. Schmelling 1993)
 580
 581   Double_t paramSum = 0;
 582   // ignore the first bin here. on purpose...
 583   for (Int_t i=1; i<fgkMaxParams; ++i)
 584     paramSum += params[i];
 585
 586   Double_t chi2 = 0;
 587   for (Int_t i=1; i<fgkMaxParams; ++i)
 588   {
 589     //Double_t tmp = params[i] / paramSum;
 590     Double_t tmp = params[i];
 591     if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
 592       chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
 593   }
 594
 595   return 10.0 + chi2;
 596 }
 597
 598 //____________________________________________________________________
 599 void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
 600 {
 601   //
 602   // fit function for minuit
 603   // does: nbd
 604   // 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);
 605   // func->SetParNames("scaling", "averagen", "k");
 606   // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
 607   // func->SetParLimits(1, 0.001, 1000);
 608   // func->SetParLimits(2, 0.001, 1000);
 609   //
 610
 611   fgNBD->SetParameters(params[0], params[1], params[2]);
 612
 613   Double_t params2[fgkMaxParams];
 614
 615   for (Int_t i=0; i<fgkMaxParams; ++i)
 616     params2[i] = fgNBD->Eval(i);
 617
 618   MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
 619
 620   printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2);
 621 }
 622
 623 //____________________________________________________________________
 624 void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
 625 {
 626   //
 627   // fit function for minuit
 628   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
 629   //
 630
 631   // d
 632   static TVectorD paramsVector(fgkMaxParams);
 633   for (Int_t i=0; i<fgkMaxParams; ++i)
 634     paramsVector[i] = params[i] * params[i];
 635
 636   // calculate penalty factor
 637   Double_t penaltyVal = 0;
 638   switch (fgRegularizationType)
 639   {
 640     case kNone:       break;
 641     case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
 642     case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
 643     case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
 644     case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
 645     case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
 646   }
 647
 648   //if (penaltyVal > 1e10)
 649   //  paramsVector2.Print();
 650
 651   penaltyVal *= fgRegularizationWeight;
 652
 653   // Ad
 654   TVectorD measGuessVector(fgkMaxInput);
 655   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
 656
 657   // Ad - m
 658   measGuessVector -= (*fgCurrentESDVector);
 659
 660   TVectorD copy(measGuessVector);
 661
 662   // (Ad - m) W
 663   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
 664   // normal way is like this:
 665   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
 666   // optimized way like this:
 667   for (Int_t i=0; i<fgkMaxParams; ++i)
 668     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
 669
 670   //measGuessVector.Print();
 671
 672   // (Ad - m) W (Ad - m)
 673   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
 674   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
 675   Double_t chi2FromFit = measGuessVector * copy * 1e6;
 676
 677   chi2 = chi2FromFit + penaltyVal;
 678
 679   static Int_t callCount = 0;
 680   if ((callCount++ % 10000) == 0)
 681     printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2);
 682 }
 683
 684 //____________________________________________________________________
 685 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
 686 {
 687   //
 688   // fills fCurrentESD, fCurrentCorrelation
 689   // resets fMultiplicityESDCorrected
 690   //
 691
 692   Bool_t display = kFALSE;
 693
 694   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 695
 696   fMultiplicityESDCorrected[correlationID]->Reset();
 697   fMultiplicityESDCorrected[correlationID]->Sumw2();
 698
 699   fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD"); //"fCurrentESD", -1, -1, "e");
 700   fCurrentESD->Sumw2();
 701
 702   // empty under/overflow bins in x, otherwise Project3D takes them into account
 703   TH3* hist = fCorrelation[correlationID];
 704   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
 705   {
 706     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
 707     {
 708       hist->SetBinContent(0, y, z, 0);
 709       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
 710     }
 711   }
 712
 713   fCurrentCorrelation = hist->Project3D("zy");
 714   //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1);
 715   //fMultiplicityESDCorrected[correlationID]->Rebin(2);
 716   fCurrentCorrelation->Sumw2();
 717
 718   if (createBigBin)
 719   {
 720     Int_t maxBin = 0;
 721     for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
 722     {
 723       if (fCurrentESD->GetBinContent(i) <= 5)
 724       {
 725         maxBin = i;
 726         break;
 727       }
 728     }
 729
 730     if (maxBin > 0)
 731     {
 732       TCanvas* canvas = 0;
 733
 734       if (display)
 735       {
 736         canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
 737         canvas->Divide(2, 2);
 738
 739         canvas->cd(1);
 740         fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
 741         fCurrentESD->SetStats(kFALSE);
 742         fCurrentESD->SetTitle(";measured multiplicity");
 743         fCurrentESD->DrawCopy();
 744         gPad->SetLogy();
 745
 746         canvas->cd(2);
 747         fCurrentCorrelation->GetXaxis()->SetRangeUser(0, 250);
 748         fCurrentCorrelation->GetYaxis()->SetRangeUser(0, 250);
 749         fCurrentCorrelation->SetTitle(";true multiplicity;measured multiplicity");
 750         fCurrentCorrelation->SetStats(kFALSE);
 751         fCurrentCorrelation->DrawCopy("COLZ");
 752       }
 753
 754       printf("Bin limit in measured spectrum is %d.\n", maxBin);
 755       fCurrentESD->SetBinContent(maxBin, fCurrentESD->Integral(maxBin, fCurrentESD->GetNbinsX()));
 756       for (Int_t i=maxBin+1; i<=fCurrentESD->GetNbinsX(); ++i)
 757       {
 758         fCurrentESD->SetBinContent(i, 0);
 759         fCurrentESD->SetBinError(i, 0);
 760       }
 761       // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
 762       fCurrentESD->SetBinError(maxBin, TMath::Sqrt(fCurrentESD->GetBinContent(maxBin)));
 763
 764       printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(maxBin), fCurrentESD->GetBinError(maxBin));
 765
 766       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
 767       {
 768         fCurrentCorrelation->SetBinContent(i, maxBin, fCurrentCorrelation->Integral(i, i, maxBin, fCurrentCorrelation->GetNbinsY()));
 769         // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
 770         fCurrentCorrelation->SetBinError(i, maxBin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, maxBin)));
 771
 772         for (Int_t j=maxBin+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
 773         {
 774           fCurrentCorrelation->SetBinContent(i, j, 0);
 775           fCurrentCorrelation->SetBinError(i, j, 0);
 776         }
 777       }
 778
 779       printf("Adjusted correlation matrix!\n");
 780
 781       if (canvas)
 782       {
 783         canvas->cd(3);
 784         fCurrentESD->DrawCopy();
 785         gPad->SetLogy();
 786
 787         canvas->cd(4);
 788         fCurrentCorrelation->DrawCopy("COLZ");
 789       }
 790     }
 791   }
 792
 793 #if 0 // does not help
 794   // null bins with one entry
 795   Int_t nNulledBins = 0;
 796   for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x)
 797     for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y)
 798     {
 799       if (fCurrentCorrelation->GetBinContent(x, y) == 1)
 800       {
 801         fCurrentCorrelation->SetBinContent(x, y, 0);
 802         fCurrentCorrelation->SetBinError(x, y, 0);
 803
 804         ++nNulledBins;
 805       }
 806     }
 807   Printf("Nulled %d bins", nNulledBins);
 808 #endif
 809
 810   fCurrentEfficiency = GetEfficiency(inputRange, eventType);
 811   //fCurrentEfficiency->Rebin(2);
 812   //fCurrentEfficiency->Scale(0.5);
 813 }
 814
 815 //____________________________________________________________________
 816 TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType)
 817 {
 818   //
 819   // calculates efficiency for given event type
 820   //
 821
 822   TH1* divisor = 0;
 823   switch (eventType)
 824   {
 825     case kTrVtx : divisor = fMultiplicityVtx[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
 826     case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
 827     case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
 828   }
 829   TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", -1, -1, "e");
 830   eff->Divide(eff, divisor, 1, 1, "B");
 831   return eff;
 832 }
 833
 834 //____________________________________________________________________
 835 void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight)
 836 {
 837   //
 838   // sets the parameters for chi2 minimization
 839   //
 840
 841   fgRegularizationType = type;
 842   fgRegularizationWeight = weight;
 843
 844   printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight);
 845 }
 846
 847 //____________________________________________________________________
 848 void AliMultiplicityCorrection::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
 849 {
 850   //
 851   // sets the parameters for Bayesian unfolding
 852   //
 853
 854   fgBayesianSmoothing = smoothing;
 855   fgBayesianIterations = nIterations;
 856
 857   printf("AliMultiplicityCorrection::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
 858 }
 859
 860 //____________________________________________________________________
 861 Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
 862 {
 863   //
 864   // implemenation of unfolding (internal function)
 865   //
 866   // unfolds <measured> using response from <correlation> and effiency <aEfficiency>
 867   // output is in <result>
 868   // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
 869   // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
 870   //
 871   // returns minuit status (0 = success), (-1 when check was set)
 872   //
 873
 874   //normalize measured
 875   measured->Scale(1.0 / measured->Integral());
 876
 877   if (!fgCorrelationMatrix)
 878     fgCorrelationMatrix = new TMatrixD(fgkMaxInput, fgkMaxParams);
 879   if (!fgCorrelationCovarianceMatrix)
 880     fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxInput, fgkMaxInput);
 881   if (!fgCurrentESDVector)
 882     fgCurrentESDVector = new TVectorD(fgkMaxInput);
 883   if (!fgEntropyAPriori)
 884     fgEntropyAPriori = new TVectorD(fgkMaxParams);
 885
 886   fgCorrelationMatrix->Zero();
 887   fgCorrelationCovarianceMatrix->Zero();
 888   fgCurrentESDVector->Zero();
 889   fgEntropyAPriori->Zero();
 890
 891   // normalize correction for given nPart
 892   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
 893   {
 894     Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
 895     if (sum <= 0)
 896       continue;
 897     Float_t maxValue = 0;
 898     Int_t maxBin = -1;
 899     for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
 900     {
 901       // find most probably value
 902       if (maxValue < correlation->GetBinContent(i, j))
 903       {
 904         maxValue = correlation->GetBinContent(i, j);
 905         maxBin = j;
 906       }
 907
 908       // npart sum to 1
 909       correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
 910       correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
 911
 912       if (i <= fgkMaxParams && j <= fgkMaxInput)
 913         (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
 914     }
 915
 916     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
 917   }
 918
 919   // Initialize TMinuit via generic fitter interface
 920   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgkMaxParams);
 921   Double_t arglist[100];
 922
 923   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
 924   arglist[0] = 0;
 925   minuit->ExecuteCommand("SET PRINT", arglist, 1);
 926
 927   // however, enable warnings
 928   //minuit->ExecuteCommand("SET WAR", arglist, 0);
 929
 930   // set minimization function
 931   minuit->SetFCN(MinuitFitFunction);
 932
 933   for (Int_t i=0; i<fgkMaxParams; i++)
 934     (*fgEntropyAPriori)[i] = 1;
 935
 936   if (initialConditions)
 937   {
 938     printf("Using different starting conditions...\n");
 939     //new TCanvas; initialConditions->DrawCopy();
 940
 941     initialConditions->Scale(1.0 / initialConditions->Integral());
 942     for (Int_t i=0; i<fgkMaxParams; i++)
 943       if (initialConditions->GetBinContent(i+1) > 0)
 944         (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
 945   }
 946   else
 947     initialConditions = measured;
 948
 949   Double_t results[fgkMaxParams+1];
 950   for (Int_t i=0; i<fgkMaxParams; ++i)
 951   {
 952     results[i] = initialConditions->GetBinContent(i+1);
 953
 954     // minimum value
 955     results[i] = TMath::Max(results[i], 1e-3);
 956
 957     Float_t stepSize = 0.1;
 958
 959     // minuit sees squared values to prevent it from going negative...
 960     results[i] = TMath::Sqrt(results[i]);
 961     //stepSize /= results[i]; // keep relative step size
 962
 963     minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
 964   }
 965   //minuit->SetParameter(fgkMaxParams, "alpha", 1e4, 1, 0, 0);
 966   // bin 0 is filled with value from bin 1 (otherwise it's 0)
 967   //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0);
 968   //results[0] = 0;
 969   //minuit->SetParameter(0, "param0", 0, 0, 0, 0);
 970
 971   for (Int_t i=0; i<fgkMaxInput; ++i)
 972   {
 973     (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
 974     if (measured->GetBinError(i+1) > 0)
 975       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
 976
 977     if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
 978       (*fgCorrelationCovarianceMatrix)(i, i) = 0;
 979   }
 980
 981   Int_t dummy;
 982   Double_t chi2;
 983   MinuitFitFunction(dummy, 0, chi2, results, 0);
 984   printf("Chi2 of initial parameters is = %f\n", chi2);
 985
 986   if (check)
 987     return -1;
 988
 989   // first param is number of iterations, second is precision....
 990   arglist[0] = 1e6;
 991   //arglist[1] = 1e-5;
 992   //minuit->ExecuteCommand("SCAN", arglist, 0);
 993   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
 994   printf("MINUIT status is %d\n", status);
 995   //minuit->ExecuteCommand("MIGRAD", arglist, 1);
 996   //minuit->ExecuteCommand("MIGRAD", arglist, 1);
 997   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
 998   //minuit->ExecuteCommand("MINOS", arglist, 0);
 999
1000   for (Int_t i=0; i<fgkMaxParams; ++i)
1001   {
1002     if (aEfficiency->GetBinContent(i+1) > 0)
1003     {
1004       result->SetBinContent(i+1, minuit->GetParameter(i) * minuit->GetParameter(i) / aEfficiency->GetBinContent(i+1));
1005       // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
1006       result->SetBinError(i+1, minuit->GetParError(i) * minuit->GetParameter(i) /  aEfficiency->GetBinContent(i+1));
1007     }
1008     else
1009     {
1010       result->SetBinContent(i+1, 0);
1011       result->SetBinError(i+1, 0);
1012     }
1013   }
1014
1015   return status;
1016 }
1017
1018 //____________________________________________________________________
1019 Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
1020 {
1021   //
1022   // correct spectrum using minuit chi2 method
1023   //
1024   // for description of parameters, see UnfoldWithMinuit
1025   //
1026
1027   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1028   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kTRUE);
1029
1030   return UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
1031 }
1032
1033 //____________________________________________________________________
1034 void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
1035 {
1036   //
1037   // correct spectrum using minuit chi2 method applying a NBD fit
1038   //
1039
1040   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1041
1042   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
1043   //normalize ESD
1044   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1045
1046   fgCorrelationMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
1047   fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
1048   fgCurrentESDVector = new TVectorD(fgkMaxParams);
1049
1050   // normalize correction for given nPart
1051   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1052   {
1053     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1054     if (sum <= 0)
1055       continue;
1056     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1057     {
1058       // npart sum to 1
1059       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1060       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1061
1062       if (i <= fgkMaxParams && j <= fgkMaxParams)
1063         (*fgCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
1064     }
1065   }
1066
1067   for (Int_t i=0; i<fgkMaxParams; ++i)
1068   {
1069     (*fgCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
1070     if (fCurrentESD->GetBinError(i+1) > 0)
1071       (*fgCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
1072   }
1073
1074   // Create NBD function
1075   if (!fgNBD)
1076   {
1077     fgNBD = 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);
1078     fgNBD->SetParNames("scaling", "averagen", "k");
1079   }
1080
1081   // Initialize TMinuit via generic fitter interface
1082   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3);
1083
1084   minuit->SetFCN(MinuitNBD);
1085   minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000);
1086   minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000);
1087   minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000);
1088
1089   Double_t arglist[100];
1090   arglist[0] = 0;
1091   minuit->ExecuteCommand("SET PRINT", arglist, 1);
1092   minuit->ExecuteCommand("MIGRAD", arglist, 0);
1093   //minuit->ExecuteCommand("MINOS", arglist, 0);
1094
1095   fgNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
1096
1097   for (Int_t i=0; i<fgkMaxParams; ++i)
1098   {
1099     printf("%d : %f\n", i, fgNBD->Eval(i));
1100     if (fgNBD->Eval(i) > 0)
1101     {
1102       fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fgNBD->Eval(i));
1103       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
1104     }
1105   }
1106 }
1107
1108 //____________________________________________________________________
1109 void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
1110 {
1111   //
1112   // normalizes a 1-d histogram to its bin width
1113   //
1114
1115   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
1116   {
1117     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
1118     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
1119   }
1120 }
1121
1122 //____________________________________________________________________
1123 void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
1124 {
1125   //
1126   // normalizes a 2-d histogram to its bin width (x width * y width)
1127   //
1128
1129   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
1130     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
1131     {
1132       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
1133       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
1134       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
1135     }
1136 }
1137
1138 //____________________________________________________________________
1139 void AliMultiplicityCorrection::DrawHistograms()
1140 {
1141   //
1142   // draws the histograms of this class
1143   //
1144
1145   printf("ESD:\n");
1146
1147   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
1148   canvas1->Divide(3, 2);
1149   for (Int_t i = 0; i < kESDHists; ++i)
1150   {
1151     canvas1->cd(i+1);
1152     fMultiplicityESD[i]->DrawCopy("COLZ");
1153     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
1154   }
1155
1156   printf("Vtx:\n");
1157
1158   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
1159   canvas2->Divide(3, 2);
1160   for (Int_t i = 0; i < kMCHists; ++i)
1161   {
1162     canvas2->cd(i+1);
1163     fMultiplicityVtx[i]->DrawCopy("COLZ");
1164     printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
1165   }
1166
1167   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
1168   canvas3->Divide(3, 3);
1169   for (Int_t i = 0; i < kCorrHists; ++i)
1170   {
1171     canvas3->cd(i+1);
1172     TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
1173     for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1174     {
1175       for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1176       {
1177         hist->SetBinContent(0, y, z, 0);
1178         hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1179       }
1180     }
1181     TH1* proj = hist->Project3D("zy");
1182     proj->DrawCopy("COLZ");
1183   }
1184 }
1185
1186 //____________________________________________________________________
1187 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple)
1188 {
1189   //mcHist->Rebin(2);
1190
1191   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1192
1193   TString tmpStr;
1194   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
1195
1196   if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
1197   {
1198     printf("ERROR. Unfolded histogram is empty\n");
1199     return;
1200   }
1201
1202   //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
1203   fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY();
1204   fCurrentESD->Sumw2();
1205   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1206
1207   // normalize unfolded result to 1
1208   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
1209
1210   //fCurrentESD->Scale(mcHist->Integral(2, 200));
1211
1212   //new TCanvas;
1213   /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1214   ratio->Divide(mcHist);
1215   ratio->Draw("HIST");
1216   ratio->Fit("pol0", "W0", "", 20, 120);
1217   Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
1218   delete ratio;
1219   mcHist->Scale(scalingFactor);*/
1220
1221   // find bin with <= 100 entries. this is used as limit for the chi2 calculation
1222   Int_t mcBinLimit = 0;
1223   for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
1224   {
1225     if (mcHist->GetBinContent(i) > 50)
1226     {
1227       mcBinLimit = i;
1228     }
1229     else
1230       break;
1231   }
1232   Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
1233
1234   // scale to 1
1235   mcHist->Sumw2();
1236   mcHist->Scale(1.0 / mcHist->Integral());
1237
1238   // calculate residual
1239
1240   // for that we convolute the response matrix with the gathered result
1241   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
1242   TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
1243   tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
1244   TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
1245   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
1246   if (convolutedProj->Integral() > 0)
1247   {
1248     convolutedProj->Scale(1.0 / convolutedProj->Integral());
1249   }
1250   else
1251     printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
1252
1253   TH1* residual = (TH1*) convolutedProj->Clone("residual");
1254   residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
1255
1256   residual->Add(fCurrentESD, -1);
1257   //residual->Divide(residual, fCurrentESD, 1, 1, "B");
1258
1259   TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
1260
1261   // TODO fix errors
1262   Float_t chi2 = 0;
1263   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
1264   {
1265     if (fCurrentESD->GetBinError(i) > 0)
1266     {
1267       Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
1268       if (i > 1)
1269         chi2 += value * value;
1270       residual->SetBinContent(i, value);
1271       residualHist->Fill(value);
1272     }
1273     else
1274     {
1275       //printf("Residual bin %d set to 0\n", i);
1276       residual->SetBinContent(i, 0);
1277     }
1278     convolutedProj->SetBinError(i, 0);
1279     residual->SetBinError(i, 0);
1280
1281     if (i == 200)
1282       fLastChi2Residuals = chi2;
1283   }
1284
1285   new TCanvas; residualHist->DrawCopy();
1286
1287   //residualHist->Fit("gaus", "N");
1288   //delete residualHist;
1289
1290   printf("Difference (Residuals) is %f for bin 2-200\n", chi2);
1291
1292   TCanvas* canvas1 = 0;
1293   if (simple)
1294   {
1295     canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
1296     canvas1->Divide(2, 1);
1297   }
1298   else
1299   {
1300     canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
1301     canvas1->Divide(2, 3);
1302   }
1303
1304   canvas1->cd(1);
1305   TH1* proj = (TH1*) mcHist->Clone("proj");
1306
1307   proj->GetXaxis()->SetRangeUser(0, 200);
1308   proj->SetTitle(";true multiplicity;Entries");
1309   proj->SetStats(kFALSE);
1310   proj->DrawCopy("HIST");
1311   gPad->SetLogy();
1312
1313   //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
1314   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
1315   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
1316
1317   TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1318   legend->AddEntry(proj, "true distribution");
1319   legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
1320   legend->SetFillColor(0);
1321   legend->Draw();
1322   // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
1323
1324   canvas1->cd(2);
1325
1326   gPad->SetLogy();
1327   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1328   //fCurrentESD->SetLineColor(2);
1329   fCurrentESD->SetTitle(";measured multiplicity;Entries");
1330   fCurrentESD->SetStats(kFALSE);
1331   fCurrentESD->DrawCopy("HIST E");
1332
1333   convolutedProj->SetLineColor(2);
1334   //proj2->SetMarkerColor(2);
1335   //proj2->SetMarkerStyle(5);
1336   convolutedProj->DrawCopy("HIST SAME");
1337
1338   legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1339   legend->AddEntry(fCurrentESD, "measured distribution");
1340   legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
1341   legend->SetFillColor(0);
1342   legend->Draw();
1343
1344   //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
1345   //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1346
1347   /*Float_t chi2 = 0;
1348   Float_t chi = 0;
1349   fLastChi2MCLimit = 0;
1350   Int_t limit = 0;
1351   for (Int_t i=2; i<=200; ++i)
1352   {
1353     if (proj->GetBinContent(i) != 0)
1354     {
1355       Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
1356       chi2 += value * value;
1357       chi += TMath::Abs(value);
1358
1359       //printf("%d %f\n", i, chi);
1360
1361       if (chi2 < 0.2)
1362         fLastChi2MCLimit = i;
1363
1364       if (chi < 3)
1365         limit = i;
1366
1367     }
1368   }*/
1369
1370   /*chi2 = 0;
1371   Float_t chi = 0;
1372   Int_t limit = 0;
1373   for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
1374   {
1375     if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
1376     {
1377       Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
1378       if (value > 1e8)
1379         value = 1e8; //prevent arithmetic exception
1380       else if (value < -1e8)
1381         value = -1e8;
1382       if (i > 1)
1383       {
1384         chi2 += value * value;
1385         chi += TMath::Abs(value);
1386       }
1387       diffMCUnfolded->SetBinContent(i, value);
1388     }
1389     else
1390     {
1391       //printf("diffMCUnfolded bin %d set to 0\n", i);
1392       diffMCUnfolded->SetBinContent(i, 0);
1393     }
1394     if (chi2 < 1000)
1395       fLastChi2MCLimit = i;
1396     if (chi < 1000)
1397       limit = i;
1398     if (i == 150)
1399       fLastChi2MC = chi2;
1400   }
1401
1402   printf("limits %d %d\n", fLastChi2MCLimit, limit);
1403   fLastChi2MCLimit = limit;
1404
1405   printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
1406
1407   if (!simple)
1408   {
1409     /*canvas1->cd(3);
1410
1411     diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
1412     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1413     diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
1414     diffMCUnfolded->DrawCopy("HIST");
1415
1416     TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
1417     for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
1418       fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
1419
1420     //new TCanvas; fluctuation->DrawCopy();
1421     delete fluctuation;*/
1422
1423     /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1424     legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1425     legend->AddEntry(fMultiplicityMC, "MC");
1426     legend->AddEntry(fMultiplicityESD, "ESD");
1427     legend->Draw();*/
1428
1429     canvas1->cd(4);
1430     //residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
1431     residual->GetXaxis()->SetRangeUser(0, 200);
1432     residual->DrawCopy();
1433
1434     canvas1->cd(5);
1435
1436     TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1437     ratio->Divide(mcHist);
1438     ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1439     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1440     ratio->GetXaxis()->SetRangeUser(0, 200);
1441     ratio->SetStats(kFALSE);
1442     ratio->Draw("HIST");
1443
1444     Double_t ratioChi2 = 0;
1445     fRatioAverage = 0;
1446     fLastChi2MCLimit = 0;
1447     for (Int_t i=2; i<=150; ++i)
1448     {
1449       Float_t value = ratio->GetBinContent(i) - 1;
1450       if (value > 1e8)
1451         value = 1e8; //prevent arithmetic exception
1452       else if (value < -1e8)
1453         value = -1e8;
1454
1455       ratioChi2 += value * value;
1456       fRatioAverage += TMath::Abs(value);
1457
1458       if (ratioChi2 < 0.1)
1459         fLastChi2MCLimit = i;
1460     }
1461     fRatioAverage /= 149;
1462
1463     printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
1464     // TODO FAKE
1465     fLastChi2MC = ratioChi2;
1466
1467     // FFT of ratio
1468     canvas1->cd(6);
1469     const Int_t kFFT = 128;
1470     Double_t fftReal[kFFT];
1471     Double_t fftImag[kFFT];
1472
1473     for (Int_t i=0; i<kFFT; ++i)
1474     {
1475       fftReal[i] = ratio->GetBinContent(i+1+10);
1476       // test: ;-)
1477       //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
1478       fftImag[i] = 0;
1479     }
1480
1481     FFT(-1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1482
1483     TH1* fftResult = (TH1*) ratio->Clone("fftResult");
1484     fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
1485     TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
1486     TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
1487     fftResult->Reset();
1488     fftResult2->Reset();
1489     fftResult3->Reset();
1490     fftResult->GetYaxis()->UnZoom();
1491     fftResult2->GetYaxis()->UnZoom();
1492     fftResult3->GetYaxis()->UnZoom();
1493
1494     Printf("AFTER FFT");
1495     for (Int_t i=0; i<kFFT; ++i)
1496     {
1497       Printf("%d: %f", i, fftReal[i]);
1498       fftResult->SetBinContent(i+1, fftReal[i]);
1499       /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
1500       {
1501         Printf("Nulled %d", i);
1502         fftReal[i] = 0;
1503       }*/
1504     }
1505
1506     fftResult->SetLineColor(1);
1507     fftResult->DrawCopy();
1508
1509
1510     Printf("IMAG");
1511     for (Int_t i=0; i<kFFT; ++i)
1512     {
1513       Printf("%d: %f", i, fftImag[i]);
1514       fftResult2->SetBinContent(i+1, fftImag[i]);
1515
1516       fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
1517     }
1518
1519     fftResult2->SetLineColor(2);
1520     fftResult2->DrawCopy("SAME");
1521
1522     fftResult3->SetLineColor(4);
1523     fftResult3->DrawCopy("SAME");
1524
1525     for (Int_t i=1; i<kFFT - 1; ++i)
1526     {
1527       if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
1528       {
1529         fftReal[i-1] = 0;
1530         fftReal[i] = 0;
1531         fftReal[i+1] = 0;
1532         fftImag[i-1] = 0;
1533         fftImag[i] = 0;
1534         fftImag[i+1] = 0;
1535         //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
1536         //fftImag[i]  = (fftImag[i-1] + fftImag[i+1]) / 2;
1537         Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
1538       }
1539     }
1540
1541     FFT(1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1542
1543     TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
1544     fftResult4->Reset();
1545
1546     Printf("AFTER BACK-TRAFO");
1547     for (Int_t i=0; i<kFFT; ++i)
1548     {
1549       Printf("%d: %f", i, fftReal[i]);
1550       fftResult4->SetBinContent(i+1+10, fftReal[i]);
1551     }
1552
1553     canvas1->cd(5);
1554     fftResult4->SetLineColor(4);
1555     fftResult4->DrawCopy("SAME");
1556
1557     // plot (MC - Unfolded) / error (MC)
1558     canvas1->cd(3);
1559
1560     TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
1561     diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1562
1563     Int_t ndfQual[kQualityRegions];
1564     for (Int_t region=0; region<kQualityRegions; ++region)
1565     {
1566       fQuality[region] = 0;
1567       ndfQual[region] = 0;
1568     }
1569
1570
1571     Double_t newChi2 = 0;
1572     Double_t newChi2Limit150 = 0;
1573     Int_t ndf = 0;
1574     for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i)
1575     {
1576       Double_t value = 0;
1577       if (proj->GetBinError(i) > 0)
1578       {
1579         value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
1580         newChi2 += value * value;
1581         if (i > 1 && i <= mcBinLimit)
1582           newChi2Limit150 += value * value;
1583         ++ndf;
1584
1585         for (Int_t region=0; region<kQualityRegions; ++region)
1586           if (diffMCUnfolded2->GetXaxis()->GetBinCenter(i) >= fgQualityRegionsB[region] - 0.1 && diffMCUnfolded2->GetXaxis()->GetBinCenter(i) <= fgQualityRegionsE[region] + 0.1) // 0.1 to avoid e.g. 3.9999 < 4 problem
1587           {
1588             fQuality[region] += TMath::Abs(value);
1589             ++ndfQual[region];
1590           }
1591       }
1592
1593       diffMCUnfolded2->SetBinContent(i, value);
1594     }
1595
1596     // normalize region to the number of entries
1597     for (Int_t region=0; region<kQualityRegions; ++region)
1598     {
1599       if (ndfQual[region] > 0)
1600         fQuality[region] /= ndfQual[region];
1601       Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
1602     }
1603
1604     if (mcBinLimit > 1)
1605     {
1606       // TODO another FAKE
1607       fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
1608       Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
1609     }
1610     else
1611       fLastChi2MC = -1;
1612
1613     Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, newChi2 / ndf);
1614
1615     diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
1616     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1617     diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
1618     diffMCUnfolded2->DrawCopy("HIST");
1619   }
1620
1621   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
1622 }
1623
1624 //____________________________________________________________________
1625 void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y)
1626 {
1627 /*-------------------------------------------------------------------------
1628    This computes an in-place complex-to-complex FFT
1629    x and y are the real and imaginary arrays of 2^m points.
1630    dir =  1 gives forward transform
1631    dir = -1 gives reverse transform
1632
1633      Formula: forward
1634                   N-1
1635                   ---
1636               1   \          - j k 2 pi n / N
1637       X(n) = ---   >   x(k) e                    = forward transform
1638               N   /                                n=0..N-1
1639                   ---
1640                   k=0
1641
1642       Formula: reverse
1643                   N-1
1644                   ---
1645                   \          j k 2 pi n / N
1646       X(n) =       >   x(k) e                    = forward transform
1647                   /                                n=0..N-1
1648                   ---
1649                   k=0
1650 */
1651
1652    Long_t   nn, i, i1, j, k, i2, l, l1, l2;
1653    Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
1654
1655    /* Calculate the number of points */
1656    nn = 1;
1657    for (i = 0; i < m; i++)
1658        nn *= 2;
1659
1660    /* Do the bit reversal */
1661    i2 = nn >> 1;
1662    j = 0;
1663    for (i= 0; i < nn - 1; i++) {
1664        if (i < j) {
1665            tx = x[i];
1666            ty = y[i];
1667            x[i] = x[j];
1668            y[i] = y[j];
1669            x[j] = tx;
1670            y[j] = ty;
1671        }
1672        k = i2;
1673        while (k <= j) {
1674            j -= k;
1675            k >>= 1;
1676        }
1677        j += k;
1678    }
1679
1680    /* Compute the FFT */
1681    c1 = -1.0;
1682    c2 = 0.0;
1683    l2 = 1;
1684    for (l = 0; l < m; l++) {
1685        l1 = l2;
1686        l2 <<= 1;
1687        u1 = 1.0;
1688        u2 = 0.0;
1689        for (j = 0;j < l1; j++) {
1690            for (i = j; i < nn; i += l2) {
1691                i1 = i + l1;
1692                t1 = u1 * x[i1] - u2 * y[i1];
1693                t2 = u1 * y[i1] + u2 * x[i1];
1694                x[i1] = x[i] - t1;
1695                y[i1] = y[i] - t2;
1696                x[i] += t1;
1697                y[i] += t2;
1698            }
1699            z =  u1 * c1 - u2 * c2;
1700            u2 = u1 * c2 + u2 * c1;
1701            u1 = z;
1702        }
1703        c2 = TMath::Sqrt((1.0 - c1) / 2.0);
1704        if (dir == 1)
1705            c2 = -c2;
1706        c1 = TMath::Sqrt((1.0 + c1) / 2.0);
1707    }
1708
1709    /* Scaling for forward transform */
1710    if (dir == 1) {
1711        for (i=0;i<nn;i++) {
1712            x[i] /= (Double_t)nn;
1713            y[i] /= (Double_t)nn;
1714        }
1715    }
1716 }
1717
1718 //____________________________________________________________________
1719 void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage) const
1720 {
1721   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1722   // These values are computed during DrawComparison, thus this function picks up the
1723   // last calculation
1724
1725   if (mc)
1726     *mc = fLastChi2MC;
1727   if (mcLimit)
1728     *mcLimit = fLastChi2MCLimit;
1729   if (residuals)
1730     *residuals = fLastChi2Residuals;
1731   if (ratioAverage)
1732     *ratioAverage = fRatioAverage;
1733 }
1734
1735 //____________________________________________________________________
1736 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1737 {
1738   //
1739   // returns the corresponding MC spectrum
1740   //
1741
1742   switch (eventType)
1743   {
1744     case kTrVtx : return fMultiplicityVtx[i]; break;
1745     case kMB : return fMultiplicityMB[i]; break;
1746     case kINEL : return fMultiplicityINEL[i]; break;
1747   }
1748
1749   return 0;
1750 }
1751
1752 //____________________________________________________________________
1753 TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
1754 {
1755   //
1756   // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
1757   // the function unfolds the spectrum using the default response matrix and several modified ones
1758   // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
1759   // these unfolded results are compared to the first result gained with the default response OR to the histogram given
1760   // in <compareTo> (optional)
1761   //
1762   // returns the error assigned to the measurement
1763   //
1764
1765   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1766
1767   // initialize seed with current time
1768   gRandom->SetSeed(0);
1769
1770   const Int_t kErrorIterations = 150;
1771
1772   TH1* maxError = 0;
1773   TH1* firstResult = 0;
1774
1775   TH1** results = new TH1*[kErrorIterations];
1776
1777   for (Int_t n=0; n<kErrorIterations; ++n)
1778   {
1779     Printf("Iteration %d of %d...", n, kErrorIterations);
1780
1781     // only done for chi2 minimization
1782     Bool_t createBigBin = kFALSE;
1783     if (methodType == kChi2Minimization)
1784       createBigBin = kTRUE;
1785
1786     SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin);
1787
1788     TH1* measured = (TH1*) fCurrentESD->Clone("measured");
1789
1790     if (n > 0)
1791     {
1792       if (randomizeResponse)
1793       {
1794         // randomize response matrix
1795         for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1796           for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1797             fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
1798       }
1799
1800       if (randomizeMeasured)
1801       {
1802         // randomize measured spectrum
1803         for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
1804         {
1805           Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1806           measured->SetBinContent(x, randomValue);
1807           measured->SetBinError(x, TMath::Sqrt(randomValue));
1808         }
1809       }
1810     }
1811
1812     // only for bayesian method we have to do it before the call to Unfold...
1813     if (methodType == kBayesian)
1814     {
1815       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1816       {
1817         // with this it is normalized to 1
1818         Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1819
1820         // with this normalized to the given efficiency
1821         if (fCurrentEfficiency->GetBinContent(i) > 0)
1822           sum /= fCurrentEfficiency->GetBinContent(i);
1823         else
1824           sum = 0;
1825
1826         for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1827         {
1828           if (sum > 0)
1829           {
1830             fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1831             fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1832           }
1833           else
1834           {
1835             fCurrentCorrelation->SetBinContent(i, j, 0);
1836             fCurrentCorrelation->SetBinError(i, j, 0);
1837           }
1838         }
1839       }
1840     }
1841
1842     TH1* result = 0;
1843     if (n == 0 && compareTo)
1844     {
1845       // in this case we just store the histogram we want to compare to
1846       result = (TH1*) compareTo->Clone("compareTo");
1847       result->Sumw2();
1848     }
1849     else
1850     {
1851       result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
1852
1853       Int_t returnCode = -1;
1854       if (methodType == kChi2Minimization)
1855       {
1856         returnCode = UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, kFALSE);
1857       }
1858       else if (methodType == kBayesian)
1859         returnCode = UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, fgBayesianSmoothing, fgBayesianIterations);
1860
1861       if (returnCode != 0)
1862         return 0;
1863     }
1864
1865     // normalize
1866     result->Scale(1.0 / result->Integral());
1867
1868     if (n == 0)
1869     {
1870       firstResult = (TH1*) result->Clone("firstResult");
1871
1872       maxError = (TH1*) result->Clone("maxError");
1873       maxError->Reset();
1874     }
1875     else
1876     {
1877       // calculate ratio
1878       TH1* ratio = (TH1*) firstResult->Clone("ratio");
1879       ratio->Divide(result);
1880
1881       // find max. deviation
1882       for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
1883         maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
1884
1885       delete ratio;
1886     }
1887
1888     results[n] = result;
1889   }
1890
1891   // find covariance matrix
1892   // results[n] is X_x
1893   // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
1894
1895   Int_t nBins = results[0]->GetNbinsX();
1896   Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
1897   Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
1898
1899   // find average, E(X)
1900   TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
1901   for (Int_t n=1; n<kErrorIterations; ++n)
1902     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1903       average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
1904   //new TCanvas; average->DrawClone();
1905
1906   // find cov. matrix
1907   TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
1908
1909   for (Int_t n=1; n<kErrorIterations; ++n)
1910     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1911       for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
1912       {
1913         // (X_x - E(X_x)) * (X_y - E(X_y)
1914         Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
1915         covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
1916       }
1917   new TCanvas; covMatrix->DrawClone("COLZ");
1918
1919 //   // fill 2D histogram that contains deviation from first
1920 //   TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
1921 //   for (Int_t n=1; n<kErrorIterations; ++n)
1922 //     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1923 //       deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
1924 //   //new TCanvas; deviations->DrawCopy("COLZ");
1925 //
1926 //   // get standard deviation "by hand"
1927 //   for (Int_t x=1; x<=nBins; x++)
1928 //   {
1929 //     if (results[0]->GetBinContent(x) > 0)
1930 //     {
1931 //       TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
1932 //       Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
1933 //       //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1934 //       Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1935 //     }
1936 //   }
1937
1938   // calculate standard deviation of (results[0] - results[n])
1939   TH1* standardDeviation = (TH1*) maxError->Clone("standardDeviation");
1940   standardDeviation->Reset();
1941
1942   for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1943   {
1944     if (results[0]->GetBinContent(x) > 0)
1945     {
1946       Double_t average = 0;
1947       for (Int_t n=1; n<kErrorIterations; ++n)
1948         average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
1949       average /= kErrorIterations-1;
1950
1951       Double_t variance = 0;
1952       for (Int_t n=1; n<kErrorIterations; ++n)
1953       {
1954         Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
1955         variance += value * value;
1956       }
1957       variance /= kErrorIterations-1;
1958
1959       Double_t standardDev = TMath::Sqrt(variance);
1960       standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1961       Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1962     }
1963   }
1964
1965   // compare maxError to RMS of profile (variable name: average)
1966   // first: calculate rms in percent of value
1967   TH1* rmsError = (TH1*) maxError->Clone("rmsError");
1968   rmsError->Reset();
1969
1970   // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
1971   average->SetErrorOption("s");
1972   for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
1973     if (average->GetBinContent(x) > 0)
1974       rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
1975
1976   // find maxError deviation from average (not from "first result"/mc as above)
1977   TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
1978   maxError2->Reset();
1979   for (Int_t n=1; n<kErrorIterations; ++n)
1980     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1981       if (average->GetBinContent(x) > 0)
1982         maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
1983
1984   //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
1985
1986   // plot difference between average and MC/first result
1987   TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
1988   averageFirstRatio->Reset();
1989   averageFirstRatio->Divide(results[0], average);
1990
1991   //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
1992   //new TCanvas; averageFirstRatio->DrawCopy();
1993
1994   // clean up
1995   for (Int_t n=0; n<kErrorIterations; ++n)
1996     delete results[n];
1997   delete[] results;
1998
1999   // fill into result histogram
2000   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2001     fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
2002
2003   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2004     fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2005
2006   return standardDeviation;
2007 }
2008
2009 //____________________________________________________________________
2010 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
2011 {
2012   //
2013   // correct spectrum using bayesian method
2014   //
2015
2016   // initialize seed with current time
2017   gRandom->SetSeed(0);
2018
2019   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2020
2021   // normalize correction for given nPart
2022   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2023   {
2024     // with this it is normalized to 1
2025     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2026
2027     // with this normalized to the given efficiency
2028     if (fCurrentEfficiency->GetBinContent(i) > 0)
2029       sum /= fCurrentEfficiency->GetBinContent(i);
2030     else
2031       sum = 0;
2032
2033     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2034     {
2035       if (sum > 0)
2036       {
2037         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2038         fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2039       }
2040       else
2041       {
2042         fCurrentCorrelation->SetBinContent(i, j, 0);
2043         fCurrentCorrelation->SetBinError(i, j, 0);
2044       }
2045     }
2046   }
2047
2048   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2049
2050   if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0)
2051     return;
2052
2053   if (!determineError)
2054   {
2055     Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
2056     return;
2057   }
2058
2059   // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
2060   // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
2061   // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
2062   // error of the unfolding method itself.
2063
2064   TH1* maxError = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("maxError");
2065   maxError->Reset();
2066
2067   TH1* normalizedResult = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("normalizedResult");
2068   normalizedResult->Scale(1.0 / normalizedResult->Integral());
2069
2070   const Int_t kErrorIterations = 20;
2071
2072   printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
2073
2074   TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
2075   TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
2076   for (Int_t n=0; n<kErrorIterations; ++n)
2077   {
2078     // randomize the content of clone following a poisson with the mean = the value of that bin
2079     for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
2080     {
2081       Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
2082       //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
2083       randomized->SetBinContent(x, randomValue);
2084       randomized->SetBinError(x, TMath::Sqrt(randomValue));
2085     }
2086
2087     result2->Reset();
2088     if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0)
2089       return;
2090
2091     result2->Scale(1.0 / result2->Integral());
2092
2093     // calculate ratio
2094     result2->Divide(normalizedResult);
2095
2096     //new TCanvas; ratio->DrawCopy("HIST");
2097
2098     // find max. deviation
2099     for (Int_t x=1; x<=result2->GetNbinsX(); x++)
2100       maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - result2->GetBinContent(x))));
2101   }
2102   delete randomized;
2103   delete result2;
2104
2105   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2106     fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2107
2108   delete maxError;
2109   delete normalizedResult;
2110 }
2111
2112 //____________________________________________________________________
2113 Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations)
2114 {
2115   //
2116   // unfolds a spectrum
2117   //
2118
2119   if (measured->Integral() <= 0)
2120   {
2121     Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
2122     return -1;
2123   }
2124
2125   const Int_t kStartBin = 0;
2126
2127   const Int_t kMaxM = fgkMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
2128   const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
2129
2130   // store information in arrays, to increase processing speed (~ factor 5)
2131   Double_t measuredCopy[kMaxM];
2132   Double_t prior[kMaxT];
2133   Double_t result[kMaxT];
2134   Double_t efficiency[kMaxT];
2135
2136   Double_t response[kMaxT][kMaxM];
2137   Double_t inverseResponse[kMaxT][kMaxM];
2138
2139   // for normalization
2140   Float_t measuredIntegral = measured->Integral();
2141   for (Int_t m=0; m<kMaxM; m++)
2142   {
2143     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
2144
2145     for (Int_t t=0; t<kMaxT; t++)
2146     {
2147       response[t][m] = correlation->GetBinContent(t+1, m+1);
2148       inverseResponse[t][m] = 0;
2149     }
2150   }
2151
2152   for (Int_t t=0; t<kMaxT; t++)
2153   {
2154     efficiency[t] = aEfficiency->GetBinContent(t+1);
2155     prior[t] = measuredCopy[t];
2156     result[t] = 0;
2157   }
2158
2159   // pick prior distribution
2160   if (initialConditions)
2161   {
2162     printf("Using different starting conditions...\n");
2163     // for normalization
2164     Float_t inputDistIntegral = initialConditions->Integral();
2165     for (Int_t i=0; i<kMaxT; i++)
2166       prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
2167   }
2168
2169   //TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
2170
2171   // unfold...
2172   for (Int_t i=0; i<nIterations; i++)
2173   {
2174     //printf(" iteration %i \n", i);
2175
2176     // calculate IR from Beyes theorem
2177     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
2178
2179     for (Int_t m=0; m<kMaxM; m++)
2180     {
2181       Float_t norm = 0;
2182       for (Int_t t = kStartBin; t<kMaxT; t++)
2183         norm += response[t][m] * prior[t];
2184
2185       if (norm > 0)
2186       {
2187         for (Int_t t = kStartBin; t<kMaxT; t++)
2188           inverseResponse[t][m] = response[t][m] * prior[t] / norm;
2189       }
2190       else
2191       {
2192         for (Int_t t = kStartBin; t<kMaxT; t++)
2193           inverseResponse[t][m] = 0;
2194       }
2195     }
2196
2197     for (Int_t t = kStartBin; t<kMaxT; t++)
2198     {
2199       Float_t value = 0;
2200       for (Int_t m=0; m<kMaxM; m++)
2201         value += inverseResponse[t][m] * measuredCopy[m];
2202
2203       if (efficiency[t] > 0)
2204         result[t] = value / efficiency[t];
2205       else
2206         result[t] = 0;
2207     }
2208
2209     // regularization (simple smoothing)
2210     for (Int_t t=kStartBin; t<kMaxT; t++)
2211     {
2212       Float_t newValue = 0;
2213       // 0 bin excluded from smoothing
2214       if (t > kStartBin+1 && t<kMaxT-1)
2215       {
2216         Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
2217
2218         // weight the average with the regularization parameter
2219         newValue = (1 - regPar) * result[t] + regPar * average;
2220       }
2221       else
2222         newValue = result[t];
2223
2224       prior[t] = newValue;
2225     }
2226
2227     // calculate chi2 (change from last iteration)
2228     //Double_t chi2 = 0;
2229     /*for (Int_t t=kStartBin; t<=hTrueSmoothed->GetNbinsX(); t++)
2230     {
2231       Float_t newValue = hTrueSmoothed->GetBinContent(t);
2232       //Float_t diff = hPrior->GetBinContent(t) - newValue;
2233       //chi2 += (Double_t) diff * diff;
2234
2235       hPrior->SetBinContent(t, newValue);
2236     }
2237     //printf("Chi2 of %d iteration = %.10f\n", i, chi2);
2238     //convergence->Fill(i+1, chi2);
2239
2240     //if (i % 5 == 0)
2241     //  DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
2242
2243     delete hTrueSmoothed;*/
2244   } // end of iterations
2245
2246   //new TCanvas;
2247   //convergence->DrawCopy();
2248   //gPad->SetLogy();
2249
2250   //delete convergence;
2251
2252   for (Int_t t=0; t<kMaxT; t++)
2253     aResult->SetBinContent(t+1, result[t]);
2254
2255   return 0;
2256
2257   // ********
2258   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
2259
2260   /*printf("Calculating covariance matrix. This may take some time...\n");
2261
2262   // TODO check if this is the right one...
2263   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2264
2265   Int_t xBins = hInverseResponseBayes->GetNbinsX();
2266   Int_t yBins = hInverseResponseBayes->GetNbinsY();
2267
2268   // calculate "unfolding matrix" Mij
2269   Float_t matrixM[251][251];
2270   for (Int_t i=1; i<=xBins; i++)
2271   {
2272     for (Int_t j=1; j<=yBins; j++)
2273     {
2274       if (fCurrentEfficiency->GetBinContent(i) > 0)
2275         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
2276       else
2277         matrixM[i-1][j-1] = 0;
2278     }
2279   }
2280
2281   Float_t* vectorn = new Float_t[yBins];
2282   for (Int_t j=1; j<=yBins; j++)
2283     vectorn[j-1] = fCurrentESD->GetBinContent(j);
2284
2285   // first part of covariance matrix, depends on input distribution n(E)
2286   Float_t cov1[251][251];
2287
2288   Float_t nEvents = fCurrentESD->Integral(); // N
2289
2290   xBins = 20;
2291   yBins = 20;
2292
2293   for (Int_t k=0; k<xBins; k++)
2294   {
2295     printf("In Cov1: %d\n", k);
2296     for (Int_t l=0; l<yBins; l++)
2297     {
2298       cov1[k][l] = 0;
2299
2300       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
2301       for (Int_t j=0; j<yBins; j++)
2302         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
2303           * (1.0 - vectorn[j] / nEvents);
2304
2305       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
2306       for (Int_t i=0; i<yBins; i++)
2307         for (Int_t j=0; j<yBins; j++)
2308         {
2309           if (i == j)
2310             continue;
2311           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
2312             * vectorn[j] / nEvents;
2313          }
2314     }
2315   }
2316
2317   printf("Cov1 finished\n");
2318
2319   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
2320   cov->Reset();
2321
2322   for (Int_t i=1; i<=xBins; i++)
2323     for (Int_t j=1; j<=yBins; j++)
2324       cov->SetBinContent(i, j, cov1[i-1][j-1]);
2325
2326   new TCanvas;
2327   cov->Draw("COLZ");
2328
2329   // second part of covariance matrix, depends on response matrix
2330   Float_t cov2[251][251];
2331
2332   // Cov[P(Er|Cu), P(Es|Cu)] term
2333   Float_t covTerm[100][100][100];
2334   for (Int_t r=0; r<yBins; r++)
2335     for (Int_t u=0; u<xBins; u++)
2336       for (Int_t s=0; s<yBins; s++)
2337       {
2338         if (r == s)
2339           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2340             * (1.0 - hResponse->GetBinContent(u+1, r+1));
2341         else
2342           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2343             * hResponse->GetBinContent(u+1, s+1);
2344       }
2345
2346   for (Int_t k=0; k<xBins; k++)
2347     for (Int_t l=0; l<yBins; l++)
2348     {
2349       cov2[k][l] = 0;
2350       printf("In Cov2: %d %d\n", k, l);
2351       for (Int_t i=0; i<yBins; i++)
2352         for (Int_t j=0; j<yBins; j++)
2353         {
2354           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
2355           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
2356           Float_t tmpCov = 0;
2357           for (Int_t r=0; r<yBins; r++)
2358             for (Int_t u=0; u<xBins; u++)
2359               for (Int_t s=0; s<yBins; s++)
2360               {
2361                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
2362                   || hResponse->GetBinContent(u+1, i+1) == 0)
2363                   continue;
2364
2365                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
2366                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
2367                   * covTerm[r][u][s];
2368               }
2369
2370           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
2371         }
2372     }
2373
2374   printf("Cov2 finished\n");
2375
2376   for (Int_t i=1; i<=xBins; i++)
2377     for (Int_t j=1; j<=yBins; j++)
2378       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
2379
2380   new TCanvas;
2381   cov->Draw("COLZ");*/
2382 }
2383
2384 //____________________________________________________________________
2385 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
2386 {
2387   //
2388   // helper function for the covariance matrix of the bayesian method
2389   //
2390
2391   Float_t result = 0;
2392
2393   if (k == u && r == i)
2394     result += 1.0 / hResponse->GetBinContent(u+1, r+1);
2395
2396   if (k == u)
2397     result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
2398
2399   if (r == i)
2400     result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
2401
2402   result *= matrixM[k][i];
2403
2404   return result;
2405 }
2406
2407 //____________________________________________________________________
2408 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
2409 {
2410   //
2411   // correct spectrum using bayesian method
2412   //
2413
2414   Float_t regPar = 0;
2415
2416   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2417   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2418
2419   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2420   //normalize ESD
2421   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2422
2423   // TODO should be taken from correlation map
2424   //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2425
2426   // normalize correction for given nPart
2427   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2428   {
2429     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2430     //Double_t sum = sumHist->GetBinContent(i);
2431     if (sum <= 0)
2432       continue;
2433     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2434     {
2435       // npart sum to 1
2436       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2437       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2438     }
2439   }
2440
2441   new TCanvas;
2442   fCurrentCorrelation->Draw("COLZ");
2443
2444   // FAKE
2445   fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
2446
2447   // pick prior distribution
2448   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
2449   Float_t norm = 1; //hPrior->Integral();
2450   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
2451     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
2452
2453   // zero distribution
2454   TH1F* zero =  (TH1F*)hPrior->Clone("zero");
2455
2456   // define temp hist
2457   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
2458   hTemp->Reset();
2459
2460   // just a shortcut
2461   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
2462
2463   // unfold...
2464   Int_t iterations = 25;
2465   for (Int_t i=0; i<iterations; i++)
2466   {
2467     //printf(" iteration %i \n", i);
2468
2469     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2470     {
2471       Float_t value = 0;
2472       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
2473         value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
2474       hTemp->SetBinContent(m, value);
2475       //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
2476     }
2477
2478     // regularization (simple smoothing)
2479     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
2480
2481     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
2482     {
2483       Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
2484                          + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
2485                          + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
2486       average *= hTrueSmoothed->GetBinWidth(t);
2487
2488       // weight the average with the regularization parameter
2489       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
2490     }
2491
2492     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2493       hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
2494
2495     // fill guess
2496     for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
2497     {
2498       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
2499       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
2500
2501       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
2502     }
2503
2504
2505     // calculate chi2 (change from last iteration)
2506     Double_t chi2 = 0;
2507
2508     // use smoothed true (normalized) as new prior
2509     Float_t norm = 1; //hTrueSmoothed->Integral();
2510
2511     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
2512     {
2513       Float_t newValue = hTemp->GetBinContent(t)/norm;
2514       Float_t diff = hPrior->GetBinContent(t) - newValue;
2515       chi2 += (Double_t) diff * diff;
2516
2517       hPrior->SetBinContent(t, newValue);
2518     }
2519
2520     printf("Chi2 of %d iteration = %.10f\n", i, chi2);
2521
2522     //if (i % 5 == 0)
2523       DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2524
2525     delete hTrueSmoothed;
2526   } // end of iterations
2527
2528   DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2529 }
2530
2531 //____________________________________________________________________
2532 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
2533 {
2534   //
2535   // correct spectrum using a simple Gaussian approach, that is model-dependent
2536   //
2537
2538   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2539   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2540
2541   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
2542   //normalize ESD
2543   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2544
2545   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
2546   correction->SetTitle("GaussianMean");
2547
2548   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
2549   correction->SetTitle("GaussianWidth");
2550
2551   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2552   {
2553     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
2554     proj->Fit("gaus", "0Q");
2555     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
2556     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
2557
2558     continue;
2559
2560     // draw for debugging
2561     new TCanvas;
2562     proj->DrawCopy();
2563     proj->GetFunction("gaus")->DrawCopy("SAME");
2564   }
2565
2566   TH1* target = fMultiplicityESDCorrected[correlationID];
2567
2568   Int_t nBins = target->GetNbinsX()*10+1;
2569   Float_t* binning = new Float_t[nBins];
2570   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2571     for (Int_t j=0; j<10; ++j)
2572       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
2573
2574   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
2575
2576   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
2577
2578   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2579   {
2580     Float_t mean = correction->GetBinContent(i);
2581     Float_t width = correctionWidth->GetBinContent(i);
2582
2583     Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
2584     Int_t fillEnd   = fineBinned->FindBin(mean + width * 5);
2585     //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
2586
2587     for (Int_t j=fillBegin; j <= fillEnd; ++j)
2588     {
2589       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
2590     }
2591   }
2592
2593   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2594   {
2595     Float_t sum = 0;
2596     for (Int_t j=1; j<=10; ++j)
2597       sum += fineBinned->GetBinContent((i-1)*10 + j);
2598     target->SetBinContent(i, sum / 10);
2599   }
2600
2601   delete[] binning;
2602
2603   DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
2604 }
2605
2606 //____________________________________________________________________
2607 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
2608 {
2609   // runs the distribution given in inputMC through the response matrix identified by
2610   // correlationMap and produces a measured distribution
2611   // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
2612   // if normalized is set, inputMC is expected to be normalized to the bin width
2613
2614   if (!inputMC)
2615     return 0;
2616
2617   if (correlationMap < 0 || correlationMap >= kCorrHists)
2618     return 0;
2619
2620   // empty under/overflow bins in x, otherwise Project3D takes them into account
2621   TH3* hist = fCorrelation[correlationMap];
2622   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
2623   {
2624     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
2625     {
2626       hist->SetBinContent(0, y, z, 0);
2627       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2628     }
2629   }
2630
2631   TH2* corr = (TH2*) hist->Project3D("zy");
2632   //corr->Rebin2D(2, 1);
2633   corr->Sumw2();
2634
2635   // normalize correction for given nPart
2636   for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
2637   {
2638     Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
2639     if (sum <= 0)
2640       continue;
2641
2642     for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
2643     {
2644       // npart sum to 1
2645       corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
2646       corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
2647     }
2648   }
2649
2650   TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
2651   target->Reset();
2652
2653   for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
2654   {
2655     Float_t measured = 0;
2656     Float_t error = 0;
2657
2658     for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
2659     {
2660       Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
2661
2662       measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
2663       error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
2664     }
2665
2666     //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
2667
2668     target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
2669     target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
2670   }
2671
2672   return target;
2673 }
2674
2675 //____________________________________________________________________
2676 void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
2677 {
2678   // uses the given function to fill the input MC histogram and generates from that
2679   // the measured histogram by applying the response matrix
2680   // this can be used to evaluate if the methods work indepedently of the input
2681   // distribution
2682   // WARNING does not respect the vertex distribution, just fills central vertex bin
2683
2684   if (!inputMC)
2685     return;
2686
2687   if (id < 0 || id >= kESDHists)
2688     return;
2689
2690   TH2F* mc = fMultiplicityVtx[id];
2691
2692   mc->Reset();
2693
2694   for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
2695   {
2696     mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
2697     mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, 0);
2698   }
2699
2700   //new TCanvas;
2701   //mc->Draw("COLZ");
2702
2703   TH1* proj = mc->ProjectionY();
2704
2705   TString tmp(fMultiplicityESD[id]->GetName());
2706   delete fMultiplicityESD[id];
2707   fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
2708   fMultiplicityESD[id]->SetName(tmp);
2709 }