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 // For each iteration, the unfolded spectrum is calculated using       //
  32 // the inverse response : the goal is to get an unfolded spectrum      //
  33 // similar (according to some criterion) to the a priori one.          //
  34 // If the difference is too big, another iteration is performed :      //
  35 // the a priori spectrum is updated to the unfolded one from the       //
  36 // previous iteration, and so on so forth, until the maximum number    //
  37 // of iterations or the similarity criterion is reached.               //
  38 //                                                                     //
  39 // Currently the similarity criterion is the Chi2 between the a priori //
  40 // and the unfolded spectrum.                                          //
  41 //                                                                     //
  42 // Currently the user has to define the max. number of iterations      //
  43 // (::SetMaxNumberOfIterations)                                        //
  44 // and the chi2 below which the procedure will stop                    //
  45 // (::SetMaxChi2 or ::SetMaxChi2PerDOF)                                //
  46 //                                                                     //
  47 // An optional possibility is to smooth the unfolded spectrum at the   //
  48 // end of each iteration, either using a fit function                  //
  49 // (only if #dimensions <=3)                                           //
  50 // or a simple averaging using the neighbouring bins values.           //
  51 // This is possible calling the function ::UseSmoothing                //
  52 // If no argument is passed to this function, then the second option   //
  53 // is used.                                                            //
  54 //                                                                     //
  55 //---------------------------------------------------------------------//
  56 // Author : renaud.vernet@cern.ch                                      //
  57 //---------------------------------------------------------------------//
  58
  59
  60 #include "AliCFUnfolding.h"
  61 #include "TMath.h"
  62 #include "TAxis.h"
  63 #include "AliLog.h"
  64 #include "TF1.h"
  65 #include "TH1D.h"
  66 #include "TH2D.h"
  67 #include "TH3D.h"
  68
  69 ClassImp(AliCFUnfolding)
  70
  71 //______________________________________________________________
  72
  73 AliCFUnfolding::AliCFUnfolding() :
  74   TNamed(),
  75   fResponse(0x0),
  76   fPrior(0x0),
  77   fEfficiency(0x0),
  78   fMeasured(0x0),
  79   fMaxNumIterations(0),
  80   fNVariables(0),
  81   fMaxChi2(0),
  82   fUseSmoothing(kFALSE),
  83   fSmoothFunction(0x0),
  84   fSmoothOption(""),
  85   fOriginalPrior(0x0),
  86   fInverseResponse(0x0),
  87   fMeasuredEstimate(0x0),
  88   fConditional(0x0),
  89   fProjResponseInT(0x0),
  90   fUnfolded(0x0),
  91   fCoordinates2N(0x0),
  92   fCoordinatesN_M(0x0),
  93   fCoordinatesN_T(0x0)
  94 {
  95   //
  96   // default constructor
  97   //
  98 }
  99
 100 //______________________________________________________________
 101
 102 AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar,
 103                                const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior) :
 104   TNamed(name,title),
 105   fResponse((THnSparse*)response->Clone()),
 106   fPrior(0x0),
 107   fEfficiency((THnSparse*)efficiency->Clone()),
 108   fMeasured((THnSparse*)measured->Clone()),
 109   fMaxNumIterations(0),
 110   fNVariables(nVar),
 111   fMaxChi2(0),
 112   fUseSmoothing(kFALSE),
 113   fSmoothFunction(0x0),
 114   fSmoothOption(""),
 115   fOriginalPrior(0x0),
 116   fInverseResponse(0x0),
 117   fMeasuredEstimate(0x0),
 118   fConditional(0x0),
 119   fProjResponseInT(0x0),
 120   fUnfolded(0x0),
 121   fCoordinates2N(0x0),
 122   fCoordinatesN_M(0x0),
 123   fCoordinatesN_T(0x0)
 124 {
 125   //
 126   // named constructor
 127   //
 128
 129   AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
 130
 131   if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
 132   else {
 133     fPrior = (THnSparse*) prior->Clone();
 134     fOriginalPrior = (THnSparse*)fPrior->Clone();
 135     if (fPrior->GetNdimensions() != fNVariables)
 136       AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
 137   }
 138
 139   if (fEfficiency->GetNdimensions() != fNVariables)
 140     AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions()));
 141   if (fMeasured->GetNdimensions() != fNVariables)
 142     AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions()));
 143   if (fResponse->GetNdimensions() != 2*fNVariables)
 144     AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions()));
 145
 146
 147   for (Int_t iVar=0; iVar<fNVariables; iVar++) {
 148     AliInfo(Form("prior      matrix has %d bins in dimension %d",fPrior     ->GetAxis(iVar)->GetNbins(),iVar));
 149     AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
 150     AliInfo(Form("measured   matrix has %d bins in dimension %d",fMeasured  ->GetAxis(iVar)->GetNbins(),iVar));
 151   }
 152   Init();
 153 }
 154
 155
 156 //______________________________________________________________
 157
 158 AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) :
 159   TNamed(c),
 160   fResponse((THnSparse*)c.fResponse->Clone()),
 161   fPrior((THnSparse*)c.fPrior->Clone()),
 162   fEfficiency((THnSparse*)c.fEfficiency->Clone()),
 163   fMeasured((THnSparse*)c.fMeasured->Clone()),
 164   fMaxNumIterations(c.fMaxNumIterations),
 165   fNVariables(c.fNVariables),
 166   fMaxChi2(c.fMaxChi2),
 167   fUseSmoothing(c.fUseSmoothing),
 168   fSmoothFunction((TF1*)c.fSmoothFunction->Clone()),
 169   fSmoothOption(fSmoothOption),
 170   fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()),
 171   fInverseResponse((THnSparse*)c.fInverseResponse->Clone()),
 172   fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()),
 173   fConditional((THnSparse*)c.fConditional->Clone()),
 174   fProjResponseInT((THnSparse*)c.fProjResponseInT->Clone()),
 175   fUnfolded((THnSparse*)c.fUnfolded->Clone()),
 176   fCoordinates2N(new Int_t(*c.fCoordinates2N)),
 177   fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)),
 178   fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T))
 179 {
 180   //
 181   // copy constructor
 182   //
 183 }
 184
 185 //______________________________________________________________
 186
 187 AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) {
 188   //
 189   // assignment operator
 190   //
 191
 192   if (this!=&c) {
 193     TNamed::operator=(c);
 194     fResponse = (THnSparse*)c.fResponse->Clone() ;
 195     fPrior = (THnSparse*)c.fPrior->Clone() ;
 196     fEfficiency = (THnSparse*)c.fEfficiency->Clone() ;
 197     fMeasured = (THnSparse*)c.fMeasured->Clone() ;
 198     fMaxNumIterations = c.fMaxNumIterations ;
 199     fNVariables = c.fNVariables ;
 200     fMaxChi2 = c.fMaxChi2 ;
 201     fUseSmoothing = c.fUseSmoothing ;
 202     fSmoothFunction = (TF1*)c.fSmoothFunction->Clone();
 203     fSmoothOption = c.fSmoothOption ;
 204     fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ;
 205     fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ;
 206     fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ;
 207     fConditional = (THnSparse*)c.fConditional->Clone() ;
 208     fProjResponseInT = (THnSparse*)c.fProjResponseInT->Clone() ;
 209     fUnfolded = (THnSparse*)c.fUnfolded->Clone() ;
 210     fCoordinates2N  = new Int_t(*c.fCoordinates2N)  ;
 211     fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ;
 212     fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ;
 213   }
 214   return *this;
 215 }
 216
 217 //______________________________________________________________
 218
 219 AliCFUnfolding::~AliCFUnfolding() {
 220   //
 221   // destructor
 222   //
 223   if (fResponse)           delete fResponse;
 224   if (fPrior)              delete fPrior;
 225   if (fEfficiency)         delete fEfficiency;
 226   if (fMeasured)           delete fMeasured;
 227   if (fSmoothFunction)     delete fSmoothFunction;
 228   if (fOriginalPrior)      delete fOriginalPrior;
 229   if (fInverseResponse)    delete fInverseResponse;
 230   if (fMeasuredEstimate)   delete fMeasuredEstimate;
 231   if (fConditional)        delete fConditional;
 232   if (fProjResponseInT)    delete fProjResponseInT;
 233   if (fCoordinates2N)      delete [] fCoordinates2N;
 234   if (fCoordinatesN_M)     delete [] fCoordinatesN_M;
 235   if (fCoordinatesN_T)     delete [] fCoordinatesN_T;
 236 }
 237
 238 //______________________________________________________________
 239
 240 void AliCFUnfolding::Init() {
 241   //
 242   // initialisation function : creates internal settings
 243   //
 244
 245   fCoordinates2N  = new Int_t[2*fNVariables];
 246   fCoordinatesN_M = new Int_t[fNVariables];
 247   fCoordinatesN_T = new Int_t[fNVariables];
 248
 249   // create the matrix of conditional probabilities P(M|T)
 250   CreateConditional();
 251
 252   // create the frame of the inverse response matrix
 253   fInverseResponse  = (THnSparse*) fResponse->Clone();
 254   // create the frame of the unfolded spectrum
 255   fUnfolded = (THnSparse*) fPrior->Clone();
 256   // create the frame of the measurement estimate spectrum
 257   fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
 258 }
 259
 260 //______________________________________________________________
 261
 262 void AliCFUnfolding::CreateEstMeasured() {
 263   //
 264   // This function creates a estimate (M) of the reconstructed spectrum
 265   // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND)
 266   //
 267   // --> P(M) = SUM   { P(M|T)    * P(T) }
 268   // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)}
 269   //
 270   // This is needed to calculate the inverse response matrix
 271   //
 272
 273
 274   // clean the measured estimate spectrum
 275   for (Long_t i=0; i<fMeasuredEstimate->GetNbins(); i++) {
 276     fMeasuredEstimate->GetBinContent(i,fCoordinatesN_M);
 277     fMeasuredEstimate->SetBinContent(fCoordinatesN_M,0.);
 278     fMeasuredEstimate->SetBinError  (fCoordinatesN_M,0.);
 279   }
 280
 281   THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
 282   priorTimesEff->Multiply(fEfficiency);
 283
 284   // fill it
 285   for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
 286     Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
 287     Double_t conditionalError = fConditional->GetBinError  (iBin);
 288     GetCoordinates();
 289     Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
 290     Double_t priorTimesEffError = priorTimesEff->GetBinError  (fCoordinatesN_T);
 291     Double_t fill = conditionalValue * priorTimesEffValue ;
 292
 293     if (fill>0.) {
 294       fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
 295
 296       // error calculation : gaussian error propagation (may be overestimated...)
 297       Double_t err2  = TMath::Power(fMeasuredEstimate->GetBinError(fCoordinatesN_M),2) ;
 298       err2 += TMath::Power(conditionalValue*priorTimesEffError,2) + TMath::Power(conditionalError*priorTimesEffValue,2) ;
 299       Double_t err = TMath::Sqrt(err2);
 300       fMeasuredEstimate->SetBinError(fCoordinatesN_M,err);
 301     }
 302   }
 303   delete priorTimesEff ;
 304 }
 305
 306 //______________________________________________________________
 307
 308 void AliCFUnfolding::CreateInvResponse() {
 309   //
 310   // Creates the inverse response matrix (INV) with Bayesian method
 311   //  : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
 312   //
 313   // --> P(T|M)   = P(M|T)    * P(T) * eff(T) / SUM   { P(M|T)    * P(T) }
 314   // --> INV(i,j) = COND(i,j) * T(j) * E(j)   / SUM_k { COND(i,k) * T(k) }
 315   //
 316
 317   THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
 318   priorTimesEff->Multiply(fEfficiency);
 319
 320   for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
 321     Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
 322     Double_t conditionalError = fConditional->GetBinError  (iBin);
 323     GetCoordinates();
 324     Double_t estMeasuredValue   = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
 325     Double_t estMeasuredError   = fMeasuredEstimate->GetBinError  (fCoordinatesN_M);
 326     Double_t priorTimesEffValue = priorTimesEff    ->GetBinContent(fCoordinatesN_T);
 327     Double_t priorTimesEffError = priorTimesEff    ->GetBinError  (fCoordinatesN_T);
 328     Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
 329     // error calculation : gaussian error propagation (may be overestimated...)
 330     Double_t err  = 0. ;
 331     if (estMeasuredValue>0.) {
 332       err = TMath::Sqrt( TMath::Power(conditionalError * priorTimesEffValue * estMeasuredValue ,2) +
 333                          TMath::Power(conditionalValue * priorTimesEffError * estMeasuredValue ,2) +
 334                          TMath::Power(conditionalValue * priorTimesEffValue * estMeasuredError ,2) )
 335         / TMath::Power(estMeasuredValue,2) ;
 336     }
 337     if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
 338       fInverseResponse->SetBinContent(fCoordinates2N,fill);
 339       fInverseResponse->SetBinError  (fCoordinates2N,err );
 340     }
 341   }
 342   delete priorTimesEff ;
 343 }
 344
 345 //______________________________________________________________
 346
 347 void AliCFUnfolding::Unfold() {
 348   //
 349   // Main routine called by the user :
 350   // it calculates the unfolded spectrum from the response matrix and the measured spectrum
 351   // several iterations are performed until a reasonable chi2 is reached
 352   //
 353
 354   Int_t iIterBayes=0 ;
 355   Double_t chi2=0 ;
 356
 357   for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
 358     CreateEstMeasured();
 359     CreateInvResponse();
 360     CreateUnfolded();
 361     chi2 = GetChi2();
 362     AliDebug(1,Form("Chi2 at iteration %d is %e",iIterBayes,chi2));
 363     if (fMaxChi2>0. && chi2<fMaxChi2) {
 364       break;
 365     }
 366     // update the prior distribution
 367     if (fUseSmoothing) {
 368       if (Smooth()) {
 369         AliError("Couldn't smooth the unfolded spectrum!!");
 370         AliInfo(Form("\n\n=======================\nFinished at iteration %d : Chi2 is %e and you required it to be < %e\n=======================\n\n",iIterBayes,chi2,fMaxChi2));
 371         return;
 372       }
 373     }
 374     fPrior = (THnSparse*)fUnfolded->Clone() ; // this should be changed (memory)
 375   }
 376   AliInfo(Form("\n\n=======================\nFinished at iteration %d : Chi2 is %e and you required it to be < %e\n=======================\n\n",iIterBayes,chi2,fMaxChi2));
 377 }
 378
 379 //______________________________________________________________
 380
 381 void AliCFUnfolding::CreateUnfolded() {
 382   //
 383   // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
 384   // We have P(T) = SUM   { P(T|M)   * P(M) }
 385   //   -->   T(i) = SUM_k { INV(i,k) * M(k) }
 386   //
 387
 388
 389   // clear the unfolded spectrum
 390   for (Long_t i=0; i<fUnfolded->GetNbins(); i++) {
 391     fUnfolded->GetBinContent(i,fCoordinatesN_T);
 392     fUnfolded->SetBinContent(fCoordinatesN_T,0.);
 393     fUnfolded->SetBinError  (fCoordinatesN_T,0.);
 394   }
 395
 396   for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
 397     Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
 398     Double_t invResponseError = fInverseResponse->GetBinError  (iBin);
 399     GetCoordinates();
 400     Double_t effValue      = fEfficiency->GetBinContent(fCoordinatesN_T);
 401     Double_t effError      = fEfficiency->GetBinError  (fCoordinatesN_T);
 402     Double_t measuredValue = fMeasured  ->GetBinContent(fCoordinatesN_M);
 403     Double_t measuredError = fMeasured  ->GetBinError  (fCoordinatesN_M);
 404     Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
 405
 406     if (fill>0.) {
 407       fUnfolded->AddBinContent(fCoordinatesN_T,fill);
 408
 409       // error calculation : gaussian error propagation (may be overestimated...)
 410       Double_t err2 = TMath::Power(fUnfolded->GetBinError(fCoordinatesN_T),2) ;
 411       err2 += TMath::Power(invResponseError * measuredValue * effValue,2) / TMath::Power(effValue,4) ;
 412       err2 += TMath::Power(invResponseValue * measuredError * effValue,2) / TMath::Power(effValue,4) ;
 413       err2 += TMath::Power(invResponseValue * measuredValue * effError,2) / TMath::Power(effValue,4) ;
 414       Double_t err = TMath::Sqrt(err2);
 415       fUnfolded->SetBinError(fCoordinatesN_T,err);
 416     }
 417   }
 418 }
 419
 420 //______________________________________________________________
 421
 422 void AliCFUnfolding::GetCoordinates() {
 423   //
 424   // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
 425   //
 426   for (Int_t i = 0; i<fNVariables ; i++) {
 427     fCoordinatesN_M[i] = fCoordinates2N[i];
 428     fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
 429   }
 430 }
 431
 432 //______________________________________________________________
 433
 434 void AliCFUnfolding::CreateConditional() {
 435   //
 436   // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
 437   //
 438   //  --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
 439   //
 440
 441   fConditional     = (THnSparse*) fResponse->Clone(); // output of this function
 442   fProjResponseInT = (THnSparse*) fPrior->Clone();    // output denominator :
 443                                                       // projection of the response matrix on the TRUE axis
 444   Int_t* dim = new Int_t [fNVariables];
 445   for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
 446   fProjResponseInT = fConditional->Projection(fNVariables,dim,"E"); //project
 447   delete [] dim;
 448
 449   // fill the conditional probability matrix
 450   for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
 451     Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
 452     Double_t responseError = fResponse->GetBinError  (iBin);
 453     GetCoordinates();
 454     Double_t projValue = fProjResponseInT->GetBinContent(fCoordinatesN_T);
 455     Double_t projError = fProjResponseInT->GetBinError  (fCoordinatesN_T);
 456
 457     Double_t fill = responseValue / projValue ;
 458     if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
 459       fConditional->SetBinContent(fCoordinates2N,fill);
 460       // gaussian error for the moment
 461       Double_t err2 = TMath::Power(responseError*projValue,2) + TMath::Power(responseValue*projError,2) ;
 462       Double_t err = TMath::Sqrt(err2);
 463       err /= TMath::Power(projValue,2) ;
 464       fConditional->SetBinError  (fCoordinates2N,err);
 465     }
 466   }
 467 }
 468
 469 //______________________________________________________________
 470
 471 Double_t AliCFUnfolding::GetChi2() {
 472   //
 473   // Returns the chi2 between unfolded and a priori spectrum
 474   //
 475
 476   Double_t chi2 = 0. ;
 477   for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
 478     Double_t priorValue = fPrior->GetBinContent(iBin);
 479 //     chi2 += (priorValue>0. ? TMath::Power(fUnfolded->GetBinContent(iBin) - priorValue,2) / priorValue : 0.) ;
 480     chi2 += (fUnfolded->GetBinError(iBin)>0. ? TMath::Power((fUnfolded->GetBinContent(iBin) - priorValue)/fUnfolded->GetBinError(iBin),2) / priorValue : 0.) ;
 481   }
 482   return chi2;
 483 }
 484
 485 //______________________________________________________________
 486
 487 void AliCFUnfolding::SetMaxChi2PerDOF(Double_t val) {
 488   //
 489   // Max. chi2 per degree of freedom : user setting
 490   //
 491
 492   Int_t nDOF = 1 ;
 493   for (Int_t iDim=0; iDim<fNVariables; iDim++) {
 494     nDOF *= fPrior->GetAxis(iDim)->GetNbins();
 495   }
 496   AliInfo(Form("Number of degrees of freedom = %d",nDOF));
 497   fMaxChi2 = val * nDOF ;
 498 }
 499
 500 //______________________________________________________________
 501
 502 Short_t AliCFUnfolding::Smooth() {
 503   //
 504   // Smoothes the unfolded spectrum
 505   //
 506   // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
 507   // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn
 508   //
 509
 510   if (fSmoothFunction) {
 511     AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
 512     return SmoothUsingFunction();
 513   }
 514   else return SmoothUsingNeighbours();
 515 }
 516
 517 //______________________________________________________________
 518
 519 Short_t AliCFUnfolding::SmoothUsingNeighbours() {
 520   //
 521   // Smoothes the unfolded spectrum using neighouring bins
 522   //
 523
 524   Int_t* numBins = new Int_t[fNVariables];
 525   for (Int_t iVar=0; iVar<fNVariables; iVar++) numBins[iVar]=fUnfolded->GetAxis(iVar)->GetNbins();
 526
 527   //need a copy because fUnfolded will be updated during the loop, and this creates problems
 528   THnSparse* copy = (THnSparse*)fUnfolded->Clone();
 529
 530   for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
 531     Double_t content = copy->GetBinContent(iBin,fCoordinatesN_T);
 532     Double_t error2  = TMath::Power(copy->GetBinError(iBin),2);
 533
 534     // skip the under/overflow bins...
 535     Bool_t isOutside = kFALSE ;
 536     for (Int_t iVar=0; iVar<fNVariables; iVar++) {
 537       if (fCoordinatesN_T[iVar]<1 || fCoordinatesN_T[iVar]>numBins[iVar]) {
 538         isOutside=kTRUE;
 539         break;
 540       }
 541     }
 542     if (isOutside) continue;
 543
 544     Int_t neighbours = 0; // number of neighbours to average with
 545
 546     for (Int_t iVar=0; iVar<fNVariables; iVar++) {
 547       if (fCoordinatesN_T[iVar] > 1) { // must not be on low edge border
 548         fCoordinatesN_T[iVar]-- ; //get lower neighbouring bin
 549         content += copy->GetBinContent(fCoordinatesN_T);
 550         error2  += TMath::Power(copy->GetBinError(fCoordinatesN_T),2);
 551         neighbours++;
 552         fCoordinatesN_T[iVar]++ ; //back to initial coordinate
 553       }
 554       if (fCoordinatesN_T[iVar] < numBins[iVar]) { // must not be on up edge border
 555         fCoordinatesN_T[iVar]++ ; //get upper neighbouring bin
 556         content += copy->GetBinContent(fCoordinatesN_T);
 557         error2  += TMath::Power(copy->GetBinError(fCoordinatesN_T),2);
 558         neighbours++;
 559         fCoordinatesN_T[iVar]-- ; //back to initial coordinate
 560       }
 561     }
 562     // make an average
 563     fUnfolded->SetBinContent(fCoordinatesN_T,content/(1.+neighbours));
 564     fUnfolded->SetBinError  (fCoordinatesN_T,TMath::Sqrt(error2)/(1.+neighbours));
 565   }
 566   delete [] numBins;
 567   delete copy;
 568   return 0;
 569 }
 570
 571 //______________________________________________________________
 572
 573 Short_t AliCFUnfolding::SmoothUsingFunction() {
 574   //
 575   // Fits the unfolded spectrum using the function fSmoothFunction
 576   //
 577
 578   AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
 579
 580   for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
 581
 582   Int_t fitResult = 0;
 583
 584   switch (fNVariables) {
 585   case 1 : fitResult = fUnfolded->Projection(0)    ->Fit(fSmoothFunction,fSmoothOption); break;
 586   case 2 : fitResult = fUnfolded->Projection(1,0)  ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
 587   case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
 588   default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
 589   }
 590
 591   if (fitResult != 0) {
 592     AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
 593     return 1;
 594   }
 595
 596   Int_t nDim = fNVariables;
 597   Int_t* bins = new Int_t[nDim]; // number of bins for each variable
 598   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
 599
 600   for (Int_t iVar=0; iVar<nDim; iVar++) {
 601     bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
 602     nBins *= bins[iVar];
 603   }
 604
 605   Int_t *bin  = new Int_t[nDim];    // bin to fill the THnSparse (holding the bin coordinates)
 606   Double_t x[3] = {0,0,0} ;         // value in bin center (max dimension is 3 (TF3))
 607
 608   // loop on the bins and update of fUnfolded
 609   // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
 610   for (Long_t iBin=0; iBin<nBins; iBin++) {
 611     Long_t bin_tmp = iBin ;
 612     for (Int_t iVar=0; iVar<nDim; iVar++) {
 613       bin[iVar] = 1 + bin_tmp % bins[iVar] ;
 614       bin_tmp /= bins[iVar] ;
 615       x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
 616     }
 617     Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
 618     fUnfolded->SetBinContent(bin,functionValue);
 619     fUnfolded->SetBinError  (bin,functionValue*fUnfolded->GetBinError(bin));
 620   }
 621   return 0;
 622 }
 623
 624 //______________________________________________________________
 625
 626 void AliCFUnfolding::CreateFlatPrior() {
 627   //
 628   // Creates a flat prior distribution
 629   //
 630
 631   AliInfo("Creating a flat a priori distribution");
 632
 633   // create the frame of the THnSparse given (for example) the one from the efficiency map
 634   fPrior = (THnSparse*) fEfficiency->Clone();
 635
 636   if (fNVariables != fPrior->GetNdimensions())
 637     AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
 638
 639   Int_t nDim = fNVariables;
 640   Int_t* bins = new Int_t[nDim]; // number of bins for each variable
 641   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
 642
 643   for (Int_t iVar=0; iVar<nDim; iVar++) {
 644     bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
 645     nBins *= bins[iVar];
 646   }
 647
 648   Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
 649
 650   // loop that sets 1 in each bin
 651   for (Long_t iBin=0; iBin<nBins; iBin++) {
 652     Long_t bin_tmp = iBin ;
 653     for (Int_t iVar=0; iVar<nDim; iVar++) {
 654       bin[iVar] = 1 + bin_tmp % bins[iVar] ;
 655       bin_tmp /= bins[iVar] ;
 656     }
 657     fPrior->SetBinContent(bin,1.); // put 1 everywhere
 658     fPrior->SetBinError  (bin,0.); // put 0 everywhere
 659   }
 660
 661   fOriginalPrior = (THnSparse*)fPrior->Clone();
 662
 663   delete [] bin;
 664   delete [] bins;
 665 }