MakeSystRatios.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/Cuts/SpillCuts.h"
#include "NuXAna/Cuts/NusCuts18.h"
#include "CAFAna/Vars/HistAxes.h"
#include "NuXAna/Vars/HistAxes.h"
#include "CAFAna/Vars/XsecTunes.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Analysis/Prod4Loaders.h"

using namespace ana;

void MakeSystRatios(TString opt) {

  bool nd = opt.Contains("neardet");

  // Here get loaders
  Prod4NomLoaders loaders(ECAFType::kNusConcat, Loaders::FluxType::kFHC);
  if (nd) loaders.SetND(true);
  loaders.SetSpillCut(kStandardSpillCuts);

  // This should surely be a vector of loaders and cuts, right?

  auto det = nd ? 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 (!nd) {
    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);
  }

  // GENIE systs

  std::vector<const ISyst*> systs = getAllXsecSysts_2018();

  // Analysis cuts and vars

  const HistAxis axis = nd ? kNCNDAxisE : kNCFDAxisE;
  const Cut sel = nd ? kNus18ND : kNus18FD;
  const Var kReweight = kXSecCVWgt2018 * kPPFXFluxCVWgt;

  // Initialise maps

  std::vector<Spectrum*> nominal(components.size());
  std::vector<std::vector<Spectrum*> > ppfx(components.size(),
    std::vector<Spectrum*>(kPPFXFluxUnivWgts.size(), nullptr));
  std::vector<std::vector<Spectrum*> > genie_up(components.size(),
    std::vector<Spectrum*>(systs.size(), nullptr));
  std::vector<std::vector<Spectrum*> > genie_down(components.size(),
    std::vector<Spectrum*>(systs.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);

    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, kXSecCVWgt2018 * kPPFXFluxUnivWgts[u]);

    for (size_t s = 0; s < systs.size(); ++s) {

      SystShifts shift_up;
      shift_up.SetShift(systs[s], 1);
      genie_up[c][s] = new Spectrum(*std::get<0>(components[c]),
        axis, sel && std::get<1>(components[c]), shift_up, kReweight);
      SystShifts shift_down;
      shift_down.SetShift(systs[s], -1);
      genie_down[c][s] = new Spectrum(*std::get<0>(components[c]),
        axis, sel && std::get<1>(components[c]), shift_down, kReweight);

    } // for genie syst

  } // for component

  loaders.Go();

  // Save em all

  TFile* f = TFile::Open("SystRatios.root", "recreate");

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

  dir = f->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, std::get<2>(components[c]).c_str());
  }

  dir = f->mkdir("genie");
  for (size_t s = 0; s < systs.size(); ++s) {
    TDirectory* sdir = dir->mkdir(systs[s]->ShortName().c_str());
    TDirectory* updir = sdir->mkdir("up");
    TDirectory* downdir = sdir->mkdir("down");
    for (size_t c = 0; c < components.size(); ++c) {
      genie_up[c][s]->SaveTo(updir, std::get<2>(components[c]).c_str());
      genie_down[c][s]->SaveTo(downdir, std::get<2>(components[c]).c_str());
    }
  }

  // loop and save things here

  delete f;

}