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/N21-04-21/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=0, double nu=10)
 
virtual ~LikelihoodCovMxExperiment ()
 
double ChiSq (osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
 
void Reset () const override
 
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)
 
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

void GetExpectedSpectrum () const
 
double GetChiSq () const
 
void InitialiseBetas () const
 
bool MaskBetas () const
 
void GetGradAndHess () const
 
double LikelihoodCovMxNewton () const
 
void EnableLMA () const
 
void ResetLambda () const
 
void DecreaseLambda () const
 
void IncreaseLambda () const
 

Protected Attributes

std::vector< covmx::SamplefSamples
 
covmx::CovarianceMatrixfCovMx
 
Eigen::ArrayXi fIndexMap
 
Eigen::ArrayXd fData
 
Eigen::ArrayXd fDataComponentwise
 
Eigen::ArrayXd fCosmic
 
Eigen::ArrayXd fMu
 
Eigen::ArrayXd fBeta
 
Eigen::ArrayXd fBetaMu
 
std::vector< bool > fBetaMask
 
Eigen::ArrayXd fSpec
 
Eigen::ArrayXd fSpecComponentwise
 
Eigen::VectorXd fGrad
 
Eigen::MatrixXd fHess
 
unsigned int fNBins
 
unsigned int fNBinsFull
 
std::vector< unsigned intfNBinsPerSample
 
std::vector< std::vector< covmx::Component > > fComps
 
bool fUseLMA
 
double fLambdaZero
 Whether to use LMA. More...
 
double fNu
 Levenberg-Marquardt starting lambda. More...
 
double fLambda
 Levenberg-Marquardt nu. More...
 
TH1D * fBetaHist
 Levenberg-Marquardt lambda. More...
 
TH1D * fMuHist
 
TH1D * fBetaMuHist
 
TH1D * fPredShiftedHist
 
TH1D * fPredUnshiftedHist
 
TH1D * fDataHist
 
Eigen::MatrixXd fMxInv
 
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 = 0,
double  nu = 10 
)
Parameters
samplesSource of concatenated oscillated MC beam predictions
covmxCovariance matrix epsilon value to add onto the diagonal to aid matrix inversion

Definition at line 52 of file LikelihoodCovMxExperiment.cxx.

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

54  : fSamples(samples), fCovMx(covmx), fNBins(0), fNBinsFull(0), fLambdaZero(lambdazero),
55  fNu(nu), fBetaHist(nullptr), fMuHist(nullptr),
56  fBetaMuHist(nullptr), fPredShiftedHist(nullptr), fPredUnshiftedHist(nullptr),
57  fDataHist(nullptr), fVerbose(false)
58  {
59  for (Sample s: fSamples) {
60  fNBins += s.GetBinning().NBins();
61  fNBinsFull += s.GetBinning().NBins() * GetComponents(s).size();
62  fNBinsPerSample.push_back(s.GetBinning().NBins());
63  fComps.push_back(GetComponents(s));
64  }
65 
66  // Initialise index map
67  fIndexMap.resize(fNBinsFull);
68  unsigned int fullOffset(0), scaledOffset(0);
69  for (size_t iS = 0; iS < fSamples.size(); ++iS) { // loop over samples
70  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
71  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) { // loop over bins
72  unsigned int iScaled = scaledOffset + iB;
73  unsigned int iFull = fullOffset + (iC*fNBinsPerSample[iS]) + iB;
74  fIndexMap(iFull) = iScaled;
75  } // for bin i
76  } // for component i
77  scaledOffset += fNBinsPerSample[iS];
78  fullOffset += fComps[iS].size() * fNBinsPerSample[iS];
79  } // for sample i
80 
81  // Fill the data vector
84  size_t i = 0;
85  for (size_t iS = 0; iS < fSamples.size(); ++iS) {
86  double pot = fSamples[iS].GetData()->POT();
87  double livetime = fSamples[iS].GetData()->Livetime();
88  TH1D* hData = fSamples[iS].GetData()->ToTH1(pot);
89  TH1D* hCos = fSamples[iS].HasCosmic() ?
90  fSamples[iS].GetCosmic()->ToTH1(livetime, kLivetime) : nullptr;
91  for (size_t iBin = 1; iBin <= fNBinsPerSample[iS]; ++iBin) {
92  fData(i) = hData->GetBinContent(iBin);
93  fCosmic(i) = hCos ? hCos->GetBinContent(iBin) : 0;
94  ++i;
95  } // for bin i
96  delete hData;
97  delete hCos;
98  } // for sample iS
99 
101  for (size_t i = 0; i < fNBinsFull; ++i)
103 
104  // Set up beta and mu
105  fMu.resize(fNBinsFull);
106  fBeta = ArrayXd::Constant(fNBinsFull, 1);
107  fBetaMask.resize(fNBinsFull);
109  fSpecComponentwise = ArrayXd::Zero(fNBinsFull);
110 
111  if (fCovMx) {
112  // Invert matrix ahead of time
113  MatrixXd mx = CovarianceMatrix::ForcePosDef(fCovMx->GetFullCovMx(), epsilon);
114  fMxInv = mx.inverse();
115  MatrixXd id = mx;
116  id.setIdentity();
117  fResidual = ((mx.transpose() * fMxInv) - id).maxCoeff();
118  if (fVerbose) cout << "Eigen matrix inversion residual is " << fResidual << endl;
119  }
120 
121  }
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()
TH1D * fBetaHist
Levenberg-Marquardt lambda.
Var Constant(double c)
Use to weight events up and down by some factor.
Definition: Var.cxx:24
const XML_Char * s
Definition: expat.h:262
#define pot
double fLambdaZero
Whether to use LMA.
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
virtual ana::LikelihoodCovMxExperiment::~LikelihoodCovMxExperiment ( )
inlinevirtual

