[u/mrichter/AliRoot.git] / ANALYSIS / AliUnfolding.h

/* $Id$ */

#ifndef ALIUNFOLDING_H
#define ALIUNFOLDING_H

//
// class that implements several unfolding methods
// I.e. chi2 minimization and bayesian unfolding
// The whole class is static and not thread-safe (due to the fact that MINUIT unfolding is not thread-safe)
//

// TMatrixD, TVectorD defined here, because it does not seem possible to predeclare these (or i do not know how)
// -->
// $ROOTSYS/include/TVectorDfwd.h:21: conflicting types for `typedef struct TVectorT<Double_t> TVectorD'
// PWG0/AliUnfolding.h:21: previous declaration as `struct TVectorD'

#include "TObject.h"
#include <TMatrixD.h>
#include <TVectorD.h>

class TH1;
class TH2;
class TF1;
class TCanvas;
class TVirtualPad;
class TAxis;

class AliUnfolding : public TObject
{
  public:
  enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature, kRatio, kPowerLaw, kLogLog };
    enum MethodType { kInvalid = -1, kChi2Minimization = 0, kBayesian = 1, kFunction = 2};

    virtual ~AliUnfolding() {};

    static void SetUnfoldingMethod(MethodType methodType);
    static void SetCreateOverflowBin(Float_t overflowBinLimit);
    static void SetSkipBinsBegin(Int_t nBins);
    static void SetNbins(Int_t nMeasured, Int_t nUnfolded);
    static void SetChi2Regularization(RegularizationType type, Float_t weight);
    static void SetMinuitStepSize(Float_t stepSize) { fgMinuitStepSize = stepSize; }
    static void SetMinuitPrecision(Float_t pres) {fgMinuitPrecision = pres;}
    static void SetMinuitMaxIterations(Int_t iter) {fgMinuitMaxIterations = iter;}
    static void SetMinuitStrategy(Double_t strat) {fgMinuitStrategy = strat;}
    static void SetMinimumInitialValue(Bool_t flag, Float_t value = -1) { fgMinimumInitialValue = flag; fgMinimumInitialValueFix = value; }
    static void SetNormalizeInput(Bool_t flag) { fgNormalizeInput = flag; }
    static void SetNotFoundEvents(Float_t notFoundEvents) { fgNotFoundEvents = notFoundEvents; }
    static void SetSkip0BinInChi2(Bool_t flag) { fgSkipBin0InChi2 = flag; }
    static void SetBayesianParameters(Float_t smoothing, Int_t nIterations);
    static void SetFunction(TF1* function);
    static void SetDebug(Bool_t flag) { fgDebug = flag; }
    static void SetPowern(Int_t n) {fgPowern = -1*n;}

    static Int_t Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check = kFALSE);
    static Int_t UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result);

    static TH1* GetPenaltyPlot(Double_t* params);
    static TH1* GetPenaltyPlot(TH1* corrected);

    static TH1* GetResidualsPlot(Double_t* params);
    static TH1* GetResidualsPlot(TH1* corrected);

    static Double_t GetChi2FromFit() {return fChi2FromFit;}
    static Double_t GetPenaltyVal()  {return fPenaltyVal;}
    static Double_t GetAvgResidual() {return fAvgResidual;}

    static void DrawResults(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TCanvas *canvas = 0, Int_t reuseHists = kFALSE,TH1 *unfolded=0);
    static void InteractiveUnfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions);
    static void RedrawInteractive();  

  protected:
    AliUnfolding() {};

    static Int_t UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check);
    static Int_t UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult);
    static Int_t UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult);

    static void CreateOverflowBin(TH2* correlation, TH1* measured); 
    static void SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency);

    static void MakePads();
    static void DrawGuess(Double_t *params, TVirtualPad *pfolded=0, TVirtualPad *pres=0, TVirtualPad *ppen=0, Int_t reuseHists = kFALSE, TH1* unfolded=0);

    static Double_t RegularizationPol0(TVectorD& params);
    static Double_t RegularizationPol1(TVectorD& params);
    static Double_t RegularizationTotalCurvature(TVectorD& params);
    static Double_t RegularizationEntropy(TVectorD& params);
    static Double_t RegularizationLog(TVectorD& params);
    static Double_t RegularizationRatio(TVectorD& params);
    static Double_t RegularizationPowerLaw(TVectorD& params);
    static Double_t RegularizationLogLog(TVectorD& params);

    static void Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
    static void TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);

    // static variable to be accessed by MINUIT
    static TMatrixD* fgCorrelationMatrix;            // contains fCurrentCorrelation in matrix form
    static TMatrixD* fgCorrelationMatrixSquared;     // contains squared fCurrentCorrelation in matrix form
    static TMatrixD* fgCorrelationCovarianceMatrix;  // contains the errors of fCurrentESD
    static TVectorD* fgCurrentESDVector;             // contains fCurrentESD
    static TVectorD* fgEntropyAPriori;               // a-priori distribution for entropy regularization
    static TVectorD* fgEfficiency;                   // efficiency
    /*
    static TVectorD* fgBinWidths;                    // bin widths to be taken into account in regularization
    static TVectorD* fgBinPos;                       // bin positions of unfolded
    */
    static TAxis *fgUnfoldedAxis;                    // bin widths and positions for unfolded
    static TAxis *fgMeasuredAxis;                    // bin widths and positions for measured

    static TF1* fgFitFunction;                       // fit function

    // --- configuration params follow ---
    static MethodType fgMethodType;                  // unfolding method to be used
    static Int_t fgMaxParams;                        // bins in unfolded histogram = number of fit params
    static Int_t fgMaxInput;                         // bins in measured histogram
    static Float_t fgOverflowBinLimit;               // to fix fluctuations at high multiplicities, all entries above the limit are summarized in one bin

    static RegularizationType fgRegularizationType;  // regularization that is used during Chi2 method
    static Float_t fgRegularizationWeight;           // factor for regularization term
    static Int_t fgSkipBinsBegin;                    // (optional) skip the given number of bins in the regularization
    static Float_t fgMinuitStepSize;                 // (usually not needed to be changed) step size in minimization
    static Float_t fgMinuitPrecision;                // precision used by minuit. default = 1e-6
    static Int_t   fgMinuitMaxIterations;            // maximum number of iterations used by minuit. default = 5000
    static Double_t fgMinuitStrategy;                // minuit strategy: 0 (low), 1 (default), 2 (high)
    static Bool_t fgMinimumInitialValue;             // set all initial values at least to the smallest value among the initial values
    static Float_t fgMinimumInitialValueFix;         // use this as the minimum initial value instead of determining it automatically
    static Bool_t fgNormalizeInput;                  // normalize input spectrum
    static Float_t fgNotFoundEvents;                 // constraint on the total number of not found events sum(guess * (1/eff -1))
    static Bool_t fgSkipBin0InChi2;                  // skip bin 0 (= 0 measured) in chi2 function

    static Float_t fgBayesianSmoothing;              // smoothing parameter (0 = no smoothing)
    static Int_t   fgBayesianIterations;             // number of iterations in Bayesian method

    static Bool_t fgDebug;                           // debug flag
    // --- end of configuration ---
    
    static Int_t fgCallCount;                        // call count to chi2 function

    static Int_t fgPowern;                           // power of power law for regularization with power law

    static Double_t fChi2FromFit;                    // Chi2 from fit at current iteration
    static Double_t fPenaltyVal;                     // Penalty value at current iteration (\beta * PU)
    static Double_t fAvgResidual;                    // Sum residuals / nbins

    static Int_t fgPrintChi2Details;                 // debug for chi2 calc

    // Pointers for interactive unfolder
    static TCanvas *fgCanvas;                        // Canvas for interactive unfolder
    static TH1 *fghUnfolded;                         // Unfolding result for interactive unfolder
    static TH2 *fghCorrelation;                      // Response matrix for interactive unfolder
    static TH1 *fghEfficiency;                       // Efficiency histo for interactive unfolder
    static TH1 *fghMeasured;                         // Measured distribution for interactive unfolder

