diff --git a/CepGen/Core/Parameters.cpp b/CepGen/Core/Parameters.cpp
index 886c597..f7a5880 100644
--- a/CepGen/Core/Parameters.cpp
+++ b/CepGen/Core/Parameters.cpp
@@ -1,269 +1,269 @@
 #include "CepGen/Parameters.h"
 
 #include "CepGen/Core/ParametersList.h"
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/TamingFunction.h"
 
 #include "CepGen/Physics/PDG.h"
 #include "CepGen/Processes/GenericProcess.h"
 #include "CepGen/Hadronisers/GenericHadroniser.h"
 
 #include "CepGen/StructureFunctions/StructureFunctions.h"
 
 #include <iomanip>
 
 namespace CepGen
 {
   Parameters::Parameters() :
     general( new ParametersList ),
     hadroniser_max_trials( 5 ),
     taming_functions( new TamingFunctionsCollection ),
     store_( false ), total_gen_time_( 0. ), num_gen_events_( 0 )
   {}
 
   Parameters::Parameters( Parameters& param ) :
     general( param.general ),
     kinematics( param.kinematics ), integrator( param.integrator ), generation( param.generation ),
     hadroniser_max_trials( param.hadroniser_max_trials ),
     taming_functions( param.taming_functions ),
     process_( std::move( param.process_ ) ),
     hadroniser_( std::move( param.hadroniser_ ) ),
     store_( false ), total_gen_time_( param.total_gen_time_ ), num_gen_events_( param.num_gen_events_ )
   {}
 
   Parameters::Parameters( const Parameters& param ) :
     general( param.general ),
     kinematics( param.kinematics ), integrator( param.integrator ), generation( param.generation ),
     hadroniser_max_trials( param.hadroniser_max_trials ),
     taming_functions( param.taming_functions ),
     store_( false ), total_gen_time_( param.total_gen_time_ ), num_gen_events_( param.num_gen_events_ )
   {}
 
   Parameters::~Parameters() // required for unique_ptr initialisation!
   {}
 
   void
   Parameters::setThetaRange( float thetamin, float thetamax )
   {
     kinematics.cuts.central.eta_single = {
       Particle::thetaToEta( thetamax ),
       Particle::thetaToEta( thetamin )
     };
 
     CG_DEBUG( "Parameters" )
       << "eta in range: " << kinematics.cuts.central.eta_single
       << " => theta(min) = " << thetamin << ", theta(max) = " << thetamax << ".";
   }
 
   void
   Parameters::clearRunStatistics()
   {
     total_gen_time_ = 0.;
     num_gen_events_ = 0;
   }
 
   void
   Parameters::addGenerationTime( double gen_time )
   {
     total_gen_time_ += gen_time;
     num_gen_events_++;
   }
 
   Process::GenericProcess*
   Parameters::process()
   {
     return process_.get();
   }
 
   std::string
   Parameters::processName() const
   {
     if ( !process_ )
       return "no process";
     return process_->name();
   }
 
   void
   Parameters::setProcess( Process::GenericProcess* proc )
   {
     if ( !proc )
       throw CG_FATAL( "Parameters" )
         << "Trying to clone an invalid process!";
     process_.reset( proc );
   }
 
   void
   Parameters::cloneProcess( const Process::GenericProcess* proc )
   {
     if ( !proc )
       throw CG_FATAL( "Parameters" )
         << "Trying to clone an invalid process!";
     process_ = std::move( proc->clone() );
   }
 
   Hadroniser::GenericHadroniser*
   Parameters::hadroniser()
   {
     return hadroniser_.get();
   }
 
   std::string
   Parameters::hadroniserName() const
   {
     if ( !hadroniser_ )
       return "";
     return hadroniser_->name();
   }
 
   void
   Parameters::setHadroniser( Hadroniser::GenericHadroniser* hadr )
   {
     hadroniser_.reset( hadr );
   }
 
   std::ostream&
   operator<<( std::ostream& os, const Parameters* p )
   {
     const bool pretty = true;
 
     const int wb = 90, wt = 40;
     os
       << "Parameters dump" << std::left << "\n\n"
       << std::setfill('_') << std::setw( wb+3 ) << "_/¯¯RUN¯INFORMATION¯¯\\_" << std::setfill( ' ' ) << "\n"
       << std::right << std::setw( wb ) << std::left << std::endl
       << std::setw( wt ) << "Process to generate";
     if ( p->process_ ) {
       os << ( pretty ? boldify( p->process_->name().c_str() ) : p->process_->name() ) << "\n"
          << std::setw( wt ) << "" << p->process_->description();
     }
     else
       os << ( pretty ? boldify( "no process!" ) : "no process!" );
     os
       << "\n"
       << std::setw( wt ) << "Events generation? "
       << ( pretty ? yesno( p->generation.enabled ) : std::to_string( p->generation.enabled ) ) << "\n"
       << std::setw( wt ) << "Number of events to generate"
       << ( pretty ? boldify( p->generation.maxgen ) : std::to_string( p->generation.maxgen ) ) << "\n";
     if ( p->generation.num_threads > 1 )
       os
         << std::setw( wt ) << "Number of threads" << p->generation.num_threads << "\n";
     os
       << std::setw( wt ) << "Number of points to try per bin" << p->generation.num_points << "\n"
       << std::setw( wt ) << "Integrand treatment"
       << ( pretty ? yesno( p->generation.treat ) : std::to_string( p->generation.treat ) ) << "\n"
       << std::setw( wt ) << "Verbosity level " << Logger::get().level << "\n";
     if ( p->hadroniser_ ) {
       os
         << "\n"
         << std::setfill( '-' ) << std::setw( wb+6 )
         << ( pretty ? boldify( " Hadronisation algorithm " ) : "Hadronisation algorithm" ) << std::setfill( ' ' ) << "\n\n"
         << std::setw( wt ) << "Name"
         << ( pretty ? boldify( p->hadroniser_->name().c_str() ) : p->hadroniser_->name() ) << "\n";
     }
     os
       << "\n"
       << std::setfill( '-' ) << std::setw( wb+6 )
       << ( pretty ? boldify( " Integration parameters " ) : "Integration parameters" ) << std::setfill( ' ' ) << "\n\n";
     std::ostringstream int_algo; int_algo << p->integrator.type;
     os
       << std::setw( wt ) << "Integration algorithm"
       << ( pretty ? boldify( int_algo.str().c_str() ) : int_algo.str() ) << "\n"
       << std::setw( wt ) << "Number of function calls" << p->integrator.ncvg << "\n"
       << std::setw( wt ) << "Random number generator seed" << p->integrator.rng_seed << "\n";
     if ( p->integrator.rng_engine )
       os
         << std::setw( wt ) << "Random number generator engine"
         << p->integrator.rng_engine->name << "\n";
     std::ostringstream proc_mode; proc_mode << p->kinematics.mode;
     os
       << "\n"
       << std::setfill('_') << std::setw( wb+3 )
       << "_/¯¯EVENTS¯KINEMATICS¯¯\\_" << std::setfill( ' ' ) << "\n\n"
       << std::setw( wt ) << "Incoming particles"
       << p->kinematics.incoming_beams.first << ",\n" << std::setw( wt ) << ""
       << p->kinematics.incoming_beams.second << "\n";
     if ( p->kinematics.mode != KinematicsMode::invalid )
       os << std::setw( wt ) << "Subprocess mode" << ( pretty ? boldify( proc_mode.str().c_str() ) : proc_mode.str() ) << "\n";
     if ( p->kinematics.mode != KinematicsMode::ElasticElastic )
-      os << std::setw( wt ) << "Structure functions" << p->kinematics.structure_functions->type << "\n";
+      os << std::setw( wt ) << "Structure functions" << *p->kinematics.structure_functions << "\n";
     os
       << "\n"
       << std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Incoming partons " ) : "Incoming partons" ) << std::setfill( ' ' ) << "\n\n";
     for ( const auto& lim : p->kinematics.cuts.initial.list() ) { // map(particles class, limits)
       if ( !lim.second.valid() )
         continue;
       os << std::setw( wt ) << lim.first << lim.second << "\n";
     }
     os
       << "\n"
       << std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Outgoing central system " ) : "Outgoing central system" ) << std::setfill( ' ' ) << "\n\n";
     for ( const auto& lim : p->kinematics.cuts.central.list() ) {
       if ( !lim.second.valid() )
         continue;
       os << std::setw( wt ) << lim.first << lim.second << "\n";
     }
     if ( p->kinematics.cuts.central_particles.size() > 0 ) {
       os << std::setw( wt ) << ( pretty ? boldify( ">>> per-particle cuts:" ) : ">>> per-particle cuts:" ) << "\n";
       for ( const auto& part_per_lim : p->kinematics.cuts.central_particles ) {
         os << " * all single " << std::setw( wt-3 ) << part_per_lim.first << "\n";
         for ( const auto& lim : part_per_lim.second.list() ) {
           if ( !lim.second.valid() )
             continue;
           os << "   - " << std::setw( wt-5 ) << lim.first << lim.second << "\n";
         }
       }
     }
     os << "\n";
     os << std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Proton / remnants " ) : "Proton / remnants" ) << std::setfill( ' ' ) << "\n\n";
     for ( const auto& lim : p->kinematics.cuts.remnants.list() )
       os << std::setw( wt ) << lim.first << lim.second << "\n";
     return os;
   }
 
   std::ostream&
   operator<<( std::ostream& os, const Parameters& p )
   {
     return os << &p;
   }
 
   //-----------------------------------------------------------------------------------------------
 
   Parameters::Integration::Integration() :
     type( Integrator::Type::Vegas ), ncvg( 500000 ),
     rng_seed( 0 ), rng_engine( (gsl_rng_type*)gsl_rng_mt19937 ),
     vegas_chisq_cut( 1.5 ),
     result( -1. ), err_result( -1. )
   {
     const size_t ndof = 10;
     {
       std::shared_ptr<gsl_monte_vegas_state> tmp_state( gsl_monte_vegas_alloc( ndof ), gsl_monte_vegas_free );
       gsl_monte_vegas_params_get( tmp_state.get(), &vegas );
     }
     {
       std::shared_ptr<gsl_monte_miser_state> tmp_state( gsl_monte_miser_alloc( ndof ), gsl_monte_miser_free );
       gsl_monte_miser_params_get( tmp_state.get(), &miser );
     }
   }
 
   Parameters::Integration::Integration( const Integration& rhs ) :
     type( rhs.type ), ncvg( rhs.ncvg ),
     rng_seed( rhs.rng_seed ), rng_engine( rhs.rng_engine ),
     vegas( rhs.vegas ), vegas_chisq_cut( rhs.vegas_chisq_cut ),
     miser( rhs.miser ),
     result( -1. ), err_result( -1. )
   {}
 
   Parameters::Integration::~Integration()
   {
     //if ( vegas.ostream && vegas.ostream != stdout && vegas.ostream != stderr )
     //  fclose( vegas.ostream );
   }
 
   //-----------------------------------------------------------------------------------------------
 
   Parameters::Generation::Generation() :
     enabled( false ), maxgen( 0 ),
     symmetrise( false ), treat( false ), ngen( 0 ), gen_print_every( 10000 ),
     num_threads( 2 ), num_points( 100 )
   {}
 }