Definition at line 23 of file LikelihoodCovMxExperiment.h.

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

23 {};

Member Function Documentation

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

Reimplemented from ana::IExperiment.

Definition at line 124 of file LikelihoodCovMxExperiment.cxx.

References calc, om::cout, osc::DowncastToSterile(), allTimeWatchdog::endl, fComps, fCovMx, fData, fMu, fNBins, fPredUnshiftedHist, fSamples, fSpec, osc::IOscCalcSterile::GetAngle(), ana::covmx::CovarianceMatrix::GetBinning(), osc::IOscCalcSterile::GetDelta(), osc::IOscCalcSterile::GetDm(), GetExpectedSpectrum(), MECModelEnuComparisons::i, calib::j, LikelihoodCovMxNewton(), ana::LogLikelihood(), ana::Binning::NBins(), pot, and ana::Spectrum::ToTH1().

Referenced by ~LikelihoodCovMxExperiment().

126  {
127  DontAddDirectory guard;
128 
129  // Get MC and data spectra
130  size_t i = 0;
131  for (size_t iS = 0; iS < fSamples.size(); ++iS) {
132  double pot = fSamples[iS].GetData()->POT();
133  SystShifts shift = fSamples[i].GetSystShifts(syst);
134  for (Component comp: fComps[iS]) {
135  Spectrum spec = fSamples[iS].GetPrediction()->PredictComponentSyst(osc,
136  shift, get<0>(comp), get<1>(comp), get<2>(comp));
137  TH1D* hMu = spec.ToTH1(pot);
138  for (size_t iBin = 1; iBin <= (size_t)hMu->GetNbinsX(); ++iBin) {
139  fMu(i) = hMu->GetBinContent(iBin);
140  ++i;
141  } // for bin
142  delete hMu;
143  } // for component
144  } // for sample
145 
146  // Catch errors in prediction, to disambiguate from fitting problems
148  if (calc && fMu.array().isNaN().any()) {
149  cout << "ERROR! NAN VALUES IN PREDICTION!" << endl;
150  cout << "---------- OSCILLATION PARAMETERS ----------" << endl
151  << "Dm21: " << calc->GetDm(2) << endl
152  << "Dm31: " << calc->GetDm(3) << endl
153  << "Dm41: " << calc->GetDm(4) << endl
154  << "Theta 13: " << calc->GetAngle(1, 3) << endl
155  << "Theta 23: " << calc->GetAngle(2, 3) << endl
156  << "Theta 14: " << calc->GetAngle(1, 4) << endl
157  << "Theta 24: " << calc->GetAngle(2, 4) << endl
158  << "Theta 34: " << calc->GetAngle(3, 4) << endl
159  << "Delta 13: " << calc->GetDelta(1, 3) << endl
160  << "Delta 24: " << calc->GetDelta(2, 4) << endl
161  << "--------------------------------------------" << endl;
162  }
163 
165 
166  // Do stats-only calculation if no matrix
167  if (!fCovMx) return ana::LogLikelihood(fSpec, fData);
168 
169  // Number of bins
170  if (fCovMx && fNBins != (size_t)fCovMx->GetBinning().NBins()) {
171  throw runtime_error(
172  "Number of bins in predictions does not match covariance matrix!");
173  }
174 
175  // Fill histogram before beta shifting
176  if (fPredUnshiftedHist) {
177  for (size_t j = 0; j < fNBins; ++j) {
178  fPredUnshiftedHist->SetBinContent(j+1, fSpec(j));
179  }
180  } // if filling unshifted spectrum
181 
182  return LikelihoodCovMxNewton();
183 
184  } // function LikelihoodCovMxExperiment::ChiSq
std::vector< std::vector< covmx::Component > > fComps
const IOscCalcSterile * DowncastToSterile(const IOscCalc *calc, bool quiet)
std::vector< double > Spectrum
Definition: Constants.h:746
osc::OscCalcDumb calc
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
std::tuple< Flavors::Flavors_t, Current::Current_t, Sign::Sign_t > Component
virtual double GetAngle(int i, int j) const =0
virtual double GetDm(int i) const =0
#define pot
const double j
Definition: BetheBloch.cxx:29
base class for sterile oscillation calculators In the context of a sterile oscillation calculator...
OStream cout
Definition: OStream.cxx:6
int NBins() const
Definition: Binning.h:29
std::vector< covmx::Sample > fSamples
virtual double GetDelta(int i, int j) const =0
void ana::LikelihoodCovMxExperiment::DecreaseLambda ( ) const
protected

