diff --git a/CepGen/Cards/PythonHandler.cpp b/CepGen/Cards/PythonHandler.cpp
index dab6e8d..29d9b9c 100644
--- a/CepGen/Cards/PythonHandler.cpp
+++ b/CepGen/Cards/PythonHandler.cpp
@@ -1,392 +1,393 @@
 #include "CepGen/Cards/PythonHandler.h"
 #include "CepGen/Core/Exception.h"
 
 #include "CepGen/Core/TamingFunction.h"
 #include "CepGen/Core/ParametersList.h"
 #include "CepGen/Core/Integrator.h"
 
 #include "CepGen/Processes/ProcessesHandler.h"
 #include "CepGen/Hadronisers/HadronisersHandler.h"
 #include "CepGen/IO/ExportHandler.h"
 #include "CepGen/StructureFunctions/StructureFunctions.h"
 
 #include "CepGen/Physics/GluonGrid.h"
 #include "CepGen/Physics/PDG.h"
 
 #include <algorithm>
 
 #if PY_MAJOR_VERSION < 3
 #  define PYTHON2
 #endif
 
 namespace cepgen
 {
   namespace card
   {
     //----- specialization for CepGen input cards
     PythonHandler::PythonHandler( const char* file )
     {
       setenv( "PYTHONPATH", ".:Cards:test:../Cards", 1 );
       setenv( "PYTHONDONTWRITEBYTECODE", "1", 1 );
       CG_DEBUG( "PythonHandler" )
         << "Python PATH: " << getenv( "PYTHONPATH" ) << ".";
       std::string filename = pythonPath( file );
       const size_t fn_len = filename.length()+1;
 
       //Py_DebugFlag = 1;
       //Py_VerboseFlag = 1;
 
 #ifdef PYTHON2
       char* sfilename = new char[fn_len];
       snprintf( sfilename, fn_len, "%s", filename.c_str() );
 #else
       wchar_t* sfilename = new wchar_t[fn_len];
       swprintf( sfilename, fn_len, L"%s", filename.c_str() );
 #endif
       if ( sfilename )
         Py_SetProgramName( sfilename );
 
       Py_InitializeEx( 1 );
 
       if ( sfilename )
         delete [] sfilename;
       if ( !Py_IsInitialized() )
         throw CG_FATAL( "PythonHandler" ) << "Failed to initialise the Python cards parser!";
 
       CG_DEBUG( "PythonHandler" )
         << "Initialised the Python cards parser\n\t"
         << "Python version: " << Py_GetVersion() << "\n\t"
         << "Platform: " << Py_GetPlatform() << ".";
 
       PyObject* cfg = PyImport_ImportModule( filename.c_str() ); // new
       if ( !cfg )
         throwPythonError( Form( "Failed to parse the configuration card %s", file ) );
 
       //--- additional particles definition
       PyObject* pextp = PyObject_GetAttrString( cfg, PDGLIST_NAME ); // new
       if ( pextp ) {
         parseExtraParticles( pextp );
         Py_CLEAR( pextp );
       }
 
       //--- process definition
       PyObject* process = PyObject_GetAttrString( cfg, PROCESS_NAME ); // new
       if ( !process )
         throwPythonError( Form( "Failed to extract a \"%s\" keyword from the configuration card %s", PROCESS_NAME, file ) );
 
       //--- list of process-specific parameters
       ParametersList proc_params;
       fillParameter( process, "processParameters", proc_params );
 
       //--- type of process to consider
       PyObject* pproc_name = element( process, MODULE_NAME ); // borrowed
       if ( !pproc_name )
         throwPythonError( Form( "Failed to extract the process name from the configuration card %s", file ) );
       const std::string proc_name = get<std::string>( pproc_name );
 
       //--- process mode
       params_.kinematics.mode = (KinematicsMode)proc_params.get<int>( "mode", (int)KinematicsMode::invalid );
       params_.setProcess( cepgen::proc::ProcessesHandler::get().build( proc_name, proc_params ) );
 
       //--- process kinematics
       PyObject* pin_kinematics = element( process, "inKinematics" ); // borrowed
       if ( pin_kinematics )
         parseIncomingKinematics( pin_kinematics );
 
       PyObject* pout_kinematics = element( process, "outKinematics" ); // borrowed
       if ( pout_kinematics )
         parseOutgoingKinematics( pout_kinematics );
 
       //--- taming functions
       PyObject* ptam = element( process, "tamingFunctions" ); // borrowed
       if ( ptam )
         for ( const auto& p : getVector<ParametersList>( ptam ) )
           params_.taming_functions->add( p.get<std::string>( "variable" ), p.get<std::string>( "expression" ) );
 
       Py_CLEAR( process );
 
       PyObject* plog = PyObject_GetAttrString( cfg, LOGGER_NAME ); // new
       if ( plog ) {
         parseLogging( plog );
         Py_CLEAR( plog );
       }
 
       //--- hadroniser parameters
       PyObject* phad = PyObject_GetAttrString( cfg, HADR_NAME ); // new
       if ( phad ) {
         parseHadroniser( phad );
         Py_CLEAR( phad );
       }
 
       //--- generation parameters
       PyObject* pint = PyObject_GetAttrString( cfg, INTEGRATOR_NAME ); // new
       if ( pint ) {
         parseIntegrator( pint );
         Py_CLEAR( pint );
       }
 
       PyObject* pgen = PyObject_GetAttrString( cfg, GENERATOR_NAME ); // new
       if ( pgen ) {
         parseGenerator( pgen );
         Py_CLEAR( pgen );
       }
 
       PyObject* pout = PyObject_GetAttrString( cfg, OUTPUT_NAME ); // new
       if ( pout ) {
         parseOutputModule( pout );
         Py_CLEAR( pout );
       }
 
       //--- finalisation
       Py_CLEAR( cfg );
     }
 
     PythonHandler::~PythonHandler()
     {
       if ( Py_IsInitialized() )
         Py_Finalize();
     }
 
     void
     PythonHandler::parseIncomingKinematics( PyObject* kin )
     {
       //--- retrieve the beams PDG ids
       std::vector<ParametersList> beams_pdg;
       fillParameter( kin, "pdgIds", beams_pdg );
       if ( !beams_pdg.empty() ) {
         if ( beams_pdg.size() != 2 )
           throwPythonError( Form( "Invalid list of PDG ids retrieved for incoming beams:\n\t2 PDG ids are expected, %d provided!", beams_pdg.size() ) );
         params_.kinematics.incoming_beams. first.pdg = (pdgid_t)beams_pdg.at( 0 ).get<int>( "pdgid" );
         params_.kinematics.incoming_beams.second.pdg = (pdgid_t)beams_pdg.at( 1 ).get<int>( "pdgid" );
       }
       //--- incoming beams kinematics
       std::vector<double> beams_pz;
       fillParameter( kin, "pz", beams_pz );
       if ( !beams_pz.empty() ) {
         if ( beams_pz.size() != 2 )
           throwPythonError( Form( "Invalid list of pz's retrieved for incoming beams:\n\t2 pz's are expected, %d provided!", beams_pz.size() ) );
         params_.kinematics.incoming_beams. first.pz = beams_pz.at( 0 );
         params_.kinematics.incoming_beams.second.pz = beams_pz.at( 1 );
       }
       double sqrt_s = -1.;
       fillParameter( kin, "cmEnergy", sqrt_s );
       if ( sqrt_s != -1. )
         params_.kinematics.setSqrtS( sqrt_s );
       //--- structure functions set for incoming beams
       PyObject* psf = element( kin, "structureFunctions" ); // borrowed
       if ( psf )
         params_.kinematics.structure_functions = strfun::StructureFunctionsHandler::get().build( get<ParametersList>( psf ) );
       //--- types of parton fluxes for kt-factorisation
       std::vector<int> kt_fluxes;
       fillParameter( kin, "ktFluxes", kt_fluxes );
       if ( !kt_fluxes.empty() ) {
         params_.kinematics.incoming_beams.first.kt_flux = (KTFlux)kt_fluxes.at( 0 );
         params_.kinematics.incoming_beams.second.kt_flux = ( kt_fluxes.size() > 1 )
           ? (KTFlux)kt_fluxes.at( 1 )
           : (KTFlux)kt_fluxes.at( 0 );
       }
       //--- specify where to look for the grid path for gluon emission
       std::string kmr_grid_path;
       fillParameter( kin, "kmrGridPath", kmr_grid_path );
       if ( !kmr_grid_path.empty() )
         kmr::GluonGrid::get( kmr_grid_path.c_str() );
       //--- parse heavy ions beams
       std::vector<int> hi_beam1, hi_beam2;
       fillParameter( kin, "heavyIonA", hi_beam1 );
       if ( hi_beam1.size() == 2 )
         params_.kinematics.incoming_beams. first.pdg = HeavyIon{ (unsigned short)hi_beam1[0], (Element)hi_beam1[1] };
       fillParameter( kin, "heavyIonB", hi_beam2 );
       if ( hi_beam2.size() == 2 )
         params_.kinematics.incoming_beams.second.pdg = HeavyIon{ (unsigned short)hi_beam2[0], (Element)hi_beam2[1] };
     }
 
     void
     PythonHandler::parseOutgoingKinematics( PyObject* kin )
     {
       std::vector<int> parts;
       fillParameter( kin, "minFinalState", parts );
       for ( const auto& pdg : parts )
         params_.kinematics.minimum_final_state.emplace_back( (pdgid_t)pdg );
 
       ParametersList part_cuts;
       fillParameter( kin, "cuts", part_cuts );
       for ( const auto& part : part_cuts.keys() ) {
         const auto pdg = (pdgid_t)stoi( part );
         const auto& cuts = part_cuts.get<ParametersList>( part );
         if ( cuts.has<Limits>( "pt" ) )
           params_.kinematics.cuts.central_particles[pdg].pt_single = cuts.get<Limits>( "pt" );
         if ( cuts.has<Limits>( "energy" ) )
           params_.kinematics.cuts.central_particles[pdg].energy_single = cuts.get<Limits>( "energy" );
         if ( cuts.has<Limits>( "eta" ) )
           params_.kinematics.cuts.central_particles[pdg].eta_single = cuts.get<Limits>( "eta" );
         if ( cuts.has<Limits>( "rapidity" ) )
           params_.kinematics.cuts.central_particles[pdg].rapidity_single = cuts.get<Limits>( "rapidity" );
       }
 
       // for LPAIR/collinear matrix elements
       fillParameter( kin, "q2", params_.kinematics.cuts.initial.q2 );
 
       // for the kT factorised matrix elements
       fillParameter( kin, "qt", params_.kinematics.cuts.initial.qt );
       fillParameter( kin, "phiqt", params_.kinematics.cuts.initial.phi_qt );
       fillParameter( kin, "ptdiff", params_.kinematics.cuts.central.pt_diff );
       fillParameter( kin, "phiptdiff", params_.kinematics.cuts.central.phi_pt_diff );
       fillParameter( kin, "rapiditydiff", params_.kinematics.cuts.central.rapidity_diff );
 
       // generic phase space limits
       fillParameter( kin, "rapidity", params_.kinematics.cuts.central.rapidity_single );
       fillParameter( kin, "eta", params_.kinematics.cuts.central.eta_single );
       fillParameter( kin, "pt", params_.kinematics.cuts.central.pt_single );
 
       fillParameter( kin, "ptsum", params_.kinematics.cuts.central.pt_sum );
       fillParameter( kin, "invmass", params_.kinematics.cuts.central.mass_sum );
 
       fillParameter( kin, "mx", params_.kinematics.cuts.remnants.mass_single );
       fillParameter( kin, "yj", params_.kinematics.cuts.remnants.rapidity_single );
 
       Limits lim_xi;
       fillParameter( kin, "xi", lim_xi );
       if ( lim_xi.valid() )
-        params_.kinematics.cuts.remnants.energy_single = ( lim_xi+(-1.) )*( -params_.kinematics.incoming_beams.first.pz );
+        //params_.kinematics.cuts.remnants.energy_single = ( lim_xi+(-1.) )*( -params_.kinematics.incoming_beams.first.pz );
+        params_.kinematics.cuts.remnants.energy_single = -( lim_xi-1. )*params_.kinematics.incoming_beams.first.pz;
     }
 
     void
     PythonHandler::parseLogging( PyObject* log )
     {
       int log_level = 0;
       fillParameter( log, "level", log_level );
       utils::Logger::get().level = (utils::Logger::Level)log_level;
       std::vector<std::string> enabled_modules;
       fillParameter( log, "enabledModules", enabled_modules );
       for ( const auto& mod : enabled_modules )
         utils::Logger::get().addExceptionRule( mod );
     }
 
     void
     PythonHandler::parseIntegrator( PyObject* integr )
     {
       if ( !PyDict_Check( integr ) )
         throwPythonError( "Integrator object should be a dictionary!" );
       PyObject* palgo = element( integr, MODULE_NAME ); // borrowed
       if ( !palgo )
         throwPythonError( "Failed to retrieve the integration algorithm name!" );
       std::string algo = get<std::string>( palgo );
       if ( algo == "plain" )
         params_.integration().type = IntegratorType::plain;
       else if ( algo == "Vegas" ) {
         params_.integration().type = IntegratorType::Vegas;
         fillParameter( integr, "alpha", (double&)params_.integration().vegas.alpha );
         fillParameter( integr, "iterations", params_.integration().vegas.iterations );
         fillParameter( integr, "mode", (int&)params_.integration().vegas.mode );
         fillParameter( integr, "verbosity", (int&)params_.integration().vegas.verbose );
         std::string vegas_logging_output = "cerr";
         fillParameter( integr, "loggingOutput", vegas_logging_output );
         if ( vegas_logging_output == "cerr" )
           // redirect all debugging information to the error stream
           params_.integration().vegas.ostream = stderr;
         else if ( vegas_logging_output == "cout" )
           // redirect all debugging information to the standard stream
           params_.integration().vegas.ostream = stdout;
         else
           params_.integration().vegas.ostream = fopen( vegas_logging_output.c_str(), "w" );
       }
       else if ( algo == "MISER" ) {
         params_.integration().type = IntegratorType::MISER;
         fillParameter( integr, "estimateFraction", (double&)params_.integration().miser.estimate_frac );
         fillParameter( integr, "minCalls", params_.integration().miser.min_calls );
         fillParameter( integr, "minCallsPerBisection", params_.integration().miser.min_calls_per_bisection );
         fillParameter( integr, "alpha", (double&)params_.integration().miser.alpha );
         fillParameter( integr, "dither", (double&)params_.integration().miser.dither );
       }
       else
         throwPythonError( Form( "Invalid integration() algorithm: %s", algo.c_str() ) );
 
       fillParameter( integr, "numFunctionCalls", params_.integration().ncvg );
       fillParameter( integr, "seed", (unsigned long&)params_.integration().rng_seed );
       unsigned int rng_engine;
       fillParameter( integr, "rngEngine", rng_engine );
       switch ( rng_engine ) {
         case 0: default: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_mt19937; break;
         case 1: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_taus2; break;
         case 2: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_gfsr4; break;
         case 3: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_ranlxs0; break;
       }
       fillParameter( integr, "chiSqCut", params_.integration().vegas_chisq_cut );
     }
 
     void
     PythonHandler::parseGenerator( PyObject* gen )
     {
       if ( !PyDict_Check( gen ) )
         throwPythonError( "Generation information object should be a dictionary!" );
       params_.generation().enabled = true;
       fillParameter( gen, "treat", params_.generation().treat );
       fillParameter( gen, "numEvents", params_.generation().maxgen );
       fillParameter( gen, "printEvery", params_.generation().gen_print_every );
       fillParameter( gen, "numThreads", params_.generation().num_threads );
       fillParameter( gen, "numPoints", params_.generation().num_points );
     }
 
     void
     PythonHandler::parseHadroniser( PyObject* hadr )
     {
       if ( !PyDict_Check( hadr ) )
         throwPythonError( "Hadroniser object should be a dictionary!" );
 
       PyObject* pname = element( hadr, MODULE_NAME ); // borrowed
       if ( !pname )
         throwPythonError( "Hadroniser name is required!" );
       std::string hadr_name = get<std::string>( pname );
 
       params_.setHadroniser( cepgen::hadr::HadronisersHandler::get().build( hadr_name, get<ParametersList>( hadr ) ) );
 
       auto h = params_.hadroniser();
       h->setParameters( params_ );
       { //--- before calling the init() method
         std::vector<std::string> config;
         fillParameter( hadr, "preConfiguration", config );
         h->readStrings( config );
       }
       h->init();
       { //--- after init() has been called
         std::vector<std::string> config;
         fillParameter( hadr, "processConfiguration", config );
         for ( const auto& block : config ) {
           std::vector<std::string> config_blk;
           fillParameter( hadr, block.c_str(), config_blk );
           h->readStrings( config_blk );
         }
       }
     }
 
     void
     PythonHandler::parseOutputModule( PyObject* pout )
     {
       if ( !is<ParametersList>( pout ) )
         throwPythonError( "Invalid type for output parameters list!" );
 
       PyObject* pname = element( pout, MODULE_NAME ); // borrowed
       if ( !pname )
         throwPythonError( "Output module name is required!" );
       params_.setOutputModule( io::ExportHandler::get().build( get<std::string>( pname ), get<ParametersList>( pout ) ) );
     }
 
     void
     PythonHandler::parseExtraParticles( PyObject* pparts )
     {
       if ( !is<ParametersList>( pparts ) )
         throwPythonError( "Extra particles definition object should be a parameters list!" );
 
       const auto& parts = get<ParametersList>( pparts );
       for ( const auto& k : parts.keys() ) {
         const auto& part = parts.get<ParticleProperties>( k );
         if ( part.pdgid == 0 || part.mass < 0. )
           continue;
         CG_DEBUG( "PythonHandler:particles" )
           << "Adding a new particle with name \"" << part.name << "\" to the PDG dictionary.";
         PDG::get().define( part );
       }
     }
   }
 }
 
