diff --git a/Cards/pptoww_cfg.py b/Cards/pptoww_cfg.py
index ae683c4..61789a4 100644
--- a/Cards/pptoww_cfg.py
+++ b/Cards/pptoww_cfg.py
@@ -1,58 +1,58 @@
 import Cards.Core as cepgen
 from Cards.integrators_cff import miser as integrator
-#from Cards.pythia8_cff import pythia8 as hadroniser
+from Cards.pythia8_cff import pythia8 as hadroniser
 
-'''hadroniser.pythiaProcessConfiguration = (
+hadroniser.pythiaProcessConfiguration = (
     # process-specific
     '13:onMode = off', # disable muon decays
     '24:onMode = off', # disable all W decays, but...
     '24:onIfAny = 11 13' # enable e-nue + mu-numu final states
-)'''
+)
 
 process = cepgen.Module('pptoww',
     mode = cepgen.ProcessMode.InelasticElastic,
     inKinematics = cepgen.Parameters(
         cmEnergy = 13.e3,
         #structureFunctions = cepgen.StructureFunctions.SuriYennie,
         #structureFunctions = cepgen.StructureFunctions.FioreBrasse,
         #structureFunctions = cepgen.StructureFunctions.SzczurekUleshchenko,
         #structureFunctions = cepgen.StructureFunctions.ALLM91,
         #structureFunctions = cepgen.StructureFunctions.ALLM97,
         structureFunctions = cepgen.StructureFunctions.LUXlike,
     ),
     outKinematics = cepgen.Parameters(
         mx = (1.07, 1000.),
         qt = (0., 1000.),
         #--- extra cuts on the pt(W+) and pt(W-) plane
         ptdiff = (0., 2000.),
         #--- distance in rapidity between W+ and W-
         #rapiditydiff = (4., 5.),
         cuts = {
             # cuts on the single W level
             24: cepgen.Parameters(
                 pt = (0.,),
                 rapidity = (-6., 6.),
             ),
             # cuts on the W decay products
             11: cepgen.Parameters(
                 pt = (20.,),
                 eta = (-2.5, 2.5),
             ),
             13: cepgen.Parameters(
                 pt = (20.,),
                 eta = (-2.5, 2.5),
             )
         },
     )
 )
 
 #integrator.numPoints = 10000
 
 #--- import the default generation parameters
 from Cards.generator_cff import generator
-generator.numEvents = 100000
-generator.printEvery = 1
+generator.numEvents = 1000
+#generator.printEvery = 1
 
 #print(process)
 #print(integrator)
 #print(hadroniser)
diff --git a/Cards/pythia8_cff.py b/Cards/pythia8_cff.py
index bc299af..a4df1e3 100644
--- a/Cards/pythia8_cff.py
+++ b/Cards/pythia8_cff.py
@@ -1,46 +1,48 @@
 import Cards.Core as cepgen
 
 pythia8 = cepgen.Module('pythia8',
     seed = 0,
     maxTrials = 1,
     pythiaPreConfiguration = (
         #'Init:showAllSettings = on',
         # disable all generation processes
         'ProcessLevel:all = off',
         #'Check:event = off',
         # printout properties
         # start by disabling some unnecessary output
         'Next:numberCount = 0',
         #'Next:numberShowInfo = 0',
         #'Next:numberShowProcess = 0',
         # parameterise the fragmentation part
         #'BeamRemnants:beamJunction = on',
         #'ProcessLevel:resonanceDecays = on',
         #'PartonLevel:Remnants = on',
         #'HadronLevel:all = on',
         'Diffraction:doHard = on',
         'HardQCD:all = on',
         # disable all Bremsstrahlung/FSR photon production
         'PartonLevel:ISR = off',
         'PartonLevel:FSR = off',
         'ParticleDecays:allowPhotonRadiation = off',
     ),
     pythiaConfiguration = (
         'Tune:preferLHAPDF = 2',
         #'Main:timesAllowErrors = 10000',
         #'Check:epTolErr = 0.01',
         'Check:abortIfVeto = off',
         'Beams:setProductionScalesFromLHEF = off',
         'SLHA:keepSM = on',
         'SLHA:minMassSM = 1000.',
         'ParticleDecays:limitTau0 = on',
         'ParticleDecays:tau0Max = 10',
         # CUEP8M1 tuning
         'Tune:pp 14', # Monash 2013 tune by Peter Skands (January 2014)
         'MultipartonInteractions:pT0Ref = 2.4024',
         'MultipartonInteractions:ecmPow = 0.25208',
         'MultipartonInteractions:expPow = 1.6',
         #'Check:epTolErr = 10000.', #FIXME FIXME FIXME!!! BAD BAD BAD
         #'Check:epTolWarn = 10.', #FIXME FIXME FIXME!!! BAD BAD BAD
+        '24:onMode = off', # disable all W decays...
+        '13:onMode = off', # disable all muon decays...
     ),
 )
diff --git a/CepGen/Cards/Handler.h b/CepGen/Cards/Handler.h
index 7a1360c..ea1f331 100644
--- a/CepGen/Cards/Handler.h
+++ b/CepGen/Cards/Handler.h
@@ -1,30 +1,36 @@
 #ifndef CepGen_Cards_Handler_h
 #define CepGen_Cards_Handler_h
 
 #include "CepGen/Parameters.h"
 
 namespace CepGen
 {
   class Parameters;
   /// Location for all steering card parsers/writers
   namespace Cards
   {
     /// Generic steering card handler
     class Handler
     {
       public:
         /// Build a configuration from an external steering card
         Handler() {}
         ~Handler() {}
 
         /// Retrieve a configuration from a parsed steering cart
         Parameters& parameters() { return params_; }
+        /// Small utility to retrieve the extension of a filename
+        ///  (naive approach)
+        static std::string getExtension( const char* filename ) {
+          const std::string file( filename );
+          return file.substr( file.find_last_of( "." )+1 );
+        }
 
       protected:
         /// List of parameters parsed from a card handler
         Parameters params_;
     };
   }
 }
 
 #endif