Definition at line 558 of file LikelihoodCovMxExperiment.cxx.

References fLambda, and fNu.

Referenced by LikelihoodCovMxNewton(), and SetVerbose().

558  {
559  fLambda /= fNu;
560  }
double fLambda
Levenberg-Marquardt nu.
double fNu
Levenberg-Marquardt starting lambda.
void ana::LikelihoodCovMxExperiment::EnableLMA ( ) const
protected
TH1D* ana::LikelihoodCovMxExperiment::GetBetaHist ( )
inline

Definition at line 37 of file LikelihoodCovMxExperiment.h.

References fBetaHist.

37 { return fBetaHist; };
TH1D * fBetaHist
Levenberg-Marquardt lambda.
TH1D* ana::LikelihoodCovMxExperiment::GetBetaMuHist ( )
inline

Definition at line 39 of file LikelihoodCovMxExperiment.h.

References fBetaMuHist.

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

Definition at line 207 of file LikelihoodCovMxExperiment.cxx.

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

Referenced by LikelihoodCovMxNewton(), and SetVerbose().

207  {
208 
209  // the statistical likelihood of the shifted spectrum...
211 
212  // ...plus the penalty the prediction picks up from its covariance
213  VectorXd vec = fBeta - 1.;
214  fSystChiSq = vec.transpose() * fMxInv * vec;
215  return fStatChiSq + fSystChiSq;
216 
217  } // 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
TH1D* ana::LikelihoodCovMxExperiment::GetDataHist ( )
inline

Definition at line 42 of file LikelihoodCovMxExperiment.h.

References fDataHist.

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

Definition at line 193 of file LikelihoodCovMxExperiment.cxx.

References fBeta, fBetaMu, fCosmic, fIndexMap, fMu, fNBinsFull, fSpec, fSpecComponentwise, and MECModelEnuComparisons::i.

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

193  {
194 
195  fBetaMu = fBeta * fMu;
196  fSpec.setZero();
197  for (size_t i = 0; i < fNBinsFull; ++i)
198  fSpec(fIndexMap(i)) += fBetaMu(i);
199  fSpec += fCosmic;
200  fSpecComponentwise.setZero();
201  for (size_t i = 0; i < fNBinsFull; ++i)
203 
204  } // function LikelihoodCovMxExperiment::GetExpectedSpectrum
void ana::LikelihoodCovMxExperiment::GetGradAndHess ( ) const
protected

Definition at line 265 of file LikelihoodCovMxExperiment.cxx.

References dg::corr, confusionMatrixTree::count, e, fBeta, fBetaMask, fDataComponentwise, fGrad, fHess, fIndexMap, fLambda, fMu, fMxInv, fNBinsFull, fSpecComponentwise, fUseLMA, MECModelEnuComparisons::i, calib::j, confusionMatrixTree::matrix, cet::sqlite::select(), tmp, submit_syst::y, and test::z.

Referenced by LikelihoodCovMxNewton(), and SetVerbose().

