Nus18Utilities.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <vector>
4 #include <tuple>
5 
8 #include "NuXAna/Analysis/NusSystsMaker.h"
10 #include "CAFAna/Core/SystShifts.h"
11 #include "CAFAna/Core/Sample.h"
19 #include "Utilities/rootlogon.C"
20 
21 #include "TLegend.h"
22 #include "TLatex.h"
23 #include "TLine.h"
24 #include "TGraph.h"
25 #include "THStack.h"
26 
27 using namespace ana;
28 
29 // Hard-coded parameters
30 const double kNus18Dm21 = 7.53e-5;
31 const double kNus18Dm32 = 2.44e-3;
32 const double kNus18Dm32Sigma = 0.08e-3;
33 
34 const double kNus18Th12 = 0.587;
35 const double kNus18Th13 = 0.145;
36 const double kNus18Th23 = 0.844;
37 const double kNus18Th23Sigma = 0.042;
38 
39 const double kNus18Delta13 = 1.21*M_PI;
40 const double kNus18Delta24 = M_PI/2;
41 
42 const int kOfficialNus18FHCColor = kAzure+2;
43 const int kOfficialNus18RHCColor = kRed+2;
44 
45 // Hard-coded loaders
46 // cosmics
47 const std::string kFDCosmicFHCProd4 = "prod_caf_R17-11-14-prod4reco.b_fd_cosmic_fhcTune_full_v1_goodruns";
48 const std::string kFDCosmicRHCProd4 = "prod_caf_R17-11-14-prod4reco.a_fd_cosmic_rhcTune_HighGain_v1_goodruns_snapshot_170116";
49 
50 // data
51 const std::string kFDDataFHC = "prod4_sumrestricteddecaf_fd_numi_fhc_full_goodruns_nus2018";
52 const std::string kFDDataRHC = "prod4_sumrestricteddecaf_fd_numi_rhc_full_goodruns_nus2018";
53 
54 // mc (in general use loaders, just for specific purposes)
55 const std::string kFDMCFHCProd4 = "prod_caf_R17-11-14-prod4reco.neutron-respin.c_fd_genie_nonswap_fhc_nova_v08_full_v1";
56 const std::string kFDMCRHCProd4 = "prod_caf_R17-11-14-prod4reco.neutron-respin.c_fd_genie_nonswap_rhc_nova_v08_full_v1";
57 
58 //===========================================================================================
59 void AddBin(std::vector<double>& v, double b)
60 { v.push_back(v.back()+b); }
61 
62 //===========================================================================================
64 {
65  std::vector<double> v = { 0. };
66  for (size_t i = 0; i < 6; ++i)
67  AddBin(v, 1);
68  for (size_t i = 0; i < 2; ++i)
69  AddBin(v, 2);
70  AddBin(v, 10);
71 
72  // AddBin(v,0.75);
73  // for (int i = 0; i < 37; ++i)
74  // AddBin(v,0.25);
75  // for (int i = 0; i < 8; ++i)
76  // AddBin(v,0.5);
77  // for (int i = 0; i < 6; ++i)
78  // AddBin(v,1);
79 
80  return Binning::Custom(v);
81 
82  /* AddBin(v, 0.75);
83  for (int i = 0; i < 37; ++i) AddBin(v, 0.25);
84  for (int i = 0; i < 2; ++i) AddBin(v, 0.5);
85  for (int i = 0; i < 2; ++i) AddBin(v, 1);
86  AddBin(v, 2);
87  AddBin(v, 5);
88  */
89 }
90 
91 //===========================================================================================
93 {
94  std::vector<double> v = { 0. };
95  AddBin(v,1.25);
96  AddBin(v,0.25);
97  AddBin(v,0.25);
98  AddBin(v,0.5);
99  AddBin(v,1);
100  AddBin(v,1.5);
101  AddBin(v,2);
102  AddBin(v,3.25);
103  AddBin(v,10);
104 
105  return Binning::Custom(v);
106 }
107 
108 //===========================================================================================
110 {
111  std::vector<double> v = { 0. };
112  AddBin(v,2.25);
113  AddBin(v,3.25);
114  AddBin(v,14.5);
115 
116  return Binning::Custom(v);
117 }
118 
119 //===========================================================================================
121 {
122  std::vector<double> v = { 0. };
123  AddBin(v,0.75);
124  for (int i = 0; i < 67; ++i)
125  AddBin(v,0.25);
126  AddBin(v,0.5);
127  AddBin(v,2.0);
128 
129  return Binning::Custom(v);
130 }
131 
132 //===========================================================================================
133 std::vector<Binning> GetNus18Binning(bool rhc, std::string detector)
134 {
135 
136  // Check detector
137 
138  bool use_nd = false;
139  bool use_fd = false;
140  if (detector == "nd") use_nd = true;
141  else if (detector == "fd") use_fd = true;
142  else if (detector == "both") {
143  use_nd = true;
144  use_fd = true;
145  }
146  else throw std::runtime_error("Detector option \"" + detector + "\" not understood.");
147 
148  // Now actually add bins
149 
150  std::vector<Binning> ret;
151 
152  if (rhc) {
153  if (use_nd) ret.push_back(RHCFDBins());
154  if (use_fd) ret.push_back(RHCFDBins());
155  } else {
156  if (use_nd) ret.push_back(FHCNDBins());
157  if (use_fd) ret.push_back(FHCFDBins());
158  }
159 
160  return ret;
161 
162 } // function GetNus18Binning
163 
164 //===========================================================================================
166 {
167  // Default vals
169 
170  // Mass splittings
171  calc->SetDm(2, kNus18Dm21);
172  calc->SetDm(3, kNus18Dm21 + kNus18Dm32);
173  calc->SetDm(4, 0.5);
174 
175  // Angles
176  calc->SetAngle(1, 3, kNus18Th13);
177  calc->SetAngle(2, 3, kNus18Th23);
178  calc->SetAngle(3, 4, 0.);
179  calc->SetAngle(2, 4, 0.);
180  calc->SetAngle(1, 4, 0.);
181 
182  // Phases
183  calc->SetDelta(1, 3, kNus18Delta13);
184  calc->SetDelta(2, 4, kNus18Delta24);
185 }
186 
187 //===========================================================================================
189 {
190  // Mass splittings
191  std::cout << "---------------------------" << std::endl;
192  std::cout << " OSCILLATION PARAMETERS " << std::endl;
193  std::cout << std::endl;
194  std::cout << " dmsq 21 - " << calc->GetDm(2) << " eV^2" << std::endl;
195  std::cout << " dmsq 31 - " << calc->GetDm(3) << " eV^2" << std::endl;
196  std::cout << " dmsq 41 - " << calc->GetDm(4) << " eV^2" << std::endl;
197  std::cout << std::endl;
198  std::cout << " th 12 - " << calc->GetAngle(1, 2) / M_PI * 180 << " deg" << std::endl;
199  std::cout << " th 13 - " << calc->GetAngle(1, 3) / M_PI * 180 << " deg" << std::endl;
200  std::cout << " th 23 - " << calc->GetAngle(2, 3) / M_PI * 180 << " deg" << std::endl;
201  std::cout << " th 14 - " << calc->GetAngle(1, 4) / M_PI * 180 << " deg" << std::endl;
202  std::cout << " th 24 - " << calc->GetAngle(2, 4) / M_PI * 180 << " deg" << std::endl;
203  std::cout << " th 34 - " << calc->GetAngle(3, 4) / M_PI * 180 << " deg" << std::endl;
204  std::cout << std::endl;
205  std::cout << " delta 13 - " << calc->GetDelta(1, 3) / M_PI * 180 << " deg" << std::endl;
206  std::cout << " delta 14 - " << calc->GetDelta(1, 4) / M_PI * 180 << " deg" << std::endl;
207  std::cout << " delta 24 - " << calc->GetDelta(2, 4) / M_PI * 180 << " deg" << std::endl;
208  std::cout << "---------------------------" << std::endl;
209 }
210 
211 //===========================================================================================
212 /// Canvas formatting utility
214 {
215  gStyle->SetCanvasDefW(1200);
216  gStyle->SetCanvasDefH(800);
217  gStyle->SetPadTopMargin(0.14);
218  gStyle->SetPadBottomMargin(0.14);
219  gStyle->SetPadLeftMargin(0.14);
220  gStyle->SetPadRightMargin(0.14);
221 }
222 
223 //===========================================================================================
224 /// Function to load prediction object
226  std::string syst_type = "all")
227 {
228  DontAddDirectory guard;
229 
230  if (detector != "neardet" && detector != "fardet")
231  throw std::runtime_error("Detector string must be either \"neardet\" or \"fardet\"!");
232 
233  std::string dname = detector == "neardet" ? "nd" : "fd";
234  std::string polarity = rhc? "rhc" : "fhc";
235 
236  // log400 hack for the time being
237  if (not rhc and detector == "neardet" and syst_type == "none") {
238 
239  std::string filename = FindCAFAnaDir() + "/data/log400/pred_nc_nd_log400.root";
240  return LoadFromFile<PredictionInterp>(filename.c_str(),
241  "pred_interp_nus_fhc_neardet_none").release();
242 
243  }
244 
245  std::string filename = "pred_nus18_" + polarity + "_" + dname
246  + "_" + syst_type + ".root";
247 
248  // std::string filename = "pred_nus_" + polarity + "_" + detector
249  // + "_" + syst_type + ".root";
250 
251  char* path = std::getenv("NUSDATA_NUS18_PRED");
252  if (!path) assert(false and
253  "NUSDATA_NUS18_PRED environment variable not set! Did you forget to set up nusdata?");
254 
255  // std::string path = FindCAFAnaDir() + "/data/joint/pred";
256  std::cout << "Loading prediction from " << path << "/" << filename << std::endl;
257 
258  std::string filepath = std::string(path) + "/" + filename;
259  return LoadFromFile<PredictionInterp>(filepath, "pred_interp").release();
260 
261 } // function LoadPrediction
262 
263 //===========================================================================================
265 {
266  DontAddDirectory guard;
267 
268  bool fhc = opt.Contains("fhc", TString::kIgnoreCase);
269  bool rhc = opt.Contains("rhc", TString::kIgnoreCase);
270  assert(fhc xor rhc);
271 
272  bool syst = opt.Contains("syst", TString::kIgnoreCase);
273 
274  std::stringstream filename;
275  filename << "pred_nus18_";
276  rhc ? filename << "rhc" : filename << "fhc";
277  filename << "_extrap";
278  if (syst) filename << "_systs";
279  filename << ".root";
280 
281  char* path = std::getenv("NUSDATA_NUS18_PRED");
282  if (!path) throw std::runtime_error(
283  "NUSDATA_NUS18_PRED environment variable not set! Did you forget to set up nusdata?");
284  else std::cout << "Loading prediction from " << path << "/" << filename.str() << std::endl;
285 
286  std::string filepath = std::string(path) + "/" + filename.str();
287 
288  //return LoadFromFile<IPrediction>(filepath, "pred_systs").release();
289  return LoadFromFile<IPrediction>(filepath, "FDNus18ExtrapPred").release();
290 
291 } // function LoadPrediction
292 
293 //===========================================================================================
294 /// Function to load covariance matrix
296 {
297  DontAddDirectory guard;
298 
299  char* path = std::getenv("NUSDATA_NUS18_COVMX");
300  if (!path) throw std::runtime_error(
301  "NUSDATA_NUS18_COVMX environment variable not set! Did you forget to set up nusdata?");
302 
303  std::string filepath = std::string(path) + "/" + file + ".root";
304  //std::string filepath = FindCAFAnaDir() + "data/joint/covmx/" + file + ".root";
305  std::cout << "Loading covariance matrix " << name << " from " << filepath << std::endl;
306 
307  return LoadFromFile<CovarianceMatrix>(filepath, name).release();
308 }
309 
310 //===========================================================================================
311 /// Function to load PPFX universe simulations
313 {
314  DontAddDirectory guard;
315 
316  char* path = std::getenv("NUSDATA_NUS18_WEIGHTS");
317  if (detector != "nd" && detector != "fd") throw std::runtime_error(
318  "Detector string must be \"nd\" or \"fd\" - " + detector + " not understood!");
319  if (!path) throw std::runtime_error(
320  "NUSDATA_NUS18_WEIGHTS environment variable not set! Did you forget to set up nusdata?");
321  std::string file = std::string(path) + "/nus18_fhc_" + detector + "_ppfx.root";
322  return LoadFromFile<PredictionInterp>(file, "universe_" + std::to_string(universe)).release();
323 }
324 
325 //===========================================================================================
326 /// Function to load the unblinded data spectrum
327 /*Spectrum* LoadData(bool rhc) {
328 
329  DontAddDirectory guard;
330 
331  if (std::getenv("PROCESS")) {
332  std::cout << "I see you're trying to run on the grid."
333  << "Now may be a good time to figure out how to read data on the grid..." << std::endl;
334  abort();
335  }
336 
337  std::string path = "/nova/app/users/wallbank/test_ndosc/NuXAna/macros/Nus18/box_opening/";
338 
339  return LoadFromFile<Spectrum>
340  (Form("%s/fout_nus18_box_opening_%s.root", path.c_str(), rhc? "rhc" : "fhc"), "spec_nus_vise_numi").release();
341 
342 }*/
343 
344 //===========================================================================================
345 /// Function to load fake data file
346 // Spectrum LoadFakeData(std::string filename, int universe)
347 // {
348 // DontAddDirectory guard;
349 
350 // char* path = std::getenv("NUSDATA_NUS18_FAKEDATA");
351 // if (!path) throw std::runtime_error(
352 // "NUSDATA_NUS18_FAKEDATA environment variable not set! Did you forget to set up nusdata?");
353 // else std::cout << "Loading fake data spectrum universe_" << universe << " from "
354 // << path << "/" << filename << ".root" << std::endl;
355 
356 // std::string filepath = std::string(path) + "/" + filename + ".root";
357 // return *(LoadFromFile<Spectrum>(filepath, "universe_" + std::to_string(universe)).release());
358 // }
359 
360 Spectrum LoadFakeData(std::string filename, int universe, double pot, double livetime)
361 {
362  DontAddDirectory guard;
363 
364  char* path = std::getenv("NUSDATA_NUS18_FAKEDATA");
365  if (!path) throw std::runtime_error(
366  "NUSDATA_NUS18_FAKEDATA environment variable not set! Did you forget to set up nusdata?");
367  else std::cout << "Loading fake data spectrum universe_" << universe << " from "
368  << path << "/" << filename << ".root" << std::endl;
369 
370  std::string filepath = std::string(path) + "/" + filename + ".root";
371  std::string histname = "universe_" + std::to_string(universe);
372 
373  TFile* f = new TFile(filepath.c_str(), "read");
374  TH1D* h = (TH1D*)f->Get(histname.c_str());
375 
376  return Spectrum(h, pot, livetime);
377 }
378 
379 //===========================================================================================
380 /// Function to load NuS data file
381 Spectrum LoadData(bool nd, bool rhc = false)
382 {
383  DontAddDirectory guard;
384 
385  std::string path = std::string(std::getenv("SRT_PRIVATE_CONTEXT"));
386  if (rhc) path = "/nova/app/users/wallbank/test_ndosc/NuXAna/macros/Nus18_neutrino/box_opening/";
387 
388  if (rhc) return *LoadFromFile<Spectrum>(path + "/fout_nus18_box_opening_rhc.root", "spec_nus_vise_numi").release();
389 
390  if (nd) return *LoadFromFile<Spectrum>(path + "/files/fhc_NDDataMCLoadOut.root","RecoEPreselPtpData").release();
391  else return *LoadFromFile<Spectrum>(path + "/files/fout_nus18_box_opening_fhc.root", "spec_nus_vise_numi").release();
392 }
393 
394 //===========================================================================================
395 /// Function to load cosmic histogram
396 TH1D* LoadCosmicHist(bool rhc = false)
397 {
398  DontAddDirectory guard;
399 
400  char* path = std::getenv("NUSDATA_NUS18_COSMIC");
401  if (!path) throw std::runtime_error(
402  "NUSDATA_NUS18_COSMIC environment variable not set! Did you forget to set up nusdata?");
403  else std::cout << "Loading cosmic data from " << path << "/nus18_cosmic_background.root" << std::endl;
404 
405  TFile f(Form("%s/nus18_cosmic_background.root", path));
406  return (TH1D*)f.Get(rhc? "rhcnumi" : "fhcnumi");
407 }
408 
409 //===========================================================================================
410 /// Function to load cosmic histogram
411 TH1D* LoadCosmicHist(TString opt)
412 {
413  DontAddDirectory guard;
414 
415  bool fhc = opt.Contains("fhc", TString::kIgnoreCase);
416  bool rhc = opt.Contains("rhc", TString::kIgnoreCase);
417  assert(fhc || rhc);
418 
419  char* path = std::getenv("NUSDATA_NUS18_COSMIC");
420  if (!path) throw std::runtime_error(
421  "NUSDATA_NUS18_COSMIC environment variable not set! Did you forget to set up nusdata?");
422  else std::cout << "Loading cosmic data from " << path << "/nus18_cosmic_background.root" << std::endl;
423 
424  TFile f(Form("%s/nus18_cosmic_background.root", path));
425  return (TH1D*)f.Get(rhc? "rhcnumi" : "fhcnumi");
426 }
427 
428 //===========================================================================================
429 /// Function to return cosmic spectrum
430 Spectrum LoadCosmicSpectrum(bool rhc = false)
431 {
432  DontAddDirectory guard;
434 }
435 
436 //===========================================================================================
437 /// Function to return cosmic spectrum
439 {
440  DontAddDirectory guard;
441 
442  bool fhc = opt.Contains("fhc", TString::kIgnoreCase);
443  bool rhc = opt.Contains("rhc", TString::kIgnoreCase);
444  assert(fhc || rhc);
445 
447 }
448 
449 /* //=========================================================================================== */
450 /* /// Function to load systematics */
451 /* NusSystematics* LoadSystematics() { */
452 
453 /* char* path = std::getenv("NUSDATA_NUS18_SYSTS"); */
454 /* if (!path) throw std::runtime_error("NUSDATA_NUS18_SYSTS environment variable not set!"); */
455 /* else std::cout << "Loading systematics from " << path << "/Nus18Systs.root" << std::endl; */
456 
457 /* TFile f(Form("%s/Nus18Systs.root", path)); */
458 /* NusSystematics* nus_systs = NusSystematics::LoadFrom(f.GetDirectory("Nus18Systs")); */
459 
460 /* return nus_systs; */
461 
462 /* } */
463 
464 //===========================================================================================
465 /// Function to load ISysts
466 std::vector<const ISyst*> LoadISysts() {
467 
468  DontAddDirectory guard;
469 
470  char* path = std::getenv("NUSDATA_NUS18_SYSTS");
471  if (!path)
472  throw std::runtime_error("NUSDATA_NUS18_SYSTS environment variable not set!");
473  else
474  std::cout << "Loading systematics from " << path << "/Nus18ISysts.root" << std::endl;
475 
476  std::vector<const ISyst*> systs;
477 
478  TFile* f = new TFile(Form("%s/Nus18ISysts.root", path), "READ");
479  TIter next(f->GetListOfKeys());
480  TKey* key;
481  while ((key = (TKey*)next()))
482  systs.push_back(NusISyst::LoadFrom(f, key->GetName())).release();
483  f->Close();
484  delete f;
485 
486  return systs;
487 
488 }
489 
490 //===========================================================================================
491 /// Function to get location of PPFX weights
493 {
494  char* ppfx_env = std::getenv("NUSDATA_NUS18_WEIGHTS");
495  if (ppfx_env == nullptr) throw std::runtime_error(
496  "NUSDATA_NUS18_WEIGHTS environment variable not set! Did you forget to set up nusdata?");
497 
498  return std::string(ppfx_env);
499 }
500 /*
501 //===========================================================================================
502 /// Function to get default Nus18 simulation
503 PredictionConcat* GetDefaultSimulation(bool rhc, std::string detector, std::string syst_type)
504 {
505  // Detector type
506  bool use_nd = false;
507  bool use_fd = false;
508  if (detector == "nd") use_nd = true;
509  else if (detector == "fd") use_fd = true;
510  else if (detector == "both") {
511  use_nd = true;
512  use_fd = true;
513  }
514  else throw std::runtime_error("Detector string must be \"nd\", \"fd\" or \"both\".");
515 
516  // Systematic option
517  if (syst_type != "none" && syst_type != "xsec" && syst_type != "nonxsec")
518  throw std::runtime_error("Syst type string must be \"none\", \"xsec\" or \"nonxsec\".");
519 
520  // Get exposures
521  double nd_pot = rhc ? kAna2018SensitivityRHCNDPOT : kAna2018SensitivityFHCNDPOT;
522  double fd_pot = rhc ? kAna2018RHCPOT : kAna2018FHCPOT;
523  double fd_livetime = rhc ? kAna2018RHCLivetime : kAna2018FHCLivetime;
524  std::string polarity = rhc? "RHC" : "FHC";
525 
526  // Vectors for predictions, POT & livetimes, systematic prefixes
527  std::vector<IPrediction*> predictions;
528  std::vector<double> pots;
529  std::vector<double> livetimes;
530  std::vector<std::string> samples;
531  std::vector<std::vector<const ISyst*>> systs;
532 
533  // Fill
534  if (use_nd) {
535  samples.push_back("Near detector");
536  predictions.push_back(LoadPrediction("neardet", rhc, syst_type));
537  pots.push_back(nd_pot);
538  livetimes.push_back(0);
539  systs.push_back(GetNus18Systs(rhc, "nd", syst_type));
540  // syst_types.push_back(polarity + "_ND");
541  }
542  if (use_fd) {
543  samples.push_back("Far detector");
544  predictions.push_back(LoadPrediction("fardet", rhc, syst_type));
545  pots.push_back(fd_pot);
546  livetimes.push_back(fd_livetime);
547  systs.push_back(GetNus18Systs(rhc, "fd", syst_type));
548  // syst_types.push_back(polarity + "_FD");
549  }
550 
551  PredictionConcat* sim = new PredictionConcat(predictions, pots, livetimes);
552  sim->SetSysts(new SystConcat(samples, {}, systs));
553  sim->SetBinning(GetNus18Binning(rhc, detector));
554 
555  return sim;
556 }
557 */
558 //===========================================================================================
559 /// Function to get default cosmics
560 std::vector<TH1D*> GetCosmics(bool rhc, std::string detector)
561 {
562  bool use_nd = false;
563  bool use_fd = false;
564  if (detector == "nd") use_nd = true;
565  else if (detector == "fd") use_fd = true;
566  else if (detector == "both") {
567  use_nd = true;
568  use_fd = true;
569  }
570  else throw std::runtime_error("Detector string \"" + detector + "\" not recognised!");
571 
572  std::vector<TH1D*> cosmics;
573  if (use_nd) cosmics.push_back(MakeTH1D(UniqueName().c_str(), "", kNus18EnergyBinning));
574  if (use_fd) cosmics.push_back(LoadCosmicHist(rhc));
575 
576  return cosmics;
577 }
578 
579 //===========================================================================================
580 /// Function to combine covariance matrices
581 void CombineMatrices(CovarianceMatrix* gen,
582  std::map<const ISyst*, CovarianceMatrix*> m, std::vector<IPrediction*> preds,
583  std::vector<double> pot, osc::IOscCalcAdjustable* calc)
584 {
585  if (m.size() == 0) throw std::runtime_error("Matrix map is empty!");
586  else if (m.size() == 1) throw std::runtime_error("Only one matrix in map; cannot add!");
587  for (auto matrix : m) gen->AddMatrix(matrix.second);
588  gen->PredictCovMx(preds, pot, calc);
589 }
590 
591 //===========================================================================================
592 /// Function to get default full covariance matrix
593 /*CovarianceMatrix* GetFullCovMx(PredictionConcat* sim,
594  osc::OscCalcSterile* calc, double pot, std::vector<bool> use_detector)
595 {
596  DontAddDirectory guard;
597 
598  std::string detector_name;
599 
600  if (use_detector[0] && use_detector[1]) detector_name = "Both";
601  else if (use_detector[1]) detector_name = "FD";
602  else if (use_detector[0]) detector_name = "ND";
603 
604  // Load all the covariance matrices
605 
606  std::map<const ISyst*, CovarianceMatrix*> m;
607  for (const ISyst* syst : GetNus18Systs(false, "base", "all")) {
608  std::string systname = syst->ShortName().substr(syst->ShortName().find_last_of("_")+1);
609  CovarianceMatrix* gen = LoadCovMx("covmx_nc_fhc", "CovMx_" + detector_name + "_" + systname);
610  m.insert(std::pair<const ISyst*, CovarianceMatrix*>(syst, gen));
611  }
612  CovarianceMatrix* ppfx = LoadCovMx("covmx_nc_fhc", "CovMx_" + detector_name + "_PPFX");
613  m.insert(std::pair<const ISyst*, CovarianceMatrix*>(nullptr, ppfx));
614 
615  std::vector<Binning> b = sim->GetBinning();
616 
617  // Generate new empty matrix
618 
619  CovarianceMatrix* matrix_combined = new CovarianceMatrix(pot, b);
620 
621  // Add all the matrices together
622 
623  CombineMatrices(matrix_combined, m, sim, calc);
624 
625  // Clean up the individual matrices
626 
627  for (auto matrix : m) delete matrix.second;
628 
629  return matrix_combined;
630 }
631 */
632 //===========================================================================================
633 /// Get fake data for a set of samples
634 
635 
636 //===========================================================================================
637 /// Get cosmics for a set of samples
638 std::vector<TH1D*> GetCosmics(std::vector<covmx::Sample> samples,
639  PredictionConcat* sim)
640 {
641  DontAddDirectory guard;
642 
643  std::map<covmx::Polarity, std::string> pname;
644  pname[covmx::Polarity::kFHC] = "fhc";
645  pname[covmx::Polarity::kRHC] = "rhc";
646 
647  std::map<covmx::Quantile, std::string> qname;
648  qname[covmx::Quantile::kQ1] = "q1";
649  qname[covmx::Quantile::kQ2] = "q2";
650  qname[covmx::Quantile::kQ3] = "q3";
651  qname[covmx::Quantile::kQ4] = "q4";
652 
653  assert(samples.size() == sim->GetNSamples() and
654  "Sample vector and concatenated prediction have different size!");
655 
656  // Open only the cosmic files we need, with a view to efficiency when running on grid
657 
658  std::string path = "/pnfs/nova/persistent/users/jhewes15/cosmic/";
659  std::vector<TFile*> f(4, nullptr);
660  for (const covmx::Sample& s : samples) {
661 
662  // NC cosmics
663 
664  if (s.analysis == covmx::Analysis::kNC and !f[0]) {
665  std::string filepath = path + "cosmics_nus.root";
666  if (RunningOnGrid()) filepath = pnfs2xrootd(filepath);
667  f[0] = TFile::Open(filepath.c_str(), "read");
668  }
669 
670  // CCNue cosmics
671 
672  if (s.analysis == covmx::Analysis::kCCNue && !f[1]) {
673  std::string filepath = path + "cosmics_nue.root";
674  if (RunningOnGrid()) filepath = pnfs2xrootd(filepath);
675  f[1] = TFile::Open(filepath.c_str(), "read");
676  }
677 
678  // CCNumu FHC cosmics
679 
680  if (s.analysis == covmx::Analysis::kCCNumu &&
681  s.polarity == covmx::Polarity::kFHC && !f[2]) {
682  std::string filepath = path + "cosmics_fhc_numu.root";
683  if (RunningOnGrid()) filepath = pnfs2xrootd(filepath);
684  f[2] = TFile::Open(filepath.c_str(), "read");
685  }
686 
687  // CCNumu RHC cosmics
688 
689  if (s.analysis == covmx::Analysis::kCCNumu &&
690  s.polarity == covmx::Polarity::kRHC && !f[3]) {
691  std::string filepath = path + "cosmics_rhc_numu.root";
692  if (RunningOnGrid()) filepath = pnfs2xrootd(filepath);
693  f[3] = TFile::Open(filepath.c_str(), "read");
694  }
695 
696  }
697 
698  // Now actually get the cosmic spectra
699 
700  std::vector<TH1D*> ret(samples.size(), nullptr);
701 
702  for (size_t i = 0; i < samples.size(); ++i) {
703 
704  // If it's a near detector spectrum, we need an empty hist for cosmics
705 
706  if (samples[i].detector == covmx::Detector::kNearDet)
707  ret[i] = MakeTH1D(UniqueName().c_str(), "", sim->GetBinningOriginal()[i]);
708 
709  // Otherwise we need to load it. The cosmics are not currently standardised
710  // between analyses, so we need to handle each differently
711 
712  else {
713 
714  // NC cosmics
715 
716  if (samples[i].analysis == covmx::Analysis::kNC) {
717  std::string hname = pname[samples[i].polarity] + "cos";
718  ret[i] = (TH1D*)f[0]->Get(hname.c_str());
719  }
720 
721  // CC nue cosmics
722 
723  else if (samples[i].analysis == covmx::Analysis::kCCNue) {
724  std::string sname = "cosmic_spect_" + pname[samples[i].polarity];
725  std::unique_ptr<Spectrum> s(Spectrum::LoadFrom(f[1], sname.c_str()));
726  ret[i] = s->ToTH1(sim->GetLivetimes()[i], kLivetime);
727  }
728 
729  // CC numu FHC cosmics
730 
731  else if (samples[i].analysis == covmx::Analysis::kCCNumu
732  && samples[i].polarity == covmx::Polarity::kFHC) {
733  std::string hname = "cosmics_" + qname[samples[i].quantile];
734  ret[i] = (TH1D*)f[2]->Get(hname.c_str());
735  }
736 
737  // CC numu RHC cosmics
738 
739  else if (samples[i].analysis == covmx::Analysis::kCCNumu
740  && samples[i].polarity == covmx::Polarity::kRHC) {
741  std::string hname = "cosmics_" + qname[samples[i].quantile];
742  ret[i] = (TH1D*)f[3]->Get(hname.c_str());
743  }
744 
745  // By this point we should have instantiated the cosmic spectrum
746 
747  assert(ret[i] && "Found uninitialised cosmic spectrum!");
748 
749  } // if far detector
750  } // for sample
751 
752  // Clean up files
753 
754  for (TFile* tmp : f) {
755  if (tmp) delete tmp;
756  }
757 
758  return ret;
759 
760 } // function GetCosmics
761 
762 //===========================================================================================
763 /// Get systematics for a given sample
764 std::vector<const ISyst*> GetSystsFromFile(covmx::Sample sample)
765 {
766  // Only supported for NC right now
767 
768  if (sample.analysis != covmx::Analysis::kNC or sample.polarity != covmx::Polarity::kFHC)
769  assert(false and "Only NC FHC is supported right now!");
770 
771  std::string detector = sample.detector == covmx::Detector::kNearDet ? "neardet" : "fardet";
772  std::string name = "systmaker_nc_" + detector;
773 
774  std::string filepath = "/pnfs/nova/persistent/users/jhewes15/nus19/isysts/isysts_" + detector + ".root";
775  if (RunningOnGrid()) filepath = pnfs2xrootd(filepath);
776 
777  std::cout << "Loading systs for sample " << name << ":" << std::endl;
778 
779  std::cout << "Loading systs for sample " << name.str() << ":" << std::endl;
780 
781  TFile* f = TFile::Open(filepath.c_str(), "read");
782  TDirectory* dir = f->GetDirectory(name.c_str());
783 
784  std::vector<const ISyst*> systs;
785 
786  TIter next(f->GetListOfKeys());
787  TKey* key;
788  while ((key = (TKey*)next())) {
789  systs.push_back(NusISyst::LoadFrom((TDirectory*)key->ReadObj()).release());
790  }
791 
792  for (auto syst : systs)
793  std::cout << " " << syst->ShortName() << std::endl;
794  std::cout << std::endl;
795 
796  return systs;
797 
798 } // function GetSystsFromFile
799 
800 //===========================================================================================
801 /// Default legend
802 TLegend* MakeLegendNus18(double x1, double y1, double x2, double y2, double textsize = 0.03)
803 {
804  TLegend* leg = new TLegend(x1, y1, x2, y2);
805  leg->SetBorderSize(0);
806  leg->SetFillColor(kWhite);
807  leg->SetFillStyle(0);
808  leg->SetFillStyle(0);
809  leg->SetTextFont(42);
810  leg->SetTextSize(textsize);
811  return leg;
812 }
813 
814 //===========================================================================================
815 /// Add text to plot with delta m^2 41 at 0.5 eV^2
816 void PlotTextNus18(const double xpos, const double start, const double step)
817 {
818  TLatex* tex = new TLatex();
819  tex->SetNDC();
820  tex->SetTextFont(42);
821  tex->SetTextSize(0.037);
822  tex->SetTextAlign(11);
823  tex->DrawLatex(xpos, start, "Neutrino Beam:");
824  tex->DrawLatex(xpos, start - 1*step, "8.0#times10^{20} ND POT,");
825  tex->DrawLatex(xpos, start - 2*step, "8.85#times10^{20} FD POT-equiv.");
826  tex->DrawLatex(xpos, start - 3*step, "#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
827  tex->DrawLatex(xpos, start - 4*step, "#theta_{13} = 8.3#circ, ^{}#theta_{23} = 48.4#circ");
828  tex->DrawLatex(xpos, start - 5*step, "#Deltam^{2}_{41} = 0.5 eV^{2}");
829 }
830 
831 //===========================================================================================
832 /// Add text to plot with delta m^2 41 at custom value
833 void PlotTextNus18Dm41(const double xpos, const double start, const double step, std::string dm)
834 {
835  TLatex* tex = new TLatex();
836  tex->SetNDC();
837  tex->SetTextFont(42);
838  tex->SetTextSize(0.037);
839  tex->SetTextAlign(11);
840  tex->DrawLatex(xpos, start, "Neutrino Beam:");
841  tex->DrawLatex(xpos, start - 1*step, "8.0#times10^{20} ND POT,");
842  tex->DrawLatex(xpos, start - 2*step, "8.85#times10^{20} FD POT-equiv.");
843  tex->DrawLatex(xpos, start - 3*step, "#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
844  tex->DrawLatex(xpos, start - 4*step, "#theta_{13} = 8.3#circ, ^{}#theta_{23} = 48.4#circ");
845  tex->DrawLatex(xpos, start - 5*step, Form("#Deltam^{2}_{41} = %s eV^{2}",dm.c_str()));
846 }
847 
848 //===========================================================================================
849 /// Add text to plot without delta m^2 41 value
850 void PlotTextNus18WithoutDm41(const double xpos, const double start, const double step)
851 {
852  TLatex* tex = new TLatex();
853  tex->SetNDC();
854  tex->SetTextFont(42);
855  tex->SetTextSize(0.037);
856  tex->SetTextAlign(11);
857  tex->DrawLatex(xpos, start, "Neutrino beam:");
858  tex->DrawLatex(xpos, start - 1*step, "8.0#times10^{20} ND POT,");
859  tex->DrawLatex(xpos, start - 2*step, "8.85#times10^{20} FD POT-equiv.");
860  tex->DrawLatex(xpos, start - 3*step, "#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
861  tex->DrawLatex(xpos, start - 4*step, "#theta_{13} = 8.3#circ, ^{}#theta_{23} = 48.4#circ");
862 }
863 
864 //===========================================================================================
865 /// Gaussian constrains on atmospheric parameters
866 std::vector<const IExperiment*> GetNus18Constraints()
867 {
868  std::vector<const IExperiment*> ret;
870  ret.push_back(new const GaussianConstraint(&kFitDmSq32Sterile, kNus18Dm32, kNus18Dm32Sigma));
871  return ret;
872 }
873 
874 //===========================================================================================
875 /// Seed values for Nus18 marginalisation parameters
877 {
878  // Sterile parameters
879  kFitDmSq41Sterile.SetValue(calc, 0.01);
880  kFitSinSqTheta24Sterile.SetValue(calc, 0.01);
881  kFitSinSqTheta34Sterile.SetValue(calc, 0.01);
883 
884  // Atmospheric parameters
887 
888  // Return seed values
889  std::vector<double> ret;
890  ret.push_back(kFitTheta23Sterile.GetValue(calc));
891  ret.push_back(kFitDmSq32Sterile.GetValue(calc));
892  ret.push_back(kFitDelta24InPiUnitsSterile.GetValue(calc));
893 
894  return ret;
895 }
896 
897 //===========================================================================================
898 std::vector<const IFitVar*> GetNus18FitVars()
899 {
900  std::vector<const IFitVar*> ret;
901  ret.push_back(&kFitTheta23Sterile);
902  ret.push_back(&kFitDmSq32Sterile);
903  ret.push_back(&kFitDelta24InPiUnitsSterile);
904  return ret;
905 }
906 
907 //===========================================================================================
908 std::vector<double> GetEdges(TH1* h)
909 {
910  int nbins = h->GetNbinsX();
911  std::vector<double> edges(nbins+1);
912  for (int i = 0; i <= nbins; ++i)
913  edges[i] = h->GetXaxis()->GetBinLowEdge(i+1);
914  return edges;
915 }
916 
917 //===========================================================================================
918 TH1D* DiagonalErrors(TH2D* m)
919 {
920  int nbins = m->GetNbinsX();
921  std::vector<double> edges = GetEdges((TH1*)m);
922  TH1D* h = new TH1D(UniqueName().c_str(),"",nbins,&edges[0]);
923  for (int i = 1; i <= nbins; ++i) {
924  h->SetBinContent(i,1);
925  double var = m->GetBinContent(i,i);
926  if (var > 0) h->SetBinError(i,sqrt(var));
927  else h->SetBinError(i,0);
928  }
929  return h;
930 }
931 
932 //===========================================================================================
933 std::vector<std::tuple<ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t>> GetComponents()
934 {
935  std::vector<std::tuple<ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t> > components;
936  std::vector<ana::Flavors::Flavors_t> flavors = {ana::Flavors::kNuEToNuE, ana::Flavors::kNuEToNuMu,
939  std::vector<ana::Sign::Sign_t> signs = {ana::Sign::kNu, ana::Sign::kAntiNu};
940  for (auto & flavor: flavors) {
941  for (auto & sign: signs) {
942  components.push_back( std::make_tuple(flavor, ana::Current::kCC, sign) );
943  }
944  }
945  components.push_back( std::make_tuple(ana::Flavors::kAll, ana::Current::kNC, ana::Sign::kBoth));
946  return components;
947 }
948 
949 //------------------------------------------------------------------------------
950 TLatex* MiscText(double x, double y, double size, TString text)
951 {
952  TLatex *l = new TLatex(x, y, text);
953  l->SetNDC();
954  l->SetTextSize(size);
955  l->SetTextColor(kBlack);
956  l->Draw();
957  return l;
958 }
959 
960 //===========================================================================================
962 {
963  if (i==0) return std::string("#nu_{e} #rightarrow #nu_{e}");
964  else if (i==1) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{e}");
965  else if (i==2) return std::string("#nu_{e} #rightarrow #nu_{#mu}");
966  else if (i==3) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{#mu}");
967  else if (i==4) return std::string("#nu_{e} #rightarrow #nu_{#tau}");
968  else if (i==5) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{#tau}");
969  else if (i==6) return std::string("#nu_{#mu} #rightarrow #nu_{e}");
970  else if (i==7) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{e}");
971  else if (i==8) return std::string("#nu_{#mu} #rightarrow #nu_{#mu}");
972  else if (i==9) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{#mu}");
973  else if (i==10) return std::string("#nu_{#mu} #rightarrow #nu_{#tau}");
974  else if (i==11) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{#tau}");
975  else if (i==12) return std::string("NC");
976  std::cerr << "What are you doing???" << std::endl;
977  return "no";
978 }
979 
980 //===========================================================================================
981 void FormatFullCovMxPlot(TH2D* h, PredictionConcat* sim)
982 {
983  double last_edge = sim->GetBinningConcatenated().Edges().back();
984  double max = 13. * last_edge;
985  h->Draw("colz");
986 
987  // Draw grey lines
988  for (int i = 0; i < 13; ++i) {
989  double pos = ((double)i + 0.5) * last_edge;
990  TLine* lx = new TLine(0, pos, max, pos);
991  lx->SetLineColor(kGray+1);
992  lx->Draw();
993  TLine* ly = new TLine(pos, 0, pos, max);
994  ly->SetLineColor(kGray+1);
995  ly->Draw();
996  }
997 
998  // Draw black lines
999  for (int i = 0; i <= 13; ++i) {
1000  double pos = i * last_edge;
1001  TLine* lx = new TLine(0, pos, max, pos);
1002  lx->Draw();
1003  TLine* ly = new TLine(pos, 0, pos, max);
1004  ly->Draw();
1005  }
1006 
1007  // Component names
1008  for (int i = 0; i < 13; ++i) {
1009  double pos = 0.16 + (0.72*(i/13.));
1010  auto l1 = MiscText(pos, 0.88, 0.04, GetNus18ComponentName(i).c_str());
1011  l1->SetTextAngle(38);
1012  auto l2 = MiscText(0.12, pos, 0.04, GetNus18ComponentName(i).c_str());
1013  l2->SetTextAlign(31);
1014  }
1015 
1016  h->SetNdivisions(0, "X");
1017  h->SetNdivisions(0, "Y");
1018 }
1019 
1020 //===========================================================================================
1021 bool CheckOption(TString opts, TString opt)
1022 {
1023  if (opts.Contains(opt, TString::ECaseCompare::kIgnoreCase))
1024  return true;
1025  else return false;
1026 }
1027 
1028 //===========================================================================================
1029 std::vector<bool> GetOptionConfig(TString opt, std::vector<std::string> options)
1030 {
1031  size_t n_opts = options.size();
1032 
1033  std::vector<bool> vec(n_opts);
1034  for (size_t i = 0; i < n_opts; ++i)
1035  if (CheckOption(opt, options[i]))
1036  vec[i] = true;
1037 
1038  if (std::none_of(vec.begin(), vec.end(), [](bool v) { return v; })) {
1039  std::string err = "Option string must contain ";
1040  for (size_t i = 0; i < n_opts-2; ++i)
1041  err += "\"" + options[i] + "\", ";
1042  err += options[n_opts-2] + " or " + options[n_opts-1] + ".";
1043  throw std::runtime_error(err);
1044  }
1045 
1046  else {
1047  std::cout << "Running for";
1048  for (size_t i = 0; i < n_opts; ++i)
1049  if (vec[i])
1050  std::cout << " " << options[i];
1051  std::cout << std::endl;
1052  }
1053 
1054  return vec;
1055 
1056 } // function GetOptionConfig
1057 
1058 //===========================================================================================
1060 {
1061  size_t val = 999;
1062 
1063  std::string optstring = std::string(opt.Data());
1064  size_t pos = optstring.find(key);
1065  if (pos != std::string::npos) {
1066  if (optstring.substr(pos+key.size(), 1) != "=") {
1067  std::ostringstream err;
1068  err << "The " << key << " option must take the form \"" << key << "=xxx\".";
1069  throw std::runtime_error(err.str());
1070  }
1071  size_t start = pos + key.size() + 1;
1072  size_t end = optstring.find("_", pos);
1073  if (end != std::string::npos)
1074  val = std::stod(optstring.substr(start, end - start));
1075  else
1076  val = std::stod(optstring.substr(start));
1077 
1078  std::cout << "Using " << key << " value " << val << std::endl;
1079  return val;
1080  }
1081  std::ostringstream err;
1082  err << "No value for " << key << " found in option string.";
1083  throw std::runtime_error(err.str());
1084 
1085 } // function GetOptionValue
1086 
1087 //===========================================================================================
1088 // Take a concatenated fake data spectrum and split it up
1089 std::vector<TH1D*> SplitFakeData(TH1D* h, std::vector<Binning> b)
1090 {
1091  // Define split-up histograms
1092 
1093  std::vector<TH1D*> h_split;
1094  size_t n_bins = 0;
1095  for (size_t i = 0; i < b.size(); ++i) {
1096  h_split.push_back(MakeTH1D(UniqueName().c_str(), "", b[i]));
1097  n_bins += b[i].NBins();
1098  }
1099 
1100  // Check binning is consistent
1101 
1102  if (h->GetNbinsX() != int(n_bins))
1103  throw std::runtime_error("There are " + std::to_string(h->GetNbinsX())
1104  + " bins in merged histogram but a total of " + std::to_string(n_bins)
1105  + " bins across the binning vector provided.");
1106 
1107  int bin_count = 0;
1108 
1109  for (size_t i = 0; i < b.size(); ++i) {
1110  for (int j = 1; j <= h_split[i]->GetNbinsX(); ++j) {
1111  double val = h->GetBinContent(bin_count+j);
1112  h_split[i]->SetBinContent(j, val);
1113  h_split[i]->SetBinError(j, sqrt(val));
1114  // std::cout << "Full hist bin " << bin_count+j
1115  // << " with value " << val << " is being assigned to bin "
1116  // << j << " of split histogram " << i+1 << std::endl;
1117  }
1118  bin_count += h_split[i]->GetNbinsX();
1119  }
1120 
1121  return h_split;
1122 
1123 } // function SplitFakeData
1124 
1125 //===========================================================================================
1126 void DrawSurfacePoint(PredictionConcat* sim,
1127  std::vector<std::pair<std::string, osc::IOscCalcAdjustable*>> calcs,
1128  TH1D* cosmic_merged,
1129  TH1D* data_merged,
1130  std::vector<double> pot,
1131  std::vector<bool> is_nd,
1132  TDirectory* dir,
1133  std::string name)
1134 {
1135  DontAddDirectory guard;
1136 
1137  std::vector<int> colours = { kAzure-2, kRed+2, kGreen+3, kMagenta+2 };
1138  int cosmic_colour = kGray+2;
1139 
1140  // Get the number of detectors and universes
1141 
1142  if (calcs.empty()) throw std::runtime_error("Spectrum vector is empty.");
1143  size_t n_universes = calcs.size();
1144  if (n_universes > colours.size())
1145  throw std::runtime_error("Only " + std::to_string(colours.size())
1146  + " universes per plot are supported!");
1147  size_t n_detectors = sim->GetNSamples();
1148  if (pot.size() != n_detectors)
1149  throw std::runtime_error(std::to_string(n_detectors) + " detectors but only "
1150  + std::to_string(pot.size()) + " elements in POT vector.");
1151  if (is_nd.size() != n_detectors)
1152  throw std::runtime_error(std::to_string(n_detectors) + " detectors but only "
1153  + std::to_string(is_nd.size()) + " elements in is_nd vector.");
1154 
1155  // Use OscCalcs and simulation to get spectra
1156 
1157  std::vector<std::vector<Spectrum>> specs;
1158  for (size_t i = 0; i < n_universes; ++i)
1159  specs.push_back(sim->GetSpectra(calcs[i].second));
1160 
1161  // Split up data histogram
1162 
1163  std::vector<TH1D*> data = SplitFakeData(data_merged, sim->GetBinning());
1164  std::vector<TH1D*> cosmic = SplitFakeData(cosmic_merged, sim->GetBinning());
1165 
1166  // Define shape of canvas
1167 
1168  TCanvas* c = new TCanvas(UniqueName().c_str(), "c", 500 * n_detectors, 800);
1169  c->Divide(n_detectors, 3);
1170 
1171  // Pad margins
1172 
1173  // std::vector<float> lm = { 0.15, 0.12 };
1174  // std::vector<float> rm = { 0.05, 0.08 };
1175  std::vector<float> margin = { 0.1, 0.05 };
1176 
1177  gStyle->SetPadLeftMargin(margin[0]);
1178  gStyle->SetPadRightMargin(margin[1]);
1179 
1180  gStyle->SetLabelFont(43, "xyz");
1181  gStyle->SetLabelSize(14, "xyz");
1182 
1183  for (size_t d = 0; d < n_detectors; ++d)
1184  {
1185  // Get pad x ranges for this detector
1186 
1187  float xmin = float(d) / float(n_detectors);
1188  float xmax = float(d+1) / float(n_detectors);
1189 
1190  // Set size of spectrum pad
1191 
1192  c->cd(d+1);
1193  gPad->SetPad(xmin, 0.45, xmax, 1);
1194 
1195  // Set size of ratio pad
1196 
1197  c->cd(n_detectors+d+1);
1198  gPad->SetPad(xmin, 0.25, xmax, 0.45);
1199 
1200  // Set size of oscillation pad
1201 
1202  c->cd((2*n_detectors)+d+1);
1203  gPad->SetPad(xmin, 0.05, xmax, 0.25);
1204 
1205  } // for detector i
1206 
1207  TLegend* l = new TLegend(0.5, 0.65, 0.85, 0.85);
1208  l->SetFillStyle(0);
1209 
1210  // Create histogram vectors and set their sizes
1211 
1212  std::vector<std::vector<TH1D*>> hists, ratio_hists;
1213 
1214  hists.resize(n_universes);
1215  for (std::vector<TH1D*>& h : hists)
1216  h.resize(n_detectors);
1217 
1218  ratio_hists.resize(n_universes);
1219  for (std::vector<TH1D*>& h : ratio_hists)
1220  h.resize(n_detectors);
1221 
1222  std::vector<TGraph*> graphs;
1223 
1224  // Y axis max for each detector's plot
1225 
1226  std::vector<double> ymax, rmax;
1227  std::vector<bool> rescaled;
1228 
1229  // Make data ratio histograms
1230 
1231  std::vector<TH1D*> data_ratio;
1232  data_ratio.resize(n_detectors);
1233  for (size_t d = 0; d < n_detectors; ++d)
1234  {
1235  data_ratio[d] = (TH1D*)data[d]->Clone();
1236  for (int b = 1; b <= data_ratio[d]->GetNbinsX(); ++b)
1237  {
1238  double val = data[d]->GetBinContent(b);
1239  double err = data[d]->GetBinError(b);
1240  data_ratio[d]->SetBinContent(b, 1);
1241 
1242  if (val == 0)
1243  data_ratio[d]->SetBinError(b, 0);
1244  else
1245  data_ratio[d]->SetBinError(b, err / val);
1246 
1247  } // for bin b
1248  } // for detector d
1249 
1250  // Loop over each detector
1251 
1252  for (size_t d = 0; d < n_detectors; ++d)
1253  {
1254  // Loop over each universe, get histograms and get ymax
1255 
1256  std::vector<double> m, r;
1257 
1258  // Rescale data histogram
1259 
1260  if (is_nd[d]) data[d]->Scale(1e-4);
1261 
1262  for (size_t u = 0; u < n_universes; ++u)
1263  {
1264  hists[u][d] = (specs[u][d].ToTH1(pot[d]));
1265 
1266  // If it's a near detector spectrum, rescale it
1267 
1268  if (is_nd[d]) hists[u][d]->Scale(1e-4);
1269  m.push_back(hists[u][d]->GetMaximum());
1270 
1271  // Clone to create ratio histogram
1272 
1273  ratio_hists[u][d] = (TH1D*)hists[u][d]->Clone();
1274  TH1D* h = ratio_hists[u][d];
1275  for (int b = 1; b <= h->GetNbinsX(); ++b)
1276  {
1277  double val = (hists[u][d]->GetBinContent(b) + cosmic[d]->GetBinContent(b)) / data[d]->GetBinContent(b);
1278  if (std::isinf(val) || std::isnan(val)) {
1279  h->SetBinContent(b, 1);
1280  h->SetBinError(b, 0);
1281  }
1282  else {
1283  h->SetBinContent(b, val);
1284  r.push_back(fabs(1-val));
1285  }
1286 
1287  } // for bin b
1288 
1289  } // for universe u
1290 
1291  // Get maximum data point
1292 
1293  int maxdata = data[d]->GetMaximumBin();
1294  m.push_back(data[d]->GetBinContent(maxdata) + data[d]->GetBinError(maxdata));
1295 
1296  // Account for data errors in ratio plot range
1297 
1298  for (int b = 1; b <= data_ratio[d]->GetNbinsX(); ++b)
1299  r.push_back(fabs(data_ratio[d]->GetBinError(b)));
1300 
1301  ymax.push_back(*std::max_element(m.begin(), m.end()));
1302  rmax.push_back(*std::max_element(r.begin(), r.end()));
1303 
1304  } // for detector d
1305 
1306  // Formatting options
1307 
1308  float ts = 0.055;
1309 
1310  // Loop over each detector to draw the cosmics
1311 
1312  for (size_t d = 0; d < n_detectors; ++d) {
1313 
1314  // Draw the cosmic spectrum
1315 
1316  c->cd(d+1);
1317  TH1D* hc = cosmic[d];
1318  hc->SetLineColor(cosmic_colour);
1319  hc->SetLineWidth(2);
1320  hc->SetMinimum(0);
1321  hc->SetMaximum(1.1 * ymax[d]);
1322  hc->Draw("hist same");
1323 
1324  }
1325 
1326  // Loop over each universe
1327 
1328  for (size_t u = 0; u < n_universes; ++u)
1329  {
1330  // Loop over each detector
1331 
1332  for (size_t d = 0; d < n_detectors; ++d)
1333  {
1334  // Format and draw the spectrum
1335 
1336  c->cd(d+1);
1337  TH1D* h = hists[u][d];
1338  h->Add(cosmic[d]);
1339  h->SetLineColor(colours[u]);
1340  h->SetLineWidth(2);
1341  h->SetMinimum(0);
1342  h->SetMaximum(1.1 * ymax[d]);
1343 
1344  if (is_nd[d])
1345  h->SetTitle(";Energy deposited in detector (GeV);10^{4} events");
1346  else
1347  h->SetTitle(";Energy deposited in detector (GeV);Events");
1348  h->GetXaxis()->SetTitleSize(ts);
1349  h->GetYaxis()->SetTitleSize(ts);
1350 
1351  h->Draw("hist same");
1352 
1353  // Now draw the ratio plot
1354 
1355  h = ratio_hists[u][d];
1356  c->cd(n_detectors+d+1);
1357  h->SetLineColor(colours[u]);
1358  h->SetLineWidth(2);
1359  h->SetMinimum(1 - (1.1 * rmax[d]));
1360  h->SetMaximum(1 + (1.1 * rmax[d]));
1361 
1362  h->SetTitle(";;Ratio");
1363  h->GetYaxis()->SetTitleSize(3*ts);
1364 
1365  h->Draw("hist same");
1366 
1367  // Now the oscillation plot
1368 
1369  c->cd((2*n_detectors)+d+1);
1370 
1371  std::vector<float> x, y;
1372  double l = calcs[u].second->GetL();
1373  if (is_nd[d]) calcs[u].second->SetL(1);
1374  else calcs[u].second->SetL(810);
1375  for (size_t k = 0; k < 200; ++k)
1376  {
1377  x.push_back(0.1 * (float(k) + 0.5));
1378  y.push_back(calcs[u].second->P(14, 0, x.back()));
1379  }
1380  calcs[u].second->SetL(l);
1381 
1382  TGraph* g = new TGraph(x.size(), &x[0], &y[0]);
1383 
1384  g->SetLineWidth(2);
1385  g->SetLineColor(colours[u]);
1386 
1387  g->SetTitle(";;P(#nu_{#mu} #rightarrow #nu_{s})");
1388  g->GetYaxis()->SetTitleSize(3*ts);
1389 
1390  g->Draw("ac");
1391  graphs.push_back(g);
1392 
1393  } // for detector d
1394 
1395  l->AddEntry(hists[u][0], calcs[u].first.c_str(), "l");
1396 
1397  } // for universe u
1398 
1399  l->AddEntry(cosmic[0], "Cosmic background", "l");
1400 
1401  if (data.size() != n_detectors)
1402  throw std::runtime_error("Number of data hists does not match number of samples!");
1403 
1404  // Loop and draw the data for each sample
1405 
1406  for (size_t d = 0; d < n_detectors; ++d)
1407  {
1408  c->cd(d+1);
1409  data[d]->SetMarkerSize(2);
1410  data[d]->Draw("e0 same");
1411 
1412  c->cd(n_detectors+d+1);
1413  data_ratio[d]->SetMarkerSize(2);
1414  data_ratio[d]->Draw("e0 same");
1415 
1416  } // for detector d
1417 
1418  l->AddEntry(data[0], "Data", "lep");
1419  c->cd(1);
1420  l->Draw();
1421 
1422  dir->WriteTObject(c, name.c_str());
1423 
1424  // std::string img = name + ".png";
1425  // c->SaveAs(img.c_str());
1426 
1427  // Now clean up
1428 
1429  delete c;
1430  delete l;
1431  for (auto g : graphs)
1432  delete g;
1433  for (auto h_vec : hists)
1434  for (auto h : h_vec)
1435  delete h;
1436  for (auto h_vec : ratio_hists)
1437  for (auto h : h_vec)
1438  delete h;
1439  for (auto h : data_ratio)
1440  delete h;
1441 
1442 } // function DrawSurfacePoint
1443 
const XML_Char * name
Definition: expat.h:151
std::vector< double > GetNus18SeedValues(osc::IOscCalcAdjustable *calc)
Seed values for Nus18 marginalisation parameters.
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
double GetDelta(int i, int j) const
IPrediction * LoadExtrapPrediction(TString opt)
Spectrum LoadCosmicSpectrum(bool rhc=false)
Function to return cosmic spectrum.
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
double stod(const std::string &val)
std::vector< std::tuple< ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t > > GetComponents()
void FormatFullCovMxPlot(TH2D *h, PredictionConcat *sim)
bool RunningOnGrid()
Is this a grid (condor) job?
Definition: Utilities.cxx:587
dictionary components
const std::string kFDMCRHCProd4
A simple Gaussian constraint on an arbitrary IFitVar.
const double kNus18Delta13
std::vector< TH1D * > SplitFakeData(TH1D *h, std::vector< Binning > b)
Float_t x1[n_points_granero]
Definition: compare.C:5
(&#39;beam &#39;)
Definition: IPrediction.h:15
void PlotTextNus18(const double xpos, const double start, const double step)
Add text to plot with delta m^2 41 at 0.5 eV^2.
TH1D * MakeTH1D(const char *name, const char *title, const Binning &bins)
Utility function to avoid need to switch on bins.IsSimple()
Definition: Utilities.cxx:215
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double kNus18Dm21
void SetDelta(int i, int j, double delta)
std::vector< const ISyst * > LoadISysts()
Function to load ISysts.
Binning RHCNDBins()
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
OStream cerr
Definition: OStream.cxx:7
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
std::string GetNus18ComponentName(int i)
string filename
Definition: shutoffs.py:106
Float_t tmp
Definition: plot.C:36
const int kOfficialNus18FHCColor
TH1D * LoadCosmicHist(bool rhc=false)
Function to load cosmic histogram.
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: Utilities.cxx:620
void PrintOscParams(osc::OscCalcSterile *calc)
const std::string kFDDataFHC
TString hists[nhists]
Definition: bdt_com.C:3
std::vector< const IExperiment * > GetNus18Constraints()
Gaussian constrains on atmospheric parameters.
void DrawSurfacePoint(PredictionConcat *sim, std::vector< std::pair< std::string, osc::IOscCalcAdjustable * >> calcs, TH1D *cosmic_merged, TH1D *data_merged, std::vector< double > pot, std::vector< bool > is_nd, TDirectory *dir, std::string name)
const double kNus18Th23Sigma
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
const double kNus18Th13
osc::OscCalcDumb calc
void FormatCanvas()
Canvas formatting utility.
std::string FindCAFAnaDir()
Definition: Utilities.cxx:429
void SetAngles(osc::OscCalcSterile *calc)
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
#define M_PI
Definition: SbMath.h:34
Binning RHCFDBins()
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const XML_Char const XML_Char * data
Definition: expat.h:268
Double_t ymax
Definition: plot.C:25
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
const int kOfficialNus18RHCColor
const XML_Char * s
Definition: expat.h:262
CovarianceMatrix * LoadCovMx(std::string file, std::string name)
Function to load covariance matrix.
const int nbins
Definition: cellShifts.C:15
Charged-current interactions.
Definition: IPrediction.h:39
const Binning kNus18EnergyBinning
Energy binnings for Nus18 for nus analysis.
Definition: Binning.cxx:321
Sum up livetimes from individual cosmic triggers.
TLegend * MakeLegendNus18(double x1, double y1, double x2, double y2, double textsize=0.03)
Default legend.
const std::string kFDDataRHC
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:607
Binning FHCNDBins()
Binning FHCFDBins()
std::string getenv(std::string const &name)
Spectrum LoadData(bool nd, bool rhc=false)
Function to load NuS data file.
std::vector< std::string > flavors
std::vector< bool > GetOptionConfig(TString opt, std::vector< std::string > options)
IPrediction * LoadPPFX(std::string detector, int universe)
Function to load PPFX universe simulations.
Detector detector
Definition: Sample.h:101
#define pot
const double kNus18Dm32
void AddBin(std::vector< double > &v, double b)
Float_t d
Definition: plot.C:236
const double kNus18Th12
const double kNus18Th23
const std::string kFDCosmicRHCProd4
const double j
Definition: BetheBloch.cxx:29
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Eigen::VectorXd vec
string pname
Definition: eplot.py:33
static Binning Custom(const std::vector< double > &edges)
Definition: Binning.cxx:83
double GetValue(const osc::IOscCalcAdjustable *osc) const override
std::vector< float > Spectrum
Definition: Constants.h:527
void PlotTextNus18WithoutDm41(const double xpos, const double start, const double step)
Add text to plot without delta m^2 41 value.
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:15
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:76
std::string GetPPFXWeightsDir()
Function to get location of PPFX weights.
bool CheckOption(TString opts, TString opt)
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
std::vector< Binning > GetNus18Binning(bool rhc, std::string detector)
(&#39; survival&#39;)
Definition: IPrediction.h:19
const std::string path
Definition: plot_BEN.C:43
TH1D * DiagonalErrors(TH2D *m)
void SetAngle(int i, int j, double th)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const std::string kFDCosmicFHCProd4
void SetDm(int i, double dm)
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
const std::string kFDMCFHCProd4
TDirectory * dir
Definition: macro.C:5
void PlotTextNus18Dm41(const double xpos, const double start, const double step, std::string dm)
Add text to plot with delta m^2 41 at custom value.
TLatex * MiscText(double x, double y, double size, TString text)
double livetime
Definition: saveFDMCHists.C:21
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
TLatex * tex
Definition: f2_nu.C:499
const hit & b
Definition: hits.cxx:21
Neutral-current interactions.
Definition: IPrediction.h:40
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
TFile * file
Definition: cellShifts.C:17
std::vector< TH1D * > GetCosmics(bool rhc, std::string detector)
Function to get default cosmics.
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
double GetDm(int i) const
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
::xsd::cxx::tree::qname< char, simple_type, uri, ncname > qname
Definition: Database.h:175
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::vector< const IFitVar * > GetNus18FitVars()
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
Prevent histograms being added to the current directory.
Definition: Utilities.h:62
Float_t e
Definition: plot.C:35
Spectrum LoadFakeData(std::string filename, int universe, double pot, double livetime)
Function to load the unblinded data spectrum.
Polarity polarity
Definition: Sample.h:100
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
std::vector< double > GetEdges(TH1 *h)
void next()
Definition: show_event.C:84
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
const double kNus18Dm32Sigma
const double kNus18Delta24
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
Definition: Exposures.h:213
def sign(x)
Definition: canMan.py:197
double GetAngle(int i, int j) const
const double kAna2018RHCLivetime
Definition: Exposures.h:214
size_t GetOptionValue(TString opt, std::string key)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void CombineMatrices(CovarianceMatrix *gen, std::map< const ISyst *, CovarianceMatrix * > m, std::vector< IPrediction * > preds, std::vector< double > pot, osc::IOscCalcAdjustable *calc)
Function to combine covariance matrices.
std::vector< const ISyst * > GetSystsFromFile(covmx::Sample sample)
Get systematics for a given sample.
const FitTheta23Sterile kFitTheta23Sterile