diff --git a/CMakeLists.txt b/CMakeLists.txt
index ae8341c..ed1d74b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,91 +1,92 @@
 cmake_minimum_required(VERSION 2.6 FATAL_ERROR)
 project(CepGen)
 set(PROJECT_VERSION 1)
 
 if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wno-long-long -pedantic-errors -std=c++11 -g")
 else()
   if(CMAKE_CXX_COMPILER_VERSION VERSION_GREATER 4.7)
     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wno-long-long -pedantic-errors -std=c++11 -g")
   else()
     message(FATAL_ERROR "gcc version >= 4.7 is required")
   endif()
 endif()
 set(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR} ${CMAKE_SOURCE_DIR}/cmake)
 set(CEPGEN_SOURCE_DIR ${PROJECT_SOURCE_DIR}/CepGen)
 set(CEPGEN_LIBRARIES CepGenCore CepGenEvent CepGenPhysics CepGenProcesses CepGenStructureFunctions CepGenExporter CepGenExternalHadronisers)
 
 #----- define all individual modules to be built beforehand
 
 set(CEPGEN_MODULES Processes Event Physics StructureFunctions Export Hadronisers)
 
 #----- enable fortran (for Pythia/Herwig/...)
 
 enable_language(Fortran)
 
 #----- include external paths
 
 include(UseEnvironment)
 set(CEPGEN_EXTERNAL_REQS "")
 if(GSL_LIB)
   message(STATUS "GSL found in ${GSL_LIB}")
   list(APPEND CEPGEN_EXTERNAL_REQS ${GSL_LIB} ${GSL_CBLAS_LIB})
 else()
   message(FATAL_ERROR "GSL not found on your system!")
 endif()
 if(MUPARSER)
   message(STATUS "muParser found in ${MUPARSER}")
   list(APPEND CEPGEN_EXTERNAL_REQS ${MUPARSER})
 endif()
 if(LIBCONFIG)
   message(STATUS "libconfig++ found in ${LIBCONFIG}")
   list(APPEND CEPGEN_EXTERNAL_REQS ${LIBCONFIG})
 endif()
 if(LHAPDF)
   message(STATUS "LHAPDF found in ${LHAPDF}")
   add_definitions(-DLIBLHAPDF)
   list(APPEND CEPGEN_EXTERNAL_REQS ${LHAPDF})
 endif()
 if(HEPMC_LIB)
   message(STATUS "HepMC found in ${HEPMC_INCLUDE}")
   include_directories(${HEPMC_INCLUDE})
