diff --git a/CepGen/Cards/LpairHandler.cpp b/CepGen/Cards/LpairHandler.cpp
index 9d6ebe2..12bbc09 100644
--- a/CepGen/Cards/LpairHandler.cpp
+++ b/CepGen/Cards/LpairHandler.cpp
@@ -1,212 +1,232 @@
 #include "CepGen/Cards/LpairHandler.h"
 #include "CepGen/Core/ParametersList.h"
 #include "CepGen/Core/Exception.h"
 
 #include "CepGen/Physics/PDG.h"
+#include "CepGen/StructureFunctions/StructureFunctionsBuilder.h"
+#include "CepGen/StructureFunctions/LHAPDF.h"
+
 #include "CepGen/Processes/GamGamLL.h"
 #include "CepGen/Processes/PPtoFF.h"
 #include "CepGen/Processes/PPtoWW.h"
 #include "CepGen/Processes/FortranProcesses.h"
 
 #include "CepGen/Hadronisers/Pythia8Hadroniser.h"
 
 #include <fstream>
 
 namespace CepGen
 {
   namespace Cards
   {
     const int LpairHandler::kInvalid = 99999;
 
     //----- specialization for LPAIR input cards
 
     LpairHandler::LpairHandler( const char* file ) :
       proc_params_( new ParametersList ),
-      hi_1_( { 0, 0 } ), hi_2_( { 0, 0 } )
+      str_fun_( 11 ), hi_1_( { 0, 0 } ), hi_2_( { 0, 0 } )
     {
       std::ifstream f( file, std::fstream::in );
       if ( !f.is_open() )
         throw CG_FATAL( "LpairHandler" ) << "Failed to parse file \"" << file << "%s\".";
 
       init( &params_ );
 
+      //--- parse all fields
       std::unordered_map<std::string, std::string> m_params;
       std::string key, value;
       std::ostringstream os;
       while ( f >> key >> value ) {
         if ( key[0] == '#' ) continue; // FIXME need to ensure there is no extra space before!
         setParameter( key, value );
         m_params.insert( { key, value } );
         if ( getDescription( key ) != "null" )
           os << "\n>> " << key << " = " << std::setw( 15 ) << getParameter( key )
              << " (" << getDescription( key ) << ")";
       }
       f.close();
 
+      //--- parse the process name
       if ( proc_name_ == "lpair" )
         params_.setProcess( new Process::GamGamLL( *proc_params_ ) );
       else if ( proc_name_ == "pptoll" || proc_name_ == "pptoff" )
         params_.setProcess( new Process::PPtoFF( *proc_params_ ) );
       else if ( proc_name_ == "pptoww" )
         params_.setProcess( new Process::PPtoWW( *proc_params_ ) );
       else {
         Process::generateFortranProcesses();
         for ( auto& proc : Process::FortranProcessesHandler::get().list() )
           if ( proc_name_ == std::string( proc.name ) )
             params_.setProcess( new Process::FortranKTProcess( *proc_params_, proc.name, proc.description, proc.method ) );
         if ( !params_.process() )
           throw CG_FATAL( "LpairHandler" ) << "Unrecognised process name: " << proc_name_ << "!";
       }
 
+      //--- parse the structure functions code
+      const unsigned long kLHAPDFCodeDec = 10000000, kLHAPDFPartDec = 1000000;
+      if ( str_fun_ / kLHAPDFCodeDec == 1 ) { // SF from parton
+        params_.kinematics.structure_functions = StructureFunctionsBuilder::get( SF::Type::LHAPDF );
+        auto sf = dynamic_cast<SF::LHAPDF*>( params_.kinematics.structure_functions.get() );
+        const unsigned long icode = str_fun_ % kLHAPDFCodeDec;
+        sf->params.pdf_code = icode % kLHAPDFPartDec;
+        sf->params.mode = (SF::LHAPDF::Parameterisation::Mode)( icode / kLHAPDFPartDec ); // 0, 1, 2
+      }
+      else
+        params_.kinematics.structure_functions = StructureFunctionsBuilder::get( (SF::Type)str_fun_ );
+
+      //--- parse the integration algorithm name
       if ( integr_type_ == "plain" )
         params_.integrator.type = Integrator::Type::plain;
       else if ( integr_type_ == "Vegas" )
         params_.integrator.type = Integrator::Type::Vegas;
       else if ( integr_type_ == "MISER" )
         params_.integrator.type = Integrator::Type::MISER;
       else if ( integr_type_ != "" )
         throw CG_FATAL( "LpairHandler" ) << "Unrecognized integrator type: " << integr_type_ << "!";
 
+      //--- parse the hadronisation algorithm name
       if ( hadr_name_ == "pythia8" )
         params_.setHadroniser( new Hadroniser::Pythia8Hadroniser( params_ ) );
 
       if ( m_params.count( "IEND" ) )
         setValue<bool>( "IEND", ( std::stoi( m_params["IEND"] ) > 1 ) );
 
+      //--- check if we are dealing with heavy ions for incoming states
       HeavyIon hi1{ hi_1_.first, (Element)hi_1_.second }, hi2{ hi_2_.first, (Element)hi_2_.second };
       if ( hi1 )
         params_.kinematics.incoming_beams.first.pdg = hi1;
       if ( hi2 )
         params_.kinematics.incoming_beams.second.pdg = hi2;
 
       CG_INFO( "LpairHandler" ) << "File '" << file << "' succesfully opened!\n\t"
         << "The following parameters are set:" << os.str();
     }
 
     void
     LpairHandler::init( Parameters* params )
     {
       //-------------------------------------------------------------------------------------------
       // Process/integration/hadronisation parameters
       //-------------------------------------------------------------------------------------------
 
       registerParameter<std::string>( "PROC", "Process name to simulate", &proc_name_ );
       registerParameter<std::string>( "ITYP", "Integration algorithm", &integr_type_ );
       registerParameter<std::string>( "HADR", "Hadronisation algorithm", &hadr_name_ );
       registerParameter<std::string>( "KMRG", "KMR grid interpolation path", &params_.kinematics.kmr_grid_path );
 
       //-------------------------------------------------------------------------------------------
       // General parameters
       //-------------------------------------------------------------------------------------------
 
       registerParameter<bool>( "IEND", "Generation type", &params->generation.enabled );
       registerParameter<bool>( "NTRT", "Smoothen the integrand", &params->generation.treat );
       registerParameter<int>( "DEBG", "Debugging verbosity", (int*)&Logger::get().level );
       registerParameter<int>( "NCVG", "Number of function calls", (int*)&params->integrator.ncvg );
       registerParameter<int>( "ITVG", "Number of integration iterations", (int*)&params->integrator.vegas.iterations );
       registerParameter<int>( "SEED", "Random generator seed", (int*)&params->integrator.rng_seed );
       registerParameter<int>( "NTHR", "Number of threads to use for events generation", (int*)&params->generation.num_threads );
       registerParameter<int>( "MODE", "Subprocess' mode", (int*)&params->kinematics.mode );
       registerParameter<int>( "NCSG", "Number of points to probe", (int*)&params->generation.num_points );
       registerParameter<int>( "NGEN", "Number of events to generate", (int*)&params->generation.maxgen );
       registerParameter<int>( "NPRN", "Number of events before printout", (int*)&params->generation.gen_print_every );
 
       //-------------------------------------------------------------------------------------------
       // Process-specific parameters
       //-------------------------------------------------------------------------------------------
 
       registerParameter<int>( "METH", "Computation method (kT-factorisation)", &proc_params_->operator[]<int>( "method" ) );
       registerParameter<int>( "IPOL", "Polarisation states to consider", &proc_params_->operator[]<int>( "polarisationStates" ) );
 
       //-------------------------------------------------------------------------------------------
       // Process kinematics parameters
       //-------------------------------------------------------------------------------------------
 
-      registerParameter<int>( "PMOD", "Outgoing primary particles' mode", (int*)&params->kinematics.structure_functions );
-      registerParameter<int>( "EMOD", "Outgoing primary particles' mode", (int*)&params->kinematics.structure_functions );
+      registerParameter<int>( "PMOD", "Outgoing primary particles' mode", &str_fun_ );
+      registerParameter<int>( "EMOD", "Outgoing primary particles' mode", &str_fun_ );
       registerParameter<int>( "PAIR", "Outgoing particles' PDG id", (int*)&proc_params_->operator[]<int>( "pair" ) );
 
       registerParameter<int>( "INA1", "Heavy ion atomic weight (1st incoming beam)", (int*)&hi_1_.first );
       registerParameter<int>( "INZ1", "Heavy ion atomic number (1st incoming beam)", (int*)&hi_1_.second );
       registerParameter<int>( "INA2", "Heavy ion atomic weight (1st incoming beam)", (int*)&hi_2_.first );
       registerParameter<int>( "INZ2", "Heavy ion atomic number (1st incoming beam)", (int*)&hi_2_.second );
       registerParameter<double>( "INP1", "Momentum (1st primary particle)", &params->kinematics.incoming_beams.first.pz );
       registerParameter<double>( "INP2", "Momentum (2nd primary particle)", &params->kinematics.incoming_beams.second.pz );
       registerParameter<double>( "INPP", "Momentum (1st primary particle)", &params->kinematics.incoming_beams.first.pz );
       registerParameter<double>( "INPE", "Momentum (2nd primary particle)", &params->kinematics.incoming_beams.second.pz );
       registerParameter<double>( "PTCT", "Minimal transverse momentum (single central outgoing particle)", &params->kinematics.cuts.central.pt_single.min() );
       registerParameter<double>( "MSCT", "Minimal central system mass", &params->kinematics.cuts.central.mass_sum.min() );
       registerParameter<double>( "ECUT", "Minimal energy (single central outgoing particle)", &params->kinematics.cuts.central.energy_single.min() );
       registerParameter<double>( "ETMN", "Minimal pseudo-rapidity (central outgoing particles)", &params->kinematics.cuts.central.eta_single.min() );
       registerParameter<double>( "ETMX", "Maximal pseudo-rapidity (central outgoing particles)", &params->kinematics.cuts.central.eta_single.max() );
       registerParameter<double>( "YMIN", "Minimal rapidity (central outgoing particles)", &params->kinematics.cuts.central.rapidity_single.min() );
       registerParameter<double>( "YMAX", "Maximal rapidity (central outgoing particles)", &params->kinematics.cuts.central.rapidity_single.max() );
       registerParameter<double>( "Q2MN", "Minimal Q² = -q² (exchanged parton)", &params->kinematics.cuts.initial.q2.min() );
       registerParameter<double>( "Q2MX", "Maximal Q² = -q² (exchanged parton)", &params->kinematics.cuts.initial.q2.max() );
       registerParameter<double>( "MXMN", "Minimal invariant mass of proton remnants", &params->kinematics.cuts.remnants.mass_single.min() );
       registerParameter<double>( "MXMX", "Maximal invariant mass of proton remnants", &params->kinematics.cuts.remnants.mass_single.max() );
     }
 
     void
     LpairHandler::store( const char* file )
     {
       std::ofstream f( file, std::fstream::out | std::fstream::trunc );
       if ( !f.is_open() ) {
         CG_ERROR( "LpairHandler" ) << "Failed to open file \"" << file << "%s\" for writing.";
       }
       for ( const auto& it : p_strings_ )
         if ( it.second.value )
           f << it.first << " = " << *it.second.value << "\n";
       for ( const auto& it : p_ints_ )
         if ( it.second.value )
           f << it.first << " = " << *it.second.value << "\n";
       for ( const auto& it : p_doubles_ )
         if ( it.second.value )
           f << it.first << " = " << *it.second.value << "\n";
       for ( const auto& it : p_bools_ )
         if ( it.second.value )
           f << it.first << " = " << *it.second.value << "\n";
       f.close();
     }
 
     void
     LpairHandler::setParameter( const std::string& key, const std::string& value )
     {
       try { setValue<double>( key.c_str(), std::stod( value ) ); } catch ( std::invalid_argument& ) {}
       try { setValue<int>( key.c_str(), std::stoi( value ) ); } catch ( std::invalid_argument& ) {}
       //setValue<bool>( key.c_str(), std::stoi( value ) );
       setValue<std::string>( key.c_str(), value );
     }
 
     std::string
     LpairHandler::getParameter( std::string key ) const
     {
       double dd = getValue<double>( key.c_str() );
       if ( dd != -999. )
         return std::to_string( dd );
 
       int ui = getValue<int>( key.c_str() );
       if ( ui != 999 )
         return std::to_string( ui );
 
       //if ( out = getValue<bool>( key.c_str() )  );
 
       return getValue<std::string>( key.c_str() );
     }
 
     std::string
     LpairHandler::getDescription( std::string key ) const
     {
       if ( p_strings_.count( key ) )
         return p_strings_.find( key )->second.description;
       if ( p_ints_.count( key ) )
         return p_ints_.find( key )->second.description;
       if ( p_doubles_.count( key ) )
         return p_doubles_.find( key )->second.description;
       if ( p_bools_.count( key ) )
         return p_bools_.find( key )->second.description;
       return "null";
     }
   }
 }
 
