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"
6 #include "CAFAna/Core/HistCache.h"
10 #include "NuXAna/nus/Nus19FHC/Utilities.h"
11 
12 using namespace ana;
13 
14 void MakeCovMx() {
15 
16  // In a utility header, we define covmx::Sample objects for the near and far
17  // detector FHC NC selections. They are instantiated as kNusFHCNearDet and
18  // kNusFHCFarDet – you can see they are heavily used throughout this macro
19  // to keep track of the two samples.
20 
21  // Get systs and predictions
22  auto nd_systs = GetSystsFromFile(kNusFHCNearDet);
23  auto fd_systs = GetSystsFromFile(kNusFHCFarDet);
24  IPrediction* nd_pred_systs = LoadPrediction(kNusFHCNearDet, "systs");
25  IPrediction* fd_pred_systs = LoadPrediction(kNusFHCFarDet, "systs");
26  //IPrediction* nd_pred_genie = LoadPrediction(kNusFHCNearDet, "genie");
27  //IPrediction* fd_pred_genie = LoadPrediction(kNusFHCFarDet, "genie");
28 
29  // Create SystConcat objects to define how systematics are correlated
30  // between different samples.
31 
32  // Start with a GENIE systematic which should be correlated between the
33  // two detectors, for which we use a SingleSyst object. The same ISyst
34  // is used across all samples, so we instantiate a SingleSyst object to
35  // deal with it.
36 /* SingleSyst ma_cc_res;
37  ma_cc_res.syst = SystRegistry::ShortNameToSyst("MaCCRES");
38  ma_cc_res.onesided = false;
39  SystConcat sc_ma_cc_res({kNusFHCNearDet,kNusFHCFarDet}, {ma_cc_res}, {});
40 */
41  // Now we consider the Cherenkov systematic, which has a separate ISyst
42  // object for each sample. Here we use the MultiSyst class to correlate
43  // this systematic between the two samples. It's also a one-sided
44  // systematic, so we need to tag it as such to only throw positive shifts
45  // when producing the matrix.
46  MultiSyst cherenkov;
47  cherenkov.name = "Cherenkov";
48  cherenkov.systs = {SystRegistry::ShortNameToSyst(kNusFHCNearDet.GetTag()+"_Cherenkov"),
49  SystRegistry::ShortNameToSyst(kNusFHCFarDet.GetTag()+"_Cherenkov")};
50  cherenkov.onesided = true;
51  SystConcat sc_cherenkov({kNusFHCNearDet,kNusFHCFarDet}, {}, {cherenkov});
52 
53  // Now we do the AbsCalib systematic, which has a separate ISyst object
54  // for each sample. However, we want to treat calibration as uncorrelated
55  // between the two detectors, so we define a separate MultiSyst for each one
56  // so they'll be thrown independently.
57  MultiSyst abs_calib_nd;
58  abs_calib_nd.name = "AbsCalib";
59  abs_calib_nd.systs = {SystRegistry::ShortNameToSyst(kNusFHCNearDet.GetTag()+"_AbsCalib"),
60  nullptr};
61  abs_calib_nd.onesided = false;
62  MultiSyst abs_calib_fd;
63  abs_calib_fd.name = "AbsCalib";
64  abs_calib_fd.systs = { nullptr,
65  SystRegistry::ShortNameToSyst(kNusFHCFarDet.GetTag()+"_AbsCalib")};
66  abs_calib_fd.onesided = false;
67  SystConcat sc_abs_calib({ kNusFHCNearDet, kNusFHCFarDet } , {}, {abs_calib_nd,abs_calib_fd});
68 
69  // Get nominal prediction vector
70  TVectorD v_nd, v_fd, v_both;
71  osc::OscCalculatorSterileTrivial calc_trivial;
72  TH1D* h_nd = nd_pred_systs->Predict(&calc_trivial).ToTH1(12e20);
73  TH1D* h_fd = fd_pred_systs->Predict(&calc_trivial).ToTH1(12e20);
74  size_t nbins_nd = h_nd->GetNbinsX(), nbins_fd = h_fd->GetNbinsX();
75  v_nd.ResizeTo(nbins_nd);
76  v_fd.ResizeTo(nbins_fd);
77  v_both.ResizeTo(nbins_nd+nbins_fd);
78  for (size_t i = 0; i < nbins_nd; ++i) {
79  v_nd(i) = h_nd->GetBinContent(i+1);
80  v_both(i) = h_nd->GetBinContent(i+1);
81  }
82  for (size_t i = 0; i < nbins_fd; ++i) {
83  v_fd(i) = h_fd->GetBinContent(i+1);
84  v_both(nbins_nd+i) = h_fd->GetBinContent(i+1);
85  }
86  HistCache::Delete(h_nd);
87  HistCache::Delete(h_fd);
88 
89  auto calc = DefaultSterileCalc(4);
90 
91  // Output files
92  TFile* f_out = TFile::Open("covmx.root", "recreate");
93  TFile* f_plot = TFile::Open("covmx_plots.root", "recreate");
94 
95  // Produce GENIE covariance matrix
96 /* CovarianceMatrix mx_ma_cc_res({nd_pred_genie,fd_pred_genie}, {12e20,12e20}, sc_ma_cc_res);
97  mx_ma_cc_res.PredictCovMx({nd_pred_genie,fd_pred_genie}, {12e20,12e20}, calc);
98  mx_ma_cc_res.SaveTo(f_out->mkdir("ma_cc_res"));
99  TH2D* h_ma_cc_res = mx_ma_cc_res.GetCovMxRelativeTH2(v_both);
100  f_plot->WriteTObject(h_ma_cc_res, "ma_cc_res");
101  HistCache::Delete(h_ma_cc_res);
102 */
103  // Produce Cherenkov covariance matrix
104  CovarianceMatrix mx_cherenkov({nd_pred_systs,fd_pred_systs}, {12e20,12e20}, sc_cherenkov);
105  mx_cherenkov.PredictCovMx({nd_pred_systs,fd_pred_systs}, {12e20,12e20}, calc);
106  mx_cherenkov.SaveTo(f_out->mkdir("cherenkov"));
107  TH2D* h_cherenkov = mx_cherenkov.GetCovMxRelativeTH2(v_both);
108  f_plot->WriteTObject(h_cherenkov, "cherenkov");
109  HistCache::Delete(h_cherenkov);
110 
111  // Produce AbsCalib covariance matrix
112  CovarianceMatrix mx_abs_calib({nd_pred_systs,fd_pred_systs}, {12e20,12e20}, sc_abs_calib);
113  mx_abs_calib.PredictCovMx({nd_pred_systs,fd_pred_systs}, {12e20,12e20}, calc);
114  mx_abs_calib.SaveTo(f_out->mkdir("abs_calib"));
115  TH2D* h_abs_calib = mx_abs_calib.GetCovMxRelativeTH2(v_both);
116  f_plot->WriteTObject(h_abs_calib, "abs_calib");
117  HistCache::Delete(h_abs_calib);
118 
119  // Clean up
120  delete nd_pred_systs;
121  delete fd_pred_systs;
122 // delete nd_pred_genie;
123 // delete fd_pred_genie;
124  delete f_out;
125  delete f_plot;
126 
127 }
128 
129 
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void MakeCovMx()
Definition: MakeCovMx.C:14
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
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
std::string GetTag() const
Definition: Sample.cxx:94
osc::OscCalcDumb calc
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
const covmx::Sample kNusFHCFarDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kFarDet)
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
std::vector< const ISyst * > GetSystsFromFile(covmx::Sample sample)
Get systematics for a given sample.