CORRFW/AliCFUnfolding.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 //---------------------------------------------------------------------//
  17 //                                                                     //
  18 // AliCFUnfolding Class                                                //
  19 // Class to handle general unfolding procedure                         //
  20 // For the moment only bayesian unfolding is supported                 //
  21 // The next steps are to add chi2 minimisation and weighting methods   //
  22 //                                                                     //
  23 //                                                                     //
  24 //                                                                     //
  25 // Use :                                                               //
  26 // -------                                                             //
  27 // The Bayesian unfolding consists of several iterations.              //
  28 // At each iteration, an inverse response matrix is calculated, given  //
  29 // the measured spectrum, the a priori (guessed) spectrum,             //
  30 // the efficiency spectrum and the response matrix.                    //
  31 //                                                                     //
  32 // Then at each iteration, the unfolded spectrum is calculated using   //
  33 // the inverse response : the goal is to get an unfolded spectrum      //
  34 // similar (according to some criterion) to the a priori one.          //
  35 // If the difference is too big, another iteration is performed :      //
  36 // the a priori spectrum is updated to the unfolded one from the       //
  37 // previous iteration, and so on so forth, until the maximum number    //
  38 // of iterations or the similarity criterion is reached.               //
  39 //                                                                     //
  40 // Chi2 calculation became absolute with the correlated error          //
  41 // calculation.                                                        //
  42 // Errors on the unfolded distribution are not known until the end     //
  43 // Use the convergence criterion instead                               //
  44 //                                                                     //
  45 // Currently the user has to define the max. number of iterations      //
  46 // (::SetMaxNumberOfIterations)                                        //
  47 // and                                                                 //
  48 //     - the chi2 below which the procedure will stop                  //
  49 // (::SetMaxChi2 or ::SetMaxChi2PerDOF)   (OBSOLETE)                   //
  50 //     - the convergence criterion below which the procedure will stop //
  51 // SetMaxConvergencePerDOF(Double_t val);                              //
  52 //                                                                     //
  53 // Correlated error calculation can be activated by using:             //
  54 // SetUseCorrelatedErrors(Bool_t b) in combination with convergence    //
  55 // criterion                                                           //
  56 // Documentation about correlated error calculation method can be      //
  57 // found in AliCFUnfolding::CalculateCorrelatedErrors()                //
  58 // Author: marta.verweij@cern.ch                                       //
  59 //                                                                     //
  60 // An optional possibility is to smooth the unfolded spectrum at the   //
  61 // end of each iteration, either using a fit function                  //
  62 // (only if #dimensions <=3)                                           //
  63 // or a simple averaging using the neighbouring bins values.           //
  64 // This is possible calling the function ::UseSmoothing                //
  65 // If no argument is passed to this function, then the second option   //
  66 // is used.                                                            //
  67 //                                                                     //
  68 // IMPORTANT:                                                          //
  69 //-----------                                                          //
  70 // With this approach, the efficiency map must be calculated           //
  71 // with *simulated* values only, otherwise the method won't work.      //
  72 //                                                                     //
  73 // ex: efficiency(bin_pt) = number_rec(bin_pt) / number_sim(bin_pt)    //
  74 //                                                                     //
  75 // the pt bin "bin_pt" must always be the same in both the efficiency  //
  76 // numerator and denominator.                                          //
  77 // This is why the efficiency map has to be created by a method        //
  78 // from which both reconstructed and simulated values are accessible   //
  79 // simultaneously.                                                     //
  80 //                                                                     //
  81 //                                                                     //
  82 //---------------------------------------------------------------------//
  83 // Author : renaud.vernet@cern.ch                                      //
  84 //---------------------------------------------------------------------//
  85
  86
  87 #include "AliCFUnfolding.h"
  88 #include "TMath.h"
  89 #include "TAxis.h"
  90 #include "TF1.h"
  91 #include "TH1D.h"
  92 #include "TH2D.h"
  93 #include "TH3D.h"
  94 #include "TRandom3.h"
  95
  96 ClassImp(AliCFUnfolding)
  97
  98 //______________________________________________________________
  99
 100 AliCFUnfolding::AliCFUnfolding() :
 101   TNamed(),
 102   fResponseOrig(0x0),
 103   fPriorOrig(0x0),
 104   fEfficiencyOrig(0x0),
 105   fMeasuredOrig(0x0),
 106   fMaxNumIterations(0),
 107   fNVariables(0),
 108   fUseSmoothing(kFALSE),
 109   fSmoothFunction(0x0),
 110   fSmoothOption("iremn"),
 111   fMaxConvergence(0),
 112   fNRandomIterations(0),
 113   fResponse(0x0),
 114   fPrior(0x0),
 115   fEfficiency(0x0),
 116   fMeasured(0x0),
 117   fInverseResponse(0x0),
 118   fMeasuredEstimate(0x0),
 119   fConditional(0x0),
 120   fUnfolded(0x0),
 121   fUnfoldedFinal(0x0),
 122   fCoordinates2N(0x0),
 123   fCoordinatesN_M(0x0),
 124   fCoordinatesN_T(0x0),
 125   fRandomResponse(0x0),
 126   fRandomEfficiency(0x0),
 127   fRandomMeasured(0x0),
 128   fRandom3(0x0),
 129   fDeltaUnfoldedP(0x0),
 130   fDeltaUnfoldedN(0x0),
 131   fNCalcCorrErrors(0),
 132   fRandomSeed(0)
 133 {
 134   //
 135   // default constructor
 136   //
 137 }
 138
 139 //______________________________________________________________
 140
 141 AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar,
 142                                const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior ,
 143                                Double_t maxConvergencePerDOF, UInt_t randomSeed, Int_t maxNumIterations
 144                                ) :
 145   TNamed(name,title),
 146   fResponseOrig((THnSparse*)response->Clone()),
 147   fPriorOrig(0x0),
 148   fEfficiencyOrig((THnSparse*)efficiency->Clone()),
 149   fMeasuredOrig((THnSparse*)measured->Clone()),
 150   fMaxNumIterations(maxNumIterations),
 151   fNVariables(nVar),
 152   fUseSmoothing(kFALSE),
 153   fSmoothFunction(0x0),
 154   fSmoothOption("iremn"),
 155   fMaxConvergence(0),
 156   fNRandomIterations(maxNumIterations),
 157   fResponse((THnSparse*)response->Clone()),
 158   fPrior(0x0),
 159   fEfficiency((THnSparse*)efficiency->Clone()),
 160   fMeasured((THnSparse*)measured->Clone()),
 161   fInverseResponse(0x0),
 162   fMeasuredEstimate(0x0),
 163   fConditional(0x0),
 164   fUnfolded(0x0),
 165   fUnfoldedFinal(0x0),
 166   fCoordinates2N(0x0),
 167   fCoordinatesN_M(0x0),
 168   fCoordinatesN_T(0x0),
 169   fRandomResponse((THnSparse*)response->Clone()),
 170   fRandomEfficiency((THnSparse*)efficiency->Clone()),
 171   fRandomMeasured((THnSparse*)measured->Clone()),
 172   fRandom3(0x0),
 173   fDeltaUnfoldedP(0x0),
 174   fDeltaUnfoldedN(0x0),
 175   fNCalcCorrErrors(0),
 176   fRandomSeed(randomSeed)
 177 {
 178   //
 179   // named constructor
 180   //
 181
 182   AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
 183
 184   if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
 185   else {
 186     fPrior = (THnSparse*) prior->Clone();
 187     fPriorOrig = (THnSparse*)fPrior->Clone();
 188     if (fPrior->GetNdimensions() != fNVariables)
 189       AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
 190   }
 191
 192   if (fEfficiency->GetNdimensions() != fNVariables)
 193     AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions()));
 194   if (fMeasured->GetNdimensions() != fNVariables)
 195     AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions()));
 196   if (fResponse->GetNdimensions() != 2*fNVariables)
 197     AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions()));
 198
 199
 200   for (Int_t iVar=0; iVar<fNVariables; iVar++) {
 201     AliInfo(Form("prior      matrix has %d bins in dimension %d",fPrior     ->GetAxis(iVar)->GetNbins(),iVar));
 202     AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
 203     AliInfo(Form("measured   matrix has %d bins in dimension %d",fMeasured  ->GetAxis(iVar)->GetNbins(),iVar));
 204   }
 205
 206   fRandomResponse  ->SetTitle("Randomized response matrix");
 207   fRandomEfficiency->SetTitle("Randomized efficiency");
 208   fRandomMeasured  ->SetTitle("Randomized measured");
 209   SetMaxConvergencePerDOF(maxConvergencePerDOF)  ;
 210   Init();
 211 }
 212
 213 //______________________________________________________________
 214
 215 AliCFUnfolding::~AliCFUnfolding() {
 216   //
 217   // destructor
 218   //
 219   if (fResponse)           delete fResponse;
 220   if (fPrior)              delete fPrior;
 221   if (fEfficiency)         delete fEfficiency;
 222   if (fEfficiencyOrig)     delete fEfficiencyOrig;
 223   if (fMeasured)           delete fMeasured;
 224   if (fMeasuredOrig)       delete fMeasuredOrig;
 225   if (fSmoothFunction)     delete fSmoothFunction;
 226   if (fPriorOrig)          delete fPriorOrig;
 227   if (fInverseResponse)    delete fInverseResponse;
 228   if (fMeasuredEstimate)   delete fMeasuredEstimate;
 229   if (fConditional)        delete fConditional;
 230   if (fUnfolded)           delete fUnfolded;
 231   if (fUnfoldedFinal)      delete fUnfoldedFinal;
 232   if (fCoordinates2N)      delete [] fCoordinates2N;
 233   if (fCoordinatesN_M)     delete [] fCoordinatesN_M;
 234   if (fCoordinatesN_T)     delete [] fCoordinatesN_T;
 235   if (fRandomResponse)     delete fRandomResponse;
 236   if (fRandomEfficiency)   delete fRandomEfficiency;
 237   if (fRandomMeasured)     delete fRandomMeasured;
 238   if (fRandom3)            delete fRandom3;
 239   if (fDeltaUnfoldedP)     delete fDeltaUnfoldedP;
 240   if (fDeltaUnfoldedN)     delete fDeltaUnfoldedN;
 241 }
 242
 243 //______________________________________________________________
 244
 245 void AliCFUnfolding::Init() {
 246   //
 247   // initialisation function : creates internal settings
 248   //
 249
 250   fRandom3 = new TRandom3(fRandomSeed);
 251
 252   fCoordinates2N  = new Int_t[2*fNVariables];
 253   fCoordinatesN_M = new Int_t[fNVariables];
 254   fCoordinatesN_T = new Int_t[fNVariables];
 255
 256   // create the matrix of conditional probabilities P(M|T)
 257   CreateConditional(); //done only once at initialization
 258
 259   // create the frame of the inverse response matrix
 260   fInverseResponse  = (THnSparse*) fResponse->Clone();
 261   // create the frame of the unfolded spectrum
 262   fUnfolded = (THnSparse*) fPrior->Clone();
 263   fUnfolded->SetTitle("Unfolded");
 264   // create the frame of the measurement estimate spectrum
 265   fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
 266
 267   // create the frame of the delta profiles
 268   fDeltaUnfoldedP = (THnSparse*)fPrior->Clone();
 269   fDeltaUnfoldedP->SetTitle("#Delta unfolded");
 270   fDeltaUnfoldedP->Reset();
 271   fDeltaUnfoldedN = (THnSparse*)fPrior->Clone();
 272   fDeltaUnfoldedN->SetTitle("");
 273   fDeltaUnfoldedN->Reset();
 274 }
 275
 276
 277 //______________________________________________________________
 278
 279 void AliCFUnfolding::CreateEstMeasured() {
 280   //
 281   // This function creates a estimate (M) of the reconstructed spectrum
 282   // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND)
 283   //
 284   // --> P(M) = SUM   { P(M|T)    * P(T) }
 285   // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)}
 286   //
 287   // This is needed to calculate the inverse response matrix
 288   //
 289
 290
 291   // clean the measured estimate spectrum
 292   fMeasuredEstimate->Reset();
 293
 294   THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
 295   priorTimesEff->Multiply(fEfficiency);
 296
 297   // fill it
 298   for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
 299     Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
 300     GetCoordinates();
 301     Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
 302     Double_t fill = conditionalValue * priorTimesEffValue ;
 303
 304     if (fill>0.) {
 305       fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
 306       fMeasuredEstimate->SetBinError(fCoordinatesN_M,0.);
 307     }
 308   }
 309   delete priorTimesEff ;
 310 }
 311
 312 //______________________________________________________________
 313
 314 void AliCFUnfolding::CreateInvResponse() {
 315   //
 316   // Creates the inverse response matrix (INV) with Bayesian method
 317   //  : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
 318   //
 319   // --> P(T|M)   = P(M|T)    * P(T) * eff(T) / SUM   { P(M|T)    * P(T) }
 320   // --> INV(i,j) = COND(i,j) * T(j) * E(j)   / SUM_k { COND(i,k) * T(k) }
 321   //
 322
 323   THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
 324   priorTimesEff->Multiply(fEfficiency);
 325
 326   for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
 327     Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
 328     GetCoordinates();
 329     Double_t estMeasuredValue   = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
 330     Double_t priorTimesEffValue = priorTimesEff    ->GetBinContent(fCoordinatesN_T);
 331     Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
 332     if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
 333       fInverseResponse->SetBinContent(fCoordinates2N,fill);
 334       fInverseResponse->SetBinError  (fCoordinates2N,0.);
 335     }
 336   }
 337   delete priorTimesEff ;
 338 }
 339
 340 //______________________________________________________________
 341
 342 void AliCFUnfolding::Unfold() {
 343   //
 344   // Main routine called by the user :
 345   // it calculates the unfolded spectrum from the response matrix, measured spectrum and efficiency
 346   // several iterations are performed until a reasonable chi2 or convergence criterion is reached
 347   //
 348
 349   Int_t iIterBayes     = 0 ;
 350   Double_t convergence = 0.;
 351
 352   for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
 353
 354     CreateEstMeasured(); // create measured estimate from prior
 355     CreateInvResponse(); // create inverse response  from prior
 356     CreateUnfolded();    // create unfoled spectrum  from measured and inverse response
 357
 358     convergence = GetConvergence();
 359     AliDebug(0,Form("convergence at iteration %d is %e",iIterBayes,convergence));
 360
 361     if (fMaxConvergence>0. && convergence<fMaxConvergence && fNCalcCorrErrors == 0) {
 362       fNRandomIterations = iIterBayes;
 363       AliDebug(0,Form("convergence is met at iteration %d",iIterBayes));
 364       break;
 365     }
 366
 367     if (fUseSmoothing) {
 368       if (Smooth()) {
 369         AliError("Couldn't smooth the unfolded spectrum!!");
 370         if (fNCalcCorrErrors>0) {
 371           AliInfo(Form("=======================\nUnfold of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
 372         }
 373         else {
 374           AliInfo(Form("\n\n=======================\nFinish at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
 375         }
 376         return;
 377       }
 378     }
 379
 380     // update the prior distribution
 381     if (fPrior) delete fPrior ;
 382     fPrior = (THnSparse*)fUnfolded->Clone() ;
 383     fPrior->SetTitle("Prior");
 384
 385   } // end bayes iteration
 386
 387   if (fNCalcCorrErrors==0) fUnfoldedFinal = (THnSparse*) fUnfolded->Clone() ;
 388
 389   //
 390   //for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) AliDebug(2,Form("%e\n",fUnfoldedFinal->GetBinError(iBin)));
 391   //
 392
 393   if (fNCalcCorrErrors == 0) {
 394     AliInfo("\n================================================\nFinished bayes iteration, now calculating errors...\n================================================\n");
 395     fNCalcCorrErrors = 1;
 396     CalculateCorrelatedErrors();
 397   }
 398
 399   if (fNCalcCorrErrors >1 ) {
 400     AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
 401   }
 402   else if(fNCalcCorrErrors>0) {
 403     AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
 404   }
 405 }
 406
 407 //______________________________________________________________
 408
 409 void AliCFUnfolding::CreateUnfolded() {
 410   //
 411   // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
 412   // We have P(T) = SUM   { P(T|M)   * P(M) }
 413   //   -->   T(i) = SUM_k { INV(i,k) * M(k) }
 414   //
 415
 416
 417   // clear the unfolded spectrum
 418   // if in the process of error calculation, the random unfolded spectrum is created
 419   // otherwise the normal unfolded spectrum is created
 420
 421   fUnfolded->Reset();
 422
 423   for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
 424     Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
 425     GetCoordinates();
 426     Double_t effValue      = fEfficiency->GetBinContent(fCoordinatesN_T);
 427     Double_t measuredValue = fMeasured  ->GetBinContent(fCoordinatesN_M);
 428     Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
 429
 430     if (fill>0.) {
 431       // set errors to zero
 432       // true errors will be filled afterwards
 433       Double_t err = 0.;
 434       fUnfolded->SetBinError  (fCoordinatesN_T,err);
 435       fUnfolded->AddBinContent(fCoordinatesN_T,fill);
 436     }
 437   }
 438 }
 439
 440 //______________________________________________________________
 441
 442 void AliCFUnfolding::CalculateCorrelatedErrors() {
 443
 444   // Step 1: Create randomized distribution (fRandomXXXX) of each bin of
 445   //         the measured spectrum to calculate correlated errors.
 446   //         Poisson statistics: mean = measured value of bin
 447   // Step 2: Unfold randomized distribution
 448   // Step 3: Store difference of unfolded spectrum from measured distribution and
 449   //         unfolded distribution from randomized distribution
 450   //         -> fDeltaUnfoldedP (TProfile with option "S")
 451   // Step 4: Repeat Step 1-3 several times (fNRandomIterations)
 452   // Step 5: The spread of fDeltaUnfoldedP for each bin is the error on the unfolded spectrum of that specific bin
 453
 454
 455   //Do fNRandomIterations = bayes iterations performed
 456   for (int i=0; i<fNRandomIterations; i++) {
 457
 458     // reset prior to original one
 459     if (fPrior) delete fPrior ;
 460     fPrior = (THnSparse*) fPriorOrig->Clone();
 461
 462     // create randomized distribution and stick measured spectrum to it
 463     CreateRandomizedDist();
 464
 465     if (fResponse) delete fResponse ;
 466     fResponse = (THnSparse*) fRandomResponse->Clone();
 467     fResponse->SetTitle("Response");
 468
 469     if (fEfficiency) delete fEfficiency ;
 470     fEfficiency = (THnSparse*) fRandomEfficiency->Clone();
 471     fEfficiency->SetTitle("Efficiency");
 472
 473     if (fMeasured)   delete fMeasured   ;
 474     fMeasured = (THnSparse*) fRandomMeasured->Clone();
 475     fMeasured->SetTitle("Measured");
 476
 477     //unfold with randomized distributions
 478     Unfold();
 479     FillDeltaUnfoldedProfile();
 480   }
 481
 482   // Get statistical errors for final unfolded spectrum
 483   // ie. spread of each pt bin in fDeltaUnfoldedP
 484   Double_t sigma       = 0.;
 485   for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) {
 486     fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M);
 487     sigma = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M);
 488     //AliDebug(2,Form("filling error %e\n",sigma));
 489     fUnfoldedFinal->SetBinError(fCoordinatesN_M,sigma);
 490   }
 491
 492   // now errors are calculated
 493   fNCalcCorrErrors = 2;
 494 }
 495
 496 //______________________________________________________________
 497 void AliCFUnfolding::CreateRandomizedDist() {
 498   //
 499   // Create randomized dist from original measured distribution
 500   // This distribution is created several times, each time with a different random number
 501   //
 502
 503   for (Long_t iBin=0; iBin<fResponseOrig->GetNbins(); iBin++) {
 504     Double_t val = fResponseOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
 505     Double_t err = fResponseOrig->GetBinError(fCoordinatesN_M);        //used as sigma
 506     Double_t ran = fRandom3->Gaus(val,err);
 507     // random        = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
 508     fRandomResponse->SetBinContent(iBin,ran);
 509   }
 510   for (Long_t iBin=0; iBin<fEfficiencyOrig->GetNbins(); iBin++) {
 511     Double_t val = fEfficiencyOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
 512     Double_t err = fEfficiencyOrig->GetBinError(fCoordinatesN_M);        //used as sigma
 513     Double_t ran = fRandom3->Gaus(val,err);
 514     // random        = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
 515     fRandomEfficiency->SetBinContent(iBin,ran);
 516   }
 517   for (Long_t iBin=0; iBin<fMeasuredOrig->GetNbins(); iBin++) {
 518     Double_t val = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
 519     Double_t err = fMeasuredOrig->GetBinError(fCoordinatesN_M);        //used as sigma
 520     Double_t ran = fRandom3->Gaus(val,err);
 521     // random        = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
 522     fRandomMeasured->SetBinContent(iBin,ran);
 523   }
 524 }
 525
 526 //______________________________________________________________
 527 void AliCFUnfolding::FillDeltaUnfoldedProfile() {
 528   //
 529   // Store difference of unfolded spectrum from measured distribution and unfolded spectrum from randomized distribution
 530   // The delta profile has been set to a THnSparse to handle N dimension
 531   // The THnSparse contains in each bin the mean value and spread of the difference
 532   // This function updates the profile wrt to its previous mean and error
 533   // The relation between iterations (n+1) and n is as follows :
 534   //  mean_{n+1} = (n*mean_n + value_{n+1}) / (n+1)
 535   // sigma_{n+1} = sqrt { 1/(n+1) * [ n*sigma_n^2 + (n^2+n)*(mean_{n+1}-mean_n)^2 ] }    (can this be optimized?)
 536
 537   for (Long_t iBin=0; iBin<fUnfolded->GetNbins(); iBin++) {
 538     Double_t deltaInBin   = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M) - fUnfolded->GetBinContent(iBin);
 539     Double_t entriesInBin = fDeltaUnfoldedN->GetBinContent(fCoordinatesN_M);
 540     //AliDebug(2,Form("%e %e ==> delta = %e\n",fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M),fUnfolded->GetBinContent(iBin),deltaInBin));
 541
 542     Double_t mean_n = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M) ;
 543     Double_t mean_nplus1 = mean_n ;
 544     mean_nplus1 *= entriesInBin ;
 545     mean_nplus1 += deltaInBin ;
 546     mean_nplus1 /= (entriesInBin+1) ;
 547
 548     Double_t sigma = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M) ;
 549     sigma *= sigma ;
 550     sigma *= entriesInBin ;
 551     sigma += ( (entriesInBin*entriesInBin+entriesInBin) * TMath::Power(mean_nplus1 - mean_n,2) ) ;
 552     sigma /= (entriesInBin+1) ;
 553     sigma = TMath::Sqrt(sigma) ;
 554
 555     //AliDebug(2,Form("sigma = %e\n",sigma));
 556
 557     fDeltaUnfoldedP->SetBinContent(fCoordinatesN_M,mean_nplus1) ;
 558     fDeltaUnfoldedP->SetBinError  (fCoordinatesN_M,sigma) ;
 559     fDeltaUnfoldedN->SetBinContent(fCoordinatesN_M,entriesInBin+1);
 560   }
 561 }
 562
 563 //______________________________________________________________
 564
 565 void AliCFUnfolding::GetCoordinates() {
 566   //
 567   // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
 568   //
 569   for (Int_t i = 0; i<fNVariables ; i++) {
 570     fCoordinatesN_M[i] = fCoordinates2N[i];
 571     fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
 572   }
 573 }
 574
 575 //______________________________________________________________
 576
 577 void AliCFUnfolding::CreateConditional() {
 578   //
 579   // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
 580   //
 581   //  --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
 582   //
 583
 584   fConditional = (THnSparse*) fResponse->Clone();  // output of this function
 585
 586   Int_t* dim = new Int_t [fNVariables];
 587   for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
 588
 589   THnSparse* responseInT = fConditional->Projection(fNVariables,dim,"E");     // output denominator :
 590                                                                               // projection of the response matrix on the TRUE axis
 591   delete [] dim;
 592
 593   // fill the conditional probability matrix
 594   for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
 595     Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
 596     GetCoordinates();
 597     Double_t projValue = responseInT->GetBinContent(fCoordinatesN_T);
 598
 599     Double_t fill = responseValue / projValue ;
 600     if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
 601       fConditional->SetBinContent(fCoordinates2N,fill);
 602       Double_t err = 0.;
 603       fConditional->SetBinError  (fCoordinates2N,err);
 604     }
 605   }
 606   delete responseInT ;
 607 }
 608 //______________________________________________________________
 609
 610 Int_t AliCFUnfolding::GetDOF() {
 611   //
 612   // number of dof = number of bins
 613   //
 614
 615   Int_t nDOF = 1 ;
 616   for (Int_t iDim=0; iDim<fNVariables; iDim++) {
 617     nDOF *= fPrior->GetAxis(iDim)->GetNbins();
 618   }
 619   AliDebug(0,Form("Number of degrees of freedom = %d",nDOF));
 620   return nDOF;
 621 }
 622
 623 //______________________________________________________________
 624
 625 Double_t AliCFUnfolding::GetChi2() {
 626   //
 627   // Returns the chi2 between unfolded and a priori spectrum
 628   // This function became absolute with the correlated error calculation.
 629   // Errors on the unfolded distribution are not known until the end
 630   // Use the convergence criterion instead
 631   //
 632
 633   Double_t chi2      = 0. ;
 634   Double_t error_unf = 0.;
 635   for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
 636     Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
 637     error_unf  = fUnfolded->GetBinError(fCoordinatesN_T);
 638     chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ;
 639   }
 640   return chi2;
 641 }
 642
 643 //______________________________________________________________
 644
 645 Double_t AliCFUnfolding::GetConvergence() {
 646   //
 647   // Returns convergence criterion = \sum_t ((U_t^{n-1}-U_t^n)/U_t^{n-1})^2
 648   // U is unfolded spectrum, t is the bin, n = current, n-1 = previous
 649   //
 650   Double_t convergence = 0.;
 651   Double_t priorValue  = 0.;
 652   Double_t currentValue = 0.;
 653   for (Long_t iBin=0; iBin < fPrior->GetNbins(); iBin++) {
 654     priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
 655     currentValue = fUnfolded->GetBinContent(fCoordinatesN_T);
 656
 657     if (priorValue > 0.)
 658       convergence += ((priorValue-currentValue)/priorValue)*((priorValue-currentValue)/priorValue);
 659     else
 660       AliWarning(Form("priorValue = %f. Adding 0 to convergence criterion.",priorValue));
 661   }
 662   return convergence;
 663 }
 664
 665 //______________________________________________________________
 666
 667 void AliCFUnfolding::SetMaxConvergencePerDOF(Double_t val) {
 668   //
 669   // Max. convergence criterion per degree of freedom : user setting
 670   // convergence criterion = DOF*val; DOF = number of bins
 671   // In Jan-Fiete's multiplicity note: Convergence criterion = DOF*0.001^2
 672   //
 673
 674   Int_t nDOF = GetDOF() ;
 675   fMaxConvergence = val * nDOF ;
 676   AliInfo(Form("MaxConvergence = %e. Number of degrees of freedom = %d",fMaxConvergence,nDOF));
 677 }
 678
 679 //______________________________________________________________
 680
 681 Short_t AliCFUnfolding::Smooth() {
 682   //
 683   // Smoothes the unfolded spectrum
 684   //
 685   // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
 686   // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn
 687   //
 688
 689   if (fSmoothFunction) {
 690     AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
 691     return SmoothUsingFunction();
 692   }
 693   else return SmoothUsingNeighbours(fUnfolded);
 694 }
 695
 696 //______________________________________________________________
 697
 698 Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) {
 699   //
 700   // Smoothes the unfolded spectrum using neighouring bins
 701   //
 702
 703   Int_t  const nDimensions = hist->GetNdimensions() ;
 704   Int_t* coordinates = new Int_t[nDimensions];
 705
 706   Int_t* numBins = new Int_t[nDimensions];
 707   for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins();
 708
 709   //need a copy because hist will be updated during the loop, and this creates problems
 710   THnSparse* copy = (THnSparse*)hist->Clone();
 711
 712   for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
 713     Double_t content = copy->GetBinContent(iBin,coordinates);
 714     Double_t error2  = TMath::Power(copy->GetBinError(iBin),2);
 715
 716     // skip the under/overflow bins...
 717     Bool_t isOutside = kFALSE ;
 718     for (Int_t iVar=0; iVar<nDimensions; iVar++) {
 719       if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) {
 720         isOutside=kTRUE;
 721         break;
 722       }
 723     }
 724     if (isOutside) continue;
 725
 726     Int_t neighbours = 0; // number of neighbours to average with
 727
 728     for (Int_t iVar=0; iVar<nDimensions; iVar++) {
 729       if (coordinates[iVar] > 1) { // must not be on low edge border
 730         coordinates[iVar]-- ; //get lower neighbouring bin
 731         content += copy->GetBinContent(coordinates);
 732         error2  += TMath::Power(copy->GetBinError(coordinates),2);
 733         neighbours++;
 734         coordinates[iVar]++ ; //back to initial coordinate
 735       }
 736       if (coordinates[iVar] < numBins[iVar]) { // must not be on up edge border
 737         coordinates[iVar]++ ; //get upper neighbouring bin
 738         content += copy->GetBinContent(coordinates);
 739         error2  += TMath::Power(copy->GetBinError(coordinates),2);
 740         neighbours++;
 741         coordinates[iVar]-- ; //back to initial coordinate
 742       }
 743     }
 744     // make an average
 745     hist->SetBinContent(coordinates,content/(1.+neighbours));
 746     hist->SetBinError  (coordinates,TMath::Sqrt(error2)/(1.+neighbours));
 747   }
 748   delete [] numBins;
 749   delete [] coordinates ;
 750   delete copy;
 751   return 0;
 752 }
 753
 754 //______________________________________________________________
 755
 756 Short_t AliCFUnfolding::SmoothUsingFunction() {
 757   //
 758   // Fits the unfolded spectrum using the function fSmoothFunction
 759   //
 760
 761   AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
 762
 763   for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
 764
 765   Int_t fitResult = 0;
 766
 767   switch (fNVariables) {
 768   case 1 : fitResult = fUnfolded->Projection(0)    ->Fit(fSmoothFunction,fSmoothOption); break;
 769   case 2 : fitResult = fUnfolded->Projection(1,0)  ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
 770   case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
 771   default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
 772   }
 773
 774   if (fitResult != 0) {
 775     AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
 776     return 1;
 777   }
 778
 779   Int_t nDim = fNVariables;
 780   Int_t* bins = new Int_t[nDim]; // number of bins for each variable
 781   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
 782
 783   for (Int_t iVar=0; iVar<nDim; iVar++) {
 784     bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
 785     nBins *= bins[iVar];
 786   }
 787
 788   Int_t *bin  = new Int_t[nDim];    // bin to fill the THnSparse (holding the bin coordinates)
 789   Double_t x[3] = {0,0,0} ;         // value in bin center (max dimension is 3 (TF3))
 790
 791   // loop on the bins and update of fUnfolded
 792   // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
 793   for (Long_t iBin=0; iBin<nBins; iBin++) {
 794     Long_t bin_tmp = iBin ;
 795     for (Int_t iVar=0; iVar<nDim; iVar++) {
 796       bin[iVar] = 1 + bin_tmp % bins[iVar] ;
 797       bin_tmp /= bins[iVar] ;
 798       x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
 799     }
 800     Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
 801     fUnfolded->SetBinError  (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin));
 802     fUnfolded->SetBinContent(bin,functionValue);
 803   }
 804   delete [] bins;
 805   delete [] bin ;
 806   return 0;
 807 }
 808
 809 //______________________________________________________________
 810
 811 void AliCFUnfolding::CreateFlatPrior() {
 812   //
 813   // Creates a flat prior distribution
 814   //
 815
 816   AliInfo("Creating a flat a priori distribution");
 817
 818   // create the frame of the THnSparse given (for example) the one from the efficiency map
 819   fPrior = (THnSparse*) fEfficiency->Clone();
 820   fPrior->SetTitle("Prior");
 821
 822   if (fNVariables != fPrior->GetNdimensions())
 823     AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
 824
 825   Int_t nDim = fNVariables;
 826   Int_t* bins = new Int_t[nDim]; // number of bins for each variable
 827   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
 828
 829   for (Int_t iVar=0; iVar<nDim; iVar++) {
 830     bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
 831     nBins *= bins[iVar];
 832   }
 833
 834   Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
 835
 836   // loop that sets 1 in each bin
 837   for (Long_t iBin=0; iBin<nBins; iBin++) {
 838     Long_t bin_tmp = iBin ;
 839     for (Int_t iVar=0; iVar<nDim; iVar++) {
 840       bin[iVar] = 1 + bin_tmp % bins[iVar] ;
 841       bin_tmp /= bins[iVar] ;
 842     }
 843     fPrior->SetBinContent(bin,1.); // put 1 everywhere
 844     fPrior->SetBinError  (bin,0.); // put 0 everywhere
 845   }
 846
 847   fPriorOrig = (THnSparse*)fPrior->Clone();
 848
 849   delete [] bin;
 850   delete [] bins;
 851 }