diff --git a/CepGen/Cards/LpairHandler.h b/CepGen/Cards/LpairHandler.h
index 78cf708..58dc8ea 100644
--- a/CepGen/Cards/LpairHandler.h
+++ b/CepGen/Cards/LpairHandler.h
@@ -1,110 +1,111 @@
 #ifndef CepGen_Cards_LpairReader_h
 #define CepGen_Cards_LpairReader_h
 
 #include "CepGen/Cards/Handler.h"
 #include <unordered_map>
 
 using std::string;
 
 namespace CepGen
 {
   class ParametersList;
   namespace Cards
   {
     /// LPAIR-like steering cards parser and writer
     class LpairHandler : public Handler
     {
       public:
         /// Read a LPAIR steering card
         explicit LpairHandler( const char* file );
 
         /// Store a configuration into a LPAIR steering card
         void store( const char* file );
 
       private:
         template<class T> struct Parameter {
           Parameter( const char* key, const char* descr, T* value ) : key( key ), description( descr ), value( value ) {}
           std::string key, description;
           T* value;
         };
         /// Register a parameter to be steered to a configuration variable
         template<class T> void registerParameter( const char* key, const char* description, T* def ) {}
         /// Set a parameter value
         template<class T> void setValue( const char* key, const T& value ) {}
         /// Retrieve a parameter value
         template<class T> T getValue( const char* key ) const {}
 
         void setParameter( const std::string& key, const std::string& value );
         std::string getParameter( std::string key ) const;
         std::string getDescription( std::string key ) const;
 
         static const int kInvalid;
 
         std::unordered_map<std::string, Parameter<std::string> > p_strings_;
         std::unordered_map<std::string, Parameter<double> > p_doubles_;
         std::unordered_map<std::string, Parameter<int> > p_ints_;
         std::unordered_map<std::string, Parameter<bool> > p_bools_;
 
         void init( Parameters* );
         std::shared_ptr<ParametersList> proc_params_;
+        int str_fun_;
         std::string proc_name_, hadr_name_, integr_type_;
         std::pair<unsigned short,unsigned short> hi_1_, hi_2_;
     };
 
     //----- specialised registerers
 
     template<> inline void LpairHandler::registerParameter<std::string>( const char* key, const char* description, std::string* def ) { p_strings_.insert( std::make_pair( key, Parameter<std::string>( key, description, def ) ) ); }
     template<> inline void LpairHandler::registerParameter<double>( const char* key, const char* description, double* def ) { p_doubles_.insert( std::make_pair( key, Parameter<double>( key, description, def ) ) ); }
     template<> inline void LpairHandler::registerParameter<int>( const char* key, const char* description, int* def ) { p_ints_.insert( std::make_pair( key, Parameter<int>( key, description, def ) ) ); }
     template<> inline void LpairHandler::registerParameter<bool>( const char* key, const char* description, bool* def ) { p_bools_.insert( std::make_pair( key, Parameter<bool>( key, description, def ) ) ); }
 
     //----- specialised setters
 
     template<> inline void LpairHandler::setValue<std::string>( const char* key, const std::string& value ) {
       auto it = p_strings_.find( key );
       if ( it != p_strings_.end() ) *it->second.value = value;
     }
     template<> inline void LpairHandler::setValue<double>( const char* key, const double& value ) {
       auto it = p_doubles_.find( key );
       if ( it != p_doubles_.end() ) *it->second.value = value;
     }
     template<> inline void LpairHandler::setValue<int>( const char* key, const int& value ) {
       auto it = p_ints_.find( key );
       if ( it != p_ints_.end() ) *it->second.value = value;
     }
     template<> inline void LpairHandler::setValue<bool>( const char* key, const bool& value ) {
       auto it = p_bools_.find( key );
       if ( it != p_bools_.end() ) *it->second.value = value;
     }
 
     //----- specialised getters
 
     /// Retrieve a string parameter value
     template<> inline std::string LpairHandler::getValue( const char* key ) const {
       const auto& it = p_strings_.find( key );
       if ( it != p_strings_.end() ) return *it->second.value;
       return "null";
     }
     /// Retrieve a floating point parameter value
     template<> inline double LpairHandler::getValue( const char* key ) const {
       const auto& it = p_doubles_.find( key );
       if ( it != p_doubles_.end() ) return *it->second.value;
       return -999.;
     }
     /// Retrieve an integer parameter value
     template<> inline int LpairHandler::getValue( const char* key ) const {
       const auto& it = p_ints_.find( key );
       if ( it != p_ints_.end() ) return *it->second.value;
       return 999;
     }
     /// Retrieve a boolean parameter value
     template<> inline bool LpairHandler::getValue( const char* key ) const {
       const auto& it = p_bools_.find( key );
       if ( it != p_bools_.end() ) return *it->second.value;
       return true;
     }
   }
 }
 
 #endif
 
