Nus18Utilities.h

Go to the documentation of this file.

#pragma once

#include <vector>
#include <tuple>

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "NuXAna/Analysis/NusSystsMaker.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Sample.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Prediction/CovarianceMatrix.h"
#include "CAFAna/Prediction/PredictionConcat.h"
#include "NuXAna/Systs/Nus18Systs.h"
#include "CAFAna/Vars/FitVarsSterile.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "Utilities/rootlogon.C"

#include "TLegend.h"
#include "TLatex.h"
#include "TLine.h"
#include "TGraph.h"
#include "THStack.h"

using namespace ana;

// Hard-coded parameters
const double kNus18Dm21 = 7.53e-5;
const double kNus18Dm32 = 2.44e-3;
const double kNus18Dm32Sigma = 0.08e-3;

const double kNus18Th12 = 0.587;
const double kNus18Th13 = 0.145;
const double kNus18Th23 = 0.844;
const double kNus18Th23Sigma = 0.042;

const double kNus18Delta13 = 1.21*M_PI;
const double kNus18Delta24 = M_PI/2;

const int kOfficialNus18FHCColor = kAzure+2;
const int kOfficialNus18RHCColor = kRed+2;

// Hard-coded loaders
// cosmics
const std::string kFDCosmicFHCProd4 = "prod_caf_R17-11-14-prod4reco.b_fd_cosmic_fhcTune_full_v1_goodruns";
const std::string kFDCosmicRHCProd4 = "prod_caf_R17-11-14-prod4reco.a_fd_cosmic_rhcTune_HighGain_v1_goodruns_snapshot_170116";

// data
const std::string kFDDataFHC = "prod4_sumrestricteddecaf_fd_numi_fhc_full_goodruns_nus2018";
const std::string kFDDataRHC = "prod4_sumrestricteddecaf_fd_numi_rhc_full_goodruns_nus2018";

// mc (in general use loaders, just for specific purposes)
const std::string kFDMCFHCProd4 = "prod_caf_R17-11-14-prod4reco.neutron-respin.c_fd_genie_nonswap_fhc_nova_v08_full_v1";
const std::string kFDMCRHCProd4 = "prod_caf_R17-11-14-prod4reco.neutron-respin.c_fd_genie_nonswap_rhc_nova_v08_full_v1";

//===========================================================================================
void AddBin(std::vector<double>& v, double b)
{ v.push_back(v.back()+b); }

//===========================================================================================
Binning FHCNDBins()
{
  std::vector<double> v = { 0. };
  for (size_t i = 0; i < 6; ++i)
    AddBin(v, 1);
  for (size_t i = 0; i < 2; ++i)
    AddBin(v, 2);
  AddBin(v, 10);

  // AddBin(v,0.75);
  // for (int i = 0; i < 37; ++i)
  //   AddBin(v,0.25);
  // for (int i = 0; i < 8; ++i)
  //   AddBin(v,0.5);
  // for (int i = 0; i < 6; ++i)
  //   AddBin(v,1);

  return Binning::Custom(v);

  /* AddBin(v, 0.75);
  for (int i = 0; i < 37; ++i) AddBin(v, 0.25);
  for (int i = 0; i < 2; ++i)  AddBin(v, 0.5);
  for (int i = 0; i < 2; ++i)  AddBin(v, 1);
  AddBin(v, 2);
  AddBin(v, 5);
  */
}

//===========================================================================================
Binning FHCFDBins()
{
  std::vector<double> v = { 0. };
  AddBin(v,1.25);
  AddBin(v,0.25);
  AddBin(v,0.25);
  AddBin(v,0.5);
  AddBin(v,1);
  AddBin(v,1.5);
  AddBin(v,2);
  AddBin(v,3.25);
  AddBin(v,10);

  return Binning::Custom(v);
}

//===========================================================================================
Binning RHCFDBins()
{
  std::vector<double> v = { 0. };
  AddBin(v,2.25);
  AddBin(v,3.25);
  AddBin(v,14.5);

  return Binning::Custom(v);
}

//===========================================================================================
Binning RHCNDBins()
{
  std::vector<double> v = { 0. };
  AddBin(v,0.75);
  for (int i = 0; i < 67; ++i)
    AddBin(v,0.25);
  AddBin(v,0.5);
  AddBin(v,2.0);

  return Binning::Custom(v);
}

//===========================================================================================
std::vector<Binning> GetNus18Binning(bool rhc, std::string detector)
{

  // Check detector

  bool use_nd = false;
  bool use_fd = false;
  if (detector == "nd") use_nd = true;
  else if (detector == "fd") use_fd = true;
  else if (detector == "both") {
    use_nd = true;
    use_fd = true;
  }
  else throw std::runtime_error("Detector option \"" + detector + "\" not understood.");

  // Now actually add bins

  std::vector<Binning> ret;

  if (rhc) {
    if (use_nd) ret.push_back(RHCFDBins());
    if (use_fd) ret.push_back(RHCFDBins());
  } else {
    if (use_nd) ret.push_back(FHCNDBins());
    if (use_fd) ret.push_back(FHCFDBins());
  }

  return ret;

} // function GetNus18Binning

//===========================================================================================
void SetAngles(osc::OscCalcSterile* calc)
{
  // Default vals
  ResetSterileCalcToDefault(calc);

  // Mass splittings
  calc->SetDm(2, kNus18Dm21);
  calc->SetDm(3, kNus18Dm21 + kNus18Dm32);
  calc->SetDm(4, 0.5);

  // Angles
  calc->SetAngle(1, 3, kNus18Th13);
  calc->SetAngle(2, 3, kNus18Th23);
  calc->SetAngle(3, 4, 0.);
  calc->SetAngle(2, 4, 0.);
  calc->SetAngle(1, 4, 0.);

  // Phases
  calc->SetDelta(1, 3, kNus18Delta13);
  calc->SetDelta(2, 4, kNus18Delta24);
}

