11 #include "TDirectory.h" 21 using std::back_inserter;
28 using std::make_tuple;
29 using std::make_unique;
30 using std::ostringstream;
32 using std::runtime_error;
34 using std::unique_ptr;
37 using Eigen::MatrixXd;
38 using Eigen::VectorXd;
48 vector<Component>
ret;
56 flavours.insert(
end(flavours), begin(fdFlavours),
end(fdFlavours));
61 for (
auto & flavour: flavours) {
62 for (
auto &
sign: signs) {
75 vector<const ISyst*>
systs,
int nUniverses)
76 : fNUniverses(nUniverses)
84 MatrixXd covMxAbsolute, vector<pair<Binning, unsigned int>>
bins)
91 vector<Sample> samples)
119 vector<const ISyst*>
systs)
124 if (!
fBins.empty()) {
125 throw runtime_error(
"Cannot reset matrix binning!");
128 Binning bins =
s.GetPrediction()->Predict(calc).GetBinnings()[0];
142 vector<double> nominalPrediction(nBinsFull);
143 vector<double> fractionalShift(nBinsFull);
147 for (
size_t iPred = 0; iPred < samples.size(); ++iPred) {
151 Spectrum spec = samples[iPred].GetPrediction()->PredictComponent(calc,
154 for (
size_t iBin = 1; iBin <= (size_t)bins[iPred].NBins(); ++iBin) {
155 nominalPrediction[
i] = hist->GetBinContent(iBin);
168 for (Long64_t iThrow = 0; iThrow <
fNUniverses; ++iThrow) {
172 for (
const ISyst* syst : systs) {
173 double shift = rnd.Gaus();
181 for (
size_t iPred = 0; iPred < samples.size(); ++iPred) {
182 SystShifts sampleShifts = samples[iPred].GetSystShifts(shifts);
186 Spectrum shiftedSpec = samples[iPred].GetPrediction()->PredictComponentSyst(calc,
188 auto shiftedHist = shiftedSpec.
ToTH1(samples[iPred].
GetPOT());
191 for (
size_t iBin = 1; iBin <= (size_t)bins[iPred].NBins(); ++iBin) {
194 double shiftedPrediction = shiftedHist->GetBinContent(iBin);
196 throw runtime_error(
"NaN value found in shifted prediction.");
199 if (nominalPrediction[i] > 0) fractionalShift[
i] = ((shiftedPrediction+0.01)/(nominalPrediction[
i]+0.01)) - 1;
200 else fractionalShift[
i] = 0;
208 for (
size_t iBin = 0; iBin < nBinsFull; ++iBin) {
209 for (
size_t jBin = 0; jBin < nBinsFull; ++jBin) {
232 vector<double> oscPrediction(nBinsFull);
235 for (
size_t iPred = 0; iPred < samples.size(); ++iPred) {
238 SystShifts sampleShifts = samples[iPred].GetSystShifts(shifts);
239 Spectrum spec = samples[iPred].GetPrediction()->PredictComponentSyst(calc,
242 for (
int iBin = 1; iBin <= spec.
GetBinnings()[0].NBins(); ++iBin) {
243 oscPrediction[
i] = hist->GetBinContent(iBin);
251 size_t iOffset(0), iFullOffset(0);
253 for (
size_t iPred = 0; iPred < samples.size(); ++iPred) {
254 size_t iBins =
fBins[iPred].first.NBins();
255 for (
size_t iComp = 0; iComp < nComps[iPred]; ++iComp) {
257 for (
size_t iBin = 0; iBin < iBins; ++iBin) {
258 size_t jOffset(0), jFullOffset(0);
259 for (
size_t jPred = 0; jPred < samples.size(); ++jPred) {
260 size_t jBins =
fBins[jPred].first.NBins();
261 for (
size_t jComp = 0; jComp < nComps[jPred]; ++jComp) {
262 for (
size_t jBin = 0; jBin < jBins; ++jBin) {
263 size_t i(iOffset+iBin), iFull(iFullOffset+iBin),
264 j(jOffset+jBin), jFull(jFullOffset+jBin);
267 jFullOffset += jBins;
272 iFullOffset += iBins;
284 MatrixXd
ret(vBetaMu.size(), vBetaMu.size());
286 size_t iOffset(0), iFullOffset(0);
288 for (
size_t iPred = 0; iPred <
fBins.size(); ++iPred) {
289 size_t iBins =
fBins[iPred].first.NBins();
290 for (
size_t iComp = 0; iComp < nComps[iPred]; ++iComp) {
292 for (
size_t iBin = 0; iBin < iBins; ++iBin) {
293 size_t jOffset(0), jFullOffset(0);
294 for (
size_t jPred = 0; jPred <
fBins.size(); ++jPred) {
295 size_t jBins =
fBins[jPred].first.NBins();
296 for (
size_t jComp = 0; jComp < nComps[jPred]; ++jComp) {
297 for (
size_t jBin = 0; jBin < jBins; ++jBin) {
298 size_t i(iOffset+iBin), iFull(iFullOffset+iBin),
299 j(jOffset+jBin), jFull(jFullOffset+jBin);
300 ret(
i,
j) +=
fFullCovMx(iFull, jFull) * vBetaMu[iFull] * vBetaMu[jFull];
302 jFullOffset += jBins;
307 iFullOffset += iBins;
322 if (nBins != nBinsComp || nBinsFull != nBinsFullComp) {
323 throw runtime_error(
"Cannot add matrices; bins do not match.");
336 for (
size_t i = 0;
i < (size_t)ret.rows(); ++
i) {
348 [](
auto const& pair){
return pair.first; });
357 vector<double> mergedBins = { 0. };
360 for (
size_t i = 1;
i <
b.Edges().size(); ++
i) {
361 double width =
b.Edges()[
i] -
b.Edges()[
i-1];
362 mergedBins.push_back(mergedBins.back() + width);
374 vector<double> mergedBins = { 0. };
376 for (
auto const&
b :
bins) {
377 vector<double> edges =
b.first.Edges();
378 for (
size_t iComp = 0; iComp <
b.second; ++iComp) {
379 for (
size_t iBin = 1; iBin < edges.size(); ++iBin) {
380 double width = edges[iBin] - edges[iBin-1];
381 mergedBins.push_back(mergedBins.back() + width);
392 vector<unsigned int>
ret;
394 [](
auto const& pair){
return pair.first.NBins(); });
402 vector<unsigned int>
ret;
404 [](
auto const& pair){
return pair.second; });
416 for (
Sample& s1 : samples) {
417 TH1D*
h1 = s1.GetPrediction()->Predict(calc).ToTH1(s1.GetPOT());
418 for (
int b1 = 1; b1 <= h1->GetNbinsX(); ++b1) {
419 double iVal = h1->GetBinContent(b1);
421 for (
Sample& s2 : samples) {
422 TH1D*
h2 = s2.GetPrediction()->Predict(calc).ToTH1(s2.GetPOT());
423 for (
int b2 = 1; b2 <= h2->GetNbinsX(); ++b2) {
424 double jVal = h2->GetBinContent(b2);
425 ret(c1, c2) /= iVal * jVal;
446 for (
int i = 0;
i < corrMx.rows(); ++
i) {
447 for (
int j = 0;
j < corrMx.cols(); ++
j) {
451 if (sigma_1 > 0 && sigma_2 > 0) {
452 corrMx(
i,
j) = cov/(
sqrt(sigma_1) *
sqrt(sigma_2));
470 vector<double> edges = bins.
Edges();
472 ";Reco E bin (components);Reco E bin (components)",
473 nBins, &edges[0], nBins, &edges[0]);
475 for (
size_t i = 0;
i <
nBins; ++
i) {
476 for (
size_t j = 0;
j <
nBins; ++
j) {
491 vector<double> edges = bins.
Edges();
493 ";Reco E bin;Reco E bin",
494 nBins, &edges[0], nBins, &edges[0]);
496 for (
size_t i = 0;
i <
nBins; ++
i) {
497 for (
size_t j = 0;
j <
nBins; ++
j) {
513 vector<double> edges = bins.
Edges();
515 ";Reco E bin;Reco E bin",
516 nBins, &edges[0], nBins, &edges[0]);
519 for (
size_t i = 0;
i <
nBins; ++
i) {
520 for (
size_t j = 0;
j <
nBins; ++
j) {
521 hist->SetBinContent(
i+1,
j+1, mx(
i,
j));
535 ";Reco E bin;Reco E bin",
539 for (
int i = 0;
i < mxBins.
NBins(); ++
i) {
540 for (
int j = 0;
j < mxBins.
NBins(); ++
j) {
544 if(sigma_1 > 0 && sigma_2 > 0)
545 hist->SetBinContent(
i+1,
j+1, cov/(
sqrt(sigma_1)*
sqrt(sigma_2)));
546 else hist->SetBinContent(
i+1,
j+1, 0);
562 ";Reco E bin (components);Reco E bin (components)",
569 for (
int i = 0;
i < mxBins.
NBins(); ++
i) {
570 for (
int j = 0;
j < mxBins.
NBins(); ++
j) {
571 hist->SetBinContent(
i+1,
j+1, corrMx(
i,
j));
583 if (!
fBins.empty()) {
584 throw runtime_error(
"Cannot reset matrix binning!");
600 TDirectory*
tmp = gDirectory;
602 dir = dir->mkdir(name.c_str());
605 TObjString(
"CovarianceMatrix").Write(
"type");
609 fullCovMx.Write(
"fullcovmx");
616 covMxAbsolute.Write(
"covmxabsolute");
619 vector<vector<double>> edges(
fBins.size());
620 vector<unsigned int> nComps(
fBins.size(), 0);
625 dir->WriteObject(&edges,
"binedges");
626 dir->WriteObject(&nComps,
"ncomps");
638 dir = dir->GetDirectory(name.c_str());
641 TObjString*
tag = (TObjString*)dir->Get(
"type");
643 assert(tag->GetString() ==
"CovarianceMatrix");
648 MatrixXd fullCovMx = Map<MatrixXd>(tmp->GetMatrixArray(),
649 tmp->GetNrows(), tmp->GetNcols());
651 tmp = (
TMatrixD*)dir->Get(
"covmxabsolute");
653 MatrixXd covMxAbsolute = Map<MatrixXd>(tmp->GetMatrixArray(),
654 tmp->GetNrows(), tmp->GetNcols());
658 vector<vector<double>>* edges;
659 vector<unsigned int>* nComps;
660 vector<pair<Binning, unsigned int>>
bins;
661 dir->GetObject(
"binedges", edges);
662 dir->GetObject(
"ncomps", nComps);
663 for (
size_t it = 0;
it < edges->size(); ++
it) {
668 unique_ptr<CovarianceMatrix>
ret(make_unique<CovarianceMatrix>
669 (fullCovMx, covMxAbsolute, bins));
TH2D * GetCovMxRelativeTH2(std::vector< Sample > &samples, osc::IOscCalcAdjustable *calc)
const std::vector< Binning > & GetBinnings() const
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
vector< Component > GetComponents(Sample sample)
fvar< T > fabs(const fvar< T > &x)
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.
Eigen::MatrixXd GetFullCovMx()
Simple record of shifts applied to systematic parameters.
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
Eigen::MatrixXd GetCovMxAbsolute()
void AddMatrix(CovarianceMatrix *gen)
bool progress
Insert implicit nodes #####.
Version of OscCalcSterile that always returns probability of 1.
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Encapsulate code to systematically shift a caf::SRProxy.
TH2D * GetCovMxAbsoluteTH2()
Representation of a spectrum in any variable, with associated POT.
int isnan(const stan::math::var &a)
Charged-current interactions.
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
std::vector< Binning > GetBinnings()
TH2D * GetFullCorrMxTH2()
Eigen::MatrixXd GetCovMxRelative(std::vector< Sample > &samples, osc::IOscCalcAdjustable *calc)
static Eigen::MatrixXd ForcePosDef(Eigen::MatrixXd mx, double epsilon=0.1)
std::vector< std::string > comps
static Binning Custom(const std::vector< double > &edges)
const std::vector< double > & Edges() const
std::vector< unsigned int > GetNComponents()
std::vector< std::pair< Binning, unsigned int > > fBins
void SetProgress(double frac)
Update the progress fraction between zero and one.
std::vector< unsigned int > GetNBins()
Neutral-current interactions.
void SetBinning(std::vector< Sample > samples)
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
Both neutrinos and antineutrinos.
double GetPOT(bool isFHC=true)
Eigen::MatrixXd GetFullCorrMx()
static std::unique_ptr< CovarianceMatrix > LoadFrom(TDirectory *dir, const std::string &name)
virtual void SaveTo(TDirectory *dir, const std::string &name) const
void BuildFullCovMx(std::vector< Sample > samples, std::vector< const ISyst * > systs)
Eigen::MatrixXd Predict(std::vector< Sample > samples, osc::IOscCalcAdjustable *calc, const SystShifts &systs=SystShifts::Nominal())
All neutrinos, any flavor.
Prevent histograms being added to the current directory.
void Done()
Call this when action is completed.
unsigned int fNUniverses
Number of universes to simulate.
CovarianceMatrix()=delete
void SetShift(const ISyst *syst, double shift, bool force=false)
std::string UniqueName()
Return a different string each time, for creating histograms.
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...