diff --git a/CepGen/Cards/PythonHandler.cpp b/CepGen/Cards/PythonHandler.cpp
index 7b00b54..96e407d 100644
--- a/CepGen/Cards/PythonHandler.cpp
+++ b/CepGen/Cards/PythonHandler.cpp
@@ -1,467 +1,466 @@
 #include "PythonHandler.h"
 #include "CepGen/Core/Exception.h"
 
 #ifdef PYTHON
 
 #include "CepGen/Core/TamingFunction.h"
 
 #include "CepGen/Processes/GamGamLL.h"
 #include "CepGen/Processes/PPtoLL.h"
 #include "CepGen/Processes/PPtoWW.h"
 
 #include "CepGen/Hadronisers/Pythia8Hadroniser.h"
 
 #include <algorithm>
 
 #ifdef PYTHIA8
 void
 feedPythia( CepGen::Hadroniser::Pythia8Hadroniser* py8, PyObject* hadr, const char* config )
 {
   PyObject* ppc = CepGen::Cards::PythonHandler::getElement( hadr, config );
   if ( !ppc )
     return;
   if ( !PyTuple_Check( ppc ) ) {
     Py_DECREF( ppc );
     return;
   }
   for ( Py_ssize_t i = 0; i < PyTuple_Size( ppc ); ++i ) {
     PyObject* pln = PyTuple_GetItem( ppc, i );
     std::string config = CepGen::Cards::PythonHandler::decode( pln );
     Py_DECREF( pln );
     py8->readString( config );
   }
   Py_DECREF( ppc );
 }
 #endif
 
 namespace CepGen
 {
   namespace Cards
   {
     //----- specialization for CepGen input cards
     PythonHandler::PythonHandler( const char* file )
     {
       setenv( "PYTHONPATH", ".:..", 1 );
       std::string filename = getPythonPath( file );
       const size_t fn_len = filename.length()+1;
 #ifdef PYTHON2
       char* sfilename = new char[fn_len];
       sprintf( sfilename, "%s", filename.c_str() );
 #else
       wchar_t* sfilename = new wchar_t[fn_len];
       swprintf( sfilename, fn_len, L"%s", filename.c_str() );
 #endif
       Py_SetProgramName( sfilename );
 
       Py_InitializeEx( 0 );
       if ( !Py_IsInitialized() )
         throw Exception( __PRETTY_FUNCTION__, "Failed to initialise the python parser!", FatalError );
 
       Debugging( Form( "Initialised the Python cards parser\n\t"
                        "Python version: %s\n\t"
                        "Platform: %s", Py_GetVersion(), Py_GetPlatform() ) );
 
       PyObject* fn = encode( filename.c_str() );
       if ( !fn )
         throwPythonError( Form( "Failed to encode the configuration filename %s", filename.c_str() ) );
 
       PyObject* cfg = PyImport_Import( fn );
       Py_DECREF( fn );
       if ( !cfg )
         throwPythonError( Form( "Failed to parse the configuration card %s", file ) );
 
       PyObject* process = PyObject_GetAttrString( cfg, "process" );
       if ( !process )
         throwPythonError( Form( "Failed to extract a \"process\" keyword from the configuration card %s", file ) );
 
       //--- type of process to consider
       PyObject* pproc_name = PyDict_GetItem( process, encode( module_name_ ) );
       if ( !pproc_name )
         throwPythonError( Form( "Failed to extract the process name from the configuration card %s", file ) );
 
       const std::string proc_name = decode( pproc_name );
       Py_DECREF( pproc_name );
       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: %s", proc_name.c_str() ) );
 
       //--- process mode
       getParameter( process, "mode", (int&)params_.kinematics.mode );
 
       //--- process kinematics
       PyObject* pin_kinematics = getElement( process, "inKinematics" );
       if ( pin_kinematics ) {
         parseIncomingKinematics( pin_kinematics );
         Py_DECREF( pin_kinematics );
       }
 
       PyObject* pout_kinematics = getElement( process, "outKinematics" );
       if ( pout_kinematics ) {
         parseOutgoingKinematics( pout_kinematics );
         Py_DECREF( pout_kinematics );
       }
 
       Py_DECREF( process );
 
       //--- hadroniser parameters
       PyObject* phad = PyObject_GetAttrString( cfg, "hadroniser" );
       if ( phad ) {
         parseHadroniser( phad );
         Py_DECREF( phad );
       }
 
       //--- generation parameters
       PyObject* pint = PyObject_GetAttrString( cfg, "integrator" );
       if ( pint ) {
         parseIntegrator( pint );
         Py_DECREF( pint );
       }
-      std::cout << "tmp=" << params_.kinematics.mode << std::endl;
 
       PyObject* pgen = PyObject_GetAttrString( cfg, "generator" );
       if ( pgen ) {
         parseGenerator( pgen );
         Py_DECREF( pgen );
       }
 
       //--- taming functions
       PyObject* ptam = PyObject_GetAttrString( cfg, "tamingFunctions" );
       if ( ptam ) {
         parseTamingFunctions( ptam );
         Py_DECREF( ptam );
       }
 
       Py_DECREF( cfg );
 #ifdef PYTHON2
       Py_Finalize();
 #else
       if ( Py_IsInitialized() )
         throw Exception( __PRETTY_FUNCTION__, "Failed to unregister the python parser!", FatalError );
       if ( Py_FinalizeEx() != 0 )
         throw Exception( __PRETTY_FUNCTION__, "Failed to unregister the python parser!", FatalError );
 #endif
       if ( sfilename ) delete sfilename;
     }
 
     void
     PythonHandler::parseIncomingKinematics( PyObject* kin )
     {
       PyObject* ppz = getElement( kin, "pz" );
       if ( ppz ) {
         if ( PyTuple_Check( ppz ) && PyTuple_Size( ppz ) == 2 ) {
           double pz0 = PyFloat_AsDouble( PyTuple_GetItem( ppz, 0 ) );
           double pz1 = PyFloat_AsDouble( PyTuple_GetItem( ppz, 1 ) );
           params_.kinematics.inp = { pz0, pz1 };
         }
         Py_DECREF( ppz );
       }
       double sqrt_s = -1.;
       getParameter( kin, "cmEnergy", sqrt_s );
       if ( sqrt_s != -1. )
         params_.kinematics.setSqrtS( sqrt_s );
       getParameter( kin, "structureFunctions", (int&)params_.kinematics.structure_functions );
     }
 
     void
     PythonHandler::parseOutgoingKinematics( PyObject* kin )
     {
       PyObject* ppair = getElement( kin, "pair" );
       if ( ppair ) {
         if ( isInteger( ppair ) ) {
           ParticleCode pair = (ParticleCode)asInteger( ppair );
           params_.kinematics.central_system = { pair, pair };
         }
         Py_DECREF( ppair );
       }
       else if ( ppair && PyTuple_Check( ppair ) ) {
         if ( PyTuple_Size( ppair ) != 2 )
           FatalError( "Invalid value for in_kinematics.pair!" );
         ParticleCode pair1 = (ParticleCode)asInteger( PyTuple_GetItem( ppair, 0 ) );
         ParticleCode pair2 = (ParticleCode)asInteger( PyTuple_GetItem( ppair, 1 ) );
         params_.kinematics.central_system = { pair1, pair2 };
         Py_DECREF( ppair );
       }
 
       PyObject* pcuts = getElement( kin, "cuts" );
       if ( pcuts && PyDict_Check( pcuts ) ) parseParticlesCuts( pcuts );
 
       // for LPAIR/collinear matrix elements
       getLimits( kin, "q2", params_.kinematics.cuts.initial[Cuts::q2] );
 
       // for the kT factorised matrix elements
       getLimits( kin, "qt", params_.kinematics.cuts.initial[Cuts::qt] );
       getLimits( kin, "ptdiff", params_.kinematics.cuts.central[Cuts::pt_diff] );
       getLimits( kin, "rapiditydiff", params_.kinematics.cuts.central[Cuts::rapidity_diff] );
 
       getLimits( kin, "mx", params_.kinematics.cuts.remnants[Cuts::mass] );
     }
 
     void
     PythonHandler::parseParticlesCuts( PyObject* cuts )
     {
       PyObject* pkey = nullptr, *pvalue = nullptr;
       Py_ssize_t pos = 0;
       while ( PyDict_Next( cuts, &pos, &pkey, &pvalue ) ) {
         ParticleCode pdg = (ParticleCode)asInteger( pkey );
         getLimits( pvalue, "pt", params_.kinematics.cuts.central_particles[pdg][Cuts::pt_single] );
         getLimits( pvalue, "energy", params_.kinematics.cuts.central_particles[pdg][Cuts::energy_single] );
         getLimits( pvalue, "eta", params_.kinematics.cuts.central_particles[pdg][Cuts::eta_single] );
         getLimits( pvalue, "rapidity", params_.kinematics.cuts.central_particles[pdg][Cuts::rapidity_single] );
       }
     }
 
     void
     PythonHandler::parseIntegrator( PyObject* integr )
     {
       if ( !PyDict_Check( integr ) )
         throwPythonError( "Integrator object should be a dictionary!" );
       PyObject* palgo = getElement( integr, module_name_ );
       if ( !palgo )
         throwPythonError( "Failed to retrieve the integration algorithm name!" );
       std::string algo = decode( palgo );
       Py_DECREF( palgo );
       if ( algo == "Plain" )
         params_.integrator.type = Integrator::Plain;
       else if ( algo == "Vegas" ) {
         params_.integrator.type = Integrator::Vegas;
         getParameter( integr, "alpha", (double&)params_.integrator.vegas.alpha );
         getParameter( integr, "iterations", (unsigned long&)params_.integrator.vegas.iterations );
         getParameter( integr, "mode", (int&)params_.integrator.vegas.mode );
         getParameter( integr, "verbosity", (int&)params_.integrator.vegas.verbose );
       }
       else if ( algo == "MISER" ) {
         params_.integrator.type = Integrator::MISER;
         getParameter( integr, "estimateFraction", (double&)params_.integrator.miser.estimate_frac );
         getParameter( integr, "minCalls", (unsigned long&)params_.integrator.miser.min_calls );
         getParameter( integr, "minCallsPerBisection", (unsigned long&)params_.integrator.miser.min_calls_per_bisection );
         getParameter( integr, "alpha", (double&)params_.integrator.miser.alpha );
         getParameter( integr, "dither", (double&)params_.integrator.miser.dither );
       }
       else
         throwPythonError( Form( "Invalid integration algorithm: %s", algo.c_str() ) );
 
       getParameter( integr, "numPoints", (int&)params_.integrator.npoints );
       getParameter( integr, "numFunctionCalls", (int&)params_.integrator.ncvg );
       getParameter( integr, "seed", (unsigned long&)params_.integrator.seed );
     }
 
     void
     PythonHandler::parseGenerator( PyObject* gen )
     {
       if ( !PyDict_Check( gen ) )
         throwPythonError( "Generation information object should be a dictionary!" );
       params_.generation.enabled = true;
       getParameter( gen, "numEvents", (int&)params_.generation.maxgen );
       getParameter( gen, "printEvery", (int&)params_.generation.gen_print_every );
     }
 
     void
     PythonHandler::parseTamingFunctions( PyObject* tf )
     {
       if ( !PyList_Check( tf ) )
         throwPythonError( "Taming functions list should be a list!" );
 
       for ( Py_ssize_t i = 0; i < PyList_Size( tf ); ++i ) {
         PyObject* pit = PyList_GetItem( tf, i );
         if ( !PyDict_Check( pit ) )
           throwPythonError( Form( "Item %d is invalid", i ) );
         PyObject* pvar = getElement( pit, "variable" );
         PyObject* pexpr = getElement( pit, "expression" );
         params_.taming_functions->add( decode( pvar ), decode( pexpr ) );
         Py_DECREF( pvar );
         Py_DECREF( pexpr );
       }
     }
 
     void
     PythonHandler::parseHadroniser( PyObject* hadr )
     {
       if ( !PyDict_Check( hadr ) )
         throwPythonError( "Hadroniser object should be a dictionary!" );
 
       PyObject* pname = getElement( hadr, module_name_ );
       if ( !pname )
         throwPythonError( "Hadroniser name is required!" );
 
       std::string hadr_name = decode( pname );
       Py_DECREF( pname );
 
       if ( hadr_name == "pythia8" ) {
 #ifdef PYTHIA8
         Hadroniser::Pythia8Hadroniser* pythia8 = new Hadroniser::Pythia8Hadroniser( params_ );
         PyObject* pseed = getElement( hadr, "seed" );
         long long seed = -1ll;
         if ( pseed ) {
           if ( isInteger( pseed ) )
             seed = PyLong_AsLongLong( pseed );
           Py_DECREF( pseed );
         }
         pythia8->setSeed( seed );
         feedPythia( pythia8, hadr, "pythiaPreConfiguration" );
         pythia8->init();
         feedPythia( pythia8, hadr, "pythiaConfiguration" );
         feedPythia( pythia8, hadr, "pythiaProcessConfiguration" );
         params_.setHadroniser( pythia8 );
 #else
         InWarning( "Pythia8 is not linked to this instance... "
                    "Ignoring this part of the configuration file." )
 #endif
       }
     }
 
     //------------------------------------------------------------------
     // Python API helpers
     //------------------------------------------------------------------
 
     std::string
     PythonHandler::getPythonPath( const char* file )
     {
       std::string s_filename = file;
       s_filename = s_filename.substr( 0, s_filename.find_last_of( "." ) );
       std::replace( s_filename.begin(), s_filename.end(), '/', '.' );
       return s_filename;
     }
 
     void
     PythonHandler::throwPythonError( const std::string& message, const ExceptionType& type )
     {
       PyObject* ptype = nullptr, *pvalue = nullptr, *ptraceback = nullptr;
       PyErr_Fetch( &ptype, &pvalue, &ptraceback );
       PyErr_Clear();
       PyErr_NormalizeException( &ptype, &pvalue, &ptraceback );
       std::ostringstream oss; oss << message;
       if ( ptype == nullptr ) {
         Py_Finalize();
         throw Exception( __PRETTY_FUNCTION__, oss.str().c_str(), type );
       }
 
       oss << "\n\tError: "
 #ifdef PYTHON2
           << PyString_AsString( PyObject_Str( pvalue ) ); // deprecated in python v3+
 #else
           << _PyUnicode_AsString( PyObject_Str( pvalue ) );
 #endif
       Py_Finalize();
       throw Exception( __PRETTY_FUNCTION__, oss.str().c_str(), type );
     }
 
     const char*
     PythonHandler::decode( PyObject* obj )
     {
 #ifdef PYTHON2
       const char* str = PyString_AsString( obj ); // deprecated in python v3+
 #else
       const char* str = _PyUnicode_AsString( obj );
 #endif
       if ( !str )
         throwPythonError( "Failed to decode a Python object!" );
       return str;
     }
 
     PyObject*
     PythonHandler::encode( const char* str )
     {
       PyObject* obj = PyUnicode_FromString( str );
       if ( !obj )
         throwPythonError( Form( "Failed to encode the following string:\n\t%s", str ) );
       return obj;
     }
 
     PyObject*
     PythonHandler::getElement( PyObject* obj, const char* key )
     {
       PyObject* nink = encode( key );
       PyObject* pout = PyDict_GetItem( obj, nink );
       Py_DECREF( nink );
       return pout;
     }
 
     void
     PythonHandler::getLimits( PyObject* obj, const char* key, Kinematics::Limits& lim )
     {
       PyObject* pobj = getElement( obj, key );
       if ( !pobj )
         return;
       if ( !PyTuple_Check( pobj ) )
         FatalError( Form( "Invalid value retrieved for %s", key ) );
       if ( PyTuple_Size( pobj ) < 1 )
         FatalError( Form( "Invalid number of values unpacked for %s!", key ) );
       double min = PyFloat_AsDouble( PyTuple_GetItem( pobj, 0 ) );
       lim.min() = min;
       if ( PyTuple_Size( pobj ) > 1 ) {
         double max = PyFloat_AsDouble( PyTuple_GetItem( pobj, 1 ) );
         if ( max != -1 )
           lim.max() = max;
       }
       Py_DECREF( pobj );
     }
 
     void
     PythonHandler::getParameter( PyObject* parent, const char* key, int& out )
     {
       PyObject* pobj = getElement( parent, key );
       if ( !pobj )
         return;
 #ifdef PYTHON2
       if ( !PyInt_Check( pobj ) )
         throwPythonError( Form( "Object \"%s\" has invalid type", key ) );
       out = _PyInt_AsInt( pobj );
 #else
       if ( !PyLong_Check( pobj ) )
         throwPythonError( Form( "Object \"%s\" has invalid type", key ) );
       out = PyLong_AsLong( pobj );
 #endif
       Py_DECREF( pobj );
     }
 
     void
     PythonHandler::getParameter( PyObject* parent, const char* key, unsigned long& out )
     {
       PyObject* pobj = getElement( parent, key );
       if ( !pobj )
         return;
       if ( !PyLong_Check( pobj ) )
         throwPythonError( Form( "Object \"%s\" has invalid type", key ) );
       out = PyLong_AsUnsignedLong( pobj );
       Py_DECREF( pobj );
     }
 
     void
     PythonHandler::getParameter( PyObject* parent, const char* key, double& out )
     {
       PyObject* pobj = getElement( parent, key );
       if ( !pobj )
         return;
       if ( !PyFloat_Check( pobj ) )
         throwPythonError( Form( "Object \"%s\" has invalid type", key ) );
       out = PyFloat_AsDouble( pobj );
       Py_DECREF( pobj );
     }
 
     bool
     PythonHandler::isInteger( PyObject* obj )
     {
 #ifdef PYTHON2
       return PyInt_Check( obj );
 #else
       return PyLong_Check( obj );
 #endif
     }
 
     int
     PythonHandler::asInteger( PyObject* obj )
     {
 #ifdef PYTHON2
       return _PyInt_AsInt( obj );
 #else
       return PyLong_AsLong( obj );
 #endif
     }
   }
 }
 
 #endif
diff --git a/CepGen/Core/Generator.cpp b/CepGen/Core/Generator.cpp
index ea1c17a..20cd79a 100644
--- a/CepGen/Core/Generator.cpp
+++ b/CepGen/Core/Generator.cpp
@@ -1,158 +1,166 @@
 #include "CepGen/Generator.h"
 #include "CepGen/Parameters.h"
 #include "CepGen/Version.h"
 
 #include "CepGen/Core/Integrator.h"
 #include "CepGen/Core/Exception.h"
 
 #include "CepGen/Processes/GenericProcess.h"
 #include "CepGen/Event/Event.h"
 
 #include <fstream>
 
 namespace CepGen
 {
   Generator::Generator() :
     parameters( std::unique_ptr<Parameters>( new Parameters ) ),
     cross_section_( -1. ), cross_section_error_( -1. )
   {
     if ( Logger::get().level > Logger::Nothing ) {
       Debugging( "Generator initialized" );
       try { printHeader(); } catch ( Exception& e ) { e.dump(); }
     }
     srand( time( 0 ) ); // Random number initialization
   }
 
   Generator::Generator( Parameters* ip ) :
     parameters( ip ),
     cross_section_( -1. ), cross_section_error_( -1. )
   {}
 
   Generator::~Generator()
   {
     if ( parameters->generation.enabled && parameters->process() && parameters->process()->numGeneratedEvents()>0 ) {
       Information( Form( "Mean generation time / event: %g ms", parameters->process()->totalGenerationTime()*1.e3/parameters->process()->numGeneratedEvents() ) );
     }
   }
 
   size_t
   Generator::numDimensions() const
   {
     if ( !parameters->process() ) return 0;
     return parameters->process()->numDimensions( parameters->kinematics.mode );
   }
 
   void
   Generator::clearRun()
   {
     integrator_.reset();
     parameters->integrator.first_run = true;
     cross_section_ = cross_section_error_ = -1.;
   }
 
   void
   Generator::setParameters( Parameters& ip )
   {
     parameters = std::unique_ptr<Parameters>( new Parameters( ip ) ); // copy constructor
   }
 
   void
   Generator::printHeader()
   {
     std::string tmp;
     std::ostringstream os; os << "version " << version() << std::endl;
     std::ifstream hf( "README" );
     if ( !hf.good() )
       throw Exception( __PRETTY_FUNCTION__, "Failed to open README file", JustWarning );
     while ( true ) {
       if ( !hf.good() ) break;
       getline( hf, tmp );
       os << "\n " << tmp;
     }
     hf.close();
     Information( os.str().c_str() );
   }
 
   double
   Generator::computePoint( double* x )
   {
     prepareFunction();
     double res = f( x, numDimensions(), (void*)parameters.get() );
     std::ostringstream os;
     for ( unsigned int i = 0; i < numDimensions(); ++i )
       os << x[i] << " ";
     Debugging( Form( "Result for x[%zu] = ( %s):\n\t%10.6f", numDimensions(), os.str().c_str(), res ) );
     return res;
   }
 
   void
   Generator::computeXsection( double& xsec, double& err )
   {
     Information( "Starting the computation of the process cross-section" );
 
     // first destroy and recreate the integrator instance
     if ( !integrator_ )
       integrator_ = std::unique_ptr<Integrator>( new Integrator( numDimensions(), f, parameters.get() ) );
     else if ( integrator_->dimensions() != numDimensions() )
       integrator_.reset( new Integrator( numDimensions(), f, parameters.get() ) );
 
     if ( Logger::get().level >= Logger::Debug ) {
       std::ostringstream topo; topo << parameters->kinematics.mode;
       Debugging( Form( "New integrator instance created\n\t"
                        "Considered topology: %s case\n\t"
                        "Will proceed with %d-dimensional integration", topo.str().c_str(), numDimensions() ) );
     }
 
     try { prepareFunction(); } catch ( Exception& e ) { e.dump(); }
 
     const int res = integrator_->integrate( cross_section_, cross_section_error_ );
     if ( res != 0 ) throw Exception( __PRETTY_FUNCTION__, Form( "Error while computing the cross-section: return value = %d", res ), FatalError );
 
     xsec = cross_section_;
     err = cross_section_error_;
 
-    Information( Form( "Total cross section: %g +/- %g pb", xsec, err ) );
+    if ( xsec < 1.e-2 ) {
+      Information( Form( "Total cross section: %g +/- %g fb", xsec*1.e3, err*1.e3 ) );
+    }
+    else if ( xsec > 5.e2 ) {
+      Information( Form( "Total cross section: %g +/- %g nb", xsec*1.e-3, err*1.e-3 ) );
+    }
+    else {
+      Information( Form( "Total cross section: %g +/- %g pb", xsec, err ) );
+    }
   }
 
   std::shared_ptr<Event>
   Generator::generateOneEvent()
   {
     if ( cross_section_ < 0. )
       computeXsection( cross_section_, cross_section_error_ );
 
     bool good = false;
     while ( !good )
       good = integrator_->generateOneEvent();
 
     return parameters->generation.last_event;
   }
 
   void
   Generator::generate( std::function<void( const Event&, unsigned int& )> callback )
   {
     Information( Form( "%d events will be generated",
                        parameters->generation.maxgen ) );
 
     unsigned int i = 0;
     while ( i < parameters->generation.maxgen )
       if ( generateOneEvent() ) {
         callback( *parameters->generation.last_event, i );
         i++;
       }
 
     Information( Form( "%d events generated", i ) );
   }
 
   void
   Generator::prepareFunction()
   {
     if ( !parameters->process() ) {
       throw Exception( __PRETTY_FUNCTION__, "No process defined!", FatalError );
     }
     Kinematics kin = parameters->kinematics;
     parameters->process()->addEventContent();
     parameters->process()->setKinematics( kin );
     Debugging( "Function prepared to be integrated!" );
   }
 }
 
