ROOT logo

#ifndef ALICFUNFOLDING_H
#define ALICFUNFOLDING_H

//--------------------------------------------------------------------//
//                                                                    //
// AliCFUnfolding Class                                               //
// Class to handle general unfolding procedure using bayesian method  //
//                                                                    //
// Author : renaud.vernet@cern.ch                                     //
//--------------------------------------------------------------------//

#include "TNamed.h"
#include "THnSparse.h"
#include "AliLog.h"

class TF1;
class TRandom3;

class AliCFUnfolding : public TNamed {

 public :

  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=0x0, 
		 Double_t maxConvergencePerDOF = 1.e-06, UInt_t randomSeed = 0,
		 Int_t maxNumIterations = 10);
  ~AliCFUnfolding();
  void UnsetCorrelatedErrors()  {AliError("===================> DEPRECATED <=====================");}
  void SetUseCorrelatedErrors() {AliError("===================> DEPRECATED <=====================");}
  void SetMaxNumberOfIterations(Int_t n = 10)  {
    AliError("===================> DEPRECATED : should be set in constructor <=====================");
    AliError("===================> DEPRECATED : will be removed soon         <=====================");
    fMaxNumIterations = n;
  }

  void SetNRandomIterations(Int_t n = 100) {fNRandomIterations = n;};

  void UseSmoothing(TF1* fcn=0x0, Option_t* opt="iremn") { // if fcn=0x0 then smooth using neighbouring bins 
    fUseSmoothing=kTRUE;                                   // this function must NOT be used if fNVariables > 3
    fSmoothFunction=fcn;                                   // the option "opt" is used if "fcn" is specified
    fSmoothOption=opt;
  } 
                                                                                                
  void Unfold();

  const THnSparse* GetResponse()             const {return fResponseOrig;}
  const THnSparse* GetEfficiency()           const {return fEfficiencyOrig;}
  const THnSparse* GetMeasured()             const {return fMeasuredOrig;}
  const THnSparse* GetOriginalPrior()        const {return fPriorOrig;}
        THnSparse* GetInverseResponse()      const {return fInverseResponse;}
        THnSparse* GetPrior()                const {return fPrior;}
	THnSparse* GetUnfolded()             const {return fUnfoldedFinal;}
        THnSparse* GetEstMeasured()          const {return fMeasuredEstimate;}
	THnSparse* GetConditional()          const {return fConditional;}
	TF1*       GetSmoothFunction()       const {return fSmoothFunction;}
	THnSparse* GetDeltaUnfoldedProfile() const {return fDeltaUnfoldedP;}
	Int_t      GetDOF();                 // Returns number of degrees of freedom

  static Short_t  SmoothUsingNeighbours(THnSparse*); // smoothes the unfolded spectrum using the neighbouring cells

 private :
  AliCFUnfolding(const AliCFUnfolding& c);
  AliCFUnfolding& operator= (const AliCFUnfolding& c);
  
  //
  // user-related settings
  //
  const THnSparse     *fResponseOrig;     // Response matrix : dimensions must be 2N = 2 x (number of variables)
                                          // dimensions 0 ->  N-1 must be filled with reconstructed values
                                          // dimensions N -> 2N-1 must be filled with generated values
  const THnSparse     *fPriorOrig;        // This is the assumed generated distribution : dimensions must be N = number of variables
                                          // it will be used at the first step 
                                          // then will be updated automatically at each iteration
  const THnSparse     *fEfficiencyOrig;   // Efficiency map : dimensions must be N = number of variables (modified)
                                          // this map must be filled only with "true" values of the variables (should not do "bin smearing")
  const THnSparse     *fMeasuredOrig;     // Measured spectrum to be unfolded : dimensions must be N = number of variables (modified)

        Int_t          fMaxNumIterations; // Maximum  number of iterations to be performed
	Int_t          fNVariables;       // Number of variables used in analysis spectra (pt, y, ...)
        Bool_t         fUseSmoothing;     // Smooth the unfolded sectrum at each iteration; default is kFALSE
	TF1           *fSmoothFunction;   // Function used to smooth the unfolded spectrum
	Option_t      *fSmoothOption;     // Option to use during the fit (with fSmoothFunction) ; default is "iremn"

  //
  // internal settings
  //
	Double_t       fMaxConvergence;    // Convergence criterion in case of correlated error calculation
        Int_t          fNRandomIterations; // Number of random distributed measured spectra
	THnSparse     *fResponse;          // Copy of the original response matrix    (modified)
	THnSparse     *fPrior;             // Copy of the original prior spectrum     (modified)
	THnSparse     *fEfficiency;        // Copy of original efficiency             (modified)
	THnSparse     *fMeasured;          // Copy of the original measureed spectrum (modified)
        THnSparse     *fInverseResponse;   // Inverse response matrix
        THnSparse     *fMeasuredEstimate;  // Estimation of the measured (M) spectrum given the a priori (T) distribution
	THnSparse     *fConditional;       // Matrix holding the conditional probabilities P(M|T)
        THnSparse     *fUnfolded;          // Unfolded spectrum (modified before and during error calculation)
	THnSparse     *fUnfoldedFinal;     // Final unfolded spectrum
        Int_t         *fCoordinates2N;     // Coordinates in 2N (measured,true) space
	Int_t         *fCoordinatesN_M;    // Coordinates in measured space
	Int_t         *fCoordinatesN_T;    // Coordinates in true space


  /* correlated error calculation */
  THnSparse     *fRandomResponse;    // Randomized distribution for each bin of the response matrix     to calculate correlated errors
  THnSparse     *fRandomEfficiency;  // Randomized distribution for each bin of the efficiency spectrum to calculate correlated errors
  THnSparse     *fRandomMeasured;    // Randomized distribution for each bin of the measured   spectrum to calculate correlated errors
  TRandom3      *fRandom3;           // Object to get random number following Poisson distribution
  THnSparse     *fDeltaUnfoldedP;    // Profile of the delta-unfolded distribution
  THnSparse     *fDeltaUnfoldedN;    // Entries of the delta-unfolded distribution (count for each bin)
  Short_t        fNCalcCorrErrors;   // Book-keeping to prevend infinite loop
  UInt_t         fRandomSeed;        // Random seed


  // functions
  void     Init();                  // initialisation of the internal settings
  void     GetCoordinates();        // gets a cell coordinates in Measured and True space
  void     CreateConditional();     // creates the conditional matrix from the response matrix
  void     CreateEstMeasured();     // creates the measured spectrum estimation from the conditional matrix and the prior distribution
  void     CreateInvResponse();     // creates the inverse response function (Bayes Theorem) from the conditional matrix and the prior distribution
  void     CreateUnfolded();        // creates the unfolded spectrum from the inverse response matrix and the measured distribution
  void     CreateFlatPrior();       // creates a flat a priori distribution in case the one given in the constructor is null
  Double_t GetChi2();               // returns the chi2 between unfolded and prior spectra
  Short_t  Smooth();                // function calling smoothing methods
  Short_t  SmoothUsingFunction();   // smoothes the unfolded spectrum using a fit function

  /* correlated error calculation */
  Double_t GetConvergence();            // Returns convergence criterion
  void     CalculateCorrelatedErrors(); // Calculates correlated errors for the final unfolded spectrum
  void     CreateRandomizedDist();      // Create randomized dist from measured distribution
  void     FillDeltaUnfoldedProfile();  // Fills the fDeltaUnfoldedP profile
  void     SetMaxConvergencePerDOF (Double_t val);

  ClassDef(AliCFUnfolding,1);
};

