GenieWeightCalc.cxx

Go to the documentation of this file.

// GenieWeightCalc.cxx
//
// Handles event weights for GENIE systematics studies
//
// Updated for merge into larsim develop branch on Feb 22 2021 by Steven Gardiner
// Heavily rewritten on Dec 9 2019
//   by Steven Gardiner <gardiner@fnal.gov>
// Updated by Marco Del Tutto on Feb 18 2017
// Ported from uboonecode to larsim on Feb 14 2017
//   by Marco Del Tutto <marco.deltutto@physics.ox.ac.uk>

// Standard library includes
#include <map>
#include <memory>
#include <set>

// Framework includes
#include "art/Framework/Principal/Event.h"
#include "cetlib_except/exception.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "larsim/EventWeight/Base/WeightCalc.h"
#include "larsim/EventWeight/Base/WeightCalcCreator.h"

#include "CLHEP/Random/RandGaussQ.h"

#include "nugen/EventGeneratorBase/GENIE/GENIE2ART.h"

#include "nusimdata/SimulationBase/GTruth.h"
#include "nusimdata/SimulationBase/MCFlux.h"
#include "nusimdata/SimulationBase/MCTruth.h"

// GENIE includes
// TODO: add legacy support for GENIE v2 (mostly changes to header file
// locations would be needed)
#include "GENIE/Framework/Conventions/KineVar.h"
#include "GENIE/Framework/EventGen/EventRecord.h"
#include "GENIE/Framework/Interaction/Interaction.h"
#include "GENIE/Framework/Interaction/Kinematics.h"
#include "GENIE/Framework/Messenger/Messenger.h"
#include "GENIE/Framework/Utils/AppInit.h"

#include "GENIE/RwCalculators/GReWeightAGKY.h"
#include "GENIE/RwCalculators/GReWeightDISNuclMod.h"
#include "GENIE/RwCalculators/GReWeightFGM.h"
#include "GENIE/RwCalculators/GReWeightFZone.h"
#include "GENIE/RwCalculators/GReWeightINuke.h"
#include "GENIE/RwCalculators/GReWeightINukeParams.h"
#include "GENIE/RwCalculators/GReWeightNonResonanceBkg.h"
#include "GENIE/RwCalculators/GReWeightNuXSecCCQE.h"
#include "GENIE/RwCalculators/GReWeightNuXSecCCQEaxial.h"
#include "GENIE/RwCalculators/GReWeightNuXSecCCQEvec.h"
#include "GENIE/RwCalculators/GReWeightNuXSecCCRES.h"
#include "GENIE/RwCalculators/GReWeightNuXSecCOH.h"
#include "GENIE/RwCalculators/GReWeightNuXSecDIS.h"
#include "GENIE/RwCalculators/GReWeightNuXSecNC.h"
#include "GENIE/RwCalculators/GReWeightNuXSecNCEL.h"
#include "GENIE/RwCalculators/GReWeightNuXSecNCRES.h"
#include "GENIE/RwCalculators/GReWeightResonanceDecay.h"
#include "GENIE/RwCalculators/GReWeightXSecEmpiricalMEC.h"
#include "GENIE/RwFramework/GReWeight.h"
#include "GENIE/RwFramework/GSyst.h"
#include "GENIE/RwFramework/GSystSet.h"

// MicroBooNE-specific reweighting tools go here. These require a special genie
// ups product.
#ifdef GENIE_UB_PATCH

// New weight calculator in GENIE v3.0.4 MicroBooNE patch 01
#include "GENIE/RwCalculators/GReWeightXSecMEC.h"

// New weight calculators in GENIE v3.0.4 MicroBooNE patch 02
#include "GENIE/RwCalculators/GReWeightDeltaradAngle.h"
#include "GENIE/RwCalculators/GReWeightNuXSecCOHuB.h"
#include "GENIE/RwCalculators/GReWeightRESBugFix.h"

#endif

namespace {

