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

Compare a single data spectrum to the MC + cosmics expectation. More...

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

Inheritance diagram for ana::LikelihoodCovMxExperiment:
ana::IExperiment

Public Member Functions

 LikelihoodCovMxExperiment (std::vector< covmx::Sample > samples, covmx::CovarianceMatrix *covmx, double epsilon=1e-5, double lambdazero=1, double nu=10)
 
virtual ~LikelihoodCovMxExperiment ()
 
virtual double ChiSq (osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const
 
double GetStatChiSq ()
 
double GetSystChiSq ()
 
double GetResidual ()
 
double GetIteration ()
 
void SaveHists (bool opt=true)
 
TH1D * GetBetaHist ()
 
TH1D * GetMuHist ()
 
TH1D * GetBetaMuHist ()
 
TH1D * GetPredShiftedHist ()
 
TH1D * GetPredUnshiftedHist ()
 
TH1D * GetDataHist ()
 
void SetVerbose (bool opt)
 
void LLUseROOT (bool val=true)
 
void LLUseSVD (bool val=true)
 
void SetShapeOnly (bool val=true)
 
TH1I * InversionTimeHist ()
 
void SetInversionTimeHist (bool timeit=true, int tmax=1000, int tmin=0)
 
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

Eigen::ArrayXd GetExpectedSpectrum () const
 
double GetChiSq (Eigen::ArrayXd e) const
 
void InitialiseBetas () const
 
bool MaskBetas () const
 
void GetGradAndHess () const
 
void GetReducedGradAndHess () const
 
double LikelihoodCovMxNewton () const
 

Protected Attributes

std::vector< covmx::SamplefSamples
 
covmx::CovarianceMatrixfCovMx
 
Eigen::ArrayXd fMu
 
Eigen::ArrayXd fBeta
 
std::vector< bool > fBetaMask
 
Eigen::ArrayXd fData
 
Eigen::ArrayXd fCosmic
 
Eigen::VectorXd fGrad
 
Eigen::MatrixXd fHess
 
Eigen::VectorXd fGradReduced
 
Eigen::MatrixXd fHessReduced
 
unsigned int fNBins
 
unsigned int fNBinsFull
 
std::vector< unsigned intfNBinsPerSample
 
std::vector< std::vector< covmx::Component > > fComps
 
double fLambdaZero
 
double fNu
 Levenberg-Marquardt starting lambda. More...
 
double fLambda
 Levenberg-Marquardt nu. More...
 
bool fLLUseROOT
 Levenberg-Marquardt lambda. More...
 
bool fLLUseSVD
 
TH1D * fBetaHist
 
TH1D * fMuHist
 
TH1D * fBetaMuHist
 
TH1D * fPredShiftedHist
 
TH1D * fPredUnshiftedHist
 
TH1D * fDataHist
 
Eigen::MatrixXd fMxInv
 
bool fShapeOnly
 
bool fVerbose
 
double fStatChiSq
 
double fSystChiSq
 
double fResidual
 
int fIteration
 

Detailed Description

Compare a single data spectrum to the MC + cosmics expectation.

Definition at line 13 of file LikelihoodCovMxExperiment.h.

Constructor & Destructor Documentation

ana::LikelihoodCovMxExperiment::LikelihoodCovMxExperiment ( std::vector< covmx::Sample samples,
covmx::CovarianceMatrix covmx,
double  epsilon = 1e-5,
double  lambdazero = 1,
double  nu = 10 
)
Parameters
samplesSource of concatenated oscillated MC beam predictions
covgenSource of predicted covariance matrix epsilon value to add onto the diagonal to aid matrix inversion

Definition at line 53 of file LikelihoodCovMxExperiment.cxx.

References om::cout, allTimeWatchdog::endl, epsilon, fBeta, fBetaMask, fComps, fCosmic, fCovMx, fData, fGrad, fHess, fMu, fMxInv, fNBins, fNBinsFull, fNBinsPerSample, fResidual, fSamples, fVerbose, GetComponents(), ana::covmx::CovarianceMatrix::GetFullCovMx(), MECModelEnuComparisons::i, ana::kLivetime, livetime, pot, and Zero().

55  : fSamples(samples), fCovMx(covmx), fNBins(0), fNBinsFull(0), fLambdaZero(lambdazero),
56  fNu(nu), fLLUseROOT(false), fLLUseSVD(false), fBetaHist(nullptr), fMuHist(nullptr),
57  fBetaMuHist(nullptr), fPredShiftedHist(nullptr), fPredUnshiftedHist(nullptr), fDataHist(nullptr),
58  /*fMxInv(nullptr), */fShapeOnly(false), fVerbose(false)
59  {
60  for (Sample s: fSamples) {
61  fNBins += s.GetBinning().NBins();
62  fNBinsFull += s.GetBinning().NBins() * GetComponents(s).size();
63  fNBinsPerSample.push_back(s.GetBinning().NBins());
64  fComps.push_back(GetComponents(s));
65  }
66 
67  if (fCovMx) {
68  // Invert matrix ahead of time
69  MatrixXd mx = CovarianceMatrix::ForcePosDef(fCovMx->GetFullCovMx(), epsilon);
70  fMxInv = mx.inverse();
71  MatrixXd id = mx;
72  id.setIdentity();
73  fResidual = ((mx.transpose() * fMxInv) - id).maxCoeff();
74  if (fVerbose) cout << "Eigen matrix inversion residual is " << fResidual << endl;
75  }
76 
77  fGrad.resize(fNBinsFull);
78  fHess.resize(fNBinsFull, fNBinsFull);
79 
80  // Set up beta and mu
81  fMu.resize(fNBinsFull);
82  fBeta.resize(fNBinsFull);
83  fBetaMask.resize(fNBinsFull);
84 
85  // Fill the data vector
88  size_t i = 0;
89  for (size_t iS = 0; iS < fSamples.size(); ++iS) {
90  double pot = fSamples[iS].GetData()->POT();
91  double livetime = fSamples[iS].GetData()->Livetime();
92  TH1D* hData = fSamples[iS].GetData()->ToTH1(pot);
93  TH1D* hCos = fSamples[iS].HasCosmic() ?
94  fSamples[iS].GetCosmic()->ToTH1(livetime, kLivetime) : nullptr;
95  for (size_t iBin = 1; iBin <= fNBinsPerSample[iS]; ++iBin) {
96  fData(i) = hData->GetBinContent(iBin);
97  fCosmic(i) = hCos ? hCos->GetBinContent(iBin) : 0;
98  ++i;
99  } // for bin i
100  delete hData;
101  delete hCos;
102  } // for sample i
103  }
std::vector< std::vector< covmx::Component > > fComps
std::vector< std::tuple< ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t > > GetComponents()
double fNu
Levenberg-Marquardt starting lambda.
Eigen::MatrixXd GetFullCovMx()
const XML_Char * s
Definition: expat.h:262
#define pot
OStream cout
Definition: OStream.cxx:6
std::vector< unsigned int > fNBinsPerSample
void Zero()
double livetime
Definition: saveFDMCHists.C:21
double epsilon
std::vector< covmx::Sample > fSamples
bool fLLUseROOT
Levenberg-Marquardt lambda.
virtual ana::LikelihoodCovMxExperiment::~LikelihoodCovMxExperiment ( )
inlinevirtual

Definition at line 23 of file LikelihoodCovMxExperiment.h.

References ChiSq(), and ana::SystShifts::Nominal().

23 {};

Member Function Documentation

double ana::LikelihoodCovMxExperiment::ChiSq ( osc::IOscCalcAdjustable osc,
const SystShifts syst = SystShifts::Nominal() 
) const
virtual

Reimplemented from ana::IExperiment.

Definition at line 106 of file LikelihoodCovMxExperiment.cxx.

References e, fBeta, fComps, fCovMx, fData, fMu, fNBins, fNBinsFull, fPredUnshiftedHist, fSamples, ana::covmx::CovarianceMatrix::GetBinning(), GetExpectedSpectrum(), MECModelEnuComparisons::i, calib::j, LikelihoodCovMxNewton(), ana::LogLikelihood(), ana::Binning::NBins(), pot, and ana::Spectrum::ToTH1().

Referenced by ~LikelihoodCovMxExperiment().

108  {
109  DontAddDirectory guard;
110 
111  // Get MC and data spectra
112  // fMu(fNBinsFull);
113  size_t i = 0;
114  for (size_t iS = 0; iS < fSamples.size(); ++iS) {
115  double pot = fSamples[iS].GetData()->POT();
116  SystShifts shift = fSamples[i].GetSystShifts(syst);
117  for (Component comp: fComps[iS]) {
118  Spectrum spec = fSamples[iS].GetPrediction()->PredictComponentSyst(osc,
119  shift, get<0>(comp), get<1>(comp), get<2>(comp));
120  TH1D* hMu = spec.ToTH1(pot);
121  for (size_t iBin = 1; iBin <= (size_t)hMu->GetNbinsX(); ++iBin) {
122  fMu(i) = hMu->GetBinContent(iBin);
123  ++i;
124  } // for bin
125  delete hMu;
126  } // for component
127  } // for sample
128  for (size_t i = 0; i < fNBinsFull; ++i) fBeta(i) = 1.;
129 
130  // Do stats-only calculation if no matrix
131  if (!fCovMx) {
132  ArrayXd e = GetExpectedSpectrum();
133  return ana::LogLikelihood(e, fData);
134  }
135 
136  // Number of bins
137  if (fCovMx && fNBins != (size_t)fCovMx->GetBinning().NBins()) {
138  throw runtime_error(
139  "Number of bins in predictions does not match covariance matrix!");
140  }
141 
142  // Fill histogram before beta shifting
143  if (fPredUnshiftedHist) {
144  ArrayXd e = GetExpectedSpectrum();
145  for (size_t j = 0; j < fNBins; ++j) {
146  fPredUnshiftedHist->SetBinContent(j+1, e(j));
147  }
148  } // if filling unshifted spectrum
149 
150  return LikelihoodCovMxNewton();
151 
152  } // function LikelihoodCovMxExperiment::ChiSq
std::vector< std::vector< covmx::Component > > fComps
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
#define pot
const double j
Definition: BetheBloch.cxx:29
std::vector< float > Spectrum
Definition: Constants.h:610
int NBins() const
Definition: Binning.h:29
std::vector< covmx::Sample > fSamples
Float_t e
Definition: plot.C:35
TH1D* ana::LikelihoodCovMxExperiment::GetBetaHist ( )
inline

Definition at line 35 of file LikelihoodCovMxExperiment.h.

References fBetaHist.

TH1D* ana::LikelihoodCovMxExperiment::GetBetaMuHist ( )
inline

Definition at line 37 of file LikelihoodCovMxExperiment.h.

References fBetaMuHist.

double ana::LikelihoodCovMxExperiment::GetChiSq ( Eigen::ArrayXd  e) const
protected

Definition at line 182 of file LikelihoodCovMxExperiment.cxx.

References fBeta, fData, fMxInv, fStatChiSq, fSystChiSq, and ana::LogLikelihood().

Referenced by LikelihoodCovMxNewton(), and SetShapeOnly().

182  {
183 
184  // the statistical likelihood of the shifted spectrum...
186 
187  // ...plus the penalty the prediction picks up from its covariance
188  VectorXd vec = fBeta - 1.;
189  fSystChiSq = vec.transpose() * fMxInv * vec;
190  return fStatChiSq + fSystChiSq;
191 
192  } // function LikelihoodCovMxExperiment::GetChiSq
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
Eigen::VectorXd vec
Float_t e
Definition: plot.C:35
TH1D* ana::LikelihoodCovMxExperiment::GetDataHist ( )
inline

Definition at line 40 of file LikelihoodCovMxExperiment.h.

References fDataHist.

ArrayXd ana::LikelihoodCovMxExperiment::GetExpectedSpectrum ( ) const
protected
Parameters
muVector of expected events
betaVector of shifted beta values

Definition at line 155 of file LikelihoodCovMxExperiment.cxx.

References fBeta, fComps, fCosmic, fMu, fNBins, fNBinsPerSample, fSamples, runNovaSAM::ret, and Zero().

Referenced by ChiSq(), LikelihoodCovMxNewton(), and SetShapeOnly().

155  {
156 
157  ArrayXd ret = ArrayXd::Zero(fNBins);
158  unsigned int fullOffset(0), scaledOffset(0);
159 
160  for (size_t iS = 0; iS < fSamples.size(); ++iS) { // loop over samples
161  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
162  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) { // loop over bins
163  unsigned int iScaled = scaledOffset + iB;
164  unsigned int iFull = fullOffset + (iC*fNBinsPerSample[iS]) + iB;
165 
166  ret(iScaled) += fMu(iFull)*fBeta(iFull);
167  } // for bin
168  } // for component
169  // Now loop one more time to add cosmics
170  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) {
171  unsigned int iScaled = scaledOffset + iB;
172  ret(iScaled) += fCosmic(iScaled);
173  } // for bin
174  scaledOffset += fNBinsPerSample[iS]; // increase offset in scaled matrix
175  fullOffset += fComps[iS].size() * fNBinsPerSample[iS];
176  } // for sample
177  return ret;
178 
179  } // function LikelihoodCovMxExperiment::GetExpectedSpectrum
std::vector< std::vector< covmx::Component > > fComps
std::vector< unsigned int > fNBinsPerSample
void Zero()
std::vector< covmx::Sample > fSamples
void ana::LikelihoodCovMxExperiment::GetGradAndHess ( ) const
protected