diff --git a/CepGen/StructureFunctions/LHAPDF.cpp b/CepGen/StructureFunctions/LHAPDF.cpp
index f60d196..042ff83 100644
--- a/CepGen/StructureFunctions/LHAPDF.cpp
+++ b/CepGen/StructureFunctions/LHAPDF.cpp
@@ -1,162 +1,170 @@
 #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_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;
     }
 
+    std::string
+    LHAPDF::description() const
+    {
+      std::ostringstream os;
+      os << "LHAPDF{" << params.pdf_set << ",m=" << params.pdf_member << ",mode=" << params.mode << "}";
+      return os.str();
+    }
+
     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
       if ( params.pdf_code != 0l )
         ::LHAPDF::initPDFSet( (int)params.pdf_code, params.pdf_member );
       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 f084ab7..a71cc01 100644
--- a/CepGen/StructureFunctions/LHAPDF.h
+++ b/CepGen/StructureFunctions/LHAPDF.h
@@ -1,59 +1,60 @@
 #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:
+        std::string description() const override;
         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
diff --git a/CepGen/StructureFunctions/Schaefer.cpp b/CepGen/StructureFunctions/Schaefer.cpp
index a56e8b3..a34a2e7 100644
--- a/CepGen/StructureFunctions/Schaefer.cpp
+++ b/CepGen/StructureFunctions/Schaefer.cpp
@@ -1,137 +1,151 @@
 #include "CepGen/StructureFunctions/Schaefer.h"
 #include "CepGen/StructureFunctions/StructureFunctionsBuilder.h"
 #include "CepGen/StructureFunctions/LHAPDF.h"
 
 #include "CepGen/Core/Exception.h"
 
 #include "CepGen/Physics/Constants.h"
 #include "CepGen/Physics/ParticleProperties.h"
 
 namespace CepGen
 {
   namespace SF
   {
     Schaefer::Parameterisation
     Schaefer::Parameterisation::mstwGrid()
     {
       Parameterisation par;
       par.q2_cut = 9.;
       par.w2_hi = 4.;
       par.w2_lo = 3.;
       par.resonances_model = StructureFunctionsBuilder::get( SF::Type::ChristyBosted );
       par.perturbative_model = StructureFunctionsBuilder::get( SF::Type::MSTWgrid );
       par.continuum_model = StructureFunctionsBuilder::get( SF::Type::GD11p );
       par.higher_twist = 0;
       return par;
     }
 
     Schaefer::Parameterisation
     Schaefer::Parameterisation::mstwParton()
     {
       Parameterisation par;
       par.q2_cut = 9.;
       par.w2_hi = 4.;
       par.w2_lo = 3.;
       par.resonances_model = StructureFunctionsBuilder::get( SF::Type::ChristyBosted );
       par.perturbative_model = std::make_shared<SF::LHAPDF>( "MSTW2008nnlo90cl" );
       par.continuum_model = StructureFunctionsBuilder::get( SF::Type::GD11p );
       par.higher_twist = 1;
       return par;
     }
 
     Schaefer::Parameterisation
     Schaefer::Parameterisation::cteq()
     {
       Parameterisation par;
       par.q2_cut = 9.;
       par.w2_hi = 4.;
       par.w2_lo = 3.;
       par.resonances_model = StructureFunctionsBuilder::get( SF::Type::ChristyBosted );
       par.perturbative_model = std::make_shared<SF::LHAPDF>( "cteq6l1" );
       par.continuum_model = StructureFunctionsBuilder::get( SF::Type::GD11p );
       par.higher_twist = 0;
       return par;
     }
 
     Schaefer::Schaefer( const Schaefer::Parameterisation& params ) :
       StructureFunctions( Type::Schaefer ),
       params( params ), initialised_( false ), inv_omega_range_( -1. )
     {}
 
+    std::string
+    Schaefer::description() const
+    {
+      std::ostringstream os;
+      os << "LUXlike{"
+         << "r=" << *params.resonances_model
+         << ",p=" << *params.perturbative_model
+         << ",c=" << *params.continuum_model;
+      if ( params.higher_twist )
+        os << ",HT";
+      os << "}";
+      return os.str();
+    }
+
     void
     Schaefer::initialise()
     {
       CG_INFO( "LUXlike" ) << "LUXlike structure functions evaluator successfully initialised.\n"
         << " * Q² cut:             " << params.q2_cut << " GeV²\n"
         << " * W² ranges:          " << params.w2_lo << " GeV² / " << params.w2_hi << " GeV²\n"
         << " * resonance model:    " << params.resonances_model->type << "\n"
         << " * perturbative model: " << params.perturbative_model->type << "\n"
         << " * continuum model:    " << params.continuum_model->type << "\n"
         << " * higher-twist?       " << std::boolalpha << params.higher_twist;
       inv_omega_range_ = 1./( params.w2_hi-params.w2_lo );
       initialised_ = true;
     }
 
     Schaefer&
     Schaefer::operator()( double xbj, double q2 )
     {
       if ( !initialised_ )
         initialise();
 
       std::pair<double,double> nv = { xbj, q2 };
       if ( nv == old_vals_ )
         return *this;
       old_vals_ = nv;
 
       const double w2 = mp2_+q2*( 1.-xbj )/xbj;
 
       StructureFunctions sel_sf;
       if ( q2 < params.q2_cut ) {
         if ( w2 < params.w2_lo )
           sel_sf = ( *params.resonances_model )( xbj, q2 );
         else if ( w2 < params.w2_hi ) {
           auto sf_r = ( *params.resonances_model )( xbj, q2 );
           auto sf_c = ( *params.continuum_model )( xbj, q2 );
           sf_r.computeFL( xbj, q2 );
           sf_c.computeFL( xbj, q2 );
           const double r = rho( w2 );
           F2 = r*sf_c.F2 + ( 1.-r )*sf_r.F2;
           FL = r*sf_c.FL + ( 1.-r )*sf_r.FL;
           return *this;
         }
         else
           sel_sf = ( *params.continuum_model )( xbj, q2 );
       }
       else {
         if ( w2 < params.w2_hi )
           sel_sf = ( *params.continuum_model )( xbj, q2 );
         else {
           auto sf_p = ( *params.perturbative_model )( xbj, q2 );
           F2 = sf_p.F2;
           sf_p.computeFL( xbj, q2 );
           FL = sel_sf.FL;
           if ( params.higher_twist )
             F2 *= ( 1.+5.5/q2 );
           return *this;
         }
       }
 
       F2 = sel_sf( xbj, q2 ).F2;
       sel_sf.computeFL( xbj, q2 );
       FL = sel_sf.FL;
 
       return *this;
     }
 
     double
     Schaefer::rho( double w2 ) const
     {
       if ( inv_omega_range_ <= 0. )
         throw CG_FATAL( "LUXlike" ) << "Invalid W² limits: "
           << params.w2_lo << " / " << params.w2_hi << " GeV²!";
       const double omega = ( w2-params.w2_lo )*inv_omega_range_;
       const double omega2 = omega*omega;
       return 2.*omega2-omega*omega;
     }
   }
 }