  // These GENIE knobs are listed in the GSyst_t enum type but are not actually implemented.
  // They will be skipped and a warning message will be printed.
  // Last updated 9 Dec 2019 by S. Gardiner
  std::set<genie::rew::GSyst_t> UNIMPLEMENTED_GENIE_KNOBS = {
    kXSecTwkDial_RnubarnuCC,  // tweak the ratio of \sigma(\bar\nu CC) / \sigma(\nu CC)
    kXSecTwkDial_NormCCQEenu, // tweak CCQE normalization (maintains dependence on neutrino energy)
    kXSecTwkDial_NormDISCC,   // tweak the inclusive DIS CC normalization
    kXSecTwkDial_DISNuclMod   // unclear intent, does anyone else know? - S. Gardiner
  };

  // Some GENIE weight calculators can work with sets of knobs that should not be used simultaneously.
  // For instance, the CCQE weight calculator can vary parameters for the dipole axial form factor model
  // (the axial mass Ma) or for the z-expansion model, but using both together is invalid. The map below
  // is used to check that all of the requested knobs from the FHiCL input are compatible with each
  // other. Assuming they pass this check, the code will then configure the weight calculators to use
  // the correct "basis" of reweighting knobs as appropriate.

  // Outer keys are names of GENIE weight calculators that use sets of incompatible knobs. The names
  // need to match those used in GenieWeightCalc::SetupWeightCalculators(), specifically in the
  // calls to genie::rew::GReWeight::AdoptWeightCalc().
  // Inner keys are integer codes used to represent (and configure) each of the mutually exclusive
  // modes (i.e., to select one of the sets of incompatible knobs to use for the weight calculation).
  // Values are GSyst_t knob labels that belong to the given mode.
  std::map<std::string, std::map<int, std::set<genie::rew::GSyst_t>>> INCOMPATIBLE_GENIE_KNOBS = {
    // CCQE (genie::rew::GReWeightNuXSecCCQE)
    {"xsec_ccqe",
     {
       {genie::rew::GReWeightNuXSecCCQE::kModeNormAndMaShape,
        {
          // Norm + shape
          kXSecTwkDial_NormCCQE, // tweak CCQE normalization (energy independent)
          kXSecTwkDial_MaCCQEshape, // tweak Ma CCQE, affects dsigma(CCQE)/dQ2 in shape only (normalized to constant integral)
          kXSecTwkDial_E0CCQEshape // tweak E0 CCQE RunningMA, affects dsigma(CCQE)/dQ2 in shape only (normalized to constant integral)
        }},
       {genie::rew::GReWeightNuXSecCCQE::kModeMa,
        {
          // Ma
          kXSecTwkDial_MaCCQE, // tweak Ma CCQE, affects dsigma(CCQE)/dQ2 both in shape and normalization
          kXSecTwkDial_E0CCQE, // tweak E0 CCQE RunningMA, affects dsigma(CCQE)/dQ2 both in shape and normalization
        }},
       {genie::rew::GReWeightNuXSecCCQE::kModeZExp,
        {
          // Z-expansion
          kXSecTwkDial_ZNormCCQE, // tweak Z-expansion CCQE normalization (energy independent)
          kXSecTwkDial_ZExpA1CCQE, // tweak Z-expansion coefficient 1, affects dsigma(CCQE)/dQ2 both in shape and normalization
          kXSecTwkDial_ZExpA2CCQE, // tweak Z-expansion coefficient 2, affects dsigma(CCQE)/dQ2 both in shape and normalization
          kXSecTwkDial_ZExpA3CCQE, // tweak Z-expansion coefficient 3, affects dsigma(CCQE)/dQ2 both in shape and normalization
          kXSecTwkDial_ZExpA4CCQE // tweak Z-expansion coefficient 4, affects dsigma(CCQE)/dQ2 both in shape and normalization
        }},
     }},

    // CCRES (genie::rew::GReWeightNuXSecCCRES)
    {"xsec_ccres",
     {{genie::rew::GReWeightNuXSecCCRES::kModeNormAndMaMvShape,
       {
         // Norm + shape
         kXSecTwkDial_NormCCRES, 
         kXSecTwkDial_MaCCRESshape, 
         kXSecTwkDial_MvCCRESshape 
       }},
      {genie::rew::GReWeightNuXSecCCRES::kModeMaMv,
       {
         // Ma + Mv
         kXSecTwkDial_MaCCRES, // tweak Ma CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
         kXSecTwkDial_MvCCRES // tweak Mv CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
       }}}},

    // NCRES (genie::rew::GReWeightNuXSecNCRES)
    {"xsec_ncres",
     {{genie::rew::GReWeightNuXSecNCRES::kModeNormAndMaMvShape,
       {
         // Norm + shape
         kXSecTwkDial_NormNCRES, 
         kXSecTwkDial_MaNCRESshape, 
         kXSecTwkDial_MvNCRESshape 
       }},
      {genie::rew::GReWeightNuXSecNCRES::kModeMaMv,
       {
         // Ma + Mv
         kXSecTwkDial_MaNCRES, // tweak Ma NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
         kXSecTwkDial_MvNCRES // tweak Mv NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
       }}}},

    // DIS (genie::rew::GReWeightNuXSecDIS)
    {"xsec_dis",
     {{genie::rew::GReWeightNuXSecDIS::kModeABCV12u,
       {
         kXSecTwkDial_AhtBY, // tweak the Bodek-Yang model parameter A_{ht} - incl. both shape and normalization effect
         kXSecTwkDial_BhtBY, // tweak the Bodek-Yang model parameter B_{ht} - incl. both shape and normalization effect
         kXSecTwkDial_CV1uBY, // tweak the Bodek-Yang model parameter CV1u - incl. both shape and normalization effect
         kXSecTwkDial_CV2uBY // tweak the Bodek-Yang model parameter CV2u - incl. both shape and normalization effect
       }},
      {genie::rew::GReWeightNuXSecDIS::kModeABCV12uShape,
       {
         kXSecTwkDial_AhtBYshape, // tweak the Bodek-Yang model parameter A_{ht} - shape only effect to d2sigma(DIS)/dxdy
         kXSecTwkDial_BhtBYshape, // tweak the Bodek-Yang model parameter B_{ht} - shape only effect to d2sigma(DIS)/dxdy
         kXSecTwkDial_CV1uBYshape, // tweak the Bodek-Yang model parameter CV1u - shape only effect to d2sigma(DIS)/dxdy
         kXSecTwkDial_CV2uBYshape // tweak the Bodek-Yang model parameter CV2u - shape only effect to d2sigma(DIS)/dxdy
       }}}}};