diff --git a/CepGen/Core/Integrand.cpp b/CepGen/Core/Integrand.cpp
index 0fae045..4109c12 100644
--- a/CepGen/Core/Integrand.cpp
+++ b/CepGen/Core/Integrand.cpp
@@ -1,202 +1,202 @@
 #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/Physics/PDG.h"
 
 #include "CepGen/Processes/GenericProcess.h"
 #include "CepGen/Hadronisers/GenericHadroniser.h"
 #include "CepGen/IO/GenericExportHandler.h"
 
 #include "CepGen/Parameters.h"
 
 #include <cmath>
 #include <sstream>
 #include <fstream>
 
 namespace cepgen
 {
   namespace integrand
   {
     utils::Logger::Level log_level;
     utils::Timer tmr;
 
     double
     eval( double* x, size_t ndim, void* func_params )
     {
       log_level = utils::Logger::get().level;
       std::shared_ptr<Event> ev;
 
       Parameters* params = nullptr;
       if ( !func_params || !( params = static_cast<Parameters*>( func_params ) ) )
         throw CG_FATAL( "Integrand" ) << "Failed to retrieve the run parameters!";
 
       proc::GenericProcess* proc = params->process();
       if ( !proc )
         throw CG_FATAL( "Integrand" ) << "Failed to retrieve the process!";
 
       //================================================================
       // start the timer
       //================================================================
 
       tmr.reset();
 
       //================================================================
       // prepare the event content prior to the process generation
       //================================================================
 
       if ( proc->hasEvent() ) // event is not empty
         ev = proc->event();
 
       params->prepareRun();
 
       //================================================================
       // specify the phase space point to probe
       //================================================================
 
       proc->setPoint( ndim, x );
 
       //================================================================
       // from this step on, the phase space point is supposed to be set
       //================================================================
 
       proc->beforeComputeWeight();
       double weight = proc->computeWeight();
 
       //================================================================
       // invalidate any unphysical behaviour
       //================================================================
 
       if ( weight <= 0. )
         return 0.;
 
       //================================================================
       // speed up the integration process if no event is to be generated
       //================================================================
 
       if ( !ev )
         return weight;
 
       if ( !params->storage()
         && !params->taming_functions
         && !params->hadroniser()
         &&  params->kinematics.cuts.central_particles.empty() )
         return weight;
 
       //================================================================
       // fill in the process' Event object
       //================================================================
 
       proc->fillKinematics();
 
       //================================================================
       // once the kinematics variables have been populated, can apply
       // the collection of taming functions
       //================================================================
 
       if ( params->taming_functions ) {
         if ( params->taming_functions->has( "m_central" )
           || params->taming_functions->has( "pt_central" ) ) {
 
           // build the kinematics of the central system
           Particle::Momentum central_system;
           for ( const auto& part : (*ev)[Particle::CentralSystem] )
             central_system += part.momentum();
 
           // tame the cross-section by the reweighting function
           if ( params->taming_functions->has( "m_central" ) )
             weight *= params->taming_functions->eval( "m_central", central_system.mass() );
           if ( params->taming_functions->has( "pt_central" ) )
             weight *= params->taming_functions->eval( "pt_central", central_system.pt() );
         }
         if ( params->taming_functions->has( "q2" ) ) {
           weight *= params->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton1 ).momentum().mass() );
           weight *= params->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton2 ).momentum().mass() );
         }
       }
 
       if ( weight <= 0. )
         return 0.;
 
       //================================================================
       // set the CepGen part of the event generation
       //================================================================
 
       if ( params->storage() )
         ev->time_generation = tmr.elapsed();
 
       //================================================================
       // event hadronisation and resonances decay
       //================================================================
 
       if ( params->hadroniser() ) {
         double br = -1.;
         if ( !params->hadroniser()->run( *ev, br, params->storage() ) || br == 0. )
           return 0.;
         weight *= br; // branching fraction for all decays
       }
 
       //================================================================
       // apply cuts on final state system (after hadronisation!)
       // (polish your cuts, as this might be very time-consuming...)
       //================================================================
 
       if ( !params->kinematics.cuts.central_particles.empty() ) {
         for ( const auto& part : (*ev)[Particle::CentralSystem] ) {
           // retrieve all cuts associated to this final state particle
           if ( params->kinematics.cuts.central_particles.count( part.pdgId() ) == 0 )
             continue;
           const auto& cuts_pdgid = params->kinematics.cuts.central_particles.at( part.pdgId() );
           // apply these cuts on the given particle
           if ( !cuts_pdgid.pt_single.passes( part.momentum().pt() ) )
             return 0.;
           if ( !cuts_pdgid.energy_single.passes( part.momentum().energy() ) )
             return 0.;
           if ( !cuts_pdgid.eta_single.passes( part.momentum().eta() ) )
             return 0.;
           if ( !cuts_pdgid.rapidity_single.passes( part.momentum().rapidity() ) )
             return 0.;
         }
       }
       for ( const auto& system : { Particle::OutgoingBeam1, Particle::OutgoingBeam2 } )
         for ( const auto& part : (*ev)[system] ) {
           if ( part.status() != Particle::Status::FinalState )
             continue;
-          if ( !params->kinematics.cuts.remnants.energy_single.passes( fabs( part.momentum().energy() ) ) )
+          if ( !params->kinematics.cuts.remnants.energy_single.passes( part.momentum().energy() ) )
             return 0.;
           if ( !params->kinematics.cuts.remnants.rapidity_single.passes( fabs( part.momentum().rapidity() ) ) )
             return 0.;
         }
 
       //================================================================
       // store the last event in parameters block for a later usage
       //================================================================
 
       if ( params->storage() ) {
         proc->last_event = ev;
         proc->last_event->time_total = tmr.elapsed();
 
         CG_DEBUG( "Integrand" )
           << "[process 0x" << std::hex << proc << std::dec << "] "
           << "Individual time (gen+hadr+cuts): " << proc->last_event->time_total*1.e3 << " ms";
       }
 
       //================================================================
       // a bit of useful debugging
       //================================================================
 
       if ( CG_LOG_MATCH( "Integrand", debugInsideLoop ) ) {
         std::ostringstream oss;
         for ( unsigned short i = 0; i < ndim; ++i )
           oss << Form( "%10.8f ", x[i] );
         CG_DEBUG( "Integrand" )
           << "f value for dim-" << ndim << " point ( " << oss.str() << "): "
           << weight;
       }
 
       return weight;
     }
   }
 }
