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
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 larsim/larsim/EventWeight/Calculators/GenieWeightCalc.cxx.

Constructor & Destructor Documentation

evwgh::GenieWeightCalc::GenieWeightCalc ( )

Definition at line 248 of file larsim/larsim/EventWeight/Calculators/GenieWeightCalc.cxx.

249 {}

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 larsim/larsim/EventWeight/Calculators/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
	)

overridevirtual

Implements evwgh::WeightCalc.

Definition at line 251 of file larsim/larsim/EventWeight/Calculators/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 larsim/larsim/EventWeight/Calculators/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 larsim/larsim/EventWeight/Calculators/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 << '\"';
     }
   }