Utilities.h
Go to the documentation of this file.
1 /// Utilities.h by J. Hewes (jhewes15@fnal.gov)
2 /// Utility functions for Nus19 covariance analysis
3 
4 #pragma once
5 
7 
9 #include "CAFAna/Core/Sample.h"
11 #include "CAFAna/Systs/SystConcat.h"
17 #include "NuXAna/Analysis/NusSystsMaker.h"
18 
19 #include "TFile.h"
20 #include "TKey.h"
21 
22 using namespace ana;
23 
24 std::string nus19_path = "/pnfs/nova/resilient/users/jhewes15/nus19";
25 
26 // Hard-coded parameters
27 const double kNus19Dm21 = 7.53e-5;
28 const double kNus19Dm32 = 2.44e-3;
29 const double kNus19Dm32Sigma = 0.08e-3;
30 
31 const double kNus19Th12 = 0.587;
32 const double kNus19Th13 = 0.145;
33 const double kNus19Th23 = 0.844;
34 const double kNus19Th23Sigma = 0.042;
35 
36 const double kNus19Delta13 = 1.21*M_PI;
37 const double kNus19Delta24 = M_PI/2;
38 
39 const double kNus19SensLivetime = 554.6519;
40 
45 
46 //===========================================================================================
47 bool CheckOption(TString opts, TString opt) {
48 
49  if (opts.Contains(opt, TString::ECaseCompare::kIgnoreCase))
50  return true;
51  else return false;
52 
53 }
54 
55 //===========================================================================================
56 std::vector<bool> GetOptionConfig(TString opt, std::vector<std::string> options) {
57 
58  size_t n_opts = options.size();
59 
60  std::vector<bool> vec(n_opts);
61  for (size_t i = 0; i < n_opts; ++i)
62  if (CheckOption(opt, options[i]))
63  vec[i] = true;
64 
65  if (std::none_of(vec.begin(), vec.end(), [](bool v) { return v; })) {
66  std::string err = "Option string must contain ";
67  for (size_t i = 0; i < n_opts-2; ++i)
68  err += "\"" + options[i] + "\", ";
69  err += options[n_opts-2] + " or " + options[n_opts-1] + ".";
70  throw std::runtime_error(err);
71  }
72 
73  else {
74  std::cout << "Running for";
75  for (size_t i = 0; i < n_opts; ++i)
76  if (vec[i])
77  std::cout << " " << options[i];
79  }
80 
81  return vec;
82 
83 } // function GetOptionConfig
84 
85 //===========================================================================================
86 /// Get prediction for a given sample
87 IPrediction* LoadPrediction(covmx::Sample sample, std::string systtype, bool cvmfs=true) {
88 
89  DontAddDirectory guard;
90 
91  // Only FHC NC for now
92  if (sample.analysis != covmx::Analysis::kNC
93  or sample.polarity != covmx::Polarity::kFHC)
94  assert(false and "Only NC and FHC samples supported for now.");
95 
96  std::string det = sample.detector == covmx::Detector::kNearDet ? "neardet" : "fardet";
97 
98  // Get path
100  if (cvmfs) {
101  // commenting out bc nusdata environment variables are temporarily hecked up
102  //const char* cpath = getenv("NUSDATA_NUS19FHC_PRED");
103  const char* cpath = getenv("NUSDATA_DIR");
104  assert (cpath and "nusdata directory not found! Aborting.");
105  path = std::string(cpath) + "/nus19/pred";
106  }
107  else {
108  path = nus19_path + "/pred";
109  }
110 
111  // Check syst type
112  std::string filename, predname;
113  if (systtype == "systs") {
114  filename = path+"/pred_nus_fhc_"+det+"_systs.root";
115  predname = "pred_interp_"+sample.GetTag()+"_systs";
116  }
117  else if (systtype == "genie") {
118  filename = path+"/pred_nus_fhc_"+det+"_genie.root";
119  predname = "pred_interp_"+sample.GetTag()+"_genie";
120  }
121  else if (systtype == "nosysts") {
122  filename = path+"/pred_nus_fhc_"+det+"_systs.root";
123  predname = "pred_"+sample.GetTag()+"_systs";
124  }
125  else assert(false and "Syst type must be \"systs\", \"genie\" or \"nosysts\".");
126 
127  // Load prediction from file
128  if (not cvmfs) filename = pnfs2xrootd(filename);
129  std::cout << "Loading " << predname << " from " << filename << std::endl;
130  TFile* f = TFile::Open(filename.c_str(), "read");
131  IPrediction* ret = nullptr;
132  if (systtype == "nosysts" and det == "neardet")
133  ret = NDPredictionSterile::LoadFrom(f, predname.c_str()).release();
134  else if (systtype == "nosysts" and det == "fardet")
135  ret = FDPredictionSterile::LoadFrom(f, predname.c_str()).release();
136  else
137  ret = PredictionInterp::LoadFrom(f, predname.c_str()).release();
138  delete f;
139  return ret;
140 
141 } // function LoadPrediction
142 
143 //===========================================================================================
144 /// Get cosmics for a given sample
145 std::unique_ptr<Spectrum> LoadCosmic(covmx::Sample sample, bool cvmfs=true) {
146 
147  DontAddDirectory guard;
148 
149  // Only FHC NC for now
150  if (sample.analysis != covmx::Analysis::kNC
151  or sample.polarity != covmx::Polarity::kFHC)
152  assert(false and "Only NC and FHC samples supported for now.");
153 
154  // No cosmics for ND
155  if (sample.detector == covmx::Detector::kNearDet)
156  return nullptr;
157 
158  // Get path
160  if (cvmfs) {
161  // commenting out bc nusdata environment variables are temporarily hecked up
162  //const char* cpath = getenv("NUSDATA_NUS19FHC_COSMIC");
163  const char* cpath = getenv("NUSDATA_DIR");
164  assert (cpath and "nusdata directory not found! Aborting.");
165  path = std::string(cpath) + "/nus19/cosmic";
166  }
167  else {
168  path = nus19_path + "/cosmic";
169  }
170 
171  std::string filename = path + "/nus19_fhc_cosmic.root";
172  if (not cvmfs) filename = pnfs2xrootd(filename);
173  TFile* f = TFile::Open(filename.c_str(), "read");
174  auto ret = LoadFromFile<Spectrum>(f, sample.GetTag());
175  f->Close();
176  return std::move(ret);
177 
178 } // function LoadCosmic
179 
180 //===========================================================================================
181 /// Get systematics for a given sample
182 std::vector<const ISyst*> GetSystsFromFile(covmx::Sample sample, bool cvmfs=true) {
183 
184  DontAddDirectory guard;
185 
186  // Only supported for NC right now
187  if (sample.analysis != covmx::Analysis::kNC or sample.polarity != covmx::Polarity::kFHC)
188  assert(false and "Only NC FHC is supported right now!");
189 
190  // Get path
192  if (cvmfs) {
193  // commenting out bc nusdata environment variables are temporarily hecked up
194  //const char* cpath = getenv("NUSDATA_NUS19FHC_ISYSTS");
195  const char* cpath = getenv("NUSDATA_DIR");
196  assert (cpath and "nusdata directory not found! Aborting.");
197  path = std::string(cpath) + "/nus19/isysts";
198  }
199  else {
200  path = nus19_path + "/isysts";
201  }
202 
203  std::string filepath = path+"/isysts_"+sample.GetDetector()+".root";
204  if (not cvmfs) filepath = pnfs2xrootd(filepath);
205 
206  // std::cout << "Loading systs for sample " << sample.GetTag() << ":" << std::endl;
207  TFile* f = TFile::Open(filepath.c_str(), "read");
208  TDirectory* dir = f->GetDirectory(sample.GetTag().c_str());
209  std::vector<const ISyst*> systs;
210  TIter next(dir->GetListOfKeys());
211  TKey* key;
212  while ((key = (TKey*)next()))
213  systs.push_back(NusISyst::LoadFrom((TDirectory*)key->ReadObj()).release());
214  for (size_t i = 0; i < systs.size(); ++i) {
215  if (systs[i]->ShortName().find("totErr") != std::string::npos) {
216  systs.erase(systs.begin()+i);
217  break;
218  }
219  }
220 
221  // for (auto syst : systs)
222  // std::cout << " " << syst->ShortName() << std::endl;
223  // std::cout << std::endl;
224 
225  f->Close();
226  return systs;
227 
228 } // function GetSystsFromFile
229 
230 //===========================================================================================
231 /// Get the last part of a systematic's ShortName
233 
234  size_t pos = name.find("_");
235  if (pos == std::string::npos) return name;
236 
237  while (pos != std::string::npos) {
238  size_t tmppos = name.find("_", pos+1);
239  if (tmppos == std::string::npos) break;
240  pos = tmppos;
241  }
242 
243  return name.substr(pos+1);
244 
245 } // function GetSystBaseName
246 
247 //===========================================================================================
248 /// Get correctly concatenated systs for ND and FD
249 SystConcat GetSystConcat(std::vector<covmx::Sample> samples,
250  std::vector<const ISyst*> systs, std::vector<IPrediction*> preds) {
251 
252  // Awkward, kinda hacky way to get the raw ISyst ShortNames. I don't like it but it works,
253  // for now at least
254  std::vector<std::string> shortnames;
255  for (auto syst : systs) {
256  shortnames.push_back(GetSystBaseName(syst->ShortName()));
257  std::cout << "Adding ShortName " << shortnames.back() << std::endl;
258  }
259 
260  std::set<std::string> uncorrelated = { "AbsCalib", "AbsCalibShape", "LightLevel" };
261  std::set<std::string> onesided = { "AbsCalibShape", "Cherenkov", "XSecTune" };
262 
263  std::vector<MultiSyst> ms;
264  for (std::string shortname : shortnames) {
265  std::vector<const ISyst*> these_systs;
266  for (auto s : samples)
267  these_systs.push_back(SystRegistry::ShortNameToSyst(s.GetTag()+"_"+shortname));
268  if (uncorrelated.count(shortname)) {
269  for (size_t i = 0; i < samples.size(); ++i) {
270  MultiSyst ms_tmp;
271  ms_tmp.name = shortname;
272  std::vector<const ISyst*> tmpsysts(preds.size(), nullptr);
273  tmpsysts[i] = these_systs[i];
274  ms_tmp.systs = tmpsysts;
275  ms_tmp.onesided = onesided.count(shortname);
276  ms.push_back(ms_tmp);
277  }
278  }
279  else {
280  MultiSyst tmp;
281  tmp.name = shortname;
282  tmp.systs = these_systs;
283  tmp.onesided = onesided.count(shortname);
284  ms.push_back(tmp);
285  }
286  }
287 
288  return SystConcat(samples, {}, ms);
289 
290 } // function GetSystConcats
291 
292 //===========================================================================================
293 /// Get correctly concatenated systs for ND and FD
294 std::vector<SystConcat> GetSystConcats(std::vector<covmx::Sample> samples,
295  std::vector<IPrediction*> preds, std::vector<const ISyst*> systs) {
296 
297  // Awkward, kinda hacky way to get the raw ISyst ShortNames. I don't like it but it works,
298  // for now at least
299  std::vector<std::string> shortnames;
300  for (auto syst : systs) {
301  size_t pos = syst->ShortName().find("_");
302  if (pos == std::string::npos) {
303  pos = 0;
304  break;
305  }
306  while (pos != std::string::npos) {
307  size_t tmppos = syst->ShortName().find("_", pos+1);
308  if (tmppos == std::string::npos) break;
309  pos = tmppos;
310  }
311  shortnames.push_back(syst->ShortName().substr(pos+1));
312  std::cout << "Adding ShortName " << shortnames.back() << std::endl;
313  }
314 
315  std::set<std::string> uncorrelated = { "AbsCalib", "AbsCalibShape", "LightLevel" };
316  std::set<std::string> onesided = { "AbsCalibShape", "Cherenkov", "XSecTune" };
317  std::vector<SystConcat> ret;
318  for (std::string shortname : shortnames) {
319  std::vector<const ISyst*> these_systs;
320  for (auto s : samples)
321  these_systs.push_back(SystRegistry::ShortNameToSyst(s.GetTag()+"_"+shortname));
322  std::vector<MultiSyst> ms;
323  if (uncorrelated.count(shortname)) {
324  for (size_t i = 0; i < samples.size(); ++i) {
325  MultiSyst ms_tmp;
326  ms_tmp.name = shortname;
327  std::vector<const ISyst*> tmpsysts(preds.size(), nullptr);
328  tmpsysts[i] = these_systs[i];
329  ms_tmp.systs = tmpsysts;
330  ms_tmp.onesided = onesided.count(shortname);
331  ms.push_back(ms_tmp);
332  }
333  }
334  else {
335  MultiSyst tmp;
336  tmp.name = shortname;
337  tmp.systs = these_systs;
338  tmp.onesided = onesided.count(shortname);
339  ms = { tmp };
340  }
341  SystConcat sc(samples, {}, ms);
342  ret.push_back(sc);
343  }
344 
345  return ret;
346 
347 } // function GetSystConcats
348 
349 //===========================================================================================
350 std::vector<CovarianceMatrix*> GetMatrices(std::vector<covmx::Sample> samples,
351  bool cvmfs=true) {
352 
353  std::string dirname;
354  if (samples.size() == 1 and samples[0] == kNusFHCNearDet) dirname = "nd";
355  else if (samples.size() == 1 and samples[0] == kNusFHCFarDet) dirname = "fd";
356  else if (samples.size() == 2 and samples[0] == kNusFHCNearDet
357  and samples[1] == kNusFHCFarDet) dirname = "both";
358  else assert (false and "Only NC FHC neardet & fardet supported for now.");
359 
361  if (cvmfs) {
362  // commenting out bc nusdata environment variables are temporarily hecked up
363  //const char* cpath = getenv("NUSDATA_NUS19FHC_COVMX");
364  const char* cpath = getenv("NUSDATA_DIR");
365  assert (cpath and "nusdata directory not found! Aborting.");
366  path = std::string(cpath) + "/nus19/covmx";
367  }
368  else {
369  path = nus19_path+"/covmx";
370  }
371  std::string filepath = path + "/covmx.root";
372  if (not cvmfs) filepath = pnfs2xrootd(filepath);
373  TFile* f = TFile::Open(filepath.c_str(), "read");
374 
375  TDirectory* dir = f->GetDirectory(dirname.c_str());
376  TIter next(dir->GetListOfKeys());
377  TKey* key;
378  std::vector<CovarianceMatrix*> ret;
379  while ((key = (TKey*)next())) {
380  std::string name(key->GetName());
381  if (name.find("totErr") != std::string::npos) continue;
382  CovarianceMatrix* mx = CovarianceMatrix::LoadFrom((TDirectory*)key->ReadObj()).release();
383  mx->SetName(name);
384  ret.push_back(mx);
385  }
386 
387  return ret;
388 
389 } // function GetMatrices
390 
391 //===========================================================================================
393 
394  // Default vals
396 
397  // Mass splittings
398  calc->SetDm(2, kNus19Dm21);
399  calc->SetDm(3, kNus19Dm21 + kNus19Dm32);
400  calc->SetDm(4, 0.5);
401 
402  // Angles
403  calc->SetAngle(1, 3, kNus19Th13);
404  calc->SetAngle(2, 3, kNus19Th23);
405  calc->SetAngle(3, 4, 0.);
406  calc->SetAngle(2, 4, 0.);
407  calc->SetAngle(1, 4, 0.);
408 
409  // Phases
410  calc->SetDelta(1, 3, kNus19Delta13);
411  calc->SetDelta(2, 4, kNus19Delta24);
412 }
413 
414 //===========================================================================================
415 /// Gaussian constrains on atmospheric parameters
416 std::vector<const IExperiment*> GetConstraints() {
417 
418  std::vector<const IExperiment*> ret;
420  ret.push_back(new const GaussianConstraint(&kFitDmSq32Sterile, kNus19Dm32, kNus19Dm32Sigma));
421  return ret;
422 
423 }
424 
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
const XML_Char * name
Definition: expat.h:151
const double kNus19Dm32
Definition: Utilities.h:28
static constexpr Double_t ms
Definition: Munits.h:192
const double kNus19Dm32Sigma
Definition: Utilities.h:29
static std::unique_ptr< FDPredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string GetSystBaseName(std::string name)
Get the last part of a systematic&#39;s ShortName.
Definition: Utilities.h:232
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
A simple Gaussian constraint on an arbitrary IFitVar.
const double kNus19Th12
Definition: Utilities.h:31
bool CheckOption(TString opts, TString opt)
Definition: Utilities.h:47
void SetDelta(int i, int j, double delta)
void SetAngles(osc::OscCalcSterile *calc)
Definition: Utilities.h:392
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
const double kNus19Th13
Definition: Utilities.h:32
std::string GetTag() const
Definition: Sample.cxx:94
string filename
Definition: shutoffs.py:106
Float_t tmp
Definition: plot.C:36
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
osc::OscCalcDumb calc
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
#define M_PI
Definition: SbMath.h:34
const XML_Char * s
Definition: expat.h:262
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
Definition: Utilities.h:416
static std::unique_ptr< NDPredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
const double kNus19Th23Sigma
Definition: Utilities.h:34
const double kNus19Delta24
Definition: Utilities.h:37
std::vector< bool > GetOptionConfig(TString opt, std::vector< std::string > options)
Definition: Utilities.h:56
std::string getenv(std::string const &name)
Detector detector
Definition: Sample.h:100
IPrediction * LoadPrediction(covmx::Sample sample, std::string systtype, bool cvmfs=true)
Get prediction for a given sample.
Definition: Utilities.h:87
Eigen::VectorXd vec
const covmx::Sample kNusFHCFarDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kFarDet)
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
Definition: Utilities.h:350
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:76
SystConcat GetSystConcat(std::vector< covmx::Sample > samples, std::vector< const ISyst * > systs, std::vector< IPrediction * > preds)
Get correctly concatenated systs for ND and FD.
Definition: Utilities.h:249
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
const double kNus19Th23
Definition: Utilities.h:33
void SetAngle(int i, int j, double th)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string GetDetector() const
Definition: Sample.cxx:58
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
void SetDm(int i, double dm)
TDirectory * dir
Definition: macro.C:5
std::string nus19_path
Definition: Utilities.h:24
std::unique_ptr< Spectrum > LoadCosmic(covmx::Sample sample, bool cvmfs=true)
Get cosmics for a given sample.
Definition: Utilities.h:145
assert(nhit_max >=nhit_nbins)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::vector< const ISyst * > GetSystsFromFile(covmx::Sample sample, bool cvmfs=true)
Get systematics for a given sample.
Definition: Utilities.h:182
const double kNus19Dm21
Definition: Utilities.h:27
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Polarity polarity
Definition: Sample.h:99
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
void next()
Definition: show_event.C:84
const double kNus19SensLivetime
Definition: Utilities.h:39
const double kNus19Delta13
Definition: Utilities.h:36
const FitTheta23Sterile kFitTheta23Sterile