diff --git a/CepGen/Event/EventBrowser.cpp b/CepGen/Event/EventBrowser.cpp
index ffe269b..365f499 100644
--- a/CepGen/Event/EventBrowser.cpp
+++ b/CepGen/Event/EventBrowser.cpp
@@ -1,86 +1,96 @@
 #include "CepGen/Event/EventBrowser.h"
 #include "CepGen/Event/Event.h"
 
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/utils.h"
 
 namespace cepgen
 {
   namespace utils
   {
     const std::regex EventBrowser::rgx_select_id_( "(\\w+)\\((\\d+)\\)" );
     const std::regex EventBrowser::rgx_select_role_( "(\\w+)\\(([a-z]+\\d?)\\)" );
 
     double
     EventBrowser::get( const Event& ev, const std::string& var ) const
     {
       std::smatch sm;
       //--- particle-level variables (indexed by integer id)
       if ( std::regex_match( var, sm, rgx_select_id_ ) ) {
         const auto& var_name = sm[1].str();
         const auto& part = ev[std::stod( sm[2].str() )];
-        return variable( part, var_name );
+        return variable( ev, part, var_name );
       }
       //--- particle-level variables (indexed by role)
       else if ( std::regex_match( var, sm, rgx_select_role_ ) ) {
         const auto& var_name = sm[1].str();
         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 INVALID_OUTPUT;
         }
         const auto& part = ev[role_str_.at( str_role )][0];
-        return variable( part, var_name );
+        return variable( ev, part, var_name );
       }
       //--- event-level variables
       else
         return variable( ev, var );
     }
 
     double
-    EventBrowser::variable( const Particle& part, const std::string& var ) const
+    EventBrowser::variable( const Event& ev, 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 == "xi"  ) {
+        const auto& moth = part.mothers();
+        if ( moth.empty() ) {
+          CG_WARNING( "EventBrowser" )
+            << "Failed to retrieve parent particle to compute xi "
+            << "for the following particle:\n"
+            << part;
+          return INVALID_OUTPUT;
+        }
+        return 1.-part.energy()/ev[*moth.begin()].energy();
+      }
       if ( var == "pdg" ) return (double)part.integerPdgId();
       if ( var == "charge" ) return part.charge();
       if ( var == "status" ) return (double)part.status();
       CG_WARNING( "EventBrowser" )
         << "Failed to retrieve variable \"" << var << "\".";
       return INVALID_OUTPUT;
     }
 
     double
     EventBrowser::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( "EventBrowser" )
         << "Failed to retrieve the event-level variable \"" << var << "\".";
       return INVALID_OUTPUT;
     }
   }
 }
 
diff --git a/CepGen/Event/EventBrowser.h b/CepGen/Event/EventBrowser.h
index ccf5e93..0366bdf 100644
--- a/CepGen/Event/EventBrowser.h
+++ b/CepGen/Event/EventBrowser.h
@@ -1,61 +1,61 @@
 #ifndef CepGen_Event_EventBrowser_h
 #define CepGen_Event_EventBrowser_h
 
 #include "CepGen/Event/Particle.h"
 #include <regex>
 
 namespace cepgen
 {
   class Event;
   namespace utils
   {
     /**
      * \brief A user-friendly browser for the Event content
      * \author Laurent Forthomme <laurent.forthomme@cern.ch>
      * \date Jul 2019
      */
     class EventBrowser
     {
       public:
         EventBrowser() = default;
         double get( const Event& ev, const std::string& var ) const;
 
       private:
         /// Retrieve a named variable from a particle
-        double variable( const Particle&, const std::string& ) const;
+        double variable( const Event&, 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.;
 
         //--- 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 }
         };
     };
   }
 }
 
 #endif
 
