//--------------------------------------------------------------------//
// //
// 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 //
+// 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);
- AliCFUnfolding(const AliCFUnfolding& c);
- AliCFUnfolding& operator= (const AliCFUnfolding& c);
+ 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 SetMaxNumberOfIterations(Int_t n) {fMaxNumIterations=n;}
- void SetMaxChi2(Double_t val) {fMaxChi2=val;}
- void SetMaxChi2PerDOF(Double_t val);
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
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;}
+ 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
- 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
- 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
- TF1 *fSmoothFunction; // Function used to smooth the unfolded spectrum
- Option_t *fSmoothOption; // Option to use during the fit (with fSmoothFunction) ; default is "iremn"
+ //
+ 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
- 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
-
+ //
+ 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 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 SmoothUsingNeighbours(); // smoothes the unfolded spectrum using the neighbouring cells
Short_t SmoothUsingFunction(); // smoothes the unfolded spectrum using a fit function
- ClassDef(AliCFUnfolding,0);
+ /* 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
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff