Var.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Var.h"
2 
3 #include "CAFAna/Core/DepMan.h"
4 
5 #include <algorithm>
6 #include <cmath>
7 #include <iostream>
8 #include <map>
9 
10 namespace ana
11 {
12  //----------------------------------------------------------------------
13  template<class T> _Var<T>::
14  _Var(const std::function<VarFunc_t>& fun)
15  : fFunc(fun), fID(fgNextID++)
16  {
17  DepMan<_Var<T>>::Instance().RegisterConstruction(this);
18  }
19 
20  //----------------------------------------------------------------------
21  template<class T> _Var<T>::
22  _Var(const std::function<VarFunc_t>& fun, int id)
23  : fFunc(fun), fID(id)
24  {
25  DepMan<_Var<T>>::Instance().RegisterConstruction(this);
26  }
27 
28  //----------------------------------------------------------------------
29  template<class T> _Var<T>::~_Var()
30  {
31  DepMan<_Var<T>>::Instance().RegisterDestruction(this);
32  }
33 
34  //----------------------------------------------------------------------
35  template<class T> _Var<T>::_Var(const _Var<T>& v)
36  : fFunc(0), fID(-1)
37  {
38  if(&v == this) return;
39 
40  if(v.fFunc){
41  fFunc = v.fFunc;
42  fID = v.fID;
43 
44  DepMan<_Var<T>>::Instance().RegisterConstruction(this);
45  }
46  else{
47  // If we are copying from a Var with NULL func, that is probably because
48  // it is all zero because it hasn't been statically constructed
49  // yet. Register our interest of getting constructed in turn once it is.
50  DepMan<_Var<T>>::Instance().RegisterDependency(&v, this);
51  }
52  }
53 
54  //----------------------------------------------------------------------
55  template<class T> _Var<T>& _Var<T>::
57  {
58  if(&v == this) return *this;
59 
60  if(v.fFunc){
61  fFunc = v.fFunc;
62  fID = v.fID;
63 
64  DepMan<_Var<T>>::Instance().RegisterConstruction(this);
65  }
66  else{
67  fFunc = 0;
68  fID = -1;
69 
70  // If we are copying from a Var with NULL func, that is probably because
71  // it is all zero because it hasn't been statically constructed
72  // yet. Register our interest of getting constructed in turn once it is.
73  DepMan<_Var<T>>::Instance().RegisterDependency(&v, this);
74  }
75 
76  return *this;
77  }
78 
79 
80  //----------------------------------------------------------------------
81  // Extract the value (trivial) and print a good error message if it's NaN
82  template<class T>
83  double ValHelper(const _Var<T>& v,
84  const std::string& op,
85  double c, const T* sr)
86  {
87  const double val = v(sr);
88  if(std::isnan(val)){
89  std::cout << "Warning: Cut compares NaN " << op << " " << c << std::endl;
90  }
91  return val;
92  }
93 
94  template<class T>
95  void ValHelper(const _Var<T>& va, const _Var<T>& vb,
96  const std::string& op,
97  const T* sr,
98  double& a, double& b)
99  {
100  a = va(sr);
101  b = vb(sr);
102  if(std::isnan(a) || std::isnan(b)){
103  std::cout << "Warning: Cut compares " << a << " " << op << " " << b << std::endl;
104  }
105  }
106 
107 
108  // These comparison operators are extremely rote. Use a macro
109 
110 #define COMPARISON(OP, OPNAME) \
111  /* Comparison of a Var with a constant */ \
112  template<class T> _Cut<T> _Var<T>:: \
113  operator OP(double c) const \
114  { \
115  struct OPNAME \
116  { \
117  _Var<T> v; double c; \
118  bool operator()(const T* sr) \
119  { \
120  return ValHelper(v, #OP, c, sr) OP c; \
121  } \
122  }; \
123  \
124  return _Cut<T>(OPNAME{*this, c}); \
125  } \
126  \
127  /* Comparison of a Var with another Var */ \
128  template<class T> _Cut<T> _Var<T>:: \
129  operator OP(const _Var<T>& v) const \
130  { \
131  struct OPNAME \
132  { \
133  _Var<T> a, b; \
134  bool operator()(const T* sr) \
135  { \
136  double va, vb; \
137  ValHelper(a, b, #OP, sr, va, vb); \
138  return va OP vb; \
139  } \
140  }; \
141  \
142  return _Cut<T>(OPNAME{*this, v}); \
143  } \
144  void dummy() // trick to require a semicolon
145 
146  COMPARISON(>, Greater);
147  COMPARISON(>=, GreaterEquals);
148  COMPARISON(<, Less);
149  COMPARISON(<=, LessEquals);
150  COMPARISON(==, Equals);
151  COMPARISON(!=, NotEquals);
152 
153 
154  // explicitly instantiate the template for the types we know we have
155  template class _Var<caf::SRProxy>;
156  template class _Var<caf::SRSpillProxy>;
157  template class _Var<caf::SRNeutrinoProxy>;
158 
159  template<class T> int _Var<T>::fgNextID = 0;
160 
161  //----------------------------------------------------------------------
162  /// Helper for \ref Var2D
163  template<class T> class Var2DFunc
164  {
165  public:
166  Var2DFunc(const _Var<T>& a, const Binning binsa,
167  const _Var<T>& b, const Binning binsb)
168  : fA(a), fBinsA(binsa),
169  fB(b), fBinsB(binsb)
170  {
171  }
172 
173  double operator()(const T* sr) const
174  {
175  // Calculate current values of the variables in StandardRecord once
176  const double va = fA(sr);
177  const double vb = fB(sr);
178 
179  // Since there are no overflow/underflow bins, check the range
180  if(va < fBinsA.Min() || vb < fBinsB.Min()) return -1;
181  if(va > fBinsA.Max() || vb > fBinsB.Max()) return fBinsA.NBins() * fBinsB.NBins();
182 
183  // FindBin uses root convention, first bin is bin 1, bin 0 is underflow
184  const int ia = fBinsA.FindBin(va) - 1;
185  const int ib = fBinsB.FindBin(vb) - 1;
186 
187  const int i = ia*fBinsB.NBins()+ib;
188 
189  return i+.5;
190  }
191 
192  protected:
193  const _Var<T> fA;
195  const _Var<T> fB;
197  };
198 
199  /// Helper for \ref Var3D
200  template<class T> class Var3DFunc
201  {
202  public:
203  Var3DFunc(const _Var<T>& a, const Binning binsa,
204  const _Var<T>& b, const Binning binsb,
205  const _Var<T>& c, const Binning binsc)
206  : fA(a), fBinsA(binsa),
207  fB(b), fBinsB(binsb),
208  fC(c), fBinsC(binsc)
209  {
210  }
211 
212  double operator()(const T* sr) const
213  {
214  /// Calculate current values of the variables in StandardRecord once
215  const double va = fA(sr);
216  const double vb = fB(sr);
217  const double vc = fC(sr);
218 
219  /// Since there are no overflow/underflow bins, check the range
220  if(va < fBinsA.Min() || vb < fBinsB.Min() || vc < fBinsC.Min()){
221  return -1.0;
222  }
223 
224  if(va > fBinsA.Max() || vb > fBinsB.Max() || vc > fBinsC.Max()){
225  return fBinsA.NBins() * fBinsB.NBins() * fBinsC.NBins();
226  }
227 
228  const int ia = fBinsA.FindBin(va) - 1;
229  const int ib = fBinsB.FindBin(vb) - 1;
230  const int ic = fBinsC.FindBin(vc) - 1;
231  const int i = ia*fBinsB.NBins()*fBinsC.NBins() + ib*fBinsC.NBins() + ic;
232 
233  return i+.5;
234  }
235 
236  protected:
237  const _Var<T> fA;
239  const _Var<T> fB;
241  const _Var<T> fC;
243  };
244 
245  //----------------------------------------------------------------------
246  template<class T> _Var<T>
247  Var2D(const _Var<T>& a, const Binning& binsa,
248  const _Var<T>& b, const Binning& binsb)
249  {
250  return _Var<T>(Var2DFunc<T>(a, binsa, b, binsb));
251  }
252 
253  //----------------------------------------------------------------------
254  template<class T> _Var<T>
255  Var2D(const _Var<T>& a, int na, double a0, double a1,
256  const _Var<T>& b, int nb, double b0, double b1)
257  {
258  return Var2D(a, Binning::Simple(na, a0, a1),
259  b, Binning::Simple(nb, b0, b1));
260  }
261 
262  // explicitly instantiate the template for the types we know we have
263  template Var Var2D(const Var&, const Binning&, const Var&, const Binning&);
264  template SpillVar Var2D(const SpillVar&, const Binning&, const SpillVar&, const Binning&);
265  template NuTruthVar Var2D(const NuTruthVar&, const Binning&, const NuTruthVar&, const Binning&);
266 
267  template Var Var2D(const Var&, int, double, double, const Var&, int, double, double);
268  template SpillVar Var2D(const SpillVar&, int, double, double, const SpillVar&, int, double, double);
269  template NuTruthVar Var2D(const NuTruthVar&, int, double, double, const NuTruthVar&, int, double, double);
270 
271  //----------------------------------------------------------------------
272  template<class T> _Var<T>
273  Var3D(const _Var<T>& a, const Binning& binsa,
274  const _Var<T>& b, const Binning& binsb,
275  const _Var<T>& c, const Binning& binsc)
276  {
277  return _Var<T>(Var3DFunc<T>(a, binsa, b, binsb, c, binsc));
278  }
279 
280  //----------------------------------------------------------------------
281  template<class T> _Var<T>
282  Var3D(const _Var<T>& a, int na, double a0, double a1,
283  const _Var<T>& b, int nb, double b0, double b1,
284  const _Var<T>& c, int nc, double c0, double c1)
285  {
286  return Var3D(a, Binning::Simple(na, a0, a1),
287  b, Binning::Simple(nb, b0, b1),
288  c, Binning::Simple(nc, c0, c1));
289  }
290 
291  // explicitly instantiate the template for the types we know we have
292  template Var Var3D(const Var&, const Binning&, const Var&, const Binning&, const Var&, const Binning&);
293  template SpillVar Var3D(const SpillVar&, const Binning&, const SpillVar&, const Binning&, const SpillVar&, const Binning&);
294  template NuTruthVar Var3D(const NuTruthVar&, const Binning&, const NuTruthVar&, const Binning&, const NuTruthVar&, const Binning&);
295 
296  template Var Var3D(const Var&, int, double, double, const Var&, int, double, double, const Var&, int, double, double);
297  template SpillVar Var3D(const SpillVar&, int, double, double, const SpillVar&, int, double, double, const SpillVar&, int, double, double);
298  template NuTruthVar Var3D(const NuTruthVar&, int, double, double, const NuTruthVar&, int, double, double, const NuTruthVar&, int, double, double);
299 
300  //----------------------------------------------------------------------
301  Var Scaled(const Var& v, double s)
302  {
303  struct Scale
304  {
305  const Var v; double s;
306  double operator()(const caf::SRProxy* sr) const {return s*v(sr);}
307  };
308 
309  return Var(Scale{v, s});
310  }
311 
312  //----------------------------------------------------------------------
313  Var Constant(double c)
314  {
315  struct Const
316  {
317  double c;
318  double operator()(const caf::SRProxy*) const {return c;}
319  };
320 
321  return Var(Const{c});
322  }
323 
324  //--------------------------------------------------------------------
325 
326  Var Sqrt(const Var& v)
327  {
328  struct Sqrter
329  {
330  const Var v;
331  double operator()(const caf::SRProxy* sr) const {return sqrt(v(sr));}
332  };
333 
334  return Var(Sqrter{v});
335  }
336 
337  //----------------------------------------------------------------------
338  template<class T> _Var<T> _Var<T>::
339  operator*(const _Var<T>& v) const
340  {
341  static std::map<std::pair<int, int>, int> ids;
342  const std::pair<int, int> key(ID(), v.ID());
343 
344  struct Times
345  {
346  const _Var<T> a, b;
347  double operator()(const T* sr) const {return a(sr) * b(sr);}
348  };
349 
350  if(ids.count(key) == 0) ids[key] = fgNextID++;
351  return _Var<T>(Times{*this, v});
352  }
353 
354  //----------------------------------------------------------------------
355  template<class T> _Var<T> _Var<T>::
356  operator/(const _Var<T>& v) const
357  {
358  static std::map<std::pair<int, int>, int> ids;
359  const std::pair<int, int> key(ID(), v.ID());
360 
361  struct Divide
362  {
363  const _Var<T> a, b;
364  double operator()(const T* sr) const
365  {
366  const double denom = b(sr);
367  if(denom != 0)
368  return a(sr) / denom;
369  else
370  return 0.0;
371  }
372  };
373 
374  if(ids.count(key) == 0) ids[key] = fgNextID++;
375  return _Var<T>(Divide{*this, v}, ids[key]);
376  }
377 
378  //----------------------------------------------------------------------
379  template<class T> _Var<T> _Var<T>::
380  operator+(const _Var<T>& v) const
381  {
382  static std::map<std::pair<int, int>, int> ids;
383  const std::pair<int, int> key(ID(), v.ID());
384 
385  struct Plus
386  {
387  const _Var<T> a, b;
388  double operator()(const T* sr) const {return a(sr) + b(sr);}
389  };
390 
391  if(ids.count(key) == 0) ids[key] = fgNextID++;
392  return _Var<T>(Plus{*this, v}, ids[key]);
393  }
394 
395  //----------------------------------------------------------------------
396  template<class T> _Var<T> _Var<T>::
397  operator-(const _Var<T>& v) const
398  {
399  static std::map<std::pair<int, int>, int> ids;
400  const std::pair<int, int> key(ID(), v.ID());
401 
402  struct Minus
403  {
404  const _Var<T> a, b;
405  double operator()(const T* sr) const {return a(sr) - b(sr);}
406  };
407 
408  if(ids.count(key) == 0) ids[key] = fgNextID++;
409  return _Var<T>(Minus{*this, v}, ids[key]);
410  }
411 } // namespace
std::function< VarFunc_t > fFunc
Definition: Var.h:62
ratio_hxv Divide(hxv, goal_hxv)
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
int ID() const
Vars with the same definition will have the same ID.
Definition: Var.h:36
double operator()(const T *sr) const
Definition: Var.cxx:212
const _Var< T > fB
Definition: Var.cxx:239
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
T sqrt(T number)
Definition: d0nt_math.hpp:156
_Var< T > Var2D(const _Var< T > &a, const Binning &binsa, const _Var< T > &b, const Binning &binsb)
Variable formed from two input variables.
Definition: Var.cxx:247
_Var< T > operator+(const _Var< T > &v) const
Definition: Var.cxx:380
Deep magic to fix static initialization order. Here Be Dragons!
Definition: DepMan.h:12
double operator()(const T *sr) const
Definition: Var.cxx:173
const Binning fBinsA
Definition: Var.cxx:238
_Var< T > operator/(const _Var< T > &v) const
Definition: Var.cxx:356
~_Var()
Definition: Var.cxx:29
static constexpr Double_t fC
Definition: Munits.h:236
_Var< T > Var3D(const _Var< T > &a, const Binning &binsa, const _Var< T > &b, const Binning &binsb, const _Var< T > &c, const Binning &binsc)
This is just like a Var2D, but useful for 3D Spectra.
Definition: Var.cxx:273
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
Var Constant(double c)
Use to weight events up and down by some factor.
Definition: Var.cxx:313
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const _Var< T > fC
Definition: Var.cxx:241
COMPARISON(>, Greater)
const XML_Char * s
Definition: expat.h:262
const Binning fBinsB
Definition: Var.cxx:240
Helper for Var3D.
Definition: Var.cxx:200
TH1F * a1
Definition: f2_nu.C:476
const double a
double operator()(const T *sr) const
Allows a variable to be called with double value = myVar(sr) syntax.
Definition: Var.h:30
caf::StandardRecord * sr
void RegisterDependency(const T *parent, T *child)
Call from copy constructor and assignment operator.
Definition: DepMan.cxx:122
const _Var< T > fA
Definition: Var.cxx:193
OStream cout
Definition: OStream.cxx:6
const _Var< T > fB
Definition: Var.cxx:195
double ValHelper(const _Var< T > &v, const std::string &op, double c, const T *sr)
Definition: Var.cxx:83
int fID
Definition: Var.h:64
_Var< T > operator-(const _Var< T > &v) const
Definition: Var.cxx:397
const Binning fBinsA
Definition: Var.cxx:194
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
void RegisterDestruction(T *x)
Call from destructor.
Definition: DepMan.cxx:97
_Var< T > operator*(const _Var< T > &v) const
Definition: Var.cxx:339
const _Var< T > fA
Definition: Var.cxx:237
const hit & b
Definition: hits.cxx:21
enum BeamMode nc
c1
Definition: demo5.py:24
_Var(const std::function< VarFunc_t > &fun)
std::function can wrap a real function, function object, or lambda
Definition: Var.cxx:14
void RegisterConstruction(T *x)
Call from constructor (in the success case)
Definition: DepMan.cxx:73
double T
Definition: Xdiff_gwt.C:5
Var2DFunc(const _Var< T > &a, const Binning binsa, const _Var< T > &b, const Binning binsb)
Definition: Var.cxx:166
Var Scaled(const Var &v, double s)
Use to rescale another variable.
Definition: Var.cxx:301
const Binning fBinsB
Definition: Var.cxx:196
static const double nb
Definition: Units.h:89
Template for Var and SpillVar.
static int fgNextID
The next ID that hasn&#39;t yet been assigned.
Definition: Var.h:66
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
simulatedPeak Scale(1/simulationNormalisationFactor)
Helper for Var2D.
Definition: Var.cxx:163
_Var & operator=(const _Var< T > &v)
Definition: Var.cxx:56
Var3DFunc(const _Var< T > &a, const Binning binsa, const _Var< T > &b, const Binning binsb, const _Var< T > &c, const Binning binsc)
Definition: Var.cxx:203
const Binning fBinsC
Definition: Var.cxx:242
enum BeamMode string