//===========================================================================================
void PrintOscParams(osc::OscCalcSterile* calc)
{
  // Mass splittings
  std::cout << "---------------------------" << std::endl;
  std::cout << "  OSCILLATION PARAMETERS   " << std::endl;
  std::cout << std::endl;
  std::cout << " dmsq 21  - " << calc->GetDm(2) << " eV^2" << std::endl;
  std::cout << " dmsq 31  - " << calc->GetDm(3) << " eV^2" << std::endl;
  std::cout << " dmsq 41  - " << calc->GetDm(4) << " eV^2" << std::endl;
  std::cout << std::endl;
  std::cout << " th 12    - " << calc->GetAngle(1, 2) / M_PI * 180 << " deg" << std::endl;
  std::cout << " th 13    - " << calc->GetAngle(1, 3) / M_PI * 180 << " deg" << std::endl;
  std::cout << " th 23    - " << calc->GetAngle(2, 3) / M_PI * 180 << " deg" << std::endl;
  std::cout << " th 14    - " << calc->GetAngle(1, 4) / M_PI * 180 << " deg" << std::endl;
  std::cout << " th 24    - " << calc->GetAngle(2, 4) / M_PI * 180 << " deg" << std::endl;
  std::cout << " th 34    - " << calc->GetAngle(3, 4) / M_PI * 180 << " deg" << std::endl;
  std::cout << std::endl;
  std::cout << " delta 13 - " << calc->GetDelta(1, 3) / M_PI * 180 << " deg" << std::endl;
  std::cout << " delta 14 - " << calc->GetDelta(1, 4) / M_PI * 180 << " deg" << std::endl;
  std::cout << " delta 24 - " << calc->GetDelta(2, 4) / M_PI * 180 << " deg" << std::endl;
  std::cout << "---------------------------" << std::endl;
}

//===========================================================================================
/// Canvas formatting utility
void FormatCanvas()
{
  gStyle->SetCanvasDefW(1200);
  gStyle->SetCanvasDefH(800);
  gStyle->SetPadTopMargin(0.14);
  gStyle->SetPadBottomMargin(0.14);
  gStyle->SetPadLeftMargin(0.14);
  gStyle->SetPadRightMargin(0.14);
}

//===========================================================================================
/// Function to load prediction object
IPrediction* LoadPrediction(std::string detector, bool rhc = false,
  std::string syst_type = "all")
{
  DontAddDirectory guard;

  if (detector != "neardet" && detector != "fardet")
    throw std::runtime_error("Detector string must be either \"neardet\" or \"fardet\"!");

  std::string dname = detector == "neardet" ? "nd" : "fd";
  std::string polarity = rhc? "rhc" : "fhc";

  // log400 hack for the time being
  if (not rhc and detector == "neardet" and syst_type == "none") {

    std::string filename = FindCAFAnaDir() + "/data/log400/pred_nc_nd_log400.root";
    return LoadFromFile<PredictionInterp>(filename.c_str(),
      "pred_interp_nus_fhc_neardet_none").release();

  }

  std::string filename = "pred_nus18_" + polarity + "_" + dname
    + "_" + syst_type + ".root";

  // std::string filename = "pred_nus_" + polarity + "_" + detector
  //   + "_" + syst_type + ".root";

  char* path = std::getenv("NUSDATA_NUS18_PRED");
  if (!path) assert(false and
    "NUSDATA_NUS18_PRED environment variable not set! Did you forget to set up nusdata?");

  // std::string path = FindCAFAnaDir() + "/data/joint/pred";
  std::cout << "Loading prediction from  " << path << "/" << filename << std::endl;

  std::string filepath = std::string(path) + "/" + filename;
  return LoadFromFile<PredictionInterp>(filepath, "pred_interp").release();

} // function LoadPrediction

//===========================================================================================
IPrediction* LoadExtrapPrediction(TString opt)
{
  DontAddDirectory guard;

  bool fhc = opt.Contains("fhc", TString::kIgnoreCase);
  bool rhc = opt.Contains("rhc", TString::kIgnoreCase);
  assert(fhc xor rhc);

  bool syst = opt.Contains("syst", TString::kIgnoreCase);

  std::stringstream filename;
  filename << "pred_nus18_";
  rhc ? filename << "rhc" : filename << "fhc";
  filename << "_extrap";
  if (syst) filename << "_systs";
  filename << ".root";

  char* path = std::getenv("NUSDATA_NUS18_PRED");
  if (!path) throw std::runtime_error(
    "NUSDATA_NUS18_PRED environment variable not set! Did you forget to set up nusdata?");
  else std::cout << "Loading prediction from  " << path << "/" << filename.str() << std::endl;

  std::string filepath = std::string(path) + "/" + filename.str();

  //return LoadFromFile<IPrediction>(filepath, "pred_systs").release();
  return LoadFromFile<IPrediction>(filepath, "FDNus18ExtrapPred").release();

} // function LoadPrediction

//===========================================================================================
/// Function to load covariance matrix
CovarianceMatrix* LoadCovMx(std::string file, std::string name)
{
  DontAddDirectory guard;

  char* path = std::getenv("NUSDATA_NUS18_COVMX");
  if (!path) throw std::runtime_error(
    "NUSDATA_NUS18_COVMX environment variable not set! Did you forget to set up nusdata?");

  std::string filepath = std::string(path) + "/" + file + ".root";
  //std::string filepath = FindCAFAnaDir() + "data/joint/covmx/" + file + ".root";
  std::cout << "Loading covariance matrix " << name << " from " << filepath << std::endl;

  return LoadFromFile<CovarianceMatrix>(filepath, name).release();
}

//===========================================================================================
/// Function to load PPFX universe simulations
IPrediction* LoadPPFX(std::string detector, int universe)
{
  DontAddDirectory guard;

  char* path = std::getenv("NUSDATA_NUS18_WEIGHTS");
  if (detector != "nd" && detector != "fd") throw std::runtime_error(
    "Detector string must be \"nd\" or \"fd\" - " + detector + " not understood!");
  if (!path) throw std::runtime_error(
    "NUSDATA_NUS18_WEIGHTS environment variable not set! Did you forget to set up nusdata?");
  std::string file = std::string(path) + "/nus18_fhc_" + detector + "_ppfx.root";
  return LoadFromFile<PredictionInterp>(file, "universe_" + std::to_string(universe)).release();
}

//===========================================================================================
/// Function to load the unblinded data spectrum
/*Spectrum* LoadData(bool rhc) {

  DontAddDirectory guard;

  if (std::getenv("PROCESS")) {
    std::cout << "I see you're trying to run on the grid."
              << "Now may be a good time to figure out how to read data on the grid..." << std::endl;
    abort();
  }

  std::string path = "/nova/app/users/wallbank/test_ndosc/NuXAna/macros/Nus18/box_opening/";

  return LoadFromFile<Spectrum>
    (Form("%s/fout_nus18_box_opening_%s.root", path.c_str(), rhc? "rhc" : "fhc"), "spec_nus_vise_numi").release();

}*/

//===========================================================================================
/// Function to load fake data file
// Spectrum LoadFakeData(std::string filename, int universe)
// {
//   DontAddDirectory guard;

//   char* path = std::getenv("NUSDATA_NUS18_FAKEDATA");
//   if (!path) throw std::runtime_error(
//     "NUSDATA_NUS18_FAKEDATA environment variable not set! Did you forget to set up nusdata?");
//   else std::cout << "Loading fake data spectrum universe_" << universe << " from "
//     << path << "/" << filename << ".root" << std::endl;

