b073c7157d58be03472aafe9f46987831825c3ed
[u/mrichter/AliRoot.git] / STAT / TKDInterpolator.h
1 #ifndef ROOT_TKDInterpolator
2 #define ROOT_TKDInterpolator
3
4 #ifndef ROOT_TKDTree
5 #include "TKDTree.h"
6 #endif
7
8 #ifndef ROOT_TVectorD
9 #include "TVectorD.h"
10 #endif
11 #ifndef ROOT_TMatrixD
12 #include "TMatrixD.h"
13 #endif
14
15 // Non parametric interpolation class based on local polinomial regression.
16 // The class can be used to approximate PDF together with TKDTree or for
17 // general regression when the data points are given.
18
19 // template <typename Value> class TVectorT;
20 // typedef struct TVectorT<Double_t> TVectorD;
21 // template <typename Value> class TMatrixT;
22 // typedef class TMatrixT<Double_t> TMatrixD;
23 // template <typename Index, typename Value> class TKDTree;
24 // typedef class TKDTree<Int_t, Float_t> TKDTreeIF;
25 class TTree;
26 class TLinearFitter;
27 class TKDInterpolator : public TKDTreeIF
28 {
29 public:
30         TKDInterpolator();
31         TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut = 0, UInt_t bsize = 100, Long64_t nentries = 1000000000, Long64_t firstentry = 0);
32         TKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data);
33         ~TKDInterpolator();
34
35                 Double_t   Eval(const Double_t *point, Double_t &result, Double_t &error, Bool_t force = kFALSE);
36                                         Float_t    GetAlpha() const {return fAlpha;}
37         inline Bool_t     GetCOGPoint(Int_t node, Float_t *coord, Float_t &val, Float_t &error) const;
38         inline Bool_t     GetDataPoint(Int_t n, Float_t *p) const;
39                 TKDTreeIF* GetHelper() {return fKDhelper;}
40                 Bool_t     GetInterpolationMethod() const {return fStatus&1;}
41                 Int_t      GetNTNodes() const {return fNTNodes;}
42                 Bool_t     GetStore() const {return fStatus&2;}
43                 Bool_t     GetWeights() const {return fStatus&4;}
44                                         
45                                         void       DrawNodes(UInt_t ax1 = 0, UInt_t ax2 = 1, Int_t depth = -1);
46                 void       DrawNode(Int_t tnode, UInt_t ax1 = 0, UInt_t ax2 = 1);
47                 void       GetStatus();
48                                         void       SetAlpha(const Float_t a){if(a>0.) fAlpha = a;}
49                 void       SetInterpolationMethod(const Bool_t on = kTRUE);
50                 void       SetStore(const Bool_t on = kTRUE);
51                 void       SetWeights(const Bool_t on = kTRUE);
52         
53 private:
54         TKDInterpolator(const TKDInterpolator &);
55         TKDInterpolator& operator=(const TKDInterpolator &);    
56                 void       Build();
57                                         
58 public:
59         class TKDNodeInfo
60         {
61         public:
62                 TKDNodeInfo(const Int_t ndim = 0);
63                 virtual  ~TKDNodeInfo();
64                 void      Build(const Int_t ndim);
65           Double_t  CookPDF(const Double_t *point, Double_t &result, Double_t &error);
66                 void      Store(const TVectorD &par, const TMatrixD &cov);
67         
68                 Int_t     fNDim;         // data dimension
69                 Float_t   *fRefPoint;    //[fNDim] node's COG
70                 Float_t   fRefValue;     // measured value for node 
71                 TMatrixD  *fCov;         // interpolator covariance matrix
72                 TVectorD  *fPar;         // interpolator parameters
73
74         private:
75                 TKDNodeInfo(const TKDNodeInfo &);
76                 TKDNodeInfo& operator=(const TKDNodeInfo &);
77                 
78                 ClassDef(TKDNodeInfo, 1) // node info for interpolator
79         };
80
81 protected:
82         Int_t     fNTNodes;        //Number of evaluation data points
83         TKDNodeInfo *fTNodes;      //[fNTNodes] interpolation node
84
85 private:
86         UChar_t   fStatus;         // status of the interpolator
87         UChar_t   fLambda;         // number of parameters in polynom
88         Short_t         fDepth;          //! depth of the KD Tree structure used
89         Float_t   fAlpha;          // alpha parameter
90         Float_t   **fRefPoints;    //! temporary storage of COG data
91         Double_t        *fBuffer;        //! working space [2*fLambda]
92         TKDTreeIF *fKDhelper;      //! kNN finder
93         TLinearFitter *fFitter;    //! linear fitter    
94
95         ClassDef(TKDInterpolator, 1)   // data interpolator based on KD tree
96 };
97
98 //__________________________________________________________________
99 Bool_t  TKDInterpolator::GetCOGPoint(Int_t node, Float_t *coord, Float_t &val, Float_t &err) const
100 {
101         if(node < 0 || node > fNTNodes) return kFALSE;
102
103         for(int idim=0; idim<fNDim; idim++) coord[idim] = fTNodes[node].fRefPoint[idim];
104         val = fTNodes[node].fRefValue;
105         err = fTNodes[node].fRefValue/TMath::Sqrt(fBucketSize);
106         return kTRUE;
107 }
108
109 //__________________________________________________________________
110 Bool_t  TKDInterpolator::GetDataPoint(Int_t n, Float_t *p) const
111 {
112         if(n < 0 || n >= fNpoints) return kFALSE;
113         if(!fData) return kFALSE;
114                 
115         for(int i=0; i<fNDim; i++) p[i] = fData[i][n];
116         return kTRUE;
117 }
118
119
120 #endif
121