SystsNDRockLoad.C

Go to the documentation of this file.

#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Cuts/Cuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "CAFAna/Extrap/ExtrapSterile.h"
#include "CAFAna/Prediction/PredictionGenerator.h"
#include "NuXAna/Prediction/PredictionGeneratorNuX.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "NuXAna/Vars/HistAxes.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "NuXAna/macros/NuSSystFunctions.h"

using namespace ana;

void SystsNDRockLoad()
{
  TH1::AddDirectory(0);

  Loaders loaders;

  loaders.SetLoaderPath(fnamenear_concat, caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(fnamefar_concat,  caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(fnameswap_concat, caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kFluxSwap);
  loaders.SetLoaderPath(fnametau_concat,  caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kTauSwap);
  loaders.SetLoaderPath(fnameneardata_concat, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);

  loaders.SetSpillCut(kOnly14DB && kStandardSpillCuts);

  // A map to set up systematics, starting with a label for the systematic,
  // and finishing with a vector for items necessary to set up the Syst object
  std::map<std::string, std::string> shiftlabels;
  shiftlabels["NDRock"] = "ND Rock Contamination";

  const Cut kNDTruthFull(
    [](const caf::SRProxy* sr)
    {
      if(sr->hdr.ismc == false) { return true; }
      if(sr->mc.nnu < 1) { return false; }
      if(sr->mc.nu[0].vtx.X() < NDL) { return false; }
      if(sr->mc.nu[0].vtx.X() > NDR) { return false; }
      if(sr->mc.nu[0].vtx.Y() < NDB) { return false; }
      if(sr->mc.nu[0].vtx.Y() > NDT) { return false; }
      if(sr->mc.nu[0].vtx.Z() < NDF) { return false; }

      double NDM = 1587.;
      if(sr->mc.nu[0].vtx.Z() > NDM) { return false; }
      if(sr->mc.nu[0].vtx.Z() > NDE) {
        double MCT = (2./3.)*(NDT - NDB) + NDB;
        if(sr->mc.nu[0].vtx.Y() > MCT) { return false; }
      }
      return true;
    });

  // A map of systematic objects, indexed by the same label as above
  std::map<std::string, Cut*> systematics;
  systematics["NDRock"] = new Cut(kNDTruthFull);

  // Define a map of samples and cuts/axes to set them up
  // Each sample, essentially, is a way to set up the analysis
  std::string labelRecoE = "Calorimetric Energy (GeV)";
  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
  std::pair<Cut*, Cut*> Ana01(new Cut(kNusND), new Cut(kNusFD));
  cut_samples["Ana01"] = std::pair<HistAxis*, std::pair<Cut*, Cut*> >(
    new HistAxis(labelRecoE, kNCDisappearanceEnergyBinning, kCaloE),
    Ana01
  );

  // Set up maps to the decompositions/predictions (each flavor component)
  // Nominal maps are indexed purely by sample label
  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
  std::map<std::string,                                     IDecomp*>     decomps_nominal;
  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
  std::map<std::string,                                     PredictionNoExtrap*>     predsNE_nominal;
  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
  std::map<std::string,                                     PredictionSterile*>      predsSt_nominal;
  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >  predsSt_shifted;

  // Set up the actual decompositions/predictions
  // sample.first = the sample label
  // shiftlabel.first = the systematic shift label
  // sample.second.first = the sample HistAxis
  // sample.second.second = the sample Cut
  // syst.first = the integer number of sigmas for the shift
  // syst.second = the systematic being shifted
  for(const auto& sample : cut_samples) {
    SterileGenerator gen(
      *sample.second.first, kNCBinsNumuCCAxis,
      *sample.second.second.second, *sample.second.second.first,
      kNumuND, kNoShift, kTuftsWeightCC
    );

    // Set up the nominal decompositions
    // Only one is needed per sample
    decomps_nominal[sample.first] = new CheatDecomp(
      loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
      *sample.second.first, *sample.second.second.first,
      kNoShift, kTuftsWeightCC
    );
    predsNE_nominal[sample.first] = new PredictionNoExtrap(
      loaders,
      *sample.second.first, *sample.second.second.second,
      kNoShift, kTuftsWeightCC
    );
    predsSt_nominal[sample.first] = (PredictionSterile*)gen.Generate(
      loaders
    ).release();

    // Set up the shifted decompositions
    // One is needed for every sample, systematic, and sigma level
    for(const auto& shiftlabel : shiftlabels) {
      decomps_shifted[sample.first][shiftlabel.first][1] = new CheatDecomp(
        loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
        *sample.second.first, *sample.second.second.first && *systematics[shiftlabel.first],
        kNoShift, kTuftsWeightCC
      );
      predsNE_shifted[sample.first][shiftlabel.first][1] = new PredictionNoExtrap(
        loaders,
        *sample.second.first, *sample.second.second.second,
        kNoShift, kTuftsWeightCC
      );

      SterileGenerator gen(
        *sample.second.first, kNCBinsNumuCCAxis,
        *sample.second.second.second,
        *sample.second.second.first && *systematics[shiftlabel.first],
        kNumuND && *systematics[shiftlabel.first],
        kNoShift, kTuftsWeightCC
      );
      predsSt_shifted[sample.first][shiftlabel.first][1] =
        (PredictionSterile*)gen.Generate(loaders).release();
    }
  }

  loaders.Go();

  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
  std::string filenm = "SystsNDRock";

  std::string fullLocation = folder + filenm + ".root";
  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");

  SaveMaps(rootF,
    decomps_nominal, decomps_shifted,
    predsNE_nominal, predsNE_shifted,
    predsSt_nominal, predsSt_shifted
  );
}