Nus18SystUtilities.h

Go to the documentation of this file.

#pragma once

// framework
#include "StandardRecord/Proxy/SRProxy.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Systs/XSecSystLists.h"

// stl
#include <vector>
#include <map>
#include <string>

using namespace ana;

// -------------------------------------------------------------------------------------
std::map<std::string, std::vector<std::string> > GetBeamTranspShiftLabels() {

  // Map to set up systematics
  std::map<std::string, std::vector<std::string> > shiftLabels;

  // First element is a label for the systematic
  shiftLabels["BeamHornCurrent"] = std::vector<std::string>();
  shiftLabels["BeamSpotSize"]    = std::vector<std::string>();
  shiftLabels["BeamposX"]        = std::vector<std::string>();
  shiftLabels["BeamposY"]        = std::vector<std::string>();
  shiftLabels["BeamH1PosX"]      = std::vector<std::string>();
  shiftLabels["BeamH1PosY"]      = std::vector<std::string>();
  shiftLabels["BeamH2PosX"]      = std::vector<std::string>();
  shiftLabels["BeamH2PosY"]      = std::vector<std::string>();
  shiftLabels["TargetPosZ"]      = std::vector<std::string>();
  shiftLabels["BeamExp"]         = std::vector<std::string>();
  shiftLabels["HornWater"]       = std::vector<std::string>();
  /* shiftLabels["totErr"]          = std::vector<std::string>(); */

  // Second element is vector for items required to set up Syst object
  shiftLabels["BeamHornCurrent"].push_back("2kA");
  shiftLabels["BeamHornCurrent"].push_back("Horn Current");
  shiftLabels["BeamSpotSize"]   .push_back("0.2mmBeamSpotSize");
  shiftLabels["BeamSpotSize"]   .push_back("Spot Size");
  shiftLabels["BeamposX"]       .push_back("1mmBeamShiftX");
  shiftLabels["BeamposX"]       .push_back("Beam Position X");
  shiftLabels["BeamposY"]       .push_back("1mmBeamShiftY");
  shiftLabels["BeamposY"]       .push_back("Beam Position Y");
  shiftLabels["BeamH1PosX"]     .push_back("3mmHorn1X");
  shiftLabels["BeamH1PosX"]     .push_back("Horn 1 X Position");
  shiftLabels["BeamH1PosY"]     .push_back("3mmHorn1Y");
  shiftLabels["BeamH1PosY"]     .push_back("Horn 1 Y Position");
  shiftLabels["BeamH2PosX"]     .push_back("3mmHorn2X");
  shiftLabels["BeamH2PosX"]     .push_back("Horn 2 X Position");
  shiftLabels["BeamH2PosY"]     .push_back("3mmHorn2Y");
  shiftLabels["BeamH2PosY"]     .push_back("Horn 2 Y Position");
  shiftLabels["TargetPosZ"]     .push_back("7mmTargetZ");
  shiftLabels["TargetPosZ"]     .push_back("Target Z position");
  shiftLabels["BeamExp"]        .push_back("Magnetic Field in Decay Pipe");
  shiftLabels["BeamExp"]        .push_back("Magnet Field in Decay Pipe");
  shiftLabels["HornWater"]      .push_back("1mmHornWater");
  shiftLabels["HornWater"]      .push_back("New Horn Geometry and 1mm water");
  /* shiftLabels["totErr"]         .push_back("totErr"); */
  /* shiftLabels["totErr"]         .push_back("Combined Beam Transport Systematics"); */

  return shiftLabels;

}

// -------------------------------------------------------------------------------------
std::vector<Var> GetPPFXFluxWeights() {

  std::vector<Var> ppfxFluxUnivWgts;

  for (unsigned int UnivIdx = 0; UnivIdx < 100; UnivIdx++) {

    const NuTruthVar tempPPFXWgt([UnivIdx](const caf::SRNeutrinoProxy* nu) {
        if(nu->rwgt.ppfx.vuniv[UnivIdx] <= 0)
          return 1.f;
        return (float)nu->rwgt.ppfx.vuniv[UnivIdx];
      });

    const Var ppfxFluxUnivWgt = VarFromNuTruthVar(tempPPFXWgt, 1);
    ppfxFluxUnivWgts.push_back(ppfxFluxUnivWgt);

  }

  return ppfxFluxUnivWgts;

}

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

}