diff --git a/CepGen/IO/TextHandler.cpp b/CepGen/IO/TextHandler.cpp
index e4031da..f0a2e6f 100644
--- a/CepGen/IO/TextHandler.cpp
+++ b/CepGen/IO/TextHandler.cpp
@@ -1,330 +1,331 @@
 #include "CepGen/IO/ExportHandler.h"
 
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/ParametersList.h"
 #include "CepGen/Core/utils.h"
 
 #include "CepGen/Event/Event.h"
 #include "CepGen/Parameters.h"
 
 #include "CepGen/Version.h"
 
 #include <gsl/gsl_histogram.h>
 
 #include <iomanip>
 #include <fstream>
 #include <regex>
 
 namespace cepgen
 {
   namespace io
   {
     /**
      * \brief Handler for the generic text file output
      * \author Laurent Forthomme <laurent.forthomme@cern.ch>
      * \date Jul 2019
      */
     class TextHandler : public GenericExportHandler
     {
       public:
         explicit TextHandler( const ParametersList& );
         ~TextHandler();
 
         void initialise( const Parameters& ) override;
         void setCrossSection( double xsec, double ) override { xsec_ = xsec; }
         void operator<<( const Event& ) override;
 
       private:
         short extractVariableProperties( const std::string& );
         std::string textHistogram( const std::string&, const gsl_histogram* ) const;
         /// Retrieve a named variable from a particle
         double variable( const Particle&, const std::string& ) const;
         /// Retrieve a named variable from the whole event
         double variable( const Event&, const std::string& ) const;
 
         static const std::regex rgx_select_id_, rgx_select_role_;
         static constexpr double INVALID_OUTPUT = -999.;
         static constexpr size_t PLOT_WIDTH = 50;
         static constexpr char PLOT_CHAR = '#';
 
         std::ofstream file_, hist_file_;
         const std::vector<std::string> variables_;
         const bool save_banner_, save_variables_;
         const bool show_hists_, save_hists_;
         const ParametersList hist_variables_;
         const std::string separator_;
 
         //--- variables definition
         std::unordered_map<short,std::string> variables_name_;
         std::unordered_map<short,bool> variable_stored_;
 
         typedef std::pair<unsigned short,std::string> IndexedVariable;
         std::unordered_map<short,std::vector<IndexedVariable> > variables_per_id_;
         std::unordered_map<Particle::Role,std::vector<IndexedVariable> > variables_per_role_;
         std::vector<IndexedVariable> variables_for_event_;
         unsigned short num_vars_;
         std::ostringstream oss_vars_;
 
         double xsec_;
 
         //--- auxiliary helper maps
         const std::unordered_map<std::string,Particle::Role> role_str_ = {
           { "ib1", Particle::Role::IncomingBeam1 }, { "ib2", Particle::Role::IncomingBeam2 },
           { "ob1", Particle::Role::OutgoingBeam1 }, { "ob2", Particle::Role::OutgoingBeam2 },
           { "pa1", Particle::Role::Parton1 }, { "pa2", Particle::Role::Parton2 },
           { "cs",  Particle::Role::CentralSystem },
           { "int", Particle::Role::Intermediate }
         };
         typedef double( Particle::Momentum::*pMethod )(void) const;
         /// Mapping of string variables to momentum getter methods
         const std::unordered_map<std::string,pMethod> m_mom_str_ = {
           { "px",  &Particle::Momentum::px },
           { "py",  &Particle::Momentum::py },
           { "pz",  &Particle::Momentum::pz },
           { "pt",  &Particle::Momentum::pt },
           { "eta", &Particle::Momentum::eta },
           { "phi", &Particle::Momentum::phi },
           { "m",   &Particle::Momentum::mass },
           { "e",   &Particle::Momentum::energy },
           { "p",   &Particle::Momentum::p },
           { "pt2", &Particle::Momentum::pt2 },
           { "th",  &Particle::Momentum::theta },
           { "y",   &Particle::Momentum::rapidity }
         };
 
         //--- kinematic variables
         double sqrts_;
         unsigned long num_evts_;
         struct gsl_histogram_deleter
         {
           void operator()( gsl_histogram* h ) {
             gsl_histogram_free( h );
           }
         };
         std::unordered_map<short,std::unique_ptr<gsl_histogram,gsl_histogram_deleter> > hists_;
     };
 
     const std::regex TextHandler::rgx_select_id_( "(\\w+)\\((\\d+)\\)" );
     const std::regex TextHandler::rgx_select_role_( "(\\w+)\\(([a-z]+\\d?)\\)" );
 
     TextHandler::TextHandler( const ParametersList& params ) :
       GenericExportHandler( "text" ),
       file_          ( params.get<std::string>( "filename", "output.txt" ) ),
       variables_     ( params.get<std::vector<std::string> >( "variables" ) ),
       save_banner_   ( params.get<bool>( "saveBanner", true ) ),
       save_variables_( params.get<bool>( "saveVariables", true ) ),
       show_hists_    ( params.get<bool>( "showHistograms", true ) ),
       save_hists_    ( params.get<bool>( "saveHistograms", false ) ),
       hist_variables_( params.get<ParametersList>( "histVariables" ) ),
       separator_     ( params.get<std::string>( "separator", "\t" ) ),
       num_vars_( 0 ), xsec_( 1. )
     {
       //--- first extract list of variables to store in output file
       oss_vars_.clear();
       std::string sep;
       for ( const auto& var : variables_ ) {
         auto id = extractVariableProperties( var );
         if ( id >= 0 ) {
           oss_vars_ << sep << var, sep = separator_;
           variable_stored_[id] = true;
         }
       }
       //--- then extract list of variables to be plotted in histogram
       for ( const auto& var : hist_variables_.keys() ) {
         auto id = extractVariableProperties( var );
         if ( id < 0 )
           continue;
         const auto& hvar = hist_variables_.get<ParametersList>( var );
         const int nbins = hvar.get<int>( "nbins", 10 );
         const double min = hvar.get<double>( "low", 0. ), max = hvar.get<double>( "high", 1. );
         hists_[id].reset( gsl_histogram_alloc( nbins ) );
         gsl_histogram_set_ranges_uniform( hists_[id].get(), min, max );
         CG_INFO( "TextHandler" )
           << "Booking a histogram with " << nbins << " bin" << utils::s( nbins )
           << " between " << min << " and " << max << " for \"" << var << "\".";
       }
       if ( save_hists_ && !hists_.empty() )
         hist_file_.open( "lastrun.hists.txt" );
     }
 
     TextHandler::~TextHandler()
     {
       //--- histograms printout
       for ( const auto& var : hist_variables_.keys() ) {
         const auto& vn = std::find_if( variables_name_.begin(), variables_name_.end(),
           [&var]( auto&& p ) { return p.second == var; } );
         if ( vn == variables_name_.end() ) {
           CG_WARNING( "TextHandler" )
             << "Failed to retrieve variable \"" << var << "\" for plotting.";
           continue;
         }
         const auto& hist = hists_.at( vn->first ).get();
         gsl_histogram_scale( hist, xsec_/( num_evts_+1 ) );
         if ( show_hists_ )
           CG_INFO( "TextHandler" )
             << textHistogram( var, hist );
         if ( save_hists_ )
           hist_file_ << "\n" << textHistogram( var, hist ) << "\n";
       }
       //--- finalisation of the output file
       file_.close();
     }
 
     void
     TextHandler::initialise( const Parameters& params )
     {
       sqrts_ = params.kinematics.sqrtS();
       num_evts_ = 0ul;
       if ( save_banner_ )
         file_ << banner( params, "#" ) << "\n";
       if ( save_variables_ )
         file_ << "# " << oss_vars_.str() << "\n";
       if ( save_hists_ && !hists_.empty() )
         hist_file_ << banner( params, "#" ) << "\n";
     }
 
     void
     TextHandler::operator<<( const Event& ev )
     {
       std::vector<double> vars( num_vars_ );
       //--- extract and order the variables to be retrieved
       //--- particle-level variables (indexed by integer id)
       for ( const auto& id_vars : variables_per_id_ ) {
         const auto& part = ev[id_vars.first];
         //--- loop over the list of variables for this particle
         for ( const auto& var : id_vars.second )
           vars[var.first] = variable( part, var.second );
       }
       //--- particle-level variables (indexed by role)
       for ( const auto& role_vars : variables_per_role_ ) {
         const auto& part = ev[role_vars.first][0];
         //--- loop over the list of variables for this particle
         for ( const auto& var : role_vars.second )
           vars[var.first] = variable( part, var.second );
       }
       //--- event-level variables
       for ( const auto& var : variables_for_event_ )
         vars[var.first] = variable( ev, var.second );
       //--- write down the variables list in the file
       std::string sep;
       unsigned short i = 0;
       for ( const auto& var : vars ) {
         if ( variable_stored_.count( i ) > 0 && variable_stored_.at( i ) )
           file_ << sep << var, sep = separator_;
         if ( hists_.count( i ) > 0 )
           gsl_histogram_increment( hists_.at( i ).get(), var );
         ++i;
       }
       file_ << "\n";
       ++num_evts_;
     }
 
     double
     TextHandler::variable( const Particle& part, const std::string& var ) const
     {
       if ( m_mom_str_.count( var ) ) {
         auto meth = m_mom_str_.at( var );
         return ( part.momentum().*meth )();
       }
       if ( var == "xi"  ) return 1.-part.momentum().energy()*2./sqrts_;
       if ( var == "pdg" ) return (double)part.integerPdgId();
       if ( var == "charge" ) return part.charge();
       if ( var == "status" ) return (double)part.status();
       CG_WARNING( "TextHandler" )
         << "Failed to retrieve variable \"" << var << "\".";
       return INVALID_OUTPUT;
     }
 
     double
     TextHandler::variable( const Event& ev, const std::string& var ) const
     {
       if ( var == "np" )
         return (double)ev.size();
       if ( var == "nev" )
         return (double)num_evts_+1;
       if ( var == "nob1" || var == "nob2" ) {
         unsigned short out = 0.;
         for ( const auto& part : ev[
           var == "nob1"
           ? Particle::Role::OutgoingBeam1
           : Particle::Role::OutgoingBeam2
         ] )
           if ( (int)part.status() > 0 )
             out++;
         return (double)out;
       }
       if ( var == "tgen" )
         return ev.time_generation;
       if ( var == "ttot" )
         return ev.time_total;
       CG_WARNING( "TextHandler" )
         << "Failed to retrieve the event-level variable \"" << var << "\".";
       return INVALID_OUTPUT;
     }
 
     short
     TextHandler::extractVariableProperties( const std::string& var )
     {
       const auto& vn = std::find_if( variables_name_.begin(), variables_name_.end(),
         [&var]( auto&& p ) { return p.second == var; } );
       if ( vn != variables_name_.end() )
         return vn->first;
       std::smatch sm;
       if ( std::regex_match( var, sm, rgx_select_id_ ) )
         variables_per_id_[std::stod( sm[2].str() )].emplace_back( std::make_pair( num_vars_, sm[1].str() ) );
       else if ( std::regex_match( var, sm, rgx_select_role_ ) ) {
         const auto& str_role = sm[2].str();
         if ( role_str_.count( str_role ) == 0 ) {
           CG_WARNING( "TextHandler" )
             << "Invalid particle role retrieved from configuration: \"" << str_role << "\".\n\t"
             << "Skipping the variable \"" << var << "\" in the output module.";
           return -1;
         }
         variables_per_role_[role_str_.at( str_role )].emplace_back( std::make_pair( num_vars_, sm[1].str() ) );
       }
       else // event-level variables
         variables_for_event_.emplace_back( std::make_pair( num_vars_, var ) );
       variables_name_[num_vars_] = var;
       return num_vars_++;
     }
 
     std::string
     TextHandler::textHistogram( const std::string& var, const gsl_histogram* hist ) const
     {
       std::ostringstream os;
       const size_t nbins = gsl_histogram_bins( hist );
       const double max_bin = gsl_histogram_max_val( hist );
       const double inv_max_bin = max_bin > 0. ? 1./max_bin : 0.;
-      const std::string sep( 15, ' ' );
+      const std::string sep( 17, ' ' );
       os
         << "plot of \"" << var << "\"\n"
-        << sep << std::string( PLOT_WIDTH-16-var.size(), ' ' )
+        << sep << std::string( PLOT_WIDTH-15-var.size(), ' ' )
         << "d(sig)/d" << var << " (pb/bin)\n"
         << sep << Form( "%-5.2f", gsl_histogram_min_val( hist ) )
-        << std::string( PLOT_WIDTH-9, ' ' )
+        << std::string( PLOT_WIDTH-8, ' ' )
         << Form( "%5.2f", gsl_histogram_max_val( hist ) ) << "\n"
         << sep << std::string( PLOT_WIDTH+2, '.' ); // abscissa axis
       for ( size_t i = 0; i < nbins; ++i ) {
         double min, max;
         gsl_histogram_get_range( hist, i, &min, &max );
         const double value = gsl_histogram_get( hist, i );
         const int val = value*PLOT_WIDTH*inv_max_bin;
         os
-          << "\n" << Form( "[%6.2f,%6.2f):", min, max )
+          << "\n" << Form( "[%7.2f,%7.2f):", min, max )
           << std::string( val, PLOT_CHAR ) << std::string( PLOT_WIDTH-val, ' ' )
           << ": " << Form( "%6.2f", value );
       }
+      const double bin_width = ( gsl_histogram_max( hist )-gsl_histogram_min( hist ) )/nbins;
       os
         << "\n"
         << Form( "%15s", var.c_str() ) << ":" << std::string( PLOT_WIDTH, '.' ) << ":\n" // 2nd abscissa axis
         << "\t("
-        << "bin width=" << ( gsl_histogram_max( hist )-gsl_histogram_min( hist ) )/nbins << ", "
+        << "bin width=" << bin_width << " unit" << utils::s( (int)bin_width ) << ", "
         << "mean=" << gsl_histogram_mean( hist ) << ", "
         << "st.dev.=" << gsl_histogram_sigma( hist )
         << ")";
       return os.str();
     }
   }
 }
 
 REGISTER_IO_MODULE( text, TextHandler )
