GenieWeights.h

Go to the documentation of this file.

/// GenieWeights.h:
///   Individual reaction cross section weights.
///   Details and references below.
///
///   If you are seeking global "tunes"
///   (combinations of weights used for specific analyses)
///   please look in XsecTunes.h.
///
///   Please don't build Vars here from CAFAna only.
///   All xsec weights should go in the NOvARwgt package
///   (so that they can be easily shared with the community as needed)
///   and be imported here using VarFromNOvARwgt().

#pragma once

#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Core/Var.h"

#include "NOvARwgt/rwgt/genie/QE/MAQEWgts.h"
#include "NOvARwgt/rwgt/genie/QE/RPAWeights.h"
#include "NOvARwgt/rwgt/genie/QE/NucOptPtlWgts.h"
#include "NOvARwgt/rwgt/genie/MEC/EmpiricalMECFixups.h"
#include "NOvARwgt/rwgt/genie/MEC/EmpiricalMECTuneSA.h"
#include "NOvARwgt/rwgt/genie/MEC/EmpiricalMECTune2017.h"
#include "NOvARwgt/rwgt/genie/MEC/EmpiricalMECTune2018.h"
#include "NOvARwgt/rwgt/genie/MEC/EmpiricalMECOtherTunes.h"
#include "NOvARwgt/rwgt/genie/DIS/HighWDISWeight.h"
#include "NOvARwgt/rwgt/genie/DIS/Nonres1piWeights.h"
#include "NOvARwgt/rwgt/genie/hN-FSI/hNFSIWeights.h"

#include "StandardRecord/NOvARwgt/NOvARwgtInterface.h"
#include "StandardRecord/SREnums.h"

namespace ana
{
  inline NuTruthVar VarFromNOvARwgt(const novarwgt::IWeightGenerator * wgtr)
  {
    return NuTruthVar([=](const caf::SRNeutrinoProxy * nu) -> double
                      {
                        auto evt = novarwgt::ConvertSRTruth(nu);

                        return wgtr->GetWeight(evt, {});
                      });
  }

  //---------------------------------------------------------------------------------

  template <typename T, typename... Args>
  inline NuTruthVar VarFromNOvARwgt(Args&&... args)
  {
    auto weighter = novarwgt::GetWeighter<T>(std::forward<Args>(args)...);

    return VarFromNOvARwgt(weighter);
  }

  //---------------------------------------------------------------------------------
  /// Adjust M_A for CCQE's central value to 1.04, per discussion following Doc 24705
  /// (GENIE 2.X default is 0.99)
  const NuTruthVar kRescaleMAQE_NT = VarFromNOvARwgt(novarwgt::kMAQEWeight_2018);
  const Var kRescaleMAQE = VarFromNuTruthVar(kRescaleMAQE_NT, 1);


  //---------------------------------------------------------------------------------
  // Nieves RPA Reweight (Minerva) - CC - q0 and q3 dependent
  // Histograms and code provided by Rik Gran
  // http://minerva-docdb.fnal.gov/cgi-bin/DocumentDatabase/  doc # 11276
  // http://arxiv.org/pdf/nucl-th/0408005.pdf
  // More documentation in NOvARwgt.

  /// Nieves RPA Reweight (Minerva) - CC QE - q0 and qmag dependent.  Second Analysis edition
  const NuTruthVar kRPAWeightCCQESA_NT = VarFromNOvARwgt(novarwgt::kRPAWeightCCQESA);
  const Var kRPAWeightCCQESA = VarFromNuTruthVar(kRPAWeightCCQESA_NT, 1);

  /// Generate a NuTruthVar that computes the RPA weight for this SRNeutino
  /// using the true Q^2.  Generally the other version (RPAWeightVar_q0q3) is to be
  /// preferred as it is more accurate, but occasionally this one is useful as well.
  const NuTruthVar RPAWeightVar_Q2 = VarFromNOvARwgt<novarwgt::RPAWeightQ2_2017>(novarwgt::kRxnCC, novarwgt::kScQuasiElastic);