diff --git a/CepGen/Core/Integrand.cpp b/CepGen/Core/Integrand.cpp
index d167723..f477493 100644
--- a/CepGen/Core/Integrand.cpp
+++ b/CepGen/Core/Integrand.cpp
@@ -1,180 +1,180 @@
 #include "CepGen/Core/Timer.h"
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/TamingFunction.h"
 #include "CepGen/Event/Event.h"
 #include "CepGen/Event/Particle.h"
 #include "CepGen/Physics/Kinematics.h"
 
 #include "CepGen/Processes/GenericProcess.h"
 #include "CepGen/Hadronisers/GenericHadroniser.h"
 
 #include "CepGen/Parameters.h"
 
 #include <sstream>
 #include <fstream>
 
 namespace CepGen
 {
   double
   f( double* x, size_t ndim, void* params )
   {
     std::ostringstream os;
 
     Parameters* p = static_cast<Parameters*>( params );
     std::shared_ptr<Event> ev = p->process()->event();
 
     if ( p->process()->hasEvent() ) {
       p->process()->clearEvent();
 
       const Particle::Momentum p1( 0., 0.,  p->kinematics.inp.first ), p2( 0., 0., -p->kinematics.inp.second );
       p->process()->setIncomingKinematics( p1, p2 ); // at some point introduce non head-on colliding beams?
 
       DebuggingInsideLoop( Form( "Function f called -- some parameters:\n\t"
                                  "  pz(p1) = %5.2f  pz(p2) = %5.2f\n\t"
                                  "  remnant mode: %d",
                                  p->kinematics.inp.first, p->kinematics.inp.second, p->kinematics.structure_functions ) );
 
       if ( p->integrator.first_run ) {
 
         if ( Logger::get().level >= Logger::Debug ) {
           std::ostringstream oss; oss << p->kinematics.mode;
           Debugging( Form( "Process mode considered: %s", oss.str().c_str() ) );
         }
 
         Particle& op1 = ev->getOneByRole( Particle::OutgoingBeam1 ),
                  &op2 = ev->getOneByRole( Particle::OutgoingBeam2 );
 
         //--- add outgoing protons or remnants
         switch ( p->kinematics.mode ) {
           case Kinematics::ElectronProton: default: { InError( Form( "Process mode %d not yet handled!", (int)p->kinematics.mode ) ); }
           case Kinematics::ElasticElastic: break; // nothing to change in the event
           case Kinematics::ElasticInelastic:
-            op2.setPdgId( DiffrProt, 1 ); break;
+            op2.setPdgId( DiffractiveProton, 1 ); break;
           case Kinematics::InelasticElastic: // set one of the outgoing protons to be fragmented
-            op1.setPdgId( DiffrProt, 1 ); break;
+            op1.setPdgId( DiffractiveProton, 1 ); break;
           case Kinematics::InelasticInelastic: // set both the outgoing protons to be fragmented
-            op1.setPdgId( DiffrProt, 1 );
-            op2.setPdgId( DiffrProt, 1 );
+            op1.setPdgId( DiffractiveProton, 1 );
+            op2.setPdgId( DiffractiveProton, 1 );
             break;
         }
 
         //--- prepare the function to be integrated
         p->process()->prepareKinematics();
 
         //--- add central system
         Particles& central_system = ev->getByRole( Particle::CentralSystem );
         if ( central_system.size() == p->kinematics.central_system.size() ) {
           unsigned short i = 0;
           for ( Particles::iterator part = central_system.begin(); part != central_system.end(); ++part ) {
             if ( p->kinematics.central_system[i] == invalidParticle ) continue;
             part->setPdgId( p->kinematics.central_system[i] );
             part->computeMass();
             i++;
           }
         }
 
         p->process()->clearRun();
         p->integrator.first_run = false;
       }
     } // event is not empty
 
     p->process()->setPoint( ndim, x );
     if ( Logger::get().level >= Logger::DebugInsideLoop ) {
       std::ostringstream oss; for ( unsigned int i = 0; i < ndim; ++i ) { oss << x[i] << " "; }
       DebuggingInsideLoop( Form( "Computing dim-%d point ( %s)", ndim, oss.str().c_str() ) );
     }
 
     //--- from this step on, the phase space point is supposed to be set by 'GenericProcess::setPoint()'
 
     p->process()->beforeComputeWeight();
 
     Timer tmr; // start the timer
 
     double integrand = p->process()->computeWeight();
 
     if ( integrand <= 0. ) return 0.;
 
     //--- fill in the process' Event object
 
     p->process()->fillKinematics();
 
     //--- once the kinematics variables have been populated,
     //    can apply the collection of taming functions
 
     if ( p->taming_functions ) {
       if ( p->taming_functions->has( "m_central" )
         || p->taming_functions->has( "pt_central" ) ) {
         Particle::Momentum central_system;
         const Particles& cm_parts = ev->getByRole( Particle::CentralSystem );
         for ( Particles::const_iterator it_p = cm_parts.begin(); it_p != cm_parts.end(); ++it_p )
           central_system += it_p->momentum();
         if ( p->taming_functions->has( "m_central" ) )
           integrand *= p->taming_functions->eval( "m_central", central_system.mass() );
         if ( p->taming_functions->has( "pt_central" ) )
           integrand *= p->taming_functions->eval( "pt_central", central_system.pt() );
       }
       if ( p->taming_functions->has( "q2" ) ) {
         integrand *= p->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton1 ).momentum().mass() );
         integrand *= p->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton2 ).momentum().mass() );
       }
     }
 
     //--- full event content (+ hadronisation) if generating events
 
     if ( p->hadroniser() ) {
       double br = 0.; // branching fraction for all decays
-      if ( !p->hadroniser()->hadronise( *ev, br, true ) )
-      //if ( !p->hadroniser()->hadronise( *ev, br, p->storage() ) )
+      //if ( !p->hadroniser()->hadronise( *ev, br, true ) )
+      if ( !p->hadroniser()->hadronise( *ev, br, p->storage() ) )
         return 0.;
       integrand *= br;
     }
 
     //--- apply cuts on final state system (after hadronisation!)
     //    (watch out your cuts, as this might be extremely time-consuming...)
 
     if ( p->kinematics.cuts.central_particles.size() > 0 ) {
       std::map<ParticleCode,std::map<Cuts::Central,Kinematics::Limits> >::const_iterator it_c;
       const Particles cs = ev->getByRole( Particle::CentralSystem );
       for ( Particles::const_iterator it_p = cs.begin(); it_p != cs.end(); ++it_p ) {
         it_c = p->kinematics.cuts.central_particles.find( it_p->pdgId() );
         if ( it_c == p->kinematics.cuts.central_particles.end() ) continue;
         // retrieve all cuts associated to this final state particle
         const std::map<Cuts::Central,Kinematics::Limits>& cm = it_c->second;
         // apply these cuts on the given particle
         if ( cm.count( Cuts::pt_single ) > 0 && !cm.at( Cuts::pt_single ).passes( it_p->momentum().pt() ) ) return 0.;
         //std::cout << it_c->first << "\t" << it_p->momentum().pt() << "\t" << cm.at( Cuts::pt_single ).passes( it_p->momentum().pt() ) << std::endl;
         if ( cm.count( Cuts::energy_single ) > 0 && !cm.at( Cuts::energy_single ).passes( it_p->momentum().energy() ) ) return 0.;
         if ( cm.count( Cuts::eta_single ) > 0 && !cm.at( Cuts::eta_single ).passes( it_p->momentum().eta() ) ) return 0.;
         if ( cm.count( Cuts::rapidity_single ) > 0 && !cm.at( Cuts::rapidity_single ).passes( it_p->momentum().rapidity() ) ) return 0.;
       }
     }
 
     if ( p->storage() ) {
 
 /*    std::ofstream dump( "dump.txt", std::ios_base::app );
     for ( unsigned short i = 0; i < ndim; ++i )
       dump << x[i] << "\t";
     dump << "\n";*/
 
       ev->time_generation = tmr.elapsed();
 
       ev->time_total = tmr.elapsed();
       p->process()->addGenerationTime( ev->time_total );
 
       Debugging( Form( "Generation time:       %5.6f sec\n\t"
                        "Total time (gen+hadr): %5.6f sec",
                        ev->time_generation,
                        ev->time_total ) );
 
       p->generation.last_event = ev;
     } // generating events
 
     if ( Logger::get().level >= Logger::DebugInsideLoop ) {
       os.str( "" ); for ( unsigned int i = 0; i < ndim; ++i ) { os << Form( "%10.8f ", x[i] ); }
       Debugging( Form( "f value for dim-%d point ( %s): %4.4e", ndim, os.str().c_str(), integrand ) );
     }
 
     return integrand;
   }
 }
 