  // Helper function that checks whether a string storing a GENIE knob name is
  // valid. Also stores the corresponding GSyst_t enum value for the knob.
  bool valid_knob_name(const std::string& knob_name, genie::rew::GSyst_t& knob)
  {
    knob = genie::rew::GSyst::FromString(knob_name);
    if (knob != kNullSystematic && knob != kNTwkDials) {
      if (UNIMPLEMENTED_GENIE_KNOBS.count(knob)) {
        MF_LOG_WARNING("GENIEWeightCalc") << "Ignoring unimplemented GENIE"
                                          << " knob " << knob_name;
        return false;
      }
    }
    else {
      throw cet::exception(__PRETTY_FUNCTION__) << "Encountered unrecognized"
                                                   "GENIE knob \""
                                                << knob_name << '\"';
      return false;
    }
    return true;
  }

  // Set of FHiCL weight calculator labels for which the tuned CV will be
  // ignored. If the name of the weight calculator doesn't appear in this set,
  // then variation weights will be thrown around the tuned CV.
  std::set<std::string> CALC_NAMES_THAT_IGNORE_TUNED_CV = {"RootinoFix"};

} // anonymous namespace

namespace evwgh {
  class GenieWeightCalc : public WeightCalc {
  public:
    GenieWeightCalc();
    void Configure(const fhicl::ParameterSet& pset, CLHEP::HepRandomEngine& engine) override;
    std::vector<std::vector<double>> GetWeight(art::Event& e) override;

  private:
    std::map<std::string, int> CheckForIncompatibleSystematics(
      const std::vector<genie::rew::GSyst_t>& knob_vec);

