SurfaceKrige.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/Kriger.h"
4 #include "CAFAna/Core/IFitVar.h"
5 #include "CAFAna/Core/Progress.h"
7 
8 #include "Utilities/func/MathUtil.h"
9 
10 #include "TH2.h"
11 
12 namespace ana
13 {
14  //----------------------------------------------------------------------
16  : fGrad(1), fGradX(1), fGradY(1)
17  {
18  }
19 
20  //----------------------------------------------------------------------
21  double KrigeKernel::operator()(double x1, double y1, double x2, double y2) const
22  {
23  // TODO - the 2D model doesn't seem to actually work
24  if(true){//fPts.size() < 10){
25  // Use a very simple function when we don't have much data
26  return fGrad*util::pythag(x1-x2, y1-y2);
27  }
28  else{
29  // Treat x and y separately later
30  return fGradX*fabs(x1-x2) + fGradY*fabs(y1-y2);
31  }
32  }
33 
34  //----------------------------------------------------------------------
35  void KrigeKernel::AddPoint(double x, double y, double z)
36  {
37  // Learn the best function of dx and dy to predict the variance
38 
39  Pt p = {x, y, z};
40  fPts.push_back(p);
41 
42  double R2 = 0, ZR = 0;
43  double X2 = 0, XY = 0, Y2 = 0, ZX = 0, ZY = 0;
44  for(unsigned int ai = 0; ai < fPts.size(); ++ai){
45  for(unsigned int bi = ai+1; bi < fPts.size(); ++bi){
46  const Pt a = fPts[ai];
47  const Pt b = fPts[bi];
48 
49  const double dz2 = util::sqr(a.z-b.z);
50 
51  const double dx = fabs(a.x-b.x);
52  const double dy = fabs(a.y-b.y);
53 
54  const double dr = util::pythag(a.x-b.x, a.y-b.y);
55 
56  R2 += dr*dr;
57  ZR += dz2*dr;
58 
59  X2 += dx*dx;
60  XY += dx*dy;
61  Y2 += dy*dy;
62  ZX += dz2 * dx;
63  ZY += dz2 * dy;
64  }
65  }
66 
67  if(R2) fGrad = ZR / R2;
68 
69  // Only compute 2D coefficients once we have some stats
70  if(fPts.size() < 10) return;
71 
72  TMatrixD M(2, 2); TVectorD v(2);
73  M(0, 0) = X2; M(0, 1) = XY; v(0) = ZX;
74  M(1, 0) = XY; M(1, 1) = Y2; v(1) = ZY;
75 
76  M.Invert();
77  TVectorD g = M*v;
78 
79  fGradX = fabs(g(0));
80  fGradY = fabs(g(1));
81  }
82 
83  //----------------------------------------------------------------------
87  const IFitVar* xvar, int nbinsx, double xmin, double xmax,
88  const IFitVar* yvar, int nbinsy, double ymin, double ymax,
89  const std::vector<const IFitVar*>& profVars,
90  const std::vector<const ISyst*>& profSysts,
91  const SeedList& seedPts,
92  const std::vector<SystShifts>& systSeedPts,
94  {
95  fLogX = fLogY = false;
96 
97  fParallel = false;
98  fFitOpts = opts;
99 
100  CreateHistograms(FitAxis(xvar, nbinsx, xmin, xmax, fLogX),
101  FitAxis(yvar, nbinsy, ymin, ymax, fLogY),
102  profVars, profSysts);
103 
104  fErr = ExpandedHistogram(";"+xvar->LatexName()+";"+yvar->LatexName(),
105  nbinsx, xmin, xmax, fLogX,
106  nbinsy, ymin, ymax, fLogY);
107 
108  for(const IFitVar* v: profVars) fSeedValues.push_back(v->GetValue( calc));
109 
110  FillSurface(expt, calc, xvar, yvar, profVars, profSysts, seedPts, systSeedPts);
111 
112  FindMinimum(expt, calc, xvar, yvar, profVars, profSysts, seedPts, systSeedPts);
113  }
114 
115  //---------------------------------------------------------------------
117  {
118  delete fErr;
119  }
120 
121  //---------------------------------------------------------------------
122  void SurfaceKrige::
125  const IFitVar* xvar, const IFitVar* yvar,
126  const std::vector<const IFitVar*>& profVars,
127  const std::vector<const ISyst*>& profSysts,
128  const SeedList& seedPts,
129  const std::vector<SystShifts>& systSeedPts)
130  {
131  KrigeKernel kern;
132  // Need this layer of indirection to capture the kernel by reference
133  Kriger kriger([&kern](double x1, double y1, double x2, double y2)
134  {
135  return kern(x1, y1, x2, y2);
136  });
137 
138  // Start off with one point in the centre
139  int targetx = fHist->GetNbinsX()/2;
140  int targety = fHist->GetNbinsY()/2;
141 
142  Progress prog(ProgressBarTitle(xvar, yvar, profVars, profSysts));
143 
144  // TODO: Arbitrary number of passes, need real stopping condition...
145  const int maxPass = 100;//fHist->GetNbinsX()*fHist->GetNbinsY()/10;
146 
147  for(int pass = 0; pass < maxPass; ++pass){
148 
149  // Evaluate the function at the target point and add the info into the
150  // Kriger
151  const double targetz = FillSurfacePoint(expt, calc,
152  xvar, BinCenterX(targetx),
153  yvar, BinCenterY(targety),
154  profVars, profSysts,
155  seedPts, systSeedPts);
156 
157  kriger.AddPoint(targetx, targety, targetz);
158 
159  kern.AddPoint(targetx, targety, targetz);
160 
161  prog.SetProgress(double(pass)/maxPass);
162 
163  // Search through the surface for the next most uncertain point
164  double maxerr = -1;
165 
166  for(int ix = 1; ix <= fHist->GetNbinsX(); ++ix){
167  for(int iy = 1; iy <= fHist->GetNbinsY(); ++iy){
168 
169  double err;
170  const double z = kriger.Predict(ix, iy, &err);
171 
172  // Keep surface and error surface up-to-date as we go
173  fHist->SetBinContent(ix, iy, z);
174  fErr->SetBinContent(ix, iy, err);
175 
176  // New candidate
177  if(err > maxerr){
178  maxerr = err;
179  targetx = ix;
180  targety = iy;
181  }
182  } // end for iy
183  } // end for ix
184 
185  } // end for pass
186  }
187 }
const std::string & LatexName() const
Definition: IFitVar.h:37
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
std::vector< Pt > fPts
Definition: SurfaceKrige.h:24
Float_t x1[n_points_granero]
Definition: compare.C:5
const char * p
Definition: xmltok.h:285
void AddPoint(double x, double y, double z)
double operator()(double x1, double y1, double x2, double y2) const
osc::OscCalcDumb calc
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Double_t X2
Definition: plot.C:264
TH2F * ExpandedHistogram(const std::string &title, int nbinsx, double xmin, double xmax, bool xlog, int nbinsy, double ymin, double ymax, bool ylog)
Internal helper for Surface and FCSurface.
Definition: Utilities.cxx:154
Double_t ymax
Definition: plot.C:25
expt
Definition: demo5.py:34
double dy[NP][NC]
double dx[NP][NC]
void FillSurface(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *xvar, const IFitVar *yvar, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts) override
const double a
https://en.wikipedia.org/wiki/Kriging
Definition: Kriger.h:12
z
Definition: test.py:28
void FindMinimum(TH1 *&h, const int minx, const int maxx)
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
Int_t nbinsx
Definition: plot.C:23
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Base class defining interface for experiments.
Definition: IExperiment.h:14
const hit & b
Definition: hits.cxx:21
A simple ascii-art progress bar.
Definition: Progress.h:9
Interface definition for fittable variables.
Definition: IFitVar.h:16
Double_t ymin
Definition: plot.C:24
SurfaceKrige(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *xvar, int nbinsx, double xmin, double xmax, const IFitVar *yvar, int nbinsy, double ymin, double ymax, const std::vector< const IFitVar * > &profVars={}, const std::vector< const ISyst * > &profSysts={}, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, MinuitFitter::FitOpts opts=MinuitFitter::kNormal)
Double_t Y2
Definition: plot.C:264
Collect information describing the axis of a fit variable.
Definition: FitAxis.h:10
Helper for SurfaceKrige.
Definition: SurfaceKrige.h:10