diff --git a/CepGen/Event/Event.cpp b/CepGen/Event/Event.cpp
index 06a2ffe..5d3d36b 100644
--- a/CepGen/Event/Event.cpp
+++ b/CepGen/Event/Event.cpp
@@ -1,326 +1,326 @@
 #include "Event.h"
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/utils.h"
 
 #include <algorithm>
 #include <math.h>
 
 namespace CepGen
 {
   Event::Event() :
     num_hadronisation_trials( 0 ),
     time_generation( -1. ), time_total( -1. )
   {}
 
   Event::Event( const Event& rhs ) :
-    particles_( rhs.particles_ ),
     num_hadronisation_trials( rhs.num_hadronisation_trials ),
     time_generation( rhs.time_generation ), time_total( rhs.time_total ),
+    particles_( rhs.particles_ ),
     evtcontent_( rhs.evtcontent_ )
   {}
 
   Event::~Event()
   {}
 
   Event&
   Event::operator=( const Event &ev_ )
   {
     particles_ = ev_.particles_;
     time_generation = ev_.time_generation;
     time_total = ev_.time_total;
     num_hadronisation_trials = ev_.num_hadronisation_trials;
     evtcontent_ = ev_.evtcontent_;
     return *this;
   }
 
   void
   Event::clear()
   {
     particles_.clear();
     time_generation = -1.;
     time_total = -1.;
   }
 
   void
   Event::freeze()
   {
     //--- store a snapshot of the primordial event block
     if ( particles_.count( Particle::CentralSystem ) > 0 )
       evtcontent_.cs = particles_[Particle::CentralSystem].size();
     if ( particles_.count( Particle::OutgoingBeam1 ) > 0 )
       evtcontent_.op1 = particles_[Particle::OutgoingBeam1].size();
     if ( particles_.count( Particle::OutgoingBeam2 ) > 0 )
       evtcontent_.op2 = particles_[Particle::OutgoingBeam2].size();
   }
 
   void
   Event::restore()
   {
     //--- remove all particles after the primordial event block
     if ( particles_.count( Particle::CentralSystem ) > 0 )
       particles_[Particle::CentralSystem].resize( evtcontent_.cs );
     if ( particles_.count( Particle::OutgoingBeam1 ) > 0 )
       particles_[Particle::OutgoingBeam1].resize( evtcontent_.op1 );
     if ( particles_.count( Particle::OutgoingBeam2 ) > 0 )
       particles_[Particle::OutgoingBeam2].resize( evtcontent_.op2 );
   }
 
   Particles&
   Event::getByRole( Particle::Role role )
   {
     //--- retrieve all particles with a given role
     return particles_[role];
   }
 
   const Particles&
   Event::getByRole( Particle::Role role ) const
   {
     //--- retrieve all particles with a given role
     return particles_.at( role );
   }
 
   ParticlesIds
   Event::getIdsByRole( Particle::Role role ) const
   {
     //--- retrieve all particles ids with a given role
     ParticlesIds out;
     Particles parts;
     try {
       parts = particles_.at( role );
     } catch ( std::out_of_range ) { return out; }
     for ( Particles::const_iterator it = parts.begin(); it != parts.end(); ++it ) {
       out.insert( it->id() );
     }
     return out;
   }
 
   Particle&
   Event::getOneByRole( Particle::Role role )
   {
     //--- retrieve the first particle a the given role
     Particles& parts_by_role = getByRole( role );
     if ( parts_by_role.size() == 0 )
       FatalError( Form( "No particle retrieved with role %d", (int)role ) );
     if ( parts_by_role.size() > 1 )
       FatalError( Form( "More than one particle with role %d: %d particles", (int)role, parts_by_role.size() ) );
     return *parts_by_role.begin();
   }
 
   const Particle&
   Event::getOneByRole( Particle::Role role ) const
   {
     //--- retrieve the first particle a the given role
     const Particles& parts_by_role = particles_.at( role );
     if ( parts_by_role.size() == 0 )
       FatalError( Form( "No particle retrieved with role %d", (int)role ) );
     if ( parts_by_role.size() > 1 )
       FatalError( Form( "More than one particle with role %d: %d particles", (int)role, parts_by_role.size() ) );
     return *parts_by_role.begin();
   }
 
   Particle&
   Event::getById( int id )
   {
     for ( ParticlesMap::iterator out = particles_.begin(); out != particles_.end(); ++out ) {
       for ( Particles::iterator part = out->second.begin(); part != out->second.end(); ++part ) {
         if ( part->id() == id ) return *part;
       }
     }
     throw Exception( __PRETTY_FUNCTION__, Form( "Failed to retrieve the particle with id=%d", id ), FatalError );
   }
 
   const Particle&
   Event::getConstById( int id ) const
   {
     for ( ParticlesMap::const_iterator out = particles_.begin(); out != particles_.end(); ++out ) {
       for ( Particles::const_iterator part = out->second.begin(); part != out->second.end(); ++part ) {
         if ( part->id() == id ) return *part;
       }
     }
     throw Exception( __PRETTY_FUNCTION__, Form( "Failed to retrieve the particle with id=%d", id ), FatalError );
   }
 
   Particles
   Event::getByIds( const ParticlesIds& ids ) const
   {
     Particles out;
     for ( ParticlesIds::const_iterator id = ids.begin(); id != ids.end(); ++id ) {
       out.emplace_back( getConstById( *id ) );
     }
     return out;
   }
 
   Particles
   Event::mothers( const Particle& part )
   {
     return getByIds( part.mothers() );
   }
 
   Particles
   Event::daughters( const Particle& part )
   {
     return getByIds( part.daughters() );
   }
 
   ParticleRoles
   Event::roles() const
   {
     ParticleRoles out;
     ParticlesMap::const_iterator it, end = particles_.end();
     for ( it = particles_.begin(); it != end; it = particles_.upper_bound( it->first ) ) {
       out.emplace_back( it->first );
     }
     return out;
   }
 
   Particle&
   Event::addParticle( Particle& part, bool replace )
   {
     DebuggingInsideLoop( Form( "Particle with PDGid = %d has role %d", part.pdgId(), part.role() ) );
     if ( part.role() <= 0 ) FatalError( Form( "Trying to add a particle with role=%d", (int)part.role() ) );
 
     //--- retrieve the list of particles with the same role
     Particles& part_with_same_role = getByRole( part.role() );
 
     //--- specify the id
     if ( part_with_same_role.empty() && part.id() < 0 ) part.setId( numParticles() ); // set the id if previously invalid/inexistent
     if ( !part_with_same_role.empty() ) {
       if ( replace ) part.setId( part_with_same_role[0].id() ); // set the previous id if replacing a particle
       else part.setId( numParticles() );
     }
 
     //--- add the particle to the collection
     if ( replace ) part_with_same_role = Particles( 1, part ); // generate a vector containing only this particle
     else part_with_same_role.emplace_back( part );
 
     return part_with_same_role.back();
   }
 
   Particle&
   Event::addParticle( Particle::Role role, bool replace )
   {
     Particle np( role );
     return addParticle( np, replace );
   }
 
   size_t
   Event::numParticles() const
   {
     size_t out = 0;
     for ( ParticlesMap::const_iterator it = particles_.begin(); it != particles_.end(); ++it ) {
       out += it->second.size();
     }
     return out;
   }
 
   const Particles
   Event::particles() const
   {
     Particles out;
     for ( ParticlesMap::const_iterator it = particles_.begin(); it != particles_.end(); ++it ) {
       out.insert( out.end(), it->second.begin(), it->second.end() );
     }
     std::sort( out.begin(), out.end() );
     return out;
   }
 
   const Particles
   Event::stableParticles() const
   {
     Particles out;
     for ( ParticlesMap::const_iterator it = particles_.begin(); it != particles_.end(); ++it ) {
       for ( Particles::const_iterator part = it->second.begin(); part != it->second.end(); ++part ) {
         if ( (short)part->status() > 0 ) out.emplace_back( *part );
       }
     }
     std::sort( out.begin(), out.end() );
     return out;
   }
 
   void
   Event::checkKinematics() const
   {
     try {
       const Particles& parts = particles();
       // check the kinematics through parentage
       for ( Particles::const_iterator p = parts.begin(); p != parts.end(); ++p ) {
         ParticlesIds daughters = p->daughters();
         if ( daughters.empty() ) continue;
         Particle::Momentum ptot;
         for ( ParticlesIds::const_iterator daugh = daughters.begin(); daugh != daughters.end(); ++daugh ) {
           const Particle& d = getConstById( *daugh );
           const ParticlesIds mothers = d.mothers();
           if ( mothers.size() > 1 ) {
             for ( ParticlesIds::const_iterator moth = mothers.begin(); moth != mothers.end(); ++moth ) {
               if ( *moth == p->id() ) continue;
               ptot -= getConstById( *moth ).momentum();
             }
           }
           ptot += d.momentum();
         }
         const double mass_diff = ( ptot-p->momentum() ).mass();
         if ( fabs( mass_diff ) > 1.e-10 ) {
           dump();
           throw Exception( __PRETTY_FUNCTION__, Form( "Error in momentum balance for particle %d: mdiff = %.5e", p->id(), mass_diff ), FatalError );
         }
       }
     } catch ( const Exception& e ) { throw Exception( __PRETTY_FUNCTION__, Form( "Event kinematics check failed:\n\t%s", e.what() ), FatalError ); }
   }
 
   void
   Event::dump( std::ostream& out, bool stable ) const
   {
     const Particles parts = ( stable ) ? stableParticles() : particles();
 
     std::ostringstream os;
 
     double pxtot = 0., pytot = 0., pztot = 0., etot = 0.;
     for ( Particles::const_iterator part_ref = parts.begin(); part_ref != parts.end(); ++part_ref ) {
       const Particle& part = *part_ref;
       const ParticlesIds mothers = part.mothers();
       {
         std::ostringstream oss;
         if ( part.pdgId() == invalidParticle && mothers.size() > 0 )
           for ( std::set<int>::const_iterator it_moth = mothers.begin(); it_moth != mothers.end(); ++it_moth )
             oss << ( ( it_moth != mothers.begin() ) ? "/" : "" ) << getConstById( *it_moth ).pdgId();
         else oss << part.pdgId();
-        os << Form( "\n %2d\t%+6d %-8s", part.id(), part.integerPdgId(), oss.str().c_str() );
+        os << Form( "\n %2d\t%-+7d %-10s", part.id(), part.integerPdgId(), oss.str().c_str() );
       }
       os << "\t";
       if ( part.charge() != 999. ) os << Form( "%6.2f ", part.charge() );
       else                         os << "\t";
       { std::ostringstream oss; oss << part.role(); os << Form( "%8s\t%6d\t", oss.str().c_str(), part.status() ); }
       if ( !mothers.empty() ) {
         std::ostringstream oss;
         for ( ParticlesIds::const_iterator moth = mothers.begin(); moth != mothers.end(); ++moth ) {
           oss << ( ( moth != mothers.begin() ) ? "+" : "" ) << *moth;
         }
         os << Form( "%6s ", oss.str().c_str() );
       }
       else os << "       ";
       const Particle::Momentum mom = part.momentum();
-      os << Form( "% 9.6e % 9.6e % 9.6e % 9.6e % 12.7f", mom.px(), mom.py(), mom.pz(), part.energy(), part.mass() );
+      os << Form( "% 9.6e % 9.6e % 9.6e % 9.6e % 12.5f", mom.px(), mom.py(), mom.pz(), part.energy(), part.mass() );
 
       // discard non-primary, decayed particles
       if ( (short)part.status() >= 0. ) {
         const int sign = ( part.status() == Particle::Undefined ) ? -1 : 1;
         pxtot += sign*mom.px();
         pytot += sign*mom.py();
         pztot += sign*mom.pz();
         etot += sign*part.energy();
       }
     }
     //--- set a threshold to the computation precision
     if ( fabs( pxtot ) < 1.e-10 ) pxtot = 0.;
     if ( fabs( pytot ) < 1.e-10 ) pytot = 0.;
     if ( fabs( pztot ) < 1.e-10 ) pztot = 0.;
     if ( fabs(  etot ) < 1.e-10 ) etot = 0.;
     //
     Information( Form( "Dump of event content:\n"
-    "Part.\tPDG id\t\tCharge\t   Role\tStatus\tMother\t\t\t\t4-Momentum (GeV)\t\tMass (GeV)\n"
-    "----\t------\t\t------\t   ----\t------\t------\t------------------------------------------------------  -----------"
+    " Id\tPDG id\tName\t\tCharge\t   Role\tStatus\tMother\tpx (GeV/c)    py (GeV/c)    pz (GeV/c)    E (GeV)\t M (GeV/c²)\n"
+    " --\t------\t----\t\t------\t   ----\t------\t------\t------------  ------------  ------------  ------------\t ----------"
     "%s\n"
-    "---------------------------------------------------------------------------------------------------------------------------\n"
-    "\t\t\t\t\t\tTotal: % 9.6e % 9.6e % 9.6e % 9.6e", os.str().c_str(), pxtot, pytot, pztot, etot ) );
+    " ----------------------------------------------------------------------------------------------------------------------------------\n"
+    "\t\t\t\t\t\t\tOut-in:% 9.6e % 9.6e % 9.6e % 9.6e", os.str().c_str(), pxtot, pytot, pztot, etot ) );
   }
 }