Definition at line 268 of file LikelihoodCovMxExperiment.cxx.

References e, fBeta, fComps, fCosmic, fData, fGrad, fHess, fLambda, fMu, fMxInv, fNBinsFull, fNBinsPerSample, fSamples, MECModelEnuComparisons::i, submit_syst::y, and test::z.

Referenced by GetReducedGradAndHess(), and SetShapeOnly().

268  {
269 
270  size_t iFullOffset(0.), iScaledOffset(0.);
271  for (size_t iS = 0; iS < fSamples.size(); ++iS) { // loop over samples
272  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
273  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) { // loop over bins
274  // Get the current scaled & full bin
275  unsigned int iScaled = iScaledOffset+iB;
276  unsigned int iFull = iFullOffset+(iC*fNBinsPerSample[iS])+iB;
277  double y = fCosmic(iScaled); // calculate y
278  for (size_t jC = 0; jC < fComps[iS].size(); ++jC) {
279  unsigned int jFull = iFullOffset+(jC*fNBinsPerSample[iS])+iB;
280  y += fMu(jFull) * fBeta(jFull);
281  } // for component j
282  if (y < 1e-40) y = 1e-40; // clamp y
283  double z(0.); // calculate z
284  for (size_t jFull = 0; jFull < fNBinsFull; ++jFull) {
285  z += (fBeta(jFull)-1)*fMxInv(iFull, jFull);
286  }
287 
288  // Calculate gradient
289  fGrad(iFull) = 2 * (fMu(iFull)*(1.0-(fData(iScaled)/y)) + z);
290 
291  // Populate Hessian matrix
292  size_t jFullOffset(0.), jScaledOffset(0.);
293  for (size_t jS = 0; jS < fSamples.size(); ++jS) { // loop over samples
294  for (size_t jC = 0; jC < fComps[jS].size(); ++jC) { // loop over components
295  for (size_t jB = 0; jB < fNBinsPerSample[jS]; ++jB) { // loop over bins
296  // Get the current scaled & full bin
297  unsigned int jScaled = jScaledOffset + jB;
298  unsigned int jFull = jFullOffset+(jC*fNBinsPerSample[jS])+jB;
299  if (iScaled != jScaled) {
300  fHess(iFull, jFull) = 2 * fMxInv(iFull, jFull);
301  } else {
302  fHess(iFull, jFull) = 2 * (fMu(iFull)*fMu(jFull)*(fData(iScaled)/(y*y)) + fMxInv(iFull, jFull));
303  }
304  } // for bin j
305  } // for component j
306  jScaledOffset += fNBinsPerSample[jS]; // increase offset in scaled matrix
307  jFullOffset += fComps[jS].size() * fNBinsPerSample[jS];
308  } // for sample j
309  } // for bin i
310  } // for component i
311  iScaledOffset += fNBinsPerSample[iS]; // increase offset in scaled matrix
312  iFullOffset += fComps[iS].size() * fNBinsPerSample[iS];
313  } // for sample i
314 
315  // Apply transformation and Levenberg-Marquardt
316  for (size_t i = 0; i < fNBinsFull; ++i) {
317  fHess(i, i) *= 1. + fLambda;
318  fGrad(i) *= -1;
319  }
320 
321  } // function LikelihoodCovMxExperiment::GetGradAndHess
double fLambda
Levenberg-Marquardt nu.
std::vector< std::vector< covmx::Component > > fComps
z
Definition: test.py:28
std::vector< unsigned int > fNBinsPerSample
std::vector< covmx::Sample > fSamples
Float_t e
Definition: plot.C:35
double ana::LikelihoodCovMxExperiment::GetIteration ( )
inline
TH1D* ana::LikelihoodCovMxExperiment::GetMuHist ( )
inline

