Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
ana::AtmConstraint Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-11-28/CAFAna/Experiment/AtmConstraint.h"

Inheritance diagram for ana::AtmConstraint:
ana::IExperiment

Public Member Functions

 AtmConstraint (const std::string &fname, const std::string &surfname_norm, const std::string &surfname_inv)
 
 AtmConstraint (const std::string &fname_norm, const std::string &surfname_norm, const std::string &fname_inv, const std::string &surfname_inv)
 
virtual ~AtmConstraint ()
 
virtual double ChiSq (osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
 
virtual stan::math::var LogLikelihood (osc::IOscCalcAdjustableStan *osc, const SystShifts &syst=kNoShift) const
 
virtual void SaveTo (TDirectory *dir, const std::string &name) const
 

Protected Member Functions

double BiCubicInterp (double x, double y) const
 
double CalculateInterp (double normX, double normY, TMatrixD coefs) const
 
int GetGridLowerLeftIdx (double coord, TAxis *axis) const
 
void GetGridMatrix (TH2 *ll, int lowerX, int lowerY, TMatrixD &grid) const
 
void GetFormMatrix (TMatrixD &form) const
 

Protected Attributes

TH2 * fLLNorm
 
TH2 * fLLInv
 
bool fIsInvHiNeg
 
bool fScaleAxes
 

Detailed Description

Constraint on the atmospheric oscillation parameters $ \sin^2\theta_{23} $ and $ \Delta m^2_{32} $ from an external likelihood surface.

Definition at line 15 of file AtmConstraint.h.

Constructor & Destructor Documentation

ana::AtmConstraint::AtmConstraint ( const std::string fname,
const std::string surfname_norm,
const std::string surfname_inv 
)

Definition at line 18 of file AtmConstraint.cxx.

References ana::assert(), plot_validation_datamc::Clone(), MakeMiniprodValidationCuts::f, stan::math::fabs(), fIsInvHiNeg, fLLInv, fLLNorm, and fScaleAxes.

Referenced by ana::MINOSAtmConstraint(), and ana::T2KAtmConstraint().

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  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
assert(nhit_max >=nhit_nbins)
ana::AtmConstraint::AtmConstraint ( const std::string fname_norm,
const std::string surfname_norm,
const std::string fname_inv,
const std::string surfname_inv 
)

Definition at line 38 of file AtmConstraint.cxx.

References ana::assert(), plot_validation_datamc::Clone(), stan::math::fabs(), fIsInvHiNeg, fLLInv, fLLNorm, and fScaleAxes.

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  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
assert(nhit_max >=nhit_nbins)
ana::AtmConstraint::~AtmConstraint ( )
virtual

Definition at line 61 of file AtmConstraint.cxx.

References fLLInv, and fLLNorm.

62  {
63  delete fLLNorm;
64  delete fLLInv;
65  }

Member Function Documentation

double ana::AtmConstraint::BiCubicInterp ( double  x,
double  y 
) const
protected

Definition at line 79 of file AtmConstraint.cxx.

References stan::math::fabs(), fIsInvHiNeg, fLLInv, fLLNorm, fScaleAxes, GetFormMatrix(), GetGridLowerLeftIdx(), GetGridMatrix(), galleryMaker::grid, MECModelEnuComparisons::i, util::ipow(), and calib::j.

Referenced by ChiSq().

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  }
void GetGridMatrix(TH2 *ll, int lowerX, int lowerY, TMatrixD &grid) const
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
int GetGridLowerLeftIdx(double coord, TAxis *axis) const
void GetFormMatrix(TMatrixD &form) const
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
double ana::AtmConstraint::CalculateInterp ( double  normX,
double  normY,
TMatrixD  coefs 
) const
protected
double ana::AtmConstraint::ChiSq ( osc::IOscCalcAdjustable osc,
const SystShifts syst = SystShifts::Nominal() 
) const
overridevirtual

Reimplemented from ana::IExperiment.

Definition at line 68 of file AtmConstraint.cxx.

References BiCubicInterp(), osc::_IOscCalcAdjustable< T >::GetDmsq32(), osc::_IOscCalcAdjustable< T >::GetTh23(), std::sin(), and util::sqr().

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  }
double BiCubicInterp(double x, double y) const
virtual T GetTh23() const
Definition: IOscCalc.h:94
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
virtual T GetDmsq32() const
Definition: IOscCalc.h:91
T sin(T number)
Definition: d0nt_math.hpp:132
void ana::AtmConstraint::GetFormMatrix ( TMatrixD form) const
protected

Definition at line 162 of file AtmConstraint.cxx.

Referenced by BiCubicInterp().

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  }
int ana::AtmConstraint::GetGridLowerLeftIdx ( double  coord,
TAxis *  axis 
) const
protected

Definition at line 119 of file AtmConstraint.cxx.

Referenced by BiCubicInterp().

121  {
122  return axis->FindBin(coord + 0.5*axis->GetBinWidth(coord)) - 2;
123  }
void ana::AtmConstraint::GetGridMatrix ( TH2 *  ll,
int  lowerX,
int  lowerY,
TMatrixD grid 
) const
protected

Definition at line 126 of file AtmConstraint.cxx.

References MakeMiniprodValidationCuts::f, galleryMaker::grid, MECModelEnuComparisons::i, and calib::j.

Referenced by BiCubicInterp().

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  }
const double j
Definition: BetheBloch.cxx:29
virtual stan::math::var ana::IExperiment::LogLikelihood ( osc::IOscCalcAdjustableStan osc,
const SystShifts syst = kNoShift 
) const
inlinevirtualinherited

Reimplemented in test::GaussQuadExperiment, ana::SingleSampleExperiment, and ana::MultiExperiment.

Definition at line 25 of file IExperiment.h.

References ana::assert(), dir, ana::IExperiment::SaveTo(), and string.

Referenced by ana::StanFitter::log_prob().

27  {
28  assert(false && "unimplemented");
29  return 0;
30  };
assert(nhit_max >=nhit_nbins)
void ana::IExperiment::SaveTo ( TDirectory *  dir,
const std::string name 
) const
virtualinherited

Reimplemented in ana::SingleSampleExperiment, ana::MultiExperiment, ana::CountingExperiment, ana::SolarConstraints, and ana::ReactorExperiment.

Definition at line 32 of file IExperiment.cxx.

References ana::assert().

Referenced by ana::IExperiment::LogLikelihood().

33  {
34  assert(0 && "Not implemented");
35  }
assert(nhit_max >=nhit_nbins)

Member Data Documentation

bool ana::AtmConstraint::fIsInvHiNeg
protected

Definition at line 39 of file AtmConstraint.h.

Referenced by AtmConstraint(), and BiCubicInterp().

TH2 * ana::AtmConstraint::fLLInv
protected

Definition at line 38 of file AtmConstraint.h.

Referenced by AtmConstraint(), BiCubicInterp(), and ~AtmConstraint().

TH2* ana::AtmConstraint::fLLNorm
protected

Definition at line 38 of file AtmConstraint.h.

Referenced by AtmConstraint(), BiCubicInterp(), and ~AtmConstraint().

bool ana::AtmConstraint::fScaleAxes
protected

Definition at line 40 of file AtmConstraint.h.

Referenced by AtmConstraint(), and BiCubicInterp().


The documentation for this class was generated from the following files: