[u/mrichter/AliRoot.git] / PWGJE / macros / PID / AliTPCPIDmathFit.h

// Class for advanced fitting methods (template fits, (weighted) LL fits, regularised fits, simultaneous fits)
//
// Based on original code of: 
// Xianguo Lu
// lu@physi.uni-heidelberg.de 
//
// Modified by:
// Benjamin Hess
// Benjamin-Andreas.Hess@Uni-Tuebingen.de

#ifndef ALITPCPIDMATHFIT_H
#define ALITPCPIDMATHFIT_H

class TString;
class TMinuit;
class TH1;

class AliTPCPIDmathFit: public TObject
{
public:
  typedef Double_t (*FitFunc_t)(const Double_t *, const Double_t *);
  
  static AliTPCPIDmathFit* Instance(const Int_t numXbinsRegularisation = -1, const Int_t numSimultaneousFits = -1,
                                    const Int_t maxDataPoints = 200000);
  
  Int_t AddRefHisto(TH1* refHisto);
  void ClearRefHistos();
  
  Int_t GetDebugLevel() const { return fDebugLevel; }
  Double_t GetEpsilon() const { return fEpsilon;  }
  Int_t GetNrefHistos() const { return fNrefHistos;  }
  TH1* GetRefHisto(Int_t index) const;
  Int_t GetIndexParametersToRegularise(Int_t index) const;
  Int_t GetMaxCalls() const { return fMaxCalls; }
  const TString GetMinimisationStrategy() const { return fMinimisationString; }; 
  Int_t GetNumParametersPerXbin() const { return fNumParametersPerXbin; };
  Int_t GetNumParametersToRegularise() const { return fNumParametersToRegularise; };
  Int_t GetNumSimultaneousFits() const { return fNumSimultaneousFits; };
  Int_t GetNumXbinsRegularisation() const { return fNumXbinsRegularisation; };
  Int_t GetXbinIndex() const { return fXbinIndex; };
  const Double_t* GetXstatisticalWeight() const { return fXstatisticalWeight; };
  const Double_t* GetXstatisticalWeightError() const { return fXstatisticalWeightError; };
  const Double_t* GetXvaluesForRegularisation() const { return fXvaluesForRegularisation; };
  Int_t GetRegularisation() const { return fRegularisation; };
  Double_t GetRegularisationFactor() const { return fRegularisationFactor; };
  Double_t GetScaleFactorError() const { return fScaleFactorError; };
  Bool_t GetUseLogLikelihood() const { return fUseLogLikelihood; }
  Bool_t GetUseRegularisation() const { return fRegularisation > 0 && fNumParametersToRegularise > 0 && fNumXbinsRegularisation > 1; };
  Bool_t GetUseWeightsForLogLikelihood() const { return fUseWeightsForLoglikelihood; }
  
  void InputData(const TH1 *hh, const Int_t indexXbinRegularisation, const Int_t indexSimultaneousFit, const Double_t leftBoundary,
                 const Double_t rightBoundary, const Double_t threshold, const Bool_t kXerr);
  Int_t MinuitFit(FitFunc_t *fn, FitFunc_t *fnError, const Int_t nPar, Double_t *par, Double_t *err, Double_t *covMatrix, 
                  Double_t &ChiSquareOrLogLikelihood, Int_t &NDF, const Double_t *stepSize = 0x0, const Double_t *lowLim = 0x0, 
                  const Double_t *hiLim = 0x0);
  
