FitVarsProduct.cxx
Go to the documentation of this file.
2 
3 #include "OscLib/IOscCalc.h"
4 
5 #include <cassert>
6 #include <cmath>
7 
8 namespace ana
9 {
10  // see the .h for why these are functions and not just const FitVarsProduct(Marg) objects
12  {
14  return fvp;
15  }
17  {
19  return fvpm;
20  }
22  {
24  return f;
25  }
26 
27  // Helper functions
28  namespace
29  {
30  // This is the definition of the coordinates alpha and beta
31  void XYToAlphaBeta(double x, double y, double& alpha, double& beta)
32  {
33  assert(x >= 0 && x <= 1);
34  assert(y >= 0 && y <= 1);
35 
36  if(x == 1 && y == 1){
37  alpha = 1;
38  beta = 1;
39  return;
40  }
41 
42  alpha = x*y;
43  beta = 2/M_PI*atan2(1-x, 1-y);
44 
45  assert(alpha >= 0 && alpha <= 1);
46  assert(beta >= 0 && beta <= 1);
47  }
48 
49  // And this is the solution back to x and y
50  void AlphaBetaToXY(double alpha, double beta, double& x, double& y)
51  {
52  assert(alpha >= 0 && alpha <= 1);
53  assert(beta >= 0 && beta <= 1);
54 
55  if(beta == 0){
56  x = 1;
57  y = alpha/x;
58  return;
59  }
60  if(beta == 1){
61  y = 1;
62  x = alpha/y;
63  return;
64  }
65 
66  const double phi = tan(beta*M_PI/2)-1;
67  const double det = sqrt(4*alpha*(phi+1)+phi*phi);
68 
69  // These are the full expressions. But for very extreme values of beta
70  // (near zero or one) there can be rounding problems. So evaluate the
71  // easy expression and then ensure that alpha truly remains fixed.
72 
73  // x = (det-phi)/2;
74  // y = (det+phi)/(2+2*phi);
75 
76  if(alpha == 0){
77  x = (det-phi)/2;
78  y = (det+phi)/(2+2*phi);
79  }
80  else{
81  if(beta < .5){
82  x = (det-phi)/2;
83  y = alpha/x;
84  }
85  else{
86  y = (det+phi)/(2+2*phi);
87  x = alpha/y;
88  }
89  }
90 
91  assert(x >= 0 && x <= 1+1e-6);
92  assert(y >= 0 && y <= 1+1e-6);
93 
94  x = std::min(x, 1.);
95  y = std::min(y, 1.);
96  }
97  }
98 
99  //----------------------------------------------------------------------
101  const IConstrainedFitVar& y)
102  : IConstrainedFitVar(x.ShortName()+"*"+y.ShortName(),
103  x.LatexName()+y.LatexName()),
104  fVarX(x), fVarY(y)
105  {
106  assert(x.LowLimit() == 0 && x.HighLimit() == 1);
107  assert(y.LowLimit() == 0 && y.HighLimit() == 1);
108  }
109 
110  //----------------------------------------------------------------------
112  {
113  // Easy, just the product of the two
114  double alpha, beta;
115  XYToAlphaBeta(fVarX.GetValue(osc), fVarY.GetValue(osc), alpha, beta);
116  return alpha;
117  }
118 
119  //----------------------------------------------------------------------
121  {
122  // Conver to the other space, hold beta (the marg variable) fixed and
123  // adjust alpha, see where that leaves us.
124  double alpha, beta;
125  XYToAlphaBeta(fVarX.GetValue(osc), fVarY.GetValue(osc), alpha, beta);
126  alpha = Clamp(val);
127 
128  double x, y;
129  AlphaBetaToXY(alpha, beta, x, y);
130  fVarX.SetValue(osc, x);
131  fVarY.SetValue(osc, y);
132  }
133 
134  //----------------------------------------------------------------------
136  const IConstrainedFitVar& y)
137  : IConstrainedFitVar(x.ShortName()+"*"+y.ShortName()+"_marg",
138  "Marginalization parameter for "+x.LatexName()+y.LatexName()),
139  fVarX(x), fVarY(y)
140  {
141  assert(x.LowLimit() == 0 && x.HighLimit() == 1);
142  assert(y.LowLimit() == 0 && y.HighLimit() == 1);
143  }
144 
145  //----------------------------------------------------------------------
147  {
148  // Just evaluate the formula for beta.
149  double alpha, beta;
150  XYToAlphaBeta(fVarX.GetValue(osc), fVarY.GetValue(osc), alpha, beta);
151  return beta;
152  }
153 
154  //----------------------------------------------------------------------
156  {
157  // Switch to the other space, hold alpha constant and adjust beta, see
158  // where that leaves us.
159  double alpha, beta;
160  XYToAlphaBeta(fVarX.GetValue(osc), fVarY.GetValue(osc), alpha, beta);
161  beta = Clamp(val);
162 
163  double x, y;
164  AlphaBetaToXY(alpha, beta, x, y);
165  fVarX.SetValue(osc, x);
166  fVarY.SetValue(osc, y);
167  }
168 } // namespace
virtual double HighLimit() const =0
FitVarsProductMarg(const IConstrainedFitVar &x, const IConstrainedFitVar &y)
const std::string & LatexName() const
Definition: IFitVar.h:38
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
const IConstrainedFitVar & fVarY
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual double LowLimit() const =0
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
Double_t beta
#define M_PI
Definition: SbMath.h:34
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
FitVarsProductMarg FitSinSqTheta23SinSq2Theta13Marg()
Fit the product of two variables (made from angles so their ranges are 0-1). See also FitVarsProductM...
FitVarsProduct FitVarSinSqTheta23SinSq2Theta13()
Return a FitVarsProduct for .
T Clamp(T val) const
Definition: IFitVar.h:61
const IConstrainedFitVar & fVarY
Oscillation probability calculators.
Definition: Calcs.h:5
FitVarsProduct(const IConstrainedFitVar &x, const IConstrainedFitVar &y)
const std::string & ShortName() const
Definition: IFitVar.h:37
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const =0
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
See the documentation for FitVarsProductMarg.
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
Fit2SinSqTheta23SinSq2Theta13 FitVar2SinSqTheta23SinSq2Theta13()
Return a Fit2SinSqTheta23SinSq2Theta13 for .
T min(const caf::Proxy< T > &a, T b)
T tan(T number)
Definition: d0nt_math.hpp:144
Float_t e
Definition: plot.C:35
const IConstrainedFitVar & fVarX
const IConstrainedFitVar & fVarX
T atan2(T number)
Definition: d0nt_math.hpp:72
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const =0