Definition at line 36 of file LikelihoodCovMxExperiment.h.

References fMuHist.

TH1D* ana::LikelihoodCovMxExperiment::GetPredShiftedHist ( )
inline

Definition at line 38 of file LikelihoodCovMxExperiment.h.

References fPredShiftedHist.

TH1D* ana::LikelihoodCovMxExperiment::GetPredUnshiftedHist ( )
inline

Definition at line 39 of file LikelihoodCovMxExperiment.h.

References fPredUnshiftedHist.

void ana::LikelihoodCovMxExperiment::GetReducedGradAndHess ( ) const
protected

Definition at line 324 of file LikelihoodCovMxExperiment.cxx.

References plot_validation_datamc::c, confusionMatrixTree::count, fBetaMask, fGrad, fGradReduced, fHess, fHessReduced, fNBinsFull, GetGradAndHess(), MECModelEnuComparisons::i, and calib::j.

Referenced by LikelihoodCovMxNewton(), and SetShapeOnly().

324  {
325 
326  // Mask off bins
327  int maskedSize = count(fBetaMask.begin(), fBetaMask.end(), true);
328 
329  GetGradAndHess();
330 
331  if (fGradReduced.rows() != maskedSize) {
332  fGradReduced.resize(maskedSize);
333  fHessReduced.resize(maskedSize, maskedSize);
334  }
335 
336  size_t c = 0;
337  for (size_t i = 0; i < fNBinsFull; ++i) {
338  if (fBetaMask[i]) {
339  fGradReduced(c) = fGrad(i);
340  ++c;
341  }
342  }
343  size_t ci = 0;
344  for (size_t i = 0; i < fNBinsFull; ++i) {
345  if (fBetaMask[i]) {
346  size_t cj = 0;
347  for (size_t j = 0; j < fNBinsFull; ++j) {
348  if (fBetaMask[j]) {
349  fHessReduced(ci, cj) = fHess(i, j);
350  ++cj;
351  } // if j unmasked
352  } // for j
353  ++ci;
354  } // if i unmasked
355  } // for i
356 
357  } // function LikelihoodCovMxExperiment::GetReducedGradAndHess
const double j
Definition: BetheBloch.cxx:29
double ana::LikelihoodCovMxExperiment::GetResidual ( )
inline