265  {
266 
267  unsigned int maskedSize = count(fBetaMask.begin(), fBetaMask.end(), true);
268  if (maskedSize != fGrad.rows()) {
269  fGrad.resize(maskedSize);
270  fHess.resize(maskedSize, maskedSize);
271  }
272  fHess.setZero();
273 
274  MatrixXd tmp = fMxInv * (fBeta-1).matrix().asDiagonal();
275  ArrayXd z = tmp.rowwise().sum();
276  ArrayXd y = (fSpecComponentwise<1e-40).select(1e-40,fSpecComponentwise);
277  ArrayXd corr = fMu * fDataComponentwise / (y*y);
278 
279  ArrayXd grad = -2 * (fMu * (1.0-(fDataComponentwise/y)) + z);
280 
281  size_t iCounter = 0;
282  for (size_t i = 0; i < fNBinsFull; ++i) {
283  if (!fBetaMask[i]) continue;
284  fGrad(iCounter) = grad(i);
285  size_t jCounter = iCounter;
286  for (size_t j = i; j < fNBinsFull; ++j) {
287  if (!fBetaMask[j]) continue;
288  fHess(iCounter, jCounter) = 2 * fMxInv(i, j);
289  if (fIndexMap(i) == fIndexMap(j)) {
290  fHess(iCounter, jCounter) += 2 * corr(i) * fMu(j);
291  }
292  ++jCounter;
293  }
294  ++iCounter;
295  } // for i
296 
297  fHess.transposeInPlace();
298 
299  // Apply transformation and Levenberg-Marquardt
300  if (fUseLMA) fHess.diagonal() *= 1. + fLambda;
301 
302  } // function LikelihoodCovMxExperiment::GetGradAndHess
double fLambda
Levenberg-Marquardt nu.
double corr
Float_t tmp
Definition: plot.C:36
const double j
Definition: BetheBloch.cxx:29
z
Definition: test.py:28
auto select(T const &...t)
Definition: select.h:146
Float_t e
Definition: plot.C:35
double ana::LikelihoodCovMxExperiment::GetIteration ( )
inline
TH1D* ana::LikelihoodCovMxExperiment::GetMuHist ( )
inline

Definition at line 38 of file LikelihoodCovMxExperiment.h.

References fMuHist.

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

Definition at line 40 of file LikelihoodCovMxExperiment.h.

References fPredShiftedHist.

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

Definition at line 41 of file LikelihoodCovMxExperiment.h.

References fPredUnshiftedHist.

double ana::LikelihoodCovMxExperiment::GetResidual ( )
inline

Definition at line 32 of file LikelihoodCovMxExperiment.h.

References fResidual.

double ana::LikelihoodCovMxExperiment::GetStatChiSq ( )
inline

Definition at line 30 of file LikelihoodCovMxExperiment.h.

References fStatChiSq.

double ana::LikelihoodCovMxExperiment::GetSystChiSq ( )
inline

Definition at line 31 of file LikelihoodCovMxExperiment.h.

References fSystChiSq.

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

Definition at line 562 of file LikelihoodCovMxExperiment.cxx.

References fLambda, fLambdaZero, fNu, and ResetLambda().

Referenced by LikelihoodCovMxNewton(), and SetVerbose().

562  {
563  fLambda *= fNu;
565  }
double fLambda
Levenberg-Marquardt nu.
double fNu
Levenberg-Marquardt starting lambda.
double fLambdaZero
Whether to use LMA.
void ana::LikelihoodCovMxExperiment::InitialiseBetas ( ) const
protected

Definition at line 220 of file LikelihoodCovMxExperiment.cxx.

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

Referenced by LikelihoodCovMxNewton(), and SetVerbose().

220  {
221  double eps(1e-40);
222  unsigned int fullOffset(0), scaledOffset(0);
223  for (size_t iS = 0; iS < fSamples.size(); ++iS) { // loop over samples
224  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) { // loop over bins
225  unsigned int iScaled = scaledOffset+iB;
226  double sum = fCosmic(iScaled);
227  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
228  unsigned int iFull = fullOffset+(iC*fNBinsPerSample[iS])+iB;
229  sum += fMu(iFull);
230  } // for component
231  double val = fData(iScaled) / (sum + eps);
232  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
233  unsigned int iFull = fullOffset+(iC*fNBinsPerSample[iS])+iB;
234  fBeta(iFull) = val;
235  } // for component
236  } // for bin
237  scaledOffset += fNBinsPerSample[iS]; // increase offset in scaled matrix
238  fullOffset += fComps[iS].size() * fNBinsPerSample[iS];
239  } // for sample
240 
241  } // 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
double ana::LikelihoodCovMxExperiment::LikelihoodCovMxNewton ( ) const
protected

Definition at line 305 of file LikelihoodCovMxExperiment.cxx.