  void SetDebugLevel(const Int_t level = 1) { fDebugLevel = level; }
  void SetEpsilon(const Double_t eps);
  void SetMaxCalls(const Int_t maxCalls);
  void SetMinimisationStrategy(TString strategy) { fMinimisationString = strategy; };
  Bool_t SetParametersToRegularise(const Int_t numParams, const Int_t numParamsPerXbin, const Int_t* indexParametersToRegularise,
                                   const Int_t* lastNotFixedIndexOfParameters, const Double_t* xValuesForRegularisation,
                                   const Double_t* xStatisticalWeight, const Double_t* xStatisticalWeightError);
  void SetRegularisation(const Int_t reg, const Double_t regFactor) { fRegularisation = reg; fRegularisationFactor = regFactor; };
  void SetScaleFactorError(const Double_t scaleFactorError) { fScaleFactorError = scaleFactorError; };
  void SetUseLogLikelihood(const Bool_t useLogLikelihood = kTRUE) { fUseLogLikelihood = useLogLikelihood; }
  void SetUseWeightsForLogLikelihood(const Bool_t useWeightsForLogLikelihood = kTRUE) 
    { fUseWeightsForLoglikelihood = useWeightsForLogLikelihood; }
  
 private:
  AliTPCPIDmathFit(const Int_t numXbinsRegularisation = 1, const Int_t numSimultaneousFits = 1, const Int_t maxDataPoints = 200000);
  virtual ~AliTPCPIDmathFit();
  
  void AddPenaltyTermForRegularisation(const Bool_t useLogLikelihood);
  inline Double_t EvalLog(Double_t x) const;
  Int_t Fit(const Double_t *inPar = 0x0, Double_t *covMatrix = 0x0, const Double_t *stepSize = 0x0, const Double_t *lowLim = 0x0,
            const Double_t *hiLim = 0x0);
  static void FitFCN(Int_t &nPar, Double_t *input, Double_t &chiSquareOrLogLikelihood, Double_t *par, Int_t flag);
  Double_t GetChiSquare();
  Double_t GetNegLogLikelihood();
  void GetCovarianceMatrix(Double_t* covMatrix) const;
  Bool_t PolynomialInterpolation(Double_t xa[], Double_t ya[], Int_t n, Double_t x, Double_t* y, Double_t* dy);

  //______________________________________________
  //______________________________________________
  
  static AliTPCPIDmathFit* fgInstance;   //  Instance of this class (singleton implementation)
  
  const Double_t fkDoubleEpsilonLimit; // Limit for comparing double with zero
  
  //------ Initialized in MinuitFit(), func, nPar, par and err passed in from outside ------
  //------ all reset after Fit() ------
  FitFunc_t *fFunc;      // Function used for fitting
  FitFunc_t *fErrorFunc; // Function used for error estimation while fitting
  Int_t fNpar;           // Number of fit parameters
  Double_t *fPar;        // Fit parameters
  Double_t *fErr;        // Error of fit parameters
  
  //chi2, fNegLogLikelihood, ndf only used internally, since if it stays, it may be confused by different MinuitFit.
  //Reset after Fit() and only passed out after a Fit()
  Double_t fChi2;             // Chi^2 of fit
  Double_t fNegLogLikelihood; // Negative log-likelihood of fit
  Int_t fNDF;                 // NDF of fit

  Int_t fNstep;               // Current fitting step
  
  //------ initialized via InputData() ------
  Int_t **fNdata;             // Array containing number of data points
  
  Double_t ***fX;             // Array containing x values of data points
  Double_t ***fY;             // Array containing y values of data points
  
  Double_t ***fXerr;          // Array containing x errors of data points
  Double_t ***fYerr;          // Array containing y errors of data points
  
  //------ Set via setter functions -----
  Int_t fMaxCalls;            // Maximum number of calls for fitting
  Double_t fEpsilon;          // Epsilon value for fit convergence
  Bool_t fUseLogLikelihood;                          // Flag indicating whether to use LL fit or chi^2 fit
  Bool_t fUseWeightsForLoglikelihood;                // Flag indicating whether to use weighted LL fit or normal one (global)
  Bool_t fUseWeightsForLoglikelihoodForCurrentFit;   // Flag indicating whether to use weighted LL fit or normal one for the current fit
  
  Double_t fScaleFactorError; // Scale errors with this factor
  
  Int_t fRegularisation;      // Do regularisation with +/- this number of bins 
  Int_t fXbinIndex;           // Index of current x bin
  
