Nus18SystsGENIELoad.C

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////
// GENIE 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/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/Systs/XSecSystLists.h"
#include "CAFAna/Systs/XSecSysts.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "NuXAna/Vars/HistAxes.h"
#include "CAFAna/Vars/PPFXWeights.h"

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

using namespace ana;

/// Get GENIE systematics labels
std::map<std::string, rwgt::ReweightLabel_t> GetGENIEShiftLabels();

/// Get cross-section systematics labels
std::map<std::string, const ISyst*> GetXSecShiftLabels();

void Nus18SystsGENIELoad(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;

  // Specific GENIE options
  bool genie_nd = true, genie_fd = true;
  bool genie_nova = true, genie_general = true;
  if (opt.Contains("nd", TString::kIgnoreCase) and
      !opt.Contains("fd", TString::kIgnoreCase))
    genie_fd = false;
  if (opt.Contains("fd", TString::kIgnoreCase) and
      !opt.Contains("nd", TString::kIgnoreCase))
    genie_nd = false;
  if (opt.Contains("NoPred", TString::kIgnoreCase)) {
    genie_nd = false;
    genie_fd = false;
  }
  if (opt.Contains("1", TString::kIgnoreCase) and
      !opt.Contains("2", TString::kIgnoreCase))
    genie_general = false;
  if (opt.Contains("2", TString::kIgnoreCase) and
      !opt.Contains("1", TString::kIgnoreCase))
    genie_nova = false;  

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

  // Vector of sigmas, levels to set each systematic
  std::vector<int> sigmas{-2, -1, +1, +2};

  // Get the labels for GENIE systematics
  std::map<std::string, rwgt::ReweightLabel_t> genieLabels;
  if (genie_nova)
    genieLabels = GetGENIEShiftLabels();

  // Get the labels for the cross-section systematics
  std::map<std::string, const ISyst*> xsecLabels;
  if (genie_general)
    xsecLabels = GetXSecShiftLabels();

  // Define map of systematic objects, indexed by same labels as previous map
  std::map<std::string, std::map<int, SystShifts*> > systematics;

  // Fill the map with all shifts
  for (const auto& sigma : sigmas) {
    for (const auto& genieLabel : genieLabels)
      systematics[genieLabel.first][sigma] = new SystShifts(GetGenieKnobSyst(genieLabel.second), sigma);
    for (const auto& xsecLabel : xsecLabels)
      systematics[xsecLabel.first][sigma] = new SystShifts(xsecLabel.second, sigma);
  }

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

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

    // Set up the shifted decompositions
    // One is needed for every sample, systematic, and sigma level
    for (const auto& systLabel : systematics) {
      for (const auto& syst : systLabel.second) {
        if (genie_nd)
          predsND_shifted[sample.first][systLabel.first][syst.first]
            = (NDPredictionSterile*)gennd.Generate(*loaders, *syst.second).release();
        if (genie_fd)
          predsFD_shifted[sample.first][systLabel.first][syst.first]
            = (FDPredictionSterile*)genfd.Generate(*loaders, *syst.second).release();
        if (extrap)
          predsExtrap_shifted[sample.first][systLabel.first][syst.first]
            = (PredictionSterile*)genExtrap->Generate(*loaders, *syst.second).release();
      }
    }

    if (extrap)
      delete genExtrap;

  }

  loaders->Go();

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

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

  // delete loaders;
  // for (const auto& syst : systematics)
  //   for (const auto& sigma : syst.second)
  //     delete sigma.second;

  return;

}