//   std::string filepath = std::string(path) + "/" + filename + ".root";
//   return *(LoadFromFile<Spectrum>(filepath, "universe_" + std::to_string(universe)).release());
// }

Spectrum LoadFakeData(std::string filename, int universe, double pot, double livetime)
{
  DontAddDirectory guard;

  char* path = std::getenv("NUSDATA_NUS18_FAKEDATA");
  if (!path) throw std::runtime_error(
    "NUSDATA_NUS18_FAKEDATA environment variable not set! Did you forget to set up nusdata?");
  else std::cout << "Loading fake data spectrum universe_" << universe << " from "
    << path << "/" << filename << ".root" << std::endl;

  std::string filepath = std::string(path) + "/" + filename + ".root";
  std::string histname = "universe_" + std::to_string(universe);

  TFile* f = new TFile(filepath.c_str(), "read");
  TH1D* h = (TH1D*)f->Get(histname.c_str());

  return Spectrum(h, pot, livetime);
}

//===========================================================================================
/// Function to load NuS data file
Spectrum LoadData(bool nd, bool rhc = false)
{
  DontAddDirectory guard;

  std::string path = std::string(std::getenv("SRT_PRIVATE_CONTEXT"));
  if (rhc) path = "/nova/app/users/wallbank/test_ndosc/NuXAna/macros/Nus18_neutrino/box_opening/";

  if (rhc) return *LoadFromFile<Spectrum>(path + "/fout_nus18_box_opening_rhc.root", "spec_nus_vise_numi").release();

  if (nd) return *LoadFromFile<Spectrum>(path + "/files/fhc_NDDataMCLoadOut.root","RecoEPreselPtpData").release();
  else return *LoadFromFile<Spectrum>(path + "/files/fout_nus18_box_opening_fhc.root", "spec_nus_vise_numi").release();
}

//===========================================================================================
/// Function to load cosmic histogram
TH1D* LoadCosmicHist(bool rhc = false)
{
  DontAddDirectory guard;

  char* path = std::getenv("NUSDATA_NUS18_COSMIC");
  if (!path) throw std::runtime_error(
    "NUSDATA_NUS18_COSMIC environment variable not set! Did you forget to set up nusdata?");
  else std::cout << "Loading cosmic data from  " << path << "/nus18_cosmic_background.root" << std::endl;

  TFile f(Form("%s/nus18_cosmic_background.root", path));
  return (TH1D*)f.Get(rhc? "rhcnumi" : "fhcnumi");
}

//===========================================================================================
/// Function to load cosmic histogram
TH1D* LoadCosmicHist(TString opt)
{
  DontAddDirectory guard;

  bool fhc = opt.Contains("fhc", TString::kIgnoreCase);
  bool rhc = opt.Contains("rhc", TString::kIgnoreCase);
  assert(fhc || rhc);

  char* path = std::getenv("NUSDATA_NUS18_COSMIC");
  if (!path) throw std::runtime_error(
    "NUSDATA_NUS18_COSMIC environment variable not set! Did you forget to set up nusdata?");
  else std::cout << "Loading cosmic data from  " << path << "/nus18_cosmic_background.root" << std::endl;

  TFile f(Form("%s/nus18_cosmic_background.root", path));
  return (TH1D*)f.Get(rhc? "rhcnumi" : "fhcnumi");
}

//===========================================================================================
/// Function to return cosmic spectrum
Spectrum LoadCosmicSpectrum(bool rhc = false)
{
  DontAddDirectory guard;
  return Spectrum(LoadCosmicHist(rhc), 0, rhc? kAna2018RHCLivetime : kAna2018FHCLivetime);
}

//===========================================================================================
/// Function to return cosmic spectrum
Spectrum LoadCosmicSpectrum(TString opt)
{
  DontAddDirectory guard;

  bool fhc = opt.Contains("fhc", TString::kIgnoreCase);
  bool rhc = opt.Contains("rhc", TString::kIgnoreCase);
  assert(fhc || rhc);

  return Spectrum(LoadCosmicHist(opt), 0, rhc? kAna2018RHCLivetime : kAna2018FHCLivetime);
}

/* //=========================================================================================== */
/* /// Function to load systematics */
/* NusSystematics* LoadSystematics() { */

/*   char* path = std::getenv("NUSDATA_NUS18_SYSTS"); */
/*   if (!path) throw std::runtime_error("NUSDATA_NUS18_SYSTS environment variable not set!"); */
/*   else std::cout << "Loading systematics from " << path << "/Nus18Systs.root" << std::endl; */

/*   TFile f(Form("%s/Nus18Systs.root", path)); */
/*   NusSystematics* nus_systs = NusSystematics::LoadFrom(f.GetDirectory("Nus18Systs")); */

/*   return nus_systs; */

/* }   */

//===========================================================================================
/// Function to load ISysts
std::vector<const ISyst*> LoadISysts() {

  DontAddDirectory guard;

  char* path = std::getenv("NUSDATA_NUS18_SYSTS");
  if (!path)
    throw std::runtime_error("NUSDATA_NUS18_SYSTS environment variable not set!");
  else
    std::cout << "Loading systematics from " << path << "/Nus18ISysts.root" << std::endl;

  std::vector<const ISyst*> systs;

  TFile* f = new TFile(Form("%s/Nus18ISysts.root", path), "READ");
  TIter next(f->GetListOfKeys());
  TKey* key;
  while ((key = (TKey*)next()))
    systs.push_back(NusISyst::LoadFrom(f, key->GetName())).release();
  f->Close();
  delete f;

  return systs;

}

