SystsGENIELoad.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/Extrap/ExtrapSterile.h"
#include "CAFAna/Prediction/PredictionGenerator.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 SystsGENIELoad()
{
  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);

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

/*
redundant
"ReweightMaCCQE"
"ReweightMaCCRESshape"
"ReweightMvCCRESshape"
"ReweightMaNCRESshape"
"ReweightMvNCRESshape"
"ReweightAhtBYshape"
"ReweightBhtBYshape"
"ReweightCV1uBYshape"
"ReweightCV2uBYshape"
"ReweightNormCCQEenu"
"ReweightNormCCRES"
"ReweightNormNCRES"
"ReweightNC"
*/

  // 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, rwgt::ReweightLabel_t> shiftlabels;
  // NCEL parameters:
  shiftlabels["ReweightMaNCEL"]              = rwgt::fReweightMaNCEL;              // Ma NCEL, affects dsigma(NCEL)/dQ2 both in shape and normalization
  shiftlabels["ReweightEtaNCEL"]             = rwgt::fReweightEtaNCEL;             // NCEL strange axial form factor eta, affects dsigma(NCEL)/dQ2 both in shape and normalization
  // CCQE parameters:
  //shiftlabels["ReweightMaCCQE"]              = rwgt::fReweightMaCCQE;            // Ma CCQE, affects dsigma(CCQE)/dQ2 both in shape and normalization
  shiftlabels["ReweightMaCCQEshape"]         = rwgt::fReweightMaCCQEshape;         // Ma CCQE, affects dsigma(CCQE)/dQ2 in shape only (normalized to constant integral)
  shiftlabels["ReweightNormCCQE"]            = rwgt::fReweightNormCCQE;            // CCQE normalization (energy independent)
  //shiftlabels["ReweightNormCCQEenu"]         = rwgt::fReweightNormCCQEenu;       // CCQE normalization (maintains dependence on neutrino energy)
  shiftlabels["ReweightCCQEPauliSupViaKF"]   = rwgt::fReweightCCQEPauliSupViaKF;   //
  shiftlabels["ReweightVecCCQEshape"]        = rwgt::fReweightVecCCQEshape;        // elastic nucleon form factors (BBA/default -> dipole) - shape only effect of dsigma(CCQE)/dQ2
  shiftlabels["ReweightCCQEMomDistroFGtoSF"] = rwgt::fReweightCCQEMomDistroFGtoSF; //
  // Resonance neutrino-production parameters:
  //shiftlabels["ReweightNormCCRES"]           = rwgt::fReweightNormCCRES;    // CCRES normalization
  //shiftlabels["ReweightNormNCRES"]           = rwgt::fReweightNormNCRES;    // NCRES normalization
  //shiftlabels["ReweightMaCCRESshape"]        = rwgt::fReweightMaCCRESshape; // Ma CCRES, affects d2sigma(CCRES)/dWdQ2 in shape only (normalized to constant integral)
  //shiftlabels["ReweightMvCCRESshape"]        = rwgt::fReweightMvCCRESshape; // Mv CCRES, affects d2sigma(CCRES)/dWdQ2 in shape only (normalized to constant integral)
  //shiftlabels["ReweightMaNCRESshape"]        = rwgt::fReweightMaNCRESshape; // Ma NCRES, affects d2sigma(NCRES)/dWdQ2 in shape only (normalized to constant integral)
  //shiftlabels["ReweightMvNCRESshape"]        = rwgt::fReweightMvNCRESshape; // Mv NCRES, affects d2sigma(NCRES)/dWdQ2 in shape only (normalized to constant integral)
  shiftlabels["ReweightMaCCRES"]             = rwgt::fReweightMaCCRES;      // Ma CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
  shiftlabels["ReweightMvCCRES"]             = rwgt::fReweightMvCCRES;      // Mv CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
  shiftlabels["ReweightMaNCRES"]             = rwgt::fReweightMaNCRES;      // Ma NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
  shiftlabels["ReweightMvNCRES"]             = rwgt::fReweightMvNCRES;      // Mv NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
  // Coherent pion production parameters:
  shiftlabels["ReweightMaCOHpi"]             = rwgt::fReweightMaCOHpi; // Ma for COH pion production
  shiftlabels["ReweightR0COHpi"]             = rwgt::fReweightR0COHpi; // R0 for COH pion production
  // Non-resonance background parameters:
  shiftlabels["ReweightRvpCC1pi"]            = rwgt::fReweightRvpCC1pi;    // the 1pi non-RES bkg in the RES region, for v+p CC
  shiftlabels["ReweightRvpCC2pi"]            = rwgt::fReweightRvpCC2pi;    // the 2pi non-RES bkg in the RES region, for v+p CC
  shiftlabels["ReweightRvnCC1pi"]            = rwgt::fReweightRvnCC1pi;    // the 1pi non-RES bkg in the RES region, for v+n CC
  shiftlabels["ReweightRvnCC2pi"]            = rwgt::fReweightRvnCC2pi;    // the 2pi non-RES bkg in the RES region, for v+n CC
  shiftlabels["ReweightRvpNC1pi"]            = rwgt::fReweightRvpNC1pi;    // the 1pi non-RES bkg in the RES region, for v+p NC
  shiftlabels["ReweightRvpNC2pi"]            = rwgt::fReweightRvpNC2pi;    // the 2pi non-RES bkg in the RES region, for v+p NC
  shiftlabels["ReweightRvnNC1pi"]            = rwgt::fReweightRvnNC1pi;    // the 1pi non-RES bkg in the RES region, for v+n NC
  shiftlabels["ReweightRvnNC2pi"]            = rwgt::fReweightRvnNC2pi;    // the 2pi non-RES bkg in the RES region, for v+n NC
  shiftlabels["ReweightRvbarpCC1pi"]         = rwgt::fReweightRvbarpCC1pi; // the 1pi non-RES bkg in the RES region, for vbar+p CC
  shiftlabels["ReweightRvbarpCC2pi"]         = rwgt::fReweightRvbarpCC2pi; // the 2pi non-RES bkg in the RES region, for vbar+p CC
  shiftlabels["ReweightRvbarnCC1pi"]         = rwgt::fReweightRvbarnCC1pi; // the 1pi non-RES bkg in the RES region, for vbar+n CC
  shiftlabels["ReweightRvbarnCC2pi"]         = rwgt::fReweightRvbarnCC2pi; // the 2pi non-RES bkg in the RES region, for vbar+n CC
  shiftlabels["ReweightRvbarpNC1pi"]         = rwgt::fReweightRvbarpNC1pi; // the 1pi non-RES bkg in the RES region, for vbar+p NC
  shiftlabels["ReweightRvbarpNC2pi"]         = rwgt::fReweightRvbarpNC2pi; // the 2pi non-RES bkg in the RES region, for vbar+p NC
  shiftlabels["ReweightRvbarnNC1pi"]         = rwgt::fReweightRvbarnNC1pi; // the 1pi non-RES bkg in the RES region, for vbar+n NC
  shiftlabels["ReweightRvbarnNC2pi"]         = rwgt::fReweightRvbarnNC2pi; // the 2pi non-RES bkg in the RES region, for vbar+n NC
  // DIS parameters - applied for DIS events with (Q2>Q2o, W>Wo),
  // typically Q2okReweight =1GeV^2, WokReweight =1.7-2.0GeV
  shiftlabels["ReweightAhtBY"]               = rwgt::fReweightAhtBY;       // the Bodek-Yang model parameter A_{ht} - incl. both shape and normalization effect
  shiftlabels["ReweightBhtBY"]               = rwgt::fReweightBhtBY;       // the Bodek-Yang model parameter B_{ht} - incl. both shape and normalization effect
  shiftlabels["ReweightCV1uBY"]              = rwgt::fReweightCV1uBY;      // the Bodek-Yang model parameter CV1u - incl. both shape and normalization effect
  shiftlabels["ReweightCV2uBY"]              = rwgt::fReweightCV2uBY;      // the Bodek-Yang model parameter CV2u - incl. both shape and normalization effect
  //shiftlabels["ReweightAhtBYshape"]          = rwgt::fReweightAhtBYshape;  // the Bodek-Yang model parameter A_{ht} - shape only effect to d2sigma(DIS)/dxdy
  //shiftlabels["ReweightBhtBYshape"]          = rwgt::fReweightBhtBYshape;  // the Bodek-Yang model parameter B_{ht} - shape only effect to d2sigma(DIS)/dxdy
  //shiftlabels["ReweightCV1uBYshape"]         = rwgt::fReweightCV1uBYshape; // the Bodek-Yang model parameter CV1u - shape only effect to d2sigma(DIS)/dxdy
  //shiftlabels["ReweightCV2uBYshape"]         = rwgt::fReweightCV2uBYshape; // the Bodek-Yang model parameter CV2u - shape only effect to d2sigma(DIS)/dxdy
  shiftlabels["ReweightNormDISCC"]           = rwgt::fReweightNormDISCC;   // the inclusive DIS CC normalization (not currently working in genie)
  shiftlabels["ReweightRnubarnuCC"]          = rwgt::fReweightRnubarnuCC;  // the ratio of \sigma(\bar\nu CC) / \sigma(\nu CC) (not currently working in genie)
  shiftlabels["ReweightDISNuclMod"]          = rwgt::fReweightDISNuclMod;  // DIS nuclear modification (shadowing, anti-shadowing, EMC). Not working in GENIE at the moment

  // Hadronization (free nucleon target)
  shiftlabels["ReweightAGKY_xF1pi"]          = rwgt::fReweightAGKY_xF1pi; // xF distribution for low multiplicity (N + pi) DIS f/s produced by AGKY
  shiftlabels["ReweightAGKY_pT1pi"]          = rwgt::fReweightAGKY_pT1pi; // pT distribution for low multiplicity (N + pi) DIS f/s produced by AGKY
  shiftlabels["ReweightFormZone"]            = rwgt::fReweightFormZone;   // Hadron formation zone
  // Resonance decays
  shiftlabels["ReweightTheta_Delta2Npi"]     = rwgt::fReweightTheta_Delta2Npi; // distort pi angular distribution in Delta -> N + pi
  shiftlabels["ReweightBR1gamma"]            = rwgt::fReweightBR1gamma;        // Resonance -> X + gamma branching ratio, eg Delta+(1232) -> p gamma
  shiftlabels["ReweightBR1eta"]              = rwgt::fReweightBR1eta;          // Resonance -> X + eta   branching ratio, eg N+(1440) -> p eta

  // 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.
  shiftlabels["ReweightMFP_N"]               = rwgt::fReweightMFP_N;       // mean free path for nucleons
  shiftlabels["ReweightFrCEx_N"]             = rwgt::fReweightFrCEx_N;     // charge exchange probability for nucleons, for given total rescattering probability
  shiftlabels["ReweightFrElas_N"]            = rwgt::fReweightFrElas_N;    // elastic    probability for nucleons, for given total rescattering probability
  shiftlabels["ReweightFrInel_N"]            = rwgt::fReweightFrInel_N;    // inelastic  probability for nucleons, for given total rescattering probability
  shiftlabels["ReweightFrAbs_N"]             = rwgt::fReweightFrAbs_N;     // absorption probability for nucleons, for given total rescattering probability
  shiftlabels["ReweightFrPiProd_N"]          = rwgt::fReweightFrPiProd_N;  // pion production probability for nucleons, for given total rescattering probability
  shiftlabels["ReweightMFP_pi"]              = rwgt::fReweightMFP_pi;      // mean free path for pions
  shiftlabels["ReweightFrCEx_pi"]            = rwgt::fReweightFrCEx_pi;    // charge exchange probability for pions, for given total rescattering probability
  shiftlabels["ReweightFrElas_pi"]           = rwgt::fReweightFrElas_pi;   // elastic   probability for pions, for given total rescattering probability
  shiftlabels["ReweightFrInel_pi"]           = rwgt::fReweightFrInel_pi;   // inelastic probability for pions, for given total rescattering probability
  shiftlabels["ReweightFrAbs_pi"]            = rwgt::fReweightFrAbs_pi;    // absorption probability for pions, for given total rescattering probability
  shiftlabels["ReweightFrPiProd_pi"]         = rwgt::fReweightFrPiProd_pi; // pion production probability for pions, for given total rescattering probability
  // Nuclear model

  //shiftlabels["ReweightNC"]                  = rwgt::fReweightNC; //

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

  for(const auto& shiftlabel : shiftlabels){
    for(const auto& sigma : sigmas) {
      systematics[shiftlabel.first][sigma] = new SystShifts(
        new GenieRwgtTableSyst(shiftlabel.second), sigma
      );
    }
  }

  // These two shifts are GENIE systs, but not in the form of reweights as above
  // Unfortunately they need a rwgt label, so give them the null one.
  shiftlabels["mecScale"] = rwgt::kReweightNull;
  shiftlabels["RPA"]      = rwgt::kReweightNull;

  for(const auto& sigma : sigmas) {
    systematics["mecScale"][sigma] = new SystShifts(new MECScaleSystSA(), sigma);
    systematics["RPA"][sigma]      = new SystShifts(new RPACCQESystSA(),  sigma);
  }

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

  // Cuts to split MC into two groups
  const Cut kMod0(
    [](const caf::SRProxy* sr)
    {
      if(sr->hdr.ismc == false) { return true; }
      int evtno = sr->hdr.evt;
      if((evtno % 2) == 0) { return true; }
      return false;
    });
  const Cut kMod1(
    [](const caf::SRProxy* sr)
    {
      if(sr->hdr.ismc == false) { return true; }
      int evtno = sr->hdr.evt;
      if((evtno % 2) == 1) { return true; }
      return false;
    });

  // 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) {
    // 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 && kMod0,
      kNoShift, kTuftsWeightCC
    );
    predsNE_nominal[sample.first] = new PredictionNoExtrap(
      loaders,
      *sample.second.first, *sample.second.second.second,
      kNoShift, kTuftsWeightCC
    );

    CheatDecomp* decompnc_n   = new CheatDecomp(
      loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
      *sample.second.first, *sample.second.second.first && kMod1,
      kNoShift, kTuftsWeightCC
    );
    CheatDecomp* decompnumu_n = new CheatDecomp(
      loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
      kNCBinsNumuCCAxis, kNumuND && kMod1,
      kNoShift, kTuftsWeightCC
    );
    ModularExtrapSterile* extrap_n = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
      loaders, *decompnc_n, *decompnumu_n,
      *sample.second.first, kNCBinsNumuCCAxis,
      *sample.second.second.second, *sample.second.second.first && kMod0, kNumuND && kMod0,
      kNoShift, kTuftsWeightCC
    ));
    predsSt_nominal[sample.first] = new PredictionSterile(extrap_n);

    // 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]) {
        decomps_shifted[sample.first][shiftlabel.first][syst.first] = new CheatDecomp(
          loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
          *sample.second.first, *sample.second.second.first && kMod0,
          *syst.second, kTuftsWeightCC
        );
        predsNE_shifted[sample.first][shiftlabel.first][syst.first] = new PredictionNoExtrap(
          loaders,
          *sample.second.first, *sample.second.second.second,
          *syst.second, kTuftsWeightCC
        );

        CheatDecomp* decompnc_s   = new CheatDecomp(
          loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
          *sample.second.first, *sample.second.second.first && kMod1,
          kNoShift, kTuftsWeightCC
        );
        CheatDecomp* decompnumu_s = new CheatDecomp(
          loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
          kNCBinsNumuCCAxis, kNumuND && kMod1,
          kNoShift, kTuftsWeightCC
        );
        ModularExtrapSterile* extrap_s = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
          loaders, *decompnc_s, *decompnumu_s,
          *sample.second.first, kNCBinsNumuCCAxis,
          *sample.second.second.second, *sample.second.second.first && kMod0, kNumuND && kMod0,
          *syst.second, kTuftsWeightCC
        ));
        predsSt_shifted[sample.first][shiftlabel.first][syst.first] = new PredictionSterile(
          extrap_s
        );
      }
    }
  }

  loaders.Go();

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

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