[u/mrichter/AliRoot.git] / CORRFW / AliCFUnfolding.h


#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"

class TF1;
class TProfile;
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);
  AliCFUnfolding(const AliCFUnfolding& c);
  AliCFUnfolding& operator= (const AliCFUnfolding& c);
  ~AliCFUnfolding();

  void SetMaxNumberOfIterations(Int_t n = 10)  {fMaxNumIterations=n; fNRandomIterations=n; }

  /* 
     The following is for correct error estimation
     thanks to Marta Verweij
  */
  void SetUseCorrelatedErrors  (Double_t maxConvergence = 1.e-06 , UInt_t randomSeed = 0)   {
    fUseCorrelatedErrors = kTRUE             ;
    fRandomSeed          = randomSeed        ;
    SetMaxConvergencePerDOF(maxConvergence)  ;
  }
  void UnsetCorrelatedErrors() {fUseCorrelatedErrors=kFALSE;}


  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();

  THnSparse* GetResponse()             const {return fResponse;}
  THnSparse* GetInverseResponse()      const {return fInverseResponse;}
  THnSparse* GetPrior()                const {return fPrior;}
  THnSparse* GetOriginalPrior()        const {return fOriginalPrior;}
  THnSparse* GetEfficiency()           const {return fEfficiency;}
  THnSparse* GetUnfolded()             const {return fUnfolded;}
  THnSparse* GetEstMeasured()          const {return fMeasuredEstimate;}
  THnSparse* GetMeasured()             const {return fMeasured;}
  THnSparse* GetConditional()          const {return fConditional;}
  TF1*       GetSmoothFunction()       const {return fSmoothFunction;}
  TProfile*  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 :
  
  // user-related settings
  THnSparse     *fResponse;           // 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
  THnSparse     *fPrior;              // 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
  THnSparse     *fEfficiency;         // Efficiency map : dimensions must be N = number of variables
                                      // this map must be filled only with "true" values of the variables (should not include resolution effects)
  THnSparse     *fMeasured;           // Measured spectrum to be unfolded : dimensions must be N = number of variables (modified)
  THnSparse     *fMeasuredOrig;       // Measured spectrum to be unfolded : dimensions must be N = number of variables (not modified)
  Int_t          fMaxNumIterations;   // Maximum  number of iterations to be performed
  Int_t          fNVariables;         // Number of variables used in analysis spectra (pt, y, ...)
/*   Double_t       fMaxChi2;            // Maximum Chi2 between unfolded and prior distributions.  */
  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"

  /* correlated error calculation */
  Double_t       fMaxConvergence;     // Convergence criterion in case of correlated error calculation
  Bool_t         fUseCorrelatedErrors;// Calculate correlated errors for the final unfolded spectrum; default is kTRUE
  Int_t          fNRandomIterations;  // Number of random distributed measured spectra

  // internal settings
  THnSparse     *fOriginalPrior;     // This is the original prior distribution : will not be 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     *fProjResponseInT;   // Projection of the response matrix on TRUE axis
  THnSparse     *fUnfolded;          // 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     *fRandomizedDist;    // Randomized distribution for each bin of the measured spectrum to calculate correlated errors
  TRandom3      *fRandom3;           // Object to get random number following Poisson distribution
  THnSparse     *fRandomUnfolded;
  TProfile      *fDeltaUnfoldedP;    // Profile of the delta-unfolded distribution
  Int_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     InitDeltaUnfoldedProfile();  // Initializes the fDeltaUnfoldedP Profiles with spread option
  void     CreateRandomizedDist();      // Create randomized dist from measured distribution
  void     FillDeltaUnfoldedProfile();  // Fills the fDeltaUnfoldedP profile
  void     SetMaxConvergencePerDOF (Double_t val);

  ClassDef(AliCFUnfolding,1);
};

#endif
Commit	Line	Data
c0b10ad4	1
	2	#ifndef ALICFUNFOLDING_H
	3	#define ALICFUNFOLDING_H
	4
	5	//--------------------------------------------------------------------//
	6	// //
	7	// AliCFUnfolding Class //