References ana::assert(), febshutoff_auto::changes, confusionMatrixTree::count, hadd_many_files::counter, om::cout, DecreaseLambda(), e, EnableLMA(), febshutoff_auto::end, allTimeWatchdog::endl, gov::fnal::cd::rms::provider::equal(), fBeta, fBetaHist, fBetaMask, fBetaMu, fBetaMuHist, fGrad, fHess, fIteration, fLambda, fMu, fMuHist, fNBins, fNBinsFull, fPredShiftedHist, fSpec, fStatChiSq, fSystChiSq, fUseLMA, fVerbose, GetChiSq(), GetExpectedSpectrum(), GetGradAndHess(), MECModelEnuComparisons::i, IncreaseLambda(), InitialiseBetas(), calib::j, MaskBetas(), recentWatchdog::now, elec2geo::pos, prev(), ResetLambda(), and febshutoff_auto::start.

Referenced by ChiSq(), and SetVerbose().

305  {
306 
308 
309  // Start with all systematic shifts unmasked
310  for (size_t i = 0; i < fNBinsFull; ++i) fBetaMask[i] = true;
311 
312  // We're trying to solve for the best expectation in each bin and each component. A good seed
313  // value is 1 (no shift).
314  // Note: beta will contain the best fit systematic shifts at the end of the process.
315  // This should be saved to show true post-fit agreement
316  const int maxIters = 5e3;
317  fIteration = 0;
318  int nFailures = 0;
319  // InitialiseBetas();
320  MaskBetas();
321  fUseLMA = false;
322  ResetLambda(); // Initialise lambda at starting value
323  bool failed = false;
324 
325  // Keep a history of masks and chi squares to help avoid infinite loops
326  bool keepChangingMask = true;
327  vector<vector<bool>> maskHistory;
328  vector<double> chisqHistory;
329 
331  double prev = GetChiSq();
332  double minChi = prev;
333  if (fVerbose) {
334  cout << "Before minimization, chisq is " << prev
335  << " (" << fStatChiSq << " stat, " << fSystChiSq << " syst)" << endl;
336  }
337 
338  while (true) {
339 
340  ++fIteration;
341  if (!fUseLMA && fIteration > 10) EnableLMA();
342  if (fIteration > maxIters) {
343  cout << "No convergence after " << maxIters << " iterations! Quitting out of the infinite loop." << endl;
344  return minChi;
345  }
346 
347  // If we have negative betas, mask them out
348  GetGradAndHess();
349 
350  VectorXd deltaBeta = fHess.llt().solve(fGrad);
351  if (deltaBeta.array().isNaN().any()) {
352  int nanCounter = 0;
353  for (int i = 0; i < deltaBeta.size(); ++i) {
354  if (isnan(deltaBeta[i])) ++nanCounter;
355  }
356  if (fUseLMA) IncreaseLambda();
357  EnableLMA();
358  if (fVerbose) {
359  cout << "LLT solution failed! There are " << nanCounter << "nan delta betas. Resetting all betas.";
360  }
361  if (failed) InitialiseBetas();
362  failed = true;
363  continue;
364  }
365 
366  // if (fVerbose) {
367  // double relativeError = (fHessReduced * deltaBeta - fGradReduced).norm() / fGradReduced.norm();
368  // cout << "Relative error of linear solution is " << relativeError << endl;
369  // }
370 
371  size_t counter = 0;
372  for (size_t i = 0; i < fNBinsFull; ++i) {
373  if (fBetaMask[i]) {
374  fBeta(i) += deltaBeta(counter);
375  if (fBeta(i) < 0) { // if we pulled negative, mask this beta off
376  fBeta(i) = 0;
377  fBetaMask[i] = false;
378  } else if (fBeta(i) > 1e10) { // if we pull too high, mask this beta off
379  fBeta(i) = 1;
380  fBetaMask[i] = false;
381  }
382  ++counter;
383  }
384  } // for bin
385 
386  // Predict collapsed spectrum
388 
389  // Update the chisq
390  double chisq = GetChiSq();
391 
392  if (isnan(chisq) || chisq > 1e20) {
393  if (fVerbose && isnan(chisq))
394  cout << "ChiSq is NaN! Resetting minimisation." << endl;
395  else if (fVerbose && chisq > 1e20) {
396  cout << "ChiSq is anomalously large! Resetting minimisation." << endl;
397  vector<double> changes(deltaBeta.size());
398  VectorXd::Map(&changes[0], deltaBeta.size()) = deltaBeta;
399  std::sort(changes.begin(), changes.end(), std::greater<double>());
400  cout << " five greatest deltas:";
401  for (int i = 0; i < 5; ++i) cout << " " << changes[i];
402  cout << endl;
403  }
404  ++nFailures;
405  ResetLambda();
406  fLambda *= (1 + ((double)nFailures/10.)); // push us out of fragile failure states
407  EnableLMA();
408  if (failed) InitialiseBetas();
409  failed = true;
410  continue;
411  }
412 
413  if (chisq < minChi) minChi = chisq;
414 
415  // Bump down the LMA lambda
416  if (fIteration > 1 && fUseLMA && chisq < prev) DecreaseLambda();
417 
418  if (fVerbose) {
419  cout << "Iteration " << fIteration << ", chisq " << chisq
420  << " (" << fStatChiSq << " stat, " << fSystChiSq << " syst)" << endl;
421  }
422 
423  if (fSystChiSq < 0) {
424  assert(false && "Negative systematic penalty!");
425  }
426 
427  // If the updates didn't change anything at all then we're done
428  double change = fabs(chisq-prev);
429 
430  // Converge to third decimal place
431  if (change/chisq < 1e-3 || change < 1e-10) {
432 
433  bool maskChange = keepChangingMask ? MaskBetas() : false;
434 
435  if (maskChange) { // re-mask betas - if the mask changed, go again
436 
437  // this block of code prevents infinte loops thru bin masks
438  maskHistory.push_back(fBetaMask);
439  chisqHistory.push_back(chisq);
440  for (size_t i = 0; i < maskHistory.size()-1; ++i) {
441  if (std::equal(maskHistory[i].begin(), maskHistory[i].end(), fBetaMask.begin())) {
442  if (fVerbose) cout << "Found a repeated mask! Setting the mask that provided the best chisq" << endl;
443  double chi = 1e100;
444  size_t pos = 999;
445  for (size_t j = 0; j < chisqHistory.size(); ++j) {
446  if (chisqHistory[j] < chi) {
447  chi = chisqHistory[j];
448  pos = j;
449  }
450  }
451  fBetaMask = maskHistory[pos];
452  for (size_t i = 0; i < fNBinsFull; ++i) if (!fBetaMask[i]) fBeta[i] = 0;
453  keepChangingMask = false;
454  break;
455  }
456  }
457 
458  ResetLambda();
459  if (fVerbose) {
460  cout << "Minimisation round finished with chisq " << chisq << ", bin mask:";
461  for (size_t i = 0; i < fNBinsFull; ++i) if (!fBetaMask[i]) cout << " " << i+1;
462  cout << endl;
463  }
464  }
465 
466  // Converge to tenth significant figure or tenth decimal place, whichever is larger
467  // Alternatively, return smaller chisq if we're flip-flopping between two beta masks
468  if (!maskChange && (change/chisq < 1e-10 || change < 1e-10)) {
469 
470  if (fVerbose) {
471  cout << "Converged to " << chisq << " (" << fStatChiSq << " stat, " << fSystChiSq
472  << " syst) after " << fIteration << " iterations and ";
474  double secs = duration_cast<seconds>(end-start).count();
475  if (secs < 1) {
476  double millisecs = duration_cast<milliseconds>(end-start).count();
477  cout << millisecs << " ms." << endl;
478  } else if (secs < 60) cout << secs << " seconds." << endl;
479  else {
480  double mins = duration_cast<minutes>(end-start).count();
481  cout << mins << " minutes." << endl;
482  }
483  }
484 
485  // Fill histograms if necessary
486  if (fBetaHist) for (size_t i = 0; i < fNBinsFull; ++i) fBetaHist->SetBinContent(i+1, fBeta(i));
487  if (fMuHist) for (size_t i = 0; i < fNBinsFull; ++i) fMuHist->SetBinContent(i+1, fMu(i));
488  if (fBetaMuHist) for (size_t i = 0; i < fNBinsFull; ++i) fBetaMuHist->SetBinContent(i+1, fBetaMu(i));
489  if (fPredShiftedHist) for (size_t i = 0; i < fNBins; ++i) fPredShiftedHist->SetBinContent(i+1, fSpec(i));
490 
491  return chisq;
492 
493  } else {
494  if (fUseLMA) DecreaseLambda(); // If we're getting close to converging, speed things up
495  }
496  }
497 
498  prev = chisq;
499 
500  } // end while
501 
502  } // function LikelihoodCovMxExperiment::LogLikelihoodCovMxNewton
double fLambda
Levenberg-Marquardt nu.
bool equal(T *first, T *second)
TH1D * fBetaHist
Levenberg-Marquardt lambda.
void prev()
Definition: show_event.C:91
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
assert(nhit_max >=nhit_nbins)
Float_t e
Definition: plot.C:35
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().

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 244 of file LikelihoodCovMxExperiment.cxx.

