SystsDecompLoad17.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/NCDecomp.h"
#include "3FlavorAna/Decomp/NueDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.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"
#include "NuXAna/Cuts/NusCuts17.h"
#include "CAFAna/Vars/PPFXWeights.h"

using namespace ana;

IDecomp* MakeDecomp(Loaders* loaders, const HistAxis& axis,
                    const Cut& cut, std::string shiftlabel,
                    const SystShifts& kMCShift = kNoShift,
                    const Var& wei = kUnweighted);
  const Var kWgtP31 = kXSecCVWgt2017*kPPFXFluxCVWgt;
void SystsDecompLoad17()
{
  TH1::AddDirectory(0);

  Loaders loaders;


  std::string fnear="defname: prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
  std::string ffar ="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1";
  std::string fswap="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1";
  std::string ftau="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_v1";
   std::string fnameneardata_con="defname: prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns";

  loaders.SetLoaderPath(fnear, caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(ffar,  caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(fswap, caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kFluxSwap);
  loaders.SetLoaderPath(ftau,  caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kTauSwap);
  loaders.SetLoaderPath(fnameneardata_con, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);

  loaders.SetSpillCut(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["NC"]      = "NC";
  shiftlabels["Numu"]    = "Numu";
  shiftlabels["Nue"]   = "Nue";

  // A map of systematic objects, indexed by the same label as above
  std::map<std::string, std::vector<IDecomp*> > systematics;

  // 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 = "Visible Energy (GeV)";
  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
  std::pair<Cut*, Cut*> Ana01(new Cut(kNus17ND), new Cut(kNus17FDAlt));
  cut_samples["Ana01"] = std::pair<HistAxis*, std::pair<Cut*, Cut*> >(
    new HistAxis(labelRecoE, kNCDisappearanceEnergyBinning, kNusNDEnergy17),
    Ana01
  );

   HistAxis* FDAxis = new HistAxis(labelRecoE, kNCDisappearanceEnergyBinning, kNusFDEnergy17);
  // 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(
      *FDAxis, kNCBinsNumuCCAxis,
      *sample.second.second.second, *sample.second.second.first,
      kNumuND, kNoShift, kWgtP31
    );

    // Set up the nominal decompositions
    // Only one is needed per sample
    decomps_nominal[sample.first] = MakeDecomp(
      &loaders, *sample.second.first, *sample.second.second.first, "",
      kNoShift, kWgtP31
    );
    predsNE_nominal[sample.first] = new PredictionNoExtrap(
      loaders,
      *FDAxis, *sample.second.second.second,
      kNoShift, kWgtP31
    );
    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] = MakeDecomp(
        &loaders, *sample.second.first, *sample.second.second.first, shiftlabel.first,
        kNoShift, kWgtP31
      );
      predsNE_shifted[sample.first][shiftlabel.first][1] = new PredictionNoExtrap(
        loaders, *FDAxis, *sample.second.second.second,
        kNoShift, kWgtP31
      );

      IDecomp* decompnc = MakeDecomp(
        &loaders, *sample.second.first, *sample.second.second.first, shiftlabel.first,
        kNoShift, kWgtP31
      );
      IDecomp* decompnm = MakeDecomp(
        &loaders, kNCBinsNumuCCAxis, kNumuND, shiftlabel.first,
        kNoShift, kWgtP31
      );
      ModularExtrapSterile* extrap = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
        loaders, *decompnc, *decompnm,
        *FDAxis, kNCBinsNumuCCAxis,
        *sample.second.second.second, *sample.second.second.first, kNumuND,
        kNoShift, kWgtP31
      ));
      predsSt_shifted[sample.first][shiftlabel.first][1] =
        new PredictionSterile(extrap);
    }
  }

  loaders.Go();

  std::string folder = "./";
  std::string filenm = "SystsDecomp";

  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
  );
}

IDecomp* MakeDecomp(Loaders* loaders, const HistAxis& axis,
                    const Cut& cut, std::string shiftlabel,
                    const SystShifts& kMCShift,
                    const Var& wei)
{
  if(shiftlabel.compare("NC") == 0) {
    return new NCDecomp(
      loaders->GetLoader(caf::kNEARDET, Loaders::kMC,   ana::kBeam, Loaders::kNonSwap),
      loaders->GetLoader(caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap),
      axis, cut, kMCShift, kNoShift, wei
    );
  }
  if(shiftlabel.compare("Numu") == 0) {
    return new NumuDecomp(
      loaders->GetLoader(caf::kNEARDET, Loaders::kMC,   ana::kBeam, Loaders::kNonSwap),
      loaders->GetLoader(caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap),
      axis, cut, kMCShift, kNoShift, wei
    );
  }
  if(shiftlabel.compare("Nue") == 0) {
    return new NueDecomp(
      loaders->GetLoader(caf::kNEARDET, Loaders::kMC,   ana::kBeam, Loaders::kNonSwap),
      loaders->GetLoader(caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap),
      axis, cut, kMCShift, kNoShift, wei
    );
  }
  std::cout << "Label not found, reverting to ProportionalDecomp." << std::endl;
  return new ProportionalDecomp(
    loaders->GetLoader(caf::kNEARDET, Loaders::kMC,   ana::kBeam, Loaders::kNonSwap),
    loaders->GetLoader(caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap),
    axis, cut, kMCShift, kNoShift, wei
  );
}