  /// 2017 RPA effect for CCQE.
  /// Unfortunately it was mistakenly applied to ALL events...
  /// this is a ~3-5% effect for numu selected events.
  /// Second argument and beyond preserve buggy 2017 behavior
  /// (so that the prediction doesn't suddenly change underneath people):
  ///   * process==kUnknownMode ensures it will apply to all events
  ///   * forceNu==true ensures the "neutrino" histogram will be used even for antineutrinos
  ///   * includeHydrogen==true additionally applies the RPA effect to Hydrogen
  ///   * forceGENIE210==true uses the weights calculated for GENIE 2.10
  const NuTruthVar kRPAWeightCCQE2017_NT = VarFromNOvARwgt(novarwgt::kRPAWeightCCQE2017);
  const Var kRPAWeightCCQE2017 = VarFromNuTruthVar(kRPAWeightCCQE2017_NT, 1);

  /// 2018 RPA effect for CCQE.
  /// Same as 2017 version except that it is correctly applied only to QE
  /// and the neutrino/antineutrino differences in weights are correctly respected.
  const NuTruthVar kRPAWeightCCQE2018_NT = VarFromNOvARwgt(novarwgt::kRPAWeightCCQE2018);
  const Var kRPAWeightCCQE2018 = VarFromNuTruthVar(kRPAWeightCCQE2018_NT, 1);

  /// 2019 RPA effect for CCQE.
  /// Same as 2018 version except that the effect is not applied to hydrogen,
  /// and contains the fixes to correctly apply to GENIE 2.12
  const NuTruthVar kRPAWeightCCQE2019_NT = VarFromNOvARwgt(novarwgt::kRPAWeightCCQE2019);
  const Var kRPAWeightCCQE2019 = VarFromNuTruthVar(kRPAWeightCCQE2019_NT, 1);

  /// Valencia group RPA (2017 edition) applied to Resonant production rather than CCQE.
  /// For 2017-2018 analyses this was used as an uncertainty,
  /// because the Valencia calculation was done only for CCQE.
  /// Doing what we're doing is sort of reckless,
  /// but it seems to reproduce the suppression effect in pion production
  /// observed elsewhere (e.g., MiniBooNE, MINOS, MINERvA).
  /// Moreover, because it significantly improves data-MC agreement
  /// in our ND distributions, from 2018 we adopt it into the central value scheme.
  /// n.b.: Late in the 2018 analysis it was discovered that the RPA effect 
  /// was incorrectly applied to Hydrogen.  includeHydrogen==true reproduces this bug.
  const NuTruthVar kRPAWeightRES2017_NT = VarFromNOvARwgt(novarwgt::kRPAWeightRES2017);
  const Var kRPAWeightRES2017 = VarFromNuTruthVar(kRPAWeightRES2017_NT, 1);

  /// 2019 RPA effect for CCRES.
  /// Same as 2017/2018 version except that the effect is not applied to hydrogen
  const NuTruthVar kRPAWeightRES2019_NT = VarFromNOvARwgt(novarwgt::kRPAWeightRES2019);
  const Var kRPAWeightRES2019 = VarFromNuTruthVar(kRPAWeightRES2019_NT, 1);

  //---------------------------------------------------------------------------------
  /// Adjust non-resonant single pion production rate based on updated fits to deuterium data.
  /// (Also note that this doesn't adjust the uncertainties, which should be made somewhat smaller.)
  /// See Doc 14971.
  const NuTruthVar kFixNonres1PiST = VarFromNOvARwgt(novarwgt::kNonres1PiWgt);
  const Var kFixNonres1Pi = VarFromNuTruthVar(kFixNonres1PiST, 1);

  //---------------------------------------------------------------------------------
  /// Older versions of GENIE (prior to 2.11)
  /// need more corrections to Dytman's "EmpiricalMEC".
  extern const NuTruthCut kIsOldGENIEMEC;

  //---------------------------------------------------------------------------------
  /// fix distribution of initial state nucleon pairs of Dytman-MEC events:
  /// in GENIE prior to 2.11, generated 80%-20% nn-np; should be 20%-80%
  const NuTruthVar kDytmanMEC_FixItlStateST = VarFromNOvARwgt(novarwgt::kDytmanMECFixItlStateWgt);
  const Var kDytmanMEC_FixItlState = VarFromNuTruthVar(kDytmanMEC_FixItlStateST, 1);