diff --git a/CepGen/StructureFunctions/Schaefer.h b/CepGen/StructureFunctions/Schaefer.h
index 8ae678b..2f81150 100644
--- a/CepGen/StructureFunctions/Schaefer.h
+++ b/CepGen/StructureFunctions/Schaefer.h
@@ -1,37 +1,38 @@
 #ifndef CepGen_StructureFunctions_Schaefer_h
 #define CepGen_StructureFunctions_Schaefer_h
 
 #include "StructureFunctions.h"
 #include <memory>
 
 namespace CepGen
 {
   namespace SF
   {
     class Schaefer : public StructureFunctions
     {
       public:
         struct Parameterisation
         {
           static Parameterisation mstwGrid();
           static Parameterisation mstwParton();
           static Parameterisation cteq();
           double q2_cut, w2_lo, w2_hi;
           std::shared_ptr<StructureFunctions> resonances_model, perturbative_model, continuum_model;
           bool higher_twist;
         };
         Schaefer( const Parameterisation& param = Parameterisation::mstwGrid() );
         Schaefer& operator()( double xbj, double q2 ) override;
 
         Parameterisation params;
 
       private:
+        std::string description() const override;
         double rho( double w2 ) const;
         void initialise();
         bool initialised_;
         double inv_omega_range_;
     };
   }
 }
 
 #endif