    void SetupWeightCalculators(genie::rew::GReWeight& rw,
                                const std::map<std::string, int>& modes_to_use);

    std::vector<genie::rew::GReWeight> reweightVector;

    std::string fGenieModuleLabel;

    bool fQuietMode;

    DECLARE_WEIGHTCALC(GenieWeightCalc)
  };

  GenieWeightCalc::GenieWeightCalc() {}

  void GenieWeightCalc::Configure(const fhicl::ParameterSet& p, CLHEP::HepRandomEngine& engine)
  {
    genie::Messenger* messenger = genie::Messenger::Instance();

    // By default, run GENIE reweighting in "quiet mode"
    // (use the logging settings in the "whisper" configuration
    // of the genie::Messenger class). The user can disable
    // quiet mode using the boolean FHiCL parameter "quiet_mode"
    fQuietMode = p.get<bool>("quiet_mode", true);
    if (fQuietMode) { genie::utils::app_init::MesgThresholds("Messenger_whisper.xml"); }
    else {
      // If quiet mode isn't enabled, then print detailed debugging
      // messages for GENIE reweighting
      messenger->SetPriorityLevel("ReW", log4cpp::Priority::DEBUG);

      MF_LOG_INFO("GENIEWeightCalc") << "Configuring GENIE weight"
                                     << " calculator " << this->GetName();
    }

    // Manually silence a couple of annoying GENIE logging messages
    // that appear a lot when running reweighting. This is done
    // whether or not "quiet mode" is enabled.
    messenger->SetPriorityLevel("TransverseEnhancementFFModel", log4cpp::Priority::WARN);
    messenger->SetPriorityLevel("Nieves", log4cpp::Priority::WARN);

    // Global Config
    fGenieModuleLabel = p.get<std::string>("genie_module_label");

    // Parameter central values from a table in the global config
    // Keys are GENIE knob names, values are knob settings that correspond to a tuned CV
    const fhicl::ParameterSet& param_CVs =
      p.get<fhicl::ParameterSet>("genie_central_values", fhicl::ParameterSet());
    std::vector<std::string> CV_knob_names = param_CVs.get_all_keys();

    // Map to store the CV knob settings
    std::map<genie::rew::GSyst_t, double> gsyst_to_cv_map;
    genie::rew::GSyst_t temp_knob;

    if (!CV_knob_names.empty() && !fQuietMode)
      MF_LOG_INFO("GENIEWeightCalc")
        << "Configuring non-default GENIE knob central values from input"
        << " FHiCL parameter set";

    for (const auto& knob_name : CV_knob_names) {
      double CV_knob_value = param_CVs.get<double>(knob_name);
      if (valid_knob_name(knob_name, temp_knob)) {
        if (gsyst_to_cv_map.count(temp_knob)) {
          throw cet::exception(__PRETTY_FUNCTION__)
            << "Duplicate central values"
            << " were configured for the " << knob_name << " GENIE knob.";
        }
        gsyst_to_cv_map[temp_knob] = CV_knob_value;
        if (!fQuietMode) {
          MF_LOG_INFO("GENIEWeightCalc")
            << "Central value for the " << knob_name << " knob was set to " << CV_knob_value;
        }
      }
    }

    // Calc Config
    const fhicl::ParameterSet& pset = p.get<fhicl::ParameterSet>(GetName());
    auto const pars =
      pset.get<std::vector<std::string>>("parameter_list", std::vector<std::string>());
    auto const par_sigmas = pset.get<std::vector<double>>("parameter_sigma", std::vector<double>());

    // Parameter limits (currently only used in minmax mode)
    // TODO: Revisit this. Could be useful for other modes.
    auto const par_mins = pset.get<std::vector<double>>("parameter_min", std::vector<double>());
    auto const par_maxes = pset.get<std::vector<double>>("parameter_max", std::vector<double>());

    auto const mode = pset.get<std::string>("mode");

    bool sigmas_ok = true;
    std::string array_name_for_exception;
    if (mode.find("central_value") == std::string::npos &&
        mode.find("minmax") == std::string::npos) {
      // For most reweighting modes, make sure that the number 1-sigma values
      // and the number of reweighting knobs match
      if (pars.size() != par_sigmas.size()) {
        sigmas_ok = false;
        array_name_for_exception = "parameter_sigma";
      }
    }
    else if (mode.find("minmax") != std::string::npos) {
      if (pars.size() != par_mins.size() || pars.size() != par_maxes.size()) {
        // For "minmax" mode, do the same for both the minimum and maximum
        // sigma values
        sigmas_ok = false;
        array_name_for_exception = "parameter_min and parameter_max";
      }
    }

    if (!sigmas_ok) {
      throw cet::exception(__PRETTY_FUNCTION__)
        << GetName() << "::Bad fcl configuration. parameter_list and " << array_name_for_exception
        << " need to have same number of parameters.";
    }

    if (!pars.empty() && !fQuietMode)
      MF_LOG_INFO("GENIEWeightCalc")
        << "Configuring"
        << " GENIE systematic knobs to be varied from the input FHiCL parameter set";

    // Convert the list of GENIE knob names from the input FHiCL configuration
    // into a vector of genie::rew::GSyst_t labels
    std::vector<genie::rew::GSyst_t> knobs_to_use;
    for (const auto& knob_name : pars) {
      if (valid_knob_name(knob_name, temp_knob)) knobs_to_use.push_back(temp_knob);
    }

    // We need to add all of the tuned CV knobs to the list when configuring
    // the weight calculators and checking for incompatibilities. Maybe all of
    // the systematic variation knobs are fine to use together, but they also
    // need to be compatible with the tuned CV. To perform the check, copy the
    // list of systematic variation knobs to use and add all the CV knobs that
    // aren't already present.
    std::vector<genie::rew::GSyst_t> all_knobs_vec = knobs_to_use;
    for (const auto& pair : gsyst_to_cv_map) {
      genie::rew::GSyst_t cv_knob = pair.first;
      auto begin = all_knobs_vec.cbegin();
      auto end = all_knobs_vec.cend();
      if (!std::count(begin, end, cv_knob)) all_knobs_vec.push_back(cv_knob);
    }

    // Check that the enabled knobs (both systematic variations and knobs used
    // for the CV tune) are all compatible with each other. The std::map
    // returned by this function call provides information needed to fine-tune
    // the configuration of the GENIE weight calculators.
    std::map<std::string, int> modes_to_use = this->CheckForIncompatibleSystematics(all_knobs_vec);

    // If we're working in "pm1sigma" or "minmax" mode, there should only be two
    // universes regardless of user input.
    size_t num_universes = 0u;
    if (mode.find("pm1sigma") != std::string::npos || mode.find("minmax") != std::string::npos) {
      num_universes = 2u;
    }
    else if (mode.find("multisim") != std::string::npos) {
      num_universes = pset.get<size_t>("number_of_multisims");

      // Since we're in multisim mode, force retrieval of the random
      // number seed. If it wasn't set, this will trigger an exception.
      // We want this check because otherwise the multisim universes
      // will not be easily reproducible.
      int dummy_seed = pset.get<int>("random_seed");
      MF_LOG_INFO("GENIEWeightCalc")
        << "GENIE weight calculator " << this->GetName() << " will generate " << num_universes
        << " multisim universes with random seed " << dummy_seed;
    }
    // If we're working in "central_value" or "default" mode, only a single
    // universe should be used
    else {
      num_universes = 1u;
    }

    // Create one default-constructed genie::rew::GReWeight object per universe
    reweightVector.resize(num_universes);

    // Set up the weight calculators for each universe
    for (auto& rwght : reweightVector) {
      this->SetupWeightCalculators(rwght, modes_to_use);
    }

    // Prepare sigmas
    size_t num_usable_knobs = knobs_to_use.size();
    std::vector<std::vector<double>> reweightingSigmas(num_usable_knobs);

    for (size_t k = 0u; k < num_usable_knobs; ++k) {
      reweightingSigmas[k].resize(num_universes);

      genie::rew::GSyst_t current_knob = knobs_to_use.at(k);

      for (size_t u = 0u; u < num_universes; ++u) {
        if (mode.find("multisim") != std::string::npos) {
          reweightingSigmas[k][u] = par_sigmas[k] * CLHEP::RandGaussQ::shoot(&engine, 0., 1.);
        }
        else if (mode.find("pm1sigma") != std::string::npos) {
          // u == 0 => +1*sigma; u == 1 => -1*sigma if pm1sigma is specified
          reweightingSigmas[k][u] = (u == 0 ? 1. : -1.) * par_sigmas.at(k);
        }
        else if (mode.find("minmax") != std::string::npos) {
          // u == 0 => max; u == 1 => min if minmax is specified
          reweightingSigmas[k][u] = (u == 0 ? par_maxes.at(k) : par_mins.at(k));
        }
        else if (mode.find("central_value") != std::string::npos) {
          // We'll correct for a modified CV below if needed
          reweightingSigmas[k][u] = 0.;
        }
        else {
          // By default, use the exact sigma value given for each knob
          reweightingSigmas[k][u] = par_sigmas[k];
        }

        if (!fQuietMode)
          MF_LOG_INFO("GENIEWeightCalc")
            << "Set sigma for the " << genie::rew::GSyst::AsString(current_knob)
            << " knob in universe #" << u << ". sigma = " << reweightingSigmas[k][u];

        // Add an offset if the central value for the current knob has been
        // configured (and is thus probably nonzero). Ignore this for minmax
        // mode (the limits should be chosen to respect a modified central
        // value)
        if (mode.find("minmax") == std::string::npos) {
          auto iter = gsyst_to_cv_map.find(current_knob);
          if (iter != gsyst_to_cv_map.end()) {
            reweightingSigmas[k][u] += iter->second;
            if (!fQuietMode)
              MF_LOG_INFO("GENIEWeightCalc")
                << "CV offset added to the " << genie::rew::GSyst::AsString(current_knob)
                << " knob. New sigma for universe #" << u << " is " << reweightingSigmas[k][u];
          }
        }
      }
    }

    // Don't adjust knobs to reflect the tuned CV if this weight calculator
    // should ignore those (as determined by whether it has one of the special
    // FHiCL names)
    if (!CALC_NAMES_THAT_IGNORE_TUNED_CV.count(this->GetName())) {

      // Add tuned CV knobs which have not been tweaked, and set them to their
      // modified central values. This ensures that weights are always thrown
      // around the modified CV.
      for (const auto& pair : gsyst_to_cv_map) {
        genie::rew::GSyst_t cv_knob = pair.first;
        double cv_value = pair.second;

        // If the current tuned CV knob is not present in the list of tweaked
        // knobs, then add it to the list with its tuned central value
        if (!std::count(knobs_to_use.cbegin(), knobs_to_use.cend(), cv_knob)) {
          ++num_usable_knobs;
          knobs_to_use.push_back(cv_knob);

          // The tuned CV knob will take the same value in every universe
          reweightingSigmas.emplace_back(std::vector<double>(num_universes, cv_value));

          if (!fQuietMode)
            MF_LOG_INFO("GENIEWeightCalc") << "Tuned knob " << genie::rew::GSyst::AsString(cv_knob)
                                           << " was not configured. Setting it to " << cv_value
                                           << " in all " << num_universes << " universes.";
        }
      }
    }

    // TODO: deal with parameters that have a priori bounds (e.g., FFCCQEVec,
    // which can vary on the interval [0,1])

    // Set up the knob values for each universe
    for (size_t u = 0; u < reweightVector.size(); ++u) {

      auto& rwght = reweightVector.at(u);
      genie::rew::GSystSet& syst = rwght.Systematics();

      for (unsigned int k = 0; k < knobs_to_use.size(); ++k) {
        genie::rew::GSyst_t knob = knobs_to_use.at(k);

        double twk_dial_value = reweightingSigmas.at(k).at(u);
        syst.Set(knob, twk_dial_value);

        if (!fQuietMode) {
          MF_LOG_INFO("GENIEWeightCalc") << "In universe #" << u << ", knob #" << k << " ("
                                         << genie::rew::GSyst::AsString(knob) << ") was set to"
                                         << " the value " << twk_dial_value;
        }
      } // loop over tweaked knobs

      rwght.Reconfigure();
      rwght.Print();
    } // loop over universes
  }