diff --git a/CepGen/Core/ParametersList.cpp b/CepGen/Core/ParametersList.cpp
index f15c953..f94cc67 100644
--- a/CepGen/Core/ParametersList.cpp
+++ b/CepGen/Core/ParametersList.cpp
@@ -1,298 +1,298 @@
 #include "CepGen/Core/ParametersList.h"
 #include "CepGen/Core/Exception.h"
 
 namespace CepGen
 {
   ParametersList&
   ParametersList::operator+=( const ParametersList& oth )
   {
     param_values_.insert( oth.param_values_.begin(), oth.param_values_.end() );
     int_values_.insert( oth.int_values_.begin(), oth.int_values_.end() );
     dbl_values_.insert( oth.dbl_values_.begin(), oth.dbl_values_.end() );
     str_values_.insert( oth.str_values_.begin(), oth.str_values_.end() );
     vec_param_values_.insert( oth.vec_param_values_.begin(), oth.vec_param_values_.end() );
     vec_int_values_.insert( oth.vec_int_values_.begin(), oth.vec_int_values_.end() );
     vec_dbl_values_.insert( oth.vec_dbl_values_.begin(), oth.vec_dbl_values_.end() );
     vec_str_values_.insert( oth.vec_str_values_.begin(), oth.vec_str_values_.end() );
     return *this;
   }
 
   std::ostream&
   operator<<( std::ostream& os, const ParametersList& params )
   {
     for ( const auto& kv : params.int_values_ )
       os << "\n" << kv.first << ": int(" << kv.second << ")";
     for ( const auto& kv : params.dbl_values_ )
       os << "\n" << kv.first << ": double(" << kv.second << ")";
     for ( const auto& kv : params.str_values_ )
       os << "\n" << kv.first << ": string(" << kv.second << ")";
     for ( const auto& kv : params.param_values_ )
       os << "\n" << kv.first << ": param({" << kv.second << "\n})";
     for ( const auto& kv : params.vec_int_values_ ) {
       os << "\n" << kv.first << ": vint(";
       bool first = true;
       for ( const auto& v : kv.second ) {
         os << ( first ? "" : ", " ) << v;
         first = false;
       }
       os << ")";
     }
     for ( const auto& kv : params.vec_dbl_values_ ) {
       os << "\n" << kv.first << ": vdouble(";
       bool first = true;
       for ( const auto& v : kv.second ) {
         os << ( first ? "" : ", " ) << v;
         first = false;
       }
       os << ")";
     }
     for ( const auto& kv : params.vec_str_values_ ) {
       os << "\n" << kv.first << ": vstring(";
       bool first = true;
       for ( const auto& v : kv.second ) {
         os << ( first ? "" : ", " ) << v;
         first = false;
       }
       os << ")";
     }
     return os;
   }
 
   //------------------------------------------------------------------
   // default template (placeholders)
   //------------------------------------------------------------------
 
   template<typename T> T
   ParametersList::get( std::string key, T def ) const
   {
     throw CG_FATAL( "ParametersList" ) << "Invalid type retrieved for key=" << key << "!";
   }
 
   template<typename T> T&
   ParametersList::operator[]( std::string key )
   {
     throw CG_FATAL( "ParametersList" ) << "Invalid type retrieved for key=" << key << "!";
   }
 
   template<typename T> void
   ParametersList::set( std::string key, const T& value )
   {
     throw CG_FATAL( "ParametersList" ) << "Invalid type to be set for key=" << key << "!";
   }
 
   //------------------------------------------------------------------
   // sub-parameters-type attributes
   //------------------------------------------------------------------
 
   template<> ParametersList
   ParametersList::get<ParametersList>( std::string key, ParametersList def ) const
   {
     for ( const auto& kv : param_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
     CG_DEBUG( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
     return def;
   }
 
   template<> ParametersList&
   ParametersList::operator[]<ParametersList>( std::string key )
   {
     for ( auto& kv : param_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
-    throw CG_FATAL( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
+    return param_values_[key];
   }
 
   template<> void
   ParametersList::set<ParametersList>( std::string key, const ParametersList& value )
   {
     param_values_[key] = value;
   }
 
   template<> std::vector<ParametersList>
   ParametersList::get<std::vector<ParametersList> >( std::string key, std::vector<ParametersList> def ) const
   {
     for ( const auto& kv : vec_param_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
     CG_DEBUG( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
     return def;
   }
 
   template<> std::vector<ParametersList>&
   ParametersList::operator[]<std::vector<ParametersList> >( std::string key )
   {
     for ( auto& kv : vec_param_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
-    throw CG_FATAL( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
+    return vec_param_values_[key];
   }
 
   template<> void
   ParametersList::set<std::vector<ParametersList> >( std::string key, const std::vector<ParametersList>& value )
   {
     vec_param_values_[key] = value;
   }
 
   //------------------------------------------------------------------
   // integer-type attributes
   //------------------------------------------------------------------
 
   template<> int
   ParametersList::get<int>( std::string key, int def ) const
   {
     for ( const auto& kv : int_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
     CG_DEBUG( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
     return def;
   }
 
   template<> int&
   ParametersList::operator[]<int>( std::string key )
   {
     for ( auto& kv : int_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
-    throw CG_FATAL( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
+    return int_values_[key];
   }
 
   template<> void
   ParametersList::set<int>( std::string key, const int& value )
   {
     int_values_[key] = value;
   }
 
   template<> std::vector<int>
   ParametersList::get<std::vector<int> >( std::string key, std::vector<int> def ) const
   {
     for ( const auto& kv : vec_int_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
     CG_DEBUG( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
     return def;
   }
 
   template<> std::vector<int>&
   ParametersList::operator[]<std::vector<int> >( std::string key )
   {
     for ( auto& kv : vec_int_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
-    throw CG_FATAL( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
+    return vec_int_values_[key];
   }
 
   template<> void
   ParametersList::set<std::vector<int> >( std::string key, const std::vector<int>& value )
   {
     vec_int_values_[key] = value;
   }
 
   //------------------------------------------------------------------
   // floating point-type attributes
   //------------------------------------------------------------------
 
   template<> double
   ParametersList::get<double>( std::string key, double def ) const
   {
     for ( const auto& kv : dbl_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
     CG_DEBUG( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
     return def;
   }
 
   template<> double&
   ParametersList::operator[]<double>( std::string key )
   {
     for ( auto& kv : dbl_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
-    throw CG_FATAL( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
+    return dbl_values_[key];
   }
 
   template<> void
   ParametersList::set<double>( std::string key, const double& value )
   {
     dbl_values_[key] = value;
   }
 
   template<> std::vector<double>
   ParametersList::get<std::vector<double> >( std::string key, std::vector<double> def ) const
   {
     for ( const auto& kv : vec_dbl_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
     CG_DEBUG( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
     return def;
   }
 
   template<> std::vector<double>&
   ParametersList::operator[]<std::vector<double> >( std::string key )
   {
     for ( auto& kv : vec_dbl_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
-    throw CG_FATAL( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
+    return vec_dbl_values_[key];
   }
 
   template<> void
   ParametersList::set<std::vector<double> >( std::string key, const std::vector<double>& value )
   {
     vec_dbl_values_[key] = value;
   }
 
   //------------------------------------------------------------------
   // string-type attributes
   //------------------------------------------------------------------
 
   template<> std::string
   ParametersList::get<std::string>( std::string key, std::string def ) const
   {
     for ( const auto& kv : str_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
     CG_DEBUG( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
     return def;
   }
 
   template<> std::string&
   ParametersList::operator[]<std::string>( std::string key )
   {
     for ( auto& kv : str_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
-    throw CG_FATAL( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
+    return str_values_[key];
   }
 
   template<> void
   ParametersList::set<std::string>( std::string key, const std::string& value )
   {
     str_values_[key] = value;
   }
 
   template<> std::vector<std::string>
   ParametersList::get<std::vector<std::string> >( std::string key, std::vector<std::string> def ) const
   {
     for ( const auto& kv : vec_str_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
     CG_DEBUG( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
     return def;
   }
 
   template<> std::vector<std::string>&
   ParametersList::operator[]<std::vector<std::string> >( std::string key )
   {
     for ( auto& kv : vec_str_values_ )
       if ( kv.first.compare( key ) == 0 )
         return kv.second;
-    throw CG_FATAL( "ParametersList" ) << "Failed to retrieve parameter with key=" << key << ".";
+    return vec_str_values_[key];
   }
 
   template<> void
   ParametersList::set<std::vector<std::string> >( std::string key, const std::vector<std::string>& value )
   {
     vec_str_values_[key] = value;
   }
 }
diff --git a/CepGen/StructureFunctions/LHAPDF.cpp b/CepGen/StructureFunctions/LHAPDF.cpp
index 9c983e4..281036c 100644
--- a/CepGen/StructureFunctions/LHAPDF.cpp
+++ b/CepGen/StructureFunctions/LHAPDF.cpp
@@ -1,150 +1,159 @@
 #include "CepGen/StructureFunctions/LHAPDF.h"
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/utils.h"
 
 namespace CepGen
 {
   namespace SF
   {
     constexpr std::array<short,6> LHAPDF::qtimes3_, LHAPDF::pdgid_;
 
     LHAPDF::Parameterisation::Parameterisation() :
-      num_flavours( 4 ), pdf_set( "cteq6" ), pdf_member( 0 ), mode( Mode::full )
+      num_flavours( 4 ), pdf_set( "cteq6" ), pdf_code( 0l ), pdf_member( 0 ), mode( Mode::full )
     {}
 
     LHAPDF::Parameterisation
     LHAPDF::Parameterisation::cteq6()
     {
       Parameterisation p;
       p.num_flavours = 4;
       p.pdf_set = "cteq6";
       p.pdf_member = 0;
       p.mode = Mode::full;
       return p;
     }
 
     LHAPDF::LHAPDF( const Parameterisation& param ) :
       StructureFunctions( Type::LHAPDF ), params( param ), initialised_( false )
     {}
 
     LHAPDF::LHAPDF( const char* set, unsigned short member, const Parameterisation::Mode& mode ) :
       StructureFunctions( Type::LHAPDF ), initialised_( false )
     {
       params.pdf_set = set;
       params.pdf_member = member;
       params.mode = mode;
     }
 
     void
     LHAPDF::initialise()
     {
       if ( initialised_ )
         return;
 #ifdef LIBLHAPDF
       std::string lhapdf_version, pdf_description, pdf_type;
 #  if defined LHAPDF_MAJOR_VERSION && LHAPDF_MAJOR_VERSION == 6
       try {
+        //--- check if PDF code is set
+        if ( params.pdf_code != 0l ) {
+          auto pdf = ::LHAPDF::lookupPDF( params.pdf_code );
+          if ( pdf.second != 0 )
+            throw CG_FATAL( "LHAPDF" ) << "Failed to retrieve PDFset with id=" << params.pdf_code << "!";
+          if ( !params.pdf_set.empty() && params.pdf_set != pdf.first )
+            CG_WARNING( "LHAPDF" ) << "PDF set name changed from \"" << params.pdf_set << "\" to \"" << pdf.first << "\".";
+          params.pdf_set = pdf.first;
+        }
         pdf_set_ = ::LHAPDF::PDFSet( params.pdf_set );
         pdf_set_.mkPDFs<std::unique_ptr<::LHAPDF::PDF> >( pdfs_ );
         lhapdf_version = ::LHAPDF::version();
         pdf_description = pdf_set_.description();
         pdf_type = pdfs_[params.pdf_member]->type();
       } catch ( const ::LHAPDF::Exception& e ) {
         throw CG_FATAL( "LHAPDF" )
           << "Caught LHAPDF exception:\n\t"
           << e.what();
       }
 #  else
       ::LHAPDF::initPDFSet( params.pdf_set, ::LHAPDF::LHGRID, params.pdf_member );
       lhapdf_version = ::LHAPDF::getVersion();
       pdf_description = ::LHAPDF::getDescription();
 #  endif
       replace_all( pdf_description, ". ", ".\n  " );
       CG_INFO( "LHAPDF" ) << "LHAPDF structure functions evaluator successfully built.\n"
         << " * LHAPDF version: " << lhapdf_version << "\n"
         << " * number of flavours: " << params.num_flavours << "\n"
         << " * PDF set: " << params.pdf_set << "\n"
         << ( pdf_description.empty() ? "" : "  "+pdf_description+"\n" )
         << " * PDF member: " << params.pdf_member << ( pdf_type.empty() ? "" : " ("+pdf_type+")" ) << "\n"
         << " * quarks mode: " << params.mode;
       initialised_ = true;
 #else
       throw CG_FATAL( "LHAPDF" ) << "LHAPDF is not liked to this instance!";
 #endif
     }
 
     LHAPDF&
     LHAPDF::operator()( double xbj, double q2 )
     {
 #ifdef LIBLHAPDF
       std::pair<double,double> nv = { xbj, q2 };
       if ( nv == old_vals_ )
         return *this;
       old_vals_ = nv;
 
       F2 = 0.;
       if ( params.num_flavours == 0 || params.num_flavours > 6 )
         return *this;
 
       if ( !initialised_ )
         initialise();
 #  if defined LHAPDF_MAJOR_VERSION && LHAPDF_MAJOR_VERSION >= 6
       auto& member = *pdfs_[params.pdf_member];
       if ( !member.inPhysicalRangeXQ2( xbj, q2 ) ) {
         CG_WARNING( "LHAPDF" ) << "(x=" << xbj << ", Q²=" << q2 << " GeV²) "
           << "not in physical range for PDF member " << params.pdf_member << ":\n\t"
           << "  min: (x=" << member.xMin() << ", Q²=" << member.q2Min() << "),\n\t"
           << "  max: (x=" << member.xMax() << ", Q²=" << member.q2Max() << ").";
         return *this;
       }
 #  else
       if ( q2 < ::LHAPDF::getQ2min( params.pdf_member ) || q2 > ::LHAPDF::getQ2max( params.pdf_member )
         || xbj < ::LHAPDF::getXmin( params.pdf_member ) || xbj > ::LHAPDF::getXmax( params.pdf_member ) ) {
         CG_WARNING( "LHAPDF" ) << "(x=" << xbj << "/Q²=" << q2 << " GeV²) "
           << "not in physical range for PDF member " << params.pdf_member << ":\n"
           << "  min: (x=" << ::LHAPDF::getXmin( params.pdf_member ) << "/Q²=" << ::LHAPDF::getQ2min( params.pdf_member ) << "),\n"
           << "  max: (x=" << ::LHAPDF::getXmax( params.pdf_member ) << "/Q²=" << ::LHAPDF::getQ2max( params.pdf_member ) << ").";
         return *this;
       }
       const double q = sqrt( q2 );
 #  endif
 
       for ( int i = 0; i < params.num_flavours; ++i ) {
         const double prefactor = 1./9.*qtimes3_[i]*qtimes3_[i];
 #  if defined LHAPDF_MAJOR_VERSION && LHAPDF_MAJOR_VERSION >= 6
         if ( !pdfs_[params.pdf_member]->hasFlavor( pdgid_[i] ) )
           throw CG_FATAL( "LHAPDF" ) << "Flavour " << pdgid_[i] << " is unsupported!";
         const double xq = member.xfxQ2( pdgid_[i], xbj, q2 );
         const double xqbar = member.xfxQ2( -pdgid_[i], xbj, q2 );
 #  else
         const double xq = ::LHAPDF::xfx( xbj, q, pdgid_[i] );
         const double xqbar = ::LHAPDF::xfx( xbj, q, -pdgid_[i] );
 #  endif
         switch ( params.mode ) {
           case Parameterisation::Mode::full:
             F2 += prefactor*( xq+xqbar ); break;
           case Parameterisation::Mode::valence:
             F2 += prefactor*( xq-xqbar ); break;
           case Parameterisation::Mode::sea:
             F2 += prefactor*( 2.*xqbar ); break;
         }
       }
 #else
       throw CG_FATAL( "LHAPDF" ) << "LHAPDF is not liked to this instance!";
 #endif
 
       return *this;
     }
   }
 
   std::ostream&
   operator<<( std::ostream& os, const SF::LHAPDF::Parameterisation::Mode& mode )
   {
     switch ( mode ) {
       case SF::LHAPDF::Parameterisation::Mode::full: return os << "all quarks";
       case SF::LHAPDF::Parameterisation::Mode::valence: return os << "valence quarks";
       case SF::LHAPDF::Parameterisation::Mode::sea: return os << "sea quarks";
     }
     return os;
   }
 }
diff --git a/CepGen/StructureFunctions/LHAPDF.h b/CepGen/StructureFunctions/LHAPDF.h
index a4e264e..f084ab7 100644
--- a/CepGen/StructureFunctions/LHAPDF.h
+++ b/CepGen/StructureFunctions/LHAPDF.h
@@ -1,58 +1,59 @@
 #ifndef CepGen_StructureFunctions_LHAPDF_h
 #define CepGen_StructureFunctions_LHAPDF_h
 
 #include "StructureFunctions.h"
 
 #ifdef LIBLHAPDF
 #include "LHAPDF/LHAPDF.h"
 #endif
 
 #include <array>
 
 namespace CepGen
 {
   namespace SF
   {
     /// Generic, tree-level import of structure functions from an external PDFs grid
     class LHAPDF : public StructureFunctions
     {
       public:
         struct Parameterisation
         {
           Parameterisation();
           static Parameterisation cteq6();
           unsigned short num_flavours;
           std::string pdf_set;
+          unsigned long pdf_code;
           unsigned short pdf_member;
           enum class Mode { full = 0, valence = 1, sea = 2 };
           Mode mode;
         };
 
         explicit LHAPDF( const Parameterisation& param = Parameterisation::cteq6() );
         explicit LHAPDF( const char* set, unsigned short member = 0, const Parameterisation::Mode& mode = Parameterisation::Mode::full );
         LHAPDF& operator()( double xbj, double q2 ) override;
 
         Parameterisation params;
 
       private:
         void initialise();
         bool initialised_;
 
 #ifdef LIBLHAPDF
 #  if defined LHAPDF_MAJOR_VERSION && LHAPDF_MAJOR_VERSION >= 6
         ::LHAPDF::PDFSet pdf_set_;
         std::vector<std::unique_ptr<::LHAPDF::PDF> > pdfs_;
 #  endif
 #endif
         static constexpr std::array<short,6> pdgid_ = { { 1, 2, 3, 4, 5, 6 } };
         static constexpr std::array<short,6> qtimes3_ = { {
           -1 /*d*/, 2 /*u*/,
           -1 /*s*/, 2 /*c*/,
           -1 /*b*/, 2 /*t*/
         } };
     };
   }
   std::ostream& operator<<( std::ostream& os, const SF::LHAPDF::Parameterisation::Mode& mode );
 }
 
 #endif