std::map<std::string, rwgt::ReweightLabel_t> GetGENIEShiftLabels() {

  std::map<std::string, rwgt::ReweightLabel_t> shiftlabels;

  // Ma NCEL, affects dsigma(NCEL)/dQ2 both in shape and normalization
  shiftlabels["ReweightMaNCEL"] = rwgt::fReweightMaNCEL;
  // NCEL strange axial form factor eta, affects dsigma(NCEL)/dQ2 both in shape and normalization
  shiftlabels["ReweightEtaNCEL"] = rwgt::fReweightEtaNCEL;

  // CCQE parameters -----------------------
  // elastic nucleon form factors (BBA/default -> dipole) - shape only effect of dsigma(CCQE)/dQ2
  shiftlabels["ReweightVecCCQEshape"] = rwgt::fReweightVecCCQEshape;
  // Ma CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
  shiftlabels["ReweightMaCCRES"] = rwgt::fReweightMaCCRES;
  // Mv CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
  shiftlabels["ReweightMvCCRES"] = rwgt::fReweightMvCCRES;
  // Ma NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
  shiftlabels["ReweightMaNCRES"] = rwgt::fReweightMaNCRES;
  // Mv NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
  shiftlabels["ReweightMvNCRES"] = rwgt::fReweightMvNCRES;

  // Coherent pion production parameters -----------------------
  // Ma for COH pion production
  shiftlabels["ReweightMaCOHpi"] = rwgt::fReweightMaCOHpi;
  // R0 for COH pion production
  shiftlabels["ReweightR0COHpi"] = rwgt::fReweightR0COHpi;
  // the Bodek-Yang model parameter A_{ht} - incl. both shape and normalization effect
  shiftlabels["ReweightAhtBY"] = rwgt::fReweightAhtBY;
  // the Bodek-Yang model parameter B_{ht} - incl. both shape and normalization effect
  shiftlabels["ReweightBhtBY"] = rwgt::fReweightBhtBY;
  // the Bodek-Yang model parameter CV1u - incl. both shape and normalization effect
  shiftlabels["ReweightCV1uBY"] = rwgt::fReweightCV1uBY;
  // the Bodek-Yang model parameter CV2u - incl. both shape and normalization effect
  shiftlabels["ReweightCV2uBY"] = rwgt::fReweightCV2uBY;
  // // the Bodek-Yang model parameter A_{ht} - shape only effect to d2sigma(DIS)/dxdy
  // shiftlabels["ReweightAhtBYshape"] = rwgt::fReweightAhtBYshape;
  // // the Bodek-Yang model parameter B_{ht} - shape only effect to d2sigma(DIS)/dxdy
  // shiftlabels["ReweightBhtBYshape"]  = rwgt::fReweightBhtBYshape;
  // // the Bodek-Yang model parameter CV1u - shape only effect to d2sigma(DIS)/dxdy
  // shiftlabels["ReweightCV1uBYshape"]  = rwgt::fReweightCV1uBYshape;
  // // the Bodek-Yang model parameter CV2u - shape only effect to d2sigma(DIS)/dxdy
  // shiftlabels["ReweightCV2uBYshape"]  = rwgt::fReweightCV2uBYshape;
  // the inclusive DIS CC normalization (not currently working in genie)
  shiftlabels["ReweightNormDISCC"] = rwgt::fReweightNormDISCC;
  // the ratio of \sigma(\bar\nu CC) / \sigma(\nu CC) (not currently working in genie)
  shiftlabels["ReweightRnubarnuCC"] = rwgt::fReweightRnubarnuCC;
  // DIS nuclear modification (shadowing, anti-shadowing, EMC). Not working in GENIE at the moment
  shiftlabels["ReweightDISNuclMod"] = rwgt::fReweightDISNuclMod;

  // NC: general -----------------------
  // NC total normalization (currently does nothing, but should retain for future)
  shiftlabels["ReweightNC"] = rwgt::fReweightNC;

  // Hadronization (free nucleon target) -----------------------
  // xF distribution for low multiplicity (N + pi) DIS f/s produced by AGKY
  shiftlabels["ReweightAGKY_xF1pi"] = rwgt::fReweightAGKY_xF1pi;
  // pT distribution for low multiplicity (N + pi) DIS f/s produced by AGKY
  shiftlabels["ReweightAGKY_pT1pi"] = rwgt::fReweightAGKY_pT1pi;
  // Hadron formation zone
  shiftlabels["ReweightFormZone"] = rwgt::fReweightFormZone;

  // Resonance decays -----------------------
  // Intranuclear rescattering systematics. There are 2 sets of parameters:
  // - parameters that control the total rescattering probability, P(total)
  // - parameters that control the fraction of each process (`fate'), given a total rescat. prob., P(fate|total)
  // These parameters are considered separately for pions and nucleons.
  // -----------------------
  // mean free path for pions
  shiftlabels["ReweightMFP_pi"] = rwgt::fReweightMFP_pi;
  // mean free path for nucleons
  shiftlabels["ReweightMFP_N"] = rwgt::fReweightMFP_N;
  // charge exchange probability for pions, for given total rescattering probability
  shiftlabels["ReweightFrCEx_pi"] = rwgt::fReweightFrCEx_pi;
  // elastic   probability for pions, for given total rescattering probability
  shiftlabels["ReweightFrElas_pi"] = rwgt::fReweightFrElas_pi;
  // inelastic probability for pions, for given total rescattering probability`
  shiftlabels["ReweightFrInel_pi"] = rwgt::fReweightFrInel_pi;
  // absorption probability for pions, for given total rescattering probability
  shiftlabels["ReweightFrAbs_pi"] = rwgt::fReweightFrAbs_pi;
  // pion production probability for pions, for given total rescattering probability
  shiftlabels["ReweightFrPiProd_pi"] = rwgt::fReweightFrPiProd_pi;
  // charge exchange probability for nucleons, for given total rescattering probability
  shiftlabels["ReweightFrCEx_N"] = rwgt::fReweightFrCEx_N;
  // elastic    probability for nucleons, for given total rescattering probability
  shiftlabels["ReweightFrElas_N"] = rwgt::fReweightFrElas_N;
  // inelastic  probability for nucleons, for given total rescattering probability
  shiftlabels["ReweightFrInel_N"] = rwgt::fReweightFrInel_N;
  // absorption probability for nucleons, for given total rescattering probability
  shiftlabels["ReweightFrAbs_N"] = rwgt::fReweightFrAbs_N;
  // pion production probability for nucleons, for given total rescattering probability
  shiftlabels["ReweightFrPiProd_N"] = rwgt::fReweightFrPiProd_N;
  shiftlabels["ReweightCCQEPauliSupViaKF"] = rwgt::fReweightCCQEPauliSupViaKF;
  shiftlabels["ReweightCCQEMomDistroFGtoSF"] = rwgt::fReweightCCQEMomDistroFGtoSF;
  // Resonance -> X + gamma branching ratio, eg Delta+(1232) -> p gamma
  shiftlabels["ReweightBR1gamma"] = rwgt::fReweightBR1gamma;
  // Resonance -> X + eta   branching ratio, eg N+(1440) -> p eta
  shiftlabels["ReweightBR1eta"] = rwgt::fReweightBR1eta;
  // distort pi angular distribution in Delta -> N + pi
  shiftlabels["ReweightTheta_Delta2Npi"] = rwgt::fReweightTheta_Delta2Npi;

  // Nuclear model -----------------------
  // shiftlabels["ReweightNC"] = rwgt::fReweightNC;

  return shiftlabels;
  
}

