ROOT logo
/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//---------------------------------------------------------------------//
//                                                                     //
// AliCFUnfolding Class                                                //
// Class to handle general unfolding procedure                         // 
// For the moment only bayesian unfolding is supported                 //
// The next steps are to add chi2 minimisation and weighting methods   //
//                                                                     //
//                                                                     //
//                                                                     //
// Use :                                                               //
// -------                                                             //
// The Bayesian unfolding consists of several iterations.              //
// At each iteration, an inverse response matrix is calculated, given  //
// the measured spectrum, the a priori (guessed) spectrum,             //
// the efficiency spectrum and the response matrix.                    //
//                                                                     //
// Then at each iteration, the unfolded spectrum is calculated using   //
// the inverse response : the goal is to get an unfolded spectrum      //
// similar (according to some criterion) to the a priori one.          //
// If the difference is too big, another iteration is performed :      //
// the a priori spectrum is updated to the unfolded one from the       //
// previous iteration, and so on so forth, until the maximum number    //
// of iterations or the similarity criterion is reached.               //
//                                                                     //
// Chi2 calculation became absolute with the correlated error          //
// calculation.                                                        //
// Errors on the unfolded distribution are not known until the end     //
// Use the convergence criterion instead                               //
//                                                                     //
// Currently the user has to define the max. number of iterations      //
// (::SetMaxNumberOfIterations)                                        //
// and                                                                 //
//     - the chi2 below which the procedure will stop                  //
// (::SetMaxChi2 or ::SetMaxChi2PerDOF)   (OBSOLETE)                   //
//     - the convergence criterion below which the procedure will stop //
// SetMaxConvergencePerDOF(Double_t val);                              //
//                                                                     //
// Correlated error calculation can be activated by using:             //
// SetUseCorrelatedErrors(Bool_t b) in combination with convergence    //
// criterion                                                           //
// Documentation about correlated error calculation method can be      //
// found in AliCFUnfolding::CalculateCorrelatedErrors()                //
// Author: marta.verweij@cern.ch                                       //
//                                                                     //
// An optional possibility is to smooth the unfolded spectrum at the   //
// end of each iteration, either using a fit function                  //
// (only if #dimensions <=3)                                           //
// or a simple averaging using the neighbouring bins values.           //
// This is possible calling the function ::UseSmoothing                //
// If no argument is passed to this function, then the second option   //
// is used.                                                            //
//                                                                     //
// IMPORTANT:                                                          //
//-----------                                                          //
// With this approach, the efficiency map must be calculated           //
// with *simulated* values only, otherwise the method won't work.      //
//                                                                     //
// ex: efficiency(bin_pt) = number_rec(bin_pt) / number_sim(bin_pt)    //
//                                                                     //
// the pt bin "bin_pt" must always be the same in both the efficiency  //
// numerator and denominator.                                          //
// This is why the efficiency map has to be created by a method        //
// from which both reconstructed and simulated values are accessible   //
// simultaneously.                                                     //
//                                                                     //
//                                                                     //
//---------------------------------------------------------------------//
// Author : renaud.vernet@cern.ch                                      //
//---------------------------------------------------------------------//


#include "AliCFUnfolding.h"
#include "TMath.h"
#include "TAxis.h"
#include "TF1.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "TRandom3.h"


ClassImp(AliCFUnfolding)

//______________________________________________________________

AliCFUnfolding::AliCFUnfolding() :
  TNamed(),
  fResponseOrig(0x0),
  fPriorOrig(0x0),
  fEfficiencyOrig(0x0),
  fMeasuredOrig(0x0),
  fMaxNumIterations(0),
  fNVariables(0),
  fUseSmoothing(kFALSE),
  fSmoothFunction(0x0),
  fSmoothOption("iremn"),
  fMaxConvergence(0),
  fNRandomIterations(0),
  fResponse(0x0),
  fPrior(0x0),
  fEfficiency(0x0),
  fMeasured(0x0),
  fInverseResponse(0x0),
  fMeasuredEstimate(0x0),
  fConditional(0x0),
  fUnfolded(0x0),
  fUnfoldedFinal(0x0),
  fCoordinates2N(0x0),
  fCoordinatesN_M(0x0),
  fCoordinatesN_T(0x0),
  fRandomResponse(0x0),
  fRandomEfficiency(0x0),
  fRandomMeasured(0x0),
  fRandom3(0x0),
  fDeltaUnfoldedP(0x0),
  fDeltaUnfoldedN(0x0),
  fNCalcCorrErrors(0),
  fRandomSeed(0)
{
  //
  // default constructor
  //
}

//______________________________________________________________

AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar, 
			       const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior ,
			       Double_t maxConvergencePerDOF, UInt_t randomSeed, Int_t maxNumIterations
			       ) :
  TNamed(name,title),
  fResponseOrig((THnSparse*)response->Clone()),
  fPriorOrig(0x0),
  fEfficiencyOrig((THnSparse*)efficiency->Clone()),
  fMeasuredOrig((THnSparse*)measured->Clone()),
  fMaxNumIterations(maxNumIterations),
  fNVariables(nVar),
  fUseSmoothing(kFALSE),
  fSmoothFunction(0x0),
  fSmoothOption("iremn"),
  fMaxConvergence(0),
  fNRandomIterations(maxNumIterations),
  fResponse((THnSparse*)response->Clone()),
  fPrior(0x0),
  fEfficiency((THnSparse*)efficiency->Clone()),
  fMeasured((THnSparse*)measured->Clone()),
  fInverseResponse(0x0),
  fMeasuredEstimate(0x0),
  fConditional(0x0),
  fUnfolded(0x0),
  fUnfoldedFinal(0x0),
  fCoordinates2N(0x0),
  fCoordinatesN_M(0x0),
  fCoordinatesN_T(0x0),
  fRandomResponse((THnSparse*)response->Clone()),
  fRandomEfficiency((THnSparse*)efficiency->Clone()),
  fRandomMeasured((THnSparse*)measured->Clone()),
  fRandom3(0x0),
  fDeltaUnfoldedP(0x0),
  fDeltaUnfoldedN(0x0),
  fNCalcCorrErrors(0),
  fRandomSeed(randomSeed)
{
  //
  // named constructor
  //

  AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));

  if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
  else {
    fPrior = (THnSparse*) prior->Clone();
    fPriorOrig = (THnSparse*)fPrior->Clone();
    if (fPrior->GetNdimensions() != fNVariables) 
      AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
  }

  if (fEfficiency->GetNdimensions() != fNVariables) 
    AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions()));
  if (fMeasured->GetNdimensions() != fNVariables) 
    AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions()));
  if (fResponse->GetNdimensions() != 2*fNVariables) 
    AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions()));
  

  for (Int_t iVar=0; iVar<fNVariables; iVar++) {
    AliInfo(Form("prior      matrix has %d bins in dimension %d",fPrior     ->GetAxis(iVar)->GetNbins(),iVar));
    AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
    AliInfo(Form("measured   matrix has %d bins in dimension %d",fMeasured  ->GetAxis(iVar)->GetNbins(),iVar));
  }

  fRandomResponse  ->SetTitle("Randomized response matrix");
  fRandomEfficiency->SetTitle("Randomized efficiency");
  fRandomMeasured  ->SetTitle("Randomized measured");
  SetMaxConvergencePerDOF(maxConvergencePerDOF)  ;
  Init();
}

//______________________________________________________________

AliCFUnfolding::~AliCFUnfolding() {
  //
  // destructor
  //
  if (fResponse)           delete fResponse;
  if (fPrior)              delete fPrior;
  if (fEfficiency)         delete fEfficiency;
  if (fEfficiencyOrig)     delete fEfficiencyOrig;
  if (fMeasured)           delete fMeasured;
  if (fMeasuredOrig)       delete fMeasuredOrig;
  if (fSmoothFunction)     delete fSmoothFunction;
  if (fPriorOrig)          delete fPriorOrig;
  if (fInverseResponse)    delete fInverseResponse;
  if (fMeasuredEstimate)   delete fMeasuredEstimate;
  if (fConditional)        delete fConditional;
  if (fUnfolded)           delete fUnfolded;
  if (fUnfoldedFinal)      delete fUnfoldedFinal;
  if (fCoordinates2N)      delete [] fCoordinates2N; 
  if (fCoordinatesN_M)     delete [] fCoordinatesN_M; 
  if (fCoordinatesN_T)     delete [] fCoordinatesN_T; 
  if (fRandomResponse)     delete fRandomResponse;
  if (fRandomEfficiency)   delete fRandomEfficiency;
  if (fRandomMeasured)     delete fRandomMeasured;
  if (fRandom3)            delete fRandom3;
  if (fDeltaUnfoldedP)     delete fDeltaUnfoldedP;
  if (fDeltaUnfoldedN)     delete fDeltaUnfoldedN;
 
}