  Int_t fNumParametersToRegularise;       // Number of fit parameter that will be regularised
  Int_t fNumParametersPerXbin;            // Number of fit parameters for each x bin
  Int_t *fIndexParametersToRegularise;    // Indicices of fit parameters that are to be regularised
  Int_t *fLastNotFixedIndexOfParameters;  // Indices of last free parameter in the corresponding x bin
  
  Double_t *fXvaluesForRegularisation;    // Array with x values that will be used for regularisation
  
  Double_t *fXstatisticalWeight;          // Statistical weights of each x bin
  Double_t *fXstatisticalWeightError;     // Errors of statistical weights of each x bin
  
  Double_t fRegularisationFactor;         // Relative weighting factor for regularisation (1 = weighted with error)
  
  TString fMinimisationString;            // Minimisation method

  //------ Initialized via constructor -----
  const Int_t fkMaxNdata;                 // Maximum number of data points
  Int_t fNumSimultaneousFits;             // Number of simultaneous fits
  Int_t fNumXbinsRegularisation;          // Number of x bins for the regularisation
  
  //------ initialized and destroyed in Fit()
  TMinuit *fMinuit;                       // Pointer to minuit object
  
  Double_t *fPolIntX;                     // Used for internal interpolation during regularisation
  Double_t *fPolIntY;                     // Used for internal interpolation during regularisation
  Double_t *fPolIntC;                     // Used for internal interpolation during regularisation
  Double_t *fPolIntD;                     // Used for internal interpolation during regularisation
  
  //------
  Int_t fDebugLevel;                      // Debug level
  
  
  // Reference histos used for the fitting
  TH1** fRefHistos;                       // Array with histos (templates) used for fitting
  Int_t fNrefHistos;                      // Number of histos (templates) for fitting
  
  AliTPCPIDmathFit(const AliTPCPIDmathFit&);            // Not implemented
  AliTPCPIDmathFit& operator=(const AliTPCPIDmathFit&); // Not implemented
  
  ClassDef(AliTPCPIDmathFit, 1);
};

