Kriger.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Kriger.h"
2 
3 namespace ana
4 {
5  //----------------------------------------------------------------------
6  void Kriger::AddPoint(double x, double y, double z)
7  {
8  fPts.emplace_back(x, y, z);
9  }
10 
11  //----------------------------------------------------------------------
12  double Kriger::Predict(double x, double y, double* err)
13  {
14  double junk;
15  if(!err) err = &junk;
16 
17  EnsureDecomp();
18 
19  Pt p(x, y, 0);
20 
21  TVectorD b(N()+1);
22  for(unsigned int i = 0; i < N(); ++i){
23  b(i) = fKern(p.x, p.y, fPts[i].x, fPts[i].y);
24  }
25  b(N()) = 1;
26 
27  bool ok;
28  TVectorD lambda = fDecomp.Solve(b, ok);
29 
30  double z = 0;
31  *err = 0;
32  for(unsigned int i = 0; i < N(); ++i){
33  z += lambda[i] * fPts[i].z;
34  *err += lambda[i] * b[i];
35  }
36  *err = sqrt(fabs(*err));
37 
38  return z;
39  }
40 
41  //----------------------------------------------------------------------
43  {
44  if(unsigned(fMat.GetNrows()) == N()+1) return;
45  fMat.ResizeTo(N()+1, N()+1);
46 
47  for(unsigned int i = 0; i < N(); ++i){
48  fMat(i, i) = 0;
49  for(unsigned int j = i+1; j < N(); ++j){
50  fMat(i, j) = fMat(j, i) = fKern(fPts[i].x, fPts[i].y, fPts[j].x, fPts[j].y);
51  }
52  }
53 
54  // Sum of coefficients is 1
55  for(unsigned int i = 0; i < N(); ++i) fMat(N(), i) = 1;
56  // Include the lagrange multiplier in each
57  for(unsigned int i = 0; i < N(); ++i) fMat(i, N()) = 1;
58  }
59 
60  //----------------------------------------------------------------------
62  {
63  EnsureMatrix();
64 
65  if(fDecomp.GetNrows() == fMat.GetNrows()) return;
66  fDecomp.SetMatrix(fMat);
67  }
68 }
void EnsureDecomp()
Definition: Kriger.cxx:61
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const char * p
Definition: xmltok.h:285
void AddPoint(double x, double y, double z)
Definition: Kriger.cxx:6
T sqrt(T number)
Definition: d0nt_math.hpp:156
double y
Definition: Kriger.h:35
TMatrixD fMat
Definition: Kriger.h:39
TDecompLU fDecomp
Definition: Kriger.h:40
double Predict(double x, double y, double *err=0)
Definition: Kriger.cxx:12
void EnsureMatrix()
Definition: Kriger.cxx:42
const double j
Definition: BetheBloch.cxx:29
z
Definition: test.py:28
double x
Definition: Kriger.h:35
Kernel_t fKern
Definition: Kriger.h:41
std::vector< Pt > fPts
Definition: Kriger.h:38
unsigned int N() const
Definition: Kriger.h:24
const hit & b
Definition: hits.cxx:21