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