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