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