std::map<std::string, const ISyst*> GetXSecShiftLabels() {

  std::map<std::string, const ISyst*> shiftlabels;

  // Get cross-section systematics
  std::vector<const ISyst*> AllSyst= getAllXsecSysts_2018();

  shiftlabels["PD3_MAQAGenieReducedSyst2018"]         = AllSyst.at(0);
  shiftlabels["PD3_RPACCQEEnhSyst2018"]               = AllSyst.at(1);
  shiftlabels["PD3_RPACCQESuppSyst2018"]              = AllSyst.at(2);
  shiftlabels["PD3_RPARESSyst2018"]                   = AllSyst.at(3);
  shiftlabels["PD3_MECEShapeSyst2018Nu"]              = AllSyst.at(4);
  // shiftlabels["PD3_MECq0ShapeSyst2018"]               = AllSyst.at(4);
  shiftlabels["PD3_MECEShapeSyst2018AntiNu"]          = AllSyst.at(5);
  shiftlabels["PD3_MECEnuShapeSyst2018Nu"]            = AllSyst.at(6);
  // shiftlabels["PD3_MECEnuShapeSyst2018"]              = AllSyst.at(6);
  shiftlabels["PD3_MECEnuShapeSyst2018AntiNu"]        = AllSyst.at(7);
  shiftlabels["PD3_MECInitStateNPFracSyst2018Nu"]     = AllSyst.at(8);
  shiftlabels["PD3_MECInitStateNPFracSyst2018AntiNu"] = AllSyst.at(9);
  shiftlabels["PD3_DISvpCC0pi"]                       = AllSyst.at(10);
  shiftlabels["PD3_DISvpCC1pi"]                       = AllSyst.at(11);
  shiftlabels["PD3_DISvpCC2pi"]                       = AllSyst.at(12);
  shiftlabels["PD3_DISvpCC3pi"]                       = AllSyst.at(13);
  shiftlabels["PD3_DISvpNC0pi"]                       = AllSyst.at(14);
  shiftlabels["PD3_DISvpNC1pi"]                       = AllSyst.at(15);
  shiftlabels["PD3_DISvpNC2pi"]                       = AllSyst.at(16);
  shiftlabels["PD3_DISvpNC3pi"]                       = AllSyst.at(17);
  shiftlabels["PD3_DISvnCC0pi"]                       = AllSyst.at(18);
  shiftlabels["PD3_DISvnCC1pi"]                       = AllSyst.at(19);
  shiftlabels["PD3_DISvnCC2pi"]                       = AllSyst.at(20);
  shiftlabels["PD3_DISvnCC3pi"]                       = AllSyst.at(21);
  shiftlabels["PD3_DISvnNC0pi"]                       = AllSyst.at(22);
  shiftlabels["PD3_DISvnNC1pi"]                       = AllSyst.at(23);
  shiftlabels["PD3_DISvnNC2pi"]                       = AllSyst.at(24);
  shiftlabels["PD3_DISvnNC3pi"]                       = AllSyst.at(25);
  shiftlabels["PD3_DISvbarpCC0pi"]                    = AllSyst.at(26);
  shiftlabels["PD3_DISvbarpCC1pi"]                    = AllSyst.at(27);
  shiftlabels["PD3_DISvbarpCC2pi"]                    = AllSyst.at(28);
  shiftlabels["PD3_DISvbarpCC3pi"]                    = AllSyst.at(29);
  shiftlabels["PD3_DISvbarpNC0pi"]                    = AllSyst.at(30);
  shiftlabels["PD3_DISvbarpNC1pi"]                    = AllSyst.at(31);
  shiftlabels["PD3_DISvbarpNC2pi"]                    = AllSyst.at(32);
  shiftlabels["PD3_DISvbarpNC3pi"]                    = AllSyst.at(33);
  shiftlabels["PD3_DISvbarnCC0pi"]                    = AllSyst.at(34);
  shiftlabels["PD3_DISvbarnCC1pi"]                    = AllSyst.at(35);
  shiftlabels["PD3_DISvbarnCC2pi"]                    = AllSyst.at(36);
  shiftlabels["PD3_DISvbarnCC3pi"]                    = AllSyst.at(37);
  shiftlabels["PD3_DISvbarnNC0pi"]                    = AllSyst.at(38);
  shiftlabels["PD3_DISvbarnNC1pi"]                    = AllSyst.at(39);
  shiftlabels["PD3_DISvbarnNC2pi"]                    = AllSyst.at(40);
  shiftlabels["PD3_DISvbarnNC3pi"]                    = AllSyst.at(41);
  shiftlabels["PD3_RadCorrNue"]                       = AllSyst.at(42);
  shiftlabels["PD3_RadCorrNuebar"]                    = AllSyst.at(43);
  shiftlabels["PD3_k2ndClassCurr"]                    = AllSyst.at(44);

  return shiftlabels;

}