  // Returns a vector of weights for each neutrino interaction in the event
  std::vector<std::vector<double>> GenieWeightCalc::GetWeight(art::Event& e)
  {
    auto const& mcTruth = e.getProduct<std::vector<simb::MCTruth>>(fGenieModuleLabel);
    auto const& gTruth = e.getProduct<std::vector<simb::GTruth>>(fGenieModuleLabel);

    size_t num_neutrinos = mcTruth.size();
    size_t num_knobs = reweightVector.size();

    // Calculate weight(s) here
    std::vector<std::vector<double>> weights(num_neutrinos);
    for (size_t v = 0u; v < num_neutrinos; ++v) {

      // Convert the MCTruth and GTruth objects from the event
      // back into the original genie::EventRecord needed to
      // compute the weights
      std::unique_ptr<genie::EventRecord> genie_event(evgb::RetrieveGHEP(mcTruth[v], gTruth[v]));

      // Set the final lepton kinetic energy and scattering cosine
      // in the owned GENIE kinematics object. This is done during
      // event generation but is not reproduced by evgb::RetrieveGHEP().
      // Several new CCMEC weight calculators developed for MicroBooNE
      // expect the variables to be set in this way (so that differential
      // cross sections can be recomputed). Failing to set them results
      // in inf and NaN weights.
      // TODO: maybe update evgb::RetrieveGHEP to handle this instead.
      genie::Interaction* interaction = genie_event->Summary();
      genie::Kinematics* kine_ptr = interaction->KinePtr();

      // Final lepton mass
      double ml = interaction->FSPrimLepton()->Mass();
      // Final lepton 4-momentum
      const TLorentzVector& p4l = kine_ptr->FSLeptonP4();
      // Final lepton kinetic energy
      double Tl = p4l.E() - ml;
      // Final lepton scattering cosine
      double ctl = p4l.CosTheta();

      kine_ptr->SetKV(kKVTl, Tl);
      kine_ptr->SetKV(kKVctl, ctl);

      // All right, the event record is fully ready. Now ask the GReWeight
      // objects to compute the weights.
      weights[v].resize(num_knobs);
      for (size_t k = 0u; k < num_knobs; ++k) {
        weights[v][k] = reweightVector.at(k).CalcWeight(*genie_event);
      }
    }
    return weights;
  }

