Classes | Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
ana::CovMxExperiment Class Reference

Compare a data spectrum to MC using a covariance matrix in the chi^2. More...

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

Inheritance diagram for ana::CovMxExperiment:
ana::SingleSampleExperiment ana::IExperiment

Classes

struct  Component
 Helper class used to specify components that are separately treated with their own covariance matrices. More...
 

Public Member Functions

 CovMxExperiment (const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, const std::map< Component, Eigen::MatrixXd > &covmx, double cosmicScaleError=0)
 
 CovMxExperiment (const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, const Eigen::MatrixXd &covmx, double cosmicScaleError=0)
 Simplified version for use with a single covariance matrix. More...
 
 CovMxExperiment (const IPrediction *pred, const Spectrum &data, const std::map< Component, Eigen::MatrixXd > &covmx)
 No-cosmics version (see first constructor for explanation of parameters) More...
 
 CovMxExperiment (const IPrediction *pred, const Spectrum &data, const Eigen::MatrixXd &covmx)
 Single-matrix variant of no-cosmics constructor. More...
 
virtual double ChiSq (osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
 
 CovMxExperiment (const SingleSampleExperiment &)=delete
 
CovMxExperimentoperator= (const SingleSampleExperiment &)=delete
 
const IPredictionGetPrediction () const
 
virtual stan::math::var LogLikelihood (osc::IOscCalcAdjustableStan *osc, const SystShifts &syst=kNoShift) const override
 
virtual void SaveTo (TDirectory *dir, const std::string &name) const override
 

Static Public Member Functions

static std::unique_ptr< SingleSampleExperimentLoadFrom (TDirectory *dir, const std::string &name)
 

Protected Member Functions

void MakeInvCovMxs (const std::map< Component, Eigen::MatrixXd > &covmx)
 
template<typename T >
Eigen::Array< T, Eigen::Dynamic, 1 > PredHistIncCosmics (osc::_IOscCalc< T > *calc, const SystShifts &syst) const
 

Protected Attributes

std::map< Component, Eigen::MatrixXd > fCovMxInv
 
const IPredictionfMC
 
Spectrum fData
 
std::optional< SpectrumfCosmic
 
double fCosmicScaleError
 
bool fPoissonError
 
std::vector< std::pair< double, double > > fBinErrors
 

Detailed Description

Compare a data spectrum to MC using a covariance matrix in the chi^2.

Definition at line 24 of file CovMxExperiment.h.

Constructor & Destructor Documentation

ana::CovMxExperiment::CovMxExperiment ( const IPrediction pred,
const Spectrum data,
const Spectrum cosmic,
const std::map< Component, Eigen::MatrixXd > &  covmx,
double  cosmicScaleError = 0 
)
inline
Parameters
predSource of oscillated MC beam predictions
dataData spectrum to compare to
cosmicCosmic ray background component
covmxBin-to-bin covariance matrix for each component you want to separately treat. See below for alternate constructor if you have just a single covariance matrix (or pass a map with the single key Component(Flavors::kAll,Current::kBoth,Sign::kBoth)).
cosmicScaleErrorfractional uncertainty on cosmic normalization

Definition at line 51 of file CovMxExperiment.h.

References MakeInvCovMxs().

Referenced by CovMxExperiment().

56  : SingleSampleExperiment(pred, data, cosmic, cosmicScaleError)
57  {
58  this->MakeInvCovMxs(covmx);
59  };
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0, bool PoissonError=false)
void MakeInvCovMxs(const std::map< Component, Eigen::MatrixXd > &covmx)
const XML_Char const XML_Char * data
Definition: expat.h:268
Sum up livetimes from individual cosmic triggers.
ana::CovMxExperiment::CovMxExperiment ( const IPrediction pred,
const Spectrum data,
const Spectrum cosmic,
const Eigen::MatrixXd &  covmx,
double  cosmicScaleError = 0 
)
inline

Simplified version for use with a single covariance matrix.

Definition at line 62 of file CovMxExperiment.h.

References ana::Flavors::kAll, ana::Current::kBoth, ana::Sign::kBoth, and MakeInvCovMxs().

67  : SingleSampleExperiment(pred, data, cosmic, cosmicScaleError)
68  {
70  };
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0, bool PoissonError=false)
void MakeInvCovMxs(const std::map< Component, Eigen::MatrixXd > &covmx)
const XML_Char const XML_Char * data
Definition: expat.h:268
Interactions of both types.
Definition: IPrediction.h:42
Sum up livetimes from individual cosmic triggers.
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
All neutrinos, any flavor.
Definition: IPrediction.h:26
ana::CovMxExperiment::CovMxExperiment ( const IPrediction pred,
const Spectrum data,
const std::map< Component, Eigen::MatrixXd > &  covmx 
)
inline