Definition at line 30 of file LikelihoodCovMxExperiment.h.

References fResidual.

double ana::LikelihoodCovMxExperiment::GetStatChiSq ( )
inline

Definition at line 28 of file LikelihoodCovMxExperiment.h.

References fStatChiSq.

double ana::LikelihoodCovMxExperiment::GetSystChiSq ( )
inline

Definition at line 29 of file LikelihoodCovMxExperiment.h.

References fSystChiSq.

void ana::LikelihoodCovMxExperiment::InitialiseBetas ( ) const
protected

Definition at line 195 of file LikelihoodCovMxExperiment.cxx.

References e, fBeta, fComps, fCosmic, fData, fMu, fNBinsPerSample, fSamples, sum, and febshutoff_auto::val.

Referenced by SetShapeOnly().

195  {
196  double eps(1e-40);
197  unsigned int fullOffset(0), scaledOffset(0);
198  for (size_t iS = 0; iS < fSamples.size(); ++iS) { // loop over samples
199  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) { // loop over bins
200  unsigned int iScaled = scaledOffset+iB;
201  double sum = fCosmic(iScaled);
202  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
203  unsigned int iFull = fullOffset+(iC*fNBinsPerSample[iS])+iB;
204  sum += fMu(iFull);
205  } // for component
206  double val = fData(iScaled) / (sum + eps);
207  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
208  unsigned int iFull = fullOffset+(iC*fNBinsPerSample[iS])+iB;
209  fBeta(iFull) = val;
210  } // for component
211  } // for bin
212  scaledOffset += fNBinsPerSample[iS]; // increase offset in scaled matrix
213  fullOffset += fComps[iS].size() * fNBinsPerSample[iS];
214  } // for sample
215 
216  } // function LikelihoodCovMxExperiment::InitialiseBetas
std::vector< std::vector< covmx::Component > > fComps
std::vector< unsigned int > fNBinsPerSample
std::vector< covmx::Sample > fSamples
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
TH1I* ana::LikelihoodCovMxExperiment::InversionTimeHist ( )