fb494025	8	// Class to handle general unfolding procedure using bayesian method //
c0b10ad4	9	// //
	10	// Author : renaud.vernet@cern.ch //
	11	//--------------------------------------------------------------------//
	12
	13	#include "TNamed.h"
	14	#include "THnSparse.h"
	15
85b6bda9	16	class TF1;
769f5114	17	class TProfile;
769f5114	18	class TRandom3;
85b6bda9	19
c0b10ad4	20	class AliCFUnfolding : public TNamed {
	21
	22	public :
	23	AliCFUnfolding();
	24	AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar,
	25	const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior=0x0);
	26	AliCFUnfolding(const AliCFUnfolding& c);
	27	AliCFUnfolding& operator= (const AliCFUnfolding& c);
	28	~AliCFUnfolding();
	29
769f5114	30	void SetMaxNumberOfIterations(Int_t n = 10) {fMaxNumIterations=n; fNRandomIterations=n; }
	31
	32	/*
	33	The following is for correct error estimation
	34	thanks to Marta Verweij
	35	*/
	36	void SetUseCorrelatedErrors (Double_t maxConvergence = 1.e-06 , UInt_t randomSeed = 0) {
	37	fUseCorrelatedErrors = kTRUE ;
	38	fRandomSeed = randomSeed ;
	39	SetMaxConvergencePerDOF(maxConvergence) ;
	40	}
489078e1	41	void UnsetCorrelatedErrors() {fUseCorrelatedErrors=kFALSE;}
769f5114	42
769f5114	43
85b6bda9	44	void UseSmoothing(TF1* fcn=0x0, Option_t* opt="iremn") { // if fcn=0x0 then smooth using neighbouring bins
	45	fUseSmoothing=kTRUE; // this function must NOT be used if fNVariables > 3
	46	fSmoothFunction=fcn; // the option "opt" is used if "fcn" is specified
	47	fSmoothOption=opt;
	48	}
	49
c0b10ad4	50	void Unfold();
c0b10ad4	51
769f5114	52	THnSparse* GetResponse() const {return fResponse;}
	53	THnSparse* GetInverseResponse() const {return fInverseResponse;}
	54	THnSparse* GetPrior() const {return fPrior;}
	55	THnSparse* GetOriginalPrior() const {return fOriginalPrior;}
	56	THnSparse* GetEfficiency() const {return fEfficiency;}
	57	THnSparse* GetUnfolded() const {return fUnfolded;}
	58	THnSparse* GetEstMeasured() const {return fMeasuredEstimate;}
	59	THnSparse* GetMeasured() const {return fMeasured;}
	60	THnSparse* GetConditional() const {return fConditional;}
	61	TF1* GetSmoothFunction() const {return fSmoothFunction;}
	62	TProfile* GetDeltaUnfoldedProfile() const {return fDeltaUnfoldedP;}
	63	Int_t GetDOF(); // Returns number of degrees of freedom
c0b10ad4	64
7036630f	65	static Short_t SmoothUsingNeighbours(THnSparse*); // smoothes the unfolded spectrum using the neighbouring cells
7036630f	66
c0b10ad4	67	private :
	68
	69	// user-related settings
	70	THnSparse *fResponse; // Response matrix : dimensions must be 2N = 2 x (number of variables)
85b6bda9	71	// dimensions 0 -> N-1 must be filled with reconstructed values
85b6bda9	72	// dimensions N -> 2N-1 must be filled with generated values
c0b10ad4	73	THnSparse *fPrior; // This is the assumed generated distribution : dimensions must be N = number of variables
	74	// it will be used at the first step
	75	// then will be updated automatically at each iteration
