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