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