Referenced by SetShapeOnly().

double ana::LikelihoodCovMxExperiment::LikelihoodCovMxNewton ( ) const
protected

Definition at line 360 of file LikelihoodCovMxExperiment.cxx.

References confusionMatrixTree::count, hadd_many_files::counter, om::cout, e, febshutoff_auto::end, allTimeWatchdog::endl, stan::math::fabs(), fBeta, fBetaHist, fBetaMask, fBetaMuHist, fGradReduced, fHessReduced, fIteration, fLambda, fMu, fMuHist, fNBins, fNBinsFull, fPredShiftedHist, fStatChiSq, fSystChiSq, fVerbose, GetChiSq(), GetExpectedSpectrum(), GetReducedGradAndHess(), MECModelEnuComparisons::i, std::isnan(), MaskBetas(), norm, recentWatchdog::now, prev(), and febshutoff_auto::start.

Referenced by ChiSq(), and SetShapeOnly().

360  {
361 
363 
364  // Start with all systematic shifts unmasked
365  for (size_t i = 0; i < fNBinsFull; ++i) fBetaMask[i] = true;
366 
367  // We're trying to solve for the best expectation in each bin and each component. A good seed
368  // value is 1 (no shift).
369  // Note: beta will contain the best fit systematic shifts at the end of the process.
370  // This should be saved to show true post-fit agreement
371  const int maxIters = 1e7;
372  fIteration = 0;
373  // InitialiseBetas();
374  MaskBetas();
375  fLambda = 0; // Initialise lambda at starting value
376 
377  ArrayXd e = GetExpectedSpectrum();
378  double prev = GetChiSq(e);
379  cout << "Before minimization, chisq is " << prev
380  << " (" << fStatChiSq << " stat, " << fSystChiSq << " syst)" << endl;
381 
382  while (true) {
383  ++fIteration;
384  if (fIteration > maxIters) {
385  std::ostringstream err;
386  err << "No convergence after " << maxIters << " iterations! Quitting out of the infinite loop.";
387  throw runtime_error(err.str());
388  }
389 
390  // If we have negative betas, mask them out
392  VectorXd deltaBeta = fHessReduced.ldlt().solve(fGradReduced);
393  if (deltaBeta.array().isNaN().any()) {
394  cout << "LDLT solution failed! Switching to SVD." << endl;
395  deltaBeta = fHessReduced.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(fGradReduced);
396  }
397 
398  if (fVerbose) {
399  double relativeError = (fHessReduced * deltaBeta - fGradReduced).norm() / fGradReduced.norm();
400  cout << "Relative error of linear solution is " << relativeError << endl;
401  }
402 
403  if (deltaBeta.array().isNaN().any()) {
404  // cout << "NaN found in delta beta vector! Setting lambda to one!" << endl;
405  // fLambda = 1;
406  // continue;
407  cout << "NaN found in delta beta vector! Dumping mus:" << endl;
408  // size_t counter = 0;
409  for (size_t i = 0; i < fNBinsFull; ++i) {
410  cout << " Mu " << i << ": " << fMu(i) << endl;
411  // if (fBetaMask[i]) {
412  // cout << " Beta for bin " << i << ": " << deltaBeta(counter) << endl;
413  // ++counter;
414  // }
415  }
416  throw runtime_error("NaN found in delta beta vector - exiting.");
417  }
418 
419  size_t counter = 0;
420  for (size_t i = 0; i < fNBinsFull; ++i) {
421  if (fBetaMask[i]) {
422  fBeta(i) += deltaBeta(counter);
423  if (fBeta(i) < 0) { // if we pulled negative, mask this beta off
424  fBeta(i) = 0;
425  fBetaMask[i] = false;
426  }
427  ++counter;
428  }
429  } // for bin
430 
431  // Predict collapsed spectrum
432  e = GetExpectedSpectrum();
433 
434  // Update the chisq
435  double chisq = GetChiSq(e);
436 
437  // Switch to LMA if first iteration led to masked betas
438  // if (fIteration == 1 && maskChanged) {
439  // cout << "Mask changed! Enabling LMA." << endl;
440  // fLambda = fLambdaZero;
441  // }
442 
443  // Bump down the LMA lambda
444  // if (fIteration > 1 && chisq < prev) fLambda /= fNu;
445 
446  if (isnan(chisq)) throw runtime_error("ChiSq value is nan! Exiting.");
447 
448  // if (fVerbose) {
449  cout << "Iteration " << fIteration << ", chisq " << chisq
450  << " (" << fStatChiSq << " stat, " << fSystChiSq << " syst)" << endl;
451  // }
452 
453  if (fSystChiSq < 0) {
454  throw runtime_error("Negative systematic penalty!");
455  }
456 
457  // If the updates didn't change anything at all then we're done
458  double change = fabs(chisq-prev);
459  // Converge to tenth significant figure or tenth decimal place, whichever is larger
460  if (change/chisq < 1e-10 || change < 1e-10) {
461 
462  // if (fVerbose) {
463  cout << "Minimisation round finished with chisq " << chisq << ", bin mask:";
464  for (size_t i = 0; i < fNBinsFull; ++i) {
465  if (!fBetaMask[i]) cout << " " << i+1;
466  }
467  cout << endl;
468  // }
469 
470  // if (MaskBetas()) { // re-mask betas - if the mask changed, go again
471  // fLambda = fLambdaZero;
472  // }
473 
474  // else {
475  if (!MaskBetas()) {
476 
477  // if (fVerbose) {
478  cout << "Converged to " << chisq << " (" << fStatChiSq << " stat, " << fSystChiSq
479  << " syst) after " << fIteration << " iterations and ";
481  double secs = duration_cast<seconds>(end-start).count();
482  if (secs < 1) {
483  double millisecs = duration_cast<milliseconds>(end-start).count();
484  cout << millisecs << " ms." << endl;
485  } else if (secs < 60) cout << secs << " seconds." << endl;
486  else {
487  double mins = duration_cast<minutes>(end-start).count();
488  cout << mins << " minutes." << endl;
489  }
490  // }
491 
492  // Fill histograms if necessary
493  if (fBetaHist) for (size_t i = 0; i < fNBinsFull; ++i) fBetaHist->SetBinContent(i+1, fBeta(i));
494  if (fMuHist) for (size_t i = 0; i < fNBinsFull; ++i) fMuHist->SetBinContent(i+1, fMu(i));
495  if (fBetaMuHist) for (size_t i = 0; i < fNBinsFull; ++i) fBetaMuHist->SetBinContent(i+1, fBeta(i)*fMu(i));
496  if (fPredShiftedHist) for (size_t i = 0; i < fNBins; ++i) fPredShiftedHist->SetBinContent(i+1, e[i]);
497 
498  return chisq;
499  }
500  }
501 
502  prev = chisq;
503 
504  } // end while
505 
506  } // function LikelihoodCovMxExperiment::LogLikelihoodCovMxNewton
double fLambda
Levenberg-Marquardt nu.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
double GetChiSq(Eigen::ArrayXd e) const
void prev()
Definition: show_event.C:91
OStream cout
Definition: OStream.cxx:6
Float_t norm
Float_t e
Definition: plot.C:35
void ana::LikelihoodCovMxExperiment::LLUseROOT ( bool  val = true)
inline