diff --git a/CepGen/StructureFunctions/StructureFunctions.cpp b/CepGen/StructureFunctions/StructureFunctions.cpp
index f98c2ed..6d6c31c 100644
--- a/CepGen/StructureFunctions/StructureFunctions.cpp
+++ b/CepGen/StructureFunctions/StructureFunctions.cpp
@@ -1,75 +1,88 @@
 #include "CepGen/StructureFunctions/StructureFunctions.h"
 
 #include "CepGen/Physics/PDG.h"
 #include "CepGen/Physics/ParticleProperties.h"
 
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/utils.h"
 
 #include <iostream>
 
 namespace CepGen
 {
   const double StructureFunctions::mp_ = ParticleProperties::mass( PDG::proton );
   const double StructureFunctions::mp2_ = StructureFunctions::mp_*StructureFunctions::mp_;
 
   double
   StructureFunctions::F1( double xbj, double q2 ) const
   {
     if ( xbj == 0. || q2 == 0. ) {
       CG_ERROR( "StructureFunctions:F1" )
         << "Invalid range for Q² = " << q2 << " or xBj = " << xbj << ".";
       return 0.;
     }
     const double F1 = 0.5*( ( 1+4.*xbj*xbj*mp2_/q2 )*F2 - FL )/xbj;
     CG_DEBUG_LOOP( "StructureFunctions:F1" )
       << "F1 for Q² = " << q2 << ", xBj = " << xbj << ": " << F1 << "\n\t"
       << "(F2 = " << F2 << ", FL = " << FL << ").";
     return F1;
   }
 
   void
   StructureFunctions::computeFL( double xbj, double q2, const SF::SigmaRatio& ratio )
   {
     double r_error = 0.;
     computeFL( xbj, q2, ratio( xbj, q2, r_error ) );
   }
 
   void
   StructureFunctions::computeFL( double xbj, double q2, double r )
   {
     const double tau = 4.*xbj*xbj*mp2_/q2;
     FL = F2 * ( 1.+tau ) * ( r/( 1.+r ) );
   }
 
+  std::string
+  StructureFunctions::description() const
+  {
+    std::ostringstream os;
+    os << type;
+    return os.str();
+  }
+
   /// Human-readable format of a structure function object
   std::ostream&
   operator<<( std::ostream& os, const StructureFunctions& sf )
   {
-    return os << "F2 = " << sf.F2 << ", FL = " << sf.FL;
+    os << sf.description();
+    if ( sf.old_vals_ != std::pair<double,double>() )
+      os << " at (" << sf.old_vals_.first << ", " << sf.old_vals_.second << "): "
+         << "F2 = " << sf.F2 << ", FL = " << sf.FL;
+    return os;
   }
+
   /// Human-readable format of a structure function type
   std::ostream&
   operator<<( std::ostream& os, const SF::Type& sf )
   {
     switch ( sf ) {
       case SF::Type::Invalid:             return os << "[INVALID]";
       case SF::Type::Electron:            return os << "electron";
       case SF::Type::ElasticProton:       return os << "elastic proton";
       case SF::Type::SuriYennie:          return os << "Suri-Yennie";
       case SF::Type::SzczurekUleshchenko: return os << "Szczurek-Uleshchenko";
       case SF::Type::FioreBrasse:         return os << "Fiore-Brasse";
       case SF::Type::ChristyBosted:       return os << "Christy-Bosted";
       case SF::Type::CLAS:                return os << "CLAS";
       case SF::Type::BlockDurandHa:       return os << "BDH";
       case SF::Type::ALLM91:              return os << "ALLM91";
       case SF::Type::ALLM97:              return os << "ALLM97";
       case SF::Type::GD07p:               return os << "GD07p";
       case SF::Type::GD11p:               return os << "GD11p";
       case SF::Type::Schaefer:            return os << "LUXlike";
       case SF::Type::MSTWgrid:            return os << "MSTW (grid)";
       case SF::Type::LHAPDF:              return os << "LHAPDF";
     }
     return os;
   }
 }