No-cosmics version (see first constructor for explanation of parameters)

Definition at line 73 of file CovMxExperiment.h.

References MakeInvCovMxs().

77  {
78  this->MakeInvCovMxs(covmx);
79  };
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0, bool PoissonError=false)
void MakeInvCovMxs(const std::map< Component, Eigen::MatrixXd > &covmx)
const XML_Char const XML_Char * data
Definition: expat.h:268
ana::CovMxExperiment::CovMxExperiment ( const IPrediction pred,
const Spectrum data,
const Eigen::MatrixXd &  covmx 
)
inline

Single-matrix variant of no-cosmics constructor.

Definition at line 82 of file CovMxExperiment.h.

References ChiSq(), CovMxExperiment(), ana::Flavors::kAll, ana::Current::kBoth, ana::Sign::kBoth, MakeInvCovMxs(), ana::SystShifts::Nominal(), and operator=().

86  {
88  };
SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, const Spectrum &cosmic, double cosmicScaleError=0, bool PoissonError=false)
void MakeInvCovMxs(const std::map< Component, Eigen::MatrixXd > &covmx)
const XML_Char const XML_Char * data
Definition: expat.h:268
Interactions of both types.
Definition: IPrediction.h:42
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
All neutrinos, any flavor.
Definition: IPrediction.h:26
ana::CovMxExperiment::CovMxExperiment ( const SingleSampleExperiment )
delete

Member Function Documentation

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

Reimplemented from ana::SingleSampleExperiment.

Definition at line 21 of file CovMxExperiment.cxx.

References bin, chi2(), ana::Chi2CovMx(), release_diff::diff, e, stan::math::fabs(), ana::SingleSampleExperiment::fCosmic, ana::SingleSampleExperiment::fCosmicScaleError, fCovMxInv, ana::SingleSampleExperiment::fData, ana::SingleSampleExperiment::fMC, ana::Spectrum::GetEigen(), ana::SystShifts::GetShift(), ana::SystShifts::IsNominal(), ana::Flavors::kAll, ana::Current::kBoth, ana::Sign::kBoth, ana::kCosmicBkgScaleSyst, pot, ana::Spectrum::POT(), ana::IPrediction::Predict(), ana::IPrediction::PredictComponent(), and scale.

Referenced by CovMxExperiment().