//===========================================================================================
/// Function to get location of PPFX weights
std::string GetPPFXWeightsDir()
{
  char* ppfx_env = std::getenv("NUSDATA_NUS18_WEIGHTS");
  if (ppfx_env == nullptr) throw std::runtime_error(
    "NUSDATA_NUS18_WEIGHTS environment variable not set! Did you forget to set up nusdata?");

  return std::string(ppfx_env);
}
/*
//===========================================================================================
/// Function to get default Nus18 simulation
PredictionConcat* GetDefaultSimulation(bool rhc, std::string detector, std::string syst_type)
{
  // Detector type
  bool use_nd = false;
  bool use_fd = false;
  if (detector == "nd") use_nd = true;
  else if (detector == "fd") use_fd = true;
  else if (detector == "both") {
    use_nd = true;
    use_fd = true;
  }
  else throw std::runtime_error("Detector string must be \"nd\", \"fd\" or \"both\".");

  // Systematic option
  if (syst_type != "none" && syst_type != "xsec" && syst_type != "nonxsec")
    throw std::runtime_error("Syst type string must be \"none\", \"xsec\" or \"nonxsec\".");

  // Get exposures
  double nd_pot = rhc ? kAna2018SensitivityRHCNDPOT : kAna2018SensitivityFHCNDPOT;
  double fd_pot = rhc ? kAna2018RHCPOT : kAna2018FHCPOT;
  double fd_livetime = rhc ? kAna2018RHCLivetime : kAna2018FHCLivetime;
  std::string polarity = rhc? "RHC" : "FHC";

  // Vectors for predictions, POT & livetimes, systematic prefixes
  std::vector<IPrediction*> predictions;
  std::vector<double> pots;
  std::vector<double> livetimes;
  std::vector<std::string> samples;
  std::vector<std::vector<const ISyst*>> systs;

  // Fill
  if (use_nd) {
    samples.push_back("Near detector");
    predictions.push_back(LoadPrediction("neardet", rhc, syst_type));
    pots.push_back(nd_pot);
    livetimes.push_back(0);
    systs.push_back(GetNus18Systs(rhc, "nd", syst_type));
    // syst_types.push_back(polarity + "_ND");
  }
  if (use_fd) {
    samples.push_back("Far detector");
    predictions.push_back(LoadPrediction("fardet", rhc, syst_type));
    pots.push_back(fd_pot);
    livetimes.push_back(fd_livetime);
    systs.push_back(GetNus18Systs(rhc, "fd", syst_type));
    // syst_types.push_back(polarity + "_FD");
  }

  PredictionConcat* sim = new PredictionConcat(predictions, pots, livetimes);
  sim->SetSysts(new SystConcat(samples, {}, systs));
  sim->SetBinning(GetNus18Binning(rhc, detector));

  return sim;
}
*/
//===========================================================================================
/// Function to get default cosmics
std::vector<TH1D*> GetCosmics(bool rhc, std::string detector)
{
  bool use_nd = false;
  bool use_fd = false;
  if (detector == "nd") use_nd = true;
  else if (detector == "fd") use_fd = true;
  else if (detector == "both") {
    use_nd = true;
    use_fd = true;
  }
  else throw std::runtime_error("Detector string \"" + detector + "\" not recognised!");

  std::vector<TH1D*> cosmics;
  if (use_nd) cosmics.push_back(MakeTH1D(UniqueName().c_str(), "", kNus18EnergyBinning));
  if (use_fd) cosmics.push_back(LoadCosmicHist(rhc));

  return cosmics;
}

//===========================================================================================
/// Function to combine covariance matrices
void CombineMatrices(CovarianceMatrix* gen,
  std::map<const ISyst*, CovarianceMatrix*> m, std::vector<IPrediction*> preds,
  std::vector<double> pot, osc::IOscCalcAdjustable* calc)
{
  if (m.size() == 0) throw std::runtime_error("Matrix map is empty!");
  else if (m.size() == 1) throw std::runtime_error("Only one matrix in map; cannot add!");
  for (auto matrix : m) gen->AddMatrix(matrix.second);
  gen->PredictCovMx(preds, pot, calc);
}

//===========================================================================================
/// Function to get default full covariance matrix
/*CovarianceMatrix* GetFullCovMx(PredictionConcat* sim,
  osc::OscCalcSterile* calc, double pot, std::vector<bool> use_detector)
{
  DontAddDirectory guard;

  std::string detector_name;

  if (use_detector[0] && use_detector[1]) detector_name = "Both";
  else if (use_detector[1]) detector_name = "FD";
  else if (use_detector[0]) detector_name = "ND";

  // Load all the covariance matrices

  std::map<const ISyst*, CovarianceMatrix*> m;
  for (const ISyst* syst : GetNus18Systs(false, "base", "all")) {
    std::string systname = syst->ShortName().substr(syst->ShortName().find_last_of("_")+1);
    CovarianceMatrix* gen = LoadCovMx("covmx_nc_fhc", "CovMx_" + detector_name + "_" + systname);
    m.insert(std::pair<const ISyst*, CovarianceMatrix*>(syst, gen));
  }
  CovarianceMatrix* ppfx = LoadCovMx("covmx_nc_fhc", "CovMx_" + detector_name + "_PPFX");
  m.insert(std::pair<const ISyst*, CovarianceMatrix*>(nullptr, ppfx));

  std::vector<Binning> b = sim->GetBinning();

  // Generate new empty matrix

  CovarianceMatrix* matrix_combined = new CovarianceMatrix(pot, b);

  // Add all the matrices together

  CombineMatrices(matrix_combined, m, sim, calc);

  // Clean up the individual matrices

  for (auto matrix : m) delete matrix.second;

  return matrix_combined;
}
*/
//===========================================================================================
/// Get fake data for a set of samples


