LikelihoodCovMxExperiment.cxx
Go to the documentation of this file.
1 
3 
6 
7 #include "OscLib/IOscCalc.h"
8 
9 #include "TDirectory.h"
10 #include "TObjString.h"
11 #include "TH1.h"
12 #include "TMatrixD.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), 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  }
104 
105  //----------------------------------------------------------------------
107  const SystShifts& syst) const
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
153 
154  //---------------------------------------------------------------------------
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
180 
181  //---------------------------------------------------------------------------
182  double LikelihoodCovMxExperiment::GetChiSq(ArrayXd e) const {
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
193 
194  //---------------------------------------------------------------------------
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
217 
218  //---------------------------------------------------------------------------
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
266 
267  //---------------------------------------------------------------------------
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
322 
323  //---------------------------------------------------------------------------
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
358 
359  //---------------------------------------------------------------------------
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
507 
508  //----------------------------------------------------------------------
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
551 
552 }
double fLambda
Levenberg-Marquardt nu.
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
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
Eigen::MatrixXd GetFullCovMx()
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Double_t beta
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
LikelihoodCovMxExperiment(std::vector< covmx::Sample > samples, covmx::CovarianceMatrix *covmx, double epsilon=1e-5, double lambdazero=1, double nu=10)
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
const std::vector< double > & Edges() const
Definition: Binning.h:34
OStream cout
Definition: OStream.cxx:6
std::vector< unsigned int > fNBinsPerSample
Float_t norm
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
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
enum BeamMode string