#endif
Commit	Line	Data
e131b05f	1	// Class for advanced fitting methods (template fits, (weighted) LL fits, regularised fits, simultaneous fits)
	2	//
	3	// Based on original code of:
	4	// Xianguo Lu
	5	// lu@physi.uni-heidelberg.de
	6	//
	7	// Modified by:
	8	// Benjamin Hess
	9	// Benjamin-Andreas.Hess@Uni-Tuebingen.de
	10
	11	#ifndef ALITPCPIDMATHFIT_H
	12	#define ALITPCPIDMATHFIT_H
	13
	14	class TString;
	15	class TMinuit;
	16	class TH1;
	17
	18	class AliTPCPIDmathFit: public TObject
	19	{
	20	public:
	21	typedef Double_t (FitFunc_t)(const Double_t , const Double_t *);
	22
	23	static AliTPCPIDmathFit* Instance(const Int_t numXbinsRegularisation = -1, const Int_t numSimultaneousFits = -1,
	24	const Int_t maxDataPoints = 200000);
	25
	26	Int_t AddRefHisto(TH1* refHisto);
	27	void ClearRefHistos();
	28
	29	Int_t GetDebugLevel() const { return fDebugLevel; }
	30	Double_t GetEpsilon() const { return fEpsilon; }
	31	Int_t GetNrefHistos() const { return fNrefHistos; }
	32	TH1* GetRefHisto(Int_t index) const;
	33	Int_t GetIndexParametersToRegularise(Int_t index) const;
	34	Int_t GetMaxCalls() const { return fMaxCalls; }
	35	const TString GetMinimisationStrategy() const { return fMinimisationString; };
	36	Int_t GetNumParametersPerXbin() const { return fNumParametersPerXbin; };
	37	Int_t GetNumParametersToRegularise() const { return fNumParametersToRegularise; };
	38	Int_t GetNumSimultaneousFits() const { return fNumSimultaneousFits; };
	39	Int_t GetNumXbinsRegularisation() const { return fNumXbinsRegularisation; };
	40	Int_t GetXbinIndex() const { return fXbinIndex; };
	41	const Double_t* GetXstatisticalWeight() const { return fXstatisticalWeight; };
	42	const Double_t* GetXstatisticalWeightError() const { return fXstatisticalWeightError; };
	43	const Double_t* GetXvaluesForRegularisation() const { return fXvaluesForRegularisation; };
	44	Int_t GetRegularisation() const { return fRegularisation; };
	45	Double_t GetRegularisationFactor() const { return fRegularisationFactor; };
	46	Double_t GetScaleFactorError() const { return fScaleFactorError; };
	47	Bool_t GetUseLogLikelihood() const { return fUseLogLikelihood; }
	48	Bool_t GetUseRegularisation() const { return fRegularisation > 0 && fNumParametersToRegularise > 0 && fNumXbinsRegularisation > 1; };
	49	Bool_t GetUseWeightsForLogLikelihood() const { return fUseWeightsForLoglikelihood; }
	50
	51	void InputData(const TH1 *hh, const Int_t indexXbinRegularisation, const Int_t indexSimultaneousFit, const Double_t leftBoundary,
	52	const Double_t rightBoundary, const Double_t threshold, const Bool_t kXerr);
	53	Int_t MinuitFit(FitFunc_t fn, FitFunc_t fnError, const Int_t nPar, Double_t par, Double_t err, Double_t *covMatrix,
	54	Double_t &ChiSquareOrLogLikelihood, Int_t &NDF, const Double_t stepSize = 0x0, const Double_t lowLim = 0x0,
	55	const Double_t *hiLim = 0x0);
	56
	57	void SetDebugLevel(const Int_t level = 1) { fDebugLevel = level; }
	58	void SetEpsilon(const Double_t eps);
	59	void SetMaxCalls(const Int_t maxCalls);
	60	void SetMinimisationStrategy(TString strategy) { fMinimisationString = strategy; };
	61	Bool_t SetParametersToRegularise(const Int_t numParams, const Int_t numParamsPerXbin, const Int_t* indexParametersToRegularise,
	62	const Int_t* lastNotFixedIndexOfParameters, const Double_t* xValuesForRegularisation,
	63	const Double_t* xStatisticalWeight, const Double_t* xStatisticalWeightError);
	64	void SetRegularisation(const Int_t reg, const Double_t regFactor) { fRegularisation = reg; fRegularisationFactor = regFactor; };
65	void SetScaleFactorError(const Double_t scaleFactorError) { fScaleFactorError = scaleFactorError; };
66	void SetUseLogLikelihood(const Bool_t useLogLikelihood = kTRUE) { fUseLogLikelihood = useLogLikelihood; }
67	void SetUseWeightsForLogLikelihood(const Bool_t useWeightsForLogLikelihood = kTRUE)
68	{ fUseWeightsForLoglikelihood = useWeightsForLogLikelihood; }
69
70	private:
71	AliTPCPIDmathFit(const Int_t numXbinsRegularisation = 1, const Int_t numSimultaneousFits = 1, const Int_t maxDataPoints = 200000);
72	virtual ~AliTPCPIDmathFit();
73
74	void AddPenaltyTermForRegularisation(const Bool_t useLogLikelihood);
75	inline Double_t EvalLog(Double_t x) const;
76	Int_t Fit(const Double_t inPar = 0x0, Double_t covMatrix = 0x0, const Double_t stepSize = 0x0, const Double_t lowLim = 0x0,
77	const Double_t *hiLim = 0x0);
78	static void FitFCN(Int_t &nPar, Double_t input, Double_t &chiSquareOrLogLikelihood, Double_t par, Int_t flag);
79	Double_t GetChiSquare();
80	Double_t GetNegLogLikelihood();
81	void GetCovarianceMatrix(Double_t* covMatrix) const;
82	Bool_t PolynomialInterpolation(Double_t xa[], Double_t ya[], Int_t n, Double_t x, Double_t* y, Double_t* dy);
83
84	//______________________________________________
85	//______________________________________________
86
87	static AliTPCPIDmathFit* fgInstance; // Instance of this class (singleton implementation)
88
89	const Double_t fkDoubleEpsilonLimit; // Limit for comparing double with zero
90
91	//------ Initialized in MinuitFit(), func, nPar, par and err passed in from outside ------
92	//------ all reset after Fit() ------
93	FitFunc_t *fFunc; // Function used for fitting
94	FitFunc_t *fErrorFunc; // Function used for error estimation while fitting
95	Int_t fNpar; // Number of fit parameters
96	Double_t *fPar; // Fit parameters
97	Double_t *fErr; // Error of fit parameters
98
99	//chi2, fNegLogLikelihood, ndf only used internally, since if it stays, it may be confused by different MinuitFit.
100	//Reset after Fit() and only passed out after a Fit()
101	Double_t fChi2; // Chi^2 of fit
102	Double_t fNegLogLikelihood; // Negative log-likelihood of fit
103	Int_t fNDF; // NDF of fit
104
105	Int_t fNstep; // Current fitting step
106
107	//------ initialized via InputData() ------
108	Int_t **fNdata; // Array containing number of data points
109
110	Double_t ***fX; // Array containing x values of data points
111	Double_t ***fY; // Array containing y values of data points
112
113	Double_t ***fXerr; // Array containing x errors of data points
114	Double_t ***fYerr; // Array containing y errors of data points
115
116	//------ Set via setter functions -----
117	Int_t fMaxCalls; // Maximum number of calls for fitting
118	Double_t fEpsilon; // Epsilon value for fit convergence
119	Bool_t fUseLogLikelihood; // Flag indicating whether to use LL fit or chi^2 fit
120	Bool_t fUseWeightsForLoglikelihood; // Flag indicating whether to use weighted LL fit or normal one (global)
121	Bool_t fUseWeightsForLoglikelihoodForCurrentFit; // Flag indicating whether to use weighted LL fit or normal one for the current fit
122
123	Double_t fScaleFactorError; // Scale errors with this factor
124
125	Int_t fRegularisation; // Do regularisation with +/- this number of bins
126	Int_t fXbinIndex; // Index of current x bin
127
128	Int_t fNumParametersToRegularise; // Number of fit parameter that will be regularised
129	Int_t fNumParametersPerXbin; // Number of fit parameters for each x bin
130	Int_t *fIndexParametersToRegularise; // Indicices of fit parameters that are to be regularised
131	Int_t *fLastNotFixedIndexOfParameters; // Indices of last free parameter in the corresponding x bin
132
133	Double_t *fXvaluesForRegularisation; // Array with x values that will be used for regularisation
134
135	Double_t *fXstatisticalWeight; // Statistical weights of each x bin
136	Double_t *fXstatisticalWeightError; // Errors of statistical weights of each x bin
137
138	Double_t fRegularisationFactor; // Relative weighting factor for regularisation (1 = weighted with error)
139
140	TString fMinimisationString; // Minimisation method
141
142	//------ Initialized via constructor -----
143	const Int_t fkMaxNdata; // Maximum number of data points
144	Int_t fNumSimultaneousFits; // Number of simultaneous fits
145	Int_t fNumXbinsRegularisation; // Number of x bins for the regularisation
146
147	//------ initialized and destroyed in Fit()
148	TMinuit *fMinuit; // Pointer to minuit object
149
150	Double_t *fPolIntX; // Used for internal interpolation during regularisation
151	Double_t *fPolIntY; // Used for internal interpolation during regularisation
152	Double_t *fPolIntC; // Used for internal interpolation during regularisation
153	Double_t *fPolIntD; // Used for internal interpolation during regularisation
154
155	//------
156	Int_t fDebugLevel; // Debug level
157
158
159	// Reference histos used for the fitting
160	TH1** fRefHistos; // Array with histos (templates) used for fitting
161	Int_t fNrefHistos; // Number of histos (templates) for fitting
162
163	AliTPCPIDmathFit(const AliTPCPIDmathFit&); // Not implemented
164	AliTPCPIDmathFit& operator=(const AliTPCPIDmathFit&); // Not implemented
165
166	ClassDef(AliTPCPIDmathFit, 1);
167	};
168
169	#endif