//===========================================================================================
/// Get cosmics for a set of samples
std::vector<TH1D*> GetCosmics(std::vector<covmx::Sample> samples,
  PredictionConcat* sim)
{
  DontAddDirectory guard;

  std::map<covmx::Polarity, std::string> pname;
  pname[covmx::Polarity::kFHC] = "fhc";
  pname[covmx::Polarity::kRHC] = "rhc";

  std::map<covmx::Quantile, std::string> qname;
  qname[covmx::Quantile::kQ1] = "q1";
  qname[covmx::Quantile::kQ2] = "q2";
  qname[covmx::Quantile::kQ3] = "q3";
  qname[covmx::Quantile::kQ4] = "q4";

  assert(samples.size() == sim->GetNSamples() and
    "Sample vector and concatenated prediction have different size!");

  // Open only the cosmic files we need, with a view to efficiency when running on grid

  std::string path = "/pnfs/nova/persistent/users/jhewes15/cosmic/";
  std::vector<TFile*> f(4, nullptr);
  for (const covmx::Sample& s : samples) {

    // NC cosmics

    if (s.analysis == covmx::Analysis::kNC and !f[0]) {
      std::string filepath = path + "cosmics_nus.root";
      if (RunningOnGrid()) filepath = pnfs2xrootd(filepath);
      f[0] = TFile::Open(filepath.c_str(), "read");
    }

    // CCNue cosmics

    if (s.analysis == covmx::Analysis::kCCNue && !f[1]) {
      std::string filepath = path + "cosmics_nue.root";
      if (RunningOnGrid()) filepath = pnfs2xrootd(filepath);
      f[1] = TFile::Open(filepath.c_str(), "read");
    }

    // CCNumu FHC cosmics

    if (s.analysis == covmx::Analysis::kCCNumu &&
        s.polarity == covmx::Polarity::kFHC && !f[2]) {
      std::string filepath = path + "cosmics_fhc_numu.root";
      if (RunningOnGrid()) filepath = pnfs2xrootd(filepath);
      f[2] = TFile::Open(filepath.c_str(), "read");
    }

    // CCNumu RHC cosmics

    if (s.analysis == covmx::Analysis::kCCNumu &&
        s.polarity == covmx::Polarity::kRHC && !f[3]) {
      std::string filepath = path + "cosmics_rhc_numu.root";
      if (RunningOnGrid()) filepath = pnfs2xrootd(filepath);
      f[3] = TFile::Open(filepath.c_str(), "read");
    }

  }

  // Now actually get the cosmic spectra

  std::vector<TH1D*> ret(samples.size(), nullptr);

  for (size_t i = 0; i < samples.size(); ++i) {

    // If it's a near detector spectrum, we need an empty hist for cosmics

    if (samples[i].detector == covmx::Detector::kNearDet)
      ret[i] = MakeTH1D(UniqueName().c_str(), "", sim->GetBinningOriginal()[i]);

    // Otherwise we need to load it. The cosmics are not currently standardised
    // between analyses, so we need to handle each differently

    else {

      // NC cosmics

      if (samples[i].analysis == covmx::Analysis::kNC) {
        std::string hname = pname[samples[i].polarity] + "cos";
        ret[i] = (TH1D*)f[0]->Get(hname.c_str());
      }

      // CC nue cosmics

      else if (samples[i].analysis == covmx::Analysis::kCCNue) {
        std::string sname = "cosmic_spect_" + pname[samples[i].polarity];
        std::unique_ptr<Spectrum> s(Spectrum::LoadFrom(f[1], sname.c_str()));
        ret[i] = s->ToTH1(sim->GetLivetimes()[i], kLivetime);
      }

      // CC numu FHC cosmics

      else if (samples[i].analysis == covmx::Analysis::kCCNumu
            && samples[i].polarity == covmx::Polarity::kFHC) {
        std::string hname = "cosmics_" + qname[samples[i].quantile];
        ret[i] = (TH1D*)f[2]->Get(hname.c_str());
      }

      // CC numu RHC cosmics

      else if (samples[i].analysis == covmx::Analysis::kCCNumu
            && samples[i].polarity == covmx::Polarity::kRHC) {
        std::string hname = "cosmics_" + qname[samples[i].quantile];
        ret[i] = (TH1D*)f[3]->Get(hname.c_str());
      }

      // By this point we should have instantiated the cosmic spectrum

      assert(ret[i] && "Found uninitialised cosmic spectrum!");

    } // if far detector
  } // for sample

  // Clean up files

  for (TFile* tmp : f) {
    if (tmp) delete tmp;
  }

  return ret;

} // function GetCosmics

//===========================================================================================
/// Get systematics for a given sample
std::vector<const ISyst*> GetSystsFromFile(covmx::Sample sample)
{
  // Only supported for NC right now

  if (sample.analysis != covmx::Analysis::kNC or sample.polarity != covmx::Polarity::kFHC)
    assert(false and "Only NC FHC is supported right now!");

  std::string detector = sample.detector == covmx::Detector::kNearDet ? "neardet" : "fardet";
  std::string name = "systmaker_nc_" + detector;

  std::string filepath = "/pnfs/nova/persistent/users/jhewes15/nus19/isysts/isysts_" + detector + ".root";
  if (RunningOnGrid()) filepath = pnfs2xrootd(filepath);

  std::cout << "Loading systs for sample " << name << ":" << std::endl;

  std::cout << "Loading systs for sample " << name.str() << ":" << std::endl;

  TFile* f = TFile::Open(filepath.c_str(), "read");
  TDirectory* dir = f->GetDirectory(name.c_str());

  std::vector<const ISyst*> systs;

  TIter next(f->GetListOfKeys());
  TKey* key;
  while ((key = (TKey*)next())) {
    systs.push_back(NusISyst::LoadFrom((TDirectory*)key->ReadObj()).release());
  }

  for (auto syst : systs)
    std::cout << "  " << syst->ShortName() << std::endl;
  std::cout << std::endl;

  return systs;

} // function GetSystsFromFile

//===========================================================================================
/// Default legend
TLegend* MakeLegendNus18(double x1, double y1, double x2, double y2, double textsize = 0.03)
{ 
  TLegend* leg = new TLegend(x1, y1, x2, y2);
  leg->SetBorderSize(0);
  leg->SetFillColor(kWhite);
  leg->SetFillStyle(0);
  leg->SetFillStyle(0);
  leg->SetTextFont(42);
  leg->SetTextSize(textsize);
  return leg;
}

//===========================================================================================
/// Add text to plot with delta m^2 41 at 0.5 eV^2
void PlotTextNus18(const double xpos, const double start, const double step)
{
  TLatex* tex = new TLatex();
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.037);
  tex->SetTextAlign(11);
  tex->DrawLatex(xpos, start, "Neutrino Beam:");
  tex->DrawLatex(xpos, start - 1*step, "8.0#times10^{20} ND POT,");
  tex->DrawLatex(xpos, start - 2*step, "8.85#times10^{20} FD POT-equiv.");
  tex->DrawLatex(xpos, start - 3*step, "#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
  tex->DrawLatex(xpos, start - 4*step, "#theta_{13} = 8.3#circ, ^{}#theta_{23} = 48.4#circ");
  tex->DrawLatex(xpos, start - 5*step, "#Deltam^{2}_{41} = 0.5 eV^{2}");
}

//===========================================================================================
/// Add text to plot with delta m^2 41 at custom value
void PlotTextNus18Dm41(const double xpos, const double start, const double step, std::string dm)
{
  TLatex* tex = new TLatex();
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.037);
  tex->SetTextAlign(11);
  tex->DrawLatex(xpos, start, "Neutrino Beam:");
  tex->DrawLatex(xpos, start - 1*step, "8.0#times10^{20} ND POT,");
  tex->DrawLatex(xpos, start - 2*step, "8.85#times10^{20} FD POT-equiv.");
  tex->DrawLatex(xpos, start - 3*step, "#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
  tex->DrawLatex(xpos, start - 4*step, "#theta_{13} = 8.3#circ, ^{}#theta_{23} = 48.4#circ");
  tex->DrawLatex(xpos, start - 5*step, Form("#Deltam^{2}_{41} = %s eV^{2}",dm.c_str()));
}

