Inheritance diagram for sbn::evwgh::GenieWeightCalc:

Public Member Functions
	GenieWeightCalc ()

void	Configure (fhicl::ParameterSet const &pset, CLHEP::HepRandomEngine &engine) override

std::vector< float >	GetWeight (art::Event &e, size_t inu) override

Public Member Functions inherited from sbn::evwgh::WeightCalc
void	SetName (std::string name)

void	SetType (std::string type)

std::string	GetName ()

std::string	GetType ()

std::string	GetFullName ()

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

std::string	fTuneName

bool	fQuietMode

Additional Inherited Members
Public Attributes inherited from sbn::evwgh::WeightCalc
EventWeightParameterSet	fParameterSet

Detailed Description

Definition at line 229 of file sbncode/sbncode/SBNEventWeight/Calculators/CrossSections/GenieWeightCalc.cxx.

Constructor & Destructor Documentation

sbn::evwgh::GenieWeightCalc::GenieWeightCalc ( )

inline

Definition at line 231 of file sbncode/sbncode/SBNEventWeight/Calculators/CrossSections/GenieWeightCalc.cxx.

231 : WeightCalc() {}

Member Function Documentation

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

private

Definition at line 564 of file sbncode/sbncode/SBNEventWeight/Calculators/CrossSections/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;//<int> code of the variable
           std::set<genie::rew::GSyst_t> knob_set = pair2.second;//incomp. knobs
 
           if ( knob_set.count(knob) ) {//<string> leads to incomp. knobs
             auto search = modes_to_use.find( calc_name );
             if ( search != modes_to_use.end() ) {//found a incomp. knob's name
               if ( search->second != mode ) {//if the knob is not connected to the variable?
                 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 sbn::evwgh::GenieWeightCalc::Configure	(	fhicl::ParameterSet const &	pset,
		CLHEP::HepRandomEngine &	engine
	)

overridevirtual

Implements sbn::evwgh::WeightCalc.

Definition at line 256 of file sbncode/sbncode/SBNEventWeight/Calculators/CrossSections/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
   //----Initialize the tune
   // Configure the appropriate GENIE tune if needed (important for v3+ only)
   // NOTE: In all normal use cases, relying on the ${GENIE_XSEC_TUNE} environment
   // variable set by the genie_xsec package should be sufficient. Only include
   // the "TuneName" FHiCL parameter for EventWeight if you really know what you're
   // doing! The same goes for the "EventGeneratorList" parameter.
   std::string genie_tune_name = p.get<std::string>("TuneName", "${GENIE_XSEC_TUNE}");
 
   // The default empty string used here will cause the subsequent call to
   // evgb::SetEventGeneratorListAndTune() to leave GENIE's current event generator
   // list name (probably "Default") alone
   std::string evgen_list_name = p.get<std::string>("EventGeneratorList", "");
 
   // Tell GENIE about the event generator list and tune
   evgb::SetEventGeneratorListAndTune( evgen_list_name, genie_tune_name );
 
 
   fGenieModuleLabel = p.get<std::string>("genie_module_label");
 
   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();
 
   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;
       }
     }
   }
 
   // Calculator config
   const fhicl::ParameterSet& pset = p.get<fhicl::ParameterSet>(GetName());
 
   auto const& pars = pset.get<std::vector<std::string> >("parameter_list");
 
   std::vector<float> parsigmas = pset.get<std::vector<float> >("parameter_sigma", std::vector<float>() );
 
 
   // 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 );
   }
 
   //check if all_knobs_vec are good
   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 );
   }
   std::map< std::string, int >
     modes_to_use = this->CheckForIncompatibleSystematics( all_knobs_vec );
 
 
   std::string mode = pset.get<std::string>("mode");
 
   int num_universes = 1; //for "central_value" or "default" mode
 
   std::string array_name_for_exception;
 
   if (pars.size() != parsigmas.size()) {
     array_name_for_exception = "parameter_sigma";
   }
   //One fParameterSet gives one set of weights; 
   //Different fParameterSet's are independent.
   if(mode.find("pmNsigma") != std::string::npos){  
     //This mode will be updated and replaced with a more useful one.
 
     num_universes = 2;
 
     for (size_t i=0; i<pars.size(); i++){ 
       fParameterSet.AddParameter(pars[i], parsigmas[i]);
     }
 
   }else if(mode.find("multisigma") != std::string::npos){
 
     //==== pars.size()==1 is supported
     array_name_for_exception = "";
     if(pars.size()!=1){
       array_name_for_exception = "multisigma";
       throw cet::exception(__PRETTY_FUNCTION__) << GetName()
         << "::Bad fcl configuration. multisigma only supports one parameter_list";
     }
 
     for (size_t i=0; i<pars.size(); i++){
       num_universes = pars.size() * parsigmas.size();
       std::vector<float> withds;
       for (size_t j=0; j<parsigmas.size(); j++){
         withds.push_back(parsigmas[j]);
       }
       fParameterSet.AddParameter(pars[i], withds);
     }
 
   }else if(mode.find("multisim") != std::string::npos){
 
     num_universes = pset.get<int>( "number_of_multisims", 1 );
 
     MF_LOG_INFO("GENIEWeightCalc") << "GENIE weight calculator "
       << this->GetName() << " will generate " << num_universes
       << " multisim universes";
 
     for (size_t i=0; i<pars.size(); i++) {
       fParameterSet.AddParameter(pars[i], parsigmas[i]);
     }
 
   }else if(mode.find("fixed") != std::string::npos){
 
     num_universes = 1;
 
     for (size_t i=0; i<pars.size(); i++) {
       fParameterSet.AddParameter(pars[i], parsigmas[i]);
     }
 
   }else { //need to identify a mode
 
     throw cet::exception(__PRETTY_FUNCTION__) << GetName()
       <<":: Mode is not recognized: "<<mode;
 
   }
 
   if( array_name_for_exception.size() > 1){
 
     throw cet::exception(__PRETTY_FUNCTION__) << GetName()
       << "::Bad fcl configuration. parameter_list and "
       << array_name_for_exception
       << " need to have same number of parameters.";
   }
 
   // Set up parameters
   fParameterSet.Configure(GetFullName(), mode, num_universes);
   fParameterSet.Sample(engine);
 
   //Initialize  genie::rew::GReWeight's
   reweightVector.resize( num_universes );// Reconfigure this later in this function.
   std::cout << GetFullName() << ": Initialize WeightCalculators with "<<num_universes<<" univeres."<<std::endl;
 
   // Initialize weight calculators for all universes;
   for ( auto& rwght : reweightVector ){ this->SetupWeightCalculators( rwght, modes_to_use );}
 
   // Set up parameters
   for ( size_t univ = 0; univ < reweightVector.size(); ++univ ) {//loop over universes
 
     auto& rwght = reweightVector.at( univ );
     genie::rew::GSystSet& syst = rwght.Systematics();
 
     rwght.Reconfigure();//Redundant, but it avoids "cushion error":
     //1638551728 FATAL ReW : [n] <GReWeightINukeParams.cxx::AddCushionTerms (583)> : There must be at least one cushion term (0 were set)
     //loop over parameters;
     for (auto const& it : fParameterSet.fParameterMap) {//loop over knobs
       //member of std::map<EventWeightParameter, std::vector<float> >
       //(see line 41 of sbnobj/Common/SBNEventWeight/EventWeightParameterSet.h)
       std::string name = it.first.fName;
       genie::rew::GSyst_t knob = genie::rew::GSyst::FromString( name );
 
       double cv_shift = 0;
       auto iter = gsyst_to_cv_map.find( knob ); //Apply CV shift, if required
       if ( iter != gsyst_to_cv_map.end() ) {
         cv_shift =  iter->second;
         if ( !fQuietMode ) MF_LOG_INFO("GENIEWeightCalc")
           << "CV offset added to the "
             << genie::rew::GSyst::AsString( knob )
             << " knob. New sigma for universe #" << univ << " is "
             << cv_shift;
       }
 
       double twk_dial_value = it.second[univ]+cv_shift;
       syst.Set( knob, twk_dial_value ); //Assign knob with new values based on variation
 
       std::cout << GetFullName() << ": univ "<<univ<<" "<<genie::rew::GSyst::AsString( knob )<<" Knob value: "<<twk_dial_value<<std::endl;
 
       if ( !fQuietMode ) {
         MF_LOG_INFO("GENIEWeightCalc") << "In universe #" << univ << ", knob name" << name // << k
           << " (" << genie::rew::GSyst::AsString( knob ) << ") was set to"
           << " the value " << twk_dial_value << " with shift "<<cv_shift;
       }
 
     }//next knob
     rwght.Reconfigure();//Apply all set knobs, i.e. use the updated syst.
     rwght.Print();
   }//next universe
 }

