evwgh::GenieWeightCalc Class Reference

Inheritance diagram for evwgh::GenieWeightCalc:

Public Member Functions
	GenieWeightCalc ()
void	Configure (const fhicl::ParameterSet &pset, CLHEP::HepRandomEngine &engine) override
std::vector< std::vector< double > >	GetWeight (art::Event &e) override
Public Member Functions inherited from evwgh::WeightCalc
virtual void	Configure (fhicl::ParameterSet const &pset, CLHEP::HepRandomEngine &)=0
void	SetName (std::string name)
std::string	GetName ()

Private Member Functions
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)

Private Attributes
std::vector< genie::rew::GReWeight >	reweightVector
std::string	fGenieModuleLabel
bool	fQuietMode

Additional Inherited Members
Static Public Member Functions inherited from evwgh::WeightCalc
static std::vector< std::vector< double > >	MultiGaussianSmearing (std::vector< double > const &centralValues, std::vector< std::vector< double >> const &inputCovarianceMatrix, int n_multisims, CLHEP::RandGaussQ &GaussRandom)
	Applies Gaussian smearing to a set of data. More...
static std::vector< double >	MultiGaussianSmearing (std::vector< double > const &centralValue, TMatrixD *const &inputCovarianceMatrix, std::vector< double > rand)
	Over load of the above function that only returns a single varied parameter set. More...
static std::vector< double >	MultiGaussianSmearing (std::vector< double > const &centralValue, TMatrixD *const &LowerTriangleCovarianceMatrix, bool isDecomposed, std::vector< double > rand)

Detailed Description

Definition at line 223 of file GenieWeightCalc.cxx.

Constructor & Destructor Documentation

evwgh::GenieWeightCalc::GenieWeightCalc

(

)

Definition at line 248 of file GenieWeightCalc.cxx.

  {}

Member Function Documentation

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

private

Definition at line 599 of file GenieWeightCalc.cxx.

  {
    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 evwgh::GenieWeightCalc::Configure ( const fhicl::ParameterSet &  pset, CLHEP::HepRandomEngine &  engine  )

override

Definition at line 251 of file GenieWeightCalc.cxx.

  {
    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

  }

std::vector< std::vector< double > > evwgh::GenieWeightCalc::GetWeight ( art::Event &  e )

overridevirtual

Implements evwgh::WeightCalc.

Definition at line 535 of file GenieWeightCalc.cxx.

  {
    // Get the MC generator information out of the event
    // These are both handles to MC information.
    art::Handle< std::vector<simb::MCTruth> > mcTruthHandle;
    art::Handle< std::vector<simb::GTruth> > gTruthHandle;

    // Actually go and get the stuff
    e.getByLabel( fGenieModuleLabel, mcTruthHandle );
    e.getByLabel( fGenieModuleLabel, gTruthHandle );

    std::vector< art::Ptr<simb::MCTruth> > mclist;
    art::fill_ptr_vector( mclist, mcTruthHandle );

    std::vector< art::Ptr<simb::GTruth > > glist;
    art::fill_ptr_vector( glist, gTruthHandle );

    size_t num_neutrinos = mclist.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(*mclist[v], *glist[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;

  }

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

private

Definition at line 631 of file GenieWeightCalc.cxx.

  {
    // 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 << '\"';
    }
  }

Member Data Documentation

std::string evwgh::GenieWeightCalc::fGenieModuleLabel

private

Definition at line 241 of file GenieWeightCalc.cxx.

bool evwgh::GenieWeightCalc::fQuietMode

private

Definition at line 243 of file GenieWeightCalc.cxx.

std::vector< genie::rew::GReWeight > evwgh::GenieWeightCalc::reweightVector

private

Definition at line 239 of file GenieWeightCalc.cxx.

---

The documentation for this class was generated from the following file:

• larsim/larsim/EventWeight/Calculators/GenieWeightCalc.cxx