MakePPFXRatios.C

Go to the documentation of this file.

#ifdef __CINT__
void MakePPFXRatios(std::string detector, int job,
  bool concat = true, bool rhc = false) {
  std::cout<< "Sorry, you gotta run in compiled mode." << std::cout;
}
#else

#include "NuXAna/Cuts/NusCuts18.h"
#include "CAFAna/Analysis/Prod4Loaders.h"
#include "NuXAna/macros/Nus18/Nus18Utilities.h"

using namespace ana;

void MakePPFXRatios(std::string detector, int job,
  bool concat = true, bool rhc = false)
{
  DontAddDirectory guard;

  bool nd;
  if (detector == "nd") nd = true;
  else if (detector == "fd") nd = false;
  else throw std::runtime_error("Detector string must be \"nd\" or \"fd\"; \""
    + detector + "\" is not understood.");

  if (job < 1 || job > 5) throw std::runtime_error("Job int must be in range 1 -> 5!");

  osc::OscCalcSterile* calc4f = DefaultSterileCalc(4);
  const HistAxis kAxisE("Energy Deposited in Detector (GeV)", kNus18EnergyBinning, kNus18Energy);

  std::string filename = detector + "_ppfx_" + std::to_string(job) + ".root";
  TFile* f = new TFile(filename.c_str(), "RECREATE");

  std::vector<Var> kPPFXFluxUnivWgts;
  for (unsigned int i = 0; i < 100; ++i) {
    const Var tempPPFXWgt([i](const caf::SRProxy* sr){
      if (!sr->hdr.ismc) return 1.f;
      if (sr->mc.nnu != 1) return 1.f;
      if (sr->mc.nu[0].rwgt.ppfx.vuniv[i] <= 0) return 1.f;
      return (float)sr->mc.nu[0].rwgt.ppfx.vuniv[i];
    });
    kPPFXFluxUnivWgts.push_back(tempPPFXWgt);
  }

  for (int n = 20 * (job-1); n < 20 * job; ++n) {

    std::string detname = nd? "near detector" : "far detector";
    std::cout << "Processing universe " << n << " for " << detname << "..." << std::endl;

    const Var kPPFXFluxUnivWgt = kPPFXFluxUnivWgts[n];
    const Var kReweight = kXSecCVWgt2018*kPPFXFluxUnivWgt;

    // Set up the loaders
    Loaders* loaders = new Prod4NomLoaders(
      concat? ECAFType::kNusConcat : ECAFType::kFullCAF,
      rhc? Loaders::FluxType::kRHC : Loaders::FluxType::kFHC);
    if (nd) loaders->SetND(true);
    loaders->SetSpillCut(kStandardSpillCuts);

    if (nd) {

      NDPredictionGenerator ppfx_nd(kAxisE, kNus18ND, kNoShift, kReweight);
      PredictionInterp universe({}, calc4f, ppfx_nd, *loaders);

      loaders->Go();

      universe.SaveTo(f, Form("universe_%i",n));

    } else {

      FDPredictionGenerator ppfx_fd(kAxisE, kNus18FD, kNoShift, kReweight);
      PredictionInterp universe({}, calc4f, ppfx_fd, *loaders);

      loaders->Go();

      universe.SaveTo(f, Form("universe_%i",n));

    }
    
  } // for n

  f->Close();
}

#endif