//______________________________________________________________

void AliCFUnfolding::Init() {
  //
  // initialisation function : creates internal settings
  //

  fRandom3 = new TRandom3(fRandomSeed);

  fCoordinates2N  = new Int_t[2*fNVariables];
  fCoordinatesN_M = new Int_t[fNVariables];
  fCoordinatesN_T = new Int_t[fNVariables];

  // create the matrix of conditional probabilities P(M|T)
  CreateConditional(); //done only once at initialization
  
  // create the frame of the inverse response matrix
  fInverseResponse  = (THnSparse*) fResponse->Clone();
  // create the frame of the unfolded spectrum
  fUnfolded = (THnSparse*) fPrior->Clone();
  fUnfolded->SetTitle("Unfolded");
  // create the frame of the measurement estimate spectrum
  fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
  
  // create the frame of the delta profiles
  fDeltaUnfoldedP = (THnSparse*)fPrior->Clone();
  fDeltaUnfoldedP->SetTitle("#Delta unfolded");
  fDeltaUnfoldedP->Reset();
  fDeltaUnfoldedN = (THnSparse*)fPrior->Clone();
  fDeltaUnfoldedN->SetTitle("");
  fDeltaUnfoldedN->Reset();


}


//______________________________________________________________

void AliCFUnfolding::CreateEstMeasured() {
  //
  // This function creates a estimate (M) of the reconstructed spectrum 
  // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND)
  //
  // --> P(M) = SUM   { P(M|T)    * P(T) }
  // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)}
  //
  // This is needed to calculate the inverse response matrix
  //


  // clean the measured estimate spectrum
  fMeasuredEstimate->Reset();

  THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
  priorTimesEff->Multiply(fEfficiency);

  // fill it
  for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
    Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
    GetCoordinates();
    Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
    Double_t fill = conditionalValue * priorTimesEffValue ;
    
    if (fill>0.) {
      fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
      fMeasuredEstimate->SetBinError(fCoordinatesN_M,0.);
    }
  }
  delete priorTimesEff ;
}

//______________________________________________________________

void AliCFUnfolding::CreateInvResponse() {
  //
  // Creates the inverse response matrix (INV) with Bayesian method
  //  : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
  //
  // --> P(T|M)   = P(M|T)    * P(T) * eff(T) / SUM   { P(M|T)    * P(T) }
  // --> INV(i,j) = COND(i,j) * T(j) * E(j)   / SUM_k { COND(i,k) * T(k) }
  //

  THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
  priorTimesEff->Multiply(fEfficiency);

  for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
    Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
    GetCoordinates();
    Double_t estMeasuredValue   = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
    Double_t priorTimesEffValue = priorTimesEff    ->GetBinContent(fCoordinatesN_T);
    Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
    if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
      fInverseResponse->SetBinContent(fCoordinates2N,fill);
      fInverseResponse->SetBinError  (fCoordinates2N,0.);
    }
  } 
  delete priorTimesEff ;
}

//______________________________________________________________

void AliCFUnfolding::Unfold() {
  //
  // Main routine called by the user : 
  // it calculates the unfolded spectrum from the response matrix, measured spectrum and efficiency
  // several iterations are performed until a reasonable chi2 or convergence criterion is reached
  //

  Int_t iIterBayes     = 0 ;
  Double_t convergence = 0.;

  for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations

    CreateEstMeasured(); // create measured estimate from prior
    CreateInvResponse(); // create inverse response  from prior
    CreateUnfolded();    // create unfoled spectrum  from measured and inverse response

    convergence = GetConvergence();
    AliDebug(0,Form("convergence at iteration %d is %e",iIterBayes,convergence));

    if (fMaxConvergence>0. && convergence<fMaxConvergence && fNCalcCorrErrors == 0) {
      fNRandomIterations = iIterBayes;
      AliDebug(0,Form("convergence is met at iteration %d",iIterBayes));
      break;
    }

    if (fUseSmoothing) {
      if (Smooth()) {
	AliError("Couldn't smooth the unfolded spectrum!!");
	if (fNCalcCorrErrors>0) {
	  AliInfo(Form("=======================\nUnfold of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
	}
	else {
	  AliInfo(Form("\n\n=======================\nFinish at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
	}
	return;
      }
    }

    // update the prior distribution
    if (fPrior) delete fPrior ;
    fPrior = (THnSparse*)fUnfolded->Clone() ;
    fPrior->SetTitle("Prior");

  } // end bayes iteration

  if (fNCalcCorrErrors==0) fUnfoldedFinal = (THnSparse*) fUnfolded->Clone() ;

  //
  //for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) AliDebug(2,Form("%e\n",fUnfoldedFinal->GetBinError(iBin)));
  //

  if (fNCalcCorrErrors == 0) {
    AliInfo("\n================================================\nFinished bayes iteration, now calculating errors...\n================================================\n");
    fNCalcCorrErrors = 1;
    CalculateCorrelatedErrors();
  }

  if (fNCalcCorrErrors >1 ) {
    AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
  }
  else if(fNCalcCorrErrors>0) {
    AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
  }
}

//______________________________________________________________

void AliCFUnfolding::CreateUnfolded() {
  //
  // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
  // We have P(T) = SUM   { P(T|M)   * P(M) } 
  //   -->   T(i) = SUM_k { INV(i,k) * M(k) }
  //


  // clear the unfolded spectrum
  // if in the process of error calculation, the random unfolded spectrum is created
  // otherwise the normal unfolded spectrum is created

  fUnfolded->Reset();
  
  for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
    Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
    GetCoordinates();
    Double_t effValue      = fEfficiency->GetBinContent(fCoordinatesN_T);
    Double_t measuredValue = fMeasured  ->GetBinContent(fCoordinatesN_M);
    Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;

    if (fill>0.) {
      // set errors to zero
      // true errors will be filled afterwards
      Double_t err = 0.;
      fUnfolded->SetBinError  (fCoordinatesN_T,err);
      fUnfolded->AddBinContent(fCoordinatesN_T,fill);
    }
  }
}

//______________________________________________________________

void AliCFUnfolding::CalculateCorrelatedErrors() {

  // Step 1: Create randomized distribution (fRandomXXXX) of each bin of 
  //         the measured spectrum to calculate correlated errors. 
  //         Poisson statistics: mean = measured value of bin
  // Step 2: Unfold randomized distribution
  // Step 3: Store difference of unfolded spectrum from measured distribution and 
  //         unfolded distribution from randomized distribution 
  //         -> fDeltaUnfoldedP (TProfile with option "S")
  // Step 4: Repeat Step 1-3 several times (fNRandomIterations)
  // Step 5: The spread of fDeltaUnfoldedP for each bin is the error on the unfolded spectrum of that specific bin


  //Do fNRandomIterations = bayes iterations performed
  for (int i=0; i<fNRandomIterations; i++) {
    
    // reset prior to original one
    if (fPrior) delete fPrior ;
    fPrior = (THnSparse*) fPriorOrig->Clone();

    // create randomized distribution and stick measured spectrum to it
    CreateRandomizedDist();

    if (fResponse) delete fResponse ;
    fResponse = (THnSparse*) fRandomResponse->Clone();
    fResponse->SetTitle("Response");

    if (fEfficiency) delete fEfficiency ;
    fEfficiency = (THnSparse*) fRandomEfficiency->Clone();
    fEfficiency->SetTitle("Efficiency");

    if (fMeasured)   delete fMeasured   ;
    fMeasured = (THnSparse*) fRandomMeasured->Clone();
    fMeasured->SetTitle("Measured");

    //unfold with randomized distributions
    Unfold();
    FillDeltaUnfoldedProfile();
  }

  // Get statistical errors for final unfolded spectrum
  // ie. spread of each pt bin in fDeltaUnfoldedP
  Double_t meanx2 = 0.;
  Double_t mean = 0.;
  Double_t checksigma = 0.;
  Double_t entriesInBin = 0.;
  for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) {
    fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M);
    mean = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M);
    meanx2 = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M);
    entriesInBin = fDeltaUnfoldedN->GetBinContent(fCoordinatesN_M);
    if(entriesInBin > 1.) checksigma = TMath::Sqrt((entriesInBin/(entriesInBin-1.))*TMath::Abs(meanx2-mean*mean));
    //printf("mean %f, meanx2 %f, sigmacheck %f, nentries %f\n",mean, meanx2, checksigma,entriesInBin);
    //AliDebug(2,Form("filling error %e\n",sigma));
    fUnfoldedFinal->SetBinError(fCoordinatesN_M,checksigma);
  }

  // now errors are calculated
  fNCalcCorrErrors = 2;
}