-  list(APPEND CEPGEN_EXTERNAL_REQS ${HEPMC_LIB} ${HEPMC_FIO_LIB})
+  list(APPEND CEPGEN_EXTERNAL_REQS ${HEPMC_LIB})
+  #list(APPEND CEPGEN_EXTERNAL_REQS ${HEPMC_FIO_LIB})
   add_definitions(-DLIBHEPMC)
 endif()
 if(PYTHIA8)
   message(STATUS "Pythia8 found in ${PYTHIA8}")
   list(APPEND CEPGEN_EXTERNAL_REQS ${PYTHIA8})
 endif()
 
 include_directories(${LHAPDF_INCLUDE})
 include_directories(${PYTHIA8_INCLUDE})
 
 #----- build all the intermediate objects
 
 include_directories(${PROJECT_SOURCE_DIR})
 add_subdirectory(${CEPGEN_SOURCE_DIR})
 foreach(_module ${CEPGEN_MODULES})
   include_directories(${CEPGEN_SOURCE_DIR}/${_module})
   add_subdirectory(${CEPGEN_SOURCE_DIR}/${_module})
 endforeach()
 
 #----- copy the input cards and other files
 
 file(GLOB input_cards RELATIVE ${PROJECT_SOURCE_DIR} Cards/*)
 foreach(_files ${input_cards})
   configure_file(${_files} ${_files} COPYONLY)
 endforeach()
 configure_file(${CEPGEN_SOURCE_DIR}/README README COPYONLY)
 
 #----- installation rules
 
 set(MODS "")
 foreach(_module ${CEPGEN_MODULES})
   list(APPEND MODS CepGen/${module})
 endforeach()
 install(DIRECTORY ${MODS} DESTINATION include/CepGen FILES_MATCHING PATTERN "*.h")
 
 #----- set the tests/utils directory
 
 add_subdirectory(test)
 
diff --git a/CepGen/Cards/LpairHandler.cpp b/CepGen/Cards/LpairHandler.cpp
index 4eb1f49..903e751 100644
--- a/CepGen/Cards/LpairHandler.cpp
+++ b/CepGen/Cards/LpairHandler.cpp
@@ -1,148 +1,148 @@
 #include "LpairHandler.h"
 #include "CepGen/Core/Exception.h"
 
 #include "CepGen/Processes/GamGamLL.h"
 #include "CepGen/Processes/PPtoLL.h"
 #include "CepGen/Processes/PPtoWW.h"
 
 #include "CepGen/Hadronisers/Pythia8Hadroniser.h"
 
 #include <fstream>
 
 namespace CepGen
 {
   namespace Cards
   {
     //----- specialization for LPAIR input cards
 
     LpairHandler::LpairHandler( const char* file ) :
       pair_( Particle::invalidParticle )
     {
       std::ifstream f( file, std::fstream::in );
       if ( !f.is_open() ) {
         FatalError( Form( "Failed to parse file \"%s\"", file ) );
         return;
       }
       init( &params_ );
 
       std::ostringstream os;
       os << Form( "File '%s' succesfully opened! The following parameters are set:\n", file );
 
       std::map<std::string, std::string> m_params;
       std::string key, value;
       while ( f >> key >> value ) {
         if ( key[0] == '#' ) continue; // FIXME need to ensure there is no extra space before!
         setParameter( key, value );
         m_params.insert( std::pair<std::string,std::string>( key, value ) );
         if ( getDescription( key ) != "null" ) os << ">> " << key << " = " << std::setw( 15 ) << getParameter( key ) << " (" << getDescription( key ) << ")" << std::endl;
       }
       f.close();
 
       if      ( proc_name_ == "lpair" )  params_.setProcess( new Process::GamGamLL() );
       else if ( proc_name_ == "pptoll" ) params_.setProcess( new Process::PPtoLL() );
       else if ( proc_name_ == "pptoww" ) params_.setProcess( new Process::PPtoWW() );
       else FatalError( Form( "Unrecognised process name: %s", proc_name_.c_str() ) );
 
       if      ( integr_type_ == "Vegas" ) params_.integrator.type = Integrator::Vegas;
       else if ( integr_type_ == "MISER" ) params_.integrator.type = Integrator::MISER;
 
 #ifdef PYTHIA8
       if ( hadr_name_ == "pythia8" ) params_.setHadroniser( new Hadroniser::Pythia8Hadroniser );
 #endif
 
       if ( m_params.count( "IEND" ) ) setValue<bool>( "IEND", ( std::stoi( m_params["IEND"] ) > 1 ) );
 
       //--- for LPAIR: specify the lepton pair to be produced
       if ( pair_ != Particle::invalidParticle ) {
         params_.kinematics.central_system = { pair_, pair_ };
       }
       Information( os.str() );
     }
 
     void
     LpairHandler::init( Parameters* params )
     {
       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<bool>( "IEND", "Generation type", &params->generation.enabled );
 
       registerParameter<unsigned int>( "DEBG", "Debugging verbosity", (unsigned int*)&Logger::get().level );
       registerParameter<unsigned int>( "NCVG", "Number of function calls", &params->integrator.ncvg );
       registerParameter<unsigned int>( "NCSG", "Number of points to probe", &params->integrator.npoints );
       registerParameter<unsigned int>( "ITVG", "Number of integration iterations", &params->integrator.itvg );
+      registerParameter<unsigned int>( "SEED", "Random generator seed", (unsigned int*)&params->integrator.seed );
       registerParameter<unsigned int>( "MODE", "Subprocess' mode", (unsigned int*)&params->kinematics.mode );
       registerParameter<unsigned int>( "PMOD", "Outgoing primary particles' mode", (unsigned int*)&params->kinematics.structure_functions );
       registerParameter<unsigned int>( "EMOD", "Outgoing primary particles' mode", (unsigned int*)&params->kinematics.structure_functions );
       registerParameter<unsigned int>( "PAIR", "Outgoing particles' PDG id", (unsigned int*)&pair_ );
       registerParameter<unsigned int>( "NGEN", "Number of events to generate", &params->generation.maxgen );
-      registerParameter<unsigned int>( "SEED", "Random generator seed", (unsigned int*)&params->vegas.seed );
 
       registerParameter<double>( "INPP", "Momentum (1st primary particle)", &params->kinematics.inp.first );
       registerParameter<double>( "INPE", "Momentum (2nd primary particle)", &params->kinematics.inp.second );
       registerParameter<double>( "PTCT", "Minimal transverse momentum (single central outgoing particle)", &params->kinematics.cuts.central[Cuts::pt_single].min() );
       registerParameter<double>( "MSCT", "Minimal central system mass", &params->kinematics.cuts.central[Cuts::mass_sum].min() );
       registerParameter<double>( "ECUT", "Minimal energy (single central outgoing particle)", &params->kinematics.cuts.central[Cuts::energy_single].min() );
       //registerParameter<double>( "THMN", "Minimal polar production angle for the central particles", &params->kinematics.eta_min );
       //registerParameter<double>( "THMX", "Maximal polar production angle for the central particles", &params->kinematics.eta_max );
       registerParameter<double>( "ETMN", "Minimal pseudo-rapidity (central outgoing particles)", &params->kinematics.cuts.central[Cuts::eta_single].min() );
       registerParameter<double>( "ETMX", "Maximal pseudo-rapidity (central outgoing particles)", &params->kinematics.cuts.central[Cuts::eta_single].max() );
       registerParameter<double>( "YMIN", "Minimal rapidity (central outgoing particles)", &params->kinematics.cuts.central[Cuts::rapidity_single].min() );
       registerParameter<double>( "YMAX", "Maximal rapidity (central outgoing particles)", &params->kinematics.cuts.central[Cuts::rapidity_single].max() );
       registerParameter<double>( "Q2MN", "Minimal Q^2 (exchanged parton)", &params->kinematics.cuts.initial[Cuts::q2].min() );
       registerParameter<double>( "Q2MX", "Maximal Q^2 (exchanged parton)", &params->kinematics.cuts.initial[Cuts::q2].max() );
       registerParameter<double>( "MXMN", "Minimal invariant mass of proton remnants", &params->kinematics.cuts.remnants[Cuts::mass].min() );
       registerParameter<double>( "MXMX", "Maximal invariant mass of proton remnants", &params->kinematics.cuts.remnants[Cuts::mass].max() );
     }
 
     void
     LpairHandler::store( const char* file )
     {
       std::ofstream f( file, std::fstream::out | std::fstream::trunc );
       if ( !f.is_open() ) {
         InError( Form( "Failed to open file \"%s\" for writing", file ) );
         return;
       }
       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( std::string key, std::string value )
     {
       try { setValue<double>( key.c_str(), std::stod( value ) ); } catch ( std::invalid_argument& ) {}
       try { setValue<unsigned 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 );
 
       unsigned int ui = getValue<unsigned 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/Core/TamingFunction.h b/CepGen/Core/TamingFunction.h
index 27a2d81..8d5c01c 100644
--- a/CepGen/Core/TamingFunction.h
+++ b/CepGen/Core/TamingFunction.h
@@ -1,38 +1,39 @@
 #ifndef CepGen_Core_TamingFunction_h
 #define CepGen_Core_TamingFunction_h
 
 #include "Functional.h"
+#include <map>
 
 namespace CepGen
 {
   /// A collection of expression/variable with the associated functional
   struct TamingFunction
   {
     TamingFunction( const std::string& var, const std::string& expr ) : variable( var ), expression( expr ), function( expr, { { variable } } ) {}
     std::string variable, expression;
     Functional<1> function;
   };
   /// A collection of taming functions evaluator with helper classes
   class TamingFunctionsCollection : public std::map<std::string, TamingFunction>
   {
     public:
       /// Insert a new variable/expression into the collection
       void add( const std::string& var, const std::string& expr ) { emplace( var, TamingFunction( var, expr ) ); }
       /// Does the collection handle a taming function for a given variable?
       bool has( const std::string& var ) const { return ( find( var ) != end() ); }
       /// Evaluate the taming function for a given variable at a given value
       double eval( const std::string& var, double x ) const {
         const auto& it = find( var );
         if ( it == end() ) return 1.0;
         return it->second.function.eval( x );
       }
       /// Dump a full list of taming functions handled
       void dump( std::ostream& os=Logger::get().outputStream ) const {
         os << "List of taming functions:\n";
         for ( const auto& it : *this )
           os << ">> \"" << it.second.expression << "\" applied on variable \"" << it.first << "\"\n";
       }
   };
 }
 
 #endif
diff --git a/test/cepgen.cpp b/test/cepgen.cpp
index 3418eee..60fdeb5 100644
--- a/test/cepgen.cpp
+++ b/test/cepgen.cpp
@@ -1,83 +1,82 @@
 #include <iostream>
 
 #include "CepGen/Generator.h"
 
 #include "CepGen/Cards/LpairHandler.h"
 #include "CepGen/Cards/ConfigHandler.h"
 #include "CepGen/Core/Logger.h"
 #include "CepGen/Core/Exception.h"
 
 #include "CepGen/Processes/GamGamLL.h"
-#include "CepGen/Hadronisers/Pythia8Hadroniser.h"
 
 #include "CepGen/StructureFunctions/StructureFunctions.h"
 
 using namespace std;
 
 /**
  * Main caller for this Monte Carlo generator. Loads the configuration files'
  * variables if set as an argument to this program, else loads a default
  * "LHC-like" configuration, then launches the cross-section computation and
  * the events generation.
  * \author Laurent Forthomme <laurent.forthomme@cern.ch>
  */
 int main( int argc, char* argv[] ) {
   CepGen::Generator mg;
   
   //CepGen::Logger::get().level = CepGen::Logger::Debug;
   //CepGen::Logger::get().level = CepGen::Logger::DebugInsideLoop;
   //CepGen::Logger::get().outputStream( ofstream( "log.txt" ) );
   
   if ( argc == 1 ) {
     Information( "No config file provided. Setting the default parameters." );
     
     mg.parameters->setProcess( new CepGen::Process::GamGamLL );
     //mg.parameters->process_mode = Kinematics::InelasticElastic;
     mg.parameters->kinematics.mode = CepGen::Kinematics::ElasticElastic;
     mg.parameters->kinematics.structure_functions = CepGen::StructureFunctions::SuriYennie;
 #ifdef PYTHIA8
     mg.parameters->setHadroniser( new CepGen::Hadroniser::Pythia8Hadroniser );
 #endif
 
     mg.parameters->kinematics.inp = { 6500., 6500. };
     mg.parameters->kinematics.central_system = { CepGen::Particle::Muon, CepGen::Particle::Muon };
     mg.parameters->kinematics.cuts.central[CepGen::Cuts::pt_single].min() = 15.;
     mg.parameters->kinematics.cuts.central[CepGen::Cuts::eta_single] = { -2.5, 2.5 };
     mg.parameters->integrator.ncvg = 5e4; //FIXME
     mg.parameters->generation.enabled = true;
     //mg.parameters->maxgen = 2;
     mg.parameters->generation.maxgen = 2e4;
   }
   else {
     Information( Form( "Reading config file stored in %s", argv[1] ) );
     //CepGen::Cards::LpairReader card( argv[1] );
     const std::string file( argv[1] ), extension = file.substr( file.find_last_of( "." )+1 );
     if ( extension == "card" ) mg.setParameters( CepGen::Cards::LpairHandler( argv[1] ).parameters() );
     else if ( extension == "cfg" ) mg.setParameters( CepGen::Cards::ConfigHandler( argv[1] ).parameters() );
   }
 
   // We might want to cross-check visually the validity of our run
   mg.parameters->dump();
 
   // Let there be cross-section...
   double xsec, err;
   mg.computeXsection( xsec, err );
 
   if ( mg.parameters->generation.enabled ) {
     // The events generation starts here !
     CepGen::Event ev;
     for ( unsigned int i=0; i<mg.parameters->generation.maxgen; i++ ) {
       ev = *mg.generateOneEvent();
       if ( i%1000==0 ) {
         Information( Form( "Generating event #%d", i ) );
         ev.dump();
       }
     }
   }
 
   // store the current configuration
   CepGen::Cards::ConfigHandler::store( mg.parameters.get(), "last_run.cfg" );
 
   return 0;
 }
 
diff --git a/test/test_physics_processes.cpp b/test/test_physics_processes.cpp
index 8e477ed..e439624 100644
--- a/test/test_physics_processes.cpp
+++ b/test/test_physics_processes.cpp
@@ -1,147 +1,148 @@
 #include "CepGen/Generator.h"
 #include "CepGen/Parameters.h"
 #include "CepGen/Core/Timer.h"
 
 #include "CepGen/Processes/GamGamLL.h"
 #include "CepGen/Processes/PPtoLL.h"
 
 #include "CepGen/StructureFunctions/StructureFunctions.h"
 
 #include "abort.h"
 
 #include <unordered_map>
 #include <assert.h>
+#include <string.h>
 
 using namespace std;
 
 int
 main( int argc, char* argv[] )
 {
   typedef vector<pair<string,pair<double,double> > > KinematicsMap;
   typedef vector<pair<float, KinematicsMap> > ValuesAtCutMap;
 
   AbortHandler ctrl_c;
 
   // values defined at pt(single lepton)>15 GeV, |eta(single lepton)|<2.5, mX<1000 GeV
   // process -> { pt cut -> { kinematics -> ( sigma, delta(sigma) ) } }
   vector<pair<string,ValuesAtCutMap> > values_map = {
     //--- LPAIR values at sqrt(s) = 13 TeV
     { "lpair", {
       { 3.0, { // pt cut
         { "elastic",    { 2.0871703e1, 3.542e-2 } },
         { "singlediss", { 1.5042536e1, 3.256e-2 } },
         { "doublediss", { 1.38835e1, 4.03018e-2 } }
       } },
       { 15.0, { // pt cut
         { "elastic",    { 4.1994803e-1, 8.328e-4 } },
         { "singlediss", { 4.8504819e-1, 1.171e-3 } },
         { "doublediss", { 6.35650e-1, 1.93968e-3 } }
       } },
     } },
     //--- PPTOLL values
     { "pptoll", {
       { 3.0, { // pt cut
         { "elastic",       { 2.0936541e1, 1.4096e-2 } },
         { "singlediss_su", { 1.4844881e1, 2.0723e-2 } }, // SU, qt<50
         { "doublediss_su", { 1.3580637e1, 2.2497e-2 } }, // SU, qt<50
       } },
       { 15.0, { // pt cut
         { "elastic",       { 4.2515888e-1, 3.0351e-4 } },
         { "singlediss_su", { 4.4903253e-1, 5.8970e-4 } }, // SU, qt<50
         { "doublediss_su", { 5.1923819e-1, 9.6549e-4 } }, // SU, qt<50
         /*{ "2_singlediss", { 4.6710287e-1, 6.4842e-4 } }, // SU, qt<500
         { "3_doublediss", { 5.6316353e-1, 1.1829e-3 } }, // SU, qt<500*/
       } },
     } },
   };
 
   const double num_sigma = 3.0;
 
   if ( argc < 3 || strcmp( argv[2], "debug" ) != 0 ) {
     CepGen::Logger::get().level = CepGen::Logger::Nothing;
   }
 
   Timer tmr;
   CepGen::Generator mg;
 
   if ( argc > 1 && strcmp( argv[1], "vegas" ) == 0 ) mg.parameters->integrator.type = CepGen::Integrator::Vegas;
   if ( argc > 1 && strcmp( argv[1], "miser" ) == 0 ) mg.parameters->integrator.type = CepGen::Integrator::MISER;
   //{ ostringstream os; os << mg.parameters->integrator.type; Information( Form( "Testing with %s", os.str().c_str() ) ); }
   { cout << "Testing with " << mg.parameters->integrator.type << endl; }
 
   mg.parameters->kinematics.setSqrtS( 13.e3 );
   mg.parameters->kinematics.cuts.central[CepGen::Cuts::eta_single].in( -2.5, 2.5 );
   mg.parameters->kinematics.cuts.remnants[CepGen::Cuts::mass].max() = 1000.;
   //mg.parameters->integrator.ncvg = 50000;
   mg.parameters->integrator.itvg = 5;
 
   Information( Form( "Initial configuration time: %.3f ms", tmr.elapsed()*1.e3 ) );
   tmr.reset();
 
   unsigned short num_tests = 0, num_tests_passed = 0;
   vector<string> failed_tests, passed_tests;
 
   try {
     for ( const auto& values_vs_generator : values_map ) { // loop over all generators
       if      ( values_vs_generator.first == "lpair"  ) mg.parameters->setProcess( new CepGen::Process::GamGamLL );
       else if ( values_vs_generator.first == "pptoll" ) {
         mg.parameters->setProcess( new CepGen::Process::PPtoLL );
         mg.parameters->kinematics.cuts.initial[CepGen::Cuts::qt].max() = 50.0;
       }
       else { InError( Form( "Unrecognized generator mode: %s", values_vs_generator.first.c_str() ) ); break; }
 
       for ( const auto& values_vs_cut : values_vs_generator.second ) { // loop over the single lepton pT cut
         mg.parameters->kinematics.cuts.central[CepGen::Cuts::pt_single].min() = values_vs_cut.first;
         for ( const auto& values_vs_kin : values_vs_cut.second ) { // loop over all possible kinematics
           if      ( values_vs_kin.first.find( "elastic"    ) != string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::ElasticElastic;
           else if ( values_vs_kin.first.find( "singlediss" ) != string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::InelasticElastic;
           else if ( values_vs_kin.first.find( "doublediss" ) != string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::InelasticInelastic;
           else { InError( Form( "Unrecognized kinematics mode: %s", values_vs_kin.first.c_str() ) ); break; }
 
           if ( values_vs_kin.first.find( "_su" ) != string::npos ) mg.parameters->kinematics.structure_functions = CepGen::StructureFunctions::SzczurekUleshchenko;
           else mg.parameters->kinematics.structure_functions = CepGen::StructureFunctions::SuriYennie;
 
           Information( Form( "Process: %s/%s\n\tConfiguration time: %.3f ms", values_vs_generator.first.c_str(), values_vs_kin.first.c_str(), tmr.elapsed()*1.e3 ) );
           tmr.reset();
 
           mg.clearRun();
           const double xsec_ref = values_vs_kin.second.first, err_xsec_ref = values_vs_kin.second.second;
           double xsec_cepgen, err_xsec_cepgen;
           mg.computeXsection( xsec_cepgen, err_xsec_cepgen );
 
           const double sigma = ( fabs( xsec_ref-xsec_cepgen ) ) / sqrt( err_xsec_cepgen*err_xsec_cepgen + err_xsec_ref*err_xsec_ref );
 
           Information( Form( "Computed cross section:\n\tRef.   = %.3e +/- %.3e\n\tCepGen = %.3e +/- %.3e\n\tPull: %.6f", xsec_ref, err_xsec_ref, xsec_cepgen, err_xsec_cepgen, sigma ) );
 
           Information( Form( "Computation time: %.3f ms", tmr.elapsed()*1.e3 ) );
           tmr.reset();
 
           string test_res = values_vs_generator.first+":"+
                                  Form( "pt-gt-%.1f", values_vs_cut.first )+":"+
                                  values_vs_kin.first+":"
                                  "ref="+Form( "%.3e", xsec_ref )+":"
                                  "got="+Form( "%.3e", xsec_cepgen )+":"
                                  "pull="+Form( "%.3f", sigma );
           if ( fabs( sigma )<num_sigma ) {
             passed_tests.emplace_back( test_res );
             num_tests_passed++;
           }
           else {
             failed_tests.emplace_back( test_res );
           }
           num_tests++;
           cout << "Test " << num_tests_passed << "/" << num_tests << " passed!" << endl;
         }
       }
     }
   } catch ( CepGen::Exception& e ) {}
   if ( failed_tests.size() != 0 ) {
     ostringstream os_failed, os_passed;
     for ( const auto& fail : failed_tests ) os_failed << " (*) " << fail << endl;
     for ( const auto& pass : passed_tests ) os_passed << " (*) " << pass << endl;
     throw CepGen::Exception( __PRETTY_FUNCTION__, Form( "Some tests failed!\n%s\nPassed tests:\n%s", os_failed.str().c_str(), os_passed.str().c_str() ), CepGen::FatalError );
   }
   else Information( "ALL TESTS PASSED!" );
 
   return 0;
 }