MakePrediction.C

Go to the documentation of this file.

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

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

#include "CAFAna/Prediction/PredictionGenerator.h"
#include "CAFAna/Prediction/PredictionInterp.h"

#include "CAFAna/Cuts/SpillCuts.h"

#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "CAFAna/Systs/XSecSystLists.h"

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Prod5Loaders.h"
#include "CAFAna/Core/Sample.h"

#include "NuXAna/Vars/NusVars.h"
#include "NuXAna/Prediction/NDPredictionSterile.h"
#include "NuXAna/Prediction/FDPredictionSterile.h"
#include "NuXAna/Prediction/PredictionGeneratorNuX.h"
#include "NuXAna/macros/Nus20/Utilities.h"

#include "OscLib/OscCalcSterile.h"

using namespace ana;

using std::begin;
using std::cout;
using std::end;
using std::endl;
using std::exception;
using std::ostringstream;
using std::pair;
using std::string;
using std::vector;

using covmx::Sample;

void MakePrediction(TString opt, string systFile="default", double res=0)
{
  try {
    DontAddDirectory guard;

    Sample sample = GetSampleFromOptString(opt);
    SetAxis(sample, res);
    SetCut(sample);

    // What type of file to run over
    vector<pair<string, bool>> loaderOpts
      = { {"fullcaf",false}, {"concat",true} };
    bool concat = ParseOption<bool>(loaderOpts, opt, true);

    // Systematic option
    vector<pair<string, string>> systOpts
      = { {"nosysts","nosysts"}, {"filesysts","file"}, {"xsecsysts","xsec"},
        {"allsysts","all"} };

    string systType = ParseOption<string>(systOpts, opt, "nosysts");

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

    // Set up systematics
    vector<const ISyst*> systs;
    if (systType == "file" || systType == "all") {
      vector<const ISyst*> fileSysts = systFile == "default" ?
        LoadSysts(sample) : LoadSystsFromFile(pnfs2xrootd(systFile), "isysts_"+sample.GetTag());
      systs.insert(end(systs), begin(fileSysts), end(fileSysts));
    }
    if (systType == "xsec" || systType == "all") {
      vector<const ISyst*> xsecSysts = getAllXsecSysts_2020();
      systs.insert(end(systs), begin(xsecSysts), end(xsecSysts));
    }
    // Add in neutron here eventually

    // Set CAF type
    enum ECAFType CAFType = ECAFType::kFullCAF;
    if (concat) {
      if      (sample.selection == covmx::kCCNue)  CAFType = ECAFType::kNueConcat;
      else if (sample.selection == covmx::kCCNumu) CAFType = ECAFType::kNumuConcat;
      else if (sample.selection == covmx::kNC)     CAFType = ECAFType::kNusConcat;
    }

    // Define the loader
    Prod5NomLoaders loaders(CAFType,
      sample.polarity == covmx::kRHC ? Loaders::FluxType::kRHC : Loaders::FluxType::kFHC);
    if (sample.detector == covmx::kNearDet) {
      loaders.SetND(true);
    }
    loaders.SetSpillCut(kStandardSpillCuts);

    // Create the prediction
    IPredictionGenerator* gen = nullptr;
    if (sample.detector == covmx::kNearDet) gen = new NDPredictionGenerator(*sample.GetAxis(), *sample.GetCut(), kNoShift, kReweight);
    else gen = new FDPredictionGenerator(*sample.GetAxis(), *sample.GetCut(), kNoShift, kReweight);
    IPrediction* pred = gen->Generate(loaders).release();
    PredictionInterp* predInterp = nullptr;
    if (systs.size() > 0) {
      predInterp = new PredictionInterp(systs, calc, *gen, loaders);
    }
    loaders.Go();

    // Initialise output file
    ostringstream outFileName;
    outFileName << "pred_" << sample.GetTag();
    if (res!=0) outFileName << "_res" << (int)res;
    outFileName << ".root";
    TFile* outFile = TFile::Open(outFileName.str().c_str(), "recreate");

    // Save prediction to file
    ostringstream dirname;
    dirname << "pred_" << sample.GetTag();
    pred->SaveTo(outFile->mkdir(dirname.str().c_str()));

    if (predInterp) {
      dirname.str("");
      dirname << "pred_interp_" << sample.GetTag();;
      predInterp->SaveTo(outFile->mkdir(dirname.str().c_str()));
    }
  } catch(const exception& e) {
    cout << endl << "Exception occurred! " << e.what() << endl;
  }

} // macro MakePrediction