  //---------------------------------------------------------------------------------
  /// fix wacky energy dependence of Dytman-MEC cross section in GENIE < 2.11 (see DocDB 14791)
  const NuTruthVar kDytmanMEC_FixXsecEdepST = VarFromNOvARwgt(novarwgt::kDytmanMECFixXsecEdepWgt);
  const Var kDytmanMEC_FixXsecEdep = VarFromNuTruthVar(kDytmanMEC_FixXsecEdepST, 1);

  //---------------------------------------------------------------------------------

  /// weight CC events to match NOvA ND data.
  /// basically reproduces Tony + Travis's prescription using Dytman-MEC.
  /// See Doc 14954 and references therein
  const NuTruthVar kTufts2p2hWgtSA_NT = VarFromNOvARwgt(novarwgt::kTufts2p2hWgtSA);
  const Var kTufts2p2hWgtSA  = VarFromNuTruthVar(kTufts2p2hWgtSA_NT, 1);


  //---------------------------------------------------------------------------------
  // alternative tunes for Empirical MEC to bound the shape for 2018 systematics

 /// Reweight Empirical MEC from fits to MC altered to be more 'qe-like' by varying existing systematic knobs, nu (reference needed)
  const NuTruthVar kMEC2018QElikeWgt_NT = VarFromNOvARwgt(novarwgt::kMEC2018_QElike_Wgt);
  const Var kMEC2018QElikeWgt = VarFromNuTruthVar(kMEC2018QElikeWgt_NT, 1);

  /// Reweight Empirical MEC from fits to MC altered to be more 'res-like' by varying existing systematic knobs, nu (reference needed)
  const NuTruthVar kMEC2018RESlikeWgt_NT = VarFromNOvARwgt(novarwgt::kMEC2018_RESlike_Wgt);
  const Var kMEC2018RESlikeWgt = VarFromNuTruthVar(kMEC2018RESlikeWgt_NT, 1);

  /// ------ MEC shape systematics for empirical MEC after Hydrogen RPA bug fixes
  /// Reweight Empirical MEC from fits to MC altered to be more 'qe-like' by varying existing systematic knobs, nu (reference needed)
  const NuTruthVar kMEC2018RPAFixQElikeWgt_NT = VarFromNOvARwgt(novarwgt::kMEC2018RPAFix_QElike_Wgt);
  const Var kMEC2018RPAFixQElikeWgt = VarFromNuTruthVar(kMEC2018RPAFixQElikeWgt_NT, 1);

  /// Reweight Empirical MEC from fits to MC altered to be more 'res-like' by varying existing systematic knobs, nu (reference needed)
  const NuTruthVar kMEC2018RPAFixRESlikeWgt_NT = VarFromNOvARwgt(novarwgt::kMEC2018RPAFix_RESlike_Wgt);
  const Var kMEC2018RPAFixRESlikeWgt = VarFromNuTruthVar(kMEC2018RPAFixRESlikeWgt_NT, 1);


  //---------------------------------------------------------------------------------
  // alternative tunes for Empirical MEC to bound the q0 response.  For systematics mostly.

  /// Reweight Empirical MEC to GENIE QE in q0, fixing q3 to data (see Doc 18741)
  const NuTruthVar kEmpiricalMECtoGENIEQEWgt_NT = VarFromNOvARwgt(novarwgt::kEmpiricalMEC_to_GENIEQE_Wgt);
  const Var kEmpiricalMECtoGENIEQEWgt = VarFromNuTruthVar(kEmpiricalMECtoGENIEQEWgt_NT, 1);

  /// Reweight Empirical MEC to GENIE RES in q0, fixing q3 to data (see Doc 18741)
  const NuTruthVar kEmpiricalMECtoGENIERESWgt_NT = VarFromNOvARwgt(novarwgt::kEmpiricalMEC_to_GENIERES_Wgt);
  const Var kEmpiricalMECtoGENIERESWgt = VarFromNuTruthVar(kEmpiricalMECtoGENIERESWgt_NT, 1);


