MakeCovMx.C
Go to the documentation of this file.
1 /*
2  * MakeCovMx.C
3  *
4  * Created on: March 22, 2018
5  * Author: J. Hewes <jhewes15@fnal.gov>
6  *
7  * Generate covariance matrix including systematic uncertainties
8  */
9 
10 #ifdef __CINT__
11 
12 #include <iostream>
13 void MakeCovMx(bool rhc = false, std::string syst_type = "all")
14 {
15  std::cerr << "Sorry, you must run in compiled mode" << std::endl;
16 }
17 
18 #else
19 
21 #include "CAFAna/Core/SystShifts.h"
22 #include "CAFAna/Core/Progress.h"
28 #include "CAFAna/Analysis/Calcs.h"
29 #include "OscLib/OscCalcSterile.h"
30 
32 
33 using namespace ana;
34 
35 void MakeCovMx(bool rhc = false, std::string syst_type = "all")
36 {
37  DontAddDirectory guard;
38 
39  // Grid options
40  int job = -1;
41  char* job_env = std::getenv("PROCESS");
42  if (!job_env) std::cout << "Running all systematics locally." << std::endl;
43  else {
44  job = std::stoi(job_env);
45  std::cout << "Grid process detected, running job " << job << "..." << std::endl;
46  }
47 
48  // Polarity flag
49  std::string polarity = rhc? "RHC" : "FHC";
50 
51  // Set the total POT and livetime
53  double fd_pot = rhc? kAna2018RHCPOT : kAna2018FHCPOT;
54  double fd_livetime = rhc? kAna2018RHCLivetime : kAna2018FHCLivetime;
55 
56  // Define oscillation calculator
58  SetAngles(calc);
59  osc::OscCalcSterileTrivial* calc_trivial
61 
62  // PPFX matrix
63  if (syst_type == "ppfx") {
64 
65  // Load predictions
66  IPrediction* nd_pred = LoadPrediction("nd", rhc, "none");
67  IPrediction* fd_pred = LoadPrediction("fd", rhc, "none");
68 
69  // Near detector prediction
70  PredictionConcat* pc_nd = new PredictionConcat({nd_pred},
71  {nd_pot}, {0});
72  pc_nd->SetSysts(new SystConcat({}, {polarity+"_ND"}));
73  pc_nd->SetNewBinEdges(GetNus18Binning(rhc, "nd"));
74 
75  // Far detector prediction
76  PredictionConcat* pc_fd = new PredictionConcat({fd_pred},
77  {fd_pot}, {fd_livetime});
78  pc_fd->SetSysts(new SystConcat({}, {polarity+"_FD"}));
79  pc_fd->SetNewBinEdges(GetNus18Binning(rhc, "fd"));
80 
81  // Both detectors prediction
82  PredictionConcat* pc_both = new PredictionConcat({nd_pred, fd_pred},
83  {nd_pot, fd_pot}, {0, fd_livetime});
84  pc_both->SetSysts(new SystConcat({}, {polarity+"_ND", polarity+"_FD"}));
85  pc_both->SetNewBinEdges(GetNus18Binning(rhc, "both"));
86 
87  // Get PPFX universes
88  std::vector<std::vector<TH1D*>> ppfx_nd;
89  std::vector<std::vector<TH1D*>> ppfx_fd;
90  std::vector<std::vector<TH1D*>> ppfx_both;
91 
92  Progress progress("Loading PPFX universes");
93 
94  for (int i = 0; i < 100; ++i) {
95 
96  progress.SetProgress((double)i/100.);
97 
98  std::vector<TH1D*> hists;
99 
100  // Load universe from file
101  std::string ppfx_dir = GetPPFXWeightsDir();
102  IPrediction* nd_ppfx_pred = LoadPPFX("nd", i);
103  IPrediction* fd_ppfx_pred = LoadPPFX("fd", i);
104 
105  assert(nd_ppfx_pred && fd_ppfx_pred);
106 
107  // ND prediction
108  PredictionConcat* pc_ppfx_nd = new PredictionConcat(
109  {nd_ppfx_pred}, {nd_pot}, {0});
110  pc_ppfx_nd->SetNewBinEdges(GetNus18Binning(rhc, "nd"));
111  for (auto comp : GetComponents()) {
112  TH1D* nom = pc_nd->PredictComponent(calc,
113  std::get<0>(comp), std::get<1>(comp), std::get<2>(comp)).ToTH1(nd_pot);
114  TH1D* shift = pc_ppfx_nd->PredictComponent(calc,
115  std::get<0>(comp), std::get<1>(comp), std::get<2>(comp)).ToTH1(nd_pot);
116  for (int b = 1; b <= nom->GetNbinsX(); ++b) {
117  if (nom->GetBinContent(b) > 0)
118  shift->SetBinContent(b, shift->GetBinContent(b) / nom->GetBinContent(b));
119  else
120  shift->SetBinContent(b, 1);
121  }
122  hists.push_back(shift);
123  delete nom;
124  }
125  ppfx_nd.push_back(hists);
126  hists.clear();
127  delete pc_ppfx_nd;
128 
129  // FD prediction
130  PredictionConcat* pc_ppfx_fd = new PredictionConcat(
131  {fd_ppfx_pred}, {fd_pot}, {fd_livetime});
132  pc_ppfx_fd->SetNewBinEdges(GetNus18Binning(rhc, "fd"));
133  for (auto comp : GetComponents()) {
134  TH1D* nom = pc_fd->PredictComponent(calc,
135  std::get<0>(comp), std::get<1>(comp), std::get<2>(comp)).ToTH1(fd_pot);
136  TH1D* shift = pc_ppfx_fd->PredictComponent(calc,
137  std::get<0>(comp), std::get<1>(comp), std::get<2>(comp)).ToTH1(fd_pot);
138  for (int b = 1; b <= nom->GetNbinsX(); ++b) {
139  if (nom->GetBinContent(b) > 0)
140  shift->SetBinContent(b, shift->GetBinContent(b) / nom->GetBinContent(b));
141  else
142  shift->SetBinContent(b, 1);
143  }
144  hists.push_back(shift);
145  delete nom;
146  }
147  ppfx_fd.push_back(hists);
148  hists.clear();
149  delete pc_ppfx_fd;
150 
151  // Both detectors prediction
152  PredictionConcat* pc_ppfx_both = new PredictionConcat(
153  {nd_ppfx_pred, fd_ppfx_pred}, {nd_pot, fd_pot}, {0, fd_livetime});
154  pc_ppfx_both->SetNewBinEdges(GetNus18Binning(rhc, "both"));
155  for (auto comp : GetComponents()) {
156  TH1D* nom = pc_both->PredictComponent(calc,
157  std::get<0>(comp), std::get<1>(comp), std::get<2>(comp)).ToTH1(fd_pot);
158  TH1D* shift = pc_ppfx_both->PredictComponent(calc,
159  std::get<0>(comp), std::get<1>(comp), std::get<2>(comp)).ToTH1(fd_pot);
160  for (int b = 1; b <= nom->GetNbinsX(); ++b) {
161  if (nom->GetBinContent(b) > 0)
162  shift->SetBinContent(b, shift->GetBinContent(b) / nom->GetBinContent(b));
163  else
164  shift->SetBinContent(b, 1);
165  }
166  hists.push_back(shift);
167  delete nom;
168  }
169  ppfx_both.push_back(hists);
170  hists.clear();
171  delete pc_ppfx_both;
172 
173  // Clean up
174  delete nd_ppfx_pred;
175  delete fd_ppfx_pred;
176  }
177 
178  progress.Done();
179 
180  int n_universes = 1e5;
181 
182  // Output file
183  std::string filename = "covmx/covmx_" + polarity + "_ppfx.root";
184  TFile* f = new TFile(filename.c_str(), "recreate");
185 
186  // Near detector matrix
187  CovarianceMatrix gen_nd(pc_nd, calc_trivial, {}, nd_pot, n_universes, ppfx_nd);
188  gen_nd.PredictCovMx(pc_nd, calc);
189  gen_nd.SaveTo(f, "CovMx_ND_PPFX");
190 
191  // Far detector matrix
192  CovarianceMatrix gen_fd(pc_fd, calc_trivial, {}, fd_pot, n_universes, ppfx_fd);
193  gen_fd.PredictCovMx(pc_fd, calc);
194  gen_fd.SaveTo(f, "CovMx_FD_PPFX");
195 
196  // Both detectors matrix
197  CovarianceMatrix gen_both(pc_both, calc_trivial, {}, fd_pot, n_universes, ppfx_both);
198  gen_both.PredictCovMx(pc_both, calc);
199  gen_both.SaveTo(f, "CovMx_Both_PPFX");
200 
201  f->Close();
202 
203  } else {
204 
205  // Load predictions
206  IPrediction* nd_pred = LoadPrediction("nd", rhc, syst_type);
207  IPrediction* fd_pred = LoadPrediction("fd", rhc, syst_type);
208 
209  // Get systematics
210  std::vector<const ISyst*> systs = GetNus18Systs(rhc, "base", syst_type);
211 
212  // Near detector simulation
213  PredictionConcat* pc_nd = new PredictionConcat({nd_pred},
214  {nd_pot}, {0.});
215  pc_nd->SetSysts(new SystConcat(systs, {polarity + "_ND"}));
216  pc_nd->SetNewBinEdges(GetNus18Binning(rhc, "nd"));
217 
218  // Far detector simulation
219  PredictionConcat* pc_fd = new PredictionConcat({fd_pred},
220  {fd_pot}, {fd_livetime});
221  pc_fd->SetSysts(new SystConcat(systs, {polarity + "_FD"}));
222  pc_fd->SetNewBinEdges(GetNus18Binning(rhc, "fd"));
223 
224  // Two-detector simulation
225  PredictionConcat* pc_both = new PredictionConcat({nd_pred, fd_pred},
226  {nd_pot, fd_pot}, {0., fd_livetime});
227  pc_both->SetSysts(new SystConcat(systs,
228  {polarity + "_ND", polarity + "_FD"}));
229  pc_both->SetNewBinEdges(GetNus18Binning(rhc, "both"));
230 
231  if (job == -1) {
232 
233  // Open output file
234  std::string filename = "covmx/covmx_" + polarity + "_" + syst_type + ".root";
235  TFile* f = new TFile(filename.c_str(),"recreate");
236 
237  // Define systematics
238  SystShifts shifts;
239  int n_done = 0;
240 
241 
242  for (auto syst : systs) {
243 
244  // Near detector
245  std::cout << "Producing " << syst->ShortName()
246  << " matrix for near detector..." << std::endl;
247  CovarianceMatrix gen_nd(pc_nd, calc_trivial, {syst}, nd_pot, 1e5);
248  gen_nd.PredictCovMx(pc_nd, calc);
249  std::string matrix_name = "CovMx_ND_" + syst->ShortName();
250  gen_nd.SaveTo(f, matrix_name.c_str());
251 
252  // Far detector
253  std::cout << "Producing " << syst->ShortName()
254  << " matrix for far detector..." << std::endl;
255  CovarianceMatrix gen_fd(pc_fd, calc_trivial, {syst}, fd_pot, 1e5);
256  gen_fd.PredictCovMx(pc_fd, calc);
257  matrix_name = "CovMx_FD_" + syst->ShortName();
258  gen_fd.SaveTo(f, matrix_name.c_str());
259 
260  // Both detectors
261  std::cout << "Producing " << syst->ShortName()
262  << " matrix for both detectors..." << std::endl;
263  CovarianceMatrix gen_both(pc_both, calc_trivial, {syst}, fd_pot, 1e5);
264  gen_both.PredictCovMx(pc_both, calc);
265  matrix_name = "CovMx_Both_" + syst->ShortName();
266  gen_both.SaveTo(f, matrix_name.c_str());
267 
268  n_done++;
269  std::cout << std::endl << "Processed " << n_done << " of " << systs.size()
270  << " systs." << std::endl << std::endl;
271  }
272 
273  // Finish up
274  f->Close();
275 
276  }
277 
278  else {
279 
280  const ISyst* syst = systs[job];
281 
282  // Open output file
283  std::string filename = "covmx_" + polarity + "_" + syst->ShortName() + ".root";
284  TFile* f = new TFile(filename.c_str(),"recreate");
285 
286  // Define systematics
287  SystShifts shifts;
288 
289  // Near detector
290  std::cout << "Producing " << syst->ShortName()
291  << " matrix for near detector..." << std::endl;
292  CovarianceMatrix gen_nd(pc_nd, calc_trivial, {syst}, nd_pot, 1e5);
293  gen_nd.PredictCovMx(pc_nd, calc);
294  std::string matrix_name = "CovMx_ND_" + syst->ShortName();
295  gen_nd.SaveTo(f, matrix_name.c_str());
296 
297  // Far detector
298  std::cout << "Producing " << syst->ShortName()
299  << " matrix for far detector..." << std::endl;
300  CovarianceMatrix gen_fd(pc_fd, calc_trivial, {syst}, fd_pot, 1e5);
301  gen_fd.PredictCovMx(pc_fd, calc);
302  matrix_name = "CovMx_FD_" + syst->ShortName();
303  gen_fd.SaveTo(f, matrix_name.c_str());
304 
305  // Both detectors
306  std::cout << "Producing " << syst->ShortName()
307  << " matrix for both detectors..." << std::endl;
308  CovarianceMatrix gen_both(pc_both, calc_trivial, {syst}, fd_pot, 1e5);
309  gen_both.PredictCovMx(pc_both, calc);
310  matrix_name = "CovMx_Both_" + syst->ShortName();
311  gen_both.SaveTo(f, matrix_name.c_str());
312 
313  // Finish up
314  f->Close();
315 
316  }
317  }
318 }
319 
320 #endif
321 
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void MakeCovMx()
Definition: MakeCovMx.C:14
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::vector< std::tuple< ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t > > GetComponents()
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Adapt the PMNS_Sterile calculator to standard interface.
OStream cerr
Definition: OStream.cxx:7
bool progress
Insert implicit nodes #####.
Version of OscCalcSterile that always returns probability of 1.
string filename
Definition: shutoffs.py:106
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
TString hists[nhists]
Definition: bdt_com.C:3
osc::OscCalcDumb calc
void SetAngles(osc::OscCalcSterile *calc)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const double kAna2018RHCPOT
Definition: Exposures.h:208
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
std::string getenv(std::string const &name)
IPrediction * LoadPPFX(std::string detector, int universe)
Function to load PPFX universe simulations.
std::string GetPPFXWeightsDir()
Function to get location of PPFX weights.
OStream cout
Definition: OStream.cxx:6
std::vector< Binning > GetNus18Binning(bool rhc, std::string detector)
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
const double kAna2018SensitivityRHCNDPOT
Definition: Exposures.h:211
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
Definition: Progress.h:9
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
std::vector< const ISyst * > GetNus18Systs(bool rhc, std::string det_type, std::string syst_type)
Definition: Nus18Systs.cxx:439
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const double kAna2018RHCLivetime
Definition: Exposures.h:214
enum BeamMode string