//===========================================================================================
/// Add text to plot without delta m^2 41 value
void PlotTextNus18WithoutDm41(const double xpos, const double start, const double step)
{
  TLatex* tex = new TLatex();
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.037);
  tex->SetTextAlign(11);
  tex->DrawLatex(xpos, start, "Neutrino beam:");
  tex->DrawLatex(xpos, start - 1*step, "8.0#times10^{20} ND POT,");
  tex->DrawLatex(xpos, start - 2*step, "8.85#times10^{20} FD POT-equiv.");
  tex->DrawLatex(xpos, start - 3*step, "#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
  tex->DrawLatex(xpos, start - 4*step, "#theta_{13} = 8.3#circ, ^{}#theta_{23} = 48.4#circ");
}

//===========================================================================================
/// Gaussian constrains on atmospheric parameters
std::vector<const IExperiment*> GetNus18Constraints()
{
  std::vector<const IExperiment*> ret;
  ret.push_back(new const GaussianConstraint(&kFitTheta23Sterile, kNus18Th23, kNus18Th23Sigma));
  ret.push_back(new const GaussianConstraint(&kFitDmSq32Sterile, kNus18Dm32, kNus18Dm32Sigma));
  return ret;
}

//===========================================================================================
/// Seed values for Nus18 marginalisation parameters
std::vector<double> GetNus18SeedValues(osc::IOscCalcAdjustable* calc)
{
  // Sterile parameters
  kFitDmSq41Sterile.SetValue(calc,  0.01);
  kFitSinSqTheta24Sterile.SetValue(calc, 0.01);
  kFitSinSqTheta34Sterile.SetValue(calc, 0.01);
  kFitDelta24InPiUnitsSterile.SetValue(calc, kNus18Delta24);

  // Atmospheric parameters
  kFitTheta23Sterile.SetValue(calc, kNus18Th23);
  kFitDmSq32Sterile.SetValue(calc, kNus18Dm32);

  // Return seed values
  std::vector<double> ret;
  ret.push_back(kFitTheta23Sterile.GetValue(calc));
  ret.push_back(kFitDmSq32Sterile.GetValue(calc));
  ret.push_back(kFitDelta24InPiUnitsSterile.GetValue(calc));

  return ret;
}

//===========================================================================================
std::vector<const IFitVar*> GetNus18FitVars()
{
  std::vector<const IFitVar*> ret;
  ret.push_back(&kFitTheta23Sterile);
  ret.push_back(&kFitDmSq32Sterile);
  ret.push_back(&kFitDelta24InPiUnitsSterile);
  return ret;
}

//===========================================================================================
std::vector<double> GetEdges(TH1* h)
{
  int nbins = h->GetNbinsX();
  std::vector<double> edges(nbins+1);
  for (int i = 0; i <= nbins; ++i)
    edges[i] = h->GetXaxis()->GetBinLowEdge(i+1);
  return edges;
}

//===========================================================================================
TH1D* DiagonalErrors(TH2D* m)
{
  int nbins = m->GetNbinsX();
  std::vector<double> edges = GetEdges((TH1*)m);
  TH1D* h = new TH1D(UniqueName().c_str(),"",nbins,&edges[0]);
  for (int i = 1; i <= nbins; ++i) {
    h->SetBinContent(i,1);
    double var = m->GetBinContent(i,i);
    if (var > 0) h->SetBinError(i,sqrt(var));
    else h->SetBinError(i,0);
  }
  return h;
}

//===========================================================================================
std::vector<std::tuple<ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t>> GetComponents()
{
  std::vector<std::tuple<ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t> > components;
  std::vector<ana::Flavors::Flavors_t> flavors = {ana::Flavors::kNuEToNuE, ana::Flavors::kNuEToNuMu,
        ana::Flavors::kNuEToNuTau, ana::Flavors::kNuMuToNuE,
        ana::Flavors::kNuMuToNuMu, ana::Flavors::kNuMuToNuTau};
  std::vector<ana::Sign::Sign_t> signs = {ana::Sign::kNu, ana::Sign::kAntiNu};
  for (auto & flavor: flavors) {
    for (auto & sign: signs) {
      components.push_back( std::make_tuple(flavor, ana::Current::kCC, sign) );
    }
  }
  components.push_back( std::make_tuple(ana::Flavors::kAll, ana::Current::kNC, ana::Sign::kBoth));
  return components;
}

//------------------------------------------------------------------------------
TLatex* MiscText(double x, double y, double size, TString text)
{
  TLatex *l = new TLatex(x, y, text);
  l->SetNDC();
  l->SetTextSize(size);
  l->SetTextColor(kBlack);
  l->Draw();
  return l;
}