References e, fBeta, fBetaMask, fBetaMu, fDataComponentwise, fMu, fMxInv, fNBinsFull, fSpecComponentwise, MECModelEnuComparisons::i, confusionMatrixTree::matrix, cet::sqlite::select(), tmp, submit_syst::y, and test::z.

Referenced by EnableLMA(), LikelihoodCovMxNewton(), and SetVerbose().

244  {
245 
246  MatrixXd tmp = fMxInv * (fBeta-1).matrix().asDiagonal();
247  tmp.diagonal() = -1 * fMxInv.diagonal();
248  ArrayXd z = tmp.rowwise().sum();
249  ArrayXd y = fSpecComponentwise - fBetaMu;
250  y = (y<1e-40).select(1e-40,y);
251  ArrayXd gradZero = 2 * (fMu * (1.0-(fDataComponentwise/y)) + z);
252 
253  Eigen::ArrayXi maskBefore(fNBinsFull);
254  for (size_t i = 0; i < fNBinsFull; ++i) maskBefore(i) = fBetaMask[i];
255  Eigen::ArrayXi maskAfter = (gradZero < 0).cast<int>();
256 
257  fBeta *= maskAfter.cast<double>();
258  for (size_t i = 0; i < fNBinsFull; ++i) fBetaMask[i] = maskAfter(i);
259 
260  return (maskBefore != maskAfter).any();
261 
262  } // function LikelihoodCovMxExperiment::MaskBetas
Float_t tmp
Definition: plot.C:36
z
Definition: test.py:28
auto select(T const &...t)
Definition: select.h:146
Float_t e
Definition: plot.C:35
void ana::LikelihoodCovMxExperiment::Reset ( ) const
overridevirtual

