CORRFW/AliCFUnfolding.h

   1
   2 #ifndef ALICFUNFOLDING_H
   3 #define ALICFUNFOLDING_H
   4
   5 //--------------------------------------------------------------------//
   6 //                                                                    //
   7 // AliCFUnfolding Class                                               //
   8 // Class to handle general unfolding procedure using bayesian method  //
   9 //                                                                    //
  10 // Author : renaud.vernet@cern.ch                                     //
  11 //--------------------------------------------------------------------//
  12
  13 #include "TNamed.h"
  14 #include "THnSparse.h"
  15
  16 class TF1;
  17 class TProfile;
  18 class TRandom3;
  19
  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
  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   }
  41   void UnsetCorrelatedErrors() {fUseCorrelatedErrors=kFALSE;}
  42
  43
  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
  50   void Unfold();
  51
  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
  64
  65   static Short_t  SmoothUsingNeighbours(THnSparse*); // smoothes the unfolded spectrum using the neighbouring cells
  66
  67  private :
  68
  69   // user-related settings
  70   THnSparse     *fResponse;           // Response matrix : dimensions must be 2N = 2 x (number of variables)
  71                                       // dimensions 0 ->  N-1 must be filled with reconstructed values
  72                                       // dimensions N -> 2N-1 must be filled with generated values
  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
  76   THnSparse     *fEfficiency;         // Efficiency map : dimensions must be N = number of variables
  77                                       // this map must be filled only with "true" values of the variables (should not include resolution effects)
  78   THnSparse     *fMeasured;           // Measured spectrum to be unfolded : dimensions must be N = number of variables (modified)
  79   THnSparse     *fMeasuredOrig;       // Measured spectrum to be unfolded : dimensions must be N = number of variables (not modified)
  80   Int_t          fMaxNumIterations;   // Maximum  number of iterations to be performed
  81   Int_t          fNVariables;         // Number of variables used in analysis spectra (pt, y, ...)
  82 /*   Double_t       fMaxChi2;            // Maximum Chi2 between unfolded and prior distributions.  */
  83   Bool_t         fUseSmoothing;       // Smooth the unfolded sectrum at each iteration; default is kFALSE
  84   TF1           *fSmoothFunction;     // Function used to smooth the unfolded spectrum
  85   Option_t      *fSmoothOption;       // Option to use during the fit (with fSmoothFunction) ; default is "iremn"
  86
  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
  92   // internal settings
  93   THnSparse     *fOriginalPrior;     // This is the original prior distribution : will not be modified
  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
 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
 112
 113   // functions
 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
 123   Short_t  SmoothUsingFunction();   // smoothes the unfolded spectrum using a fit function
 124
 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);
 134 };
 135
 136 #endif