22  {
23  // the covariance matrix replaces the systematic shift approach.
24  if (!syst.IsNominal())
25  throw std::runtime_error("ExperimentWithCovMx::ChiSq(): Systematic shifts cannot be specified when fitting against covariance matrices");
26 
27  const double pot = fData.POT();
28  Eigen::ArrayXd hdata = fData.GetEigen(pot);
29 
30  std::map<Component, Eigen::ArrayXd> predictions;
31  for (const auto & mxPair : this->fCovMxInv)
32  {
33  const auto & comp = mxPair.first;
34  if (comp.flavor == Flavors::kAll && comp.current == Current::kBoth && comp.sign == Sign::kBoth)
35  {
36  // if you specified one matrix as applying to the whole prediction,
37  // why did you give me more of them??
38  if (this->fCovMxInv.size() > 1)
39  throw std::runtime_error("CovMxExperiment::ChiSq(): one component corresponds to the whole prediction but you also provided others?!");
40 
41  predictions.emplace(comp, fMC->Predict(calc).GetEigen(pot));
42  }
43  else
44  predictions.emplace(mxPair.first, fMC->PredictComponent(calc, comp.flavor, comp.current, comp.sign).GetEigen(pot));
45  }
46 
47  // will want the total prediction in a moment,
48  // but since we assume the components add to the total
49  // in the following calculations, might as well ensure it here.
50  Eigen::ArrayXd totalPred = fMC->Predict(calc).GetEigen(pot);
51  Eigen::ArrayXd diff = totalPred;
52  for (const auto & predPair : predictions)
53  diff -= predPair.second;
54  for (int bin = 0; bin < totalPred.size(); bin++)
55  {
56  if ( std::fabs(diff[bin]) > 1e-4 )
57  throw std::runtime_error("CovMxExperiment::ChiSq(): components don't sum to total prediction!");
58  }
59 
60  // add the cosmics into the prediction.
61  if(fCosmic)
62  {
63  if(fCosmicScaleError != 0)
64  {
65  const double scale = 1 + syst.GetShift(&kCosmicBkgScaleSyst) * fCosmicScaleError;
66  totalPred += fCosmic->GetEigen(pot) * scale;
67  }
68  else
69  totalPred += fCosmic->GetEigen(pot);
70 
71  }
72 
73  // first, compute the statistical uncertainty covariance matrix
74  // (which is diagonal). we'll need this repeatedly below.
75  // use the *total* prediction to get the sigma^2 here because
76  // that's the expected fluctuation size
77  Eigen::MatrixXd stat_mx(totalPred.size()-2, totalPred.size()-2);
78  for (int binNum = 1; binNum <= totalPred.size()-2; binNum++)
79  stat_mx(binNum-1, binNum-1) = totalPred[binNum];
80 
81 
82  // We need to include the statistical covariance matrix
83  // to do the chi^2 calculation or it'll be way off.
84  // we need to recalculate the sum for *every* new configuration of
85  // the IPrediction, unfortunately--
86  // and since we have to then invert it, this could be very costly.
87  // StackExchange to the rescue: somebody there wanted to solve a similar problem
88  // (where we already have the inverse of a matrix A
89  // and want to compute (A+B)^{-1] with some other invertible B):
90  // http://math.stackexchange.com/a/17780
91  // method has apparently been discovered in various forms by multiple people,
92  // but this version is due to K. Miller, Math.Mag. 54, 67 (1981)
93  //
94 
95 
96  // Strictly speaking, it's not possible to divide the prediction up
97  // into pieces and then use the covariance matrices for those pieces
98  // to compute a chi^2 because it's not possible to know which was which
99  // in the data. However, we'll cheat a little and do the following:
100  // for each component, we'll subtract the rest of the prediction
101  // and compare background-subtracted data to that component.
102  // Because we have to include the statistical covariance matrix every time,
103  // this over-counts the statistical uncertainty by a factor of
104  // roughly (N-1), where N is the number of components.
105  // Probably this isn't exactly the right statistical thing to do,
106  // but for a first approximation, we'll take an average of the chi^2s
107  // weighted by how many events they contribute to the prediction.
108 
109  double chi2 = 0;
110  for (auto it_comp = this->fCovMxInv.begin(); it_comp != this->fCovMxInv.end(); ++it_comp)
111  {
112  const auto & comp = it_comp->first;
113  const auto & compCovMxInv = it_comp->second;
114  const auto & predComp = (predictions.at(comp));
115  Eigen::ArrayXd hdata_compSub = hdata - ( totalPred - predComp );
116 
117 
118  // now recursively update the matrix
119  // by calculating the inverse of (A+B_k),
120  // where B_k is one of the N rank-1 matrices
121  // (when N is the rank of B) that B can be decomposed into.
122  // our case is even easier because we know B is diagonal
123  // (it's the statistical uncertainty covariance matrix).
124  Eigen::MatrixXd Bk(stat_mx.rows(), stat_mx.cols());
125  Eigen::MatrixXd CkInv = compCovMxInv; // here we're setting what Miller calls C_1 = A
126 
127  // the rank of the stat mx is just the number of rows:
128  // it's diagonal and the bins are orthogonal.
129  // note that we're counting from zero instead of 1
130  // for convenience in the TMatrixDs
131  for (int k = 0; k < stat_mx.rows(); k++)
132  {
133  // reset the last iteration so we can re-use
134  if (k > 0)
135  Bk(k-1, k-1) = 0;
136 
137  // set current iteration B_k
138  Bk(k, k) = stat_mx(k, k);
139 
140  // can short-circuit the trace in Miller's algorithm
141  // since our Bk come from a diagonal matrix:
142  // they're single entries on the diagonal
143  double gk = 1. / (1 + CkInv(k, k)*Bk(k, k));
144 
145  CkInv -= gk * CkInv * Bk * CkInv;
146  }
147 
148  // event-weighted chi^2
149  chi2 += predComp.sum() * Chi2CovMx(predictions.at(comp), hdata_compSub, CkInv);
150 
151  }
152 
153  // divide out the total number of events to the get the weighted chi^2
154  chi2 /= totalPred.sum();
155 
156  return chi2;
157  }
double Chi2CovMx(const Eigen::ArrayXd &e, const Eigen::ArrayXd &o, const Eigen::MatrixXd &covmxinv)
Definition: EigenUtils.cxx:52
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
std::optional< Spectrum > fCosmic
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:189
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
osc::OscCalcDumb calc
double chi2()
Double_t scale
Definition: plot.C:25
Interactions of both types.
Definition: IPrediction.h:42
std::map< Component, Eigen::MatrixXd > fCovMxInv
#define pot
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:227
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
const IPrediction* ana::SingleSampleExperiment::GetPrediction ( ) const
inlineinherited
std::unique_ptr< SingleSampleExperiment > ana::SingleSampleExperiment::LoadFrom ( TDirectory *  dir,
const std::string name 
)
staticinherited