diff --git a/CepGen/Physics/Limits.cpp b/CepGen/Physics/Limits.cpp
index 5ea5805..af80887 100644
--- a/CepGen/Physics/Limits.cpp
+++ b/CepGen/Physics/Limits.cpp
@@ -1,132 +1,154 @@
 #include "CepGen/Physics/Limits.h"
 
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/utils.h"
 
 namespace cepgen
 {
   Limits::Limits( double min, double max ) :
     std::pair<double,double>( min, max )
   {}
 
   Limits::Limits( const Limits& rhs ) :
     std::pair<double,double>( rhs.first, rhs.second )
   {}
 
+  Limits
+  Limits::operator-() const
+  {
+    return Limits( -second, -first );
+  }
+
   Limits&
   Limits::operator+=( double c )
   {
     first += c;
     second += c;
     return *this;
   }
 
   Limits&
+  Limits::operator-=( double c )
+  {
+    first -= c;
+    second -= c;
+    return *this;
+  }
+
+  Limits&
   Limits::operator*=( double c )
   {
     first *= c;
     second *= c;
     if ( c < 0. ) {
       const double tmp = first; // will be optimised by compiler anyway...
       first = second;
       second = tmp;
     }
     return *this;
   }
 
   void
   Limits::in( double low, double up )
   {
     first = low;
     second = up;
   }
 
   double
   Limits::range() const
   {
     if ( !hasMin() || !hasMax() )
       return 0.;
     return second-first;
   }
 
   bool
   Limits::hasMin() const
   {
     return first != INVALID;
   }
 
   bool
   Limits::hasMax() const
   {
     return second != INVALID;
   }
 
   bool
   Limits::passes( double val ) const
   {
     if ( hasMin() && val < min() )
       return false;
     if ( hasMax() && val > max() )
       return false;
     return true;
   }
 
   bool
   Limits::valid() const
   {
     return hasMin() || hasMax();
   }
 
   void
   Limits::save( bool& on, double& lmin, double& lmax ) const
   {
     on = false; lmin = lmax = 0.;
     if ( !valid() )
       return;
     on = true;
     if ( hasMin() )
       lmin = min();
     if ( hasMax() )
       lmax = max();
     if ( lmin == lmax )
       on = false;
   }
 
   double
   Limits::x( double v ) const
   {
     if ( v < 0. || v > 1. )
       throw CG_ERROR( "Limits:shoot" )
         << "x must be comprised between 0 and 1; x value = " << v << ".";
     if ( !valid() )
       return INVALID;
 
     return first + ( second-first ) * v;
   }
 
   std::ostream&
   operator<<( std::ostream& os, const Limits& lim )
   {
     if ( !lim.hasMin() && !lim.hasMax() )
       return os << "no cuts";
     if ( !lim.hasMin() )
       return os << Form( "below %g", lim.max() );
     if ( !lim.hasMax() )
       return os << Form( "above %g", lim.min() );
     return os << Form( "%g to %g", lim.min(), lim.max() );
   }
+
   Limits
   operator+( Limits lim, double c )
   {
     lim += c;
     return lim;
   }
 
   Limits
+  operator-( Limits lim, double c )
+  {
+    lim -= c;
+    return lim;
+  }
+
+  Limits
   operator*( Limits lim, double c )
   {
     lim *= c;
     return lim;
   }
 
 }
