MakeCovMx.C
Go to the documentation of this file.
1 /// MakeCovMx.C by J. Hewes (jhewes15@fnal.gov)
2 // Macro to make covariance matrices
3 
4 #include "CAFAna/Core/Sample.h"
5 #include "CAFAna/Core/ISyst.h"
10 
11 using namespace ana;
12 
13 void MakeCovMx() {
14 
15  // Get samples
20 
21  // Get systs and predictions
22  IPrediction* nd_pred = LoadPrediction(nd, "systs");
23  IPrediction* fd_pred = LoadPrediction(fd, "systs");
24  auto sc = GetSystConcats({nd,fd}, {nd_pred,fd_pred}, GetSystsFromFile(kNusFHCNearDet));
25 
26  // Get nominal prediction vector
27  TVectorD v_nd, v_fd, v_both;
29  TH1D* h_nd = nd_pred->Predict(&calc_trivial).ToTH1(12e20);
30  TH1D* h_fd = fd_pred->Predict(&calc_trivial).ToTH1(12e20);
31  size_t nbins_nd = h_nd->GetNbinsX(), nbins_fd = h_fd->GetNbinsX();
32  v_nd.ResizeTo(nbins_nd);
33  v_fd.ResizeTo(nbins_fd);
34  v_both.ResizeTo(nbins_nd+nbins_fd);
35  for (size_t i = 0; i < nbins_nd; ++i) {
36  v_nd(i) = h_nd->GetBinContent(i+1);
37  v_both(i) = h_nd->GetBinContent(i+1);
38  }
39  for (size_t i = 0; i < nbins_fd; ++i) {
40  v_fd(i) = h_fd->GetBinContent(i+1);
41  v_both(nbins_nd+i) = h_fd->GetBinContent(i+1);
42  }
43  HistCache::Delete(h_nd);
44  HistCache::Delete(h_fd);
45 
46  auto calc = DefaultSterileCalc(4);
47 
48  // Output files
49  TFile* f_out = TFile::Open("covmx.root", "recreate");
50  TFile* f_plot = TFile::Open("covmx_plots.root", "recreate");
51 
52  // One-sided systs
53  std::vector<std::string> onesided = { "AbsCalibShape", "Cherenkov", "XSecTune" };
54 
55  // ND matrices
56  TDirectory* dir = f_out->mkdir("nd");
57  TDirectory* plotdir = f_plot->mkdir("nd");
58  for (SystConcat it_sc : sc) {
59  CovarianceMatrix mx({nd_pred}, {nd}, {12e20}, it_sc.GetSysts(0), onesided);
60  mx.PredictCovMx({nd_pred}, {12e20}, calc);
61  auto names = it_sc.GetIndividualSystNames();
62  mx.SaveTo(dir->mkdir(names[0].c_str()));
63  TH2D* h = mx.GetCovMxRelativeTH2(v_nd);
64  plotdir->WriteTObject(h, names[0].c_str());
66  }
67 
68  // FD matrices
69  dir = f_out->mkdir("fd");
70  plotdir = f_plot->mkdir("fd");
71  for (SystConcat it_sc : sc) {
72  CovarianceMatrix mx({fd_pred}, {fd}, {12e20}, it_sc.GetSysts(1), onesided);
73  mx.PredictCovMx({fd_pred}, {12e20}, calc);
74  auto names = it_sc.GetIndividualSystNames();
75  mx.SaveTo(dir->mkdir(names[0].c_str()));
76  TH2D* h = mx.GetCovMxRelativeTH2(v_fd);
77  plotdir->WriteTObject(h, names[0].c_str());
79  }
80 
81  // Two-detector matrices
82  dir = f_out->mkdir("both");
83  plotdir = f_plot->mkdir("both");
84  for (SystConcat it_sc : sc) {
85  CovarianceMatrix mx({nd_pred,fd_pred}, {12e20,12e20}, it_sc);
86  mx.PredictCovMx({nd_pred,fd_pred}, {12e20,12e20}, calc);
87  auto names = it_sc.GetIndividualSystNames();
88  mx.SaveTo(dir->mkdir(names[0].c_str()));
89  TH2D* h = mx.GetCovMxRelativeTH2(v_both);
90  plotdir->WriteTObject(h, names[0].c_str());
92  }
93 
94  delete nd_pred;
95  delete fd_pred;
96 
97  // Now GENIE systematics
98  nd_pred = LoadPrediction(nd, "genie");
99  fd_pred = LoadPrediction(fd, "genie");
100 
101  dir = f_out->GetDirectory("nd");
102  plotdir = f_plot->GetDirectory("nd");
103  for (const ISyst* syst : getAllXsecSysts_2018()) {
104  CovarianceMatrix mx({nd_pred}, {nd}, {12e20}, {syst}, {});
105  mx.PredictCovMx({nd_pred}, {12e20}, calc);
106  mx.SaveTo(dir->mkdir(syst->ShortName().c_str()));
107  TH2D* h = mx.GetCovMxRelativeTH2(v_nd);
108  plotdir->WriteTObject(h, syst->ShortName().c_str());
110  }
111 
112  dir = f_out->GetDirectory("fd");
113  plotdir = f_plot->GetDirectory("fd");
114  for (const ISyst* syst : getAllXsecSysts_2018()) {
115  CovarianceMatrix mx({fd_pred}, {fd}, {12e20}, {syst}, {});
116  mx.PredictCovMx({fd_pred}, {12e20}, calc);
117  mx.SaveTo(dir->mkdir(syst->ShortName().c_str()));
118  TH2D* h = mx.GetCovMxRelativeTH2(v_fd);
119  plotdir->WriteTObject(h, syst->ShortName().c_str());
121  }
122 
123  dir = f_out->GetDirectory("both");
124  plotdir = f_plot->GetDirectory("both");
125  for (const ISyst* syst : getAllXsecSysts_2018()) {
126  CovarianceMatrix mx({nd_pred,fd_pred}, {nd,fd}, {12e20,12e20}, {syst}, {});
127  mx.PredictCovMx({nd_pred,fd_pred}, {12e20,12e20}, calc);
128  mx.SaveTo(dir->mkdir(syst->ShortName().c_str()));
129  TH2D* h = mx.GetCovMxRelativeTH2(v_both);
130  plotdir->WriteTObject(h, syst->ShortName().c_str());
132  }
133 
134  delete f_out;
135  delete f_plot;
136 
137 }
138 
139 
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
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:553
TMatrixD PredictCovMx(std::vector< IPrediction * > preds, std::vector< double > pot, osc::IOscCalculatorAdjustable *calc, SystConcat *sc=nullptr, const SystShifts &systs=SystShifts::Nominal())
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
static void Delete(TH1D *&h)
Definition: HistCache.cxx:215
Version of OscCalculatorSterile that always returns probability of 1.
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
void MakeCovMx(bool rhc=false, std::string syst_type="all")
Definition: MakeCovMx.C:35
TDirectory * dir
Definition: macro.C:5
string syst
Definition: plotSysts.py:176
std::vector< SystConcat > GetSystConcats(std::vector< covmx::Sample > samples, std::vector< IPrediction * > preds, std::vector< const ISyst * > systs)
Get correctly concatenated systs for ND and FD.
Definition: Utilities.h:294
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
osc::OscCalculatorSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
h
Definition: demo3.py:41
std::vector< const ISyst * > GetSystsFromFile(covmx::Sample sample)
Get systematics for a given sample.