//===========================================================================================
std::string GetNus18ComponentName(int i)
{
  if      (i==0) return std::string("#nu_{e} #rightarrow #nu_{e}");
  else if (i==1) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{e}");
  else if (i==2) return std::string("#nu_{e} #rightarrow #nu_{#mu}");
  else if (i==3) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{#mu}");
  else if (i==4) return std::string("#nu_{e} #rightarrow #nu_{#tau}");
  else if (i==5) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{#tau}");
  else if (i==6) return std::string("#nu_{#mu} #rightarrow #nu_{e}");
  else if (i==7) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{e}");
  else if (i==8) return std::string("#nu_{#mu} #rightarrow #nu_{#mu}");
  else if (i==9) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{#mu}");
  else if (i==10) return std::string("#nu_{#mu} #rightarrow #nu_{#tau}");
  else if (i==11) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{#tau}");
  else if (i==12) return std::string("NC");
  std::cerr << "What are you doing???" << std::endl;
  return "no";
}

//===========================================================================================
void FormatFullCovMxPlot(TH2D* h, PredictionConcat* sim)
{
  double last_edge = sim->GetBinningConcatenated().Edges().back();
  double max = 13. * last_edge;
  h->Draw("colz");

  // Draw grey lines
  for (int i = 0; i < 13; ++i) {
    double pos = ((double)i + 0.5) * last_edge;
    TLine* lx = new TLine(0, pos, max, pos);
    lx->SetLineColor(kGray+1);
    lx->Draw();
    TLine* ly = new TLine(pos, 0, pos, max);
    ly->SetLineColor(kGray+1);
    ly->Draw();
  }

  // Draw black lines
  for (int i = 0; i <= 13; ++i) {
    double pos = i * last_edge;
    TLine* lx = new TLine(0, pos, max, pos);
    lx->Draw();
    TLine* ly = new TLine(pos, 0, pos, max);
    ly->Draw();
  }

  // Component names
  for (int i = 0; i < 13; ++i) {
    double pos = 0.16 + (0.72*(i/13.));
    auto l1 = MiscText(pos, 0.88, 0.04, GetNus18ComponentName(i).c_str());
    l1->SetTextAngle(38);
    auto l2 = MiscText(0.12, pos, 0.04, GetNus18ComponentName(i).c_str());
    l2->SetTextAlign(31);
  }

  h->SetNdivisions(0, "X");
  h->SetNdivisions(0, "Y");
}

//===========================================================================================
bool CheckOption(TString opts, TString opt)
{
  if (opts.Contains(opt, TString::ECaseCompare::kIgnoreCase))
    return true;
  else return false;
}

//===========================================================================================
std::vector<bool> GetOptionConfig(TString opt, std::vector<std::string> options)
{
  size_t n_opts = options.size();

  std::vector<bool> vec(n_opts);
  for (size_t i = 0; i < n_opts; ++i)
    if (CheckOption(opt, options[i]))
      vec[i] = true;

  if (std::none_of(vec.begin(), vec.end(), [](bool v) { return v; })) {
    std::string err = "Option string must contain ";
    for (size_t i = 0; i < n_opts-2; ++i)
      err += "\"" + options[i] + "\", ";
    err += options[n_opts-2] + " or " + options[n_opts-1] + ".";
    throw std::runtime_error(err);
  }

  else {
    std::cout << "Running for";
    for (size_t i = 0; i < n_opts; ++i)
      if (vec[i])
        std::cout << " " << options[i];
    std::cout << std::endl;
  }

  return vec;

} // function GetOptionConfig

//===========================================================================================
size_t GetOptionValue(TString opt, std::string key)
{
  size_t val = 999;

  std::string optstring = std::string(opt.Data());
  size_t pos = optstring.find(key);
  if (pos != std::string::npos) {
    if (optstring.substr(pos+key.size(), 1) != "=") {
      std::ostringstream err;
      err << "The " << key << " option must take the form \"" << key << "=xxx\".";
      throw std::runtime_error(err.str());
    }
    size_t start = pos + key.size() + 1;
    size_t end = optstring.find("_", pos);
    if (end != std::string::npos)
      val = std::stod(optstring.substr(start, end - start));
    else
      val = std::stod(optstring.substr(start));

    std::cout << "Using " << key << " value " << val << std::endl;
    return val;
  }
  std::ostringstream err;
  err << "No value for " << key << " found in option string.";
  throw std::runtime_error(err.str());

} // function GetOptionValue

//===========================================================================================
// Take a concatenated fake data spectrum and split it up
std::vector<TH1D*> SplitFakeData(TH1D* h, std::vector<Binning> b)
{
  // Define split-up histograms

  std::vector<TH1D*> h_split;
  size_t n_bins = 0;
  for (size_t i = 0; i < b.size(); ++i) {
    h_split.push_back(MakeTH1D(UniqueName().c_str(), "", b[i]));
    n_bins += b[i].NBins();
  }

  // Check binning is consistent

  if (h->GetNbinsX() != int(n_bins))
    throw std::runtime_error("There are " + std::to_string(h->GetNbinsX())
      + " bins in merged histogram but a total of " + std::to_string(n_bins)
      + " bins across the binning vector provided.");

  int bin_count = 0;

  for (size_t i = 0; i < b.size(); ++i) {
    for (int j = 1; j <= h_split[i]->GetNbinsX(); ++j) {
      double val = h->GetBinContent(bin_count+j);
      h_split[i]->SetBinContent(j, val);
      h_split[i]->SetBinError(j, sqrt(val));
      // std::cout << "Full hist bin " << bin_count+j
      //   << " with value " << val << " is being assigned to bin "
      //   << j << " of split histogram " << i+1 << std::endl;
    }
    bin_count += h_split[i]->GetNbinsX();
  }

  return h_split;

} // function SplitFakeData

//===========================================================================================
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)
{
  DontAddDirectory guard;

  std::vector<int> colours = { kAzure-2, kRed+2, kGreen+3, kMagenta+2 };
  int cosmic_colour = kGray+2;

  // Get the number of detectors and universes

  if (calcs.empty()) throw std::runtime_error("Spectrum vector is empty.");
  size_t n_universes = calcs.size();
  if (n_universes > colours.size())
    throw std::runtime_error("Only " + std::to_string(colours.size())
      + " universes per plot are supported!");
  size_t n_detectors = sim->GetNSamples();
  if (pot.size() != n_detectors)
    throw std::runtime_error(std::to_string(n_detectors) + " detectors but only "
      + std::to_string(pot.size()) + " elements in POT vector.");
  if (is_nd.size() != n_detectors)
    throw std::runtime_error(std::to_string(n_detectors) + " detectors but only "
      + std::to_string(is_nd.size()) + " elements in is_nd vector.");

  // Use OscCalcs and simulation to get spectra

  std::vector<std::vector<Spectrum>> specs;
  for (size_t i = 0; i < n_universes; ++i)
    specs.push_back(sim->GetSpectra(calcs[i].second));

  // Split up data histogram

  std::vector<TH1D*> data   = SplitFakeData(data_merged,   sim->GetBinning());
  std::vector<TH1D*> cosmic = SplitFakeData(cosmic_merged, sim->GetBinning());

  // Define shape of canvas

  TCanvas* c = new TCanvas(UniqueName().c_str(), "c", 500 * n_detectors, 800);
  c->Divide(n_detectors, 3);

  // Pad margins

  // std::vector<float> lm = { 0.15, 0.12 };
  // std::vector<float> rm = { 0.05, 0.08 };
  std::vector<float> margin = { 0.1, 0.05 };

  gStyle->SetPadLeftMargin(margin[0]);
  gStyle->SetPadRightMargin(margin[1]);

  gStyle->SetLabelFont(43, "xyz");
  gStyle->SetLabelSize(14, "xyz");

  for (size_t d = 0; d < n_detectors; ++d)
  {
    // Get pad x ranges for this detector

    float xmin = float(d) / float(n_detectors);
    float xmax = float(d+1) / float(n_detectors);

    // Set size of spectrum pad

    c->cd(d+1);
    gPad->SetPad(xmin, 0.45, xmax, 1);

    // Set size of ratio pad

    c->cd(n_detectors+d+1);
    gPad->SetPad(xmin, 0.25, xmax, 0.45);

    // Set size of oscillation pad

    c->cd((2*n_detectors)+d+1);
    gPad->SetPad(xmin, 0.05, xmax, 0.25);

  } // for detector i

  TLegend* l = new TLegend(0.5, 0.65, 0.85, 0.85);
  l->SetFillStyle(0);

  // Create histogram vectors and set their sizes

  std::vector<std::vector<TH1D*>> hists, ratio_hists;

  hists.resize(n_universes);
  for (std::vector<TH1D*>& h : hists)
    h.resize(n_detectors);

  ratio_hists.resize(n_universes);
  for (std::vector<TH1D*>& h : ratio_hists)
    h.resize(n_detectors);

  std::vector<TGraph*> graphs;

  // Y axis max for each detector's plot

  std::vector<double> ymax, rmax;
  std::vector<bool> rescaled;

  // Make data ratio histograms

  std::vector<TH1D*> data_ratio;
  data_ratio.resize(n_detectors);
  for (size_t d = 0; d < n_detectors; ++d)
  {
    data_ratio[d] = (TH1D*)data[d]->Clone();
    for (int b = 1; b <= data_ratio[d]->GetNbinsX(); ++b)
    {
      double val = data[d]->GetBinContent(b);
      double err = data[d]->GetBinError(b);
      data_ratio[d]->SetBinContent(b, 1);

      if (val == 0)
        data_ratio[d]->SetBinError(b, 0);
      else
        data_ratio[d]->SetBinError(b, err / val);

    } // for bin b
  } // for detector d

  // Loop over each detector

  for (size_t d = 0; d < n_detectors; ++d)
  {
    // Loop over each universe, get histograms and get ymax

    std::vector<double> m, r;

    // Rescale data histogram

    if (is_nd[d]) data[d]->Scale(1e-4);

    for (size_t u = 0; u < n_universes; ++u)
    {
      hists[u][d] = (specs[u][d].ToTH1(pot[d]));

      // If it's a near detector spectrum, rescale it

      if (is_nd[d]) hists[u][d]->Scale(1e-4);
      m.push_back(hists[u][d]->GetMaximum());

      // Clone to create ratio histogram

      ratio_hists[u][d] = (TH1D*)hists[u][d]->Clone();
      TH1D* h = ratio_hists[u][d];
      for (int b = 1; b <= h->GetNbinsX(); ++b)
      {
        double val = (hists[u][d]->GetBinContent(b) + cosmic[d]->GetBinContent(b)) / data[d]->GetBinContent(b);
        if (std::isinf(val) || std::isnan(val)) {
          h->SetBinContent(b, 1);
          h->SetBinError(b, 0);
        }
        else {
          h->SetBinContent(b, val);
          r.push_back(fabs(1-val));
        }

      } // for bin b

    } // for universe u

    // Get maximum data point

    int maxdata = data[d]->GetMaximumBin();
    m.push_back(data[d]->GetBinContent(maxdata) + data[d]->GetBinError(maxdata));

    // Account for data errors in ratio plot range

    for (int b = 1; b <= data_ratio[d]->GetNbinsX(); ++b)
      r.push_back(fabs(data_ratio[d]->GetBinError(b)));

    ymax.push_back(*std::max_element(m.begin(), m.end()));
    rmax.push_back(*std::max_element(r.begin(), r.end()));

  } // for detector d

  // Formatting options

  float ts = 0.055;

  // Loop over each detector to draw the cosmics

  for (size_t d = 0; d < n_detectors; ++d) {

    // Draw the cosmic spectrum

    c->cd(d+1);
    TH1D* hc = cosmic[d];
    hc->SetLineColor(cosmic_colour);
    hc->SetLineWidth(2);
    hc->SetMinimum(0);
    hc->SetMaximum(1.1 * ymax[d]);
    hc->Draw("hist same");

  }

  // Loop over each universe

  for (size_t u = 0; u < n_universes; ++u)
  {
    // Loop over each detector

    for (size_t d = 0; d < n_detectors; ++d)
    {
      // Format and draw the spectrum

      c->cd(d+1);
      TH1D* h = hists[u][d];
      h->Add(cosmic[d]);
      h->SetLineColor(colours[u]);
      h->SetLineWidth(2);
      h->SetMinimum(0);
      h->SetMaximum(1.1 * ymax[d]);

      if (is_nd[d])
        h->SetTitle(";Energy deposited in detector (GeV);10^{4} events");
      else
        h->SetTitle(";Energy deposited in detector (GeV);Events");
      h->GetXaxis()->SetTitleSize(ts);
      h->GetYaxis()->SetTitleSize(ts);

      h->Draw("hist same");

      // Now draw the ratio plot

      h = ratio_hists[u][d];
      c->cd(n_detectors+d+1);
      h->SetLineColor(colours[u]);
      h->SetLineWidth(2);
      h->SetMinimum(1 - (1.1 * rmax[d]));
      h->SetMaximum(1 + (1.1 * rmax[d]));

      h->SetTitle(";;Ratio");
      h->GetYaxis()->SetTitleSize(3*ts);

      h->Draw("hist same");

      // Now the oscillation plot

      c->cd((2*n_detectors)+d+1);

      std::vector<float> x, y;
      double l = calcs[u].second->GetL();
      if (is_nd[d]) calcs[u].second->SetL(1);
      else calcs[u].second->SetL(810);
      for (size_t k = 0; k < 200; ++k)
      {
        x.push_back(0.1 * (float(k) + 0.5));
        y.push_back(calcs[u].second->P(14, 0, x.back()));
      }
      calcs[u].second->SetL(l);

      TGraph* g = new TGraph(x.size(), &x[0], &y[0]);

      g->SetLineWidth(2);
      g->SetLineColor(colours[u]);

      g->SetTitle(";;P(#nu_{#mu} #rightarrow #nu_{s})");
      g->GetYaxis()->SetTitleSize(3*ts);

      g->Draw("ac");
      graphs.push_back(g);

    } // for detector d

    l->AddEntry(hists[u][0], calcs[u].first.c_str(), "l");

  } // for universe u

  l->AddEntry(cosmic[0], "Cosmic background", "l");

  if (data.size() != n_detectors)
    throw std::runtime_error("Number of data hists does not match number of samples!");

  // Loop and draw the data for each sample

  for (size_t d = 0; d < n_detectors; ++d)
  {
    c->cd(d+1);
    data[d]->SetMarkerSize(2);
    data[d]->Draw("e0 same");

    c->cd(n_detectors+d+1);
    data_ratio[d]->SetMarkerSize(2);
    data_ratio[d]->Draw("e0 same");

  } // for detector d

  l->AddEntry(data[0], "Data", "lep");
  c->cd(1);
  l->Draw();

  dir->WriteTObject(c, name.c_str());

  // std::string img = name + ".png";
  // c->SaveAs(img.c_str());

  // Now clean up

  delete c;
  delete l;
  for (auto g : graphs)
    delete g;
  for (auto h_vec : hists)
    for (auto h : h_vec)
      delete h;
  for (auto h_vec : ratio_hists)
    for (auto h : h_vec)
      delete h;
  for (auto h : data_ratio)
    delete h;

} // function DrawSurfacePoint