diff --git a/CepGen/Event/Particle.cpp b/CepGen/Event/Particle.cpp
index 4dd99a0..ca67386 100644
--- a/CepGen/Event/Particle.cpp
+++ b/CepGen/Event/Particle.cpp
@@ -1,286 +1,284 @@
 #include "CepGen/Event/Particle.h"
 
 #include "CepGen/Physics/PDG.h"
 
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/utils.h"
 
 #include <iomanip>
 
 namespace cepgen
 {
   Particle::Particle() :
     id_( -1 ), charge_sign_( 1 ),
     mass_( -1. ), helicity_( 0. ),
     role_( UnknownRole ), status_( Status::Undefined ), pdg_id_( PDG::invalid )
   {}
 
   Particle::Particle( Role role, pdgid_t pdgId, Status st ) :
     id_( -1 ), charge_sign_( 1 ),
     mass_( -1. ), helicity_( 0. ),
     role_( role ), status_( st ),
     pdg_id_( pdgId )
   {
     try {
       phys_prop_ = PDG::get()( pdg_id_ );
     } catch ( const Exception& ) {}
     if ( pdg_id_ != PDG::invalid )
       computeMass();
   }
 
   Particle::Particle( const Particle& part ) :
     id_( part.id_ ), charge_sign_( part.charge_sign_ ),
     momentum_( part.momentum_ ), mass_( part.mass_ ), helicity_( part.helicity_ ),
     role_( part.role_ ), status_( part.status_ ),
     mothers_( part.mothers_ ), daughters_( part.daughters_ ),
     pdg_id_( part.pdg_id_ )
   {
     try {
       phys_prop_ = PDG::get()( pdg_id_ );
     } catch ( const Exception& ) {}
   }
 
   bool
   Particle::operator<( const Particle& rhs ) const
   {
     return id_ >= 0
         && rhs.id_ > 0
         && id_ < rhs.id_;
   }
 
   double
   Particle::thetaToEta( double theta )
   {
     return -log( tan( 0.5*theta*M_PI/180. ) );
   }
 
   double
   Particle::etaToTheta( double eta )
   {
     return 2.*atan( exp( -eta ) )*180.*M_1_PI;
   }
 
   bool
   Particle::valid()
   {
     if ( pdg_id_ == PDG::invalid )
       return false;
     if ( momentum_.p() == 0. && mass_ == 0. )
       return false;
     return true;
   }
 
   float
   Particle::charge() const
   {
     return charge_sign_ * phys_prop_.charge/3.;
   }
 
   void
   Particle::computeMass( bool off_shell )
   {
     if ( !off_shell && pdg_id_ != PDG::invalid ) { // retrieve the mass from the on-shell particle's properties
       mass_ = phys_prop_.mass;
     }
     else if ( momentum_.energy() >= 0. ) {
       mass_ = sqrt( energy2() - momentum_.p2() );
     }
 
     //--- finish by setting the energy accordingly
     if ( momentum_.energy() < 0. ) { // invalid energy
       momentum_.setEnergy( sqrt( momentum_.p2() + mass2() ) );
     }
   }
 
   void
   Particle::setMass( double m )
   {
     if ( m >= 0. )
       mass_ = m;
     else
       computeMass();
   }
 
   void
   Particle::addMother( Particle& part )
   {
     mothers_.insert( part.id() );
 
     CG_DEBUG_LOOP( "Particle" )
       <<  "Particle " << id() << " (pdgId=" << part.integerPdgId() << ") "
       << "is the new mother of " << id_ << " (pdgId=" << (int)pdg_id_ << ").";
 
     part.addDaughter( *this );
   }
 
   void
   Particle::addDaughter( Particle& part )
   {
     const auto ret = daughters_.insert( part.id() );
 
     if ( CG_LOG_MATCH( "Particle", debugInsideLoop ) ) {
       std::ostringstream os;
       for ( const auto& daugh : daughters_ )
         os << Form( "\n\t * id=%d", daugh );
       CG_DEBUG_LOOP( "Particle" )
         << "Particle " << role_ << " (pdgId=" << (int)pdg_id_ << ") "
         << "has now " << daughters_.size() << " daughter(s):"
         << os.str();
     }
 
     if ( ret.second ) {
       CG_DEBUG_LOOP( "Particle" )
         << "Particle " << part.role() << " (pdgId=" << part.integerPdgId() << ") "
         << "is a new daughter of " << role_ << " (pdgId=" << pdg_id_ << ").";
 
       if ( part.mothers().find( id_ ) == part.mothers().end() )
         part.addMother( *this );
     }
   }
 
   void
   Particle::setMomentum( const Momentum& mom, bool offshell )
   {
     momentum_ = mom;
     if ( !offshell && mom.mass() > 0. )
       mass_ = momentum_.mass();
     else
       computeMass();
   }
 
   void
   Particle::setMomentum( double px, double py, double pz, double e )
   {
     momentum_.setP( px, py, pz );
     setEnergy( e );
     if ( fabs( e-momentum_.energy() ) > 1.e-6 ) // more than 1 eV difference
       CG_WARNING( "Particle" ) << "Energy difference: " << e-momentum_.energy();
   }
 
   double
   Particle::energy() const
   {
     return ( momentum_.energy() < 0.
       ? std::hypot( mass_, momentum_.p() )
       : momentum_.energy() );
   }
 
   void
   Particle::setEnergy( double e )
   {
     if ( e < 0. && mass_ >= 0. )
       e = std::hypot( mass_, momentum_.p() );
     momentum_.setEnergy( e );
   }
 
   pdgid_t
   Particle::pdgId() const
   {
     return pdg_id_;
   }
 
   void
   Particle::setPdgId( short pdg )
   {
     pdg_id_ = abs( pdg );
     try {
       phys_prop_ = PDG::get()( pdg_id_ );
     } catch ( const Exception& ) {}
     switch ( pdg_id_ ) {
       case PDG::electron: case PDG::muon: case PDG::tau:
         charge_sign_ = -pdg/abs( pdg ); break;
       default:
         charge_sign_ = pdg/abs( pdg ); break;
     }
   }
 
   void
   Particle::setPdgId( pdgid_t pdg, short ch )
   {
     pdg_id_ = pdg;
     switch ( pdg_id_ ) {
       case PDG::electron: case PDG::muon: case PDG::tau:
         charge_sign_ = -ch; break;
       default:
         charge_sign_ = ch; break;
     }
   }
 
   int
   Particle::integerPdgId() const
   {
     const float ch = phys_prop_.charge/3.;
     if ( ch == 0 )
       return static_cast<int>( pdg_id_ );
     return static_cast<int>( pdg_id_ ) * charge_sign_ * ( ch/fabs( ch ) );
   }
 
-  void
-  Particle::dump() const
+  double
+  Particle::etaToY( double eta_, double m_, double pt_ )
+  {
+    const double m2 = m_*m_, mt = std::hypot( m_, pt_ );
+    return asinh( sqrt( ( ( mt*mt-m2 )*cosh( 2.*eta_ )+m2 )/ mt*mt - 1. )*M_SQRT1_2 );
+  }
+
+  std::ostream&
+  operator<<( std::ostream& os, const Particle& part )
   {
-    std::ostringstream os;
     os << std::resetiosflags( std::ios::showbase )
-      << "Particle[" << id_ << "]{role=" << role_ << ", status=" << (int)status_ << ", "
-      << "pdg=" << pdg_id_ << ", p4=" << momentum_ << " GeV, m=" << mass_ << " GeV, "
-      << "p⟂=" << momentum_.pt() << " GeV, eta=" << momentum_.eta() << ", phi=" << momentum_.phi();
-    if ( primary() )
+      << "Particle[" << part.id_ << "]{role=" << part.role_ << ", status=" << (int)part.status_ << ", "
+      << "pdg=" << part.pdg_id_ << ", p4=" << part.momentum_ << " GeV, m=" << part.mass_ << " GeV, "
+      << "p⟂=" << part.momentum_.pt() << " GeV, eta=" << part.momentum_.eta() << ", phi=" << part.momentum_.phi();
+    if ( part.primary() )
       os << ", primary";
     else {
-      os << ", mother" << utils::s( mothers_.size() ) << "=";
+      os << ", " << utils::s( "mother", part.mothers_.size() ) << "=";
       std::string delim;
-      for ( const auto& moth : mothers_ )
+      for ( const auto& moth : part.mothers_ )
         os << delim << moth, delim = ",";
     }
-    const auto& daughters_list = daughters();
+    const auto& daughters_list = part.daughters();
     if ( !daughters_list.empty() ) {
-      os << ", daughter" << utils::s( daughters_list.size() ) << "=";
+      os << ", " << utils::s( "daughter", daughters_list.size() ) << "=";
       std::string delim;
       for ( const auto& daugh : daughters_list )
         os << delim << daugh, delim = ",";
     }
-    os << "}";
-    CG_INFO( "Particle" ) << os.str();
-  }
-
-  double
-  Particle::etaToY( double eta_, double m_, double pt_ )
-  {
-    const double m2 = m_*m_, mt = std::hypot( m_, pt_ );
-    return asinh( sqrt( ( ( mt*mt-m2 )*cosh( 2.*eta_ )+m2 )/ mt*mt - 1. )*M_SQRT1_2 );
+    return os << "}";
   }
 
   std::ostream&
   operator<<( std::ostream& os, const Particle::Role& rl )
   {
     switch ( rl ) {
       case Particle::UnknownRole:   return os << "unknown";
       case Particle::IncomingBeam1: return os << "i.beam 1";
       case Particle::IncomingBeam2: return os << "i.beam 2";
       case Particle::OutgoingBeam1: return os << "o.beam 1";
       case Particle::OutgoingBeam2: return os << "o.beam 2";
       case Particle::Parton1:       return os << "parton 1";
       case Particle::Parton2:       return os << "parton 2";
       case Particle::Intermediate:  return os << "hard pr.";
       case Particle::CentralSystem: return os << "central";
     }
     return os;
   }
 
   double
   CMEnergy( const Particle& p1, const Particle& p2 )
   {
     if ( p1.mass()*p2.mass() < 0.
       || p1.energy()*p2.energy() < 0. )
       return 0.;
     return sqrt( p1.mass2()+p2.mass2() + 2.*p1.energy()*p2.energy() - 2.*( p1.momentum()*p2.momentum() ) );
   }
 
   double
   CMEnergy( const Particle::Momentum& m1, const Particle::Momentum& m2 )
   {
     if ( m1.mass()*m2.mass() < 0.
       || m1.energy()*m2.energy() < 0. )
       return 0.;
     return sqrt( m1.mass2()+m2.mass2() + 2.*m1.energy()*m2.energy() - 2.*( m1*m2 ) );
   }
 }
