PWG0/AliUnfolding.h

   1 /* $Id$ */
   2
   3 #ifndef ALIUNFOLDING_H
   4 #define ALIUNFOLDING_H
   5
   6 //
   7 // class that implements several unfolding methods
   8 // E.g. chi2 minimization and bayesian unfolding
   9 //
  10
  11 // TMatrixD, TVectorD defined here, because it does not seem possible to predeclare these (or i do not know how)
  12 // -->
  13 // $ROOTSYS/include/TVectorDfwd.h:21: conflicting types for `typedef struct TVectorT<Double_t> TVectorD'
  14 // PWG0/AliUnfolding.h:21: previous declaration as `struct TVectorD'
  15
  16 #include "TObject.h"
  17 #include <TMatrixD.h>
  18 #include <TVectorD.h>
  19
  20 class TH1;
  21 class TH2;
  22 class TH1F;
  23 class TH2F;
  24 class TH3F;
  25 class TF1;
  26 class TCollection;
  27
  28 class AliUnfolding : public TObject
  29 {
  30   public:
  31     enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature };
  32     enum MethodType { kChi2Minimization = 0, kBayesian = 1 };
  33
  34     AliUnfolding();
  35     virtual ~AliUnfolding();
  36
  37     void SetInput(TH2* correlationMatrix, TH1* efficiency, TH1* measured) { fCurrentCorrelation = correlationMatrix; fCurrentEfficiency = efficiency; fCurrentESD = measured; }
  38     void SetInitialConditions(TH1* initialConditions) { fInitialConditions = initialConditions; }
  39     const TH1* GetResult() const { return fResult; }
  40
  41     static void SetParameters(Int_t measuredBins, Int_t unfoldedBins, Bool_t bigbin) { fMaxInput = measuredBins; fMaxParams = unfoldedBins; fgCreateBigBin = bigbin; }
  42     static void SetChi2MinimizationParameters(RegularizationType type, Float_t weight) { fgRegularizationType = type; fgRegularizationWeight = weight; }
  43     static void SetRegularizationRange(Int_t start, Int_t end) { fgRegularizationRangeStart = start; fgRegularizationRangeEnd = end; }
  44     static void SetBayesianParameters(Float_t smoothing, Int_t nIterations) { fgBayesianSmoothing = smoothing; fgBayesianIterations = nIterations; }
  45
  46     Int_t ApplyMinuitFit(Bool_t check = kFALSE);
  47     Int_t ApplyBayesianMethod(Bool_t determineError = kTRUE);
  48     Int_t ApplyNBDFit();
  49     Int_t ApplyLaszloMethod();
  50
  51     TH1* StatisticalUncertainty(MethodType methodType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo = 0);
  52
  53   protected:
  54     static Double_t RegularizationPol0(TVectorD& params);
  55     static Double_t RegularizationPol1(TVectorD& params);
  56     static Double_t RegularizationTotalCurvature(TVectorD& params);
  57     static Double_t RegularizationEntropy(TVectorD& params);
  58     static Double_t RegularizationLog(TVectorD& params);
  59
  60     static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
  61     static void MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
  62
  63     void SetupCurrentHists();
  64
  65     Int_t UnfoldWithBayesian(* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations);
  66     Int_t UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check);
  67
  68     Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u);
  69
  70     TH1* fCurrentESD;         //! measured spectrum
  71     TH2* fCurrentCorrelation; //! correlation matrix
  72     TH1* fCurrentEfficiency;  //! efficiency
  73     TH1* fInitialConditions;  //! initial conditions
  74     TH1* fResult;             //! unfolding result
  75
  76     // static variable to be accessed by MINUIT
  77     static TMatrixD* fgCorrelationMatrix;            //! contains fCurrentCorrelation in matrix form
  78     static TMatrixD* fgCorrelationCovarianceMatrix;  //! contains the errors of fCurrentESD
  79     static TVectorD* fgCurrentESDVector;             //! contains fCurrentESD
  80     static TVectorD* fgEntropyAPriori;               //! a-priori distribution for entropy regularization
  81
  82     static TF1* fgNBD;   //! negative binomial distribution
  83
  84     static Int_t fgMaxParams;  //! bins in unfolded histogram = number of fit params
  85     static Int_t fgMaxInput;   //! bins in measured histogram
  86
  87     // configuration params follow
  88     static RegularizationType fgRegularizationType; //! regularization that is used during Chi2 method
  89     static Float_t fgRegularizationWeight;          //! factor for regularization term
  90     static Int_t fgRegularizationRangeStart;        //! first bin where regularization is applied
  91     static Int_t fgRegularizationRangeEnd;          //! last bin + 1 where regularization is applied
  92     static Bool_t  fgCreateBigBin;                  //! to fix fluctuations at high multiplicities, all entries above a certain limit are summarized in one bin
  93
  94     static Float_t fgBayesianSmoothing;             //! smoothing parameter (0 = no smoothing)
  95     static Int_t   fgBayesianIterations;            //! number of iterations in Bayesian method
  96     // end of configuration
  97
  98  private:
  99     AliUnfolding(const AliUnfolding&);
 100     AliUnfolding& operator=(const AliUnfolding&);
 101
 102   ClassDef(AliUnfolding, 0);
 103 };
 104
 105 #endif
 106