Undoing some unwanted method renames
[u/mrichter/AliRoot.git] / PWGJE / macros / PID / AliTPCPIDmathFit.h
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