diff --git a/CepGen/Physics/Limits.h b/CepGen/Physics/Limits.h
index 7c982b7..1ec246a 100644
--- a/CepGen/Physics/Limits.h
+++ b/CepGen/Physics/Limits.h
@@ -1,56 +1,59 @@
 #ifndef CepGen_Physics_Limits_h
 #define CepGen_Physics_Limits_h
 
 #include <utility>
 #include <ostream>
 
 namespace cepgen
 {
   /// Validity interval for a variable
   class Limits : private std::pair<double,double>
   {
     public:
       /// Define lower and upper limits on a quantity
       Limits( double min = INVALID, double max = INVALID );
       Limits( const Limits& );
 
+      Limits operator-() const;
       Limits& operator+=( double c );
+      Limits& operator-=( double c );
       Limits& operator*=( double c );
       friend Limits operator+( Limits lim, double c );
+      friend Limits operator-( Limits lim, double c );
       friend Limits operator*( Limits lim, double c );
 
       /// Lower limit to apply on the variable
       double min() const { return first; }
       /// Lower limit to apply on the variable
       double& min() { return first; }
       /// Upper limit to apply on the variable
       double max() const { return second; }
       /// Upper limit to apply on the variable
       double& max() { return second; }
       /// Export the limits into external variables
       void save( bool& on, double& lmin, double& lmax ) const;
       /// Find the [0,1] value scaled between minimum and maximum
       double x( double v ) const;
       /// Specify the lower and upper limits on the variable
       void in( double low, double up );
       /// Full variable range allowed
       double range() const;
       /// Have a lower limit?
       bool hasMin() const;
       /// Have an upper limit?
       bool hasMax() const;
       /// Check if the value is inside limits' boundaries
       bool passes( double val ) const;
       /// Is there a lower and upper limit?
       bool valid() const;
 
       /// Human-readable expression of the limits
       friend std::ostream& operator<<( std::ostream&, const Limits& );
 
     private:
       static constexpr double INVALID = -999.999;
   };
 }
 
 #endif