c0b10ad4	76	THnSparse *fEfficiency; // Efficiency map : dimensions must be N = number of variables
c0b10ad4	77	// this map must be filled only with "true" values of the variables (should not include resolution effects)
769f5114	78	THnSparse *fMeasured; // Measured spectrum to be unfolded : dimensions must be N = number of variables (modified)
769f5114	79	THnSparse *fMeasuredOrig; // Measured spectrum to be unfolded : dimensions must be N = number of variables (not modified)
c0b10ad4	80	Int_t fMaxNumIterations; // Maximum number of iterations to be performed
c0b10ad4	81	Int_t fNVariables; // Number of variables used in analysis spectra (pt, y, ...)
769f5114	82	/* Double_t fMaxChi2; // Maximum Chi2 between unfolded and prior distributions. */
769f5114	83	Bool_t fUseSmoothing; // Smooth the unfolded sectrum at each iteration; default is kFALSE
85b6bda9	84	TF1 *fSmoothFunction; // Function used to smooth the unfolded spectrum
85b6bda9	85	Option_t *fSmoothOption; // Option to use during the fit (with fSmoothFunction) ; default is "iremn"
c0b10ad4	86
769f5114	87	/* correlated error calculation */
	88	Double_t fMaxConvergence; // Convergence criterion in case of correlated error calculation
	89	Bool_t fUseCorrelatedErrors;// Calculate correlated errors for the final unfolded spectrum; default is kTRUE
	90	Int_t fNRandomIterations; // Number of random distributed measured spectra
	91
c0b10ad4	92	// internal settings
85b6bda9	93	THnSparse *fOriginalPrior; // This is the original prior distribution : will not be modified
c0b10ad4	94	THnSparse *fInverseResponse; // Inverse response matrix
	95	THnSparse *fMeasuredEstimate; // Estimation of the measured (M) spectrum given the a priori (T) distribution
	96	THnSparse *fConditional; // Matrix holding the conditional probabilities P(M\|T)
	97	THnSparse *fProjResponseInT; // Projection of the response matrix on TRUE axis
	98	THnSparse *fUnfolded; // Unfolded spectrum
	99	Int_t *fCoordinates2N; // Coordinates in 2N (measured,true) space
	100	Int_t *fCoordinatesN_M; // Coordinates in measured space
	101	Int_t *fCoordinatesN_T; // Coordinates in true space
769f5114	102
	103
	104	/* correlated error calculation */
	105	THnSparse *fRandomizedDist; // Randomized distribution for each bin of the measured spectrum to calculate correlated errors
	106	TRandom3 *fRandom3; // Object to get random number following Poisson distribution
	107	THnSparse *fRandomUnfolded;
	108	TProfile *fDeltaUnfoldedP; // Profile of the delta-unfolded distribution
	109	Int_t fNCalcCorrErrors; // book keeping to prevend infinite loop
	110	UInt_t fRandomSeed; // Random seed
	111
c0b10ad4	112
c0b10ad4	113	// functions
85b6bda9	114	void Init(); // initialisation of the internal settings
	115	void GetCoordinates(); // gets a cell coordinates in Measured and True space
	116	void CreateConditional(); // creates the conditional matrix from the response matrix
	117	void CreateEstMeasured(); // creates the measured spectrum estimation from the conditional matrix and the prior distribution
	118	void CreateInvResponse(); // creates the inverse response function (Bayes Theorem) from the conditional matrix and the prior distribution
	119	void CreateUnfolded(); // creates the unfolded spectrum from the inverse response matrix and the measured distribution
	120	void CreateFlatPrior(); // creates a flat a priori distribution in case the one given in the constructor is null
	121	Double_t GetChi2(); // returns the chi2 between unfolded and prior spectra
	122	Short_t Smooth(); // function calling smoothing methods
85b6bda9	123	Short_t SmoothUsingFunction(); // smoothes the unfolded spectrum using a fit function
c0b10ad4	124
769f5114	125	/* correlated error calculation */
	126	Double_t GetConvergence(); // Returns convergence criterion
	127	void CalculateCorrelatedErrors(); // Calculates correlated errors for the final unfolded spectrum
	128	void InitDeltaUnfoldedProfile(); // Initializes the fDeltaUnfoldedP Profiles with spread option
	129	void CreateRandomizedDist(); // Create randomized dist from measured distribution
	130	void FillDeltaUnfoldedProfile(); // Fills the fDeltaUnfoldedP profile
	131	void SetMaxConvergencePerDOF (Double_t val);
	132
	133	ClassDef(AliCFUnfolding,1);
c0b10ad4	134	};
	135
	136	#endif