MakeGENIEMatrix.C
Go to the documentation of this file.
1 
2 #include "CAFAna/Core/Sample.h"
9 
10 using namespace ana;
11 
13 
14  // Run for NC FHC, ND & FD
15 
16  std::vector<covmx::Sample> samples =
23  };
24 
25  // Get nominal concatenated prediction
26 
27  PredictionConcat* pred_nom = CovMxManager::GetPrediction(samples);
28  pred_nom->SetBinning({ FHCNDBins(), FHCFDBins() }); // new binning schemes defined in Nus18Utilities.h
29  // reeeeeally need to finalise binning schemes for each detector
30  // so we can stop rebinning downstream which is super annoying
31 
32  // Get the nominal prediction as a vector (so we can draw the matrix)
33 
34  auto calc = DefaultSterileCalc(4);
35  SetAngles(calc);
36  TVectorD v_pred(pred_nom->NBins());
37  TH1D* h_pred = pred_nom->Predict(calc).ToTH1(pred_nom->GetPOTs().back());
38  for (size_t i = 0; i < pred_nom->NBins(); ++i)
39  v_pred(i) = h_pred->GetBinContent(i+1);
40  delete h_pred;
41 
42  // Get all the GENIE systematics (PredictionInterps for these are already made!)
43 
44  std::vector<const ISyst*> genie_systs = getAllXsecSysts_2018();
45 
46  // Instantiate an empty matrix
47 
48  CovarianceMatrix* mx = new CovarianceMatrix(pred_nom->GetPOTs().back(),
49  pred_nom->GetBinning());
50 
51  // Output file!
52 
53  TFile* fout = TFile::Open("covmx.root", "recreate");
54  TDirectory* hist_dir = fout->mkdir("hists");
55 
56  // Now we loop over each systematic, make a matrix and add it to the
57  // originally empty matrix we defined earlier
58 
59  // Actually there are a ton of these so let's just loop over the first
60  // ten to save some time
61 
62  for (size_t i = 0; i < 10; ++i) {
63 
64  const ISyst* syst = genie_systs[i];
65 
66  // Get the prediction for this systematic
67 
68  PredictionConcat* pred_syst = CovMxManager::GetPrediction(samples, syst->ShortName());
69  pred_syst->SetBinning({ FHCNDBins(), FHCFDBins() });
70  std::vector<std::vector<const ISyst*>> corrsysts(2, std::vector<const ISyst*>());
71  pred_syst->SetSysts(new SystConcat({"nc_fhc_nd", "nc_fhc_fd"}, {syst}, corrsysts));
72 
73  // Actually make the (not empty) matrix for this syst
74 
75  CovarianceMatrix* mx_syst = new CovarianceMatrix(pred_syst, { syst }, pred_syst->GetPOTs().back());
76  mx_syst->PredictCovMx(pred_nom, calc);
77  TH2* h_mx = mx_syst->GetCovMxRelativeTH2(v_pred);
78  std::string histname = "mxhist_" + syst->ShortName();
79  hist_dir->WriteTObject(h_mx, histname.c_str());
80 
81  // Add it to our total matrix, save it, and then throw it out
82 
83  mx->AddMatrix(mx_syst);
84  std::string mxname = "covmx_" + syst->ShortName();
85  mx_syst->SaveTo(fout, mxname.c_str());
86  delete mx_syst;
87 
88  } // for syst
89 
90  // And now save our combined matrix too!
91 
92  mx->PredictCovMx(pred_nom, calc);
93  mx->SaveTo(fout, "covmx_all");
94 
95  TH2* h_mx = mx->GetCovMxRelativeWithStatsTH2(v_pred);
96  hist_dir->WriteTObject(h_mx, "mxhist_all");
97 
98  delete mx;
99  delete fout;
100 
101  // bye
102 
103 }
104 
105 
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
IPrediction * GetPrediction(SystPredType predType, const HistAxis *axis, const HistAxis *numuAxis, const Cut *FDCut, const Cut *NDCut, const Cut *numuCut, const SystShifts *shiftData, const SystShifts *shiftMC, const Var *weight, Loaders *loaders, Loaders *loadersND)
Definition: SystMaker.cxx:44
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
osc::OscCalcDumb calc
void SetAngles(osc::OscCalcSterile *calc)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
Binning FHCNDBins()
Binning FHCFDBins()
void MakeGENIEMatrix()
enum BeamMode string