MakePPFXRatios.C

Go to the documentation of this file.

#include "TString.h"
#include "TFile.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include "CAFAna/Core/ISyst.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Sample.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Vars/HistAxes.h"
#include "CAFAna/Vars/XsecTunes.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Analysis/Prod4Loaders.h"

#include "NuXAna/Cuts/NusCuts18.h"
#include "NuXAna/macros/Nus20/Utilities.h"

using namespace ana;

void MakePPFXRatios(TString opt) {

  covmx::Sample sample = opt.Contains("neardet", TString::kIgnoreCase) ?
    kNusFHCNearDet : kNusFHCFarDet;

  // Here get loaders
  Prod4NomLoaders loaders(ECAFType::kNusConcat, Loaders::FluxType::kFHC);
  if (sample.detector == covmx::kNearDet) loaders.SetND(true);
  loaders.SetSpillCut(kStandardSpillCuts);

  // This should surely be a vector of loaders and cuts, right?
  auto det = (sample.detector == covmx::kNearDet) ? caf::kNEARDET : caf::kFARDET;

  std::vector<std::tuple<SpectrumLoaderBase*, const Cut, std::string>> components;

  SpectrumLoaderBase* nonswap = &loaders.GetLoader(det,
    Loaders::kMC, ana::kBeam, Loaders::kNonSwap);

  // Unoscillated components (for both detectors)
  components.push_back(std::make_tuple(nonswap, kIsBeamNue && !kIsAntiNu, "cc_nue_to_nue"        ));
  components.push_back(std::make_tuple(nonswap, kIsBeamNue && kIsAntiNu,  "cc_nuebar_to_nuebar"  ));
  components.push_back(std::make_tuple(nonswap, kIsNumuCC  && !kIsAntiNu, "cc_numu_to_numu"      ));
  components.push_back(std::make_tuple(nonswap, kIsNumuCC  && kIsAntiNu,  "cc_numubar_to_numubar"));
  components.push_back(std::make_tuple(nonswap, kIsNC,                    "nc"                   ));

  // CC appearance channels (for far detector only)
  if (sample.detector == covmx::kFarDet) {
    SpectrumLoaderBase* swap = &loaders.GetLoader(det,
      Loaders::kMC, ana::kBeam, Loaders::kFluxSwap);
    SpectrumLoaderBase* tau = &loaders.GetLoader(det,
      Loaders::kMC, ana::kBeam, Loaders::kTauSwap);
    components.push_back(std::make_tuple(swap, kIsSig       && !kIsAntiNu, "cc_numu_to_nue"        ));
    components.push_back(std::make_tuple(swap, kIsSig       && kIsAntiNu,  "cc_numubar_to_nuebar"  ));
    components.push_back(std::make_tuple(swap, kIsNumuApp   && !kIsAntiNu, "cc_nue_to_numu"        ));
    components.push_back(std::make_tuple(swap, kIsNumuApp   && kIsAntiNu,  "cc_nuebar_to_numubar"  ));
    components.push_back(std::make_tuple(tau,  kIsTauFromE  && !kIsAntiNu, "cc_nue_to_nutau"       ));
    components.push_back(std::make_tuple(tau,  kIsTauFromE  && kIsAntiNu,  "cc_nuebar_to_nutaubar" ));
    components.push_back(std::make_tuple(tau,  kIsTauFromMu && !kIsAntiNu, "cc_numu_to_nutau"      ));
    components.push_back(std::make_tuple(tau,  kIsTauFromMu && kIsAntiNu,  "cc_numubar_to_nutaubar"));
  }

  // PPFX weight vars
  std::vector<Var> kPPFXFluxUnivWgts;
  for (size_t 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);
  }

  // Analysis cuts and vars
  const HistAxis axis = sample.detector == covmx::kNearDet ? kNCNDAxisE : kNCFDAxisE;
  const Cut sel = sample.detector == covmx::kNearDet ? kNus18ND : kNus18FD;
  const Var kReweight = kXSecCVWgt2018RPAFix * kPPFXFluxCVWgt;

  // Initialise maps
  std::vector<Spectrum*> nominal(components.size());
  std::vector<std::vector<Spectrum*> > ppfx(components.size(),
    std::vector<Spectrum*>(kPPFXFluxUnivWgts.size(), nullptr));

  // Loop over each component
  for (size_t c = 0; c < components.size(); ++c) {

    // Nominal spectra
    nominal[c] = new Spectrum(*std::get<0>(components[c]),
      axis, sel && std::get<1>(components[c]), kNoShift, kReweight);

    // PPFX universe spectra
    for (size_t u = 0; u < 100; ++u)
      ppfx[c][u] = new Spectrum(*std::get<0>(components[c]),
        axis, sel && std::get<1>(components[c]), kNoShift, kXSecCVWgt2018RPAFix * kPPFXFluxUnivWgts[u]);

  } // for component

  loaders.Go();

  // Save em all
  std::string fileName = "ppfx_ratios_" + sample.GetTag() + ".root";
  TFile* f = TFile::Open(fileName.c_str(), "recreate");
  TDirectory* outDir = f->mkdir(sample.GetTag().c_str());

  TDirectory* dir = outDir->mkdir("nominal");
  for (size_t c = 0; c < components.size(); ++c) {
    nominal[c]->SaveTo(dir->mkdir(std::get<2>(components[c]).c_str()));
  }

  dir = outDir->mkdir("ppfx");
  for (size_t u = 0; u < 100; ++u) {
    TDirectory* udir = dir->mkdir(Form("u%zu",u+1));
    for (size_t c = 0; c < components.size(); ++c) {
      ppfx[c][u]->SaveTo(udir->mkdir(std::get<2>(components[c]).c_str()));
    }
  }

  delete f;

}