Functions
MakeCovMx.C File Reference
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Prediction/PredictionConcat.h"
#include "CAFAna/Prediction/CovarianceMatrix.h"
#include "NuXAna/Systs/Nus18Systs.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Calcs.h"
#include "OscLib/OscCalcSterile.h"
#include "NuXAna/macros/Nus18/Nus18Utilities.h"

Go to the source code of this file.

Functions

void MakeCovMx (bool rhc=false, std::string syst_type="all")
 

Function Documentation

void MakeCovMx ( bool  rhc = false,
std::string  syst_type = "all" 
)

Definition at line 35 of file MakeCovMx.C.

References ana::assert(), b, calc, om::cout, ana::DefaultSterileCalc(), ana::Progress::Done(), allTimeWatchdog::endl, MakeMiniprodValidationCuts::f, shutoffs::filename, GetComponents(), cet::getenv(), GetNus18Binning(), ana::GetNus18Systs(), GetPPFXWeightsDir(), hists, MECModelEnuComparisons::i, site_stats_from_log::job, ana::kAna2018FHCLivetime, ana::kAna2018FHCPOT, ana::kAna2018RHCLivetime, ana::kAna2018RHCPOT, ana::kAna2018SensitivityFHCNDPOT, ana::kAna2018SensitivityRHCNDPOT, LoadPPFX(), LoadPrediction(), file_size_ana::progress, SetAngles(), ana::Progress::SetProgress(), ana::ISyst::ShortName(), string, and systs.

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
52  double nd_pot = rhc? kAna2018SensitivityRHCNDPOT : kAna2018SensitivityFHCNDPOT;
53  double fd_pot = rhc? kAna2018RHCPOT : kAna2018FHCPOT;
54  double fd_livetime = rhc? kAna2018RHCLivetime : kAna2018FHCLivetime;
55 
56  // Define oscillation calculator
57  osc::OscCalcSterile* calc = DefaultSterileCalc(4);
58  SetAngles(calc);
59  osc::OscCalcSterileTrivial* calc_trivial
60  = new osc::OscCalcSterileTrivial(*DefaultSterileCalc(4));
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 }
systs
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
GetComponents
std::vector< std::tuple< ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t > > GetComponents()
Definition: Nus18Utilities.h:933
ana::ISyst::ShortName
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
ana::SystShifts
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
osc::OscCalcSterile
Adapt the PMNS_Sterile calculator to standard interface.
Definition: OscCalcSterile.h:23
file_size_ana.progress
bool progress
Insert implicit nodes #####.
Definition: file_size_ana.py:101
osc::OscCalcSterileTrivial
Version of OscCalcSterile that always returns probability of 1.
Definition: OscCalcSterile.h:76
shutoffs.filename
string filename
Definition: shutoffs.py:106
ana::kAna2018SensitivityFHCNDPOT
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
hists
TString hists[nhists]
Definition: bdt_com.C:3
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
SetAngles
void SetAngles(osc::OscCalcSterile *calc)
Definition: Nus18Utilities.h:165
ana::ISyst
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
ana::kAna2018RHCPOT
const double kAna2018RHCPOT
Definition: Exposures.h:208
ana::DefaultSterileCalc
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
cet::getenv
std::string getenv(std::string const &name)
LoadPPFX
IPrediction * LoadPPFX(std::string detector, int universe)
Function to load PPFX universe simulations.
Definition: Nus18Utilities.h:312
GetPPFXWeightsDir
std::string GetPPFXWeightsDir()
Function to get location of PPFX weights.
Definition: Nus18Utilities.h:492
om::cout
OStream cout
Definition: OStream.cxx:6
GetNus18Binning
std::vector< Binning > GetNus18Binning(bool rhc, std::string detector)
Definition: Nus18Utilities.h:133
LoadPrediction
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
Definition: Nus18Utilities.h:225
site_stats_from_log.job
job
Definition: site_stats_from_log.py:84
ana::kAna2018SensitivityRHCNDPOT
const double kAna2018SensitivityRHCNDPOT
Definition: Exposures.h:211
b
const hit & b
Definition: hits.cxx:21
ana::assert
assert(nhit_max >=nhit_nbins)
ana::Progress
A simple ascii-art progress bar.
Definition: Progress.h:9
ana::IPrediction
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
ana::kAna2018FHCPOT
const double kAna2018FHCPOT
Definition: Exposures.h:207
ana::DontAddDirectory
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
ana::GetNus18Systs
std::vector< const ISyst * > GetNus18Systs(bool rhc, std::string det_type, std::string syst_type)
Definition: Nus18Systs.cxx:439
ana::kAna2018FHCLivetime
const double kAna2018FHCLivetime
Definition: Exposures.h:213
ana::kAna2018RHCLivetime
const double kAna2018RHCLivetime
Definition: Exposures.h:214
string
enum BeamMode string
Generated on Thu Jan 21 2021 14:16:39 for NOvA by   doxygen 1.8.11