//______________________________________________________________
void AliCFUnfolding::CreateRandomizedDist() {
  //
  // Create randomized dist from original measured distribution
  // This distribution is created several times, each time with a different random number
  //

  for (Long_t iBin=0; iBin<fResponseOrig->GetNbins(); iBin++) {
    Double_t val = fResponseOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
    Double_t err = fResponseOrig->GetBinError(fCoordinatesN_M);        //used as sigma
    Double_t ran = fRandom3->Gaus(val,err);
    // random        = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
    fRandomResponse->SetBinContent(iBin,ran);
  }
  for (Long_t iBin=0; iBin<fEfficiencyOrig->GetNbins(); iBin++) {
    Double_t val = fEfficiencyOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
    Double_t err = fEfficiencyOrig->GetBinError(fCoordinatesN_M);        //used as sigma
    Double_t ran = fRandom3->Gaus(val,err);
    // random        = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
    fRandomEfficiency->SetBinContent(iBin,ran);
  }
  for (Long_t iBin=0; iBin<fMeasuredOrig->GetNbins(); iBin++) {
    Double_t val = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
    Double_t err = fMeasuredOrig->GetBinError(fCoordinatesN_M);        //used as sigma
    Double_t ran = fRandom3->Gaus(val,err);
    // random        = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
    fRandomMeasured->SetBinContent(iBin,ran);
  }
}

//______________________________________________________________
void AliCFUnfolding::FillDeltaUnfoldedProfile() {
  //
  // Store difference of unfolded spectrum from measured distribution and unfolded spectrum from randomized distribution
  // The delta profile has been set to a THnSparse to handle N dimension
  // The THnSparse contains in each bin the mean value and spread of the difference 
  // This function updates the profile wrt to its previous mean and error
  // The relation between iterations (n+1) and n is as follows :
  //  mean_{n+1} = (n*mean_n + value_{n+1}) / (n+1)
  // sigma_{n+1} = sqrt { 1/(n+1) * [ n*sigma_n^2 + (n^2+n)*(mean_{n+1}-mean_n)^2 ] }    (can this be optimized?)

  for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) {
    Double_t deltaInBin   = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M) - fUnfolded->GetBinContent(fCoordinatesN_M);
    Double_t entriesInBin = fDeltaUnfoldedN->GetBinContent(fCoordinatesN_M);
    //AliDebug(2,Form("%e %e ==> delta = %e\n",fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M),fUnfolded->GetBinContent(iBin),deltaInBin));

    //printf("deltaInBin %f\n",deltaInBin);
    //printf("pt %f\n",ptaxis->GetBinCenter(iBin+1));
    
    Double_t mean_n = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M) ;
    Double_t mean_nplus1 = mean_n ;
    mean_nplus1 *= entriesInBin ;
    mean_nplus1 += deltaInBin ;
    mean_nplus1 /= (entriesInBin+1) ;

    Double_t meanx2_n = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M) ;
    Double_t meanx2_nplus1 = meanx2_n ;
    meanx2_nplus1 *= entriesInBin ;
    meanx2_nplus1 += (deltaInBin*deltaInBin) ;
    meanx2_nplus1 /= (entriesInBin+1) ;

    //AliDebug(2,Form("sigma = %e\n",sigma));

    fDeltaUnfoldedP->SetBinError(fCoordinatesN_M,meanx2_nplus1) ;
    fDeltaUnfoldedP->SetBinContent(fCoordinatesN_M,mean_nplus1) ;
    fDeltaUnfoldedN->SetBinContent(fCoordinatesN_M,entriesInBin+1);
  }
}

//______________________________________________________________

void AliCFUnfolding::GetCoordinates() {
  //
  // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
  //
  for (Int_t i = 0; i<fNVariables ; i++) {
    fCoordinatesN_M[i] = fCoordinates2N[i];
    fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
  }
}

//______________________________________________________________

void AliCFUnfolding::CreateConditional() {
  //
  // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
  //
  //  --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
  //

  fConditional = (THnSparse*) fResponse->Clone();  // output of this function

  Int_t* dim = new Int_t [fNVariables];
  for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)

  THnSparse* responseInT = fConditional->Projection(fNVariables,dim,"E");     // output denominator : 
                                                                              // projection of the response matrix on the TRUE axis
  delete [] dim; 

  // fill the conditional probability matrix
  for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
    Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
    GetCoordinates();
    Double_t projValue = responseInT->GetBinContent(fCoordinatesN_T);
   
    Double_t fill = responseValue / projValue ;
    if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
      fConditional->SetBinContent(fCoordinates2N,fill);
      Double_t err = 0.;
      fConditional->SetBinError  (fCoordinates2N,err);
    }
  }
  delete responseInT ;
}
//______________________________________________________________

Int_t AliCFUnfolding::GetDOF() {
  //
  // number of dof = number of bins
  //

  Int_t nDOF = 1 ;
  for (Int_t iDim=0; iDim<fNVariables; iDim++) {
    nDOF *= fPrior->GetAxis(iDim)->GetNbins();
  }
  AliDebug(0,Form("Number of degrees of freedom = %d",nDOF));
  return nDOF;
}

//______________________________________________________________

Double_t AliCFUnfolding::GetChi2() {
  //
  // Returns the chi2 between unfolded and a priori spectrum
  // This function became absolute with the correlated error calculation.
  // Errors on the unfolded distribution are not known until the end
  // Use the convergence criterion instead
  //

  Double_t chi2      = 0. ;
  Double_t error_unf = 0.;
  for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
    Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
    error_unf  = fUnfolded->GetBinError(fCoordinatesN_T);
    chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ;
  }
  return chi2;
}

//______________________________________________________________

Double_t AliCFUnfolding::GetConvergence() {
  //
  // Returns convergence criterion = \sum_t ((U_t^{n-1}-U_t^n)/U_t^{n-1})^2
  // U is unfolded spectrum, t is the bin, n = current, n-1 = previous
  //
  Double_t convergence = 0.;
  Double_t priorValue  = 0.;
  Double_t currentValue = 0.;
  for (Long_t iBin=0; iBin < fPrior->GetNbins(); iBin++) {
    priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
    currentValue = fUnfolded->GetBinContent(fCoordinatesN_T);

    if (priorValue > 0.)
      convergence += ((priorValue-currentValue)/priorValue)*((priorValue-currentValue)/priorValue);
    else 
      AliWarning(Form("priorValue = %f. Adding 0 to convergence criterion.",priorValue)); 
  }
  return convergence;
}

//______________________________________________________________

void AliCFUnfolding::SetMaxConvergencePerDOF(Double_t val) {
  //
  // Max. convergence criterion per degree of freedom : user setting
  // convergence criterion = DOF*val; DOF = number of bins
  // In Jan-Fiete's multiplicity note: Convergence criterion = DOF*0.001^2
  //

  Int_t nDOF = GetDOF() ;
  fMaxConvergence = val * nDOF ;
  AliInfo(Form("MaxConvergence = %e. Number of degrees of freedom = %d",fMaxConvergence,nDOF));
}

//______________________________________________________________

Short_t AliCFUnfolding::Smooth() {
  //
  // Smoothes the unfolded spectrum
  //
  // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
  // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn 
  //
  
  if (fSmoothFunction) {
    AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
    return SmoothUsingFunction();
  }
  else return SmoothUsingNeighbours(fUnfolded);
}