  std::map<std::string, int> GenieWeightCalc::CheckForIncompatibleSystematics(
    const std::vector<genie::rew::GSyst_t>& knob_vec)
  {
    std::map<std::string, int> modes_to_use;

    for (const auto& knob : knob_vec) {
      for (const auto& pair1 : INCOMPATIBLE_GENIE_KNOBS) {
        std::string calc_name = pair1.first;
        const auto& mode_map = pair1.second;
        for (const auto& pair2 : mode_map) {
          int mode = pair2.first;
          std::set<genie::rew::GSyst_t> knob_set = pair2.second;

          if (knob_set.count(knob)) {
            auto search = modes_to_use.find(calc_name);
            if (search != modes_to_use.end()) {
              if (search->second != mode) {
                auto knob_str = genie::rew::GSyst::AsString(knob);
                throw cet::exception(__PRETTY_FUNCTION__)
                  << this->GetName() << ": the GENIE knob " << knob_str << " is incompatible"
                  << " with others that are already configured";
              }
            }
            else
              modes_to_use[calc_name] = mode;
          }
        }
      }
    }

    return modes_to_use;
  }

  void GenieWeightCalc::SetupWeightCalculators(genie::rew::GReWeight& rw,
                                               const std::map<std::string, int>& modes_to_use)
  {
    // Based on the list from the GENIE command-line tool grwght1p
    rw.AdoptWghtCalc("xsec_ncel", new GReWeightNuXSecNCEL);
    rw.AdoptWghtCalc("xsec_ccqe", new GReWeightNuXSecCCQE);
    rw.AdoptWghtCalc("xsec_ccqe_axial", new GReWeightNuXSecCCQEaxial);
    rw.AdoptWghtCalc("xsec_ccqe_vec", new GReWeightNuXSecCCQEvec);
    rw.AdoptWghtCalc("xsec_ccres", new GReWeightNuXSecCCRES);
    rw.AdoptWghtCalc("xsec_ncres", new GReWeightNuXSecNCRES);
    rw.AdoptWghtCalc("xsec_nonresbkg", new GReWeightNonResonanceBkg);
    rw.AdoptWghtCalc("xsec_coh", new GReWeightNuXSecCOH);
    rw.AdoptWghtCalc("xsec_dis", new GReWeightNuXSecDIS);
    rw.AdoptWghtCalc("nuclear_qe", new GReWeightFGM);
    rw.AdoptWghtCalc("hadro_res_decay", new GReWeightResonanceDecay);
    rw.AdoptWghtCalc("hadro_fzone", new GReWeightFZone);
    rw.AdoptWghtCalc("hadro_intranuke", new GReWeightINuke);
    rw.AdoptWghtCalc("hadro_agky", new GReWeightAGKY);
    rw.AdoptWghtCalc("xsec_nc", new GReWeightNuXSecNC);
    rw.AdoptWghtCalc("res_dk", new GReWeightResonanceDecay);
    rw.AdoptWghtCalc("xsec_empmec", new GReWeightXSecEmpiricalMEC);
    // GReWeightDISNuclMod::CalcWeight() is not implemented, so we won't
    // bother to use it here. - S. Gardiner, 9 Dec 2019
    //rw.AdoptWghtCalc( "nuclear_dis",   new GReWeightDISNuclMod );

#ifdef GENIE_UB_PATCH

    // New weight calculator in GENIE v3.0.4 MicroBooNE patch 01
    rw.AdoptWghtCalc("xsec_mec", new GReWeightXSecMEC);

    // New weight calculators in GENIE v3.0.4 MicroBooNE patch 02
    rw.AdoptWghtCalc("deltarad_angle", new GReWeightDeltaradAngle);
    rw.AdoptWghtCalc("xsec_coh_ub", new GReWeightNuXSecCOHuB);
    rw.AdoptWghtCalc("res_bug_fix", new GReWeightRESBugFix);

#endif

    // Set the modes for the weight calculators that need them to be specified
    for (const auto& pair : modes_to_use) {
      std::string calc_name = pair.first;
      int mode = pair.second;

      genie::rew::GReWeightI* calc = rw.WghtCalc(calc_name);
      if (!calc)
        throw cet::exception(__PRETTY_FUNCTION__)
          << "Failed to retrieve the GENIE weight calculator labeled \"" << calc_name << '\"';

      // The GReWeightI base class doesn't have a SetMode(int) function,
      // so we'll just try dynamic casting until we get the right one.
      // If none work, then throw an exception.
      // TODO: Add a virtual function GReWeightI::SetMode( int ) in GENIE's
      // Reweight framework. Then we can avoid the hacky dynamic casts here.
      auto* calc_ccqe = dynamic_cast<genie::rew::GReWeightNuXSecCCQE*>(calc);
      auto* calc_ccres = dynamic_cast<genie::rew::GReWeightNuXSecCCRES*>(calc);
      auto* calc_ncres = dynamic_cast<genie::rew::GReWeightNuXSecNCRES*>(calc);
      auto* calc_dis = dynamic_cast<genie::rew::GReWeightNuXSecDIS*>(calc);
      if (calc_ccqe)
        calc_ccqe->SetMode(mode);
      else if (calc_ccres)
        calc_ccres->SetMode(mode);
      else if (calc_ncres)
        calc_ncres->SetMode(mode);
      else if (calc_dis)
        calc_dis->SetMode(mode);
      else
        throw cet::exception(__PRETTY_FUNCTION__)
          << "Request to set the mode of an unrecognized GENIE weight calculator \"" << calc_name
          << '\"';
    }
  }

  REGISTER_WEIGHTCALC(GenieWeightCalc)
}