diff --git a/CepGen/Event/Particle.h b/CepGen/Event/Particle.h
index b7188dc..75c4ecb 100644
--- a/CepGen/Event/Particle.h
+++ b/CepGen/Event/Particle.h
@@ -1,355 +1,355 @@
 #ifndef CepGen_Event_Particle_h
 #define CepGen_Event_Particle_h
 
 #include "CepGen/Core/Hasher.h"
 
 #include "CepGen/Physics/Constants.h"
 #include "CepGen/Physics/ParticleProperties.h"
 
 #include <set>
 #include <unordered_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 class Status {
         PrimordialIncoming = -9, ///< Incoming beam particle
         DebugResonance = -5, ///< Intermediate resonance (for processes developers)
         Resonance = -4, ///< Already decayed intermediate resonance
         Fragmented = -3, ///< Already fragmented outgoing beam
         Propagator = -2, ///< Generic propagator
         Incoming = -1, ///< Incoming parton
         Undefined = 0, ///< Undefined particle
         FinalState = 1, ///< Stable, final state particle
         Undecayed = 2, ///< Particle to be decayed externally
         Unfragmented = 3 ///< Particle to be hadronised externally
       };
       /// Role of the particle in the process
       enum Role {
         UnknownRole = -1, ///< Undefined role
         IncomingBeam1 = 1, ///< \f$z>0\f$ incoming beam particle
         IncomingBeam2 = 2, ///< \f$z<0\f$ incoming beam particle
         OutgoingBeam1 = 3, ///< \f$z<0\f$ outgoing beam state/particle
         OutgoingBeam2 = 5, ///< \f$z>0\f$ outgoing beam state/particle
         CentralSystem = 6, ///< Central particles system
         Intermediate = 4, ///< Intermediate two-parton system
         Parton1 = 41, ///< \f$z>0\f$ beam incoming parton
         Parton2 = 42 ///< \f$z<0\f$ beam incoming parton
       };
       /**
        * 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. );
           /// Build a 4-momentum using its 3-momentum coordinates and its energy
           Momentum( double* p );
 
           // --- 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 );
           /// Build a 4-momentum from its three momentum coordinates and mass
           static Momentum fromPxPyPzM( double px, double py, double pz, double m );
           /// Build a 4-momentum from its transverse momentum, rapidity and mass
           static Momentum fromPxPyYM( double px, double py, double rap, double m );
 
           // --- 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;
           /// Vector product of the 3-momentum with another 3-momentum
           double crossProduct( const Momentum& ) const;
           /// Compute the 4-vector sum of two 4-momenta
           Momentum operator+( const Momentum& ) const;
           /// Add a 4-momentum through a 4-vector sum
           Momentum& operator+=( const Momentum& );
           /// Unary inverse operator
           Momentum operator-() const;
           /// Compute the inverse per-coordinate 4-vector
           Momentum operator-( const Momentum& ) const;
           /// Subtract a 4-momentum through a 4-vector sum
           Momentum& operator-=( const Momentum& );
           /// Scalar product of two 3-momenta
           double operator*( const Momentum& ) const;
           /// Scalar product of the 3-momentum with another 3-momentum
           double operator*=( const Momentum& );
           /// Multiply all components of a 4-momentum by a scalar
           Momentum operator*( double c ) const;
           /// Multiply all 4-momentum coordinates by a scalar
           Momentum& operator*=( double c );
           /// Left-multiply all 4-momentum coordinates by a scalar
           friend Momentum operator*( double, const Momentum& );
           /// Equality operator
           bool operator==( const Momentum& ) const;
           /// Human-readable format for a particle's momentum
           friend std::ostream& operator<<( std::ostream&, const Momentum& );
 
           Momentum& betaGammaBoost( double gamma, double betagamma );
           /// Forward Lorentz boost
           Momentum& lorentzBoost( const Momentum& p );
 
           // --- setters and getters
 
           /// Set all the components of the 4-momentum (in GeV)
           Momentum& setP( double px, double py, double pz, double e );
           /// Set all the components of the 3-momentum (in GeV)
           Momentum& setP( double px, double py, double pz );
           /// Set the energy (in GeV)
           inline Momentum& setEnergy( double e ) { energy_ = e; return *this; }
           /// Compute the energy from the mass
           inline Momentum& setMass( double m ) { return setMass2( m*m ); }
           /// Compute the energy from the mass
           Momentum& 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$^2\f$)
           double pt2() const;
           /// 5-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$^2\f$)
           inline double p2() const { return p_*p_; }
           /// Energy (in GeV)
           inline double energy() const { return energy_; }
           /// Squared energy (in GeV\f$^2\f$)
           inline double energy2() const { return energy_*energy_; }
           /// Squared mass (in GeV\f$^2\f$) 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;
           Momentum& truncate( double tolerance = 1.e-10 );
           /// 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
           Momentum& 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 role in the event
       friend std::ostream& operator<<( std::ostream& os, const Role& rl );
 
       //----- 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] role Role of the particle in the process
       /// \param[in] id PDG identifier
       /// \param[in] st Current status
       Particle( Role role, pdgid_t id, Status st = Status::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 ); }
 
       // --- 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;
       /// 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 PDG identifier
       /// \param[in] ch Electric charge (0, 1, or -1)
       void setPdgId( pdgid_t pdg, short ch = 0 );
       /// Set the PDG identifier (along with the particle's electric charge)
       /// \param[in] pdg_id PDG identifier (incl. electric charge in e)
       void setPdgId( short pdg_id );
       /// Retrieve the objectified PDG identifier
       pdgid_t pdgId() const;
       /// 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; }
       /// Particle mass in GeV/c\f$^2\f$
       /// \return Particle's mass
       inline double mass() const { return mass_; };
       /// Compute the particle mass
       /// \param[in] off_shell Allow the particle to be produced off-shell?
       /// \note This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority)
       void computeMass( bool off_shell = false );
       /// Set the particle mass, in GeV/c\f$^2\f$
       /// \param m Mass in GeV/c\f$^2\f$
       /// \note This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority)
       void setMass( double m = -1. );
       /// Particle squared mass, in GeV\f$^2\f$/c\f$^4\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- or 4-momentum associated to the particle
        * \param[in] px Momentum along the \f$x\f$-axis, in GeV/c
        * \param[in] py Momentum along the \f$y\f$-axis, in GeV/c
        * \param[in] pz Momentum along the \f$z\f$-axis, in GeV/c
        * \param[in] e Energy, in GeV
        */
       void setMomentum( double px, double py, double pz, double e = -1. );
       /// 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] ); }
       /// 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 GeV\f$^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(); }
       /// Set the mother particle
       /// \param[in] part A Particle object containing all the information on the mother particle
       void addMother( Particle& part );
       /// Get 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(); };
       /// 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;
+      friend std::ostream& operator<<( std::ostream&, const Particle& );
 
     protected:
       /// Unique identifier in an event
       int id_;
       /// Electric charge (+-1 or 0)
       short charge_sign_;
       /// Momentum properties handler
       Momentum momentum_;
       /// Mass, in GeV/c\f$^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
       pdgid_t pdg_id_;
       /// Collection of standard, bare-level physical properties
       ParticleProperties phys_prop_;
   };
 
   /// 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::unordered_map<Particle::Role,Particles,utils::EnumHash<Particle::Role> > ParticlesMap;
 }
 
 #endif