//______________________________________________________________

Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) {
  //
  // Smoothes the unfolded spectrum using neighouring bins
  //

  Int_t  const nDimensions = hist->GetNdimensions() ;
  Int_t* coordinates = new Int_t[nDimensions];

  Int_t* numBins = new Int_t[nDimensions];
  for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins();
  
  //need a copy because hist will be updated during the loop, and this creates problems
  THnSparse* copy = (THnSparse*)hist->Clone();

  for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
    Double_t content = copy->GetBinContent(iBin,coordinates);
    Double_t error2  = TMath::Power(copy->GetBinError(iBin),2);

    // skip the under/overflow bins...
    Bool_t isOutside = kFALSE ;
    for (Int_t iVar=0; iVar<nDimensions; iVar++) {
      if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) {
	isOutside=kTRUE;
	break;
      }
    }
    if (isOutside) continue;
    
    Int_t neighbours = 0; // number of neighbours to average with

    for (Int_t iVar=0; iVar<nDimensions; iVar++) {
      if (coordinates[iVar] > 1) { // must not be on low edge border
	coordinates[iVar]-- ; //get lower neighbouring bin 
	content += copy->GetBinContent(coordinates);
	error2  += TMath::Power(copy->GetBinError(coordinates),2);
	neighbours++;
	coordinates[iVar]++ ; //back to initial coordinate
      }
      if (coordinates[iVar] < numBins[iVar]) { // must not be on up edge border
	coordinates[iVar]++ ; //get upper neighbouring bin
	content += copy->GetBinContent(coordinates);
	error2  += TMath::Power(copy->GetBinError(coordinates),2);
	neighbours++;
	coordinates[iVar]-- ; //back to initial coordinate
      }
    }
    // make an average
    hist->SetBinContent(coordinates,content/(1.+neighbours));
    hist->SetBinError  (coordinates,TMath::Sqrt(error2)/(1.+neighbours));
  }
  delete [] numBins;
  delete [] coordinates ;
  delete copy;
  return 0;
}

//______________________________________________________________

Short_t AliCFUnfolding::SmoothUsingFunction() {
  //
  // Fits the unfolded spectrum using the function fSmoothFunction
  //

  AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));

  for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));

  Int_t fitResult = 0;

  switch (fNVariables) {
  case 1 : fitResult = fUnfolded->Projection(0)    ->Fit(fSmoothFunction,fSmoothOption); break;
  case 2 : fitResult = fUnfolded->Projection(1,0)  ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
  case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
  default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
  }

  if (fitResult != 0) {
    AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
    return 1;
  }

  Int_t nDim = fNVariables;
  Int_t* bins = new Int_t[nDim]; // number of bins for each variable
  Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse

  for (Int_t iVar=0; iVar<nDim; iVar++) {
    bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
    nBins *= bins[iVar];
  }

  Int_t *bin  = new Int_t[nDim];    // bin to fill the THnSparse (holding the bin coordinates)
  Double_t x[3] = {0,0,0} ;         // value in bin center (max dimension is 3 (TF3))

  // loop on the bins and update of fUnfolded
  // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
  for (Long_t iBin=0; iBin<nBins; iBin++) {
    Long_t bin_tmp = iBin ;
    for (Int_t iVar=0; iVar<nDim; iVar++) {
      bin[iVar] = 1 + bin_tmp % bins[iVar] ;
      bin_tmp /= bins[iVar] ;
      x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
    }
    Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
    fUnfolded->SetBinError  (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin));
    fUnfolded->SetBinContent(bin,functionValue);
  }
  delete [] bins;
  delete [] bin ;
  return 0;
}

//______________________________________________________________