Definition at line 154 of file SingleSampleExperiment.cxx.

References ana::assert(), dir, ana::Spectrum::LoadFrom(), ana::LoadFrom< IPrediction >(), mc, runNovaSAM::release, and runNovaSAM::ret.

155  {
156  dir = dir->GetDirectory(name.c_str()); // switch to subdir
157  assert(dir);
158 
159  TObjString* ptag = (TObjString*)dir->Get("type");
160  assert(ptag);
161  assert(ptag->GetString() == "SingleSampleExperiment");
162  delete ptag;
163 
164  const IPrediction* mc = ana::LoadFrom<IPrediction>(dir, "mc").release();
165  const std::unique_ptr<Spectrum> data = Spectrum::LoadFrom(dir, "data");
166 
167  auto ret = std::make_unique<SingleSampleExperiment>(mc, *data);
168 
169  if(dir->Get("cosmic")){
170  ret->fCosmic = *ana::LoadFrom<Spectrum>(dir, "cosmic");
171  }
172 
173  delete dir;
174 
175  return ret;
176  }
const XML_Char * name
Definition: expat.h:151
const XML_Char const XML_Char * data
Definition: expat.h:268
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:18
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
stan::math::var ana::SingleSampleExperiment::LogLikelihood ( osc::IOscCalcAdjustableStan osc,
const SystShifts syst = kNoShift 
) const
overridevirtualinherited

Reimplemented from ana::IExperiment.

Definition at line 122 of file SingleSampleExperiment.cxx.

References ana::SingleSampleExperiment::fData, ana::Spectrum::GetEigen(), ana::LogLikelihood(), ana::Spectrum::POT(), plot_validation_datamc::pred, and ana::SingleSampleExperiment::PredHistIncCosmics().

Referenced by test_stanfit_systpulls().

124  {
126  Eigen::ArrayXd data = fData.GetEigen(fData.POT());
127 
128  // fully-qualified so that we get the one in StanUtils.h
129  return ana::LogLikelihood(pred, data) / -2.; // LogLikelihood(), confusingly, returns chi2=-2*LL
130  }
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:189
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
const XML_Char const XML_Char * data
Definition: expat.h:268
double POT() const
Definition: Spectrum.h:227
Eigen::Array< stan::math::var, Eigen::Dynamic, 1 > ArrayXstan
Definition: StanUtils.h:7
Eigen::Array< T, Eigen::Dynamic, 1 > PredHistIncCosmics(osc::_IOscCalc< T > *calc, const SystShifts &syst) const
void ana::CovMxExperiment::MakeInvCovMxs ( const std::map< Component, Eigen::MatrixXd > &  covmx)
inlineprotected

Definition at line 98 of file CovMxExperiment.h.

References fCovMxInv, and ana::SymmMxInverse().

Referenced by CovMxExperiment().

99  {
100  for (const auto & mxPair : covmx)
101  this->fCovMxInv.emplace( mxPair.first, SymmMxInverse(mxPair.second) );
102 
103  }
std::map< Component, Eigen::MatrixXd > fCovMxInv
Eigen::MatrixXd SymmMxInverse(const Eigen::MatrixXd &mx)
Definition: EigenUtils.cxx:64
CovMxExperiment& ana::CovMxExperiment::operator= ( const SingleSampleExperiment )
delete

Referenced by CovMxExperiment().

template<typename T >
template Eigen::ArrayXstan ana::SingleSampleExperiment::PredHistIncCosmics ( osc::_IOscCalc< T > *  calc,
const SystShifts syst 
) const
protectedinherited

Definition at line 64 of file SingleSampleExperiment.cxx.