Definition at line 44 of file LikelihoodCovMxExperiment.h.

References fLLUseROOT, and febshutoff_auto::val.

44 { fLLUseROOT = val; };
bool fLLUseROOT
Levenberg-Marquardt lambda.
void ana::LikelihoodCovMxExperiment::LLUseSVD ( bool  val = true)
inline
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)
bool ana::LikelihoodCovMxExperiment::MaskBetas ( ) const
protected

Definition at line 219 of file LikelihoodCovMxExperiment.cxx.

References beta, e, fBeta, fBetaMask, fComps, fCosmic, fData, fMu, fMxInv, fNBinsFull, fNBinsPerSample, fSamples, submit_syst::y, and test::z.

Referenced by LikelihoodCovMxNewton(), and SetShapeOnly().

219  {
220 
221  bool maskChange = false;
222 
223  size_t iFullOffset(0.), iScaledOffset(0.);
224  for (size_t iS = 0; iS < fSamples.size(); ++iS) { // loop over samples
225  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
226  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) { // loop over bins
227  // Get the current scaled & full bin
228  unsigned int iScaled = iScaledOffset+iB;
229  unsigned int iFull = iFullOffset+(iC*fNBinsPerSample[iS])+iB;
230 
231  double y = fCosmic(iScaled);
232  for (size_t jC = 0; jC < fComps[iS].size(); ++jC) {
233  if (iC != jC) { // beta mu from this bin is 0
234  unsigned int jFull = iFullOffset+(jC*fNBinsPerSample[iS])+iB;
235  y += fMu(jFull) * fBeta(jFull);
236  }
237  } // for component j
238  if (y < 1e-40) y = 1e-40; // clamp y
239 
240  double z(0.); // calculate z
241  for (size_t jFull = 0; jFull < fNBinsFull; ++jFull) {
242  double beta = (jFull == iFull) ? 0 : fBeta[jFull]; // set this beta to 0
243  z += (beta-1)*fMxInv(iFull, jFull);
244  }
245 
246  // Calculate gradient
247  double gradZero = 2 * (fMu(iFull)*(1.0-(fData(iScaled)/y)) + z);
248  if (gradZero > 0 && fBetaMask[iFull]) { // if we're masking a bin off
249  fBetaMask[iFull] = false;
250  maskChange = true;
251  }
252  if (gradZero <= 0 && !fBetaMask[iFull]) { // if we're masking a bin on
253  fBetaMask[iFull] = true;
254  fBeta(iFull) = 0;
255  maskChange = true;
256  }
257  } // for bin i
258  } // for component i
259  iScaledOffset += fNBinsPerSample[iS]; // increase offset in scaled matrix
260  iFullOffset += fComps[iS].size() * fNBinsPerSample[iS];
261  } // for sample i
262 
263  return maskChange;
264 
265  } // function LikelihoodCovMxExperiment::MaskBetas
std::vector< std::vector< covmx::Component > > fComps
Double_t beta
z
Definition: test.py:28
std::vector< unsigned int > fNBinsPerSample
std::vector< covmx::Sample > fSamples
Float_t e
Definition: plot.C:35
void ana::LikelihoodCovMxExperiment::SaveHists ( bool  opt = true)

