LikelihoodCovMxExperiment.cxx
Go to the documentation of this file.
1 
3 
6 
7 #include "OscLib/IOscCalc.h"
8 #include "OscLib/OscCalcSterile.h" // for debugging
9 
10 #include "TDirectory.h"
11 #include "TObjString.h"
12 #include "TH1.h"
13 
14 #include <iomanip>
15 #include <iostream>
16 #include <chrono>
17 #include <cmath> // for isnan
18 
19 #include <Eigen/Dense>
20 
21 namespace ana
22 {
23  // Standard C++ namespace
24  using std::count;
25  using std::cout;
26  using std::endl;
27  using std::fixed;
28  using std::get;
29  using std::isnan;
30  using std::runtime_error;
31  using std::setprecision;
32  using std::string;
33  using std::vector;
34 
35  // Standard C++ chrono namespace
36  using std::chrono::high_resolution_clock;
37  using std::chrono::duration_cast;
38  using std::chrono::milliseconds;
39  using std::chrono::seconds;
40  using std::chrono::minutes;
41 
42  // Eigen tools
43  using Eigen::ArrayXd;
44  using Eigen::MatrixXd;
45  using Eigen::VectorXd;
46 
47  // CAFAna covariance matrix namespace
48  using covmx::Component;
49  using covmx::Sample;
50  using covmx::CovarianceMatrix;
51 
52  //----------------------------------------------------------------------
54  CovarianceMatrix* covmx, double epsilon, double lambdazero, double nu)
55  : fSamples(samples), fCovMx(covmx), fNBins(0), fNBinsFull(0), fLambdaZero(lambdazero),
56  fNu(nu), fBetaHist(nullptr), fMuHist(nullptr),
57  fBetaMuHist(nullptr), fPredShiftedHist(nullptr), fPredUnshiftedHist(nullptr),
58  fDataHist(nullptr), fResetShifts(true), 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  fGrad.resize(fNBinsFull);
68  fHess.resize(fNBinsFull, fNBinsFull);
69 
70  // Set up beta and mu
71  fMu.resize(fNBinsFull);
72  fBeta.resize(fNBinsFull);
73  fBetaMask.resize(fNBinsFull);
74  for (size_t i = 0; i < fNBinsFull; ++i) {
75  fBeta(i) = 1.;
76  }
77 
78  // Fill the data vector
81  size_t i = 0;
82  for (size_t iS = 0; iS < fSamples.size(); ++iS) {
83  double pot = fSamples[iS].GetData()->POT();
84  double livetime = fSamples[iS].GetData()->Livetime();
85  TH1D* hData = fSamples[iS].GetData()->ToTH1(pot);
86  TH1D* hCos = fSamples[iS].HasCosmic() ?
87  fSamples[iS].GetCosmic()->ToTH1(livetime, kLivetime) : nullptr;
88  for (size_t iBin = 1; iBin <= fNBinsPerSample[iS]; ++iBin) {
89  fData(i) = hData->GetBinContent(iBin);
90  fCosmic(i) = hCos ? hCos->GetBinContent(iBin) : 0;
91  ++i;
92  } // for bin i
93  delete hData;
94  delete hCos;
95  } // for sample iS
96 
97  if (fCovMx) {
98  // Invert matrix ahead of time
99  MatrixXd mx = CovarianceMatrix::ForcePosDef(fCovMx->GetFullCovMx(), epsilon);
100  fMxInv = mx.inverse();
101  MatrixXd id = mx;
102  id.setIdentity();
103  fResidual = ((mx.transpose() * fMxInv) - id).maxCoeff();
104  if (fVerbose) cout << "Eigen matrix inversion residual is " << fResidual << endl;
105  }
106 
107  }
108 
109  //----------------------------------------------------------------------
111  const SystShifts& syst) const
112  {
113  DontAddDirectory guard;
114 
115  // Get MC and data spectra
116  size_t i = 0;
117  for (size_t iS = 0; iS < fSamples.size(); ++iS) {
118  double pot = fSamples[iS].GetData()->POT();
119  SystShifts shift = fSamples[i].GetSystShifts(syst);
120  for (Component comp: fComps[iS]) {
121  Spectrum spec = fSamples[iS].GetPrediction()->PredictComponentSyst(osc,
122  shift, get<0>(comp), get<1>(comp), get<2>(comp));
123  TH1D* hMu = spec.ToTH1(pot);
124  for (size_t iBin = 1; iBin <= (size_t)hMu->GetNbinsX(); ++iBin) {
125  fMu(i) = hMu->GetBinContent(iBin);
126  ++i;
127  } // for bin
128  delete hMu;
129  } // for component
130  } // for sample
131 
132  // Catch errors in prediction, to disambiguate from fitting problems
133  osc::OscCalcSterile* calc = dynamic_cast<osc::OscCalcSterile*>(osc);
134  if (calc && fMu.array().isNaN().any()) {
135  cout << "ERROR! NAN VALUES IN PREDICTION!" << endl;
136  cout << "---------- OSCILLATION PARAMETERS ----------" << endl
137  << "Dm21: " << calc->GetDm(2) << endl
138  << "Dm31: " << calc->GetDm(3) << endl
139  << "Dm41: " << calc->GetDm(4) << endl
140  << "Theta 13: " << calc->GetAngle(1, 3) << endl
141  << "Theta 23: " << calc->GetAngle(2, 3) << endl
142  << "Theta 14: " << calc->GetAngle(1, 4) << endl
143  << "Theta 24: " << calc->GetAngle(2, 4) << endl
144  << "Theta 34: " << calc->GetAngle(3, 4) << endl
145  << "Delta 13: " << calc->GetDelta(1, 3) << endl
146  << "Delta 24: " << calc->GetDelta(2, 4) << endl
147  << "--------------------------------------------" << endl;
148  }
149 
150  if (fResetShifts) {
151  for (size_t i = 0; i < fNBinsFull; ++i) fBeta(i) = 1.;
152  }
153 
154  // Do stats-only calculation if no matrix
155  if (!fCovMx) {
156  ArrayXd e = GetExpectedSpectrum();
157  return ana::LogLikelihood(e, fData);
158  }
159 
160  // Number of bins
161  if (fCovMx && fNBins != (size_t)fCovMx->GetBinning().NBins()) {
162  throw runtime_error(
163  "Number of bins in predictions does not match covariance matrix!");
164  }
165 
166  // Fill histogram before beta shifting
167  if (fPredUnshiftedHist) {
168  ArrayXd e = GetExpectedSpectrum();
169  for (size_t j = 0; j < fNBins; ++j) {
170  fPredUnshiftedHist->SetBinContent(j+1, e(j));
171  }
172  } // if filling unshifted spectrum
173 
174  return LikelihoodCovMxNewton();
175 
176  } // function LikelihoodCovMxExperiment::ChiSq
177 
178  //---------------------------------------------------------------------------
180 
181  ArrayXd ret = ArrayXd::Zero(fNBins);
182  unsigned int fullOffset(0), scaledOffset(0);
183 
184  for (size_t iS = 0; iS < fSamples.size(); ++iS) { // loop over samples
185  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
186  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) { // loop over bins
187  unsigned int iScaled = scaledOffset + iB;
188  unsigned int iFull = fullOffset + (iC*fNBinsPerSample[iS]) + iB;
189 
190  ret(iScaled) += fMu(iFull)*fBeta(iFull);
191  } // for bin
192  } // for component
193  // Now loop one more time to add cosmics
194  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) {
195  unsigned int iScaled = scaledOffset + iB;
196  ret(iScaled) += fCosmic(iScaled);
197  } // for bin
198  scaledOffset += fNBinsPerSample[iS]; // increase offset in scaled matrix
199  fullOffset += fComps[iS].size() * fNBinsPerSample[iS];
200  } // for sample
201  return ret;
202 
203  } // function LikelihoodCovMxExperiment::GetExpectedSpectrum
204 
205  //---------------------------------------------------------------------------
206  double LikelihoodCovMxExperiment::GetChiSq(ArrayXd e) const {
207 
208  // the statistical likelihood of the shifted spectrum...
210 
211  // ...plus the penalty the prediction picks up from its covariance
212  VectorXd vec = fBeta - 1.;
213  fSystChiSq = vec.transpose() * fMxInv * vec;
214  return fStatChiSq + fSystChiSq;
215 
216  } // function LikelihoodCovMxExperiment::GetChiSq
217 
218  //---------------------------------------------------------------------------
220  double eps(1e-40);
221  unsigned int fullOffset(0), scaledOffset(0);
222  for (size_t iS = 0; iS < fSamples.size(); ++iS) { // loop over samples
223  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) { // loop over bins
224  unsigned int iScaled = scaledOffset+iB;
225  double sum = fCosmic(iScaled);
226  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
227  unsigned int iFull = fullOffset+(iC*fNBinsPerSample[iS])+iB;
228  sum += fMu(iFull);
229  } // for component
230  double val = fData(iScaled) / (sum + eps);
231  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
232  unsigned int iFull = fullOffset+(iC*fNBinsPerSample[iS])+iB;
233  fBeta(iFull) = val;
234  } // for component
235  } // for bin
236  scaledOffset += fNBinsPerSample[iS]; // increase offset in scaled matrix
237  fullOffset += fComps[iS].size() * fNBinsPerSample[iS];
238  } // for sample
239 
240  } // function LikelihoodCovMxExperiment::InitialiseBetas
241 
242  //---------------------------------------------------------------------------
244 
245  bool maskChange = false;
246 
247  size_t iFullOffset(0.), iScaledOffset(0.);
248  for (size_t iS = 0; iS < fSamples.size(); ++iS) { // loop over samples
249  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
250  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) { // loop over bins
251  // Get the current scaled & full bin
252  unsigned int iScaled = iScaledOffset+iB;
253  unsigned int iFull = iFullOffset+(iC*fNBinsPerSample[iS])+iB;
254 
255  double y = fCosmic(iScaled);
256  for (size_t jC = 0; jC < fComps[iS].size(); ++jC) {
257  if (iC != jC) { // beta mu from this bin is 0
258  unsigned int jFull = iFullOffset+(jC*fNBinsPerSample[iS])+iB;
259  y += fMu(jFull) * fBeta(jFull);
260  }
261  } // for component j
262  if (y < 1e-40) y = 1e-40; // clamp y
263 
264  double z(0.); // calculate z
265  for (size_t jFull = 0; jFull < fNBinsFull; ++jFull) {
266  double beta = (jFull == iFull) ? 0 : fBeta[jFull]; // set this beta to 0
267  z += (beta-1)*fMxInv(iFull, jFull);
268  }
269 
270  // Calculate gradient
271  double gradZero = 2 * (fMu(iFull)*(1.0-(fData(iScaled)/y)) + z);
272  if (gradZero > 0 && fBetaMask[iFull]) { // if we're masking a bin off
273  fBetaMask[iFull] = false;
274  maskChange = true;
275  }
276  if (gradZero <= 0 && !fBetaMask[iFull]) { // if we're masking a bin on
277  fBetaMask[iFull] = true;
278  fBeta(iFull) = 0;
279  maskChange = true;
280  }
281  } // for bin i
282  } // for component i
283  iScaledOffset += fNBinsPerSample[iS]; // increase offset in scaled matrix
284  iFullOffset += fComps[iS].size() * fNBinsPerSample[iS];
285  } // for sample i
286 
287  return maskChange;
288 
289  } // function LikelihoodCovMxExperiment::MaskBetas
290 
291  //---------------------------------------------------------------------------
293 
294  size_t iFullOffset(0.), iScaledOffset(0.);
295  for (size_t iS = 0; iS < fSamples.size(); ++iS) { // loop over samples
296  for (size_t iC = 0; iC < fComps[iS].size(); ++iC) { // loop over components
297  for (size_t iB = 0; iB < fNBinsPerSample[iS]; ++iB) { // loop over bins
298  // Get the current scaled & full bin
299  unsigned int iScaled = iScaledOffset+iB;
300  unsigned int iFull = iFullOffset+(iC*fNBinsPerSample[iS])+iB;
301  double y = fCosmic(iScaled); // calculate y
302  for (size_t jC = 0; jC < fComps[iS].size(); ++jC) {
303  unsigned int jFull = iFullOffset+(jC*fNBinsPerSample[iS])+iB;
304  y += fMu(jFull) * fBeta(jFull);
305  } // for component j
306  if (y < 1e-40) y = 1e-40; // clamp y
307  double z(0.); // calculate z
308  for (size_t jFull = 0; jFull < fNBinsFull; ++jFull) {
309  z += (fBeta(jFull)-1)*fMxInv(iFull, jFull);
310  }
311 
312  // Calculate gradient
313  fGrad(iFull) = 2 * (fMu(iFull)*(1.0-(fData(iScaled)/y)) + z);
314 
315  // Populate Hessian matrix
316  size_t jFullOffset(0.), jScaledOffset(0.);
317  for (size_t jS = 0; jS < fSamples.size(); ++jS) { // loop over samples
318  for (size_t jC = 0; jC < fComps[jS].size(); ++jC) { // loop over components
319  for (size_t jB = 0; jB < fNBinsPerSample[jS]; ++jB) { // loop over bins
320  // Get the current scaled & full bin
321  unsigned int jScaled = jScaledOffset + jB;
322  unsigned int jFull = jFullOffset+(jC*fNBinsPerSample[jS])+jB;
323  if (iScaled != jScaled) {
324  fHess(iFull, jFull) = 2 * fMxInv(iFull, jFull);
325  } else {
326  fHess(iFull, jFull) = 2 * (fMu(iFull)*fMu(jFull)*(fData(iScaled)/(y*y))
327  + fMxInv(iFull, jFull));
328  }
329  } // for bin j
330  } // for component j
331  jScaledOffset += fNBinsPerSample[jS]; // increase offset in scaled matrix
332  jFullOffset += fComps[jS].size() * fNBinsPerSample[jS];
333  } // for sample j
334  } // for bin i
335  } // for component i
336  iScaledOffset += fNBinsPerSample[iS]; // increase offset in scaled matrix
337  iFullOffset += fComps[iS].size() * fNBinsPerSample[iS];
338  } // for sample i
339 
340  // Apply transformation and Levenberg-Marquardt
341  for (size_t i = 0; i < fNBinsFull; ++i) {
342  if (fUseLMA) fHess(i, i) *= 1. + fLambda;
343  fGrad(i) *= -1;
344  }
345 
346  } // function LikelihoodCovMxExperiment::GetGradAndHess
347 
348  //---------------------------------------------------------------------------
350 
351  // Mask off bins
352  int maskedSize = count(fBetaMask.begin(), fBetaMask.end(), true);
353 
354  GetGradAndHess();
355 
356  if (fGradReduced.rows() != maskedSize) {
357  fGradReduced.resize(maskedSize);
358  fHessReduced.resize(maskedSize, maskedSize);
359  }
360 
361  size_t c = 0;
362  for (size_t i = 0; i < fNBinsFull; ++i) {
363  if (fBetaMask[i]) {
364  fGradReduced(c) = fGrad(i);
365  ++c;
366  }
367  }
368  size_t ci = 0;
369  for (size_t i = 0; i < fNBinsFull; ++i) {
370  if (fBetaMask[i]) {
371  size_t cj = 0;
372  for (size_t j = 0; j < fNBinsFull; ++j) {
373  if (fBetaMask[j]) {
374  fHessReduced(ci, cj) = fHess(i, j);
375  ++cj;
376  } // if j unmasked
377  } // for j
378  ++ci;
379  } // if i unmasked
380  } // for i
381 
382  } // function LikelihoodCovMxExperiment::GetReducedGradAndHess
383 
384  //---------------------------------------------------------------------------
386 
388 
389  // Start with all systematic shifts unmasked
390  for (size_t i = 0; i < fNBinsFull; ++i) fBetaMask[i] = true;
391 
392  // We're trying to solve for the best expectation in each bin and each component. A good seed
393  // value is 1 (no shift).
394  // Note: beta will contain the best fit systematic shifts at the end of the process.
395  // This should be saved to show true post-fit agreement
396  const int maxIters = 5e3;
397  fIteration = 0;
398  int nFailures = 0;
399  // InitialiseBetas();
400  MaskBetas();
401  fUseLMA = false;
402  ResetLambda(); // Initialise lambda at starting value
403  bool failed = false;
404 
405  // Keep a history of masks and chi squares to help avoid infinite loops
406  bool keepChangingMask = true;
407  vector<vector<bool>> maskHistory;
408  vector<double> chisqHistory;
409 
410  ArrayXd e = GetExpectedSpectrum();
411  double prev = GetChiSq(e);
412  double minChi = prev;
413  if (fVerbose) {
414  cout << "Before minimization, chisq is " << prev
415  << " (" << fStatChiSq << " stat, " << fSystChiSq << " syst)" << endl;
416  }
417 
418  while (true) {
419  ++fIteration;
420  if (!fUseLMA && fIteration > 10) EnableLMA();
421  if (fIteration > maxIters) {
422  cout << "No convergence after " << maxIters << " iterations! Quitting out of the infinite loop." << endl;
423  return minChi;
424  }
425 
426  // If we have negative betas, mask them out
428  VectorXd deltaBeta = fHessReduced.ldlt().solve(fGradReduced);
429  if (deltaBeta.array().isNaN().any()) {
430  int nanCounter = 0;
431  for (int i = 0; i < deltaBeta.size(); ++i) {
432  if (isnan(deltaBeta[i])) ++nanCounter;
433  }
434  if (fUseLMA) IncreaseLambda();
435  EnableLMA();
436  if (fVerbose) {
437  cout << "LDLT solution failed! There are " << nanCounter << "nan delta betas. Resetting all betas.";
438  }
439  if (failed) InitialiseBetas();
440  failed = true;
441  continue;
442  }
443 
444  // if (fVerbose) {
445  // double relativeError = (fHessReduced * deltaBeta - fGradReduced).norm() / fGradReduced.norm();
446  // cout << "Relative error of linear solution is " << relativeError << endl;
447  // }
448 
449  size_t counter = 0;
450  for (size_t i = 0; i < fNBinsFull; ++i) {
451  if (fBetaMask[i]) {
452  fBeta(i) += deltaBeta(counter);
453  if (fBeta(i) < 0) { // if we pulled negative, mask this beta off
454  fBeta(i) = 0;
455  fBetaMask[i] = false;
456  } else if (fBeta(i) > 1e10) { // if we pull too high, mask this beta off
457  fBeta(i) = 1;
458  fBetaMask[i] = false;
459  }
460  ++counter;
461  }
462  } // for bin
463 
464  // Predict collapsed spectrum
465  e = GetExpectedSpectrum();
466 
467  // Update the chisq
468  double chisq = GetChiSq(e);
469 
470  if (isnan(chisq) || chisq > 1e20) {
471  if (fVerbose && isnan(chisq))
472  cout << "ChiSq is NaN! Resetting minimisation." << endl;
473  else if (fVerbose && chisq > 1e20) {
474  cout << "ChiSq is anomalously large! Resetting minimisation." << endl;
475  vector<double> changes(deltaBeta.size());
476  VectorXd::Map(&changes[0], deltaBeta.size()) = deltaBeta;
477  std::sort(changes.begin(), changes.end(), std::greater<double>());
478  cout << " five greatest deltas:";
479  for (int i = 0; i < 5; ++i) cout << " " << changes[i];
480  cout << endl;
481  }
482  ++nFailures;
483  ResetLambda();
484  fLambda *= (1 + ((double)nFailures/10.)); // push us out of fragile failure states
485  EnableLMA();
486  if (failed) InitialiseBetas();
487  failed = true;
488  continue;
489  }
490 
491  if (chisq < minChi) minChi = chisq;
492 
493  // Bump down the LMA lambda
494  if (fIteration > 1 && fUseLMA && chisq < prev) DecreaseLambda();
495 
496  if (fVerbose) {
497  cout << "Iteration " << fIteration << ", chisq " << chisq
498  << " (" << fStatChiSq << " stat, " << fSystChiSq << " syst)" << endl;
499  }
500 
501  if (fSystChiSq < 0) {
502  assert(false && "Negative systematic penalty!");
503  }
504 
505  // If the updates didn't change anything at all then we're done
506  double change = fabs(chisq-prev);
507 
508  // Converge to third decimal place
509  if (change/chisq < 1e-3 || change < 1e-10) {
510 
511  bool maskChange = keepChangingMask ? MaskBetas() : false;
512 
513  if (maskChange) { // re-mask betas - if the mask changed, go again
514 
515  // this block of code prevents infinte loops thru bin masks
516  maskHistory.push_back(fBetaMask);
517  chisqHistory.push_back(chisq);
518  for (size_t i = 0; i < maskHistory.size()-1; ++i) {
519  if (std::equal(maskHistory[i].begin(), maskHistory[i].end(), fBetaMask.begin())) {
520  if (fVerbose) cout << "Found a repeated mask! Setting the mask that provided the best chisq" << endl;
521  double chi = 1e100;
522  size_t pos = 999;
523  for (size_t j = 0; j < chisqHistory.size(); ++j) {
524  if (chisqHistory[j] < chi) {
525  chi = chisqHistory[j];
526  pos = j;
527  }
528  }
529  fBetaMask = maskHistory[pos];
530  for (size_t i = 0; i < fNBinsFull; ++i) if (!fBetaMask[i]) fBeta[i] = 0;
531  keepChangingMask = false;
532  break;
533  }
534  }
535 
536  ResetLambda();
537  if (fVerbose) {
538  cout << "Minimisation round finished with chisq " << chisq << ", bin mask:";
539  for (size_t i = 0; i < fNBinsFull; ++i) if (!fBetaMask[i]) cout << " " << i+1;
540  cout << endl;
541  }
542  }
543 
544  // Converge to tenth significant figure or tenth decimal place, whichever is larger
545  // Alternatively, return smaller chisq if we're flip-flopping between two beta masks
546  if (!maskChange && (change/chisq < 1e-10 || change < 1e-10)) {
547 
548  if (fVerbose) {
549  cout << "Converged to " << chisq << " (" << fStatChiSq << " stat, " << fSystChiSq
550  << " syst) after " << fIteration << " iterations and ";
552  double secs = duration_cast<seconds>(end-start).count();
553  if (secs < 1) {
554  double millisecs = duration_cast<milliseconds>(end-start).count();
555  cout << millisecs << " ms." << endl;
556  } else if (secs < 60) cout << secs << " seconds." << endl;
557  else {
558  double mins = duration_cast<minutes>(end-start).count();
559  cout << mins << " minutes." << endl;
560  }
561  }
562 
563  // Fill histograms if necessary
564  if (fBetaHist) for (size_t i = 0; i < fNBinsFull; ++i) fBetaHist->SetBinContent(i+1, fBeta(i));
565  if (fMuHist) for (size_t i = 0; i < fNBinsFull; ++i) fMuHist->SetBinContent(i+1, fMu(i));
566  if (fBetaMuHist) for (size_t i = 0; i < fNBinsFull; ++i) fBetaMuHist->SetBinContent(i+1, fBeta(i)*fMu(i));
567  if (fPredShiftedHist) for (size_t i = 0; i < fNBins; ++i) fPredShiftedHist->SetBinContent(i+1, e[i]);
568 
569  return chisq;
570 
571  } else {
572  if (fUseLMA) DecreaseLambda(); // If we're getting close to converging, speed things up
573  }
574  }
575 
576  prev = chisq;
577 
578  } // end while
579 
580  } // function LikelihoodCovMxExperiment::LogLikelihoodCovMxNewton
581 
582  //----------------------------------------------------------------------
584  if (!fCovMx) return;
585  // Start by clearing up existing beta histogram
586  if (fBetaHist) {
587  delete fBetaHist;
588  fBetaHist = nullptr;
589  }
590  if (fMuHist) {
591  delete fMuHist;
592  fMuHist = nullptr;
593  }
594  if (fBetaMuHist) {
595  delete fBetaMuHist;
596  fBetaMuHist = nullptr;
597  }
598  if (fPredShiftedHist) {
599  delete fPredShiftedHist;
600  fPredShiftedHist = nullptr;
601  }
602  if (fPredUnshiftedHist) {
603  delete fPredUnshiftedHist;
604  fPredUnshiftedHist = nullptr;
605  }
606  if (fDataHist) {
607  delete fDataHist;
608  fDataHist = nullptr;
609  }
610  // Now re-instantiate it if necessary
611  if (opt) {
612  Binning fullBins = fCovMx->GetFullBinning();
613  Binning scaledBins = fCovMx->GetBinning();
614  fBetaHist = new TH1D(UniqueName().c_str(), ";;", fullBins.NBins(), &fullBins.Edges()[0]);
615  fMuHist = new TH1D(UniqueName().c_str(), ";;", fullBins.NBins(), &fullBins.Edges()[0]);
616  fBetaMuHist = new TH1D(UniqueName().c_str(), ";;", fullBins.NBins(), &fullBins.Edges()[0]);
617  fPredShiftedHist = new TH1D(UniqueName().c_str(), ";;", scaledBins.NBins(), &scaledBins.Edges()[0]);
618  fPredUnshiftedHist = new TH1D(UniqueName().c_str(), ";;", scaledBins.NBins(), &scaledBins.Edges()[0]);
619  fDataHist = new TH1D(UniqueName().c_str(), ";;", scaledBins.NBins(), &scaledBins.Edges()[0]);
620  for (size_t i = 0; i < fNBins; ++i) {
621  fDataHist->SetBinContent(i+1, fData(i)); // Populate data hist
622  }
623  }
624  } // function LikelihoodCovMxExperiment::SaveHists
625 
627  fUseLMA = true;
628  for (size_t i = 0; i < fNBinsFull; ++i) fBeta(i) = 1.;
629  MaskBetas();
630  }
631 
634  }
635 
637  fLambda /= fNu;
638  }
639 
641  fLambda *= fNu;
643  }
644 
645 }
double fLambda
Levenberg-Marquardt nu.
double GetDelta(int i, int j) const
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< std::vector< covmx::Component > > fComps
LikelihoodCovMxExperiment(std::vector< covmx::Sample > samples, covmx::CovarianceMatrix *covmx, double epsilon=1e-5, double lambdazero=0, double nu=10)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::vector< std::tuple< ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t > > GetComponents()
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
double fNu
Levenberg-Marquardt starting lambda.
Eigen::MatrixXd GetFullCovMx()
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
bool equal(T *first, T *second)
Adapt the PMNS_Sterile calculator to standard interface.
Double_t beta
TH1D * fBetaHist
Levenberg-Marquardt lambda.
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
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const XML_Char * s
Definition: expat.h:262
double GetChiSq(Eigen::ArrayXd e) const
void prev()
Definition: show_event.C:91
#define pot
const double j
Definition: BetheBloch.cxx:29
Eigen::VectorXd vec
z
Definition: test.py:28
Oscillation probability calculators.
Definition: Calcs.h:5
double fLambdaZero
Whether to use LMA.
const std::vector< double > & Edges() const
Definition: Binning.h:34
OStream cout
Definition: OStream.cxx:6
std::vector< unsigned int > fNBinsPerSample
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const
void Zero()
double livetime
Definition: saveFDMCHists.C:21
int NBins() const
Definition: Binning.h:29
double epsilon
std::vector< covmx::Sample > fSamples
assert(nhit_max >=nhit_nbins)
double GetDm(int i) const
Double_t sum
Definition: plot.C:31
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
double GetAngle(int i, int j) const
enum BeamMode string