diff --git a/CepGen/StructureFunctions/StructureFunctions.h b/CepGen/StructureFunctions/StructureFunctions.h
index 0d16122..c313fd7 100644
--- a/CepGen/StructureFunctions/StructureFunctions.h
+++ b/CepGen/StructureFunctions/StructureFunctions.h
@@ -1,71 +1,72 @@
 #ifndef CepGen_StructureFunctions_StructureFunctions_h
 #define CepGen_StructureFunctions_StructureFunctions_h
 
 #include "CepGen/StructureFunctions/SigmaRatio.h"
 
 #include <iostream>
 #include <map>
 
 namespace CepGen
 {
   namespace SF
   {
     /// Proton structure function to be used in the outgoing state description
     /// \note Values correspond to the LPAIR legacy steering card values
     enum struct Type {
       Invalid             = 0,
       Electron            = 1,
       ElasticProton       = 2,
       SuriYennie          = 11,
       SzczurekUleshchenko = 12,
       BlockDurandHa       = 13,
       FioreBrasse         = 101,
       ChristyBosted       = 102,
       CLAS                = 103,
       ALLM91              = 201,
       ALLM97              = 202,
       GD07p               = 203,
       GD11p               = 204,
       MSTWgrid            = 205,
       Schaefer            = 301,
       LHAPDF              = 401,
     };
   }
   std::ostream& operator<<( std::ostream&, const SF::Type& );
 
   class StructureFunctionsFactory;
   class StructureFunctions
   {
     public:
       StructureFunctions( const StructureFunctions& sf ) :
         type( sf.type ), F2( sf.F2 ), FL( sf.FL ), old_vals_( sf.old_vals_ ) {}
       StructureFunctions( const SF::Type& type = SF::Type::Invalid, double f2 = 0., double fl = 0. ) :
         type( type ), F2( f2 ), FL( fl ), old_vals_({ 0., 0. }) {}
       ~StructureFunctions() {}
 
+      friend std::ostream& operator<<( std::ostream&, const StructureFunctions& );
       StructureFunctions& operator=( const StructureFunctions& sf ) {
         type = sf.type, F2 = sf.F2, FL = sf.FL, old_vals_ = sf.old_vals_;
         return *this;
       }
 
       static StructureFunctions builder( const SF::Type& );
 
       virtual StructureFunctions& operator()( double xbj, double q2 ) { return *this; }
       virtual void computeFL( double xbj, double q2, const SF::SigmaRatio& ratio = SF::E143Ratio() );
       virtual void computeFL( double xbj, double q2, double r );
       double F1( double xbj, double q2 ) const;
 
       SF::Type type;
       double F2, FL;
 
     protected:
+      virtual std::string description() const;
       static const double mp_, mp2_;
       std::pair<double,double> old_vals_;
 
     private:
       std::string name_;
   };
-  std::ostream& operator<<( std::ostream&, const StructureFunctions& );
 }
 
 #endif