diff --git a/CepGen/Event/Particle.h b/CepGen/Event/Particle.h
index 80665a1..6a7de08 100644
--- a/CepGen/Event/Particle.h
+++ b/CepGen/Event/Particle.h
@@ -1,361 +1,362 @@
 #ifndef CepGen_Event_Particle_h
 #define CepGen_Event_Particle_h
 
 #include "CepGen/Physics/ParticleProperties.h"
 
 #include <set>
 #include <map>
 #include <vector>
 
 namespace CepGen
 {
 
   /// A set of integer-type particle identifiers
   typedef std::set<int> ParticlesIds;
 
   /// Kinematic information for one particle
   class Particle {
     public:
       /// Internal status code for a particle
       enum Status {
         PrimordialIncoming = -9,
-        DebugResonance = -4,
-        Resonance = -3,
+        DebugResonance = -5,
+        Resonance = -4,
+        Fragmented = -3,
         Propagator = -2,
         Incoming = -1,
         Undefined = 0,
         FinalState = 1,
         Undecayed = 2, Unfragmented = 3
       };
       /// Role of the particle in the process
       enum Role {
         UnknownRole = -1,
         IncomingBeam1 = 1, IncomingBeam2 = 2,
         OutgoingBeam1 = 3, OutgoingBeam2 = 5,
         CentralSystem = 6,
         Intermediate = 4,
         Parton1 = 41, Parton2 = 42, Parton3 = 43
       };
       /**
        * Container for a particle's 4-momentum, along with useful methods to ease the development of any matrix element level generator
        * \brief 4-momentum for a particle
        * \date Dec 2015
        * \author Laurent Forthomme <laurent.forthomme@cern.ch>
        */
       class Momentum {
         public:
           /// Build a 4-momentum at rest with an invalid energy (no mass information known)
           Momentum();
           /// Build a 4-momentum using its 3-momentum coordinates and its energy
           Momentum( double x, double y, double z, double t = -1. );
 
           // --- static definitions
 
           /// Build a 3-momentum from its three pseudo-cylindric coordinates
           static Momentum fromPtEtaPhi( double pt, double eta, double phi, double e = -1. );
           /// Build a 4-momentum from its scalar momentum, and its polar and azimuthal angles
           static Momentum fromPThetaPhi( double p, double theta, double phi, double e = -1. );
           /// Build a 4-momentum from its four momentum and energy coordinates
           static Momentum fromPxPyPzE( double px, double py, double pz, double e );
 
           // --- vector and scalar operators
 
           /// Scalar product of the 3-momentum with another 3-momentum
           double threeProduct( const Momentum& ) const;
           /// Scalar product of the 4-momentum with another 4-momentum
           double fourProduct( const Momentum& ) const;
           /// Add a 4-momentum through a 4-vector sum
           Momentum& operator+=( const Momentum& );
           /// Subtract a 4-momentum through a 4-vector sum
           Momentum& operator-=( const Momentum& );
           /// Scalar product of the 3-momentum with another 3-momentum
           double operator*=( const Momentum& );
           /// Multiply all 4-momentum coordinates by a scalar
           Momentum& operator*=( double c );
           /// Equality operator
           bool operator==( const Momentum& ) const;
           /// Human-readable format for a particle's momentum
           friend std::ostream& operator<<( std::ostream& os, const Particle::Momentum& mom );
 
           Momentum& betaGammaBoost( double gamma, double betagamma );
           /// Forward Lorentz boost
           Momentum& lorentzBoost( const Particle::Momentum& p );
 
           // --- setters and getters
 
           /// Set all the components of the 4-momentum (in GeV)
           void setP( double px, double py, double pz, double e );
           /// Set all the components of the 3-momentum (in GeV)
           void setP( double px, double py, double pz );
           /// Set an individual component of the 4-momentum (in GeV)
           void setP( unsigned int i, double p );
           /// Set the energy (in GeV)
           inline void setEnergy( double e ) { energy_ = e; }
           /// Compute the energy from the mass
           inline void setMass( double m ) { setMass2( m*m ); }
           /// Compute the energy from the mass
           void setMass2( double m2 );
           /// Get one component of the 4-momentum (in GeV)
           double operator[]( const unsigned int i ) const;
           /// Get one component of the 4-momentum (in GeV)
           double& operator[]( const unsigned int i );
           /// Momentum along the \f$x\f$-axis (in GeV)
           inline double px() const { return px_; }
           /// Momentum along the \f$y\f$-axis (in GeV)
           inline double py() const { return py_; }
           /// Longitudinal momentum (in GeV)
           inline double pz() const { return pz_; }
           /// Transverse momentum (in GeV)
           double pt() const;
           /// Squared transverse momentum (in GeV\f$^\textrm{2}\f$)
           inline double pt2() const { return ( px()*px()+py()*py() ); }
           /// 4-vector of double precision floats (in GeV)
           const std::vector<double> pVector() const;
           /// 3-momentum norm (in GeV)
           inline double p() const { return p_; }
           /// Squared 3-momentum norm (in GeV\f$^\textrm{2}\f$)
           inline double p2() const { return p_*p_; }
           /// Energy (in GeV)
           inline double energy() const { return energy_; }
           /// Squared energy (in GeV^2)
           inline double energy2() const { return energy_*energy_; }
           /// Squared mass (in GeV^2) as computed from its energy and momentum
           inline double mass2() const { return energy2()-p2(); }
           /// Mass (in GeV) as computed from its energy and momentum
           /// \note Returns \f$-\sqrt{|E^2-\mathbf{p}^2|}<0\f$ if \f$\mathbf{p}^2>E^2\f$
           double mass() const;
           /// Polar angle (angle with respect to the longitudinal direction)
           double theta() const;
           /// Azimutal angle (angle in the transverse plane)
           double phi() const;
           /// Pseudo-rapidity
           double eta() const;
           /// Rapidity
           double rapidity() const;
           /// Rotate the transverse components by an angle phi (and reflect the y coordinate)
           Momentum& rotatePhi( double phi, double sign );
           /// Rotate the particle's momentum by a polar/azimuthal angle
           Momentum& rotateThetaPhi( double theta_, double phi_ );
           /// Apply a \f$ z\rightarrow -z\f$ transformation
           inline Momentum& mirrorZ() { pz_ = -pz_; return *this; }
         private:
           /// Compute the 3-momentum's norm
           void computeP();
           /// Momentum along the \f$x\f$-axis
           double px_;
           /// Momentum along the \f$y\f$-axis
           double py_;
           /// Momentum along the \f$z\f$-axis
           double pz_;
           /// 3-momentum's norm (in GeV/c)
           double p_;
           /// Energy (in GeV)
           double energy_;
       };
       /// Human-readable format for a particle's PDG code
       friend std::ostream& operator<<( std::ostream& os, const ParticleCode& pc );
       /// Human-readable format for a particle's role in the event
       friend std::ostream& operator<<( std::ostream& os, const Particle::Role& rl );
       /// Compute the 4-vector sum of two 4-momenta
       friend Particle::Momentum operator+( const Particle::Momentum& mom1, const Particle::Momentum& mom2 );
       /// Compute the 4-vector difference of two 4-momenta
       friend Particle::Momentum operator-( const Particle::Momentum& mom1, const Particle::Momentum& mom2 );
       /// Scalar product of two 3-momenta
       friend double operator*( const Particle::Momentum& mom1, const Particle::Momentum& mom2 );
       /// Multiply all components of a 4-momentum by a scalar
       friend Particle::Momentum operator*( const Particle::Momentum& mom, double c );
       /// Multiply all components of a 4-momentum by a scalar
       friend Particle::Momentum operator*( double c, const Particle::Momentum& mom );
 
       //----- static getters
 
       /// Convert a polar angle to a pseudo-rapidity
       static double thetaToEta( double theta );
       /// Convert a pseudo-rapidity to a polar angle
       static double etaToTheta( double eta );
       /// Convert a pseudo-rapidity to a rapidity
       static double etaToY( double eta_, double m_, double pt_ );
 
       Particle();
       /// Build using the role of the particle in the process and its PDG id
       /// \param[in] pdgId ParticleCode (PDG ID)
       /// \param[in] role Role of the particle in the process
       /// \param[in] st Current status
       Particle( Role role, ParticleCode pdgId = invalidParticle, Status st = Undefined );
       /// Copy constructor
       Particle( const Particle& );
       inline ~Particle() {}
       /// Comparison operator (from unique identifier)
       bool operator<( const Particle& rhs ) const;
       /// Comparison operator (from their reference's unique identifier)
       //bool operator<( Particle *rhs ) const { return ( id < rhs->id ); }
       Particle& lorentzBoost( double m_, const Momentum& mom_ );
       /// Lorentz boost (shamelessly stolen from ROOT)
       std::vector<double> lorentzBoost( const Momentum& mom_ );
 
       // --- general particle properties
 
       /// Unique identifier (in a Event object context)
       int id() const { return id_; }
       //void setId( int id ) { id_ = id; }
       /// Set the particle unique identifier in an event
       void setId( int id ) { id_ = id; }
       /// Electric charge (given as a float number, for the quarks and bound states)
       float charge() const { return charge_sign_ * ParticleProperties::charge( pdg_id_ ); }
       /// Set the electric charge sign (+-1 for charged or 0 for neutral particles)
       void setChargeSign( int sign ) { charge_sign_ = sign; }
       /// Role in the considered process
       Role role() const { return role_; }
       /// Set the particle role in the process
       void setRole( const Role& role ) { role_ = role; }
       /**
        * Codes 1-10 correspond to currently existing partons/particles, and larger codes contain partons/particles which no longer exist, or other kinds of event information
        * \brief Particle status
        */
       Status status() const { return status_; }
       /// Set the particle decay/stability status
       void setStatus( Status status ) { status_ = status; }
 
       /// Set the PDG identifier (along with the particle's electric charge)
       /// \param[in] pdg ParticleCode (PDG ID)
       /// \param[in] ch Electric charge (0, 1, or -1)
       void setPdgId( const ParticleCode& pdg, short ch = 0 );
       /// Retrieve the objectified PDG identifier
       inline ParticleCode pdgId() const { return pdg_id_; }
       /// Retrieve the integer value of the PDG identifier
       int integerPdgId() const;
       /// Particle's helicity
       float helicity() const { return helicity_; }
       /// Set the helicity of the particle
       void setHelicity( float heli ) { helicity_ = heli; }
       /**
        * Gets the particle's mass in \f$\textrm{GeV}/c^{2}\f$.
        * \brief Gets the particle's mass
        * \return The particle's mass
        */
       inline double mass() const { return mass_; };
       /**
        * Set the mass of the particle in \f$\textrm{GeV}/c^{2}\f$ while ensuring that the kinematics is properly set (the mass is set according to the energy and the momentum in priority)
        * \brief Compute the particle's mass in \f$\textrm{GeV}/c^{2}\f$
        */
       void computeMass( bool off_shell=false );
       /**
        * Set the mass of the particle in \f$\textrm{GeV}/c^{2}\f$ according to a value given as an argument. This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority)
        * \param m The mass in \f$\textrm{GeV}/c^{2}\f$ to set
        * \brief Set the particle's mass in \f$\textrm{GeV}/c^{2}\f$
        */
       void setMass( double m=-1. );
       /// Get the particle's squared mass (in \f$\textrm{GeV}^\textrm{2}\f$)
       inline double mass2() const { return mass_*mass_; };
       /// Retrieve the momentum object associated with this particle
       inline Momentum& momentum() { return momentum_; }
       /// Retrieve the momentum object associated with this particle
       inline Momentum momentum() const { return momentum_; }
       /// Associate a momentum object to this particle
       void setMomentum( const Momentum& mom, bool offshell = false );
       /**
        * \brief Set the 3-momentum associated to the particle
        * \param[in] px Momentum along the \f$x\f$-axis, in \f$\textrm{GeV}/c\f$
        * \param[in] py Momentum along the \f$y\f$-axis, in \f$\textrm{GeV}/c\f$
        * \param[in] pz Momentum along the \f$z\f$-axis, in \f$\textrm{GeV}/c\f$
        */
       void setMomentum( double px, double py, double pz );
       /**
        * \brief Set the 4-momentum associated to the particle
        * \param[in] px Momentum along the \f$x\f$-axis, in \f$\textrm{GeV}/c\f$
        * \param[in] py Momentum along the \f$y\f$-axis, in \f$\textrm{GeV}/c\f$
        * \param[in] pz Momentum along the \f$z\f$-axis, in \f$\textrm{GeV}/c\f$
        * \param[in] e Energy, in GeV
        */
       void setMomentum( double px, double py, double pz, double e );
       /**
        * \brief Set the 4-momentum associated to the particle
        * \param[in] p 4-momentum
        */
       inline void setMomentum( double p[4] ) { setMomentum( p[0], p[1], p[2], p[3] ); }
       /**
        * \brief Set the particle's energy
        * \param[in] e Energy, in GeV
        */
       void setEnergy( double e=-1. );
       /// Get the particle's energy (in GeV)
       double energy() const;
       /// Get the particle's squared energy (in \f$\textrm{GeV}^\textrm{2}\f$)
       inline double energy2() const { return energy()*energy(); };
       /// Is this particle a valid particle which can be used for kinematic computations?
       bool valid();
 
       // --- particle relations
 
       /// Is this particle a primary particle?
       inline bool primary() const { return mothers_.empty(); }
       /**
        * \brief Set the mother particle
        * \param[in] part A Particle object containing all the information on the mother particle
        */
       void addMother( Particle& part );
       /**
        * \brief Gets the unique identifier to the mother particle from which this particle arises
        * \return An integer representing the unique identifier to the mother of this particle in the event
        */
       inline ParticlesIds mothers() const { return mothers_; }
       /**
        * \brief Add a decay product
        * \param[in] part The Particle object in which this particle will desintegrate or convert
        * \return A boolean stating if the particle has been added to the daughters list or if it was already present before
        */
       void addDaughter( Particle& part );
       /// Gets the number of daughter particles
       inline unsigned int numDaughters() const { return daughters_.size(); };
       /**
        * \brief Get an identifiers list all daughter particles
        * \return An integer vector containing all the daughters' unique identifier in the event
        */
       inline ParticlesIds daughters() const { return daughters_; }
 
       // --- global particle information extraction
 
       /// Dump all the information on this particle into the standard output stream
       void dump() const;    
 
     private:
       /// Unique identifier in an event
       int id_;
       /// Electric charge (+-1 or 0)
       short charge_sign_;
       /// Momentum properties handler
       Momentum momentum_;
       /// Mass in \f$\textrm{GeV}/c^2\f$
       double mass_;
       /// Helicity
       float helicity_;
       /// Role in the process
       Role role_;
       /// Decay/stability status
       Status status_;
       /// List of mother particles
       ParticlesIds mothers_;
       /// List of daughter particles
       ParticlesIds daughters_;
       /// PDG id
       ParticleCode pdg_id_;
   };
 
   /// Compute the centre of mass energy of two particles (incoming or outgoing states)
   double CMEnergy( const Particle& p1, const Particle& p2 );
   /// Compute the centre of mass energy of two particles (incoming or outgoing states)
   double CMEnergy( const Particle::Momentum& m1, const Particle::Momentum& m2 );
 
   //bool operator<( const Particle& a, const Particle& b ) { return a.id<b.id; }
 
   // --- particle containers
 
   /// List of Particle objects
   typedef std::vector<Particle> Particles;
   /// List of particles' roles
   typedef std::vector<Particle::Role> ParticleRoles;
   /// Map between a particle's role and its associated Particle object
   typedef std::map<Particle::Role,Particles> ParticlesMap;
 }
 
 #endif