  //---------------------------------------------------------------------------------
  // weight to turn Empirical MEC into Valencia MEC (or as close as we can get) in (q0,q3)
  const NuTruthVar kEmpiricalMECtoValenciaMECWgt_NT = VarFromNOvARwgt(novarwgt::kEmpiricalMEC_to_Valencia_Wgt);
  const Var kEmpiricalMECtoValenciaMECWgt = VarFromNuTruthVar(kEmpiricalMECtoValenciaMECWgt_NT, 1);

  //---------------------------------------------------------------------------------
  //----------------------------------------------------------------------------
  /// Implementation of various reweight schemes from MINERvA
  /// which attempt to fit to the excess in their data under different MEC hypotheses.
  /// There are four choices, each of which much be applied under different circumstances:
  ///   * "QE" --> weights applied only to QE.  MEC is weighted to 0 in this version.
  ///   * "MEC" --> weights applied to MEC.
  ///   * "MECNP" --> weights applied only to MEC NP component; MEC {NN,NP} component will be left alone.
  ///   * "MECPP" --> weights applied only to MEC NN,NP component; MEC NN component will be left alone.
  ///
  /// These are as publicized by R. Gran at NuInt 2017, https://meetings.triumf.ca/indico/event/6/session/4/contribution/77
  /// Citation forthcoming.  (Weight functions from K. McFarland in private communication to J. Wolcott.)
  /// They live here, rather than in NOvARwgt, because we don't have permission to distribute them.
  class MINERvA_MEC_Wgt
  {
      const std::string histFile = FindCAFAnaDir() + "/data/xs/MINERvA_MEC_wgts.root";
    public:
      MINERvA_MEC_Wgt(const std::string & whichWgt);

      double operator()(const caf::SRNeutrinoProxy * nu) const {return GetWeight(nu);}
    private:
      /// Draws the weight from the histogram
      double GetWeight(const caf::SRNeutrinoProxy * nu) const;

      const novarwgt::HistWrapper<TH2> fHist_nu;
      const novarwgt::HistWrapper<TH2> fHist_nubar;
      caf::mode_type_ fRxn;

      /// which nucleon pair to adjust.
      /// >0 means 'adjust pair with this pdg code';
      /// <0 means 'adjust anything but pairs with this pdg code;
      /// 0 means 'do nothing'
      int fHitNuc;
  };

  /// "1p1h" fit (QE was fit, MEC disabled)
  extern const NuTruthVar kMINERvA_Wgt_QE_NT;
  const Var kMINERvA_Wgt_QE = VarFromNuTruthVar(kMINERvA_Wgt_QE_NT, 1);

  /// "2p2h" fit (all MEC was fit)
  extern const NuTruthVar kMINERvA_Wgt_MEC_NT;
  const Var kMINERvA_Wgt_MEC = VarFromNuTruthVar(kMINERvA_Wgt_MEC_NT, 1);

  /// "2p2h-np" fit (only 'np' component of MEC was fit; 'nn' or 'pp' (depending on nu vs antinu) left alone)
  extern const NuTruthVar kMINERvA_Wgt_MECNP_NT;
  const Var kMINERvA_Wgt_MECNP = VarFromNuTruthVar(kMINERvA_Wgt_MECNP_NT, 1);

  /// "2p2h-pp" fit (either 'pp' or 'nn' component of MEC (depending on nu vs antinu) was fit; 'np' left alone)
  extern const NuTruthVar kMINERvA_Wgt_MECPP_NT;
  const Var kMINERvA_Wgt_MECPP = VarFromNuTruthVar(kMINERvA_Wgt_MECPP_NT, 1);

  //---------------------------------------------------------------------------------

  /// 2017 "best guess" for the MEC cross section (spill truth version)
  const NuTruthVar kEmpiricalMECWgt2017ST = VarFromNOvARwgt(novarwgt::kEmpiricalMECWgt2017);
  /// See kEmpiricalMECWgt2017ST
  const Var kEmpiricalMECWgt2017 = VarFromNuTruthVar(kEmpiricalMECWgt2017ST, 1);

  //---------------------------------------------------------------------------------

