Utilities.h
Go to the documentation of this file.
1 /// Utilities.h by J. Hewes (jhewes15@fnal.gov)
2 /// Utility functions for Nus20 covariance analysis
3 
4 #pragma once
5 
7 
9 #include "CAFAna/Core/Sample.h"
10 #include "CAFAna/Core/Registry.h"
17 
22 
23 #include "NuXAna/Vars/HistAxes.h"
24 #include "NuXAna/Cuts/NusCuts20.h"
29 
30 #include "TFile.h"
31 #include "TKey.h"
32 
33 using namespace ana;
34 
35 // Hard-coded parameters
36 const double kNux20Dm21 = 7.53e-5;
37 const double kNux20Dm32 = 2.44e-3;
38 const double kNux20Dm32Sigma = 0.08e-3;
39 
40 const double kNux20Th12 = 0.587;
41 const double kNux20Th13 = 0.145;
42 const double kNux20Th23 = 0.844;
43 const double kNux20Th23Sigma = 0.042;
44 
45 const double kNux20Delta13 = 1.21*M_PI;
46 const double kNux20Delta24 = M_PI/2;
47 
48 const double kNux20FHCNDPOT = 11e20;
50 const double kNux20RHCNDPOT = 11.8e20;
54 
65 
70 
71 //===========================================================================================
72 bool CheckOption(TString opts, TString opt) {
73 
74  if (opts.Contains(opt, TString::ECaseCompare::kIgnoreCase))
75  return true;
76  else return false;
77 
78 }
79 
80 //===========================================================================================
81 template <class T>
82 T ParseOption(std::vector<std::pair<std::string, T>> opts, TString optString, T def) {
83 
84  T ret = def;
85  size_t nOpts = 0;
86  for (auto opt : opts) {
87  if (optString.Contains(opt.first.c_str(), TString::kIgnoreCase)) {
88  ++nOpts;
89  ret = opt.second;
90  }
91  }
92  if (nOpts > 1) {
93  std::ostringstream err;
94  err << "Option string must include one of: ";
95  for (auto opt : opts) err << "\"" << opt.first << "\", ";
96  std::cerr << err.str().substr(0, err.str().size()-2) << std::endl;
97  assert(false);
98  }
99  return ret;
100 
101 } // template function ParseOption
102 
103 //===========================================================================================
104 template <class T>
105 std::vector<T> ParseOptions(std::vector<std::pair<std::string, T>> opts, TString optString) {
106 
107  std::vector<T> ret;
108  for (auto opt : opts) {
109  if (optString.Contains(opt.first.c_str(), TString::kIgnoreCase)) {
110  ret.push_back(opt.second);
111  }
112  }
113  return ret;
114 
115 } // template function ParseOption
116 
117 //===========================================================================================
118 std::vector<covmx::Sample> GetNumuNDSamples() {
119  std::vector<covmx::Sample> ret;
122  ret.push_back(covmx::Sample(covmx::kCCNumu, p, covmx::kNearDet, q));
123  }
124  }
125  return ret;
126 }
127 
128 //===========================================================================================
129 std::vector<covmx::Sample> GetNumuFDSamples() {
130  std::vector<covmx::Sample> ret;
133  ret.push_back(covmx::Sample(covmx::kCCNumu, p, covmx::kFarDet, q));
134  }
135  }
136  return ret;
137 }
138 
139 //===========================================================================================
140 std::vector<covmx::Sample> GetNumuSamples() {
141  std::vector<covmx::Sample> ret;
142  for (auto tmp : { GetNumuNDSamples(), GetNumuFDSamples() }) {
143  ret.insert(ret.end(), tmp.begin(), tmp.end());
144  }
145  return ret;
146 }
147 
148 //===========================================================================================
149 std::vector<covmx::Sample> GetNCNDSamples() {
150  std::vector<covmx::Sample> ret;
152  ret.push_back(covmx::Sample(covmx::kNC, p, covmx::kNearDet));
153  }
154  return ret;
155 }
156 
157 //===========================================================================================
158 std::vector<covmx::Sample> GetNCFDSamples() {
159  std::vector<covmx::Sample> ret;
161  ret.push_back(covmx::Sample(covmx::kNC, p, covmx::kFarDet));
162  }
163  return ret;
164 }
165 
166 //===========================================================================================
167 std::vector<covmx::Sample> GetNCSamples() {
168  std::vector<covmx::Sample> ret;
169  for (auto tmp : { GetNCNDSamples(), GetNCFDSamples() }) {
170  ret.insert(ret.end(), tmp.begin(), tmp.end());
171  }
172  return ret;
173 }
174 
175 //===========================================================================================
176 std::vector<covmx::Sample> GetAllSamples() {
177  std::vector<covmx::Sample> ret;
178  for (auto tmp : { GetNumuSamples(), GetNCSamples() }) {
179  ret.insert(ret.end(), tmp.begin(), tmp.end());
180  }
181  return ret;
182 }
183 
184 //===========================================================================================
185 double GetDefaultPOT(covmx::Sample sample) {
186  if (sample.detector == covmx::kNearDet) {
187  if (sample.polarity == covmx::kFHC) return kNux20FHCNDPOT;
188  else return kNux20RHCNDPOT;
189  } else {
190  if (sample.polarity == covmx::kFHC) return kNux20FHCFDPOT;
191  else return kNux20RHCFDPOT;
192  }
193 }
194 
195 //===========================================================================================
196 void SetPOT(covmx::Sample& sample) {
197  sample.SetPOT(GetDefaultPOT(sample));
198 }
199 
200 //===========================================================================================
202  if (sample.detector == covmx::kNearDet) return -1; // no livetime at ND
203  else {
204  if (sample.polarity == covmx::kFHC) return kNux20FHCLivetime;
205  else return kNux20RHCLivetime;
206  }
207 }
208 
209 //===========================================================================================
210 void SetLivetime(covmx::Sample& sample) {
211  sample.SetLivetime(GetDefaultLivetime(sample));
212 }
213 
214 //===========================================================================================
215 const HistAxis* GetDefaultAxis(covmx::Sample sample, double res=0) {
216 
217  if (sample.selection == covmx::kNC) {
218  if (res != 0) { // Resolution-based binning
219  std::vector<double> bins = { 0.5 };
220  while (true) {
221  double next = (1 + (res/100)) * bins.back();
222  double nextnext = (1 + (res/100)) * next;
223  if (nextnext > 20.) {
224  bins.push_back(20);
225  break;
226  } else {
227  bins.push_back(next);
228  }
229  } // while true
230  std::cout << "Created dynamic binning scheme with resolution " << res << ", and edges:";
231  for (double bin : bins) std::cout << " " << bin;
232  std::cout << endl;
233  return new HistAxis("Energy Deposited in Scintillator (GeV)",
235  } else {
236  if ( sample.detector == covmx::kNearDet ) return &kNC20NDAxisE;
237  else if ( sample.detector == covmx::kFarDet ) return &kNC20FDAxisE;
238  }
239  }
240  else if (sample.selection == covmx::kCCNumu){
241  return &kNumuCCOptimisedAxis2020;
242  }
243  else if (sample.selection == covmx::kCCNue){
244  return &kNue2020Axis;
245  }
246  else{
247  std::cout << "No default binning for " << sample.GetName() << std::endl;
248  assert(false);
249  }
250  return &kNC20NDAxisE; // moot return to avoid -Wreturn-type non-void warning
251 
252 }
253 
254 //===========================================================================================
255 void SetAxis(covmx::Sample& sample, double res=0) {
256  sample.SetAxis(GetDefaultAxis(sample, res));
257 }
258 
259 //===========================================================================================
260 
261 void SetCut(covmx::Sample& sample) {
262 
263  // Neutral current cuts
264  if (sample.selection == covmx::kNC) {
265  if (sample.detector == covmx::kNearDet) sample.SetCut(&kNus20NDCuts);
266  else sample.SetCut(&kNus20FDCuts_ML);
267  }
268 
269  else if (sample.selection == covmx::kCCNue) {
270  if (sample.detector == covmx::kNearDet) sample.SetCut(&kNue2020ND);
271  else sample.SetCut(&kNue2020FDAllSamples_ML);
272  }
273 
274  else if (sample.selection == covmx::kCCNumu) {
275  if (sample.quantile == covmx::kNoQ) {
276  if (sample.detector == covmx::kNearDet) sample.SetCut(&kNumu2020ND);
277  else sample.SetCut(&kNumu2020FD_ML);
278  } else {
279  std::string quantDir = std::getenv("NUMUDATA_DIR");
280  std::string infileQuant = quantDir + "/lib/ana2020/Quantiles/quantiles_" + sample.GetPolarity() + "_full_numu2020.root";
281  TFile* infile = TFile::Open( pnfs2xrootd(infileQuant).c_str() );
282  TH2* spec2DFD = (TH2*)infile->FindObjectAny("FDSpec2D");
283  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(spec2DFD, kNumuCCOptimisedAxis2020, kHadEFracAxis, 4);
284  infile->Close();
285  if (sample.detector == covmx::kNearDet) {
286  sample.SetCut(new const Cut(kNumu2020ND && HadEFracQuantCuts[sample.quantile]));
287  } else {
288  sample.SetCut(new const Cut(kNumu2020FD_ML && HadEFracQuantCuts[sample.quantile]));
289  }
290  } // if using quantiles
291  } // numu selection
292 
293 } // function SetCut
294 
295 //===========================================================================================
296 void SetCuts(std::vector<covmx::Sample>& samples) {
297 
298  // Instantiate this once so we can reuse it if necessary without reloading
299  std::vector<Cut> HadEFracQuantCuts;
300 
301  for (covmx::Sample& sample : samples) {
302 
303  // Neutral current cuts
304 
305  if (sample.selection == covmx::kNC) {
306  if (sample.detector == covmx::kNearDet) sample.SetCut(&kNus20NDCuts);
307  else sample.SetCut(&kNus20FDCuts_ML);
308  }
309 
310  else if (sample.selection == covmx::kCCNue) {
311  if (sample.detector == covmx::kNearDet) sample.SetCut(&kNue2020ND);
312  else sample.SetCut(&kNue2020FDAllSamples_ML);
313  }
314 
315  else if (sample.selection == covmx::kCCNumu) {
316  if (sample.quantile == covmx::kNoQ) {
317  if (sample.detector == covmx::kNearDet) sample.SetCut(&kNumu2020ND);
318  else sample.SetCut(&kNumu2020FD_ML);
319  } else {
320  if (HadEFracQuantCuts.size() == 0) { // Load quantiles
321  std::string quantDir = std::getenv("NUMUDATA_DIR");
322  std::string infileQuant = quantDir + "/lib/ana2020/Quantiles/quantiles_" + sample.GetPolarity() + "_full_numu2020.root";
323  TFile* infile = TFile::Open( pnfs2xrootd(infileQuant).c_str() );
324  TH2* spec2DFD = (TH2*)infile->FindObjectAny("FDSpec2D");
325  HadEFracQuantCuts = QuantileCutsFromTH2(spec2DFD, kNumuCCOptimisedAxis2020, kHadEFracAxis, 4);
326  infile->Close();
327  }
328 
329  if (sample.detector == covmx::kNearDet) {
330  sample.SetCut(new const Cut(kNumu2020ND && HadEFracQuantCuts[sample.quantile]));
331  } else {
332  sample.SetCut(new const Cut(kNumu2020FD_ML && HadEFracQuantCuts[sample.quantile]));
333  }
334  } // if using quantiles
335  } // numu selection
336  } // for sample
337 
338 } // function SetCuts
339 
340 //===========================================================================================
341 /// Function to take an option TString and return a single associated covmx::Sample
343 
344  // Get selection
345  std::vector<std::pair<std::string, covmx::Selection>> selOpts
346  = { {"nuesel", covmx::kCCNue}, {"numusel", covmx::kCCNumu}, {"ncsel", covmx::kNC} };
347  covmx::Selection sel = ParseOption<covmx::Selection>(selOpts, optString, covmx::kNoSel);
348  // Get polarity
349  std::vector<std::pair<std::string, covmx::Polarity>> polOpts
350  = { {"fhc", covmx::kFHC}, {"rhc", covmx::kRHC } };
351  covmx::Polarity pol = ParseOption<covmx::Polarity>(polOpts, optString, covmx::kNoPol);
352  // Get detector
353  std::vector<std::pair<std::string, covmx::Detector>> detOpts
354  = { {"neardet", covmx::kNearDet}, {"fardet", covmx::kFarDet} };
355  covmx::Detector det = ParseOption<covmx::Detector>(detOpts, optString, covmx::kNoDet);
356  // Define sample
357  covmx::Sample ret(sel, pol, det);
358  // Get quantile if applicable
359  if (sel == covmx::kCCNumu) {
360  std::vector<std::pair<std::string, covmx::Quantile>> qOpts
361  = { {"q0", covmx::kQ1}, {"q1", covmx::kQ2}, {"q2", covmx::kQ3}, {"q3", covmx::kQ4} };
362  covmx::Quantile q = ParseOption<covmx::Quantile>(qOpts, optString, covmx::kNoQ);
363  ret.quantile = q;
364  }
365  return ret;
366 
367 } // function GetSampleFromOptString
368 
369 //===========================================================================================
370 /// Function to take an option TString and return a vector of associated covmx::Samples
371 std::vector<covmx::Sample> GetSamplesFromOptString(TString optString) {
372 
373  // Get selections
374  std::vector<std::pair<std::string, covmx::Selection>> selOpts
375  = { {"nuesel", covmx::kCCNue}, {"numusel", covmx::kCCNumu}, {"ncsel", covmx::kNC} };
376  std::vector<covmx::Selection> sels
377  = ParseOptions<covmx::Selection>(selOpts, optString);
378  // Get polarities
379  std::vector<std::pair<std::string, covmx::Polarity>> polOpts
380  = { {"fhc", covmx::kFHC}, {"rhc", covmx::kRHC } };
381  std::vector<covmx::Polarity> pols
382  = ParseOptions<covmx::Polarity>(polOpts, optString);
383  // Get detectors
384  std::vector<std::pair<std::string, covmx::Detector>> detOpts
385  = { {"neardet", covmx::kNearDet}, {"fardet", covmx::kFarDet} };
386  std::vector<covmx::Detector> dets
387  = ParseOptions<covmx::Detector>(detOpts, optString);
388  bool quantiles = optString.Contains("quantiles");
389  // Check for optional no-quantile config
390  // std::vector<std::pair<std::string, covmx::Quantile>> qOpts
391  // = { {"noquantiles", covmx::kNoQ}, {"q0", covmx::kQ1}, {"q1", covmx::kQ2},
392  // {"q2", covmx::kQ3}, {"q3", covmx::kQ4} };
393  // std::vector<covmx::Quantile> quantiles
394  // = ParseOptions<covmx::Quantile>(qOpts, optString);
395 
396  // Loop over samples and add to vector
397  std::vector<covmx::Sample> ret;
398  for (covmx::Selection sel : sels) {
399  for (covmx::Polarity pol : pols) {
400  for (covmx::Detector det : dets) {
401  // If it's numu, add all the quantiles unless configured for no quantiles
402  if (sel == covmx::kCCNumu && quantiles) {
403  // if (quantiles.empty()) {
404  // ret.push_back(covmx::Sample(sel, pol, det, covmx::kNoQ));
405  // } else {
407  ret.push_back(covmx::Sample(sel, pol, det, q));
408  }
409  // }
410  }
411  else {
412  ret.push_back(covmx::Sample(sel, pol, det));
413  }
414  } // for detector
415  } // for polarity
416  } // for selection
417 
418  std::cout << "Samples: " << std::endl;
419  for(auto& s : ret){
420  SetAxis(s);
421  SetPOT(s);
422  SetLivetime(s);
423  std::cout << s.GetTag() << " (" << s.GetBinning().NBins() << " bins)" << std::endl;
424  }
425  return ret;
426 
427 } // function GetSampleFromOptString
428 
429 //===========================================================================================
430 std::vector<bool> GetOptionConfig(TString opt, std::vector<std::string> options) {
431 
432  size_t n_opts = options.size();
433 
434  std::vector<bool> vec(n_opts);
435  for (size_t i = 0; i < n_opts; ++i)
436  if (CheckOption(opt, options[i]))
437  vec[i] = true;
438 
439  if (std::none_of(vec.begin(), vec.end(), [](bool v) { return v; })) {
440  std::string err = "Option string must contain ";
441  for (size_t i = 0; i < n_opts-2; ++i)
442  err += "\"" + options[i] + "\", ";
443  err += options[n_opts-2] + " or " + options[n_opts-1] + ".";
444  throw std::runtime_error(err);
445  }
446 
447  else {
448  std::cout << "Running for";
449  for (size_t i = 0; i < n_opts; ++i)
450  if (vec[i])
451  std::cout << " " << options[i];
452  std::cout << std::endl;
453  }
454 
455  return vec;
456 
457 } // function GetOptionConfig
458 
459 
460 //===========================================================================================
463  or sample.polarity != covmx::Polarity::kFHC) {
464  assert(false and "Nus20 currently only supports up to FHC NC + Numu.");
465  }
466 }
467 
468 //===========================================================================================
469 /// Get cosmics for a given sample
470 std::unique_ptr<Spectrum> LoadCosmic(covmx::Sample sample) {
471 
472  DontAddDirectory guard;
473 
474  // No cosmics for ND
475  if (sample.detector == covmx::kNearDet)
476  return nullptr;
477 
478  std::string filename = InputPath() + "/cosmics/cosmics.root";
479  TFile* f = TFile::Open(filename.c_str(), "read");
480  auto ret = Spectrum::LoadFrom(f, sample.GetTag());
481  delete f;
482  return ret;
483 
484 } // function LoadCosmic
485 
486 //===========================================================================================
487 /// Set cosmics in a given sample
488 void SetCosmic(covmx::Sample& sample) {
489 
490  sample.SetCosmic(LoadCosmic(sample).release());
491 
492 } // function LoadCosmic
493 
494 //===========================================================================================
495 void SetInfo(covmx::Sample& sample, double res=0) {
496  SetAxis(sample, res);
497  SetPOT(sample);
498  SetCut(sample);
499  SetLivetime(sample);
500 }
501 
502 //===========================================================================================
503 void SetPrediction(covmx::Sample& sample, bool systs=true) {
504 
505  string filename = InputPath() + "/pred/pred_" + sample.GetTag() + ".root";
506 
507  // Load prediction from file
508  TFile* f = TFile::Open(filename.c_str(), "read");
509  IPrediction* p = nullptr;
510  if (systs) {
511  getNusAna2020AllSysts(sample);
512  p = PredictionInterp::LoadFrom(f, (sample.GetTag()+"/pred_interp").c_str()).release();
513  }
514  else {
515  p = (sample.detector == covmx::kNearDet) ?
516  (IPrediction*)NDPredictionSterile::LoadFrom(f, (sample.GetTag()+"/pred").c_str()).release() :
517  (IPrediction*)FDPredictionSterile::LoadFrom(f, (sample.GetTag()+"/pred").c_str()).release();
518  }
519  sample.SetPrediction(p);
520  delete f;
521 
522 }
523 
524 //===========================================================================================
525 void SetEverything(covmx::Sample& sample, bool systs=false, double res=0) {
526  SetInfo(sample, res);
527  getNusAna2020AllSysts(sample);
528  SetPrediction(sample, systs);
529 }
530 
531 //===========================================================================================
532 /// Get the last part of a systematic's ShortName
534 
535  size_t pos = name.find("_");
536  if (pos == std::string::npos) return name;
537 
538  while (pos != std::string::npos) {
539  size_t tmppos = name.find("_", pos+1);
540  if (tmppos == std::string::npos) break;
541  pos = tmppos;
542  }
543 
544  return name.substr(pos+1);
545 
546 } // function GetSystBaseName
547 
548 //===========================================================================================
549 std::vector<covmx::CovarianceMatrix*> GetMatrices(std::vector<covmx::Sample> samples,
550  bool cvmfs=true) {
551  for(auto sample : samples) AssertSupported(sample);
552 
553  DontAddDirectory guard;
554 
555  std::string dirname = "both";
556  if (samples.size() == 1 and samples[0] == kNusFHCNearDet) dirname = "nd";
557  else if (samples.size() == 1 and samples[0] == kNusFHCFarDet) dirname = "fd";
558  //else if (samples.size() == 2 and samples[0] == kNusFHCNearDet
559  // and samples[1] == kNusFHCFarDet) dirname = "both";
560  //else if (samples.size() == 10) dirname = "both";
561  //else assert (false and "Only NC FHC neardet & fardet supported for now.");
562 
564  if (cvmfs) {
565  // commenting out bc nusdata environment variables are temporarily hecked up
566  //const char* cpath = getenv("NUSDATA_NUS19FHC_COVMX");
567  const char* cpath = getenv("NUSDATA_DIR");
568  assert (cpath and "nusdata directory not found! Aborting.");
569  path = std::string(cpath) + "/nus19/covmx";
570  }
571  else {
572  path = kNus20Path+"/covmx";
573  }
574  std::string filepath = path + "/covmx.root";
575 
576  // Temporary hack so I can cautiously preserve the existing covmx file
577  for(auto sample : samples){
578  if(sample.selection == covmx::kCCNumu){
579  filepath = path + "/covmx_nus_numu_fhc.root";
580  break;
581  }
582  }
583 
584  //if (not cvmfs) filepath = pnfs2xrootd(filepath); will this work without xrootd?
585  TFile* f = TFile::Open(filepath.c_str(), "read");
586 
587  TDirectory* dir = f->GetDirectory(dirname.c_str());
588  TIter next(dir->GetListOfKeys());
589  TKey* key;
590  std::vector<covmx::CovarianceMatrix*> ret;
591  while ((key = (TKey*)next())) {
592  std::string name(key->GetName());
593  if (name.find("totErr") != std::string::npos) continue;
594  covmx::CovarianceMatrix* mx = covmx::CovarianceMatrix::LoadFrom((TDirectory*)key->ReadObj(), "").release();
595  mx->SetName(name);
596  ret.push_back(mx);
597  }
598 
599  delete f;
600  return ret;
601 
602 } // function GetMatrices
603 
604 //===========================================================================================
606 
607  covmx::Sample ndSample(kNusFHCNearDet);
608  SetAxis(ndSample);
609  covmx::Sample fdSample(kNusFHCFarDet);
610  SetAxis(fdSample);
611  covmx::CovarianceMatrix* mx = new covmx::CovarianceMatrix({ndSample,fdSample});
612  for (auto syst : GetMatrices({ndSample,fdSample}, cvmfs))
613  mx->AddMatrix(syst);
614  return mx;
615 
616 } // function GetFullMatrix
617 
618 //===========================================================================================
620 
621  if(type=="3flav"){
622  // Mass splittings
623  std::cout<<" 3-Flavor signal"<<std::endl;
624  calc->SetDm(2, kNux20Dm21);
625  calc->SetDm(3, kNux20Dm21 + kNux20Dm32);
626  calc->SetDm(4, 0.);
627 
628  // Angles
629  calc->SetAngle(1, 3, kNux20Th13);
630  calc->SetAngle(2, 3, kNux20Th23);
631  calc->SetAngle(3, 4, 0.);
632  calc->SetAngle(2, 4, 0.);
633  calc->SetAngle(1, 4, 0.);
634 
635  // Phases
636  calc->SetDelta(1, 3, kNux20Delta13);
637  calc->SetDelta(2, 4, kNux20Delta24);
638 
639  }
640  else if (type == "lsnd") { // fake signal 1: LSND
641  // Mass splittings
642  std::cout<<" LSND Like signal"<<std::endl;
643  calc->SetDm(2, kNux20Dm21);
644  calc->SetDm(3, kNux20Dm21 + kNux20Dm32);
645  calc->SetDm(4, 1.2);
646  // Angles
647  Double_t theta_24 = TMath::ASin(TMath::Sqrt(0.0166));
648  Double_t theta_14 = TMath::ASin(TMath::Sqrt(0.18))/2;
649  calc->SetAngle(1, 3, kNux20Th13);
650  calc->SetAngle(2, 3, kNux20Th23);
651  calc->SetAngle(3, 4, 0.);
652  calc->SetAngle(2, 4, theta_24);
653  calc->SetAngle(1, 4, theta_14);
654 
655  // Phases
656  calc->SetDelta(1, 3, kNux20Delta13);
657  calc->SetDelta(2, 4, kNux20Delta24);
658 
659  }
660  else if (type=="miniboone"){
661  std::cout<<" MiniBooNE Like signal"<<std::endl;
662  // Mass splittings
663  calc->SetDm(2, kNux20Dm21);
664  calc->SetDm(3, kNux20Dm21 + kNux20Dm32);
665  calc->SetDm(4, 0.041);
666 
667  // Angles
668  Double_t theta_24 = TMath::ASin(TMath::Sqrt(0.95));
669  calc->SetAngle(1, 3, kNux20Th13);
670  calc->SetAngle(2, 3, kNux20Th23);
671  calc->SetAngle(3, 4, 0.);
672  calc->SetAngle(2, 4, theta_24);
673  calc->SetAngle(1, 4, 0.785398);
674 
675  // Phases
676  calc->SetDelta(1, 3, kNux20Delta13);
677  calc->SetDelta(2, 4, kNux20Delta24);
678  }
679  else assert(false && ("Oscillation parameter set \"" + type + "\" not recognised.").c_str());
680 
681 }
682 
683 //===========================================================================================
685 
686  SetNus20Params(calc); // Start with nominal parameters always
687 
688  if (type == 1 || type == 2) { // fake signal 1 & 2: MiniBooNE sterile parameters
689  Double_t dm41 = 0.041;
690  Double_t theta_24 = TMath::ASin(TMath::Sqrt(0.9580));
691  Double_t theta_34 = TMath::ASin(TMath::Sqrt(0.9580));
692  calc->SetDm(4, dm41);
693  calc->SetAngle(2, 4, theta_24);
694  calc->SetAngle(3, 4, theta_34);
695  } else if (type == 3) {
696  return; // Use nominal parameters for this set!
697  } else { // Throw an error if fake signal type not recognised.
698  std::ostringstream err;
699  err << "Error in SetFakeSignalParams! Invalid fake signal set "
700  << type << " was requested!";
701  throw std::runtime_error(err.str());
702  }
703 
704 } // function SetFakeSignalParams
705 
706 //===========================================================================================
708 
709  std::cout << "---------- OSCILLATION PARAMETERS ----------" << std::endl
710  << "Dm21: " << calc->GetDm(2) << std::endl
711  << "Dm31: " << calc->GetDm(3) << std::endl
712  << "Dm41: " << calc->GetDm(4) << std::endl
713  << "Theta 13: " << calc->GetAngle(1, 3) << std::endl
714  << "Theta 23: " << calc->GetAngle(2, 3) << std::endl
715  << "Theta 14: " << calc->GetAngle(1, 4) << std::endl
716  << "Theta 24: " << calc->GetAngle(2, 4) << std::endl
717  << "Theta 34: " << calc->GetAngle(3, 4) << std::endl
718  << "Delta 13: " << calc->GetDelta(1, 3) << std::endl
719  << "Delta 24: " << calc->GetDelta(2, 4) << std::endl
720  << "--------------------------------------------" << std::endl;
721 }
722 
723 //===========================================================================================
724 /// Gaussian constrains on atmospheric parameters
725 std::vector<const IExperiment*> GetConstraints() {
726 
727  std::vector<const IExperiment*> ret;
729  ret.push_back(new const GaussianConstraint(&kFitDmSq32Sterile, kNux20Dm32, kNux20Dm32Sigma));
730  return ret;
731 
732 }
733 
734 //===========================================================================================
735 /// Load fake data from file
736 std::vector<Spectrum*> LoadFakeData(std::vector<covmx::Sample> samples, size_t universe,
737  bool systs)
738 {
739  DontAddDirectory guard;
740  std::vector<Spectrum*> ret(samples.size());
741 
742  std::string filepath = InputPath() + "/fakedata/FakeData" + (systs ? "Systs" : "Stats") + ".root";
743  //if (not cvmfs) filepath = pnfs2xrootd(filepath);
744  TFile* f = TFile::Open(filepath.c_str(), "read");
745 
746  // Get universe and spectra
747  TDirectory* dir = f->GetDirectory(Form("universe_%zu", universe));
748  for (size_t i = 0; i < ret.size(); ++i) {
749  ret[i] = Spectrum::LoadFrom(dir, samples[i].GetTag()).release();
750  }
751  return ret;
752 
753 } // function LoadFakeData
754 
755 // //===========================================================================================
756 // /// Get decorrelated near detector errors
757 // std::vector<double> GetDecorrelatedNDErrors() {
758 
759 // // Get covariance matrix and predictions
760 // covmx::CovarianceMatrix* mx = GetFullMatrix(kUseCVMFS);
761 // auto predND = LoadPrediction(kNusFHCNearDet, "nosysts", kUseCVMFS);
762 // auto predFD = LoadPrediction(kNusFHCFarDet, "nosysts", kUseCVMFS);
763 // covmx::Sample ndSample(kNusFHCNearDet);
764 // ndSample.SetPrediction(predND);
765 // covmx::Sample fdSample(kNusFHCFarDet);
766 // fdSample.SetPrediction(predFD);
767 // auto calc = DefaultSterileCalc(4);
768 // SetNus20Params(calc);
769 // TMatrixD matrix = mx->Predict({ndSample,fdSample}, calc);
770 
771 // size_t ndbins = kNCNDBinning.NBins();
772 // size_t fdbins = kNCFDBinning.NBins();
773 
774 // // Create and populate submatrices
775 // TMatrixD s11(ndbins, ndbins);
776 // TMatrixD s12(fdbins, ndbins);
777 // TMatrixD s21(ndbins, fdbins);
778 // TMatrixD s22(fdbins, fdbins);
779 
780 // for (size_t i = 0; i < ndbins; ++i) {
781 // for (size_t j = 0; j < ndbins; ++j) {
782 // s11(i,j) = matrix(i,j);
783 // }
784 // }
785 
786 // for (size_t i = 0; i < ndbins; ++i) {
787 // for (size_t j = 0; j < fdbins; ++j) {
788 // s12(j,i) = matrix(i,ndbins+j);
789 // s21(i,j) = matrix(i,ndbins+j);
790 // }
791 // }
792 
793 // for (size_t i = 0; i < fdbins; ++i) {
794 // for (size_t j = 0; j < fdbins; ++j) {
795 // s22(i,j) = matrix(ndbins+i, ndbins+j);
796 // }
797 // }
798 
799 // // Decorrelate
800 // TMatrixD tmp(ndbins, ndbins);
801 // s22.Invert();
802 // s21 *= s22;
803 // tmp.Mult(s21, s12);
804 // s11 -= tmp;
805 
806 // std::vector<double> ret(ndbins);
807 // for (size_t i = 0; i < ndbins; ++i) {
808 // ret[i] = sqrt(s11(i, i));
809 // }
810 
811 // // Clean up
812 // delete mx;
813 // delete predND;
814 // delete predFD;
815 // delete calc;
816 
817 // return ret;
818 
819 // } // function GetDecorrelatedNDErrors
820 
821 // //===========================================================================================
822 // /// Get decorrelated far detector errors
823 // std::vector<double> GetDecorrelatedFDErrors() {
824 
825 // // Get covariance matrix and predictions
826 // covmx::CovarianceMatrix* mx = GetFullMatrix(kUseCVMFS);
827 // auto predND = LoadPrediction(kNusFHCNearDet, "nosysts", kUseCVMFS);
828 // auto predFD = LoadPrediction(kNusFHCFarDet, "nosysts", kUseCVMFS);
829 // covmx::Sample ndSample(kNusFHCNearDet);
830 // ndSample.SetPrediction(predND);
831 // covmx::Sample fdSample(kNusFHCFarDet);
832 // fdSample.SetPrediction(predFD);
833 // auto calc = DefaultSterileCalc(4);
834 // SetNus20Params(calc);
835 // TMatrixD matrix = mx->Predict({ndSample,fdSample}, calc);
836 
837 // size_t ndbins = kNCNDBinning.NBins();
838 // size_t fdbins = kNCFDBinning.NBins();
839 
840 // TMatrixD s11(fdbins, fdbins);
841 // TMatrixD s12(ndbins, fdbins);
842 // TMatrixD s21(fdbins, ndbins);
843 // TMatrixD s22(ndbins, ndbins);
844 
845 // for (size_t i = 0; i < ndbins; ++i) {
846 // for (size_t j = 0; j < ndbins; ++j) {
847 // s22(i,j) = matrix(i,j);
848 // }
849 // }
850 
851 // for (size_t i = 0; i < ndbins; ++i) {
852 // for (size_t j = 0; j < fdbins; ++j) {
853 // s12(i,j) = matrix(i,ndbins+j);
854 // s21(j,i) = matrix(i,ndbins+j);
855 // }
856 // }
857 
858 // for (size_t i = 0; i < fdbins; ++i) {
859 // for (size_t j = 0; j < fdbins; ++j) {
860 // s11(i,j) = matrix(ndbins+i, ndbins+j);
861 // }
862 // }
863 
864 // TMatrixD tmp(fdbins, fdbins);
865 // s22.Invert();
866 // s21 *= s22;
867 // tmp.Mult(s21, s12);
868 // s11 -= tmp;
869 
870 // std::vector<double> ret(fdbins);
871 // for (size_t i = 0; i < fdbins; ++i) {
872 // ret[i] = sqrt(s11(i, i));
873 // }
874 
875 // // Clean up
876 // delete mx;
877 // delete predND;
878 // delete predFD;
879 // delete calc;
880 
881 // return ret;
882 
883 // } // function GetDecorrelatedFDErrors
884 
885 //===========================================================================================
886 // Function to pack covariance matrix to remove empty ND components prior to inversion
888 
889  size_t mxSize = mx.GetNrows();
890  size_t nComps = 13; // if this is a hardcoded hack, may as well go all in
891  std::vector<Binning> bins = { kNCNDBinning, kNCFDBinning };
892  size_t nSamples = bins.size();
893  size_t nBinsPerComp = kNCNDBinning.NBins() + kNCFDBinning.NBins();
894  std::vector<size_t> ndComps = { 0, 1, 8, 9, 12 };
895  size_t mxPackedSize = mxSize - ((8*bins[0].NBins())+1);
896  TMatrixD ret(mxPackedSize, mxPackedSize);
897  ret.Zero();
898  size_t iPack = 0;
899  // Loop over all components, samples, bins
900  for (size_t iComp = 0; iComp < nComps; ++iComp) {
901  for (size_t iSample = 0; iSample < nSamples; ++iSample) {
902  // If we're in an ND flux swap component, throw it away
903  if (iSample == 0 && std::find(ndComps.begin(), ndComps.end(), iComp) == ndComps.end()) {
904  // std::cout << "Skipping component " << iComp << ", sample " << iSample << std::endl;
905  continue;
906  }
907  for (size_t iBin = 0; iBin < (size_t)bins[iSample].NBins(); ++iBin) {
908  size_t i = (iComp*nBinsPerComp) + (iSample*kNCNDBinning.NBins()) + iBin;
909  if (i == 124) continue; // awful hack to remove problem bin
910  size_t jPack = 0;
911  for (size_t jComp = 0; jComp < nComps; ++jComp) {
912  for (size_t jSample = 0; jSample < nSamples; ++jSample) {
913  // If we're in an ND flux swap component, throw it away
914  if (jSample == 0 && std::find(ndComps.begin(), ndComps.end(), jComp) == ndComps.end()) {
915  // std::cout << "Skipping component " << jComp << ", sample " << jSample << std::endl;
916  continue;
917  }
918  for (size_t jBin = 0; jBin < (size_t)bins[jSample].NBins(); ++jBin) {
919  size_t j = (jComp*nBinsPerComp) + (jSample*kNCNDBinning.NBins()) + jBin;
920  if (j == 124) continue; // awful hack to remove problem bin
921  ret(iPack,jPack) = mx(i,j);
922  ++jPack;
923  } // for bin j
924  } // for sample j
925  } // for component j
926  ++iPack;
927  } // for bin i
928  } // for sample i
929  } // for component i
930 
931  return ret;
932 
933 } // function PackMatrix
934 
935 //===========================================================================================
936 // Unpack matrix to restore empty ND components
938 
939  size_t nComps = 13; // if this is a hardcoded hack, may as well go all in
940  std::vector<Binning> bins = { kNCNDBinning, kNCFDBinning };
941  size_t nSamples = bins.size();
942  size_t nBinsPerComp = kNCNDBinning.NBins() + kNCFDBinning.NBins();
943  std::vector<size_t> ndComps = { 0, 1, 8, 9, 12 };
944  size_t mxSize = nBinsPerComp * nComps;
945  TMatrixD ret(mxSize, mxSize);
946  ret.Zero();
947 
948  // Now if we packed the matrix, we need to unpack it!
949  size_t iPack = 0;
950  // Loop over all components, samples, bins
951  for (size_t iComp = 0; iComp < nComps; ++iComp) {
952  for (size_t iSample = 0; iSample < nSamples; ++iSample) {
953  // If we're in an ND flux swap component, throw it away
954  if (iSample == 0 && std::find(ndComps.begin(), ndComps.end(), iComp) == ndComps.end()) {
955  // std::cout << "Skipping component " << iComp << ", sample " << iSample << std::endl;
956  continue;
957  }
958  for (size_t iBin = 0; iBin < (size_t)bins[iSample].NBins(); ++iBin) {
959  size_t i = (iComp*nBinsPerComp) + (iSample*kNCNDBinning.NBins()) + iBin;
960  if (i == 124) continue; // awful hack to remove problem bin
961  size_t jPack = 0;
962  for (size_t jComp = 0; jComp < nComps; ++jComp) {
963  for (size_t jSample = 0; jSample < nSamples; ++jSample) {
964  // If we're in an ND flux swap component, throw it away
965  if (jSample == 0 && std::find(ndComps.begin(), ndComps.end(), jComp) == ndComps.end()) {
966  // std::cout << "Skipping component " << jComp << ", sample " << jSample << std::endl;
967  continue;
968  }
969  for (size_t jBin = 0; jBin < (size_t)bins[jSample].NBins(); ++jBin) {
970  size_t j = (jComp*nBinsPerComp) + (jSample*kNCNDBinning.NBins()) + jBin;
971  if (j == 124) continue; // awful hack to remove problem bin
972  ret(i,j) = mx(iPack,jPack);
973  ++jPack;
974  } // for bin j
975  } // for sample j
976  } // for component j
977  ++iPack;
978  } // for bin i
979  } // for sample i
980  } // for component i
981 
982  return ret;
983 
984 } // function UnpackMatrix
std::vector< covmx::Sample > GetNumuNDSamples()
Definition: Utilities.h:118
const double kNux20Th23
Definition: Utilities.h:42
const covmx::Sample kNumuRHCFarDet(covmx::kCCNumu, covmx::kRHC, covmx::kFarDet)
const XML_Char * name
Definition: expat.h:151
TMatrixD UnpackMatrix(TMatrixD mx)
Definition: Utilities.h:937
std::vector< double > quantiles(TH1D *h)
Definition: absCal.cxx:528
static std::unique_ptr< FDPredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
double GetDelta(int i, int j) const
const double kNux20RHCLivetime
Definition: Utilities.h:53
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
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
const double kNux20RHCNDPOT
Definition: Utilities.h:50
void SetNus20Params(osc::OscCalcSterile *calc, std::string type="3flav")
Definition: Utilities.h:619
const covmx::Sample kNumuQ1FHCFarDet(covmx::kCCNumu, covmx::kFHC, covmx::kFarDet, covmx::kQ1)
void SetEverything(covmx::Sample &sample, bool systs=false, double res=0)
Definition: Utilities.h:525
std::vector< covmx::Sample > GetNCFDSamples()
Definition: Utilities.h:158
const double kNux20Th12
Definition: Utilities.h:40
A simple Gaussian constraint on an arbitrary IFitVar.
vector< const ISyst * > getNusAna2020AllSysts(covmx::Sample &s)
std::vector< covmx::Sample > GetNumuSamples()
Definition: Utilities.h:140
const char * p
Definition: xmltok.h:285
void SetPrediction(covmx::Sample &sample, bool systs=true)
Definition: Utilities.h:503
void SetPOT(covmx::Sample &sample)
Definition: Utilities.h:196
const Var kNus20Energy([](const caf::SRProxy *sr){if(sr->hdr.det!=caf::kNEARDET &&sr->hdr.det!=caf::kFARDET) return(double) sr->slc.calE;double pars[4][6]={{1.049, 0.795, 0.8409, 0.17, 0.82,-1.00},{1.025, 0.797, 0.9162, 0.53,-0.26,-1.48},{1.000, 1.000, 1.0000, 0.00, 0.00, 0.00},{1.000, 1.000, 1.0000, 0.00, 0.00, 0.00}};int detCur=(sr->hdr.det==caf::kFARDET)+((sr->spill.isRHC)<< 1);double e_EM=ana::kEME_2020(sr);double e_Had=ana::kHADE_2020(sr);double e_Cal=sr->slc.calE;return(e_EM/pars[detCur][0]+e_Had/pars[detCur][1])/(pars[detCur][2]+pars[detCur][3]*std::pow(e_Cal+pars[detCur][4], 2)*std::exp(pars[detCur][5]*e_Cal));})
Definition: NusVars.h:64
bool CheckOption(TString opts, TString opt)
Definition: Utilities.h:47
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
void SetDelta(int i, int j, double delta)
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
void SetLivetime(double l)
Definition: Sample.h:53
const covmx::Sample kNumuRHCNearDet(covmx::kCCNumu, covmx::kRHC, covmx::kNearDet)
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
OStream cerr
Definition: OStream.cxx:7
std::vector< Spectrum * > LoadFakeData(std::vector< covmx::Sample > samples, size_t universe, bool systs)
Load fake data from file.
Definition: Utilities.h:736
const double kNux20Th13
Definition: Utilities.h:41
void SetInfo(covmx::Sample &sample, double res=0)
Definition: Utilities.h:495
void AddMatrix(CovarianceMatrix *gen)
std::string GetTag() const
Definition: Sample.cxx:94
string filename
Definition: shutoffs.py:106
Float_t tmp
Definition: plot.C:36
void SetCut(covmx::Sample &sample)
Definition: Utilities.h:261
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
const double kNux20RHCFDPOT
Definition: Utilities.h:51
const HistAxis kNumuCCOptimisedAxis2020("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kNumuE2020)
Definition: HistAxes.h:25
osc::OscCalcDumb calc
const covmx::Sample kNumuQ2FHCFarDet(covmx::kCCNumu, covmx::kFHC, covmx::kFarDet, covmx::kQ2)
covmx::Sample GetSampleFromOptString(TString optString)
Function to take an option TString and return a single associated covmx::Sample.
Definition: Utilities.h:342
const Cut kNue2020ND
Definition: NueCuts2020.h:178
#define M_PI
Definition: SbMath.h:34
const double kNux20FHCFDPOT
Definition: Utilities.h:49
const covmx::Sample kNumuQ4FHCFarDet(covmx::kCCNumu, covmx::kFHC, covmx::kFarDet, covmx::kQ4)
void SetCut(const Cut *c)
Definition: Sample.h:51
const Cut kNus20NDCuts
Definition: NusCuts20.h:102
const HistAxis kNC20FDAxisE("Energy Deposited in Scintillator (GeV)", kNCFDBinning, kNus20Energy)
Definition: HistAxes.h:21
const XML_Char * s
Definition: expat.h:262
std::string GetPolarity() const
Definition: Sample.cxx:52
void SetFakeSignalParams(osc::OscCalcSterile *calc, size_t type)
Definition: Utilities.h:684
double GetDefaultPOT(covmx::Sample sample)
Definition: Utilities.h:185
const double kNux20FHCLivetime
Definition: Utilities.h:52
void SetAxis(const HistAxis *a)
Definition: Sample.h:50
T ParseOption(std::vector< std::pair< std::string, T >> opts, TString optString, T def)
Definition: Utilities.h:82
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
Definition: Utilities.h:416
string infile
static std::unique_ptr< NDPredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
const Cut sels[kNumSels]
Definition: vars.h:44
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
std::vector< bool > GetOptionConfig(TString opt, std::vector< std::string > options)
Definition: Utilities.h:56
const std::string kNus20Path
void SetPrediction(IPrediction *p)
Definition: Sample.cxx:103
std::string getenv(std::string const &name)
Selection selection
Definition: Sample.h:95
const HistAxis kNue2020Axis("NuE Energy / Analysis Bin", kNue2020Binning, kNue2020AnaBin)
Use this Axis for Ana2020, official Axis.
Definition: NueCuts2020.h:195
Detector detector
Definition: Sample.h:97
void SetCosmic(covmx::Sample &sample)
Set cosmics in a given sample.
Definition: Utilities.h:488
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
const double kAna2020FHCLivetime
Definition: Exposures.h:237
const double kNux20Delta24
Definition: Utilities.h:46
void PrintOscParams(osc::OscCalcSterile *calc)
Definition: Utilities.h:707
void SetCosmic(Spectrum *s)
Definition: Sample.cxx:108
const std::vector< std::string > ndComps
Definition: TestCovMxNew.C:12
const HistAxis * GetDefaultAxis(covmx::Sample sample, double res=0)
Definition: Utilities.h:215
const double j
Definition: BetheBloch.cxx:29
Eigen::VectorXd vec
std::vector< covmx::Sample > GetNCSamples()
Definition: Utilities.h:167
double GetDefaultLivetime(covmx::Sample sample)
Definition: Utilities.h:201
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
const double kNux20Dm32Sigma
Definition: Utilities.h:38
const double kAna2020FHCPOT
Definition: Exposures.h:233
static Binning Custom(const std::vector< double > &edges)
Definition: Binning.cxx:145
float bin[41]
Definition: plottest35.C:14
void SetLivetime(covmx::Sample &sample)
Definition: Utilities.h:210
const double kAna2020RHCPOT
Definition: Exposures.h:235
std::vector< covmx::Sample > GetAllSamples()
Definition: Utilities.h:176
OStream cout
Definition: OStream.cxx:6
const covmx::Sample kNumuQ2FHCNearDet(covmx::kCCNumu, covmx::kFHC, covmx::kNearDet, covmx::kQ2)
string InputPath()
TMatrixD PackMatrix(TMatrixD mx)
Definition: Utilities.h:887
const std::string path
Definition: plot_BEN.C:43
const Binning bins
Definition: NumuCC_CPiBin.h:8
void SetAngle(int i, int j, double th)
const HistAxis kNC20NDAxisE("Energy Deposited in Scintillator (GeV)", kNCNDBinning, kNus20Energy)
NC covariance hist axes.
Definition: HistAxes.h:20
const Binning kNCNDBinning
NC custom binning.
Definition: Binning.cxx:104
const covmx::Sample kNumuQ1FHCNearDet(covmx::kCCNumu, covmx::kFHC, covmx::kNearDet, covmx::kQ1)
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
void SetName(std::string name)
void SetDm(int i, double dm)
TDirectory * dir
Definition: macro.C:5
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
Quantile quantile
Definition: Sample.h:98
const double kAna2020RHCLivetime
Definition: Exposures.h:238
int NBins() const
Definition: Binning.h:29
const covmx::Sample kNumuQ3FHCNearDet(covmx::kCCNumu, covmx::kFHC, covmx::kNearDet, covmx::kQ3)
std::vector< T > ParseOptions(std::vector< std::pair< std::string, T >> opts, TString optString)
Definition: Utilities.h:105
std::string GetName() const
Definition: Sample.cxx:72
const Binning kNCFDBinning
Definition: Binning.cxx:105
std::vector< covmx::Sample > GetSamplesFromOptString(TString optString)
Function to take an option TString and return a vector of associated covmx::Samples.
Definition: Utilities.h:371
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
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)
const double kNux20Th23Sigma
Definition: Utilities.h:43
double GetDm(int i) const
void SetPOT(double d)
Definition: Sample.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
static std::unique_ptr< CovarianceMatrix > LoadFrom(TDirectory *dir, const std::string &name)
const covmx::Sample kNumuQ3FHCFarDet(covmx::kCCNumu, covmx::kFHC, covmx::kFarDet, covmx::kQ3)
const double kNux20Dm32
Definition: Utilities.h:37
covmx::CovarianceMatrix * GetFullMatrix(bool cvmfs=true)
Definition: Utilities.h:605
double T
Definition: Xdiff_gwt.C:5
const Cut kNue2020FDAllSamples_ML
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
const Cut kNus20FDCuts_ML
Polarity polarity
Definition: Sample.h:96
std::vector< covmx::Sample > GetNCNDSamples()
Definition: Utilities.h:149
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
const Cut kNumu2020FD_ML
void next()
Definition: show_event.C:84
const double kNux20Delta13
Definition: Utilities.h:45
const double kNux20FHCNDPOT
Definition: Utilities.h:48
void AssertSupported(covmx::Sample sample)
Definition: Utilities.h:461
void SetCuts(std::vector< covmx::Sample > &samples)
Definition: Utilities.h:296
const covmx::Sample kNumuFHCNearDet(covmx::kCCNumu, covmx::kFHC, covmx::kNearDet)
double GetAngle(int i, int j) const
const double kNux20Dm21
Definition: Utilities.h:36
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...
const covmx::Sample kNumuFHCFarDet(covmx::kCCNumu, covmx::kFHC, covmx::kFarDet)
const covmx::Sample kNumuQ4FHCNearDet(covmx::kCCNumu, covmx::kFHC, covmx::kNearDet, covmx::kQ4)
const FitTheta23Sterile kFitTheta23Sterile
void SetAxis(covmx::Sample &sample, double res=0)
Definition: Utilities.h:255
std::vector< covmx::Sample > GetNumuFDSamples()
Definition: Utilities.h:129
enum BeamMode string