private:
    AliUnfolding(const AliUnfolding&);
    AliUnfolding& operator=(const AliUnfolding&);

  ClassDef(AliUnfolding, 0);
};

#endif
Commit	Line	Data
e6339065	1	/* $Id$ */
	2
	3	#ifndef ALIUNFOLDING_H
	4	#define ALIUNFOLDING_H
	5
	6	//
	7	// class that implements several unfolding methods
19442b86	8	// I.e. chi2 minimization and bayesian unfolding
19442b86	9	// The whole class is static and not thread-safe (due to the fact that MINUIT unfolding is not thread-safe)
e6339065	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;
e6339065	23	class TF1;
9e065ad2	24	class TCanvas;
	25	class TVirtualPad;
	26	class TAxis;
e6339065	27
	28	class AliUnfolding : public TObject
	29	{
	30	public:
9e065ad2	31	enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature, kRatio, kPowerLaw, kLogLog };
19442b86	32	enum MethodType { kInvalid = -1, kChi2Minimization = 0, kBayesian = 1, kFunction = 2};
e6339065	33
19442b86	34	virtual ~AliUnfolding() {};
e6339065	35
19442b86	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; }
9e065ad2	42	static void SetMinuitPrecision(Float_t pres) {fgMinuitPrecision = pres;}
d22beb7f	43	static void SetMinuitMaxIterations(Int_t iter) {fgMinuitMaxIterations = iter;}
9f070ef5	44	static void SetMinuitStrategy(Double_t strat) {fgMinuitStrategy = strat;}
95e970ca	45	static void SetMinimumInitialValue(Bool_t flag, Float_t value = -1) { fgMinimumInitialValue = flag; fgMinimumInitialValueFix = value; }
19442b86	46	static void SetNormalizeInput(Bool_t flag) { fgNormalizeInput = flag; }
95e970ca	47	static void SetNotFoundEvents(Float_t notFoundEvents) { fgNotFoundEvents = notFoundEvents; }
95e970ca	48	static void SetSkip0BinInChi2(Bool_t flag) { fgSkipBin0InChi2 = flag; }
19442b86	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; }
9e065ad2	52	static void SetPowern(Int_t n) {fgPowern = -1*n;}
e6339065	53
19442b86	54	static Int_t Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check = kFALSE);
95e970ca	55	static Int_t UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result);
e6339065	56
95e970ca	57	static TH1* GetPenaltyPlot(Double_t* params);
95e970ca	58	static TH1* GetPenaltyPlot(TH1* corrected);
9e065ad2	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
19442b86	71	protected:
19442b86	72	AliUnfolding() {};
e6339065	73
19442b86	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
19442b86	78	static void CreateOverflowBin(TH2* correlation, TH1* measured);
95e970ca	79	static void SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency);
e6339065	80
9e065ad2	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
e6339065	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);
95e970ca	89	static Double_t RegularizationRatio(TVectorD& params);
9e065ad2	90	static Double_t RegularizationPowerLaw(TVectorD& params);
9e065ad2	91	static Double_t RegularizationLogLog(TVectorD& params);
e6339065	92
19442b86	93	static void Chi2Function(Int_t&, Double_t, Double_t& chi2, Double_t params, Int_t);
19442b86	94	static void TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
e6339065	95
e6339065	96	// static variable to be accessed by MINUIT
19442b86	97	static TMatrixD* fgCorrelationMatrix; // contains fCurrentCorrelation in matrix form
95e970ca	98	static TMatrixD* fgCorrelationMatrixSquared; // contains squared fCurrentCorrelation in matrix form
19442b86	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
95e970ca	102	static TVectorD* fgEfficiency; // efficiency
9e065ad2	103	/*
95e970ca	104	static TVectorD* fgBinWidths; // bin widths to be taken into account in regularization
9e065ad2	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
19442b86	109
19442b86	110	static TF1* fgFitFunction; // fit function
e6339065	111
19442b86	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
e6339065	117
19442b86	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
9e065ad2	122	static Float_t fgMinuitPrecision; // precision used by minuit. default = 1e-6
d22beb7f	123	static Int_t fgMinuitMaxIterations; // maximum number of iterations used by minuit. default = 5000
9f070ef5	124	static Double_t fgMinuitStrategy; // minuit strategy: 0 (low), 1 (default), 2 (high)
95e970ca	125	static Bool_t fgMinimumInitialValue; // set all initial values at least to the smallest value among the initial values
95e970ca	126	static Float_t fgMinimumInitialValueFix; // use this as the minimum initial value instead of determining it automatically
19442b86	127	static Bool_t fgNormalizeInput; // normalize input spectrum
95e970ca	128	static Float_t fgNotFoundEvents; // constraint on the total number of not found events sum(guess * (1/eff -1))
95e970ca	129	static Bool_t fgSkipBin0InChi2; // skip bin 0 (= 0 measured) in chi2 function
e6339065	130
19442b86	131	static Float_t fgBayesianSmoothing; // smoothing parameter (0 = no smoothing)
19442b86	132	static Int_t fgBayesianIterations; // number of iterations in Bayesian method
e6339065	133
19442b86	134	static Bool_t fgDebug; // debug flag
19442b86	135	// --- end of configuration ---
95e970ca	136
95e970ca	137	static Int_t fgCallCount; // call count to chi2 function
e6339065	138
9e065ad2	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
19442b86	154	private:
e6339065	155	AliUnfolding(const AliUnfolding&);
	156	AliUnfolding& operator=(const AliUnfolding&);
	157
	158	ClassDef(AliUnfolding, 0);
	159	};
	160
	161	#endif
	162