ANALYSIS/AliUnfolding.h

   1 /* $Id$ */
   2
   3 #ifndef ALIUNFOLDING_H
   4 #define ALIUNFOLDING_H
   5
   6 //
   7 // class that implements several unfolding methods
   8 // I.e. chi2 minimization and bayesian unfolding
   9 // The whole class is static and not thread-safe (due to the fact that MINUIT unfolding is not thread-safe)
  10 //
  11
  12 // TMatrixD, TVectorD defined here, because it does not seem possible to predeclare these (or i do not know how)
  13 // -->
  14 // $ROOTSYS/include/TVectorDfwd.h:21: conflicting types for `typedef struct TVectorT<Double_t> TVectorD'
  15 // PWG0/AliUnfolding.h:21: previous declaration as `struct TVectorD'
  16
  17 #include "TObject.h"
  18 #include <TMatrixD.h>
  19 #include <TVectorD.h>
  20
  21 class TH1;
  22 class TH2;
  23 class TF1;
  24 class TCanvas;
  25 class TVirtualPad;
  26 class TAxis;
  27
  28 class AliUnfolding : public TObject
  29 {
  30   public:
  31   enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature, kRatio, kPowerLaw, kLogLog };
  32     enum MethodType { kInvalid = -1, kChi2Minimization = 0, kBayesian = 1, kFunction = 2};
  33
  34     virtual ~AliUnfolding() {};
  35
  36     static void SetUnfoldingMethod(MethodType methodType);
  37     static void SetCreateOverflowBin(Float_t overflowBinLimit);
  38     static void SetSkipBinsBegin(Int_t nBins);
  39     static void SetNbins(Int_t nMeasured, Int_t nUnfolded);
  40     static void SetChi2Regularization(RegularizationType type, Float_t weight);
  41     static void SetMinuitStepSize(Float_t stepSize) { fgMinuitStepSize = stepSize; }
  42     static void SetMinuitPrecision(Float_t pres) {fgMinuitPrecision = pres;}
  43     static void SetMinuitMaxIterations(Int_t iter) {fgMinuitMaxIterations = iter;}
  44     static void SetMinuitStrategy(Double_t strat) {fgMinuitStrategy = strat;}
  45     static void SetMinimumInitialValue(Bool_t flag, Float_t value = -1) { fgMinimumInitialValue = flag; fgMinimumInitialValueFix = value; }
  46     static void SetNormalizeInput(Bool_t flag) { fgNormalizeInput = flag; }
  47     static void SetNotFoundEvents(Float_t notFoundEvents) { fgNotFoundEvents = notFoundEvents; }
  48     static void SetSkip0BinInChi2(Bool_t flag) { fgSkipBin0InChi2 = flag; }
  49     static void SetBayesianParameters(Float_t smoothing, Int_t nIterations);
  50     static void SetFunction(TF1* function);
  51     static void SetDebug(Bool_t flag) { fgDebug = flag; }
  52     static void SetPowern(Int_t n) {fgPowern = -1*n;}
  53
  54     static Int_t Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check = kFALSE);
  55     static Int_t UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result);
  56
  57     static TH1* GetPenaltyPlot(Double_t* params);
  58     static TH1* GetPenaltyPlot(TH1* corrected);
  59
  60     static TH1* GetResidualsPlot(Double_t* params);
  61     static TH1* GetResidualsPlot(TH1* corrected);
  62
  63     static Double_t GetChi2FromFit() {return fChi2FromFit;}
  64     static Double_t GetPenaltyVal()  {return fPenaltyVal;}
  65     static Double_t GetAvgResidual() {return fAvgResidual;}
  66
  67     static void DrawResults(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TCanvas *canvas = 0, Int_t reuseHists = kFALSE,TH1 *unfolded=0);
  68     static void InteractiveUnfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions);
  69     static void RedrawInteractive();
  70
  71   protected:
  72     AliUnfolding() {};
  73
  74     static Int_t UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check);
  75     static Int_t UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult);
  76     static Int_t UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult);
  77
  78     static void CreateOverflowBin(TH2* correlation, TH1* measured);
  79     static void SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency);
  80
  81     static void MakePads();
  82     static void DrawGuess(Double_t *params, TVirtualPad *pfolded=0, TVirtualPad *pres=0, TVirtualPad *ppen=0, Int_t reuseHists = kFALSE, TH1* unfolded=0);
  83
  84     static Double_t RegularizationPol0(TVectorD& params);
  85     static Double_t RegularizationPol1(TVectorD& params);
  86     static Double_t RegularizationTotalCurvature(TVectorD& params);
  87     static Double_t RegularizationEntropy(TVectorD& params);
  88     static Double_t RegularizationLog(TVectorD& params);
  89     static Double_t RegularizationRatio(TVectorD& params);
  90     static Double_t RegularizationPowerLaw(TVectorD& params);
  91     static Double_t RegularizationLogLog(TVectorD& params);
  92
  93     static void Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
  94     static void TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
  95
  96     // static variable to be accessed by MINUIT
  97     static TMatrixD* fgCorrelationMatrix;            // contains fCurrentCorrelation in matrix form
  98     static TMatrixD* fgCorrelationMatrixSquared;     // contains squared fCurrentCorrelation in matrix form
  99     static TMatrixD* fgCorrelationCovarianceMatrix;  // contains the errors of fCurrentESD
 100     static TVectorD* fgCurrentESDVector;             // contains fCurrentESD
 101     static TVectorD* fgEntropyAPriori;               // a-priori distribution for entropy regularization
 102     static TVectorD* fgEfficiency;                   // efficiency
 103     /*
 104     static TVectorD* fgBinWidths;                    // bin widths to be taken into account in regularization
 105     static TVectorD* fgBinPos;                       // bin positions of unfolded
 106     */
 107     static TAxis *fgUnfoldedAxis;                    // bin widths and positions for unfolded
 108     static TAxis *fgMeasuredAxis;                    // bin widths and positions for measured
 109
 110     static TF1* fgFitFunction;                       // fit function
 111
 112     // --- configuration params follow ---
 113     static MethodType fgMethodType;                  // unfolding method to be used
 114     static Int_t fgMaxParams;                        // bins in unfolded histogram = number of fit params
 115     static Int_t fgMaxInput;                         // bins in measured histogram
 116     static Float_t fgOverflowBinLimit;               // to fix fluctuations at high multiplicities, all entries above the limit are summarized in one bin
 117
 118     static RegularizationType fgRegularizationType;  // regularization that is used during Chi2 method
 119     static Float_t fgRegularizationWeight;           // factor for regularization term
 120     static Int_t fgSkipBinsBegin;                    // (optional) skip the given number of bins in the regularization
 121     static Float_t fgMinuitStepSize;                 // (usually not needed to be changed) step size in minimization
 122     static Float_t fgMinuitPrecision;                // precision used by minuit. default = 1e-6
 123     static Int_t   fgMinuitMaxIterations;            // maximum number of iterations used by minuit. default = 5000
 124     static Double_t fgMinuitStrategy;                // minuit strategy: 0 (low), 1 (default), 2 (high)
 125     static Bool_t fgMinimumInitialValue;             // set all initial values at least to the smallest value among the initial values
 126     static Float_t fgMinimumInitialValueFix;         // use this as the minimum initial value instead of determining it automatically
 127     static Bool_t fgNormalizeInput;                  // normalize input spectrum
 128     static Float_t fgNotFoundEvents;                 // constraint on the total number of not found events sum(guess * (1/eff -1))
 129     static Bool_t fgSkipBin0InChi2;                  // skip bin 0 (= 0 measured) in chi2 function
 130
 131     static Float_t fgBayesianSmoothing;              // smoothing parameter (0 = no smoothing)
 132     static Int_t   fgBayesianIterations;             // number of iterations in Bayesian method
 133
 134     static Bool_t fgDebug;                           // debug flag
 135     // --- end of configuration ---
 136
 137     static Int_t fgCallCount;                        // call count to chi2 function
 138
 139     static Int_t fgPowern;                           // power of power law for regularization with power law
 140
 141     static Double_t fChi2FromFit;                    // Chi2 from fit at current iteration
 142     static Double_t fPenaltyVal;                     // Penalty value at current iteration (\beta * PU)
 143     static Double_t fAvgResidual;                    // Sum residuals / nbins
 144
 145     static Int_t fgPrintChi2Details;                 // debug for chi2 calc
 146
 147     // Pointers for interactive unfolder
 148     static TCanvas *fgCanvas;                        // Canvas for interactive unfolder
 149     static TH1 *fghUnfolded;                         // Unfolding result for interactive unfolder
 150     static TH2 *fghCorrelation;                      // Response matrix for interactive unfolder
 151     static TH1 *fghEfficiency;                       // Efficiency histo for interactive unfolder
 152     static TH1 *fghMeasured;                         // Measured distribution for interactive unfolder
 153
 154 private:
 155     AliUnfolding(const AliUnfolding&);
 156     AliUnfolding& operator=(const AliUnfolding&);
 157
 158   ClassDef(AliUnfolding, 0);
 159 };
 160
 161 #endif
 162