void AliCFUnfolding::CreateFlatPrior() {
  //
  // Creates a flat prior distribution
  // 

  AliInfo("Creating a flat a priori distribution");
  
  // create the frame of the THnSparse given (for example) the one from the efficiency map
  fPrior = (THnSparse*) fEfficiency->Clone();
  fPrior->SetTitle("Prior");

  if (fNVariables != fPrior->GetNdimensions()) 
    AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));

  Int_t nDim = fNVariables;
  Int_t* bins = new Int_t[nDim]; // number of bins for each variable
  Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse

  for (Int_t iVar=0; iVar<nDim; iVar++) {
    bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
    nBins *= bins[iVar];
  }

  Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)

  // loop that sets 1 in each bin
  for (Long_t iBin=0; iBin<nBins; iBin++) {
    Long_t bin_tmp = iBin ;
    for (Int_t iVar=0; iVar<nDim; iVar++) {
      bin[iVar] = 1 + bin_tmp % bins[iVar] ;
      bin_tmp /= bins[iVar] ;
    }
    fPrior->SetBinContent(bin,1.); // put 1 everywhere
    fPrior->SetBinError  (bin,0.); // put 0 everywhere
  }
  
  fPriorOrig = (THnSparse*)fPrior->Clone();

  delete [] bin;
  delete [] bins;
}
 AliCFUnfolding.cxx:1
 AliCFUnfolding.cxx:2
 AliCFUnfolding.cxx:3
 AliCFUnfolding.cxx:4
 AliCFUnfolding.cxx:5
 AliCFUnfolding.cxx:6
 AliCFUnfolding.cxx:7
 AliCFUnfolding.cxx:8
 AliCFUnfolding.cxx:9
 AliCFUnfolding.cxx:10
 AliCFUnfolding.cxx:11
 AliCFUnfolding.cxx:12
 AliCFUnfolding.cxx:13
 AliCFUnfolding.cxx:14
 AliCFUnfolding.cxx:15
 AliCFUnfolding.cxx:16
 AliCFUnfolding.cxx:17
 AliCFUnfolding.cxx:18
 AliCFUnfolding.cxx:19
 AliCFUnfolding.cxx:20
 AliCFUnfolding.cxx:21
 AliCFUnfolding.cxx:22
 AliCFUnfolding.cxx:23
 AliCFUnfolding.cxx:24
 AliCFUnfolding.cxx:25
 AliCFUnfolding.cxx:26
 AliCFUnfolding.cxx:27
 AliCFUnfolding.cxx:28
 AliCFUnfolding.cxx:29
 AliCFUnfolding.cxx:30
 AliCFUnfolding.cxx:31
 AliCFUnfolding.cxx:32
 AliCFUnfolding.cxx:33
 AliCFUnfolding.cxx:34
 AliCFUnfolding.cxx:35
 AliCFUnfolding.cxx:36
 AliCFUnfolding.cxx:37
 AliCFUnfolding.cxx:38
 AliCFUnfolding.cxx:39
 AliCFUnfolding.cxx:40
 AliCFUnfolding.cxx:41
 AliCFUnfolding.cxx:42
 AliCFUnfolding.cxx:43
 AliCFUnfolding.cxx:44
 AliCFUnfolding.cxx:45
 AliCFUnfolding.cxx:46
 AliCFUnfolding.cxx:47
 AliCFUnfolding.cxx:48
 AliCFUnfolding.cxx:49
 AliCFUnfolding.cxx:50
 AliCFUnfolding.cxx:51
 AliCFUnfolding.cxx:52
 AliCFUnfolding.cxx:53
 AliCFUnfolding.cxx:54
 AliCFUnfolding.cxx:55
 AliCFUnfolding.cxx:56
 AliCFUnfolding.cxx:57
 AliCFUnfolding.cxx:58
 AliCFUnfolding.cxx:59
 AliCFUnfolding.cxx:60
 AliCFUnfolding.cxx:61
 AliCFUnfolding.cxx:62
 AliCFUnfolding.cxx:63
 AliCFUnfolding.cxx:64
 AliCFUnfolding.cxx:65
 AliCFUnfolding.cxx:66
 AliCFUnfolding.cxx:67
 AliCFUnfolding.cxx:68
 AliCFUnfolding.cxx:69
 AliCFUnfolding.cxx:70
 AliCFUnfolding.cxx:71
 AliCFUnfolding.cxx:72
 AliCFUnfolding.cxx:73
 AliCFUnfolding.cxx:74
 AliCFUnfolding.cxx:75
 AliCFUnfolding.cxx:76
 AliCFUnfolding.cxx:77
 AliCFUnfolding.cxx:78
 AliCFUnfolding.cxx:79
 AliCFUnfolding.cxx:80
 AliCFUnfolding.cxx:81
 AliCFUnfolding.cxx:82
 AliCFUnfolding.cxx:83
 AliCFUnfolding.cxx:84
 AliCFUnfolding.cxx:85
 AliCFUnfolding.cxx:86
 AliCFUnfolding.cxx:87
 AliCFUnfolding.cxx:88
 AliCFUnfolding.cxx:89
 AliCFUnfolding.cxx:90
 AliCFUnfolding.cxx:91
 AliCFUnfolding.cxx:92
 AliCFUnfolding.cxx:93
 AliCFUnfolding.cxx:94
 AliCFUnfolding.cxx:95
 AliCFUnfolding.cxx:96
 AliCFUnfolding.cxx:97
 AliCFUnfolding.cxx:98
 AliCFUnfolding.cxx:99
 AliCFUnfolding.cxx:100
 AliCFUnfolding.cxx:101
 AliCFUnfolding.cxx:102
 AliCFUnfolding.cxx:103
 AliCFUnfolding.cxx:104
 AliCFUnfolding.cxx:105
 AliCFUnfolding.cxx:106
 AliCFUnfolding.cxx:107
 AliCFUnfolding.cxx:108
 AliCFUnfolding.cxx:109
 AliCFUnfolding.cxx:110
 AliCFUnfolding.cxx:111
 AliCFUnfolding.cxx:112
 AliCFUnfolding.cxx:113
 AliCFUnfolding.cxx:114
 AliCFUnfolding.cxx:115
 AliCFUnfolding.cxx:116
 AliCFUnfolding.cxx:117
 AliCFUnfolding.cxx:118
 AliCFUnfolding.cxx:119
 AliCFUnfolding.cxx:120
 AliCFUnfolding.cxx:121
 AliCFUnfolding.cxx:122
 AliCFUnfolding.cxx:123
 AliCFUnfolding.cxx:124
 AliCFUnfolding.cxx:125
 AliCFUnfolding.cxx:126
 AliCFUnfolding.cxx:127
 AliCFUnfolding.cxx:128
 AliCFUnfolding.cxx:129
 AliCFUnfolding.cxx:130
 AliCFUnfolding.cxx:131
 AliCFUnfolding.cxx:132
 AliCFUnfolding.cxx:133
 AliCFUnfolding.cxx:134
 AliCFUnfolding.cxx:135
 AliCFUnfolding.cxx:136
 AliCFUnfolding.cxx:137
 AliCFUnfolding.cxx:138
 AliCFUnfolding.cxx:139
 AliCFUnfolding.cxx:140
 AliCFUnfolding.cxx:141
 AliCFUnfolding.cxx:142
 AliCFUnfolding.cxx:143
 AliCFUnfolding.cxx:144
 AliCFUnfolding.cxx:145
 AliCFUnfolding.cxx:146
 AliCFUnfolding.cxx:147
 AliCFUnfolding.cxx:148
 AliCFUnfolding.cxx:149
 AliCFUnfolding.cxx:150
 AliCFUnfolding.cxx:151
 AliCFUnfolding.cxx:152
 AliCFUnfolding.cxx:153
 AliCFUnfolding.cxx:154
 AliCFUnfolding.cxx:155
 AliCFUnfolding.cxx:156
 AliCFUnfolding.cxx:157
 AliCFUnfolding.cxx:158
 AliCFUnfolding.cxx:159
 AliCFUnfolding.cxx:160
 AliCFUnfolding.cxx:161
 AliCFUnfolding.cxx:162
 AliCFUnfolding.cxx:163
 AliCFUnfolding.cxx:164
 AliCFUnfolding.cxx:165
 AliCFUnfolding.cxx:166
 AliCFUnfolding.cxx:167
 AliCFUnfolding.cxx:168
 AliCFUnfolding.cxx:169
 AliCFUnfolding.cxx:170
 AliCFUnfolding.cxx:171
 AliCFUnfolding.cxx:172
 AliCFUnfolding.cxx:173
 AliCFUnfolding.cxx:174
 AliCFUnfolding.cxx:175
 AliCFUnfolding.cxx:176
 AliCFUnfolding.cxx:177
 AliCFUnfolding.cxx:178
 AliCFUnfolding.cxx:179
 AliCFUnfolding.cxx:180
 AliCFUnfolding.cxx:181
 AliCFUnfolding.cxx:182
 AliCFUnfolding.cxx:183
 AliCFUnfolding.cxx:184
 AliCFUnfolding.cxx:185
 AliCFUnfolding.cxx:186
 AliCFUnfolding.cxx:187
 AliCFUnfolding.cxx:188
 AliCFUnfolding.cxx:189
 AliCFUnfolding.cxx:190
 AliCFUnfolding.cxx:191
 AliCFUnfolding.cxx:192
 AliCFUnfolding.cxx:193
 AliCFUnfolding.cxx:194
 AliCFUnfolding.cxx:195
 AliCFUnfolding.cxx:196
 AliCFUnfolding.cxx:197
 AliCFUnfolding.cxx:198
 AliCFUnfolding.cxx:199
 AliCFUnfolding.cxx:200
 AliCFUnfolding.cxx:201
 AliCFUnfolding.cxx:202
 AliCFUnfolding.cxx:203
 AliCFUnfolding.cxx:204
 AliCFUnfolding.cxx:205
 AliCFUnfolding.cxx:206
 AliCFUnfolding.cxx:207
 AliCFUnfolding.cxx:208
 AliCFUnfolding.cxx:209
 AliCFUnfolding.cxx:210
 AliCFUnfolding.cxx:211
 AliCFUnfolding.cxx:212
 AliCFUnfolding.cxx:213
 AliCFUnfolding.cxx:214
 AliCFUnfolding.cxx:215
 AliCFUnfolding.cxx:216
 AliCFUnfolding.cxx:217
 AliCFUnfolding.cxx:218
 AliCFUnfolding.cxx:219
 AliCFUnfolding.cxx:220
 AliCFUnfolding.cxx:221
 AliCFUnfolding.cxx:222
 AliCFUnfolding.cxx:223
 AliCFUnfolding.cxx:224
 AliCFUnfolding.cxx:225
 AliCFUnfolding.cxx:226
 AliCFUnfolding.cxx:227
 AliCFUnfolding.cxx:228
 AliCFUnfolding.cxx:229
 AliCFUnfolding.cxx:230
 AliCFUnfolding.cxx:231
 AliCFUnfolding.cxx:232
 AliCFUnfolding.cxx:233
 AliCFUnfolding.cxx:234
 AliCFUnfolding.cxx:235
 AliCFUnfolding.cxx:236
 AliCFUnfolding.cxx:237
 AliCFUnfolding.cxx:238
 AliCFUnfolding.cxx:239
 AliCFUnfolding.cxx:240
 AliCFUnfolding.cxx:241
 AliCFUnfolding.cxx:242
 AliCFUnfolding.cxx:243
 AliCFUnfolding.cxx:244
 AliCFUnfolding.cxx:245
 AliCFUnfolding.cxx:246
 AliCFUnfolding.cxx:247
 AliCFUnfolding.cxx:248
 AliCFUnfolding.cxx:249
 AliCFUnfolding.cxx:250
 AliCFUnfolding.cxx:251
 AliCFUnfolding.cxx:252
 AliCFUnfolding.cxx:253
 AliCFUnfolding.cxx:254
 AliCFUnfolding.cxx:255
 AliCFUnfolding.cxx:256
 AliCFUnfolding.cxx:257
 AliCFUnfolding.cxx:258
 AliCFUnfolding.cxx:259
 AliCFUnfolding.cxx:260
 AliCFUnfolding.cxx:261
 AliCFUnfolding.cxx:262
 AliCFUnfolding.cxx:263
 AliCFUnfolding.cxx:264
 AliCFUnfolding.cxx:265
 AliCFUnfolding.cxx:266
 AliCFUnfolding.cxx:267
 AliCFUnfolding.cxx:268
 AliCFUnfolding.cxx:269
 AliCFUnfolding.cxx:270
 AliCFUnfolding.cxx:271
 AliCFUnfolding.cxx:272
 AliCFUnfolding.cxx:273
 AliCFUnfolding.cxx:274
 AliCFUnfolding.cxx:275
 AliCFUnfolding.cxx:276
 AliCFUnfolding.cxx:277
 AliCFUnfolding.cxx:278
 AliCFUnfolding.cxx:279
 AliCFUnfolding.cxx:280
 AliCFUnfolding.cxx:281
 AliCFUnfolding.cxx:282
 AliCFUnfolding.cxx:283
 AliCFUnfolding.cxx:284
 AliCFUnfolding.cxx:285
 AliCFUnfolding.cxx:286
 AliCFUnfolding.cxx:287
 AliCFUnfolding.cxx:288
 AliCFUnfolding.cxx:289
 AliCFUnfolding.cxx:290
 AliCFUnfolding.cxx:291
 AliCFUnfolding.cxx:292
 AliCFUnfolding.cxx:293
 AliCFUnfolding.cxx:294
 AliCFUnfolding.cxx:295
 AliCFUnfolding.cxx:296
 AliCFUnfolding.cxx:297
 AliCFUnfolding.cxx:298
 AliCFUnfolding.cxx:299
 AliCFUnfolding.cxx:300
 AliCFUnfolding.cxx:301
 AliCFUnfolding.cxx:302
 AliCFUnfolding.cxx:303
 AliCFUnfolding.cxx:304
 AliCFUnfolding.cxx:305
 AliCFUnfolding.cxx:306
 AliCFUnfolding.cxx:307
 AliCFUnfolding.cxx:308
 AliCFUnfolding.cxx:309
 AliCFUnfolding.cxx:310
 AliCFUnfolding.cxx:311
 AliCFUnfolding.cxx:312
 AliCFUnfolding.cxx:313
 AliCFUnfolding.cxx:314
 AliCFUnfolding.cxx:315
 AliCFUnfolding.cxx:316
 AliCFUnfolding.cxx:317
 AliCFUnfolding.cxx:318
 AliCFUnfolding.cxx:319
 AliCFUnfolding.cxx:320
 AliCFUnfolding.cxx:321
 AliCFUnfolding.cxx:322
 AliCFUnfolding.cxx:323
 AliCFUnfolding.cxx:324
 AliCFUnfolding.cxx:325
 AliCFUnfolding.cxx:326
 AliCFUnfolding.cxx:327
 AliCFUnfolding.cxx:328
 AliCFUnfolding.cxx:329
 AliCFUnfolding.cxx:330
 AliCFUnfolding.cxx:331
 AliCFUnfolding.cxx:332
 AliCFUnfolding.cxx:333
 AliCFUnfolding.cxx:334
 AliCFUnfolding.cxx:335
 AliCFUnfolding.cxx:336
 AliCFUnfolding.cxx:337
 AliCFUnfolding.cxx:338
 AliCFUnfolding.cxx:339
 AliCFUnfolding.cxx:340
 AliCFUnfolding.cxx:341
 AliCFUnfolding.cxx:342
 AliCFUnfolding.cxx:343
 AliCFUnfolding.cxx:344
 AliCFUnfolding.cxx:345
 AliCFUnfolding.cxx:346
 AliCFUnfolding.cxx:347
 AliCFUnfolding.cxx:348
 AliCFUnfolding.cxx:349
 AliCFUnfolding.cxx:350
 AliCFUnfolding.cxx:351
 AliCFUnfolding.cxx:352
 AliCFUnfolding.cxx:353
 AliCFUnfolding.cxx:354
 AliCFUnfolding.cxx:355
 AliCFUnfolding.cxx:356
 AliCFUnfolding.cxx:357
 AliCFUnfolding.cxx:358
 AliCFUnfolding.cxx:359
 AliCFUnfolding.cxx:360
 AliCFUnfolding.cxx:361
 AliCFUnfolding.cxx:362
 AliCFUnfolding.cxx:363
 AliCFUnfolding.cxx:364
 AliCFUnfolding.cxx:365
 AliCFUnfolding.cxx:366
 AliCFUnfolding.cxx:367
 AliCFUnfolding.cxx:368
 AliCFUnfolding.cxx:369
 AliCFUnfolding.cxx:370
 AliCFUnfolding.cxx:371
 AliCFUnfolding.cxx:372
 AliCFUnfolding.cxx:373
 AliCFUnfolding.cxx:374
 AliCFUnfolding.cxx:375
 AliCFUnfolding.cxx:376
 AliCFUnfolding.cxx:377
 AliCFUnfolding.cxx:378
 AliCFUnfolding.cxx:379
 AliCFUnfolding.cxx:380
 AliCFUnfolding.cxx:381
 AliCFUnfolding.cxx:382
 AliCFUnfolding.cxx:383
 AliCFUnfolding.cxx:384
 AliCFUnfolding.cxx:385
 AliCFUnfolding.cxx:386
 AliCFUnfolding.cxx:387
 AliCFUnfolding.cxx:388
 AliCFUnfolding.cxx:389
 AliCFUnfolding.cxx:390
 AliCFUnfolding.cxx:391
 AliCFUnfolding.cxx:392
 AliCFUnfolding.cxx:393
 AliCFUnfolding.cxx:394
 AliCFUnfolding.cxx:395
 AliCFUnfolding.cxx:396
 AliCFUnfolding.cxx:397
 AliCFUnfolding.cxx:398
 AliCFUnfolding.cxx:399
 AliCFUnfolding.cxx:400
 AliCFUnfolding.cxx:401
 AliCFUnfolding.cxx:402
 AliCFUnfolding.cxx:403
 AliCFUnfolding.cxx:404
 AliCFUnfolding.cxx:405
 AliCFUnfolding.cxx:406
 AliCFUnfolding.cxx:407
 AliCFUnfolding.cxx:408
 AliCFUnfolding.cxx:409
 AliCFUnfolding.cxx:410
 AliCFUnfolding.cxx:411
 AliCFUnfolding.cxx:412
 AliCFUnfolding.cxx:413
 AliCFUnfolding.cxx:414
 AliCFUnfolding.cxx:415
 AliCFUnfolding.cxx:416
 AliCFUnfolding.cxx:417
 AliCFUnfolding.cxx:418
 AliCFUnfolding.cxx:419
 AliCFUnfolding.cxx:420
 AliCFUnfolding.cxx:421
 AliCFUnfolding.cxx:422
 AliCFUnfolding.cxx:423
 AliCFUnfolding.cxx:424
 AliCFUnfolding.cxx:425
 AliCFUnfolding.cxx:426
 AliCFUnfolding.cxx:427
 AliCFUnfolding.cxx:428
 AliCFUnfolding.cxx:429
 AliCFUnfolding.cxx:430
 AliCFUnfolding.cxx:431
 AliCFUnfolding.cxx:432
 AliCFUnfolding.cxx:433
 AliCFUnfolding.cxx:434
 AliCFUnfolding.cxx:435
 AliCFUnfolding.cxx:436
 AliCFUnfolding.cxx:437
 AliCFUnfolding.cxx:438
 AliCFUnfolding.cxx:439
 AliCFUnfolding.cxx:440
 AliCFUnfolding.cxx:441
 AliCFUnfolding.cxx:442
 AliCFUnfolding.cxx:443
 AliCFUnfolding.cxx:444
 AliCFUnfolding.cxx:445
 AliCFUnfolding.cxx:446
 AliCFUnfolding.cxx:447
 AliCFUnfolding.cxx:448
 AliCFUnfolding.cxx:449
 AliCFUnfolding.cxx:450
 AliCFUnfolding.cxx:451
 AliCFUnfolding.cxx:452
 AliCFUnfolding.cxx:453
 AliCFUnfolding.cxx:454
 AliCFUnfolding.cxx:455
 AliCFUnfolding.cxx:456
 AliCFUnfolding.cxx:457
 AliCFUnfolding.cxx:458
 AliCFUnfolding.cxx:459
 AliCFUnfolding.cxx:460
 AliCFUnfolding.cxx:461
 AliCFUnfolding.cxx:462
 AliCFUnfolding.cxx:463
 AliCFUnfolding.cxx:464
 AliCFUnfolding.cxx:465
 AliCFUnfolding.cxx:466
 AliCFUnfolding.cxx:467
 AliCFUnfolding.cxx:468
 AliCFUnfolding.cxx:469
 AliCFUnfolding.cxx:470
 AliCFUnfolding.cxx:471
 AliCFUnfolding.cxx:472
 AliCFUnfolding.cxx:473
 AliCFUnfolding.cxx:474
 AliCFUnfolding.cxx:475
 AliCFUnfolding.cxx:476
 AliCFUnfolding.cxx:477
 AliCFUnfolding.cxx:478
 AliCFUnfolding.cxx:479
 AliCFUnfolding.cxx:480
 AliCFUnfolding.cxx:481
 AliCFUnfolding.cxx:482
 AliCFUnfolding.cxx:483
 AliCFUnfolding.cxx:484
 AliCFUnfolding.cxx:485
 AliCFUnfolding.cxx:486
 AliCFUnfolding.cxx:487
 AliCFUnfolding.cxx:488
 AliCFUnfolding.cxx:489
 AliCFUnfolding.cxx:490
 AliCFUnfolding.cxx:491
 AliCFUnfolding.cxx:492
 AliCFUnfolding.cxx:493
 AliCFUnfolding.cxx:494
 AliCFUnfolding.cxx:495
 AliCFUnfolding.cxx:496
 AliCFUnfolding.cxx:497
 AliCFUnfolding.cxx:498
 AliCFUnfolding.cxx:499
 AliCFUnfolding.cxx:500
 AliCFUnfolding.cxx:501
 AliCFUnfolding.cxx:502
 AliCFUnfolding.cxx:503
 AliCFUnfolding.cxx:504
 AliCFUnfolding.cxx:505
 AliCFUnfolding.cxx:506
 AliCFUnfolding.cxx:507
 AliCFUnfolding.cxx:508
 AliCFUnfolding.cxx:509
 AliCFUnfolding.cxx:510
 AliCFUnfolding.cxx:511
 AliCFUnfolding.cxx:512
 AliCFUnfolding.cxx:513
 AliCFUnfolding.cxx:514
 AliCFUnfolding.cxx:515
 AliCFUnfolding.cxx:516
 AliCFUnfolding.cxx:517
 AliCFUnfolding.cxx:518
 AliCFUnfolding.cxx:519
 AliCFUnfolding.cxx:520
 AliCFUnfolding.cxx:521
 AliCFUnfolding.cxx:522
 AliCFUnfolding.cxx:523
 AliCFUnfolding.cxx:524
 AliCFUnfolding.cxx:525
 AliCFUnfolding.cxx:526
 AliCFUnfolding.cxx:527
 AliCFUnfolding.cxx:528
 AliCFUnfolding.cxx:529
 AliCFUnfolding.cxx:530
 AliCFUnfolding.cxx:531
 AliCFUnfolding.cxx:532
 AliCFUnfolding.cxx:533
 AliCFUnfolding.cxx:534
 AliCFUnfolding.cxx:535
 AliCFUnfolding.cxx:536
 AliCFUnfolding.cxx:537
 AliCFUnfolding.cxx:538
 AliCFUnfolding.cxx:539
 AliCFUnfolding.cxx:540
 AliCFUnfolding.cxx:541
 AliCFUnfolding.cxx:542
 AliCFUnfolding.cxx:543
 AliCFUnfolding.cxx:544
 AliCFUnfolding.cxx:545
 AliCFUnfolding.cxx:546
 AliCFUnfolding.cxx:547
 AliCFUnfolding.cxx:548
 AliCFUnfolding.cxx:549
 AliCFUnfolding.cxx:550
 AliCFUnfolding.cxx:551
 AliCFUnfolding.cxx:552
 AliCFUnfolding.cxx:553
 AliCFUnfolding.cxx:554
 AliCFUnfolding.cxx:555
 AliCFUnfolding.cxx:556
 AliCFUnfolding.cxx:557
 AliCFUnfolding.cxx:558
 AliCFUnfolding.cxx:559
 AliCFUnfolding.cxx:560
 AliCFUnfolding.cxx:561
 AliCFUnfolding.cxx:562
 AliCFUnfolding.cxx:563
 AliCFUnfolding.cxx:564
 AliCFUnfolding.cxx:565
 AliCFUnfolding.cxx:566
 AliCFUnfolding.cxx:567
 AliCFUnfolding.cxx:568
 AliCFUnfolding.cxx:569
 AliCFUnfolding.cxx:570
 AliCFUnfolding.cxx:571
 AliCFUnfolding.cxx:572
 AliCFUnfolding.cxx:573
 AliCFUnfolding.cxx:574
 AliCFUnfolding.cxx:575
 AliCFUnfolding.cxx:576
 AliCFUnfolding.cxx:577
 AliCFUnfolding.cxx:578
 AliCFUnfolding.cxx:579
 AliCFUnfolding.cxx:580
 AliCFUnfolding.cxx:581
 AliCFUnfolding.cxx:582
 AliCFUnfolding.cxx:583
 AliCFUnfolding.cxx:584
 AliCFUnfolding.cxx:585
 AliCFUnfolding.cxx:586
 AliCFUnfolding.cxx:587
 AliCFUnfolding.cxx:588
 AliCFUnfolding.cxx:589
 AliCFUnfolding.cxx:590
 AliCFUnfolding.cxx:591
 AliCFUnfolding.cxx:592
 AliCFUnfolding.cxx:593
 AliCFUnfolding.cxx:594
 AliCFUnfolding.cxx:595
 AliCFUnfolding.cxx:596
 AliCFUnfolding.cxx:597
 AliCFUnfolding.cxx:598
 AliCFUnfolding.cxx:599
 AliCFUnfolding.cxx:600
 AliCFUnfolding.cxx:601
 AliCFUnfolding.cxx:602
 AliCFUnfolding.cxx:603
 AliCFUnfolding.cxx:604
 AliCFUnfolding.cxx:605
 AliCFUnfolding.cxx:606
 AliCFUnfolding.cxx:607
 AliCFUnfolding.cxx:608
 AliCFUnfolding.cxx:609
 AliCFUnfolding.cxx:610
 AliCFUnfolding.cxx:611
 AliCFUnfolding.cxx:612
 AliCFUnfolding.cxx:613
 AliCFUnfolding.cxx:614
 AliCFUnfolding.cxx:615
 AliCFUnfolding.cxx:616
 AliCFUnfolding.cxx:617
 AliCFUnfolding.cxx:618
 AliCFUnfolding.cxx:619
 AliCFUnfolding.cxx:620
 AliCFUnfolding.cxx:621
 AliCFUnfolding.cxx:622
 AliCFUnfolding.cxx:623
 AliCFUnfolding.cxx:624
 AliCFUnfolding.cxx:625
 AliCFUnfolding.cxx:626
 AliCFUnfolding.cxx:627
 AliCFUnfolding.cxx:628
 AliCFUnfolding.cxx:629
 AliCFUnfolding.cxx:630
 AliCFUnfolding.cxx:631
 AliCFUnfolding.cxx:632
 AliCFUnfolding.cxx:633
 AliCFUnfolding.cxx:634
 AliCFUnfolding.cxx:635
 AliCFUnfolding.cxx:636
 AliCFUnfolding.cxx:637
 AliCFUnfolding.cxx:638
 AliCFUnfolding.cxx:639
 AliCFUnfolding.cxx:640
 AliCFUnfolding.cxx:641
 AliCFUnfolding.cxx:642
 AliCFUnfolding.cxx:643
 AliCFUnfolding.cxx:644
 AliCFUnfolding.cxx:645
 AliCFUnfolding.cxx:646
 AliCFUnfolding.cxx:647
 AliCFUnfolding.cxx:648
 AliCFUnfolding.cxx:649
 AliCFUnfolding.cxx:650
 AliCFUnfolding.cxx:651
 AliCFUnfolding.cxx:652
 AliCFUnfolding.cxx:653
 AliCFUnfolding.cxx:654
 AliCFUnfolding.cxx:655
 AliCFUnfolding.cxx:656
 AliCFUnfolding.cxx:657
 AliCFUnfolding.cxx:658
 AliCFUnfolding.cxx:659
 AliCFUnfolding.cxx:660
 AliCFUnfolding.cxx:661
 AliCFUnfolding.cxx:662
 AliCFUnfolding.cxx:663
 AliCFUnfolding.cxx:664
 AliCFUnfolding.cxx:665
 AliCFUnfolding.cxx:666
 AliCFUnfolding.cxx:667
 AliCFUnfolding.cxx:668
 AliCFUnfolding.cxx:669
 AliCFUnfolding.cxx:670
 AliCFUnfolding.cxx:671
 AliCFUnfolding.cxx:672
 AliCFUnfolding.cxx:673
 AliCFUnfolding.cxx:674
 AliCFUnfolding.cxx:675
 AliCFUnfolding.cxx:676
 AliCFUnfolding.cxx:677
 AliCFUnfolding.cxx:678
 AliCFUnfolding.cxx:679
 AliCFUnfolding.cxx:680
 AliCFUnfolding.cxx:681
 AliCFUnfolding.cxx:682
 AliCFUnfolding.cxx:683
 AliCFUnfolding.cxx:684
 AliCFUnfolding.cxx:685
 AliCFUnfolding.cxx:686
 AliCFUnfolding.cxx:687
 AliCFUnfolding.cxx:688
 AliCFUnfolding.cxx:689
 AliCFUnfolding.cxx:690
 AliCFUnfolding.cxx:691
 AliCFUnfolding.cxx:692
 AliCFUnfolding.cxx:693
 AliCFUnfolding.cxx:694
 AliCFUnfolding.cxx:695
 AliCFUnfolding.cxx:696
 AliCFUnfolding.cxx:697
 AliCFUnfolding.cxx:698
 AliCFUnfolding.cxx:699
 AliCFUnfolding.cxx:700
 AliCFUnfolding.cxx:701
 AliCFUnfolding.cxx:702
 AliCFUnfolding.cxx:703
 AliCFUnfolding.cxx:704
 AliCFUnfolding.cxx:705
 AliCFUnfolding.cxx:706
 AliCFUnfolding.cxx:707
 AliCFUnfolding.cxx:708
 AliCFUnfolding.cxx:709
 AliCFUnfolding.cxx:710
 AliCFUnfolding.cxx:711
 AliCFUnfolding.cxx:712
 AliCFUnfolding.cxx:713
 AliCFUnfolding.cxx:714
 AliCFUnfolding.cxx:715
 AliCFUnfolding.cxx:716
 AliCFUnfolding.cxx:717
 AliCFUnfolding.cxx:718
 AliCFUnfolding.cxx:719
 AliCFUnfolding.cxx:720
 AliCFUnfolding.cxx:721
 AliCFUnfolding.cxx:722
 AliCFUnfolding.cxx:723
 AliCFUnfolding.cxx:724
 AliCFUnfolding.cxx:725
 AliCFUnfolding.cxx:726
 AliCFUnfolding.cxx:727
 AliCFUnfolding.cxx:728
 AliCFUnfolding.cxx:729
 AliCFUnfolding.cxx:730
 AliCFUnfolding.cxx:731
 AliCFUnfolding.cxx:732
 AliCFUnfolding.cxx:733
 AliCFUnfolding.cxx:734
 AliCFUnfolding.cxx:735
 AliCFUnfolding.cxx:736
 AliCFUnfolding.cxx:737
 AliCFUnfolding.cxx:738
 AliCFUnfolding.cxx:739
 AliCFUnfolding.cxx:740
 AliCFUnfolding.cxx:741
 AliCFUnfolding.cxx:742
 AliCFUnfolding.cxx:743
 AliCFUnfolding.cxx:744
 AliCFUnfolding.cxx:745
 AliCFUnfolding.cxx:746
 AliCFUnfolding.cxx:747
 AliCFUnfolding.cxx:748
 AliCFUnfolding.cxx:749
 AliCFUnfolding.cxx:750
 AliCFUnfolding.cxx:751
 AliCFUnfolding.cxx:752
 AliCFUnfolding.cxx:753
 AliCFUnfolding.cxx:754
 AliCFUnfolding.cxx:755
 AliCFUnfolding.cxx:756
 AliCFUnfolding.cxx:757
 AliCFUnfolding.cxx:758
 AliCFUnfolding.cxx:759
 AliCFUnfolding.cxx:760
 AliCFUnfolding.cxx:761
 AliCFUnfolding.cxx:762
 AliCFUnfolding.cxx:763
 AliCFUnfolding.cxx:764
 AliCFUnfolding.cxx:765
 AliCFUnfolding.cxx:766
 AliCFUnfolding.cxx:767
 AliCFUnfolding.cxx:768
 AliCFUnfolding.cxx:769
 AliCFUnfolding.cxx:770
 AliCFUnfolding.cxx:771
 AliCFUnfolding.cxx:772
 AliCFUnfolding.cxx:773
 AliCFUnfolding.cxx:774
 AliCFUnfolding.cxx:775
 AliCFUnfolding.cxx:776
 AliCFUnfolding.cxx:777
 AliCFUnfolding.cxx:778
 AliCFUnfolding.cxx:779
 AliCFUnfolding.cxx:780
 AliCFUnfolding.cxx:781
 AliCFUnfolding.cxx:782
 AliCFUnfolding.cxx:783
 AliCFUnfolding.cxx:784
 AliCFUnfolding.cxx:785
 AliCFUnfolding.cxx:786
 AliCFUnfolding.cxx:787
 AliCFUnfolding.cxx:788
 AliCFUnfolding.cxx:789
 AliCFUnfolding.cxx:790
 AliCFUnfolding.cxx:791
 AliCFUnfolding.cxx:792
 AliCFUnfolding.cxx:793
 AliCFUnfolding.cxx:794
 AliCFUnfolding.cxx:795
 AliCFUnfolding.cxx:796
 AliCFUnfolding.cxx:797
 AliCFUnfolding.cxx:798
 AliCFUnfolding.cxx:799
 AliCFUnfolding.cxx:800
 AliCFUnfolding.cxx:801
 AliCFUnfolding.cxx:802
 AliCFUnfolding.cxx:803
 AliCFUnfolding.cxx:804
 AliCFUnfolding.cxx:805
 AliCFUnfolding.cxx:806
 AliCFUnfolding.cxx:807
 AliCFUnfolding.cxx:808
 AliCFUnfolding.cxx:809
 AliCFUnfolding.cxx:810
 AliCFUnfolding.cxx:811
 AliCFUnfolding.cxx:812
 AliCFUnfolding.cxx:813
 AliCFUnfolding.cxx:814
 AliCFUnfolding.cxx:815
 AliCFUnfolding.cxx:816
 AliCFUnfolding.cxx:817
 AliCFUnfolding.cxx:818
 AliCFUnfolding.cxx:819
 AliCFUnfolding.cxx:820
 AliCFUnfolding.cxx:821
 AliCFUnfolding.cxx:822
 AliCFUnfolding.cxx:823
 AliCFUnfolding.cxx:824
 AliCFUnfolding.cxx:825
 AliCFUnfolding.cxx:826
 AliCFUnfolding.cxx:827
 AliCFUnfolding.cxx:828
 AliCFUnfolding.cxx:829
 AliCFUnfolding.cxx:830
 AliCFUnfolding.cxx:831
 AliCFUnfolding.cxx:832
 AliCFUnfolding.cxx:833
 AliCFUnfolding.cxx:834
 AliCFUnfolding.cxx:835
 AliCFUnfolding.cxx:836
 AliCFUnfolding.cxx:837
 AliCFUnfolding.cxx:838
 AliCFUnfolding.cxx:839
 AliCFUnfolding.cxx:840
 AliCFUnfolding.cxx:841
 AliCFUnfolding.cxx:842
 AliCFUnfolding.cxx:843
 AliCFUnfolding.cxx:844
 AliCFUnfolding.cxx:845
 AliCFUnfolding.cxx:846
 AliCFUnfolding.cxx:847
 AliCFUnfolding.cxx:848
 AliCFUnfolding.cxx:849
 AliCFUnfolding.cxx:850
 AliCFUnfolding.cxx:851
 AliCFUnfolding.cxx:852
 AliCFUnfolding.cxx:853
 AliCFUnfolding.cxx:854
 AliCFUnfolding.cxx:855
 AliCFUnfolding.cxx:856
 AliCFUnfolding.cxx:857
 AliCFUnfolding.cxx:858
 AliCFUnfolding.cxx:859
 AliCFUnfolding.cxx:860
 AliCFUnfolding.cxx:861
 AliCFUnfolding.cxx:862
 AliCFUnfolding.cxx:863
 AliCFUnfolding.cxx:864