Definition at line 509 of file LikelihoodCovMxExperiment.cxx.

References ana::Binning::Edges(), fBetaHist, fBetaMuHist, fCovMx, fData, fDataHist, fMuHist, fNBins, fPredShiftedHist, fPredUnshiftedHist, ana::covmx::CovarianceMatrix::GetBinning(), ana::covmx::CovarianceMatrix::GetFullBinning(), MECModelEnuComparisons::i, ana::Binning::NBins(), and ana::UniqueName().

Referenced by GetIteration().

509  {
510  if (!fCovMx) return;
511  // Start by clearing up existing beta histogram
512  if (fBetaHist) {
513  delete fBetaHist;
514  fBetaHist = nullptr;
515  }
516  if (fMuHist) {
517  delete fMuHist;
518  fMuHist = nullptr;
519  }
520  if (fBetaMuHist) {
521  delete fBetaMuHist;
522  fBetaMuHist = nullptr;
523  }
524  if (fPredShiftedHist) {
525  delete fPredShiftedHist;
526  fPredShiftedHist = nullptr;
527  }
528  if (fPredUnshiftedHist) {
529  delete fPredUnshiftedHist;
530  fPredUnshiftedHist = nullptr;
531  }
532  if (fDataHist) {
533  delete fDataHist;
534  fDataHist = nullptr;
535  }
536  // Now re-instantiate it if necessary
537  if (opt) {
538  Binning fullBins = fCovMx->GetFullBinning();
539  Binning scaledBins = fCovMx->GetBinning();
540  fBetaHist = new TH1D(UniqueName().c_str(), ";;", fullBins.NBins(), &fullBins.Edges()[0]);
541  fMuHist = new TH1D(UniqueName().c_str(), ";;", fullBins.NBins(), &fullBins.Edges()[0]);
542  fBetaMuHist = new TH1D(UniqueName().c_str(), ";;", fullBins.NBins(), &fullBins.Edges()[0]);
543  fPredShiftedHist = new TH1D(UniqueName().c_str(), ";;", scaledBins.NBins(), &scaledBins.Edges()[0]);
544  fPredUnshiftedHist = new TH1D(UniqueName().c_str(), ";;", scaledBins.NBins(), &scaledBins.Edges()[0]);
545  fDataHist = new TH1D(UniqueName().c_str(), ";;", scaledBins.NBins(), &scaledBins.Edges()[0]);
546  for (size_t i = 0; i < fNBins; ++i) {
547  fDataHist->SetBinContent(i+1, fData(i)); // Populate data hist
548  }
549  }
550  } // function LikelihoodCovMxExperiment::SaveHists
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
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)
void ana::LikelihoodCovMxExperiment::SetInversionTimeHist ( bool  timeit = true,
int  tmax = 1000,
int  tmin = 0 
)
Parameters
timeitWhether or not to save inversion time in a histogram
tmaxUpper boundary on time hist in milliseconds
tminLower boundary on time hist in milliseconds

Referenced by SetShapeOnly().

void ana::LikelihoodCovMxExperiment::SetShapeOnly ( bool  val = true)
inline
void ana::LikelihoodCovMxExperiment::SetVerbose ( bool  opt)
inline

Member Data Documentation

Eigen::ArrayXd ana::LikelihoodCovMxExperiment::fBeta
mutableprotected
TH1D* ana::LikelihoodCovMxExperiment::fBetaHist
protected

Definition at line 94 of file LikelihoodCovMxExperiment.h.

Referenced by GetBetaHist(), LikelihoodCovMxNewton(), and SaveHists().