#endif
 AliCFUnfolding.h:1
 AliCFUnfolding.h:2
 AliCFUnfolding.h:3
 AliCFUnfolding.h:4
 AliCFUnfolding.h:5
 AliCFUnfolding.h:6
 AliCFUnfolding.h:7
 AliCFUnfolding.h:8
 AliCFUnfolding.h:9
 AliCFUnfolding.h:10
 AliCFUnfolding.h:11
 AliCFUnfolding.h:12
 AliCFUnfolding.h:13
 AliCFUnfolding.h:14
 AliCFUnfolding.h:15
 AliCFUnfolding.h:16
 AliCFUnfolding.h:17
 AliCFUnfolding.h:18
 AliCFUnfolding.h:19
 AliCFUnfolding.h:20
 AliCFUnfolding.h:21
 AliCFUnfolding.h:22
 AliCFUnfolding.h:23
 AliCFUnfolding.h:24
 AliCFUnfolding.h:25
 AliCFUnfolding.h:26
 AliCFUnfolding.h:27
 AliCFUnfolding.h:28
 AliCFUnfolding.h:29
 AliCFUnfolding.h:30
 AliCFUnfolding.h:31
 AliCFUnfolding.h:32
 AliCFUnfolding.h:33
 AliCFUnfolding.h:34
 AliCFUnfolding.h:35
 AliCFUnfolding.h:36
 AliCFUnfolding.h:37
 AliCFUnfolding.h:38
 AliCFUnfolding.h:39
 AliCFUnfolding.h:40
 AliCFUnfolding.h:41
 AliCFUnfolding.h:42
 AliCFUnfolding.h:43
 AliCFUnfolding.h:44
 AliCFUnfolding.h:45
 AliCFUnfolding.h:46
 AliCFUnfolding.h:47
 AliCFUnfolding.h:48
 AliCFUnfolding.h:49
 AliCFUnfolding.h:50
 AliCFUnfolding.h:51
 AliCFUnfolding.h:52
 AliCFUnfolding.h:53
 AliCFUnfolding.h:54
 AliCFUnfolding.h:55
 AliCFUnfolding.h:56
 AliCFUnfolding.h:57
 AliCFUnfolding.h:58
 AliCFUnfolding.h:59
 AliCFUnfolding.h:60
 AliCFUnfolding.h:61
 AliCFUnfolding.h:62
 AliCFUnfolding.h:63
 AliCFUnfolding.h:64
 AliCFUnfolding.h:65
 AliCFUnfolding.h:66
 AliCFUnfolding.h:67
 AliCFUnfolding.h:68
 AliCFUnfolding.h:69
 AliCFUnfolding.h:70
 AliCFUnfolding.h:71
 AliCFUnfolding.h:72
 AliCFUnfolding.h:73
 AliCFUnfolding.h:74
 AliCFUnfolding.h:75
 AliCFUnfolding.h:76
 AliCFUnfolding.h:77
 AliCFUnfolding.h:78
 AliCFUnfolding.h:79
 AliCFUnfolding.h:80
 AliCFUnfolding.h:81
 AliCFUnfolding.h:82
 AliCFUnfolding.h:83
 AliCFUnfolding.h:84
 AliCFUnfolding.h:85
 AliCFUnfolding.h:86
 AliCFUnfolding.h:87
 AliCFUnfolding.h:88
 AliCFUnfolding.h:89
 AliCFUnfolding.h:90
 AliCFUnfolding.h:91
 AliCFUnfolding.h:92
 AliCFUnfolding.h:93
 AliCFUnfolding.h:94
 AliCFUnfolding.h:95
 AliCFUnfolding.h:96
 AliCFUnfolding.h:97
 AliCFUnfolding.h:98
 AliCFUnfolding.h:99
 AliCFUnfolding.h:100
 AliCFUnfolding.h:101
 AliCFUnfolding.h:102
 AliCFUnfolding.h:103
 AliCFUnfolding.h:104
 AliCFUnfolding.h:105
 AliCFUnfolding.h:106
 AliCFUnfolding.h:107
 AliCFUnfolding.h:108
 AliCFUnfolding.h:109
 AliCFUnfolding.h:110
 AliCFUnfolding.h:111
 AliCFUnfolding.h:112
 AliCFUnfolding.h:113
 AliCFUnfolding.h:114
 AliCFUnfolding.h:115
 AliCFUnfolding.h:116
 AliCFUnfolding.h:117
 AliCFUnfolding.h:118
 AliCFUnfolding.h:119
 AliCFUnfolding.h:120
 AliCFUnfolding.h:121
 AliCFUnfolding.h:122
 AliCFUnfolding.h:123
 AliCFUnfolding.h:124
 AliCFUnfolding.h:125
 AliCFUnfolding.h:126
 AliCFUnfolding.h:127
 AliCFUnfolding.h:128
 AliCFUnfolding.h:129
 AliCFUnfolding.h:130
 AliCFUnfolding.h:131
 AliCFUnfolding.h:132
 AliCFUnfolding.h:133
 AliCFUnfolding.h:134
 AliCFUnfolding.h:135
 AliCFUnfolding.h:136
 AliCFUnfolding.h:137
 AliCFUnfolding.h:138