References calc, ana::SingleSampleExperiment::fBinErrors, ana::SingleSampleExperiment::fCosmic, ana::SingleSampleExperiment::fCosmicScaleError, ana::SingleSampleExperiment::fData, ana::SingleSampleExperiment::fMC, ana::SingleSampleExperiment::fPoissonError, ana::Spectrum::GetEigen(), ana::Spectrum::GetEigenStan(), ana::SystShifts::GetShift(), ana::Spectrum::HasStan(), MECModelEnuComparisons::i, ana::kLivetime, ana::Spectrum::Livetime(), ana::Spectrum::POT(), ana::IPrediction::PredictSyst(), ana::SystShifts::RemoveShift(), scale, sigma(), and T.

Referenced by ana::SingleSampleExperiment::ChiSq(), and ana::SingleSampleExperiment::LogLikelihood().

66  {
67  using EigenArray = Eigen::Array<T, Eigen::Dynamic, 1>;
68 
69  SystShifts systNoCosmic = syst;
70  systNoCosmic.RemoveShift(&kCosmicBkgScaleSyst);
71 
72  Spectrum specPred = fMC->PredictSyst(calc, systNoCosmic);
73  EigenArray apred;
74  if (specPred.HasStan())
75  {
76  if constexpr (std::is_same_v<T, stan::math::var>)
77  apred = specPred.GetEigenStan(fData.POT());
78  }
79  else
80  apred = specPred.GetEigen(fData.POT());
81 
82  if(fCosmic.has_value()){
83  // note the conversion here from double to stan::math::var if T is stan::math::var
84  EigenArray acosm = fCosmic->GetEigen(fData.Livetime(), kLivetime);
85  if(fCosmicScaleError != 0 && !fPoissonError){
86  const T scale = 1 + syst.GetShift(&kCosmicBkgScaleSyst) * fCosmicScaleError;
87  acosm *= scale;
88  }
89  if(fPoissonError){
90  T sigma = syst.GetShift(&kCosmicBkgScaleSyst);
91  for(int i = 1; i<= acosm.size()-2; i++){
92  if (sigma>0) acosm[i] *= (1+sigma*fBinErrors[i-1].first);
93  if (sigma<0) acosm[i] *= (1+sigma*fBinErrors[i-1].second);
94  }
95  }
96 
97  apred += acosm;
98  }
99 
100  return apred;
101  }
std::optional< Spectrum > fCosmic
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
Double_t scale
Definition: plot.C:25
if(dump)
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:49
std::vector< std::pair< double, double > > fBinErrors
std::vector< float > Spectrum
Definition: Constants.h:610
double sigma(TH1F *hist, double percentile)
std::vector< double > POT
double T
Definition: Xdiff_gwt.C:5
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
void ana::SingleSampleExperiment::SaveTo ( TDirectory *  dir,
const std::string name 
) const
overridevirtualinherited

Reimplemented from ana::IExperiment.

Definition at line 133 of file SingleSampleExperiment.cxx.

References dir, ana::SingleSampleExperiment::fCosmic, ana::SingleSampleExperiment::fData, ana::SingleSampleExperiment::fMC, ana::IPrediction::SaveTo(), ana::Spectrum::SaveTo(), and tmp.

134  {
135  TDirectory* tmp = dir;
136 
137  dir = dir->mkdir(name.c_str()); // switch to subdir
138  dir->cd();
139 
140  TObjString("SingleSampleExperiment").Write("type");
141 
142  fMC->SaveTo(dir, "mc");
143  fData.SaveTo(dir, "data");
144 
145  if(fCosmic.has_value()) fCosmic->SaveTo(dir, "cosmic");
146 
147  dir->Write();
148  delete dir;
149 
150  tmp->cd();
151  }
const XML_Char * name
Definition: expat.h:151
std::optional< Spectrum > fCosmic
Float_t tmp
Definition: plot.C:36
virtual void SaveTo(TDirectory *dir, const std::string &name) const
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
TDirectory * dir
Definition: macro.C:5

Member Data Documentation

std::vector< std::pair <double, double> > ana::SingleSampleExperiment::fBinErrors
protectedinherited
std::optional<Spectrum> ana::SingleSampleExperiment::fCosmic
protectedinherited
double ana::SingleSampleExperiment::fCosmicScaleError
protectedinherited
std::map<Component, Eigen::MatrixXd> ana::CovMxExperiment::fCovMxInv
protected

Definition at line 105 of file CovMxExperiment.h.

Referenced by ChiSq(), and MakeInvCovMxs().

Spectrum ana::SingleSampleExperiment::fData
protectedinherited
const IPrediction* ana::SingleSampleExperiment::fMC
protectedinherited
bool ana::SingleSampleExperiment::fPoissonError
protectedinherited

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