std::vector<bool> ana::LikelihoodCovMxExperiment::fBetaMask
mutableprotected
TH1D* ana::LikelihoodCovMxExperiment::fBetaMuHist
protected

Definition at line 96 of file LikelihoodCovMxExperiment.h.

Referenced by GetBetaMuHist(), LikelihoodCovMxNewton(), and SaveHists().

std::vector<std::vector<covmx::Component> > ana::LikelihoodCovMxExperiment::fComps
protected
Eigen::ArrayXd ana::LikelihoodCovMxExperiment::fCosmic
protected
covmx::CovarianceMatrix* ana::LikelihoodCovMxExperiment::fCovMx
protected

Definition at line 71 of file LikelihoodCovMxExperiment.h.

Referenced by ChiSq(), LikelihoodCovMxExperiment(), and SaveHists().

Eigen::ArrayXd ana::LikelihoodCovMxExperiment::fData
protected
TH1D* ana::LikelihoodCovMxExperiment::fDataHist
protected

Definition at line 99 of file LikelihoodCovMxExperiment.h.

Referenced by GetDataHist(), and SaveHists().

Eigen::VectorXd ana::LikelihoodCovMxExperiment::fGrad
mutableprotected
Eigen::VectorXd ana::LikelihoodCovMxExperiment::fGradReduced
mutableprotected

Definition at line 80 of file LikelihoodCovMxExperiment.h.

Referenced by GetReducedGradAndHess(), and LikelihoodCovMxNewton().

Eigen::MatrixXd ana::LikelihoodCovMxExperiment::fHess
mutableprotected
Eigen::MatrixXd ana::LikelihoodCovMxExperiment::fHessReduced
mutableprotected

Definition at line 81 of file LikelihoodCovMxExperiment.h.

Referenced by GetReducedGradAndHess(), and LikelihoodCovMxNewton().

int ana::LikelihoodCovMxExperiment::fIteration
mutableprotected

Definition at line 108 of file LikelihoodCovMxExperiment.h.

Referenced by GetIteration(), and LikelihoodCovMxNewton().

double ana::LikelihoodCovMxExperiment::fLambda
mutableprotected

Levenberg-Marquardt nu.

Definition at line 90 of file LikelihoodCovMxExperiment.h.

Referenced by GetGradAndHess(), and LikelihoodCovMxNewton().

double ana::LikelihoodCovMxExperiment::fLambdaZero
protected

Definition at line 88 of file LikelihoodCovMxExperiment.h.

bool ana::LikelihoodCovMxExperiment::fLLUseROOT
protected

Levenberg-Marquardt lambda.

Definition at line 92 of file LikelihoodCovMxExperiment.h.

Referenced by LLUseROOT().

bool ana::LikelihoodCovMxExperiment::fLLUseSVD
protected

Definition at line 93 of file LikelihoodCovMxExperiment.h.

Referenced by LLUseSVD().

Eigen::ArrayXd ana::LikelihoodCovMxExperiment::fMu
mutableprotected
TH1D* ana::LikelihoodCovMxExperiment::fMuHist
protected

Definition at line 95 of file LikelihoodCovMxExperiment.h.

Referenced by GetMuHist(), LikelihoodCovMxNewton(), and SaveHists().

Eigen::MatrixXd ana::LikelihoodCovMxExperiment::fMxInv
protected
unsigned int ana::LikelihoodCovMxExperiment::fNBins
protected
unsigned int ana::LikelihoodCovMxExperiment::fNBinsFull
protected
std::vector<unsigned int> ana::LikelihoodCovMxExperiment::fNBinsPerSample
protected
double ana::LikelihoodCovMxExperiment::fNu
protected

Levenberg-Marquardt starting lambda.

Definition at line 89 of file LikelihoodCovMxExperiment.h.

TH1D* ana::LikelihoodCovMxExperiment::fPredShiftedHist
protected
TH1D* ana::LikelihoodCovMxExperiment::fPredUnshiftedHist
protected

Definition at line 98 of file LikelihoodCovMxExperiment.h.

Referenced by ChiSq(), GetPredUnshiftedHist(), and SaveHists().

double ana::LikelihoodCovMxExperiment::fResidual
mutableprotected

Definition at line 107 of file LikelihoodCovMxExperiment.h.

Referenced by GetResidual(), and LikelihoodCovMxExperiment().

std::vector<covmx::Sample> ana::LikelihoodCovMxExperiment::fSamples
protected
bool ana::LikelihoodCovMxExperiment::fShapeOnly
protected

Definition at line 101 of file LikelihoodCovMxExperiment.h.

Referenced by SetShapeOnly().

double ana::LikelihoodCovMxExperiment::fStatChiSq
mutableprotected

Definition at line 105 of file LikelihoodCovMxExperiment.h.

Referenced by GetChiSq(), GetStatChiSq(), and LikelihoodCovMxNewton().

double ana::LikelihoodCovMxExperiment::fSystChiSq
mutableprotected

Definition at line 106 of file LikelihoodCovMxExperiment.h.

Referenced by GetChiSq(), GetSystChiSq(), and LikelihoodCovMxNewton().

bool ana::LikelihoodCovMxExperiment::fVerbose
protected

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