Nus18SystsKaonLoad.C

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////
// Kaon systematics for Nus18
//
// Author: Sijith Edayath
// Modified by Mike Wallbank (wallbank@fnal.gov)
/////////////////////////////////////////////////////////////////////

#include "CAFAna/Analysis/Prod4Loaders.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Cuts/Cuts.h"
#include "NuXAna/Cuts/NusCuts18.h"
#include "3FlavorAna/Cuts/NumuCuts2018.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/NDPredictionSterile.h"
#include "NuXAna/Prediction/FDPredictionSterile.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "CAFAna/Systs/BeamSysts.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "NuXAna/Vars/HistAxes.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "NuXAna/macros/Nus18/Nus18Systs/HelperSyst18/NuSSystFunctions18.h"

using namespace ana;

/// Run systematics for Kaons
void Nus18SystsKaonLoad(TString opt) {

  TH1::AddDirectory(0);
  const Var kReweight = kXSecCVWgt2018*kPPFXFluxCVWgt;

  // Get options
  ECAFType cafType;
  Loaders::FluxType fluxType;
  bool extrap = false;
  if (opt.Contains("caf", TString::kIgnoreCase))
    cafType = ECAFType::kFullCAF;
  else if (opt.Contains("concat", TString::kIgnoreCase))
    cafType = ECAFType::kNusConcat;
  else {
    std::cout << "ECAFType required: caf or concat" << std::endl;
    abort();
  }
  if (opt.Contains("fhc", TString::kIgnoreCase))
    fluxType = Loaders::FluxType::kFHC;
  else if (opt.Contains("rhc", TString::kIgnoreCase))
    fluxType = Loaders::FluxType::kRHC;
  else {
    std::cout << "FluxType required: fhc or rhc" << std::endl;
    abort();
  }
  if (opt.Contains("extrap", TString::kIgnoreCase))
    extrap = true;

  // Set up the loaders
  Loaders* loaders = new Prod4NomLoaders(cafType, fluxType);
  loaders->SetSpillCut(kStandardSpillCuts);

  // Get the labels for Kaon systematics
  std::map<std::string, std::string> shiftLabels;
  shiftLabels["Kaon"] = "Kaon";
  const Var kKaonP = kReweight*kNus18BeamWeightP;
  const Var kKaonM = kReweight*kNus18BeamWeightM;

  // Define map of systematic objects, indexed by same labels as previous map
  std::map<std::string, std::map<int, const Var*> > systematics;
  systematics["Kaon"][-1] = &kKaonM;
  systematics["Kaon"][+1] = &kKaonP;

  // Define a map of samples and cuts/axes to set them up
  // Each sample, essentially, is a way to set up the analysis
  HistAxis AxisRecoE = kNus18AxisE;
  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
  cut_samples["Nus18"] = std::make_pair(&AxisRecoE, std::make_pair(new Cut(kNus18ND), new Cut(kNus18FD)));

  // 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, NDPredictionSterile*> predsND_nominal;
  std::map<std::string, FDPredictionSterile*> predsFD_nominal;
  std::map<std::string, PredictionSterile*>   predsExtrap_nominal;
  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > > predsND_shifted;
  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > > predsFD_shifted;
  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >   predsExtrap_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) {

    // const NDPredictionGenerator gennd(*sample.second.first, *sample.second.second.first, kNoShift, kReweight);
    // const FDPredictionGenerator genfd(*sample.second.first, *sample.second.second.second, kNoShift, kReweight);
    SterileGenerator* genExtrap;
    // if (extrap) {
      genExtrap = new SterileGenerator(*sample.second.first, kNus18BinsNumuCCAxis,
                                       *sample.second.second.second, *sample.second.second.first,    
                                       kNumuCutND2018, kNoShift, kReweight);
      // ProportionalDecomp* decompNC_nominal = new ProportionalDecomp
      //        (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
      //         loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kData),
      //         *sample.second.first, *sample.second.second.first,
      //         kNoShift, kNoShift, kReweight);
      // ProportionalDecomp* decompNuMu_nominal = new ProportionalDecomp
      //        (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
      //         loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kData),
      //         kNus18BinsNumuCCAxis, kNumuCutND2018,
      //         kNoShift, kNoShift, kReweight);
      // ModularExtrapSterile* extrap_nominal = new ModularExtrapSterile
      //        (ModularExtrapSterile::NCDisappearance
      //         (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
      //          loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kFluxSwap),
      //          loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kNonSwap),
      //          loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kTauSwap),
      //          *decompNC_nominal, *decompNuMu_nominal,
      //          *sample.second.first, kNus18BinsNumuCCAxis,
      //          *sample.second.second.second, *sample.second.second.first, kNumuCutND2018,
      //          kNoShift, kReweight));
      // predsExtrap_nominal[sample.first] = new PredictionSterile(extrap_nominal);
    // }

    // Set up the nominal decompositions
    // Only one is needed per sample
    // predsND_nominal[sample.first] = (NDPredictionSterile*)gennd.Generate(*loaders).release();
    // predsFD_nominal[sample.first] = (FDPredictionSterile*)genfd.Generate(*loaders).release();
    // if (extrap) {
      predsExtrap_nominal[sample.first] = (PredictionSterile*)genExtrap->Generate(*loaders).release();
    // }

    // if (extrap)
    //   delete genExtrap;

    // Set up the shifted decompositions
    // One is needed for every sample, systematic, and sigma level
    for(const auto& shiftLabel : shiftLabels) {
      for(const auto& syst : systematics[shiftLabel.first]) {

        // // Create new prediction generators with correct reweights
        // const NDPredictionGenerator gennd_kaon
        //   (*sample.second.first, *sample.second.second.first, kNoShift, *syst.second);
        // const FDPredictionGenerator genfd_kaon
        //   (*sample.second.first, *sample.second.second.second, kNoShift, *syst.second);
        SterileGenerator* genExtrap_kaon;
        // if (extrap) {
          genExtrap_kaon = new SterileGenerator(*sample.second.first, kNus18BinsNumuCCAxis,
                                                *sample.second.second.second, *sample.second.second.first,
                                                kNumuCutND2018, kNoShift, *syst.second);
          // ProportionalDecomp* decompNC_shifted = new ProportionalDecomp
          //   (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
          //    loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kData),
          //    *sample.second.first, *sample.second.second.first,
          //    kNoShift, kNoShift, *syst.second);
          // ProportionalDecomp* decompNuMu_shifted = new ProportionalDecomp
          //   (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
          //    loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kData),
          //    kNus18BinsNumuCCAxis, kNumuCutND2018,
          //    kNoShift, kNoShift, *syst.second);
          // ModularExtrapSterile* extrap_shifted = new ModularExtrapSterile
          //   (ModularExtrapSterile::NCDisappearance
          //    (loaders->GetLoader(caf::Det_t::kNEARDET, Loaders::DataMC::kMC),
          //     loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kFluxSwap),
          //     loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kNonSwap),
          //     loaders->GetLoader(caf::Det_t::kFARDET, Loaders::DataMC::kMC, DataSource::kBeam, Loaders::SwappingConfig::kTauSwap),
          //     *decompNC_shifted, *decompNuMu_shifted,
          //     *sample.second.first, kNus18BinsNumuCCAxis,
          //     *sample.second.second.second, *sample.second.second.first, kNumuCutND2018,
          //     kNoShift, *syst.second));
          // predsExtrap_shifted[sample.first][shiftLabel.first][syst.first] = new PredictionSterile(extrap_shifted);
        // }

        // Shift the spectra
        // predsND_shifted[sample.first][shiftLabel.first][syst.first]
        //   = (NDPredictionSterile*)gennd_kaon.Generate(*loaders).release();
        // predsFD_shifted[sample.first][shiftLabel.first][syst.first]
        //   = (FDPredictionSterile*)genfd_kaon.Generate(*loaders).release();
        // if (extrap) {
          predsExtrap_shifted[sample.first][shiftLabel.first][syst.first]
            = (PredictionSterile*)genExtrap_kaon->Generate(*loaders).release();
        // }

        // if (extrap)
        //   delete genExtrap_kaon;

      }
    }

  }

  loaders->Go();

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

  std::string fullLocation = folder + filenm + ".root";
  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");
  // SaveMaps(rootF, predsND_nominal, predsND_shifted);
  // SaveMaps(rootF, predsFD_nominal, predsFD_shifted);
  // if (extrap)
    SaveMaps(rootF, predsExtrap_nominal, predsExtrap_shifted);
  rootF->Close();

  // delete loaders;

  return;

}