Reimplemented from ana::IExperiment.

Definition at line 187 of file LikelihoodCovMxExperiment.cxx.

References ana::Constant(), om::cout, allTimeWatchdog::endl, fBeta, and fNBinsFull.

Referenced by ~LikelihoodCovMxExperiment().

187  {
188  cout << "Resetting syst shifts due to new seed." << endl;
190  }
Var Constant(double c)
Use to weight events up and down by some factor.
Definition: Var.cxx:24
OStream cout
Definition: OStream.cxx:6
void ana::LikelihoodCovMxExperiment::ResetLambda ( ) const
protected

Definition at line 554 of file LikelihoodCovMxExperiment.cxx.

References fLambda, and fLambdaZero.

Referenced by IncreaseLambda(), LikelihoodCovMxNewton(), and SetVerbose().

554  {
556  }
double fLambda
Levenberg-Marquardt nu.
double fLambdaZero
Whether to use LMA.
void ana::LikelihoodCovMxExperiment::SaveHists ( bool  opt = true)

Definition at line 505 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().

505  {
506  if (!fCovMx) return;
507  // Start by clearing up existing beta histogram
508  if (fBetaHist) {
509  delete fBetaHist;
510  fBetaHist = nullptr;
511  }
512  if (fMuHist) {
513  delete fMuHist;
514  fMuHist = nullptr;
515  }
516  if (fBetaMuHist) {
517  delete fBetaMuHist;
518  fBetaMuHist = nullptr;
519  }
520  if (fPredShiftedHist) {
521  delete fPredShiftedHist;
522  fPredShiftedHist = nullptr;
523  }
524  if (fPredUnshiftedHist) {
525  delete fPredUnshiftedHist;
526  fPredUnshiftedHist = nullptr;
527  }
528  if (fDataHist) {
529  delete fDataHist;
530  fDataHist = nullptr;
531  }
532  // Now re-instantiate it if necessary
533  if (opt) {
534  Binning fullBins = fCovMx->GetFullBinning();
535  Binning scaledBins = fCovMx->GetBinning();
536  fBetaHist = new TH1D(UniqueName().c_str(), ";;", fullBins.NBins(), &fullBins.Edges()[0]);
537  fMuHist = new TH1D(UniqueName().c_str(), ";;", fullBins.NBins(), &fullBins.Edges()[0]);
538  fBetaMuHist = new TH1D(UniqueName().c_str(), ";;", fullBins.NBins(), &fullBins.Edges()[0]);
539  fPredShiftedHist = new TH1D(UniqueName().c_str(), ";;", scaledBins.NBins(), &scaledBins.Edges()[0]);
540  fPredUnshiftedHist = new TH1D(UniqueName().c_str(), ";;", scaledBins.NBins(), &scaledBins.Edges()[0]);
541  fDataHist = new TH1D(UniqueName().c_str(), ";;", scaledBins.NBins(), &scaledBins.Edges()[0]);
542  for (size_t i = 0; i < fNBins; ++i) {
543  fDataHist->SetBinContent(i+1, fData(i)); // Populate data hist
544  }
545  }
546  } // function LikelihoodCovMxExperiment::SaveHists
TH1D * fBetaHist
Levenberg-Marquardt lambda.
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:28
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::Reset().