  /// 2018 tune for the MEC component of the cross section (NuTruth version)
  const NuTruthVar kEmpiricalMECWgt2018_NT = VarFromNOvARwgt(novarwgt::kEmpiricalMECWgt2018);
  /// See kEmpiricalMECWgt2018_NT
  const Var kEmpiricalMECWgt2018 = VarFromNuTruthVar(kEmpiricalMECWgt2018_NT, 1);

  //---------------------------------------------------------------------------------

  /// 2018 tune for the MEC component of the cross section after Hydrogen RPA bug fix
  const NuTruthVar kEmpiricalMECWgt2018RPAFix_NT = VarFromNOvARwgt(novarwgt::kEmpiricalMECWgt2018RPAFix);
  /// See kEmpiricalMECWgt2018_NT
  const Var kEmpiricalMECWgt2018RPAFix = VarFromNuTruthVar(kEmpiricalMECWgt2018RPAFix_NT, 1);

  //---------------------------------------------------------------------------------
  // weights obtained using miniprod5
  class EmpiricalMECWgt2020 : public novarwgt::MECq0q3NuNubarHistTuneWgt
  {
    public:
      explicit EmpiricalMECWgt2020(const novarwgt::IRegisterable::ClassID<EmpiricalMECWgt2020>& clID)
          : novarwgt::MECq0q3NuNubarHistTuneWgt(clID, "EmpiricalMECWgt2020",
                                                {novarwgt::StoredGenSupportCfg(novarwgt::GenCfg::kGENIE_Prod5Only)},
                                                FindCAFAnaDir()+"/data/xs/mec_fit_hists_fhc_bins_mp5.root", "mec_weights_smoothed",
                                                FindCAFAnaDir()+"/data/xs/mec_fit_hists_rhc_bins_mp5.root", "mec_weights_smoothed")
      {}
  };

  extern const EmpiricalMECWgt2020 * kEmpiricalMECWgt2020_Provisional;

  class ValenciaMECWgt2020_Provisional : public novarwgt::DoubleGaussMECWgt
  {
    public:
      explicit ValenciaMECWgt2020_Provisional(const novarwgt::IRegisterable::ClassID<ValenciaMECWgt2020_Provisional>& clID)
          : DoubleGaussMECWgt(clID, "ValenciaMECWgt2020_Provisional",
                              {novarwgt::StoredGenSupportCfg(novarwgt::GenCfg::kGENIE_Prod5Only)},
                              ValenciaMECWgt2020_Provisional::sParams)
      {}

    private:
      static constexpr DoubleGaussParameters sParams
          {
              12.78,     // "norm"
              0.33,   // "mean_q0"
              0.80,   // "mean_q3"
              0.10,   // "sigma_q0"
              0.29,   // "sigma_q3"
              0.91,       // "corr"

              38.8,     // "norm_2"
              0.036, // "mean_q0_2"
              0.45,  // "mean_q3_2"
              0.042, // "sigma_q0_2"
              0.233, // "sigma_q3_2"
              0.71,      // "corr_2"

              1        // "baseline"
                  };

  };

  extern const ValenciaMECWgt2020_Provisional * kValenciaMECWgt2020_Provisional;

  //---------------------------------------------------------------------------------

  /// Reweight DIS events at higher W; addresses the "short track anomaly"
  /// (see Doc 18027 and follow-ups)
  const NuTruthVar kRescaleHighWDIS_NT =  VarFromNOvARwgt(novarwgt::kHighWDISWgt_2018);
  const Var kRescaleHighWDIS = VarFromNuTruthVar(kRescaleHighWDIS_NT, 1);

  //---------------------------------------------------------------------------------

  /// Reweight GENIE 3.0.6 nominal FSI model to tuned CV (see Doc 42870)
  const NuTruthVar khNCVWgt2020_NT = VarFromNOvARwgt(novarwgt::khNCVWgt2020);
  const Var khNCVWgt2020 = VarFromNuTruthVar(khNCVWgt2020_NT, 1);

 


  //---------------------------------------------------------------------------------
   ///Reweight for Nuclear Optical Potential (see Doc ....) 
   const NuTruthVar kNucOptPtlWgt_NT = VarFromNOvARwgt(novarwgt::kNucOptPtl_wgt);
   const Var kNucOptPtlWgt = VarFromNuTruthVar(kNucOptPtlWgt_NT, 1);

} // namespace