MakeCovarSim.C

Go to the documentation of this file.

/// MakeCovarSim.C by J. Hewes <jhewes15@fnal.gov>

#include "CAFAna/Extrap/NDExtrap.h"
#include "CAFAna/Extrap/FDExtrap.h"

#include "CAFAna/Prediction/PredictionGenerator.h"
#include "NuXAna/Prediction/PredictionGeneratorNuX.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "NuXAna/Prediction/NDPredictionSterile.h"
#include "NuXAna/Prediction/FDPredictionSterile.h"

#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "NuXAna/Cuts/NusCuts18.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"

#include "NuXAna/Vars/NusVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "NuXAna/Vars/HistAxes.h"

#include "CAFAna/Systs/XSecSystLists.h"
#include "NuXAna/Systs/Nus18Systs.h"

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Prod4Loaders.h"
#include "3FlavorAna/Systs/JointAna2018Systs.h"
#include "CAFAna/Core/Sample.h"

#include "NuXAna/macros/Nus19FHC/Utilities.h"

#include "OscLib/OscCalcSterile.h"

using namespace ana;

void MakeCovarSim(TString opt)
{
  DontAddDirectory guard;

  // Set detector and analysis type

  bool nd;
  covmx::Detector c_det;
  if (CheckOption(opt, "neardet")) {
    nd = true;
    c_det = covmx::Detector::kNearDet;
  }
  else if (CheckOption(opt, "fardet")) {
    nd = false;
    c_det = covmx::Detector::kFarDet;
  }
  else assert(false and "Option string must include detector: \"neardet\" or \"fardet\".");
  std::string detector = nd ? "neardet" : "fardet";

  std::string ana;
  covmx::Analysis c_ana;
  if (CheckOption(opt, "nue")) {
    ana = "nue";
    c_ana = covmx::Analysis::kCCNue;
  }
  else if (CheckOption(opt, "numu")) {
    ana = "numu";
    c_ana = covmx::Analysis::kCCNumu;
  }
  else if (CheckOption(opt, "nus")) {
    ana = "nus";
    c_ana = covmx::Analysis::kNC;
  }
  else assert(false and "Option string must include analysis: \"nue\", \"numu\" or \"nus\".");

  bool rhc;
  covmx::Polarity c_pol;
  if (CheckOption(opt, "fhc")) {
    rhc = false;
    c_pol = covmx::Polarity::kFHC;
  }
  else if (CheckOption(opt, "rhc")) {
    rhc = true;
    c_pol = covmx::Polarity::kRHC;
  }
  else assert(false and "Option string must include beam polarity: \"fhc\" or \"rhc\".");
  std::string polarity = rhc ? "rhc" : "fhc";

  covmx::Sample sample(c_ana, c_pol, c_det);

  bool concat;
  if (CheckOption(opt, "fullcaf")) concat = false;
  else if (CheckOption(opt, "concat")) concat = true;
  else assert(false and "Option string must include file type: \"concat\" or \"fullcaf\".");

  std::string syst_type;
  if (CheckOption(opt, "nosysts")) syst_type = "nosysts";
  else if (CheckOption(opt, "filesysts")) syst_type = "systs";
  else if (CheckOption(opt, "geniesysts")) syst_type = "genie";
  else assert(false and "Option string must include syst type: \"nosysts\", \"filesysts\" or \"geniesysts\".");

  std::string optstring = std::string(opt.Data());

  size_t quantile = 999;
  if (ana == "numu") {
    size_t qpos = optstring.find("quantile");
    if (qpos != std::string::npos) {
      if (optstring.substr(qpos+8, 1) != "=")
        assert(false and "The quantile option must take the form \"quantile=xxx\".");
      size_t start = qpos+9;
      size_t end = optstring.find("_", qpos);
      if (end != std::string::npos)
        quantile = std::stod(optstring.substr(start, end-start));
      else
        quantile = std::stod(optstring.substr(start));
      std::cout << "Using quantile value " << quantile << std::endl;
    }
  }

  size_t usesyst = 999;
  size_t spos = optstring.find("usesyst");
  if (spos != std::string::npos) {
    if (optstring.substr(spos+7, 1) != "=")
      assert(false and "The usesyst option must take the form \"usesyst=xxx\".");
    size_t start = spos+8;
    size_t end = optstring.find("_", spos);
    if (end != std::string::npos)
      usesyst = std::stod(optstring.substr(start, end-start));
    else
      usesyst = std::stod(optstring.substr(start));
    std::cout << "Using systematic in position " << usesyst << std::endl;
  }
  // else assert(false && "You must provide a \"usesyst\" option!");

  // Important variables
  const Var kReweight = kXSecCVWgt2018 * kPPFXFluxCVWgt;
  osc::OscCalcSterile* calc = DefaultSterileCalc(4);

  // Set up systematics
  std::vector<const ISyst*> all_systs, systs;
  if (syst_type == "systs")
    all_systs = GetSystsFromFile(sample);
  if (syst_type == "genie")
    all_systs = getAllXsecSysts_2018();

  if (all_systs.size() > 0) {

    if (usesyst == 999) {
      std::cout << "No \"usesyst\" option specified, so running for all systs simultaneously!" << std::endl;
      systs = all_systs;
    }

    else if (usesyst >= all_systs.size()) {
      std::ostringstream err;
      err << "You asked for the systematic in position "
          << usesyst << ", but there are only "
          << all_systs.size() << " systematics in the vector provided.";
      assert(false and err.str().c_str());
    }

    else {
      const ISyst* syst = all_systs[usesyst];
      std::cout << "Running for systematic " << syst->ShortName() << std::endl;
      systs = { syst };
    }

  } // if allsysts not empty

  // Set CAF type
  enum ECAFType caf_type;
  if (!concat) caf_type = ECAFType::kFullCAF;
  else if (ana == "nue")  caf_type = ECAFType::kNueConcat;
  else if (ana == "numu") caf_type = ECAFType::kNumuConcat;
  else if (ana == "nus")  caf_type = ECAFType::kNusConcat;
  else assert(false and "Couldn't find valid caf type.");

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

  // Set axes and syst type
  const HistAxis* axis = nullptr;
  const Cut* cut = nullptr;

  // Nue axis and cuts
  if (ana == "nue") {
    axis = new const HistAxis(kNue2018AxisMergedPeripheral);
    if (nd) cut = new const Cut(kNue2018NDCVNSsb);
    else    cut = new const Cut(kNue2018FDAllSamples);
  }

  // Numu axis and cuts
  else if (ana == "numu") {
    axis = new const HistAxis(kNumuCCOptimisedAxis);
    if (quantile > 3)
      assert(false and "For a numu prediction you must define a quantile between 0 and 3!");
    std::string quant_dir = "/pnfs/nova/persistent/users/jhewes15/quantiles";
    std::string infile_quant = quant_dir + "/quantiles__" + polarity + "_full__numu2018.root";
    std::cout << "\n\n --- get quantile cuts from file ---" << std::endl;
    TFile* infile = TFile::Open( pnfs2xrootd(infile_quant).c_str() );
    TH2* FDSpec2D = (TH2*)infile->FindObjectAny("FDSpec2D");
    std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
    infile->Close();
    delete infile;
    if (nd) cut = new const Cut(kNumuCutND2018 && HadEFracQuantCuts[quantile]);
    else    cut = new const Cut(kNumuCutFD2018 && HadEFracQuantCuts[quantile]);
  }

  // Neutral current axis and cuts

  else if (ana == "nus") {
    if (nd) axis = new const HistAxis(kNCNDAxisE);
    else    axis = new const HistAxis(kNCFDAxisE);
    if (nd) cut = new const Cut(kNus18ND);
    else    cut = new const Cut(kNus18FD);
  }

  // Define the loader
  Prod4NomLoaders loaders(caf_type,
    rhc ? Loaders::FluxType::kRHC : Loaders::FluxType::kFHC);
  if (nd) loaders.SetND(true);
  loaders.SetSpillCut(kStandardSpillCuts);

  // Create the prediction
  IPredictionGenerator* gen = nullptr;
  if (nd) gen = new NDPredictionGenerator(*axis, *cut, kNoShift, kReweight);
  else    gen = new FDPredictionGenerator(*axis, *cut, kNoShift, kReweight);
  IPrediction* pred = gen->Generate(loaders).release();
  PredictionInterp* pred_interp = nullptr;
  if (systs.size() > 0)
    pred_interp = new PredictionInterp(systs, calc, *gen, loaders);
  loaders.Go();

  // Initialise output file
  std::string systname;
  if (systs.size() == 1)
    systname = systs[0]->ShortName();
  else if (systs.size() == 0)
    systname = "nosysts";
  else
    systname = syst_type;

  std::ostringstream f_out_name;
  f_out_name << "pred_" << ana << "_" << polarity << "_" << detector;
  if (ana == "numu") f_out_name << "_q" << quantile;
  f_out_name << "_" << systname << ".root";
  TFile* f_out = TFile::Open(f_out_name.str().c_str(), "recreate");

  // Save prediction to file
  std::ostringstream dirname;
  dirname << "pred_" << ana << "_" << polarity << "_" << detector;
  if (ana == "numu") dirname << "_q" << quantile;
  dirname << "_" << systname;

  pred->SaveTo(f_out, dirname.str().c_str());

  if (pred_interp) {
    dirname.str("");
    dirname << "pred_interp_" << ana << "_" << polarity << "_" << detector;
    if (ana == "numu") dirname << "_q" << quantile;
    dirname << "_" << systname;
    pred_interp->SaveTo(f_out, dirname.str().c_str());
 }

} // function MakeCovarSim