33  {
34  assert(0 && "Not implemented");
35  }
assert(nhit_max >=nhit_nbins)
void ana::LikelihoodCovMxExperiment::SetVerbose ( bool  opt)
inline

Member Data Documentation

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

Levenberg-Marquardt lambda.

Definition at line 87 of file LikelihoodCovMxExperiment.h.

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

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

Definition at line 89 of file LikelihoodCovMxExperiment.h.

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

std::vector<std::vector<covmx::Component> > ana::LikelihoodCovMxExperiment::fComps
protected

Definition at line 80 of file LikelihoodCovMxExperiment.h.

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

Eigen::ArrayXd ana::LikelihoodCovMxExperiment::fCosmic
protected
covmx::CovarianceMatrix* ana::LikelihoodCovMxExperiment::fCovMx
protected

Definition at line 62 of file LikelihoodCovMxExperiment.h.

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

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

Definition at line 92 of file LikelihoodCovMxExperiment.h.

Referenced by GetDataHist(), and SaveHists().

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

Definition at line 74 of file LikelihoodCovMxExperiment.h.

Referenced by GetGradAndHess(), and LikelihoodCovMxNewton().

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

Definition at line 75 of file LikelihoodCovMxExperiment.h.

Referenced by GetGradAndHess(), and LikelihoodCovMxNewton().

Eigen::ArrayXi ana::LikelihoodCovMxExperiment::fIndexMap
protected
int ana::LikelihoodCovMxExperiment::fIteration
mutableprotected

Definition at line 99 of file LikelihoodCovMxExperiment.h.

Referenced by GetIteration(), and LikelihoodCovMxNewton().

double ana::LikelihoodCovMxExperiment::fLambda
mutableprotected

Levenberg-Marquardt nu.

Definition at line 85 of file LikelihoodCovMxExperiment.h.

Referenced by DecreaseLambda(), GetGradAndHess(), IncreaseLambda(), LikelihoodCovMxNewton(), and ResetLambda().

double ana::LikelihoodCovMxExperiment::fLambdaZero
protected

Whether to use LMA.

Definition at line 83 of file LikelihoodCovMxExperiment.h.

Referenced by IncreaseLambda(), and ResetLambda().

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

Definition at line 88 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

Definition at line 79 of file LikelihoodCovMxExperiment.h.

Referenced by InitialiseBetas(), and LikelihoodCovMxExperiment().

double ana::LikelihoodCovMxExperiment::fNu
protected

Levenberg-Marquardt starting lambda.

Definition at line 84 of file LikelihoodCovMxExperiment.h.

Referenced by DecreaseLambda(), and IncreaseLambda().

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

Definition at line 91 of file LikelihoodCovMxExperiment.h.

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

double ana::LikelihoodCovMxExperiment::fResidual
mutableprotected

Definition at line 98 of file LikelihoodCovMxExperiment.h.

Referenced by GetResidual(), and LikelihoodCovMxExperiment().

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

Definition at line 61 of file LikelihoodCovMxExperiment.h.

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

Eigen::ArrayXd ana::LikelihoodCovMxExperiment::fSpec
mutableprotected
Eigen::ArrayXd ana::LikelihoodCovMxExperiment::fSpecComponentwise
mutableprotected
double ana::LikelihoodCovMxExperiment::fStatChiSq
mutableprotected

Definition at line 96 of file LikelihoodCovMxExperiment.h.

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

double ana::LikelihoodCovMxExperiment::fSystChiSq
mutableprotected

Definition at line 97 of file LikelihoodCovMxExperiment.h.

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

bool ana::LikelihoodCovMxExperiment::fUseLMA
mutableprotected

Definition at line 82 of file LikelihoodCovMxExperiment.h.

Referenced by EnableLMA(), GetGradAndHess(), and LikelihoodCovMxNewton().

bool ana::LikelihoodCovMxExperiment::fVerbose
protected

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