diff --git a/CepGen/Hadronisers/Pythia8Hadroniser.cpp b/CepGen/Hadronisers/Pythia8Hadroniser.cpp
index 9ec9c13..9c422aa 100644
--- a/CepGen/Hadronisers/Pythia8Hadroniser.cpp
+++ b/CepGen/Hadronisers/Pythia8Hadroniser.cpp
@@ -1,267 +1,270 @@
 #ifdef PYTHIA8
 #include "Pythia8Hadroniser.h"
 
 #include "CepGen/Parameters.h"
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/utils.h"
 #include "CepGen/Event/Event.h"
 #include "CepGen/Event/Particle.h"
 #include "CepGen/Version.h"
 
 namespace CepGen
 {
   namespace Hadroniser
   {
     Pythia8Hadroniser::Pythia8Hadroniser( const Parameters& params ) :
       GenericHadroniser( "pythia8" ), max_attempts_( params.hadroniser_max_trials )
     {
 #ifdef PYTHIA8
       pythia_.reset( new Pythia8::Pythia );
-      std::cout << params.kinematics.inpdg.first << "\t" << params.kinematics.inpdg.second << std::endl;
       readString( Form( "Beams:idA = %d", params.kinematics.inpdg.first ) );
       readString( Form( "Beams:idB = %d", params.kinematics.inpdg.second ) );
       readString( Form( "Beams:eCM = %.2f", params.kinematics.sqrtS() ) );
 #endif
     }
 
     Pythia8Hadroniser::~Pythia8Hadroniser()
     {}
 
     bool
     Pythia8Hadroniser::init()
     {
       bool res = pythia_->init();
       if ( !res )
         FatalError( "Failed to initialise the Pythia8 core!\n\t"
                     "See the message above for more details." );
       return res;
     }
 
     void
     Pythia8Hadroniser::readString( const char* param )
     {
       if ( !pythia_->readString( param ) )
         FatalError( Form( "The Pythia8 core failed to parse the following setting:\n\t%s", param ) );
     }
 
     void
     Pythia8Hadroniser::setSeed( long long seed )
     {
 #ifdef PYTHIA8
       if ( seed == -1ll ) {
         pythia_->settings.flag( "Random:setSeed", false );
         return;
       }
       pythia_->settings.flag( "Random:setSeed", true );
       pythia_->settings.parm( "Random:seed", seed );
 #endif
     }
 
     bool
     Pythia8Hadroniser::hadronise( Event& ev, double& weight, bool proton_fragment )
     {
       weight = 1.;
 #ifdef PYTHIA8
       //--- start by cleaning up the previous runs leftovers
       const double mp = ParticleProperties::mass( Proton ), mp2 = mp*mp;
 
       pythia_->event.reset();
       /*const double sqrt_s = ev.cmEnergy();
       pythia_->event[0].e( sqrt_s );
       pythia_->event[0].m( sqrt_s );*/
+      // handle slight energy conservation violation from the process...
       Particle::Momentum p_out = ev.getOneByRole( Particle::OutgoingBeam1 ).momentum()
                                + ev.getOneByRole( Particle::OutgoingBeam2 ).momentum();
       for ( const auto& p : ev.getByRole( Particle::CentralSystem ) )
         p_out += p.momentum();
       //std::cout << p_out.pz() << "|" << p_out.mass() << "|" << p_out.energy() << std::endl;
       pythia_->event[0].e( p_out.energy() );
       pythia_->event[0].px( p_out.px() );
       pythia_->event[0].py( p_out.py() );
       pythia_->event[0].pz( p_out.pz() );
       pythia_->event[0].m( p_out.mass() );
 
       const unsigned short num_before = ev.numParticles();
       std::map<short,short> py_cg_corresp, cg_py_corresp;
 
       //--- loop to add the particles
 
       unsigned short idx_remn1 = invalid_idx_, idx_remn2 = invalid_idx_;
       for ( unsigned short i = 0; i < num_before; ++i ) {
         const Particle& part = ev.getConstById( i );
 
         const Particle::Momentum mom = part.momentum();
         Pythia8::Particle py8part( part.integerPdgId(), 0, 0, 0, 0, 0, 0, 0, mom.px(), mom.py(), mom.pz(), mom.energy(), part.mass() );
         unsigned short py_id = invalid_idx_;
         switch ( part.role() ) {
           case Particle::IncomingBeam1:
           case Particle::IncomingBeam2:
           case Particle::Parton1:
           case Particle::Parton2:
           case Particle::Parton3:
           case Particle::Intermediate:
           case Particle::UnknownRole:
             continue;
           case Particle::CentralSystem: {
-            if ( pythia_->particleData.canDecay( py8part.id() )
-              || pythia_->particleData.mayDecay( py8part.id() ) )
-              py8part.status( 93 );
-            else
-              py8part.status( 23 ); // outgoing particles of the hardest subprocess
+            py8part.status( 23 ); // outgoing particles of the hardest subprocess
             py_id = pythia_->event.append( py8part );
           } break;
           case Particle::OutgoingBeam1:
           case Particle::OutgoingBeam2: {
-            py8part.id( 2212 );
-            py8part.status( ( proton_fragment && part.status() == Particle::Unfragmented )
-              ? 15
-              : 14 // final state proton
+            py8part.status(
+              ( proton_fragment && part.status() == Particle::Unfragmented )
+                ? 15
+                : 14 // final state proton
             );
             py_id = pythia_->event.append( py8part );
 
             if ( proton_fragment && part.status() == Particle::Unfragmented ) {
-              if ( part.role() == Particle::OutgoingBeam1 ) idx_remn1 = py_id;
-              if ( part.role() == Particle::OutgoingBeam2 ) idx_remn2 = py_id;
+              if ( part.role() == Particle::OutgoingBeam1 )
+                idx_remn1 = py_id;
+              if ( part.role() == Particle::OutgoingBeam2 )
+                idx_remn2 = py_id;
             }
           } break;
         }
         // populate the CepGen id <-> Pythia id map
         if ( py_id != invalid_idx_ ) {
           cg_py_corresp[part.id()] = py_id;
           py_cg_corresp[py_id] = part.id();
         }
       }
       if ( proton_fragment ) {
         if ( idx_remn1 != invalid_idx_ ) {
           const Particle::Momentum& p0 = ev.getOneByRole( Particle::IncomingBeam1 ).momentum();
           const Particle::Momentum& p = ev.getOneByRole( Particle::OutgoingBeam1 ).momentum();
           const double q2 = -( p-p0 ).mass2();
           const double mx2 = p.mass2();
           const double xbj = q2/( q2+mx2-mp2 );
           fragmentState( idx_remn1, xbj );
         }
         if ( idx_remn2 != invalid_idx_ ) {
           const Particle::Momentum& p0 = ev.getOneByRole( Particle::IncomingBeam2 ).momentum();
           const Particle::Momentum& p = ev.getOneByRole( Particle::OutgoingBeam2 ).momentum();
           const double q2 = -( p-p0 ).mass2();
           const double my2 = p.mass2();
           const double xbj = q2/( q2+my2-mp2 );
           fragmentState( idx_remn2, xbj );
         }
       }
 
       const unsigned short num_py_parts = pythia_->event.size();
 //pythia_->event.list(true,true);
       //ev.dump();
 //      exit(0);
 
       //std::cout << ":::::" << pythia_->settings.parm("Check:mTolErr") << std::endl;
 
       bool success = false;
       ev.num_hadronisation_trials = 0;
       while ( !success ) {
-        success = pythia_->next();
-        if ( proton_fragment ){
-          std::cout << success << "attempt " << ev.num_hadronisation_trials << " / " << max_attempts_ << std::endl;
+        //success = pythia_->next();
+        pythia_->next(); success = pythia_->event.size() != num_py_parts; //FIXME discards any pythia error!
+        /*if ( proton_fragment ) {
+          //std::cout << success << "attempt " << ev.num_hadronisation_trials << " / " << max_attempts_ << std::endl;
           pythia_->event.list(true,true);
-      }
+        }*/
         if ( ++ev.num_hadronisation_trials > max_attempts_ )
           return false;
       }
       /*for ( unsigned short i = 0; i < pythia_->event.size(); ++i ) {
         const double err = abs( pythia_->event[i].mCalc()-pythia_->event[i].m() )/std::max( 1.0, pythia_->event[i].e());
         std::cout << ">> " << pythia_->event[i].id() << "::" << err << "::" << (err>pythia_->settings.parm("Check:mTolErr")) << std::endl;
       }*/
       //std::cout << "bad" << std::endl;
       //InWarning( "Pythia8 failed to process the event." );
       //pythia_->event.list( true, true );
 
       // check if something happened in the event processing by Pythia
       // if not, return the event as it is...
       if ( pythia_->event.size() == num_py_parts )
         return true;
 
       for ( unsigned short i = 1; i < pythia_->event.size(); ++i ) { // skip the central system
         const Pythia8::Particle& p = pythia_->event[i];
         std::map<short,short>::const_iterator it = py_cg_corresp.find( i );
         if ( it != py_cg_corresp.end() ) { // the particle is already in the event content
-          if ( p.daughterList().size() > 0 && i != idx_remn1 && i != idx_remn2 ) {
-            weight *= p.particleDataEntry().pickChannel().bRatio();
-            ev.getById( it->second ).setStatus( Particle::Resonance );
+          if ( p.daughterList().size() > 0 ) {
+            if ( i != idx_remn1 && i != idx_remn2 ) {
+              weight *= p.particleDataEntry().pickChannel().bRatio();
+              ev.getById( it->second ).setStatus( Particle::Resonance );
+            }
+            else
+              ev.getById( it->second ).setStatus( Particle::Fragmented );
           }
         }
         else { // the particle was not yet included in the CepGen event
           const std::vector<int> mothers = p.motherList();
           if ( mothers.size() == 0 ) continue; // isolated particle
 
           Particle::Role role = Particle::CentralSystem;
           std::map<short,short>::const_iterator it_m = py_cg_corresp.find( mothers[0] );
           if ( it_m != py_cg_corresp.end() ) {
             const Particle& moth = ev.getById( it_m->second );
             if ( mothers[0] == idx_remn1 || moth.role() == Particle::OutgoingBeam1 )
               role = Particle::OutgoingBeam1;
             else if ( mothers[0] == idx_remn2 || moth.role() == Particle::OutgoingBeam2 )
               role = Particle::OutgoingBeam2;
           }
 
           Particle& op = ev.addParticle( role );
           cg_py_corresp[op.id()] = i;
           py_cg_corresp[i] = op.id();
 
           op.setPdgId( static_cast<ParticleCode>( abs( p.id() ) ), p.charge() );
           op.setStatus( p.isFinal()
             ? Particle::FinalState
             : Particle::Propagator
           );
           op.setMomentum( Particle::Momentum( p.px(), p.py(), p.pz(), p.e() ) );
           for ( std::vector<int>::const_iterator moth = mothers.begin(); moth != mothers.end(); ++moth ) {
             if ( *moth != 0 && py_cg_corresp.count( *moth ) == 0 )
               FatalError( Form( "Particle with id=%d was not found in the event content!", *moth ) );
             op.addMother( ev.getById( py_cg_corresp[*moth] ) );
           }
         }
       }
 #else
       FatalError( "Pythia8 is not linked to this instance!" );
 #endif
       return true;
     }
 
     void
     Pythia8Hadroniser::fragmentState( unsigned short idx, double xbj )
     {
       const Pythia8::Particle& remn = pythia_->event[idx];
       // specify the quark/diquark flavours
       //FIXME naive approach (weighted by e_q^2/1-e_dq^2); to be improved
       const double rnd = 1./RAND_MAX * rand();
       unsigned short pdg_q = 0, pdg_dq = 0;
       if      ( rnd < 1./9. ) { pdg_q = 1; pdg_dq = 2203; }
       else if ( rnd < 5./9. ) { pdg_q = 2; pdg_dq = 2101; }
       else                    { pdg_q = 2; pdg_dq = 2103; }
       // then assign the quark/diquark a 4-momentum
       const double px_x = remn.px(), px_y = remn.py(), px_z = remn.pz(), ex = remn.e();
       const double xdq = 1.-xbj;
       //const double m_q = pythia_->particleData.m0( pdg_q );
       //const double m_dq = pythia_->particleData.m0( pdg_dq );
       // fractional momenta of the two partons:
       // -> x * p_X for the quark
       // -> ( 1-x ) * p_X for the diquark
       Pythia8::Particle diquark( pdg_dq, 63, idx, 0, 0, 0, 0, 100+idx, px_x*xdq, px_y*xdq, px_z*xdq, ex*xdq );
       Pythia8::Particle   quark( pdg_q,  63, idx, 0, 0, 0, 100+idx, 0, px_x*xbj, px_y*xbj, px_z*xbj, ex*xbj );
       //diquark.e( diquark.eCalc() );
       //quark.e( quark.eCalc() );
       diquark.m( diquark.mCalc() );
       quark.m( quark.mCalc() );
       //std::cout << "> " << xdq << "|" << xbj << std::endl;
       const unsigned short id_dq = pythia_->event.append( diquark );
       const unsigned short id_q = pythia_->event.append( quark );
       // keep up with the particles parentage
       pythia_->event[idx].daughter1( id_dq );
       pythia_->event[idx].daughter2( id_q );
       //std::cout << "-->" << diquark.m() << "|" << quark.m() << std::endl;
       // set the quark/diquark to be hadronised through a string
       pythia_->event[idx].status( -15 );
     }
   }
 }
 
 #endif
diff --git a/CepGen/Physics/ParticleProperties.h b/CepGen/Physics/ParticleProperties.h
index 45ebc07..f59140f 100644
--- a/CepGen/Physics/ParticleProperties.h
+++ b/CepGen/Physics/ParticleProperties.h
@@ -1,144 +1,144 @@
 #ifndef CepGen_Physics_ParticleProperties_h
 #define CepGen_Physics_ParticleProperties_h
 
 #include <iostream>
 
 namespace CepGen
 {
   /** Unique identifier for a particle type. From \cite Beringer:1900zz :
    * `The Monte Carlo particle numbering scheme [...] is intended to facilitate interfacing between event generators, detector simulators, and analysis packages used in particle physics.`
    * \brief PDG ids of all known particles
    */
   enum ParticleCode {
     invalidParticle = 0,
     //--- fundamental particles
     TopQuark = 6,
     Electron = 11, ElectronNeutrino = 12,
     Muon = 13, MuonNeutrino = 14,
     Tau = 15, TauNeutrino = 16,
     Gluon = 21, Photon = 22, Z = 23, W = 24,
     //--- composite particles
     PiPlus = 211, PiZero = 111,
+    KPlus = 321, DPlus = 411,
     Rho770_0 = 113, Rho1450_0 = 100113, Rho1700_0 = 30113,
     Eta = 221, Omega782 = 223,
     h1380_1 = 10333,
     JPsi= 443,
     Phi1680 = 100333,
     Upsilon1S = 553, Upsilon2S = 100553, Upsilon3S = 200553,
     Proton = 2212, Neutron = 2112,
     Pomeron = 990, Reggeon = 110,
-    DiffrProt = 9902210
+    DiffractiveProton = 9902210
   };
   namespace ParticleProperties
   {
     /**
      * \brief Gets the mass of a particle
      * \param pdgId ParticleCode (PDG ID)
      * \return Mass of the particle in \f$\textrm{GeV}/c^2\f$
      */
     inline double mass( const ParticleCode& pdgId ) {
       switch ( pdgId ) {
         case Electron:     return 0.510998928e-3;
         case Muon:         return 0.1056583715;
         case Tau:          return 1.77682;
         case TopQuark:     return 172.44;
         case ElectronNeutrino: case MuonNeutrino: case TauNeutrino: return 0.;
         case Gluon: case Photon: return 0.;
         case Z:            return 91.1876;
         case W:            return 80.385;
         case PiPlus:       return 0.13957018;
         case PiZero:       return 0.1349766;
+        case KPlus:        return 0.49368;
+        case DPlus:        return 1.86962;
         case JPsi:         return 20.;            //FIXME FIXME FIXME
         case Proton:       return 0.938272046;
         case Neutron:      return 0.939565346;
         case Upsilon1S:    return 9.46030;
         case Upsilon2S:    return 10.02326;
         case Upsilon3S:    return 10.3552;
         case Rho770_0:     return 0.77526;
         case Rho1450_0:    return 1.465;
         case Rho1700_0:    return 1.720;
         case h1380_1:      return 1.38619;
         case Eta:          return 0.547862;
         case invalidParticle:
         default:           return -1.;
       }
     }
-    /**
-     * \brief Gets the electric charge of a particle
-     * \param pdgId ParticleCode (PDG ID)
-     * \return Charge of the particle in \f$e\f$
-     */
-    inline double charge( const ParticleCode& pdgId ) {
-      switch ( pdgId ) {
-        case Proton: case DiffrProt: return +1.;
+    /// Electric charge of a particle, in \f$e\f$
+    /// \param[in] pdg_id PDG id
+    inline double charge( const ParticleCode& pdg_id ) {
+      switch ( pdg_id ) {
+        case Proton: case DiffractiveProton: return +1.;
         case Electron: case Muon: case Tau: return -1.;
-        case ElectronNeutrino: case MuonNeutrino: case TauNeutrino: return 0.;
-        case Gluon: case Z: case Photon: return 0.;
         case TopQuark: return +2./3;
         case W: return +1.;
-        case PiPlus: return +1.;
-        case PiZero: return 0.;
-        case Neutron: return 0.;
-        case Eta: return 0.;
+        case PiPlus: case KPlus: case DPlus: return +1.;
         default: return 0.;
       }
     }
-    /**
-     * \brief Total decay width of one unstable particle
-     * \param[in] pdgId ParticleCode (PDG ID)
-     * \return Decay width in GeV
-     */
-    inline double width( const ParticleCode& pdgId ) {
-      switch ( pdgId ) {
+    /// Electric charge of a particle, in \f$e\f$
+    /// \param[in] id integer PDG id
+    inline double charge( int id ) {
+      const short sign = id / abs( id );
+      return sign * charge( (ParticleCode)abs( id ) );
+    }
+    /// Total decay width of an unstable particle, in GeV
+    /// \param[in] pdg_id ParticleCode (PDG ID)
+    inline double width( const ParticleCode& pdg_id ) {
+      switch ( pdg_id ) {
         case JPsi:      return 5.; //FIXME
         case Z:         return 2.4952;
         case W:         return 2.085;
         case Upsilon1S: return 54.02e-6;
         case Upsilon2S: return 31.98e-6;
         case Upsilon3S: return 20.32e-6;
         case Rho770_0:  return 0.150; // PDG
         case Rho1450_0: return 0.400; // PDG
         case Rho1700_0: return 0.250; // PDG
         default:        return -1.;
       }
     }
   }
 
   inline std::ostream& operator<<( std::ostream& os, const ParticleCode& pc ) {
     switch ( pc ) {
-      case Electron:         return os << "e±";
-      case ElectronNeutrino: return os << "ν_e";
-      case Muon:             return os << "µ±";
-      case MuonNeutrino:     return os << "ν_µ";
-      case Tau:              return os << "τ±";
-      case TauNeutrino:      return os << "ν_τ";
+      case Electron:         return os << "e± ";
+      case ElectronNeutrino: return os << "ν_e ";
+      case Muon:             return os << "µ±  ";
+      case MuonNeutrino:     return os << "ν_µ  ";
+      case Tau:              return os << "τ±  ";
+      case TauNeutrino:      return os << "ν_τ  ";
       case Gluon:            return os << "gluon";
-      case Photon:           return os << "ɣ";
+      case Photon:           return os << "ɣ ";
       case Z:                return os << "Z";
-      case W:                return os << "W±";
-      case PiPlus:           return os << "π±";
-      case PiZero:           return os << "π⁰";
-      case Rho770_0:         return os << "ρ(770)₀";
-      case Rho1450_0:        return os << "ρ(1450)₀";
-      case Rho1700_0:        return os << "ρ(1700)₀";
-      case h1380_1:          return os << "h(1380)₁";
+      case W:                return os << "W± ";
+      case PiPlus:           return os << "π±  ";
+      case PiZero:           return os << "π⁰   ";
+      case KPlus:            return os << "K± ";
+      case DPlus:            return os << "D± ";
+      case Rho770_0:         return os << "ρ(770)₀  ";
+      case Rho1450_0:        return os << "ρ(1450)₀  ";
+      case Rho1700_0:        return os << "ρ(1700)₀  ";
+      case h1380_1:          return os << "h(1380)₁ ";
       case Eta:              return os << "η meson";
-      case Omega782:         return os << "ω(782)";
-      case JPsi:             return os << "J/ψ";
-      case Phi1680:          return os << "ɸ(1680)";
-      case Upsilon1S:        return os << "Υ(1S)";
-      case Upsilon2S:        return os << "Υ(2S)";
-      case Upsilon3S:        return os << "Υ(3S)";;
+      case Omega782:         return os << "ω(782) ";
+      case JPsi:             return os << "J/ψ ";
+      case Phi1680:          return os << "ɸ(1680) ";
+      case Upsilon1S:        return os << "Υ(1S) ";
+      case Upsilon2S:        return os << "Υ(2S) ";
+      case Upsilon3S:        return os << "Υ(3S) ";;
       case Proton:           return os << "proton";
-      case DiffrProt:        return os << "diff.prot.";
+      case DiffractiveProton:return os << "diffr.proton";
       case Neutron:          return os << "neutron";
       case Pomeron:          return os << "pomeron";
       case Reggeon:          return os << "reggeon";
       case TopQuark:         return os << "t";
       case invalidParticle:  return os << "[...]";
     }
     return os;
   }
 }
 
 #endif
diff --git a/CepGen/StructureFunctions/SuriYennie.cpp b/CepGen/StructureFunctions/SuriYennie.cpp
index 214b8a6..d0ce107 100644
--- a/CepGen/StructureFunctions/SuriYennie.cpp
+++ b/CepGen/StructureFunctions/SuriYennie.cpp
@@ -1,57 +1,62 @@
 #include "SuriYennie.h"
 #include "CepGen/Physics/ParticleProperties.h"
 #include <math.h>
 
 namespace CepGen
 {
   namespace SF
   {
     SuriYennie::Parameterisation
     SuriYennie::Parameterisation::standard()
     {
       Parameterisation p;
       p.C1 = 0.86926;
       p.C2 = 2.23422;
       p.D1 = 0.12549;
       p.rho2 = 0.585;
       p.Cp = 0.96;
       p.Bp = 0.63;
       return p;
     }
 
     SuriYennie::Parameterisation
     SuriYennie::Parameterisation::alternative()
     {
       Parameterisation p;
       p.C1 = 0.6303;
       p.C2 = 2.3049;
       p.D1 = 0.04681;
       p.rho2 = 1.05;
       p.Cp = 1.23;
       p.Bp = 0.61;
       return p;
     }
 
     SuriYennie::SuriYennie( const SuriYennie::Parameterisation& param ) :
-      FE( 0. ), FM( 0. ), params_( param )
+      F1( 0. ), FE( 0. ), FM( 0. ), params_( param )
     {}
 
     SuriYennie
     SuriYennie::operator()( double q2, double xbj ) const
     {
       const double mp = ParticleProperties::mass( Proton ), mp2 = mp*mp;
-      const double mx2 = q2 * ( 1.-xbj )/xbj + mp2, // [GeV^2]
-                   nu = 0.5 * ( q2 + mx2 - mp2 ) / mp; // [GeV]
-      const double dm2 = mx2-mp2, Xpr = q2/( q2+mx2 ), En = dm2+q2, Tau = 0.25 * q2/mp2, MQ = params_.rho2+q2;
+      const double inv_q2 = 1./q2;
+      const double dm2 = q2*( 1.-xbj )/xbj;
+      const double mx2 = mp2 + dm2; // [GeV^2]
+      const double en = dm2+q2, en2 = en*en; // [GeV^2]
+      const double nu = 0.5 * en / mp, Xpr = q2/( q2+mx2 ), tau = 0.25 * q2/mp2;
+      const double mq = params_.rho2+q2;
 
       SuriYennie sy;
-      sy.FM = ( params_.C1*dm2*pow( params_.rho2/MQ, 2 ) + ( params_.C2*mp2*pow( 1.-Xpr, 4 ) ) / ( 1.+Xpr*( Xpr*params_.Cp-2.*params_.Bp ) ) )/q2;
-      sy.FE = ( Tau*sy.FM + params_.D1*dm2*q2*params_.rho2/mp2*pow( dm2/MQ/En, 2 ) ) / ( 1.+0.25*En*En/mp2/q2 );
+      sy.FM = ( params_.C1*dm2*pow( params_.rho2/mq, 2 )
+            + ( params_.C2*mp2*pow( 1.-Xpr, 4 ) ) / ( 1.+Xpr*( Xpr*params_.Cp-2.*params_.Bp ) ) ) * inv_q2;
+      sy.FE = ( tau*sy.FM + params_.D1*dm2*q2*params_.rho2/mp2*pow( dm2/mq, 2 )/en2 ) / ( 1.+0.25*en2/mp2*inv_q2 );
 
-      const double w2 = 2.*mp*sy.FE/*, w1 = 0.5 * sy.FM*q2/mp*/;
+      //const double w2 = 2.*mp*sy.FE;
 
-      sy.F2 = nu/mp*w2;
+      sy.F1 = 0.5*sy.FM*q2/mp;
+      sy.F2 = 2.*nu*sy.FE;
       return sy;
     }
   }
 }
diff --git a/CepGen/StructureFunctions/SuriYennie.h b/CepGen/StructureFunctions/SuriYennie.h
index 2c443ad..97d29a8 100644
--- a/CepGen/StructureFunctions/SuriYennie.h
+++ b/CepGen/StructureFunctions/SuriYennie.h
@@ -1,31 +1,31 @@
 #ifndef CepGen_StructureFunctions_SuriYennie_h
 #define CepGen_StructureFunctions_SuriYennie_h
 
 #include "StructureFunctions.h"
 
 namespace CepGen
 {
   namespace SF
   {
     class SuriYennie : public StructureFunctions
     {
       public:
         struct Parameterisation {
           // values extracted from experimental fits
           static Parameterisation standard();
           static Parameterisation alternative();
 
           double C1, C2, D1, rho2, Cp, Bp;
         };
 
         explicit SuriYennie( const SuriYennie::Parameterisation& param = SuriYennie::Parameterisation::standard() );
         SuriYennie operator()( double q2, double xbj ) const;
 
-        double FE, FM;
+        double F1, FE, FM;
       private:
         Parameterisation params_;
     };
   }
 }
 
 #endif
diff --git a/test/cepgen-lhe.cpp b/test/cepgen-lhe.cpp
index 1f8d4db..eef9a72 100644
--- a/test/cepgen-lhe.cpp
+++ b/test/cepgen-lhe.cpp
@@ -1,49 +1,49 @@
-#include <iostream>
-
+#include "CepGen/Cards/PythonHandler.h"
 #include "CepGen/Generator.h"
-#include "CepGen/Cards/PythiaHandler.h"
 #include "CepGen/Export/LHEFHandler.h"
 
 #include "HepMC/Version.h"
 
+#include <iostream>
+
 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;
   
   if ( argc == 1 ) FatalError( "No config file provided." );
 
   Debugging( Form( "Reading config file stored in %s", argv[1] ) );
-  CepGen::Cards::PythiaHandler card( argv[1] );
+  CepGen::Cards::PythonHandler card( argv[1] );
   mg.setParameters( card.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 );
 
   CepGen::OutputHandler::LHEFHandler writer( "example.dat" );
   writer.setCrossSection( xsec, err );
   writer.initialise( *mg.parameters );
 
   // The events generation starts here !
   for ( unsigned int i = 0; i < mg.parameters->generation.maxgen; ++i ) {
     if ( i%1000 == 0 )
       cout << "Generating event #" << i+1 << endl;
     try {
       writer << mg.generateOneEvent().get();
     } catch ( CepGen::Exception& e ) { e.dump(); }
   }
 
   return 0;
 }
 
diff --git a/test/cepgen.cpp b/test/cepgen.cpp
index 09cd624..3ea9aed 100644
--- a/test/cepgen.cpp
+++ b/test/cepgen.cpp
@@ -1,79 +1,79 @@
 #include "CepGen/Cards/PythonHandler.h"
 #include "CepGen/Cards/LpairHandler.h"
 
 #include "CepGen/Generator.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"
 
 #include <iostream>
 
 using namespace std;
 
 void printEvent( const CepGen::Event& ev, unsigned int& ev_id )
 {
   if ( ev_id % 5000 == 0 ) {
     Information( Form( "Generating event #%d", ev_id ) );
     ev.dump();
   }
 }
 
 /**
  * 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;
     mg.parameters->kinematics.inp = { 6500., 6500. };
     mg.parameters->kinematics.central_system = { CepGen::Muon, CepGen::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 );
+    const std::string extension = CepGen::Cards::Handler::getExtension( argv[1] );
     if ( extension == "card" ) mg.setParameters( CepGen::Cards::LpairHandler( argv[1] ).parameters() );
 #ifdef PYTHON
     else if ( extension == "py" ) mg.setParameters( CepGen::Cards::PythonHandler( argv[1] ).parameters() );
 #endif
   }
 
   // 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 !
     mg.generate( printEvent ); // use a callback function!
 
   return 0;
 }
 
diff --git a/test/generate_root_output.cxx b/test/generate_root_output.cxx
index dc8352c..38451fc 100644
--- a/test/generate_root_output.cxx
+++ b/test/generate_root_output.cxx
@@ -1,104 +1,108 @@
-#include <iostream>
+#include "CepGen/Cards/PythonHandler.h"
+#include "CepGen/Cards/LpairHandler.h"
+#include "CepGen/Generator.h"
+#include "CepGen/Event/Event.h"
 
 // ROOT includes
 #include "TFile.h"
 #include "TTree.h"
 #include "TLorentzVector.h"
 
 #include "TreeEvent.h"
 
-#include "CepGen/Generator.h"
-#include "CepGen/Cards/LpairHandler.h"
-#include "CepGen/Cards/PythiaHandler.h"
+#include <iostream>
 
 using namespace std;
 
 /**
  * Generation of events and storage in a ROOT format
  * @author Laurent Forthomme <laurent.forthomme@cern.ch>
  * @date 27 jan 2014
  */
 int main( int argc, char* argv[] ) {
   CepGen::Generator mg;
 
   if ( argc < 2 ) {
     InError( Form( "Usage: %s <input card> [output .root filename]", argv[0] ) );
     return -1;
   }
-  const std::string incard( argv[1] ), extension = incard.substr( incard.find_last_of( "." )+1 );
-  if ( extension == "card" ) mg.setParameters( CepGen::Cards::LpairHandler( argv[1] ).parameters() );
-  else if ( extension == "py" ) mg.setParameters( CepGen::Cards::PythiaHandler( argv[1] ).parameters() );
+  const std::string extension = CepGen::Cards::Handler::getExtension( argv[1] );
+  if ( extension == "card" )
+    mg.setParameters( CepGen::Cards::LpairHandler( argv[1] ).parameters() );
+  else if ( extension == "py" )
+    mg.setParameters( CepGen::Cards::PythonHandler( argv[1] ).parameters() );
 
   mg.parameters->generation.enabled = true;
   mg.parameters->dump();
 
   //----- open the output root file
 
   const TString filename = ( argc > 2 ) ? argv[2] : "events.root";
   auto file = TFile::Open( filename, "recreate" );
   if ( !file ) {
-    cout << "ERROR while trying to create the output file!" << endl;
+    cerr << "ERROR while trying to create the output file!" << endl;
+    return -1;
   }
 
   //----- start by computing the cross section for the list of parameters applied
   double xsec, err;
   mg.computeXsection( xsec, err );
 
   //----- then generate the events and the container tree structure
 
   auto ev_tree = new TTree( "events", "A TTree containing information from the events produced from CepGen" );
 
   CepGen::TreeRun run;
   run.create();
   run.xsect = xsec;
   run.errxsect = err;
   run.litigious_events = 0;
   run.fill();
 
   CepGen::TreeEvent ev;
   ev.create( ev_tree );
 
   for ( unsigned int i = 0; i < mg.parameters->generation.maxgen; ++i ) {
     const auto event = mg.generateOneEvent();
     if ( !event ) FatalError( "Failed to generate the event!" );
 
     ev.clear();
     if ( i % 10000 == 0 ) {
       cout << ">> event " << i << " generated" << endl;
       //event->dump();
     }
 
     ev.gen_time = event->time_generation;
     ev.tot_time = event->time_total;
     ev.np = 0;
     for ( const auto& p : event->particles() ) {
       const CepGen::Particle::Momentum m = p.momentum();
 
       //ev.kinematics[ev.np].SetXYZM( m.px(), m.py(), m.pz(), m.mass() );
       ev.rapidity[ev.np] = m.rapidity();
       ev.pt[ev.np] = m.pt();
       ev.eta[ev.np] = m.eta();
       ev.phi[ev.np] = m.phi();
       ev.E[ev.np] = p.energy();
       ev.m[ev.np] = p.mass();
       ev.pdg_id[ev.np] = p.integerPdgId();
       ev.parent1[ev.np] = ( p.mothers().size() > 0 ) ? *p.mothers().begin() : -1;
       ev.parent2[ev.np] = ( p.mothers().size() > 1 ) ? *p.mothers().rbegin() : -1;
       ev.status[ev.np] = p.status();
       ev.stable[ev.np] = ( (short)p.status() > 0 );
       ev.charge[ev.np] = p.charge();
       ev.role[ev.np] = p.role();
 
       ev.np++;
     }
 
     ev_tree->Fill();
   }
   //cout << "Number of litigious events = " << run.litigious_events << " -> fraction = " << ( run.litigious_events*100./ngen ) << "%" << endl;
 
   file->Write();
   delete file;
 
   return 0;
 }
 
diff --git a/test/produce_distributions.cxx b/test/produce_distributions.cxx
index bd04843..cff0aef 100644
--- a/test/produce_distributions.cxx
+++ b/test/produce_distributions.cxx
@@ -1,59 +1,58 @@
+#include "CepGen/Cards/PythonHandler.h"
 #include "CepGen/Generator.h"
-#include "CepGen/Cards/PythiaHandler.h"
+#include "CepGen/Event/Event.h"
 
 #include "Canvas.h"
 #include "TH1.h"
 
 #include <sstream>
 
 void produce_plot( const char* name, TH1* hist )
 {
   CepGen::Canvas c( name, "CepGen Simulation" );
   hist->Draw( "hist" );
   c.Prettify( hist );
   c.SetLogy();
   c.Save( "pdf" );
 }
 
 int main( int argc, char* argv[] )
 {
   CepGen::Generator mg;
   //CepGen::Logger::get().level = CepGen::Logger::Debug;
 
   if ( argc < 2 ) {
     InError( Form( "Usage: %s [input card]", argv[0] ) );
     return -1;
   }
-  mg.setParameters( CepGen::Cards::PythiaHandler( argv[1] ).parameters() );
+  mg.setParameters( CepGen::Cards::PythonHandler( argv[1] ).parameters() );
 
   TH1D h_mass( "invm", "Dilepton invariant mass\\d#sigma/dM\\GeV?.2f", 1000, 0., 500. ),
        h_ptpair( "ptpair", "Dilepton p_{T}\\d#sigma/dp_{T}\\GeV?.2f", 500, 0., 50. ),
        h_ptsingle( "pt_single", "Single lepton p_{T}\\d#sigma/dp_{T}\\?.2f", 100, 0., 100. ),
        h_etasingle( "eta_single", "Single lepton #eta\\d#sigma/d#eta\\?.2f", 60, -3., 3. );
 
-  std::ostringstream gen_name;
-  gen_name << mg.parameters->process()->name();
-  Information( Form( "Process name: %s", gen_name.str().c_str() ) );
+  Information( Form( "Process name: %s", mg.parameters->processName().c_str() ) );
   //mg.parameters->taming_functions.dump();
 
   for ( unsigned int i = 0; i < mg.parameters->generation.maxgen; ++i ) {
     const CepGen::Event& ev = *mg.generateOneEvent();
     if ( i%100==0 ) Information( Form( "Produced event #%d", i ) );
     const auto central_system = ev.getByRole( CepGen::Particle::CentralSystem );
     const auto pl1 = central_system[0].momentum(), pl2 = central_system[1].momentum();
     h_mass.Fill( ( pl1+pl2 ).mass() );
     h_ptpair.Fill( ( pl1+pl2 ).pt() );
     h_ptsingle.Fill( pl1.pt() );
     h_etasingle.Fill( pl1.eta() );
   }
   const double weight = mg.crossSection()/mg.parameters->generation.maxgen;
   h_mass.Scale( weight );
   h_ptpair.Scale( weight );
 
   produce_plot( "dilepton_invm", &h_mass );
   produce_plot( "dilepton_ptpair", &h_ptpair );
   produce_plot( "singlelepton_pt", &h_ptsingle );
   produce_plot( "singlelepton_eta", &h_etasingle );
 
   return 0;
 }
diff --git a/test/test_cross_section_scan.cpp b/test/test_cross_section_scan.cpp
index f1b538b..924318e 100644
--- a/test/test_cross_section_scan.cpp
+++ b/test/test_cross_section_scan.cpp
@@ -1,47 +1,51 @@
 #include "CepGen/Generator.h"
+#include "CepGen/Parameters.h"
+#include "CepGen/Core/Exception.h"
 #include "CepGen/Processes/GamGamLL.h"
 
+#include <fstream>
+
 using namespace std;
 
 int main( int argc, char* argv[] )
 {
   if ( argc<5 ) {
     InError( Form( "Usage: %s <process mode=1..4> <num points> <min value> <max value> [output file=xsect.dat]", argv[0] ) );
     return -1;
   }
   const unsigned int proc_mode = atoi( argv[1] ),
                      npoints = atoi( argv[2] );
   const float min_value = atof( argv[3] ),
               max_value = atof( argv[4] );
   const char* output_file = ( argc>5 ) ? argv[5] : "xsect.dat";
 
   CepGen::Generator mg;
 
   CepGen::Logger::get().level = CepGen::Logger::Error;
 
   CepGen::Parameters* par = mg.parameters.get();
   par->kinematics.inp = { 6500., 6500. };
   par->kinematics.cuts.central[CepGen::Cuts::eta_single] = { -2.5, 2.5 };
   par->kinematics.cuts.remnants[CepGen::Cuts::mass].max() = 1000.0;
   par->setProcess( new CepGen::Process::GamGamLL );
   par->kinematics.mode = static_cast<CepGen::Kinematics::ProcessMode>( proc_mode );
   par->dump();
 
   double xsect, err_xsect;
 
   ofstream xsect_file( output_file );
   if ( !xsect_file.is_open() ) {
     InError( Form( "Output file \"%s\" cannot be opened!", output_file ) );
     return -2;
   }
 
   for ( unsigned int i=0; i<npoints; i++ ) {
     par->kinematics.cuts.central[CepGen::Cuts::pt_single].min() = min_value + (max_value-min_value)*i/npoints;
     par->kinematics.dump();
     mg.computeXsection( xsect, err_xsect );
     xsect_file << Form( "%.2f\t%.5f\t%.5f\n", par->kinematics.cuts.central[CepGen::Cuts::pt_single].min(), xsect, err_xsect );
     xsect_file.flush();
   }
 
   return 0;
 }
diff --git a/test/test_function_parser.cpp b/test/test_function_parser.cpp
index 82eb2de..bef7933 100644
--- a/test/test_function_parser.cpp
+++ b/test/test_function_parser.cpp
@@ -1,44 +1,45 @@
 #include "CepGen/Core/Functional.h"
 #include "CepGen/Core/Exception.h"
 
+#include <math.h>
 #include <string>
 #include <iostream>
 
 using namespace std;
 
 int
 main()
 {
   const double epsilon = 1.e-9; // tolerance
 
   { // test with a 1-variable function
     const double exp_result_test1 = 6.795704571;
     CepGen::Functional<1> test( "2.5*exp(0.1*x)", { { "x" } } );
     if ( fabs( test.eval( 10. ) - exp_result_test1 ) > epsilon  ) throw CepGen::Exception( __PRETTY_FUNCTION__, "Test 1.1 failed!", CepGen::FatalError );
     if ( fabs( test.eval( { { 10. } } ) - exp_result_test1 ) > epsilon  ) throw CepGen::Exception( __PRETTY_FUNCTION__, "Test 1.2 failed!", CepGen::FatalError );
     cout << "Test 1 passed!" << endl;
   }
   { // test with an invalid function
     bool passed = false;
     try {
       CepGen::Functional<1> test( "sqrt(x+x**3-log(10)", { { "x" } } );
       test.eval( 10 );
     } catch ( CepGen::Exception& e ) { cout << "Test 2 passed!" /*<< e.what()*/ << endl; passed = true; }
     if ( !passed ) throw CepGen::Exception( __PRETTY_FUNCTION__, "Test 2 failed!", CepGen::FatalError );
   }
   { // test with a 2-variables function
     CepGen::Functional<2> test( "sqrt(a^2+b^2)", { { "a", "b" } } );
     if ( fabs( test.eval( { { 3, 4 } } ) - 5.0 ) > epsilon  ) throw CepGen::Exception( __PRETTY_FUNCTION__, "Test 3 failed!", CepGen::FatalError );
     cout << "Test 3 passed!" << endl;
   }
   { // test with an invalid function
     CepGen::Functional<1> test( "a**2", { { "a" } } );
     bool passed = true;
     try { test.eval( 10 ); passed = false; } catch ( CepGen::Exception& e ) {}
     try { test.eval( { { 10 } } ); passed = false; } catch ( CepGen::Exception& e ) {}
     if ( !passed  ) throw CepGen::Exception( __PRETTY_FUNCTION__, "Test 4 failed!", CepGen::FatalError );
     cout << "Test 4 passed!" << endl;
   }
 
   return 0;
 }