std::vector< float > sbn::evwgh::GenieWeightCalc::GetWeight	(	art::Event &	e,
		size_t	inu
	)

overridevirtual

Implements sbn::evwgh::WeightCalc.

Definition at line 489 of file sbncode/sbncode/SBNEventWeight/Calculators/CrossSections/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_knobs = reweightVector.size();
 
     // Calculate weight(s) here
     std::vector< float > weights( num_knobs );
 //    std::vector< std::vector<double> > weights( num_neutrinos );
 
       // 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[inu], *glist[inu]) );
 
       // 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.
 
     for (size_t k = 0u; k < num_knobs; ++k ) {//one "knob" for one universe;
       //Add exception to avoid "FATAL KineLimits", 
       //  see https://github.com/GENIE-MC/Reweight/issues/12
       std::vector< genie::rew::GSyst_t > ak = reweightVector.at( k ).Systematics().AllIncluded();//all konbs used
 
       //NOTE: this line is to skip reweighting stange-CCQE events 
       //Among 7500 events, 207 of them are strange events; ~3%
       if(glist[inu]->fIsStrange && std::find(ak.begin(), ak.end(),  kXSecTwkDial_MaCCQE) != ak.end()){
         weights[k] = 1;
       } else{
         weights[k] = reweightVector.at( k ).CalcWeight( *genie_event );
         reweightVector.at(k).Print();
       }
 //  std::cout<<"CHECK "<<k<<" universe weight "<<weights[k]<<std::endl;
     }
 //  std::cout<<"Finish Event # "<<e.event()<<std::endl;
 
     return weights;
 }

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

private

Definition at line 596 of file sbncode/sbncode/SBNEventWeight/Calculators/CrossSections/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 << '\"';
   }
 
 }