AtmConstraint.cxx
Go to the documentation of this file.
2 
4 
5 #include "Utilities/func/MathUtil.h"
6 
7 #include "TAxis.h"
8 #include "TMatrixD.h"
9 #include "TGraph.h"
10 #include "TH2.h"
11 #include "TMatrixD.h"
12 
13 #include "OscLib/IOscCalc.h"
14 
15 namespace ana
16 {
17  //----------------------------------------------------------------------
19  const std::string& surfname_norm,
20  const std::string& surfname_inv)
21  {
22  TFile f(fname.c_str());
23  assert(!f.IsZombie());
24 
25  fLLNorm = (TH2*)f.Get(surfname_norm.c_str())->Clone();
26  assert(fLLNorm);
27  fLLNorm->SetDirectory(0);
28 
29  fLLInv = (TH2*)f.Get(surfname_inv.c_str())->Clone();
30  assert(fLLInv);
31  fLLInv->SetDirectory(0);
32 
33  fIsInvHiNeg = (fLLInv->GetYaxis()->GetXmin() < 0);
34  fScaleAxes = fabs(fLLNorm->GetYaxis()->GetBinCenter(1)) > 1;
35  }
36 
37  //----------------------------------------------------------------------
39  const std::string& surfname_norm,
40  const std::string& fname_inv,
41  const std::string& surfname_inv)
42  {
43  TFile fNorm(fname_norm.c_str());
44  assert(!fNorm.IsZombie());
45  TFile fInv(fname_inv.c_str());
46  assert(!fInv.IsZombie());
47 
48  fLLNorm = (TH2*)fNorm.Get(surfname_norm.c_str())->Clone();
49  assert(fLLNorm);
50  fLLNorm->SetDirectory(0);
51 
52  fLLInv = (TH2*)fInv.Get(surfname_inv.c_str())->Clone();
53  assert(fLLInv);
54  fLLInv->SetDirectory(0);
55 
56  fIsInvHiNeg = (fLLInv->GetYaxis()->GetXmin() < 0);
57  fScaleAxes = fabs(fLLNorm->GetYaxis()->GetBinCenter(1)) > 1;
58  }
59 
60  //----------------------------------------------------------------------
62  {
63  delete fLLNorm;
64  delete fLLInv;
65  }
66 
67  //----------------------------------------------------------------------
69  const ana::SystShifts& syst) const
70  {
71  double Th23 = osc->GetTh23();
72  double SsqTh23 = util::sqr(std::sin(Th23));
73  double Dmsq23 = osc->GetDmsq32();
74 
75  return BiCubicInterp(SsqTh23, Dmsq23);
76  }
77 
78  //----------------------------------------------------------------------
79  double AtmConstraint::BiCubicInterp(double x, double y) const
80  {
81  TH2* ll = (y < 0) ? fLLInv : fLLNorm;
82  y = (fIsInvHiNeg) ? y : fabs(y);
83  if (fScaleAxes) y *= 1000; // NOvA puts Dmsq32 in units of 1e-3eV^2
84 
85  int xLeftIdx = GetGridLowerLeftIdx(x, ll->GetXaxis());
86  int yLowerIdx = GetGridLowerLeftIdx(y, ll->GetYaxis());
87  // If you're out of bounds, return a large number
88  if (xLeftIdx < 1 || yLowerIdx < 1) return 1e3;
89  if (xLeftIdx+3 > ll->GetNbinsX()) return 1e3;
90  if (yLowerIdx+3 > ll->GetNbinsY()) return 1e3;
91 
92  TMatrixD grid(4,4);
93  GetGridMatrix(ll, xLeftIdx, yLowerIdx, grid);
94  TMatrixD leftForm(4,4);
95  GetFormMatrix(leftForm);
96  TMatrixD rightForm(4,4);
97  GetFormMatrix(rightForm);
98  rightForm.T();
99 
100  TMatrixD intermediate(4,4);
101  intermediate.Mult(leftForm, grid);
102  TMatrixD coef(4,4);
103  coef.Mult(intermediate, rightForm);
104 
105  // coordinates need to be mapped to unit square [0,1]x[0,1]
106  double normX = (x-ll->GetXaxis()->GetBinCenter(xLeftIdx+1)) /
107  (ll->GetXaxis()->GetBinWidth(xLeftIdx+1));
108  double normY = (y-ll->GetYaxis()->GetBinCenter(yLowerIdx+1)) /
109  (ll->GetYaxis()->GetBinWidth(yLowerIdx+1));
110 
111  double interp = 0;
112  for (int i = 0; i < 4; i++)
113  for (int j = 0; j < 4; j++)
114  interp += coef(i,j) * util::ipow(normX,i) * util::ipow(normY,j);
115  return interp;
116  }
117 
118  //----------------------------------------------------------------------
120  TAxis* axis) const
121  {
122  return axis->FindBin(coord + 0.5*axis->GetBinWidth(coord)) - 2;
123  }
124 
125  //----------------------------------------------------------------------
127  int leftX, int lowerY,
128  TMatrixD &grid) const
129  {
130  double f[4][4];
131 
132  for (int i = 0; i < 4; i++)
133  for (int j = 0; j < 4; j++)
134  f[i][j] = ll->GetBinContent(leftX+i, lowerY+j);
135 
136  // Block of hist values
137  grid(0,0) = f[1][1];
138  grid(0,1) = f[1][2];
139  grid(1,0) = f[2][1];
140  grid(1,1) = f[2][2];
141 
142  // Block of Y-derivatives
143  grid(0,2) = (f[1][2] - f[1][0]) / 2;
144  grid(0,3) = (f[1][3] - f[1][1]) / 2;
145  grid(1,2) = (f[2][2] - f[2][0]) / 2;
146  grid(1,3) = (f[2][3] - f[2][1]) / 2;
147 
148  // Block of X-derivatives
149  grid(2,0) = (f[2][1] - f[0][1]) / 2;
150  grid(2,1) = (f[2][2] - f[0][2]) / 2;
151  grid(3,0) = (f[3][1] - f[1][1]) / 2;
152  grid(3,1) = (f[3][2] - f[1][2]) / 2;
153 
154  // Block of XY-derivatives
155  grid(2,2) = (f[2][2] + f[0][0] - f[2][0] - f[0][2]) / 4;
156  grid(2,3) = (f[2][3] + f[0][1] - f[2][1] - f[0][3]) / 4;
157  grid(3,2) = (f[3][2] + f[1][0] - f[3][0] - f[1][2]) / 4;
158  grid(3,3) = (f[3][3] + f[1][1] - f[3][1] - f[1][3]) / 4;
159  }
160 
161  //----------------------------------------------------------------------
163  {
164  // From Wikipedia, a magic matrix that lets you do the interp without a
165  // matrix inversion ==> it's way faster to implement
166  form(0,0) = 1;
167  form(0,1) = 0;
168  form(0,2) = 0;
169  form(0,3) = 0;
170 
171  form(1,0) = 0;
172  form(1,1) = 0;
173  form(1,2) = 1;
174  form(1,3) = 0;
175 
176  form(2,0) = -3;
177  form(2,1) = 3;
178  form(2,2) = -2;
179  form(2,3) = -1;
180 
181  form(3,0) = 2;
182  form(3,1) = -2;
183  form(3,2) = 1;
184  form(3,3) = 1;
185  }
186 
187  //----------------------------------------------------------------------
189  {
190  const std::string fname = FindCAFAnaDir() + "/data/expt/T2K-numuDisappearanceData-Run1to4-2014.root";
191  return new AtmConstraint(fname,
192  "t2krun1to4_normal_2D_n2dlnLsurface",
193  "t2krun1to4_inverted_2D_n2dlnLsurface");
194  }
195 
196  //----------------------------------------------------------------------
198  {
199  const std::string fname = FindCAFAnaDir() + "/data/expt/MINOS_ThreeFlavour_Surfaces_PRL_112_191801.root";
200  return new AtmConstraint(fname,
201  "hMINOS_NumuAndNue_dChi2_NH",
202  "hMINOS_NumuAndNue_dChi2_IH");
203  }
204 }
double BiCubicInterp(double x, double y) const
void GetGridMatrix(TH2 *ll, int lowerX, int lowerY, TMatrixD &grid) const
virtual ~AtmConstraint()
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const AtmConstraint * MINOSAtmConstraint()
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
int GetGridLowerLeftIdx(double coord, TAxis *axis) const
virtual T GetTh23() const
Definition: IOscCalc.h:89
std::string FindCAFAnaDir()
Definition: Utilities.cxx:429
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
void GetFormMatrix(TMatrixD &form) const
const AtmConstraint * T2KAtmConstraint()
Constraint from T2K&#39;s Oct 31, 2014 paper, arXiv 1410.1532v1.
const double j
Definition: BetheBloch.cxx:29
T ipow(T x, unsigned int n)
More efficient exponentiation function than pow(x,n) for small n.
Definition: MathUtil.h:29
virtual T GetDmsq32() const
Definition: IOscCalc.h:86
Oscillation probability calculators.
Definition: Calcs.h:5
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T sin(T number)
Definition: d0nt_math.hpp:132
assert(nhit_max >=nhit_nbins)
AtmConstraint(const std::string &fname, const std::string &surfname_norm, const std::string &surfname_inv)