diff --git a/Cards/pythia8.cfi b/Cards/pythia8.cfi
index f6c7267..c70632a 100644
--- a/Cards/pythia8.cfi
+++ b/Cards/pythia8.cfi
@@ -1,4 +1,28 @@
 hadroniser = {
   name = "pythia8";
-  decayWto = ( 11, 13 );
+  seed = 0L;
+  pythiaConfiguration = (
+    "Random:setSeed = on", // random engine
+    "Tune:preferLHAPDF = 2",
+    "Main:timesAllowErrors = 10000",
+    "Check:epTolErr = 0.01",
+    "Beams:setProductionScalesFromLHEF = off",
+    "SLHA:keepSM = on",
+    "SLHA:minMassSM = 1000.",
+    "ParticleDecays:limitTau0 = on",
+    "ParticleDecays:tau0Max = 10",
+    "ParticleDecays:allowPhotonRadiation = on",
+    //"ParticleDecays:allowPhotonRadiation = off", //FIXME
+    // CUEP8M1 tuning
+    "Tune:pp 14",
+    "Tune:ee 7",
+    "MultipartonInteractions:pT0Ref = 2.4024",
+    "MultipartonInteractions:ecmPow = 0.25208",
+    "MultipartonInteractions:expPow = 1.6",
+    // process-specific
+    "ProcessLevel:all = off",
+    "ProcessLevel:resonanceDecays = on",
+    "24:onMode = off", // disable all W decays, but...
+    "24:onIfAny = 11 13" // enable e-nue + mu-numu final states
+  );
 };
diff --git a/CepGen/Cards/ConfigHandler.cpp b/CepGen/Cards/ConfigHandler.cpp
index 20cdf6a..3d704bf 100644
--- a/CepGen/Cards/ConfigHandler.cpp
+++ b/CepGen/Cards/ConfigHandler.cpp
@@ -1,319 +1,304 @@
 #include "ConfigHandler.h"
 #include "CepGen/Core/Exception.h"
 
 #include "CepGen/Processes/GamGamLL.h"
 #include "CepGen/Processes/PPtoLL.h"
 #include "CepGen/Processes/PPtoWW.h"
 
 #include "CepGen/Hadronisers/Pythia8Hadroniser.h"
 
 #ifdef LIBCONFIG
 
 namespace CepGen
 {
   namespace Cards
   {
     //----- specialization for CepGen input cards
     ConfigHandler::ConfigHandler( const char* file )
     {
 #ifdef LIBCONFIG
       libconfig::Config cfg;
       try { cfg.readFile( file ); } catch ( const libconfig::ParseException& pe ) {
         FatalError( Form( "Failed to parse the configuration card \"%s\".\n\tParser error: %s (at line %d)", file, pe.getError(), pe.getLine() ) );
       }
       try {
         const libconfig::Setting& root = cfg.getRoot();
         const libconfig::Setting& proc = root["process"];
 
         //--- type of process to consider
         const char* proc_name = proc["name"]; const std::string str_proc_name = proc_name;
         if ( str_proc_name == "lpair" ) params_.setProcess( new Process::GamGamLL );
         else if ( str_proc_name == "pptoll" ) params_.setProcess( new Process::PPtoLL );
         else if ( str_proc_name == "pptoww" ) params_.setProcess( new Process::PPtoWW );
         else FatalError( Form( "Unrecognised process: %s", proc_name ) );
 
         //--- process mode
         int int_mode; const char* mode;
         if ( proc.lookupValue( "mode", int_mode ) ) {
           params_.kinematics.mode = (Kinematics::ProcessMode)int_mode;
         }
         else if ( proc.lookupValue( "mode", mode ) ) {
           const std::string str_mode = mode;
           if ( str_mode == "elastic/elastic" ) params_.kinematics.mode = Kinematics::ProcessMode::ElasticElastic;
           else if ( str_mode == "elastic/inelastic" ) params_.kinematics.mode = Kinematics::ProcessMode::ElasticInelastic;
           else if ( str_mode == "inelastic/elastic" ) params_.kinematics.mode = Kinematics::ProcessMode::InelasticElastic;
           else if ( str_mode == "inelastic/inelastic" ) params_.kinematics.mode = Kinematics::ProcessMode::InelasticInelastic;
           else FatalError( Form( "Unrecognised interaction mode: %s", mode ) );
         }
 
         //--- process kinematics
         if ( proc.exists( "in_kinematics" ) ) parseIncomingKinematics( proc["in_kinematics"] );
         if ( proc.exists( "out_kinematics" ) ) parseOutgoingKinematics( proc["out_kinematics"] );
 
         //--- hadroniser parameters
         if ( root.exists( "hadroniser" ) ) parseHadroniser( root["hadroniser"] );
 
         //--- generation parameters
         if ( root.exists( "integrator" ) ) parseIntegrator( root["integrator"] );
         if ( root.exists( "vegas" ) ) parseIntegrator( root["vegas"] ); //FIXME for backward compatibility
         if ( root.exists( "generator" ) ) parseGenerator( root["generator"] );
 
         //--- taming functions
         if ( proc.exists( "taming_functions" ) ) parseTamingFunctions( proc["taming_functions"] );
 
       } catch ( const libconfig::SettingNotFoundException& nfe ) {
         FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
       } catch ( const libconfig::SettingTypeException& te ) {
         FatalError( Form( "Field \"%s\" has wrong type.", te.getPath() ) );
       }
 #else
       InWarning( "libconfig++ is not present on this machine" );
 #endif
     }
 
 #ifdef LIBCONFIG
     void
     ConfigHandler::parseIncomingKinematics( const libconfig::Setting& kin )
     {
       try {
         if ( kin.exists( "beam1_pz" ) ) params_.kinematics.inp.first = (double)kin["beam1_pz"];
         if ( kin.exists( "beam2_pz" ) ) params_.kinematics.inp.second = (double)kin["beam2_pz"];
         if ( kin.exists( "structure_functions" ) ) {
           const char* sf = kin["structure_functions" ]; const std::string sf_str = sf;
           if ( sf_str == "electron" ) params_.kinematics.structure_functions = StructureFunctions::Electron;
           else if ( sf_str == "elastic proton" ) params_.kinematics.structure_functions = StructureFunctions::ElasticProton;
           else if ( sf_str == "Suri-Yennie" ) params_.kinematics.structure_functions = StructureFunctions::SuriYennie;
           else if ( sf_str == "Szczurek-Uleshchenko" ) params_.kinematics.structure_functions = StructureFunctions::SzczurekUleshchenko;
           else if ( sf_str == "Fiore" ) params_.kinematics.structure_functions = StructureFunctions::FioreBrasse;
           else if ( sf_str == "ALLM" ) params_.kinematics.structure_functions = StructureFunctions::ALLM97;
           else if ( sf_str == "ALLM;91" ) params_.kinematics.structure_functions = StructureFunctions::ALLM91;
           else if ( sf_str == "ALLM;97" ) params_.kinematics.structure_functions = StructureFunctions::ALLM97;
           //else if ( sf_str == "ALLM;HHT" ) params_.kinematics.structure_functions = StructureFunctions::ALLM_HHT;
           //else if ( sf_str == "ALLM;HHT-FT" ) params_.kinematics.structure_functions = StructureFunctions::ALLM_HHT_FT;
           else if ( sf_str == "ALLM;GD07p" ) params_.kinematics.structure_functions = StructureFunctions::GD07p;
           else if ( sf_str == "ALLM;GD11p" ) params_.kinematics.structure_functions = StructureFunctions::GD11p;
           else FatalError( Form( "Invalid structure functions mode: %s", sf ) );
         }
       } catch ( const libconfig::SettingNotFoundException& nfe ) {
         FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
       } catch ( const libconfig::SettingTypeException& te ) {
         FatalError( Form( "Field \"%s\" has wrong type.", te.getPath() ) );
       }
     }
 
     void
     ConfigHandler::parseOutgoingKinematics( const libconfig::Setting& kin )
     {
       try {
         if ( kin.exists( "pair" ) ) {
           Particle::ParticleCode pair = (Particle::ParticleCode)(int)kin["pair"];
           params_.kinematics.central_system = { pair, pair };
         }
         if ( kin.exists( "min_pt" ) ) params_.kinematics.cuts.central[Cuts::pt_single].min() = (double)kin["min_pt"];
         if ( kin.exists( "max_pt" ) ) params_.kinematics.cuts.central[Cuts::pt_single].max() = (double)kin["max_pt"];
         if ( kin.exists( "min_ptdiff" ) ) params_.kinematics.cuts.central[Cuts::pt_diff].min() = (double)kin["min_ptdiff"];
         if ( kin.exists( "max_ptdiff" ) ) params_.kinematics.cuts.central[Cuts::pt_diff].max() = (double)kin["max_ptdiff"];
         if ( kin.exists( "min_rapiditydiff" ) ) params_.kinematics.cuts.central[Cuts::rapidity_diff].min() = (double)kin["min_rapiditydiff"];
         if ( kin.exists( "max_rapiditydiff" ) ) params_.kinematics.cuts.central[Cuts::rapidity_diff].max() = (double)kin["max_rapiditydiff"];
         if ( kin.exists( "min_energy" ) ) params_.kinematics.cuts.central[Cuts::energy_single].min() = (double)kin["min_energy"];
         if ( kin.exists( "max_energy" ) ) params_.kinematics.cuts.central[Cuts::energy_single].max() = (double)kin["max_energy"];
         if ( kin.exists( "min_eta" ) ) params_.kinematics.cuts.central[Cuts::eta_single].min() = (double)kin["min_eta"];
         if ( kin.exists( "max_eta" ) ) params_.kinematics.cuts.central[Cuts::eta_single].max() = (double)kin["max_eta"];
         if ( kin.exists( "min_rapidity" ) ) params_.kinematics.cuts.central[Cuts::rapidity_single].min() = (double)kin["min_rapidity"];
         if ( kin.exists( "max_rapidity" ) ) params_.kinematics.cuts.central[Cuts::rapidity_single].max() = (double)kin["max_rapidity"];
         if ( kin.exists( "min_mx" ) ) params_.kinematics.cuts.remnants[Cuts::mass].min() = (double)kin["min_mx"];
         if ( kin.exists( "max_mx" ) ) params_.kinematics.cuts.remnants[Cuts::mass].max() = (double)kin["max_mx"];
       } catch ( const libconfig::SettingNotFoundException& nfe ) {
         FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
       } catch ( const libconfig::SettingTypeException& te ) {
         FatalError( Form( "Field \"%s\" has wrong type.", te.getPath() ) );
       }
     }
 
     void
     ConfigHandler::parseIntegrator( const libconfig::Setting& integr )
     {
       try {
         if ( integr.exists( "algorithm" ) ) {
           const std::string algo = integr["algorithm"];
           if ( algo == "Vegas" ) params_.integrator.type = Integrator::Vegas;
           else if ( algo == "MISER" ) params_.integrator.type = Integrator::MISER;
         }
         if ( integr.exists( "num_points" ) ) params_.integrator.npoints = (int)integr["num_points"];
         if ( integr.exists( "num_integration_calls" ) ) params_.integrator.ncvg = (int)integr["num_integration_calls"];
         if ( integr.exists( "num_integration_iterations" ) ) params_.integrator.itvg = (int)integr["num_integration_iterations"];
         if ( integr.exists( "seed" ) ) params_.integrator.seed = (unsigned long)integr["seed"];
       } catch ( const libconfig::SettingNotFoundException& nfe ) {
         FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
       }
     }
 
     void
     ConfigHandler::parseGenerator( const libconfig::Setting& gen )
     {
       params_.generation.enabled = true;
       try {
         if ( gen.exists( "num_events" ) ) params_.generation.maxgen = (int)gen["num_events"];
         if ( gen.exists( "print_every" ) ) params_.generation.gen_print_every = (int)gen["print_every"];
       } catch ( const libconfig::SettingNotFoundException& nfe ) {
         FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
       }
     }
 
     void
     ConfigHandler::parseTamingFunctions( const libconfig::Setting& tf )
     {
       if ( !tf.isList() ) FatalError( "The taming functions definition must be wrapped within a list!" );
       for ( unsigned short i = 0; i < tf.getLength(); ++i ) {
         params_.taming_functions.add( tf[i]["variable"], tf[i]["expression"] );
       }
     }
 
     void
     ConfigHandler::parseHadroniser( const libconfig::Setting& hadr )
     {
       try {
         std::string hadr_name = hadr["name"];
         if ( hadr_name == "pythia8" ) {
 #ifdef PYTHIA8
           Hadroniser::Pythia8Hadroniser* pythia8 = new Hadroniser::Pythia8Hadroniser;
+          long long seed = ( hadr.exists( "seed" ) ) ? hadr["seed"] : -1ll;
+          pythia8->setSeed( seed );
           pythia8->readString( Form( "Beams:eCM = %.2f", params_.kinematics.sqrtS() ) );
-          if ( hadr.exists( "decayWto" ) ) {
-            libconfig::Setting& pdg_ids = hadr["decayWto"];
-            if ( !pdg_ids.isList() || pdg_ids.getLength() != 2 ) throw libconfig::SettingTypeException( pdg_ids );
-            std::ostringstream os;
-            pythia8->readString( "24:onMode = off" ); // disable all W decays, but...
-            for ( unsigned short i = 0; i < pdg_ids.getLength(); ++i ) {
-              int pdg_id = pdg_ids[i];
-              switch ( abs( pdg_id ) ) {
-                case 1: case 2: // u/d quark
-                  os << " 1 2"; break;
-                case 3: case 4: // c/s quark
-                  os << " 3 4"; break;
-                case 5: case 6: // t/b quark
-                  os << " 5 6"; break;
-                case 11: // e/nu_e
-                  os << " 11 12"; break;
-                case 13: // mu/nu_mu
-                  os << " 13 14"; break;
-                case 15: // tau/nu_tau
-                  os << " 15 16"; break;
-                default: { FatalError( Form( "Unsupported W decay mode: %d", pdg_id ) ); } break;
-              }
+          if ( hadr.exists( "pythiaConfiguration" ) ) {
+            libconfig::Setting& configs = hadr["pythiaConfiguration"];
+            if ( !configs.isList() ) throw libconfig::SettingTypeException( configs );
+            for ( unsigned short i = 0; i < configs.getLength(); ++i ) {
+              std::string config = configs[i];
+              pythia8->readString( config );
             }
-            pythia8->readString( Form( "24:onIfAny = %s", os.str().c_str() ) );
           }
           pythia8->init();
           params_.setHadroniser( pythia8 );
 #endif
         }
       } catch ( const libconfig::SettingTypeException& nfe ) {
         FatalError( Form( "Wrong type to the field \"%s\".", nfe.getPath() ) );
       } catch ( const libconfig::SettingNotFoundException& nfe ) {
         FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) );
       }
     }
 
     void
     ConfigHandler::writeProcess( const Parameters* params, libconfig::Setting& root )
     {
       libconfig::Setting& proc = root.add( "process", libconfig::Setting::TypeGroup );
       proc.add( "name", libconfig::Setting::TypeString ) = params->processName();
       std::ostringstream os; os << params->kinematics.mode;
       proc.add( "mode", libconfig::Setting::TypeString ) = os.str();
     }
 
     void
     ConfigHandler::writeIncomingKinematics( const Parameters* params, libconfig::Setting& root )
     {
       libconfig::Setting& kin = root.add( "in_kinematics", libconfig::Setting::TypeGroup );
       kin.add( "beam1_pz", libconfig::Setting::TypeFloat ) = params->kinematics.inp.first;
       kin.add( "beam2_pz", libconfig::Setting::TypeFloat ) = params->kinematics.inp.second;
       std::ostringstream os; os << params->kinematics.structure_functions;
       kin.add( "structure_function", libconfig::Setting::TypeString ) = os.str();
     }
 
     void
     ConfigHandler::writeOutgoingKinematics( const Parameters* params, libconfig::Setting& root )
     {
       libconfig::Setting& kin = root.add( "out_kinematics", libconfig::Setting::TypeGroup );
       if ( params->kinematics.central_system.size() > 0 )
         kin.add( "pair", libconfig::Setting::TypeInt ) = (int)params->kinematics.central_system[0];
       if ( params->kinematics.cuts.central.count( Cuts::pt_single ) ) {
         kin.add( "min_pt", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::pt_single ).min();
         kin.add( "max_pt", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::pt_single ).max();
       }
       if ( params->kinematics.cuts.central.count( Cuts::pt_diff ) ) {
         kin.add( "min_ptdiff", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::pt_diff ).min();
         kin.add( "max_ptdiff", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::pt_diff ).max();
       }
       if ( params->kinematics.cuts.central.count( Cuts::rapidity_diff ) ) {
         kin.add( "min_dely", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::rapidity_diff ).min();
         kin.add( "max_dely", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::rapidity_diff ).max();
       }
       if ( params->kinematics.cuts.central.count( Cuts::energy_single ) ) {
         kin.add( "min_energy", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::energy_single ).min();
         kin.add( "max_energy", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::energy_single ).max();
       }
       if ( params->kinematics.cuts.central.count( Cuts::eta_single ) ) {
         kin.add( "min_eta", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::eta_single ).min();
         kin.add( "max_eta", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::eta_single ).max();
       }
       if ( params->kinematics.cuts.remnants.count( Cuts::mass ) ) {
         kin.add( "min_mx", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.remnants.at( Cuts::mass ).min();
         kin.add( "max_mx", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.remnants.at( Cuts::mass ).max();
       }
     }
 
     void
     ConfigHandler::writeTamingFunctions( const Parameters* params, libconfig::Setting& root )
     {
       libconfig::Setting& tf = root.add( "taming_functions", libconfig::Setting::TypeList );
       for ( std::map<std::string,TamingFunction>::const_iterator it = params->taming_functions.begin(); it != params->taming_functions.end(); ++it ) {
         libconfig::Setting& fun = tf.add( libconfig::Setting::TypeGroup );
         fun.add( "variable", libconfig::Setting::TypeString ) = it->first;
         fun.add( "expression", libconfig::Setting::TypeString ) = it->second.expression;
       }
     }
 
     void
     ConfigHandler::writeIntegrator( const Parameters* params, libconfig::Setting& root )
     {
       libconfig::Setting& integr = root.add( "integrator", libconfig::Setting::TypeGroup );
       std::ostringstream os; os << params->integrator.type;
       integr.add( "algorithm", libconfig::Setting::TypeString ) = os.str();
       integr.add( "num_points", libconfig::Setting::TypeInt ) = (int)params->integrator.npoints;
       integr.add( "num_integration_calls", libconfig::Setting::TypeInt ) = (int)params->integrator.ncvg;
       integr.add( "num_integration_iterations", libconfig::Setting::TypeInt ) = (int)params->integrator.itvg;
       integr.add( "seed", libconfig::Setting::TypeInt64 ) = (long)params->integrator.seed;
     }
 
     void
     ConfigHandler::writeGenerator( const Parameters* params, libconfig::Setting& root )
     {
       if ( !params->generation.enabled ) return;
       libconfig::Setting& gen = root.add( "generator", libconfig::Setting::TypeGroup );
       gen.add( "num_events", libconfig::Setting::TypeInt ) = (int)params->generation.maxgen;
       gen.add( "print_every", libconfig::Setting::TypeInt ) = (int)params->generation.gen_print_every;
     }
 #endif
 
     void
     ConfigHandler::store( const Parameters* params, const char* file )
     {
 #ifdef LIBCONFIG
       libconfig::Config cfg;
       libconfig::Setting& root = cfg.getRoot();
       writeProcess( params, root );
       writeIncomingKinematics( params, root["process"] );
       writeOutgoingKinematics( params, root["process"] );
       writeTamingFunctions( params, root["process"] );
       writeIntegrator( params, root );
       writeGenerator( params, root );
       cfg.writeFile( file );
 #endif
     }
   }
 }
 
 #endif
diff --git a/CepGen/Core/Integrand.cpp b/CepGen/Core/Integrand.cpp
index 97a9f0f..b223c01 100644
--- a/CepGen/Core/Integrand.cpp
+++ b/CepGen/Core/Integrand.cpp
@@ -1,141 +1,143 @@
 #include "CepGen/Core/Timer.h"
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Event/Event.h"
 #include "CepGen/Event/Particle.h"
 #include "CepGen/Physics/Kinematics.h"
 #include "CepGen/Processes/GenericProcess.h"
 #include "CepGen/Parameters.h"
 
 #include <sstream>
 
 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();
 
+      std::cout << "event cleared!!! " << ev->getByRole( Particle::CentralSystem ).size() << std::endl;
+
       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( Particle::uQuark, 1 ); break;
           case Kinematics::InelasticElastic: // set one of the outgoing protons to be fragmented
             op1.setPdgId( Particle::uQuark, 1 ); break;
           case Kinematics::InelasticInelastic: // set both the outgoing protons to be fragmented
             op1.setPdgId( Particle::uQuark, 1 );
             op2.setPdgId( Particle::uQuark, 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] == Particle::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.;
 
     //--- only fill in the process' Event object if storage is requested
     //    or if taming functions are to be applied
 
     p->process()->fillKinematics();
 
     //--- once the kinematics variables have been populated,
     //    can apply the collection of taming functions
 
     double taming = 1.0;
     if ( p->taming_functions.has( "m_central" ) || p->taming_functions.has( "pt_central" ) ) {
       const Particle::Momentum central_system( ev->getByRole( Particle::CentralSystem )[0].momentum() + ev->getByRole( Particle::CentralSystem )[1].momentum() );
       taming *= p->taming_functions.eval( "m_central", central_system.mass() );
       taming *= p->taming_functions.eval( "pt_central", central_system.pt() );
     }
     if ( p->taming_functions.has( "q2" ) ) {
       taming *= p->taming_functions.eval( "q2", ev->getOneByRole( Particle::Parton1 ).momentum().mass() );
       taming *= p->taming_functions.eval( "q2", ev->getOneByRole( Particle::Parton2 ).momentum().mass() );
     }
     integrand *= taming;
 
     //--- full event content (+ hadronisation) if generating events
 
     if ( p->hadroniser() ) {
       p->hadroniser()->hadronise( ev.get() );
     }
 
     if ( p->storage() ) {
 
       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 f732c84..89c2da1 100644
--- a/CepGen/Event/Event.cpp
+++ b/CepGen/Event/Event.cpp
@@ -1,264 +1,294 @@
 #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()
   {}
 
   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;
     return *this;
   }
 
   void
   Event::clear()
   {
     particles_.clear();
     time_generation = -1.;
     time_total = -1.;
   }
 
   void
-  Event::init()
+  Event::freeze()
   {
-    last_particle_ = particles_.end();
+    //--- store a snapshot of the primordial event block
+    if ( particles_.count( Particle::CentralSystem ) != 0 )
+      last_cs_particle_ = particles_[Particle::CentralSystem].end();
+    if ( particles_.count( Particle::OutgoingBeam1 ) != 0 )
+      last_op1_particle_ = particles_[Particle::OutgoingBeam1].end();
+    if ( particles_.count( Particle::OutgoingBeam2 ) != 0 )
+      last_op2_particle_ = particles_[Particle::OutgoingBeam2].end();
   }
 
   void
   Event::restore()
   {
     //--- remove all particles after the primordial event block
-    particles_.erase( last_particle_, particles_.end() );
+    if ( particles_.count( Particle::CentralSystem ) != 0 )
+      particles_[Particle::CentralSystem].erase( last_cs_particle_, particles_[Particle::CentralSystem].end() );
+    if ( particles_.count( Particle::OutgoingBeam1 ) != 0 )
+      particles_[Particle::OutgoingBeam1].erase( last_op1_particle_, particles_[Particle::OutgoingBeam1].end() );
+    if ( particles_.count( Particle::OutgoingBeam2 ) != 0 )
+      particles_[Particle::OutgoingBeam2].erase( last_op2_particle_, particles_[Particle::OutgoingBeam2].end() );
   }
 
   Particles&
   Event::getByRole( const Particle::Role& role )
   {
     //--- retrieve all particles with a given role
     return particles_[role];
   }
 
+  ParticlesIds
+  Event::getIdsByRole( const 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( const Particle::Role& role )
   {
     //--- retrieve the first particle a the given role
     Particles& parts_by_role = getByRole( 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;
     for ( it=particles_.begin(), end=particles_.end(); it!=end; it=particles_.upper_bound( it->first ) ) {
       out.emplace_back( it->first );
     }
     return out;
   }
 
-  void
-  Event::addParticle( Particle part, bool replace )
+  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();
   }
 
-  void
+  Particle&
   Event::addParticle( const Particle::Role& role, bool replace )
   {
-    addParticle( Particle( role ), 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 ( part->status() == Particle::Undefined || part->status() == Particle::FinalState ) {
           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
   {
     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;
       { std::ostringstream oss; oss << part.pdgId(); os << Form( "\n %2d\t%+6d%8s", 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() ); }
       const ParticlesIds mothers = part.mothers();
       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 % 9.5e", mom.px(), mom.py(), mom.pz(), part.energy(), part.mass() );
+      os << Form( "% 9.6e % 9.6e % 9.6e % 9.6e % 12.7f", mom.px(), mom.py(), mom.pz(), part.energy(), part.mass() );
       if ( part.status() == Particle::Undefined
         || part.status() == Particle::Undecayed
+        || part.status() == Particle::Unfragmented
         || part.status() == Particle::FinalState ) {
         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------------------------------------------------------  -----------"
     "%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 ) );
   }
 }
diff --git a/CepGen/Event/Event.h b/CepGen/Event/Event.h
index 9fef87c..5ec67de 100644
--- a/CepGen/Event/Event.h
+++ b/CepGen/Event/Event.h
@@ -1,119 +1,120 @@
 #ifndef CepGen_Event_Event_h
 #define CepGen_Event_Event_h
 
 #include "Particle.h"
 #include "CepGen/Core/Logger.h"
 
 namespace CepGen
 {
   /**
    * Class containing all the information on the in- and outgoing particles' kinematics
    * \brief Kinematic information on the particles in the event
    */
   class Event {
     public:
       Event();
       ~Event();
       /**
        * \brief Copies all the relevant quantities from one Event object to another
        */
       Event& operator=( const Event& );
       /**
        * \brief Empties the whole event content
        */
       void clear();
       /// Initialize an "empty" event collection
-      void init();
+      void freeze();
       /// Restore the event to its "empty" state
       void restore();
 
       /// Dump all the known information on every Particle object contained in this Event container in the output stream
       /// \param[out] os Output stream where to dump the information
       /// \param[in] stable_ Do we only show the stable particles in this event?
       void dump( std::ostream& os = Logger::get().outputStream, bool stable_=false ) const;
 
       //----- particles adders
 
       /// Set the information on one particle in the process
       /// \param[in] part The Particle object to insert or modify in the event
       /// \param[in] replace Do we replace the particle if already present in the event or do we append another particle with the same role ?
-      void addParticle( Particle part, bool replace=false );
+      Particle& addParticle( Particle& part, bool replace=false );
       /// \brief Create a new particle in the event, with no kinematic information but the role it has to play in the process
       /// \param[in] role The role the particle will play in the process
       /// \param[in] replace Do we replace the particle if already present in the event or do we append another particle with the same role ?
-      void addParticle( const Particle::Role& role, bool replace=false );
+      Particle& addParticle( const Particle::Role& role, bool replace=false );
 
       //----- particles retrievers
 
       /// Number of particles in the event
       size_t numParticles() const;
       /// \brief Vector of all particles in the event
       const Particles particles() const;
       /// \brief Vector of all stable particles in the event
       const Particles stableParticles() const;
       /**
        * Returns the list of Particle objects corresponding to a certain role in the process kinematics
        * \brief Gets a list of particles by their role in the event
        * \param[in] role The role the particles have to play in the process
        * \return A vector of references to the requested Particle objects
        */
       Particles& getByRole( const Particle::Role& role );
+      ParticlesIds getIdsByRole( const Particle::Role& role ) const;
       /**
        * Returns the first Particle object in the particles list whose role corresponds to the given argument
        * \param[in] role The role the particle has to play in the event
        * \return A Particle object corresponding to the first particle with the role
        */
       Particle& getOneByRole( const Particle::Role& role );
       /**
        * Returns the reference to the Particle object corresponding to a unique identifier in the event
        * \brief Gets one particle by its unique identifier in the event
        * \param[in] id_ The unique identifier to this particle in the event
        * \return A reference to the requested Particle object
        */
       Particle& getById( int id_ );
       /// Get a const Particle object using its unique identifier
       /// \param[in] id_ Unique identifier of the particle in the event
       /// \return Constant object to be retrieved
       const Particle& getConstById( int id_ ) const;
       /**
        * Returns the references to the Particle objects corresponding to the unique identifiers in the event
        * \brief Gets a vector of particles by their unique identifier in the event
        * \param[in] ids_ The unique identifiers to the particles to be selected in the event
        * \return A vector of references to the requested Particle objects
        */
       Particles getByIds( const ParticlesIds& ids_ ) const;
       /**
        * Returns the list of mother particles of any given Particle object in this event
        * \param[in] part The reference to the Particle object from which we want to extract the mother particles
        * \return A list of parenting Particle object
        */
 
       //----- general particles information retriever
 
       Particles mothers( const Particle& part );
       /// Get a vector containing all the daughters from a particle
       /// \param[in] part The particle for which the daughter particles have to be retrieved
       /// \return Vector of Particle objects containing all the daughters' kinematic information
       Particles daughters( const Particle& part );
       /// Get a list of roles for the given event (really process-dependant for the central system)
       /// \return Vector of integers corresponding to all the roles the particles can play in the event
       ParticleRoles roles() const;
 
       /// Number of trials before the event was "correctly" hadronised
       int num_hadronisation_trials;
       /// Time needed to generate the event at parton level (in seconds)
       float time_generation;
       /// Time needed to generate the hadronised (if needed) event (in seconds)
       float time_total;
 
     private:
       /// Check if the event kinematics is properly defined
       void checkKinematics() const;
       /// List of particles in the event, mapped to their role in the process
       ParticlesMap particles_;
       /// Last particle in an "empty" event
-      ParticlesMap::iterator last_particle_;
+      Particles::iterator last_cs_particle_, last_op1_particle_, last_op2_particle_;
   };
 }
 
 #endif
diff --git a/CepGen/Event/Momentum.cpp b/CepGen/Event/Momentum.cpp
index 235d6a6..ce55b33 100644
--- a/CepGen/Event/Momentum.cpp
+++ b/CepGen/Event/Momentum.cpp
@@ -1,315 +1,321 @@
 #include "Particle.h"
 
 #include "CepGen/Core/Exception.h"
 #include "CepGen/Core/utils.h"
 
 #include <math.h>
 
 namespace CepGen
 {
   using Momentum = Particle::Momentum;
 
   //----- Particle momentum methods
 
   Momentum::Momentum() :
     px_( 0. ), py_( 0. ), pz_( 0. ), p_( 0. ), energy_( 0. )
   {}
 
   Momentum::Momentum( double x, double y, double z, double t ) :
     px_( x ), py_( y ), pz_( z ), energy_( t )
   {
     computeP();
   }
 
   //--- static constructors
 
   Momentum
   Momentum::fromPtEtaPhi( double pt, double eta, double phi, double e )
   {
     const double px = pt*cos( phi ),
                  py = pt*sin( phi ),
                  pz = pt*sinh( eta );
     return Momentum( px, py, pz, e );
   }
 
   Momentum
   Momentum::fromPThetaPhi( double p, double theta, double phi, double e )
   {
     const double px = p*sin( theta )*cos( phi ),
                  py = p*sin( theta )*sin( phi ),
                  pz = p*cos( theta );
     return Momentum( px, py, pz, e );
   }
 
   Momentum
   Momentum::fromPxPyPzE( double px, double py, double pz, double e )
   {
     return Momentum( px, py, pz, e );
   }
 
   //--- arithmetic operators
 
   Momentum&
   Momentum::operator+=( const Momentum& mom )
   {
     px_ += mom.px_;
     py_ += mom.py_;
     pz_ += mom.pz_;
     energy_ += mom.energy_; //FIXME not supposed to be this way!
     computeP();
     return *this;
   }
 
   Momentum&
   Momentum::operator-=( const Momentum& mom )
   {
     px_ -= mom.px_;
     py_ -= mom.py_;
     pz_ -= mom.pz_;
     energy_ -= mom.energy_; //FIXME not supposed to be this way!
     computeP();
     return *this;
   }
 
+  bool
+  Momentum::operator==( const Momentum& mom ) const
+  {
+    return ( px_ == mom.px_ && py_ == mom.py_ && pz_ == mom.pz_ && energy_ == mom.energy_ );
+  }
+
   double
   Momentum::threeProduct( const Momentum& mom ) const
   {
     DebuggingInsideLoop( Form( "  (%f, %f, %f, %f)\n\t* (%f, %f, %f, %f)\n\t= %f",
       px_, py_, pz_, energy_,
       mom.px_, mom.py_, mom.pz_, mom.energy_,
       px_*mom.px_+py_*mom.py_+pz_*mom.pz_ ) );
     return px_*mom.px_+py_*mom.py_+pz_*mom.pz_;
   }
 
   double
   Momentum::fourProduct( const Momentum& mom ) const
   {
     DebuggingInsideLoop( Form( "  (%f, %f, %f, %f)\n\t* (%f, %f, %f, %f)\n\t= %f",
       px_, py_, pz_, energy_,
       mom.px_, mom.py_, mom.pz_, mom.energy_,
       px_*mom.px_+py_*mom.py_+pz_*mom.pz_ ) );
     return energy_*mom.energy_-threeProduct(mom);
   }
 
   double
   Momentum::operator*=( const Momentum& mom )
   {
     return threeProduct( mom );
   }
 
   Momentum&
   Momentum::operator*=( double c )
   {
     px_ *= c;
     py_ *= c;
     pz_ *= c;
     computeP();
     return *this;
   }
 
   Momentum
   operator*( const Momentum& mom, double c )
   {
     Momentum out = mom;
     out *= c;
     return out;
   }
 
   Momentum
   operator*( double c, const Momentum& mom )
   {
     Momentum out = mom;
     out *= c;
     return out;
   }
 
   Momentum
   operator+( const Momentum& mom1, const Momentum& mom2 )
   {
     Momentum out = mom1;
     out += mom2;
     return out;
   }
 
   Momentum
   operator-( const Momentum& mom1, const Momentum& mom2 )
   {
     Momentum out = mom1;
     out -= mom2;
     return out;
   }
 
   double
   operator*( const Momentum& mom1, const Momentum& mom2 )
   {
     Momentum tmp = mom1;
     return tmp.threeProduct( mom2 );
   }
 
   //--- various setters
 
   void
   Momentum::setMass2( double m2 )
   {
     energy_ = sqrt( p2()+m2 );
   }
 
   void
   Momentum::setP( double px, double py, double pz, double e )
   {
     setP( px, py, pz );
     setEnergy( e );
   }
 
   void
   Momentum::setP( double px, double py, double pz )
   {
     px_ = px;
     py_ = py;
     pz_ = pz;
     computeP();
   }
 
   void
   Momentum::setP( unsigned int i, double p )
   {
     switch ( i ) {
       case 0: px_ = p; break;
       case 1: py_ = p; break;
       case 2: pz_ = p; break;
       case 3: energy_ = p; break;
       default: return;
     }
     computeP();
   }
 
   void
   Momentum::computeP()
   {
     p_ = sqrt( px_*px_ + py_*py_ + pz_*pz_ );
   }
 
   //--- various getters
 
   double
   Momentum::operator[]( const unsigned int i ) const {
     switch ( i ) {
       case 0: return px_;
       case 1: return py_;
       case 2: return pz_;
       case 3: return energy_;
       default: return -1.;
     }
   }
 
   const std::vector<double>
   Momentum::pVector() const
   {
     return std::vector<double>( { px(), py(), pz(), energy(), mass() } );
   }
 
   double
   Momentum::mass() const
   {
     if ( mass2() >= 0. ) return sqrt( mass2() );
     return -sqrt( -mass2() );
   }
 
   double
   Momentum::theta() const
   {
     return atan2( pt(), pz() );
   }
 
   double
   Momentum::phi() const
   {
     return atan2( py(), px() );
   }
 
   double
   Momentum::pt() const
   {
     return sqrt( pt2() );
   }
 
   double
   Momentum::eta() const
   {
     const int sign = ( pz()/fabs( pz() ) );
     return ( pt() != 0. )
       ? log( ( p()+fabs( pz() ) )/pt() )*sign
       : 9999.*sign;
   }
 
   double
   Momentum::rapidity() const
   {
     const int sign = ( pz()/fabs( pz() ) );
     return ( energy() >= 0. )
       ? log( ( energy()+pz() )/( energy()-pz() ) )*0.5
       : 999.*sign;
   }
 
   //--- boosts/rotations
 
   Momentum&
   Momentum::betaGammaBoost( double gamma, double betagamma )
   {
     if ( gamma == 1. && betagamma == 0. ) return *this; // trivial case
     const double pz = pz_, e = energy_;
     pz_ = gamma*pz+betagamma*e;
     energy_  = gamma*e +betagamma*pz;
     computeP();
     return *this;
   }
 
   Momentum&
   Momentum::lorentzBoost( const Momentum& p )
   {
     const double m = p.mass();
     if ( m == p.energy() ) return *this;
 
     const double pf4 = ( ( *this )*p ) / m,
                  fn = ( pf4+energy() )/( p.energy()+m );
     *this -= p*fn;
     setEnergy( pf4 );
     return *this;
   }
 
   Momentum&
   Momentum::rotatePhi( double phi, double sign )
   {
     const double px = px_*cos( phi )+py_*sin( phi )*sign,
                  py =-px_*sin( phi )+py_*cos( phi )*sign;
     px_ = px;
     py_ = py;
     return *this;
   }
 
   Momentum&
   Momentum::rotateThetaPhi( double theta, double phi )
   {
     double rotmtx[3][3], mom[3]; //FIXME check this! cos(phi)->-sin(phi) & sin(phi)->cos(phi) --> phi->phi+pi/2 ?
     rotmtx[0][0] = -sin( phi ); rotmtx[0][1] = -cos( theta )*cos( phi ); rotmtx[0][2] =  sin( theta )*cos( phi );
     rotmtx[1][0] =  cos( phi ); rotmtx[1][1] = -cos( theta )*sin( phi ); rotmtx[1][2] =  sin( theta )*sin( phi );
     rotmtx[2][0] =  0.;         rotmtx[2][1] =  sin( theta );            rotmtx[2][2] =  cos( theta );
 
     for (int i=0; i<3; i++) {
       mom[i] = 0.;
       for (int j=0; j<3; j++) {
         mom[i] += rotmtx[i][j]*operator[]( j );
       }
     }
     setP( mom[0], mom[1], mom[2] );
     return *this;
   }
 
   //--- printout
 
   std::ostream&
   operator<<( std::ostream& os, const Momentum& mom )
   {
     return os << "(E ; p) = (" << mom.energy_ << " ; " << mom.px_ << ", " << mom.py_ << ", " << mom.pz_ << ")";
   }
 }
diff --git a/CepGen/Event/Particle.h b/CepGen/Event/Particle.h
index b561a3e..724948a 100644
--- a/CepGen/Event/Particle.h
+++ b/CepGen/Event/Particle.h
@@ -1,402 +1,404 @@
 #ifndef CepGen_Event_Particle_h
 #define CepGen_Event_Particle_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:
       /** 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
         dQuark = 1, uQuark = 2,
         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,
         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,
         ud0Diquark = 2101, ud1Diquark = 2103, uu1Diquark = 2203,
         Proton = 2212, Neutron = 2112,
         Pomeron = 990, Reggeon = 110
       };
       /// Internal status code for a particle
       enum Status {
         PrimordialIncoming = -9,
-        Undecayed = -3,
+        Unfragmented = -4, Undecayed = -3,
         Propagator = -2,
         Incoming = -1,
         Undefined = 0,
         FinalState = 1,
         Resonance = 2,
         DebugResonance = 3,
         PythiaHIncoming = 21
       };
       /// 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. );
           inline ~Momentum() {}
 
           // --- 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;
           /// 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 Particle::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
 
       /**
        * \brief Gets the mass of a particle
        * \param pdgId ParticleCode (PDG ID)
        * \return Mass of the particle in \f$\textrm{GeV}/c^2\f$
        */
       static double massFromPDGId( const Particle::ParticleCode& pdgId );
       /**
        * \brief Gets the electric charge of a particle
        * \param pdgId ParticleCode (PDG ID)
        * \return Charge of the particle in \f$e\f$
        */
       static double chargeFromPDGId( const Particle::ParticleCode& pdgId );
       /**
        * \brief Total decay width of one unstable particle
        * \param[in] pdgId ParticleCode (PDG ID)
        * \return Decay width in GeV
        */
       static double widthFromPDGId( const Particle::ParticleCode& pdgId );
 
       /// 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_ * chargeFromPDGId( 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( const Status& status ) { status_ = 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 1d1414f..0d06ff5 100644
--- a/CepGen/Hadronisers/Pythia8Hadroniser.cpp
+++ b/CepGen/Hadronisers/Pythia8Hadroniser.cpp
@@ -1,55 +1,116 @@
 #ifdef PYTHIA8
 #include "Pythia8Hadroniser.h"
 
 namespace CepGen
 {
   namespace Hadroniser
   {
     Pythia8Hadroniser::Pythia8Hadroniser():
       GenericHadroniser( "pythia8" ),
       pythia_( std::unique_ptr<Pythia8::Pythia>( new Pythia8::Pythia ) )
     {}
 
     Pythia8Hadroniser::~Pythia8Hadroniser()
     {}
 
+    void
+    Pythia8Hadroniser::setSeed( long long seed )
+    {
+      if ( seed == -1ll ) {
+        pythia_->readString( "Random:setSeed = off" );
+        return;
+      }
+      pythia_->readString( "Random:setSeed = on" );
+      pythia_->readString( Form( "Random:seed = %llu", seed ) );
+    }
+
     bool
     Pythia8Hadroniser::hadronise( const Particle* part )
     {
       return true;
     }
 
     bool
     Pythia8Hadroniser::hadronise( Event* ev )
     {
+      //--- first start by cleaning up the previous runs leftovers
+
+      pythia_->event.reset();
+      std::cout << ">>>>> before hadronisation" << std::endl;
+      ev->dump();
+
       //--- outgoing proton remnants fragmentation
 
       Particles& op1 = ev->getByRole( Particle::OutgoingBeam1 ), &op2 = ev->getByRole( Particle::OutgoingBeam2 );
       for ( Particles::const_iterator p_it = op1.begin(); p_it != op1.end(); ++p_it ) {
         //... excited proton fragmentation
       }
       for ( Particles::const_iterator p_it = op2.begin(); p_it != op2.end(); ++p_it ) {
         //... excited proton fragmentation
       }
 
       //--- central system hadronisation/decay/...
 
-      Particles& cs = ev->getByRole( Particle::CentralSystem );
-      for ( Particles::const_iterator p_it = cs.begin(); p_it != cs.end(); ++p_it ) {
-        if ( p_it->status() != Particle::Undecayed ) continue;
-        p_it->dump();
+      const ParticlesIds cs = ev->getIdsByRole( Particle::CentralSystem );
+      for ( ParticlesIds::iterator p_it = cs.begin(); p_it != cs.end(); ++p_it ) {
+        const Particle& part = ev->getById( *p_it );
+
+        // check if the particle can be decayed
+        if ( part.status() != Particle::Undecayed ) continue;
+        if ( !pythia_->particleData.canDecay( part.pdgId() ) ) continue;
+
+        // check if the particle has a valid kinematics
+        const Particle::Momentum mom = part.momentum();
+        if ( mom == Particle::Momentum() ) continue; // FIXME probably during warm-up...
+
+        pythia_->event.reset();
+        Pythia8::Particle py8part( part.integerPdgId(), 93, 0, 0, 0, 0, 0, 0, mom.px(), mom.py(), mom.pz(), mom.energy(), part.mass() );
+        py8part.tau( pythia_->particleData.tau0( part.pdgId() ) );
+        pythia_->event.append( py8part );
+
+        const unsigned short num_before = pythia_->event.size(); // number of particles before the decay
+        pythia_->next(); // launch the decay
+
+        // map { "pythia id" -> "cepgen id" }
+        std::map<unsigned short, unsigned short> py_cg_corresp = { { 1, part.id() } };
+
+        const unsigned short num_after = pythia_->event.size(); // number of particles after the decay
+        if ( num_before == num_after ) continue; // nothing happened... so no decay!
+
+        for ( unsigned short i = 2; i < pythia_->event.size(); ++i ) { // skip the initial system
+          const Pythia8::Particle p = pythia_->event[i];
+          const std::vector<int> mothers = p.motherList();
+          if ( mothers.size() == 0 ) continue;
+          Particle& op = ev->addParticle( Particle::CentralSystem );
+          py_cg_corresp[i] = op.id();
+
+//std::cout << "dumping ids correspondence (" << *p_it << ")" << std::endl; for ( const auto& p : py_cg_corresp ) { std::cout << p.first << "--->" << p.second << std::endl; } std::cout << "-------enddump-------" << std::endl;
+          op.setPdgId( ( Particle::ParticleCode )abs( p.id() ), p.charge() );
+          if ( p.isFinal() ) op.setStatus( Particle::FinalState );
+          else op.setStatus( 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 ( 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] ) );
+          }
+        }
+        ev->getById( *p_it ).setStatus( Particle::Resonance );
+        //std::cout << part.id() << "->" << part.integerPdgId() << "->" << part.status() << std::endl;
       }
 
+      ev->dump();
+
       return true;
     }
 
     bool
     Pythia8Hadroniser::prepareHadronisation( Event *ev )
     {
-      pythia_->init();
+      init();
       return true;
     }
   }
 }
 
 #endif
diff --git a/CepGen/Hadronisers/Pythia8Hadroniser.h b/CepGen/Hadronisers/Pythia8Hadroniser.h
index ac6a22f..6f76797 100644
--- a/CepGen/Hadronisers/Pythia8Hadroniser.h
+++ b/CepGen/Hadronisers/Pythia8Hadroniser.h
@@ -1,36 +1,37 @@
 #ifndef Pythia8Hadroniser_h
 #define Pythia8Hadroniser_h
 
 #include <Pythia8/Pythia.h>
 #include <memory>
 
 #include "GenericHadroniser.h"
 
 namespace CepGen
 {
   namespace Hadroniser
   {
     /**
      * Full interface to the Pythia8 hadronisation algorithm. It can be used in a single particle decay mode as well as a full event hadronisation using the string model, as in Jetset.
      * \brief Pythia8 hadronisation algorithm
      */
     class Pythia8Hadroniser : public GenericHadroniser
     {
       public:
         Pythia8Hadroniser();
         ~Pythia8Hadroniser();
         bool hadronise( const Particle* part ) override;
         bool hadronise( Event* ev ) override;
 
+        void setSeed( long long seed );
         void readString( const char* param ) { pythia_->readString( param ); }
         void readString( const std::string& param ) { pythia_->readString( param ); }
         void init() { pythia_->init(); }
 
       private:
         bool prepareHadronisation( Event *ev );
         std::unique_ptr<Pythia8::Pythia> pythia_;
     };
   }
 }
 
 #endif
diff --git a/CepGen/Processes/GamGamLL.cpp b/CepGen/Processes/GamGamLL.cpp
index dab5dfa..4d1fc92 100644
--- a/CepGen/Processes/GamGamLL.cpp
+++ b/CepGen/Processes/GamGamLL.cpp
@@ -1,877 +1,877 @@
 #include "GamGamLL.h"
 #include "CepGen/Core/Exception.h"
 
 using namespace CepGen::Process;
 
 GamGamLL::GamGamLL( int nopt ) :
   GenericProcess( "lpair", "pp -> p(*) (gamma gamma -> l+ l-) p(*)" ),
   n_opt_( nopt ),
   MX2_( 0. ), MY2_( 0. ), Ml2_( 0. ),
   ep1_( 0. ), ep2_( 0. ), p_cm_( 0. ),
   w12_( 0. ), w31_( 0. ), dw31_( 0. ), w52_( 0. ), dw52_( 0. ),
   ec4_( 0. ), pc4_( 0. ), mc4_( 0. ), w4_( 0. ),
   p12_( 0. ), p1k2_( 0. ), p2k1_( 0. ),
   p13_( 0. ), p14_( 0. ), p25_( 0. ),
   q1dq_( 0. ), q1dq2_( 0. ),
   s1_( 0. ), s2_( 0. ),
   epsi_( 0. ),
   g5_( 0. ), g6_( 0. ), a5_( 0. ), a6_( 0. ), bb_( 0. ),
   gram_( 0. ),
   dd1_( 0. ), dd2_( 0. ), dd3_( 0. ), dd4_( 0. ), dd5_( 0. ),
   delta_( 0. ),
   g4_( 0. ), sa1_( 0. ), sa2_( 0. ),
   sl1_( 0. ),
   cos_theta4_( 0. ), sin_theta4_( 0. ),
   al4_( 0. ), be4_( 0. ), de3_( 0. ), de5_( 0. ),
   pt4_( 0. ),
   jacobian_( 0. )
 {}
 
 void
 GamGamLL::addEventContent()
 {
   GenericProcess::setEventContent( {
     { Particle::IncomingBeam1, Particle::Proton },
     { Particle::IncomingBeam2, Particle::Proton },
     { Particle::Parton1, Particle::Photon },
     { Particle::Parton2, Particle::Photon }
   }, {
     { Particle::OutgoingBeam1, { Particle::Proton } },
     { Particle::OutgoingBeam2, { Particle::Proton } },
     { Particle::CentralSystem, { Particle::Muon, Particle::Muon } }
   } );
 }
 
 unsigned int
 GamGamLL::numDimensions( const Kinematics::ProcessMode& process_mode ) const
 {
   switch ( process_mode ) {
     case Kinematics::ElectronProton:     { InError( "Not supported yet!" ); }
     case Kinematics::ElasticElastic:
     default:                             return 7;
     case Kinematics::ElasticInelastic:
     case Kinematics::InelasticElastic:   return 8;
     case Kinematics::InelasticInelastic: return 9;
   }
 }
 
 bool
 GamGamLL::pickin()
 {
   DebuggingInsideLoop( Form( "Optimised mode? %i", n_opt_ ) );
 
   jacobian_ = 0.;
 
   w4_ = mc4_*mc4_;
 
   // sig1 = sigma and sig2 = sigma' in [1]
   const double sig = mc4_+MY_;
   double sig1 = sig*sig,
          sig2 = sig1;
 
   DebuggingInsideLoop( Form( "mc4 = %f\n\t"
                              "sig1 = %f\n\t"
                              "sig2 = %f", mc4_, sig1, sig2 ) );
 
   // Mass difference between the first outgoing particle and the first incoming
   // particle
   w31_ = MX2_-w1_;
   // Mass difference between the second outgoing particle and the second
   // incoming particle
   w52_ = MY2_-w2_;
   // Mass difference between the two incoming particles
   w12_ = w1_-w2_;
   // Mass difference between the central two-photons system and the second
   // outgoing particle
   const double d6 = w4_-MY2_;
 
   DebuggingInsideLoop( Form( "w1 = %f\n\t"
                              "w2 = %f\n\t"
                              "w3 = %f\n\t"
                              "w4 = %f\n\t"
                              "w5 = %f",
                              w1_, w2_, MX2_, w4_, MY2_ ) );
 
   DebuggingInsideLoop( Form( "w31 = %f\n\tw52 = %f\n\tw12 = %f", w31_, w52_, w12_ ) );
 
   const double ss = s_+w12_;
 
   const double rl1 = ss*ss-4.*w1_*s_; // lambda(s, m1**2, m2**2)
   if ( rl1 <= 0. ) { InWarning( Form( "rl1 = %f <= 0", rl1 ) ); return false; }
   sl1_ = sqrt( rl1 );
 
   s2_ = 0.;
   double ds2 = 0.;
   if ( n_opt_ == 0 ) {
     const double smax = s_+MX2_-2.*MX_*sqs_;
     Map( x(2), sig1, smax, s2_, ds2, "s2" );
     sig1 = s2_; //FIXME!!!!!!!!!!!!!!!!!!!!
   }
 
   DebuggingInsideLoop( Form( "s2 = %f", s2_ ) );
 
   //std::cout << "s=" << _s << ", w3=" << MX2_ << ", sig1=" << sig1 << std::endl;
   const double sp = s_+MX2_-sig1,
                d3 = sig1-w2_;
 
   const double rl2 = sp*sp-4.*s_*MX2_; // lambda(s, m3**2, sigma)
   if ( rl2 <= 0. ) { InWarning( Form( "rl2 = %f <= 0", rl2 ) ); return false; }
   const double sl2 = sqrt( rl2 );
 
   //std::cout << "ss=" << ss << ", sp=" << sp << ", sl1=" << sl1_ << ", sl2=" << sl2 << std::endl;
   double t1_max = w1_+MX2_-( ss*sp+sl1_*sl2 )/( 2.*s_ ), // definition from eq. (A.4) in [1]
          t1_min = ( w31_*d3+( d3-w31_ )*( d3*w1_-w31_*w2_ )/s_ )/t1_max; // definition from eq. (A.5) in [1]
 
   // FIXME dropped in CDF version
   const Kinematics::Limits q2_limits = cuts_.cuts.initial[Cuts::q2];
   if ( t1_max > -q2_limits.min() ) { InWarning( Form( "t1max = %f > -q2min = %f", t1_max, -q2_limits.min() ) ); return false; }
   if ( t1_min < -q2_limits.max() && q2_limits.hasMax() ) { Debugging( Form( "t1min = %f < -q2max = %f", t1_min, -q2_limits.max() ) ); return false; }
   if ( t1_max < -q2_limits.max() && q2_limits.hasMax() ) t1_max = -q2_limits.max();
   if ( t1_min > -q2_limits.min() )                       t1_min = -q2_limits.min();
   /////
 
   // t1, the first photon propagator, is defined here
   t1_ = 0.;
   double dt1 = 0.;
   Map( x(0), t1_min, t1_max, t1_, dt1, "t1" );
   // changes wrt mapt1 : dx->-dx
   dt1 = -dt1;
 
   DebuggingInsideLoop( Form( "Definition of t1 = %f according to\n\t"
                              "(t1min, t1max) = (%f, %f)", t1_, t1_min, t1_max ) );
 
   dd4_ = w4_-t1_;
 
   const double d8 = t1_-w2_,
 	       t13 = t1_-w1_-MX2_;
 
   sa1_ = -pow( t1_-w31_, 2 )/4.+w1_*t1_;
   if ( sa1_ >= 0. ) { InWarning( Form( "sa1_ = %f >= 0", sa1_ ) ); return false; }
 
   const double sl3 = sqrt( -sa1_ );
 
   // one computes splus and (s2x=s2max)
   double splus, s2max;
   if ( w1_ != 0. ) {
     const double sb =( s_*( t1_-w31_ )+w12_*t13 )/( 2.*w1_ )+MX2_,
                  sd = sl1_*sl3/w1_,
                  se =( s_*( t1_*( s_+t13-w2_ )-w2_*w31_ )+MX2_*( w12_*d8+w2_*MX2_ ) )/w1_;
 
     if ( fabs( ( sb-sd )/sd )>=1. ) { splus = sb-sd; s2max = se/splus; }
     else                            { s2max = sb+sd; splus = se/s2max; }
   }
   else { // 3
     s2max = ( s_*(t1_*(s_+d8-MX2_)-w2_*MX2_ )+w2_*MX2_*( w2_+MX2_-t1_ ) )/( ss*t13 );
     splus = sig2;
   }
   // 4
   double s2x = s2max;
 
   DebuggingInsideLoop( Form( "s2x = s2max = %f", s2x ) );
 
   if ( n_opt_ < 0 ) { // 5
     if ( splus > sig2 ) {
       sig2 = splus;
       DebuggingInsideLoop( Form( "sig2 truncated to splus = %f", splus ) );
     }
     if ( n_opt_ < -1 ) { Map( x(2), sig2, s2max, s2_, ds2, "s2" ); }
     else               { Mapla( t1_, w2_, x(2), sig2, s2max, s2_, ds2 ); } // n_opt_==-1
     s2x = s2_;
   }
   else if ( n_opt_ == 0 ) { s2x = s2_; } // 6
 
   DebuggingInsideLoop( Form( "s2x = %f", s2x ) );
 
   // 7
   const double r1 = s2x-d8,
                r2 = s2x-d6;
 
   const double rl4 = ( r1*r1-4.*w2_*s2x )*( r2*r2-4.*MY2_*s2x );
   if ( rl4 <= 0. ) {
     DebuggingInsideLoop( Form( "rl4 = %f <= 0", rl4 ) );
     return false;
   }
   const double sl4 = sqrt( rl4 );
 
   // t2max, t2min definitions from eq. (A.12) and (A.13) in [1]
   const double t2_max = w2_+MY2_-( r1*r2+sl4 )/s2x * 0.5,
                t2_min = ( w52_*dd4_+( dd4_-w52_ )*( dd4_*w2_-w52_*t1_ )/s2x )/t2_max;
 
   // t2, the second photon propagator, is defined here
   t2_ = 0.;
   double dt2 = 0.;
   Map( x(1), t2_min, t2_max, t2_, dt2, "t2" );
   // changes wrt mapt2 : dx->-dx
 
   dt2 = -dt2;
 
   // \f$\delta_6=m_4^2-m_5^2\f$ as defined in Vermaseren's paper
   const double tau = t1_-t2_,
                r3 = dd4_-t2_,
                r4 = w52_-t2_;
 
   DebuggingInsideLoop( Form( "r1 = %f\n\tr2 = %f\n\tr3 = %f\n\tr4 = %f", r1, r2, r3, r4 ) );
 
   const double b = r3*r4-2.*( t1_+w2_ )*t2_,
                c = t2_*d6*d8+( d6-d8 )*( d6*w2_-d8*MY2_ );
 
   const double t25 = t2_-w2_-MY2_;
 
   sa2_ = -r4*r4/4.+w2_*t2_;
   if ( sa2_ >= 0. ) { InWarning( Form( "sa2_ = %f >= 0", sa2_ ) ); return false; }
 
   const double sl6 = 2.*sqrt( -sa2_ );
 
   g4_ = -r3*r3/4.+t1_*t2_;
   if ( g4_ >= 0. )  { InWarning( Form( "g4_ = %f >= 0", g4_ ) );   return false; }
 
   const double sl7 = 2.*sqrt( -g4_ ),
                sl5 = sl6*sl7;
 
   double s2p, s2min;
   if ( fabs( ( sl5-b )/sl5 ) >= 1. ) {
     s2p = ( sl5-b )/t2_ * 0.5;
     s2min = c/( t2_*s2p );
   }
   else { // 8
     s2min = ( -sl5-b )/t2_ * 0.5;
     s2p = c/( t2_*s2min );
   }
   // 9
   if ( n_opt_ > 1 )       Map( x( 2 ), s2min, s2max, s2_, ds2, "s2" );
   else if ( n_opt_ == 1 ) Mapla( t1_, w2_, x( 2 ), s2min, s2max, s2_, ds2 );
 
   const double ap = -0.25*pow( s2_+d8, 2 )+s2_*t1_;
 
   DebuggingInsideLoop( Form( "s2 = %f, s2max = %f, splus = %f", s2_, s2max, splus ) );
 
   if ( w1_!=0. ) dd1_ = -0.25 * ( s2_-s2max ) * ( s2_-splus ) * w1_; // 10
   else           dd1_ =  0.25 * ( s2_-s2max ) * ss * t13;
   // 11
   dd2_ = -t2_*( s2_-s2p )*( s2_-s2min ) * 0.25;
 
   DebuggingInsideLoop( Form( "t2 =%f\n\ts2 =%f\n\ts2p=%f\n\ts2min=%f\n\tdd2=%f", t2_, s2_, s2p, s2min, dd2_ ) );
 
   const double yy4 = cos( M_PI*x( 3 ) );
   const double dd = dd1_*dd2_;
   p12_ = ( s_-w1_-w2_ )*0.5;
   const double st = s2_-t1_-w2_;
   const double delb = ( 2.*w2_*r3+r4*st )*( 4.*p12_*t1_-( t1_-w31_ )*st )/( 16.*ap );
 
   if ( dd <= 0. ) {
     DebuggingInsideLoop( Form( "dd = %e <= 0\n\tdd1 = %e\tdd2 = %e", dd, dd1_, dd2_ ) );
     return false;
   }
 
   delta_ = delb - yy4*st*sqrt( dd )/ ap * 0.5;
   s1_ = t2_+w1_+( 2.*p12_*r3-4.*delta_ )/st;
 
   if ( ap >= 0. ) {
     DebuggingInsideLoop( Form( "ap = %f >= 0", ap ) );
     return false;
   }
 
   jacobian_ = ds2
             * dt1
             * dt2
             * M_PI*M_PI/( 8.*sl1_*sqrt( -ap ) );
 
   DebuggingInsideLoop( Form( "Jacobian = %e", jacobian_ ) );
 
   gram_ = ( 1.-yy4*yy4 )*dd/ap;
 
   p13_ = -t13 * 0.5;
   p14_ = (tau+s1_-MX2_) * 0.5;
   p25_ = -t25 * 0.5;
 
   /*const double //p15 = (s_+t2_-s1_-w2_)/2.,
                //p23 = (s_+t1_-s2_-w1_)/2.,
                //p24 = (s2_-tau-MY2_)/2.,
                //p34 = (s1_-MX2_-w4_)/2.,
                //p35 = (s_+w4_-s1_-s2_)/2.,
                //p45 = (s2_-w4_-MY2_)/2.;*/
 
   p1k2_ = ( s1_-t2_-w1_ ) * 0.5;
   p2k1_ = st * 0.5;
 
   double s1p, s1m;
   if ( w2_ != 0. ) {
     const double sbb = ( s_*( t2_-w52_ )-w12_*t25 )/w2_ * 0.5 + MY2_,
                  sdd = sl1_*sl6/w2_ * 0.5,
                  see = ( s_*( t2_*( s_+t25-w1_ )-w1_*w52_ )+MY2_*( w1_*MY2_-w12_*( t2_-w1_ ) ) )/w2_;
     if ( sbb/sdd >= 0. ) { s1p = sbb+sdd; s1m = see/s1p; }
     else                 { s1m = sbb-sdd; s1p = see/s1m; } // 12
     dd3_ = -w2_*( s1p-s1_ )*( s1m-s1_ ) * 0.25; // 13
   }
   else { // 14
     s1p = ( s_*( t2_*( s_-MY2_+t2_-w1_ )-w1_*MY2_ )+w1_*MY2_*( w1_+MY2_-t2_ ) )/( t25*( s_-w12_ ) );
     dd3_ = -t25*( s_-w12_ )*( s1p-s1_ ) * 0.25;
   }
   // 15
   //const double acc3 = (s1p-s1_)/(s1p+s1_);
 
   const double ssb = t2_+w1_-r3*( w31_-t1_ )/t1_ * 0.5,
                ssd = sl3*sl7/t1_,
                sse = ( t2_-w1_ )*( w4_-MX2_ )+( t2_-w4_+w31_ )*( ( t2_-w1_)*MX2_-(w4_-MX2_)*w1_)/t1_;
 
   double s1pp, s1pm;
   if (ssb/ssd>=0.) { s1pp = ssb+ssd; s1pm = sse/s1pp; }
   else             { s1pm = ssb-ssd; s1pp = sse/s1pm; } // 16
   // 17
   dd4_ = -t1_*( s1_-s1pp )*( s1_-s1pm ) * 0.25;
   //const double acc4 = ( s1_-s1pm )/( s1_+s1pm );
   dd5_ = dd1_+dd3_+( ( p12_*( t1_-w31_ )*0.5-w1_*p2k1_ )*( p2k1_*( t2_-w52_ )-w2_*r3 )-delta_*( 2.*p12_*p2k1_-w2_*( t1_-w31_ ) ) ) / p2k1_;
 
   return true;
 }
 
 bool
 GamGamLL::orient()
 {
   if ( !pickin() or jacobian_ == 0. ) { DebuggingInsideLoop( Form( "Pickin failed! Jacobian = %f", jacobian_ ) ); return false; }
 
   const double re = 0.5 / sqs_;
   ep1_ = re*( s_+w12_ );
   ep2_ = re*( s_-w12_ );
 
   DebuggingInsideLoop( Form( " re = %e\n\tw12_ = %e", re, w12_ ) );
   DebuggingInsideLoop( Form( "Incoming particles' energy = %f, %f", ep1_, ep2_ ) );
 
   p_cm_ = re*sl1_;
 
   de3_ = re*(s2_-MX2_+w12_);
   de5_ = re*(s1_-MY2_-w12_);
 
   // Final state energies
   const double ep3 = ep1_-de3_,
                ep5 = ep2_-de5_;
   ec4_ = de3_+de5_;
 
   if ( ec4_ < mc4_ ) { InWarning( Form( "ec4_ = %f < mc4_ = %f\n\t==> de3 = %f, de5 = %f", ec4_, mc4_, de3_, de5_ ) ); return false; }
 
   // What if the protons' momenta are not along the z-axis?
   pc4_ = sqrt( ec4_*ec4_-mc4_*mc4_ );
 
   if ( pc4_==0. ) {
     InWarning( "pzc4==0" );
     return false;
   }
 
   const double pp3 = sqrt( ep3*ep3-MX2_ ), pt3 = sqrt( dd1_/s_ )/p_cm_,
                pp5 = sqrt( ep5*ep5-MY2_ ), pt5 = sqrt( dd3_/s_ )/p_cm_;
 
   DebuggingInsideLoop( Form( "Central system's energy: E4 = %f\n\t"
                              "                 momentum: p4 = %f\n\t"
                              "                 invariant mass: m4 = %f\n\t"
                              "Outgoing particles' energy: E3 = %f\n\t"
                              "                            E5 = %f",
                              ec4_, pc4_, mc4_, ep3, ep5 ) );
 
   const double sin_theta3 = pt3/pp3,
                sin_theta5 = pt5/pp5;
 
   DebuggingInsideLoop( Form( "sin_theta3 = %e\n\tsin_theta5 = %e", sin_theta3, sin_theta5 ) );
 
   if (sin_theta3>1.) { InWarning( Form( "sin_theta3 = %e > 1", sin_theta3 ) ); return false; }
   if (sin_theta5>1.) { InWarning( Form( "sin_theta5 = %e > 1", sin_theta5 ) ); return false; }
 
   double ct3 = sqrt( 1.-sin_theta3*sin_theta3 ),
          ct5 = sqrt( 1.-sin_theta5*sin_theta5 );
 
   if ( ep1_*ep3 < p13_ ) ct3 *= -1.;
   if ( ep2_*ep5 > p25_ ) ct5 *= -1.;
 
   DebuggingInsideLoop( Form( "ct3 = %e\n\tct5 = %e", ct3, ct5 ) );
 
   if ( dd5_ < 0. ) { InWarning( Form( "dd5 = %f < 0", dd5_ ) ); return false; }
 
   // Centre of mass system kinematics (theta4 and phi4)
   pt4_ = sqrt( dd5_/s_ )/p_cm_;
   sin_theta4_ = pt4_/pc4_;
 
   if ( sin_theta4_ > 1. ) { InWarning( Form( "st4 = %f > 1", sin_theta4_ ) ); return false; }
 
   cos_theta4_ = sqrt( 1.-sin_theta4_*sin_theta4_ );
   if ( ep1_*ec4_ < p14_ ) cos_theta4_ *= -1.;
 
   al4_ = 1.-cos_theta4_;
   be4_ = 1.+cos_theta4_;
 
   if (cos_theta4_<0.) be4_ = sin_theta4_*sin_theta4_/al4_;
   else                al4_ = sin_theta4_*sin_theta4_/be4_;
 
   DebuggingInsideLoop( Form( "ct4 = %f\n\tal4 = %f, be4 = %f", cos_theta4_, al4_, be4_ ) );
 
   const double rr  = sqrt( -gram_/s_ )/( p_cm_*pt4_ );
   const double sin_phi3 =  rr / pt3,
                sin_phi5 = -rr / pt5;
 
   if ( fabs( sin_phi3 ) > 1. ) { InWarning( Form( "sin(phi_3) = %e while it must be in [-1 ; 1]", sin_phi3 ) ); return false; }
   if ( fabs( sin_phi5 ) > 1. ) { InWarning( Form( "sin(phi_5) = %e while it must be in [-1 ; 1]", sin_phi5 ) ); return false; }
 
   const double cos_phi3 = -sqrt( 1.-sin_phi3*sin_phi3 ),
                cos_phi5 = -sqrt( 1.-sin_phi5*sin_phi5 );
 
   p3_lab_ = Particle::Momentum( pp3*sin_theta3*cos_phi3, pp3*sin_theta3*sin_phi3, pp3*ct3, ep3 );
   p5_lab_ = Particle::Momentum( pp5*sin_theta5*cos_phi5, pp5*sin_theta5*sin_phi5, pp5*ct5, ep5 );
 
   const double a1 = p3_lab_.px()-p5_lab_.px();
 
   DebuggingInsideLoop( Form( "Kinematic quantities\n\t"
                              "cos(theta3) = %1.4f\tsin(theta3) = %1.4f\tcos( phi3 ) = %1.4f\tsin( phi3 ) = %1.4f\n\t"
                              "cos(theta4) = %1.4f\tsin(theta4) = %1.4f\n\t"
                              "cos(theta5) = %1.4f\tsin(theta5) = %1.4f\tcos( phi5 ) = %1.4f\tsin( phi5 ) = %1.4f\n\t"
                              "a1 = %f",
                              ct3, sin_theta3, cos_phi3, sin_phi3,
                              cos_theta4_, sin_theta4_,
                              ct5, sin_theta5, cos_phi5, sin_phi5, a1 ) );
 
   if ( fabs( pt4_+p3_lab_.px()+p5_lab_.px() ) < fabs( fabs( a1 )-pt4_ ) ) {
     DebuggingInsideLoop( Form( "|pp4+pp3*cos(phi3)+pp5*cos(phi5)| < | |a1|-pp4 |\n\t"
                                "pp4 = %f\tpp5 = %f\n\t"
                                "cos(phi3) = %f\tcos(phi5) = %f"
                                "a1 = %f",
                                pt4_, pt5, cos_phi3, cos_phi5, a1 ) );
     return true;
   }
   if ( a1 < 0. ) p5_lab_.setP( 0, -p5_lab_.px() );
   else           p3_lab_.setP( 0, -p3_lab_.px() );
   return true;
 }
 
 double
 GamGamLL::computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dw )
 {
   const double mx0 = Particle::massFromPDGId( Particle::Proton )+Particle::massFromPDGId( Particle::PiPlus ); // 1.07
   const Kinematics::Limits mx_limits = cuts_.cuts.remnants[Cuts::mass];
   const double wx2min = pow( std::max( mx0, mx_limits.min() ), 2 ),
                wx2max = pow( std::min( sqs_-outmass-2.*lepmass, mx_limits.max() ), 2 );
 
   double mx2 = 0., dmx2 = 0.;
   Map( x, wx2min, wx2max, mx2, dmx2, "mx2" );
 
   DebuggingInsideLoop( Form( "mX^2 in range (%f, %f), x = %f\n\t"
                              "mX^2 = %f, d(mX^2) = %f\n\t"
                              "mX = %f, d(mX) = %f", wx2min, wx2max, x, mx2, dmx2, sqrt( mx2 ), sqrt( dmx2 ) ) );
 
   dw = sqrt( dmx2 );
   return sqrt( mx2 );
 }
 
 void
 GamGamLL::beforeComputeWeight()
 {
   if ( !GenericProcess::is_point_set_ ) return;
 
   const Particle& p1 = event_->getOneByRole( Particle::IncomingBeam1 ),
                  &p2 = event_->getOneByRole( Particle::IncomingBeam2 );
 
   ep1_ = p1.energy();
   ep2_ = p2.energy();
 
   Ml2_ = event_->getByRole( Particle::CentralSystem )[0].mass2();
 
   switch ( cuts_.mode ) {
     case Kinematics::ElectronProton: default:
       { InError( "Case not yet supported!" ); } break;
     case Kinematics::ElasticElastic:
       dw31_ = dw52_ = 0.; break;
     case Kinematics::InelasticElastic: {
       const double m = computeOutgoingPrimaryParticlesMasses( x( 7 ), p1.mass(), sqrt( Ml2_ ), dw31_ );
       event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( m );
       event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( Particle::massFromPDGId( p2.pdgId() ) );
     } break;
     case Kinematics::ElasticInelastic: {
       const double m = computeOutgoingPrimaryParticlesMasses( x( 7 ), p2.mass(), sqrt( Ml2_ ), dw52_ );
       event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( Particle::massFromPDGId( p1.pdgId() ) );
       event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( m );
     } break;
     case Kinematics::InelasticInelastic: {
       const double mx = computeOutgoingPrimaryParticlesMasses( x( 7 ), p2.mass(), sqrt( Ml2_ ), dw31_ );
       event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( mx );
       const double my = computeOutgoingPrimaryParticlesMasses( x( 8 ), p1.mass(), sqrt( Ml2_ ), dw52_ );
       event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( my );
     } break;
   }
   MX_ = event_->getOneByRole( Particle::OutgoingBeam1 ).mass();
   MY_ = event_->getOneByRole( Particle::OutgoingBeam2 ).mass();
   MX2_ = MX_*MX_;
   MY2_ = MY_*MY_;
 }
 
 double
 GamGamLL::computeWeight()
 {
   if ( !is_outgoing_state_set_ ) { InWarning( "Output state not set!" ); return 0.; }
 
   DebuggingInsideLoop( Form( "sqrt(s)=%f\n\tm(X1)=%f\tm(X2)=%f", sqs_, MX_, MY_ ) );
 
   Kinematics::Limits& w_limits = cuts_.cuts.initial[Cuts::w];
   if ( !w_limits.hasMax() ) w_limits.max() = s_;
 
   // The minimal energy for the central system is its outgoing leptons' mass energy (or wmin_ if specified)
   const double wmin = std::max( 4.*Ml2_, w_limits.min() );
 
   // The maximal energy for the central system is its CM energy with the outgoing particles' mass energy substracted (or _wmax if specified)
   const double wmax = std::min( pow( sqs_-MX_-MY_, 2 ), w_limits.max() );
 
   DebuggingInsideLoop( Form( "wmin = %f\n\twmax = %f\n\twmax/wmin = %f", wmin, wmax, wmax/wmin ) );
 
   // compute the two-photon energy for this point
   w4_ = 0.;
   double dw4 = 0.;
   Map( x( 4 ), wmin, wmax, w4_, dw4, "w4" );
   Map( x( 4 ), wmin, wmax, w4_, dw4, "w4" );
   mc4_ = sqrt( w4_ );
 
   DebuggingInsideLoop( Form( "Computed value for w4 = %f -> mc4 = %f", w4_, mc4_ ) );
 
   if ( !orient() ) return 0.;
 
   if ( jacobian_ == 0. ) { InWarning( Form( "dj = %f", jacobian_ ) ); return 0.; }
 
   if ( t1_>0. )  { InWarning( Form( "t1 = %f > 0", t1_ ) ); return 0.; }
   if ( t2_>0. )  { InWarning( Form( "t2 = %f > 0", t2_ ) ); return 0.; }
 
   const double ecm6 = w4_ / ( 2.*mc4_ ),
                pp6cm = sqrt( ecm6*ecm6-Ml2_ );
 
   jacobian_ *= dw4*pp6cm/( mc4_*Constants::sconstb*s_ );
 
   // Let the most obscure part of this code begin...
 
   const double e1mp1 = w1_ / ( ep1_+p_cm_ ),
                e3mp3 = MX2_ / ( p3_lab_.energy()+p3_lab_.p() );
 
   const double al3 = pow( sin( p3_lab_.theta() ), 2 )/( 1.+( p3_lab_.theta() ) );
 
   // 2-photon system kinematics ?!
   const double eg = ( w4_+t1_-t2_ )/( 2.*mc4_ );
   double pg = sqrt( eg*eg-t1_ );
 
   const double pgx = -p3_lab_.px()*cos_theta4_-sin_theta4_*( de3_-e1mp1 + e3mp3 + p3_lab_.p()*al3 ),
                pgy = -p3_lab_.py(),
                pgz = mc4_*de3_/( ec4_+pc4_ )-ec4_*de3_*al4_/mc4_-p3_lab_.px()*ec4_*sin_theta4_/mc4_+ec4_*cos_theta4_/mc4_*( p3_lab_.p()*al3+e3mp3-e1mp1 );
 
   DebuggingInsideLoop( Form( "pg3 = (%f, %f, %f)\n\t"
                              "pg3^2 = %f",
                              pgx, pgy, pgz,
                              sqrt( pgx*pgx+pgy*pgy+pgz*pgz ) ) );
 
   const double pgp = sqrt( pgx*pgx + pgy*pgy ), // outgoing proton (3)'s transverse momentum
                pgg = sqrt( pgp*pgp + pgz*pgz ); // outgoing proton (3)'s momentum
   if ( pgg > pgp*0.9 && pgg > pg ) { pg = pgg; } //FIXME ???
 
   // Phi angle for the 2-photon system ?!
   const double cpg = pgx/pgp,
                spg = pgy/pgp;
 
   // Theta angle for the 2-photon system ?!
   const double stg = pgp/pg;
 
   const int theta_sign = ( pgz>0. ) ? 1 : -1;
   const double ctg = theta_sign*sqrt( 1.-stg*stg );
 
   double xx6 = x( 5 );
 
   const double amap = 0.5 * (w4_-t1_-t2_),
                bmap = 0.5 * sqrt( ( pow( w4_-t1_-t2_, 2 )-4.*t1_*t2_ )*( 1.-4.*Ml2_/w4_ ) ),
                ymap = ( amap+bmap )/( amap-bmap ),
                beta = pow( ymap, 2.*xx6-1. );
   xx6 = 0.5 * ( 1. + amap/bmap*( beta-1. )/( beta+1. ) );
   xx6 = std::max( 0., std::min( xx6, 1. ) ); // xx6 in [0., 1.]
 
   DebuggingInsideLoop( Form( "amap = %f\n\tbmap = %f\n\tymap = %f\n\tbeta = %f", amap, bmap, ymap, beta ) );
 
   // 3D rotation of the first outgoing lepton wrt the CM system
   const double theta6cm = acos( 1.-2.*xx6 );
 
   // match the Jacobian
   jacobian_ *= ( amap+bmap*cos( theta6cm ) );
   jacobian_ *= ( amap-bmap*cos( theta6cm ) );
   jacobian_ /= amap;
   jacobian_ /= bmap;
   jacobian_ *= log( ymap );
   jacobian_ *= 0.5;
 
   DebuggingInsideLoop( Form( "Jacobian = %e", jacobian_ ) );
 
   DebuggingInsideLoop( Form( "ctcm6 = %f\n\tstcm6 = %f", cos( theta6cm ), sin( theta6cm ) ) );
 
   const double phi6cm = 2.*M_PI*x( 6 );
 
   // First outgoing lepton's 3-momentum in the centre of mass system
   Particle::Momentum p6cm = Particle::Momentum::fromPThetaPhi( pp6cm, theta6cm, phi6cm );
 
   DebuggingInsideLoop( Form( "p3cm6 = (%f, %f, %f)", p6cm.px(), p6cm.py(), p6cm.pz() ) );
 
   const double h1 = stg*p6cm.pz()+ctg*p6cm.px();
   const double pc6z = ctg*p6cm.pz()-stg*p6cm.px(), pc6x = cpg*h1-spg*p6cm.py();
 
   const double qcx = 2.*pc6x, qcz = 2.*pc6z;
   // qcy == QCY is never defined
 
   const double el6 = ( ec4_*ecm6+pc4_*pc6z ) / mc4_,
                h2  = ( ec4_*pc6z+pc4_*ecm6 ) / mc4_;
 
   DebuggingInsideLoop( Form( "h1 = %f\n\th2 = %f", h1, h2 ) );
 
   // first outgoing lepton's 3-momentum
   const double p6x = cos_theta4_*pc6x+sin_theta4_*h2,
                p6y = cpg*p6cm.py()+spg*h1,
                p6z = cos_theta4_*h2-sin_theta4_*pc6x;
 
   // first outgoing lepton's kinematics
   p6_cm_ = Particle::Momentum( p6x, p6y, p6z, el6 );
   DebuggingInsideLoop( Form( "E6(cm) = %f\n\tP6(cm) = (%f, %f, %f)", el6, p6x, p6y, p6z ) );
 
   const double hq = ec4_*qcz/mc4_;
 
   const Particle::Momentum qve(
     cos_theta4_*qcx+sin_theta4_*hq,
     2.*p6y,
     cos_theta4_*hq-sin_theta4_*qcx,
     pc4_*qcz/mc4_ // energy
   );
 
   // Available energy for the second lepton is the 2-photon system's energy with the first lepton's energy removed
   const double el7 = ec4_-el6;
 
   DebuggingInsideLoop( Form( "Outgoing kinematics\n\t"
                              " first outgoing lepton: p = %f, E = %f\n\t"
                              "second outgoing lepton: p = %f, E = %f",
                              p6_cm_.p(), p6_cm_.energy(), p7_cm_.p(), p6_cm_.energy() ) );
 
   // Second outgoing lepton's 3-momentum
   const double p7x = -p6x + pt4_,
                p7y = -p6y,
                p7z = -p6z + pc4_*cos_theta4_;
 
   // second outgoing lepton's kinematics
   p7_cm_ = Particle::Momentum( p7x, p7y, p7z, el7 );
 
   //p6_cm_ = Particle::Momentum(pl6*st6*cp6, pl6*st6*sp6, pl6*ct6, el6);
   //p7_cm_ = Particle::Momentum(pl7*st7*cp7, pl7*st7*sp7, pl7*ct7, el7);
 
   q1dq_ = eg*( 2.*ecm6-mc4_ )-2.*pg*p6cm.pz();
   q1dq2_ = ( w4_-t1_-t2_ ) * 0.5;
 
   const double phi3 = p3_lab_.phi(), cos_phi3 = cos( phi3 ), sin_phi3 = sin( phi3 ),
                phi5 = p5_lab_.phi(), cos_phi5 = cos( phi5 ), sin_phi5 = sin( phi5 );
 
   bb_ = t1_*t2_+( w4_*pow( sin( theta6cm ), 2 ) + 4.*Ml2_*pow( cos( theta6cm ), 2 ) )*pg*pg;
 
   const double c1 = p3_lab_.pt() * ( qve.px()*sin_phi3  - qve.py()*cos_phi3   ),
                c2 = p3_lab_.pt() * ( qve.pz()*ep1_ - qve.energy() *p_cm_ ),
                c3 = ( w31_*ep1_*ep1_ + 2.*w1_*de3_*ep1_ - w1_*de3_*de3_ + p3_lab_.pt2()*ep1_*ep1_ ) / ( p3_lab_.energy()*p_cm_ + p3_lab_.pz()*ep1_ );
 
   const double b1 = p5_lab_.pt() * ( qve.px()*sin_phi5  - qve.py()*cos_phi5   ),
                b2 = p5_lab_.pt() * ( qve.pz()*ep2_ + qve.energy() *p_cm_ ),
                b3 = ( w52_*ep2_*ep2_ + 2.*w2_*de5_*ep2_ - w2_*de5_*de5_ + p5_lab_.pt2()*ep2_*ep2_ ) / ( ep2_*p5_lab_.pz() - p5_lab_.energy()*p_cm_ );
 
   const double r12 =  c2*sin_phi3 + qve.py()*c3,
                r13 = -c2*cos_phi3 - qve.px()*c3;
 
   const double r22 =  b2*sin_phi5 + qve.py()*b3,
                r23 = -b2*cos_phi5 - qve.px()*b3;
 
   epsi_ = p12_*c1*b1 + r12*r22 + r13*r23;
 
   g5_ = w1_*c1*c1 + r12*r12 + r13*r13;
   g6_ = w2_*b1*b1 + r22*r22 + r23*r23;
 
   a5_ = -( qve.px()*cos_phi3 + qve.py()*sin_phi3 )*p3_lab_.pt()*p1k2_-( ep1_*qve.energy()-p_cm_*qve.pz() )*( cos_phi3*cos_phi5 + sin_phi3*sin_phi5 )*p3_lab_.pt()*p5_lab_.pt() + ( de5_*qve.pz()+qve.energy()*( p_cm_+p5_lab_.pz() ) )*c3;
   a6_ = -( qve.px()*cos_phi5 + qve.py()*sin_phi5 )*p5_lab_.pt()*p2k1_-( ep2_*qve.energy()+p_cm_*qve.pz() )*( cos_phi3*cos_phi5 + sin_phi3*sin_phi5 )*p3_lab_.pt()*p5_lab_.pt() + ( de3_*qve.pz()-qve.energy()*( p_cm_-p3_lab_.pz() ) )*b3;
 
   DebuggingInsideLoop( Form( "a5 = %f\n\ta6 = %f", a5_, a6_ ) );
 
   ////////////////////////////////////////////////////////////////
   // END of GAMGAMLL subroutine in the FORTRAN version
   ////////////////////////////////////////////////////////////////
 
   const Particle::Momentum cm = event_->getOneByRole( Particle::IncomingBeam1 ).momentum() + event_->getOneByRole( Particle::IncomingBeam2 ).momentum();
 
   ////////////////////////////////////////////////////////////////
   // INFO from f.f
   ////////////////////////////////////////////////////////////////
 
   const double gamma = cm.energy() / sqs_, betgam = cm.pz() / sqs_;
 
   //--- kinematics computation for both leptons
 
   p6_cm_.betaGammaBoost( gamma, betgam );
   p7_cm_.betaGammaBoost( gamma, betgam );
 
   //--- cut on mass of final hadronic system (MX/Y)
 
   const Kinematics::Limits mx_limits = cuts_.cuts.remnants[Cuts::mass];
   if ( cuts_.mode == Kinematics::InelasticElastic || cuts_.mode == Kinematics::InelasticInelastic ) {
     if ( mx_limits.hasMin() && MX_ < mx_limits.min() ) return 0.;
     if ( mx_limits.hasMax() && MX_ > mx_limits.max() ) return 0.;
   }
   if ( cuts_.mode == Kinematics::ElasticInelastic || cuts_.mode == Kinematics::InelasticInelastic ) {
     if ( mx_limits.hasMin() && MY_ < mx_limits.min() ) return 0.;
     if ( mx_limits.hasMax() && MY_ > mx_limits.max() ) return 0.;
   }
 
   //--- cut on the proton's Q2 (first photon propagator T1)
 
   const Kinematics::Limits q2_limits = cuts_.cuts.initial[Cuts::q2];
   if ( q2_limits.hasMax() && t1_<-q2_limits.max() ) return 0;
   if ( q2_limits.hasMin() && t1_>-q2_limits.min() ) return 0.;
 
   //--- cuts on outgoing leptons' kinematics
 
   const Kinematics::Limits m_limits = cuts_.cuts.central[Cuts::mass_sum];
   if ( m_limits.hasMin() && ( p6_cm_+p7_cm_ ).mass() < m_limits.min() ) return 0.;
   if ( m_limits.hasMax() && ( p6_cm_+p7_cm_ ).mass() > m_limits.max() ) return 0.;
 
   //----- cuts on the individual leptons
 
   const Kinematics::Limits pt_limits = cuts_.cuts.central[Cuts::pt_single];
   if ( pt_limits.hasMin() && ( p6_cm_.pt() < pt_limits.min() || p7_cm_.pt() < pt_limits.min() ) ) return 0.;
   if ( pt_limits.hasMax() && ( p6_cm_.pt() > pt_limits.max() || p7_cm_.pt() > pt_limits.max() ) ) return 0.;
 
   const Kinematics::Limits energy_limits = cuts_.cuts.central[Cuts::energy_single];
   if ( energy_limits.hasMin() && ( p6_cm_.energy() < energy_limits.min() || p7_cm_.energy() < energy_limits.min() ) ) return 0.;
   if ( energy_limits.hasMax() && ( p6_cm_.energy() > energy_limits.max() || p7_cm_.energy() > energy_limits.max() ) ) return 0.;
 
   const Kinematics::Limits eta_limits = cuts_.cuts.central[Cuts::eta_single];
   if ( eta_limits.hasMin() && ( p6_cm_.eta() < eta_limits.min() || p7_cm_.eta() < eta_limits.min() ) ) return 0.;
   if ( eta_limits.hasMax() && ( p6_cm_.eta() > eta_limits.max() || p7_cm_.eta() > eta_limits.max() ) ) return 0.;
 
   //--- compute the structure functions factors
 
   switch ( cuts_.mode ) { // inherited from CDF version
     case Kinematics::ElectronProton: default: jacobian_ *= periPP( 1, 2 ); break; // ep case
     case Kinematics::ElasticElastic:          jacobian_ *= periPP( 2, 2 ); break; // elastic case
     case Kinematics::InelasticElastic:        jacobian_ *= periPP( 3, 2 )*( dw31_*dw31_ ); break;
     case Kinematics::ElasticInelastic:        jacobian_ *= periPP( 3, 2 )*( dw52_*dw52_ ); break; // single-dissociative case
     case Kinematics::InelasticInelastic:      jacobian_ *= periPP( 3, 3 )*( dw31_*dw31_ )*( dw52_*dw52_ ); break; // double-dissociative case
   }
 
   //--- compute the event weight using the Jacobian
 
   return Constants::GeV2toBarn*jacobian_;
 }
 
 void
 GamGamLL::fillKinematics( bool )
 {
   const Particle::Momentum cm = event_->getOneByRole( Particle::IncomingBeam1 ).momentum() + event_->getOneByRole( Particle::IncomingBeam2 ).momentum();
 
   const double gamma  = cm.energy()/sqs_, betgam = cm.pz()/sqs_;
 
   Particle::Momentum plab_ip1 = Particle::Momentum( 0., 0.,  p_cm_, ep1_ ).betaGammaBoost( gamma, betgam );
   Particle::Momentum plab_ip2 = Particle::Momentum( 0., 0., -p_cm_, ep2_ ).betaGammaBoost( gamma, betgam );
   p3_lab_.betaGammaBoost( gamma, betgam );
   p5_lab_.betaGammaBoost( gamma, betgam );
 
   //----- needed to parametrise a random rotation around z-axis
   const int rany = ( rand() >= RAND_MAX*0.5 ) ? 1 : -1,
             ransign = ( rand() >= RAND_MAX*0.5 ) ? 1 : -1;
   const double ranphi = ( rand()/RAND_MAX*0.5 )*2.*M_PI;
 
   Particle::Momentum plab_ph1 = ( plab_ip1-p3_lab_ ).rotatePhi( ranphi, rany );
   Particle::Momentum plab_ph2 = ( plab_ip2-p5_lab_ ).rotatePhi( ranphi, rany );
 
   p3_lab_.rotatePhi( ranphi, rany );
   p5_lab_.rotatePhi( ranphi, rany );
   p6_cm_.rotatePhi( ranphi, rany );
   p7_cm_.rotatePhi( ranphi, rany );
 
   /*if ( symmetrise_ && rand() >= .5*RAND_MAX ) {
     p6_cm_.mirrorZ();
     p7_cm_.mirrorZ();
   }*/
 
   //----- incoming protons
   event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( plab_ip1 );
   event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( plab_ip2 );
 
   //----- first outgoing proton
   Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 );
 
   op1.setMomentum( p3_lab_ );
   switch ( cuts_.mode ) {
     case Kinematics::ElasticElastic:
     case Kinematics::ElasticInelastic:
     default:
       op1.setStatus( Particle::FinalState ); // stable proton
       break;
     case Kinematics::InelasticElastic:
     case Kinematics::InelasticInelastic:
-      op1.setStatus( Particle::Undecayed ); // fragmenting remnants
+      op1.setStatus( Particle::Unfragmented ); // fragmenting remnants
       op1.setMass( MX_ );
       break;
   }
 
   //----- second outgoing proton
   Particle& op2 = event_->getOneByRole( Particle::OutgoingBeam2 );
   op2.setMomentum( p5_lab_ );
   switch ( cuts_.mode ) {
     case Kinematics::ElasticElastic:
     case Kinematics::InelasticElastic:
     default:
       op2.setStatus( Particle::FinalState ); // stable proton
       break;
     case Kinematics::ElasticInelastic:
     case Kinematics::InelasticInelastic:
-      op2.setStatus( Particle::Undecayed ); // fragmenting remnants
+      op2.setStatus( Particle::Unfragmented ); // fragmenting remnants
       op2.setMass( MY_ );
       break;
   }
 
   //----- first incoming photon
   Particle& ph1 = event_->getOneByRole( Particle::Parton1 );
   ph1.setMomentum( plab_ph1 );
   ph1.setStatus( Particle::Incoming );
 
   //----- second incoming photon
   Particle& ph2 = event_->getOneByRole( Particle::Parton2 );
   ph2.setMomentum( plab_ph2 );
   ph2.setStatus( Particle::Incoming );
 
   Particles& central_system = event_->getByRole( Particle::CentralSystem );
 
   //----- first outgoing lepton
   Particle& ol1 = central_system[0];
   ol1.setPdgId( ol1.pdgId(), ransign );
   ol1.setMomentum( p6_cm_ );
   ol1.setStatus( Particle::FinalState );
 
   //----- second outgoing lepton
   Particle& ol2 = central_system[1];
   ol2.setPdgId( ol2.pdgId(), -ransign );
   ol2.setMomentum( p7_cm_ );
   ol2.setStatus( Particle::FinalState );
 
   //----- intermediate two-lepton system
   event_->getOneByRole( Particle::Intermediate ).setMomentum( p6_cm_+p7_cm_ );
 }
 
 double
 GamGamLL::periPP( int nup_, int ndown_ )
 {
   DebuggingInsideLoop( Form( " Nup  = %d\n\tNdown = %d", nup_, ndown_ ) );
 
   FormFactors fp1, fp2;
   GenericProcess::formFactors( -t1_, -t2_, fp1, fp2 );
 
   DebuggingInsideLoop( Form( "u1 = %f\n\tu2 = %f\n\tv1 = %f\n\tv2 = %f", fp1.FM, fp1.FE, fp2.FM, fp2.FE ) );
 
   const double qqq = q1dq_*q1dq_,
                qdq = 4.*Ml2_-w4_;
   const double t11 = 64. *(  bb_*( qqq-g4_-qdq*( t1_+t2_+2.*Ml2_ ) )-2.*( t1_+2.*Ml2_ )*( t2_+2.*Ml2_ )*qqq ) * t1_*t2_, // magnetic-magnetic
                t12 = 128.*( -bb_*( dd2_+g6_ )-2.*( t1_+2.*Ml2_ )*( sa2_*qqq+a6_*a6_ ) ) * t1_, // electric-magnetic
                t21 = 128.*( -bb_*( dd4_+g5_ )-2.*( t2_+2.*Ml2_ )*( sa1_*qqq+a5_*a5_ ) ) * t2_, // magnetic-electric
                t22 = 512.*(  bb_*( delta_*delta_-gram_ )-pow( epsi_-delta_*( qdq+q1dq2_ ), 2 )-sa1_*a6_*a6_-sa2_*a5_*a5_-sa1_*sa2_*qqq ); // electric-electric
 
   const double peripp = ( fp1.FM*fp2.FM*t11 + fp1.FE*fp2.FM*t21 + fp1.FM*fp2.FE*t12 + fp1.FE*fp2.FE*t22 ) / pow( 2.*t1_*t2_*bb_, 2 );
 
   DebuggingInsideLoop( Form( "t11 = %5.2f\tt12 = %5.2f\n\t"
                              "t21 = %5.2f\tt22 = %5.2f\n\t"
                              "=> PeriPP = %e", t11, t12, t21, t22, peripp ) );
 
   return peripp;
 }
diff --git a/CepGen/Processes/GenericKTProcess.cpp b/CepGen/Processes/GenericKTProcess.cpp
index f186176..82e35a2 100644
--- a/CepGen/Processes/GenericKTProcess.cpp
+++ b/CepGen/Processes/GenericKTProcess.cpp
@@ -1,269 +1,269 @@
 #include "GenericKTProcess.h"
 #include "CepGen/StructureFunctions/StructureFunctionsBuilder.h"
 #include "CepGen/Core/Exception.h"
 
 namespace CepGen
 {
   namespace Process
   {
     GenericKTProcess::GenericKTProcess( const std::string& name,
                                         const std::string& description,
                                         const unsigned int& num_user_dimensions,
                                         const std::array<Particle::ParticleCode,2>& partons,
                                         const std::vector<Particle::ParticleCode>& central ) :
       GenericProcess( name, description+" (kT-factorisation approach)" ),
       kNumUserDimensions( num_user_dimensions ),
       kIntermediateParts( partons ), kProducedParts( central )
     {}
 
     GenericKTProcess::~GenericKTProcess()
     {}
 
     void
     GenericKTProcess::prepareKTKinematics()
     {
       DebuggingInsideLoop("Dummy kinematics prepared!");
     }
 
     void
     GenericKTProcess::addEventContent()
     {
       GenericProcess::setEventContent(
         { // incoming state
           { Particle::IncomingBeam1, Particle::Proton },
           { Particle::IncomingBeam2, Particle::Proton },
           { Particle::Parton1, kIntermediateParts[0] },
           { Particle::Parton2, kIntermediateParts[1] }
         },
         { // outgoing state
           { Particle::OutgoingBeam1, { Particle::Proton } },
           { Particle::OutgoingBeam2, { Particle::Proton } },
           { Particle::CentralSystem, kProducedParts }
         }
       );
       setExtraContent();
     }
 
     unsigned int
     GenericKTProcess::numDimensions( const Kinematics::ProcessMode& process_mode ) const
     {
       switch ( process_mode ) {
         default:
         case Kinematics::ElasticElastic:     return kNumRequiredDimensions+kNumUserDimensions;
         case Kinematics::ElasticInelastic:
         case Kinematics::InelasticElastic:   return kNumRequiredDimensions+kNumUserDimensions+1;
         case Kinematics::InelasticInelastic: return kNumRequiredDimensions+kNumUserDimensions+2;
       }
     }
 
     void
     GenericKTProcess::addPartonContent()
     {
       // Incoming partons
       qt1_ = exp( log_qmin_+( log_qmax_-log_qmin_ )*x( 0 ) );
       qt2_ = exp( log_qmin_+( log_qmax_-log_qmin_ )*x( 1 ) );
       phi_qt1_ = 2.*M_PI*x( 2 );
       phi_qt2_ = 2.*M_PI*x( 3 );
       DebuggingInsideLoop( Form( "photons transverse virtualities (qt):\n\t"
                                  "  mag = %f / %f (%.2f < log(qt) < %.2f)\n\t"
                                  "  phi = %f / %f",
                                  qt1_, qt2_, log_qmin_, log_qmax_, phi_qt1_, phi_qt2_ ) );
     }
 
     void
     GenericKTProcess::setKinematics( const Kinematics& kin )
     {
       cuts_ = kin;
       log_qmin_ = -10.; // FIXME //log_qmin_ = std::log( std::sqrt( cuts_.q2min ) );
       log_qmax_ = log( cuts_.cuts.initial[Cuts::qt].max() );
     }
 
     double
     GenericKTProcess::computeWeight()
     {
       addPartonContent();
       prepareKTKinematics();
       computeOutgoingPrimaryParticlesMasses();
 
       const double jac = computeJacobian(),
                    integrand = computeKTFactorisedMatrixElement(),
                    weight = jac*integrand;
       DebuggingInsideLoop( Form( "Jacobian = %f\n\tIntegrand = %f\n\tdW = %f", jac, integrand, weight ) );
 
       return weight;
     }
 
     void
     GenericKTProcess::computeOutgoingPrimaryParticlesMasses()
     {
       const unsigned int op_index = kNumRequiredDimensions+kNumUserDimensions;
       const Kinematics::Limits remn_mx_cuts = cuts_.cuts.remnants[Cuts::mass];
       switch ( cuts_.mode ) {
         case Kinematics::ElectronProton: default: {
           InError( "This kT factorisation process is intended for p-on-p collisions! Aborting!" );
           exit( 0 ); } break;
         case Kinematics::ElasticElastic: {
           MX_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass();
           MY_ = event_->getOneByRole( Particle::IncomingBeam2 ).mass();
         } break;
         case Kinematics::ElasticInelastic: {
           const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range();
           MX_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass();
           MY_ = mx_min + mx_range*x( op_index );
         } break;
         case Kinematics::InelasticElastic: {
           const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range();
           MX_ = mx_min + mx_range*x( op_index );
           MY_ = event_->getOneByRole( Particle::IncomingBeam2 ).mass();
         } break;
         case Kinematics::InelasticInelastic: {
           const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range();
           MX_ = mx_min + mx_range*x( op_index );
           MY_ = mx_min + mx_range*x( op_index+1 );
         } break;
       }
       DebuggingInsideLoop( Form( "outgoing remnants invariant mass: %f / %f (%.2f < M(X/Y) < %.2f)", MX_, MY_, remn_mx_cuts.min(), remn_mx_cuts.max() ) );
     }
 
     void
     GenericKTProcess::computeIncomingFluxes( double x1, double q1t2, double x2, double q2t2 )
     {
       flux1_ = flux2_ = 0.;
       switch ( cuts_.mode ) {
         case Kinematics::ElasticElastic:
           flux1_ = elasticFlux( x1, q1t2 );
           flux2_ = elasticFlux( x2, q2t2 );
           break;
         case Kinematics::ElasticInelastic:
           flux1_ = elasticFlux( x1, q1t2 );
           flux2_ = inelasticFlux( x2, q2t2, MY_ );
           break;
         case Kinematics::InelasticElastic:
           flux1_ = inelasticFlux( x1, q1t2, MX_ );
           flux2_ = elasticFlux( x2, q2t2 );
           break;
         case Kinematics::InelasticInelastic:
           flux1_ = inelasticFlux( x1, q1t2, MX_ );
           flux2_ = inelasticFlux( x2, q2t2, MY_ );
           break;
         default: return;
       }
       flux1_ = std::max( flux1_, 1.e-20 );
       flux2_ = std::max( flux2_, 1.e-20 );
       DebuggingInsideLoop( Form( "Form factors: %e / %e", flux1_, flux2_ ) );
     }
 
     void
     GenericKTProcess::fillKinematics( bool )
     {
       fillPrimaryParticlesKinematics();
       fillCentralParticlesKinematics();
     }
 
     void
     GenericKTProcess::fillPrimaryParticlesKinematics()
     {
       //=================================================================
       //     outgoing protons
       //=================================================================
       Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 ),
                &op2 = event_->getOneByRole( Particle::OutgoingBeam2 );
 
       op1.setMomentum( PX_ );
       op2.setMomentum( PY_ );
 
       switch ( cuts_.mode ) {
         case Kinematics::ElasticElastic:
           op1.setStatus( Particle::FinalState );
           op2.setStatus( Particle::FinalState );
           break;
         case Kinematics::ElasticInelastic:
           op1.setStatus( Particle::FinalState );
-          op2.setStatus( Particle::Undecayed ); op2.setMass( MY_ );
+          op2.setStatus( Particle::Unfragmented ); op2.setMass( MY_ );
           break;
         case Kinematics::InelasticElastic:
-          op1.setStatus( Particle::Undecayed ); op1.setMass( MX_ );
+          op1.setStatus( Particle::Unfragmented ); op1.setMass( MX_ );
           op2.setStatus( Particle::FinalState );
           break;
         case Kinematics::InelasticInelastic:
-          op1.setStatus( Particle::Undecayed ); op1.setMass( MX_ );
-          op2.setStatus( Particle::Undecayed ); op2.setMass( MY_ );
+          op1.setStatus( Particle::Unfragmented ); op1.setMass( MX_ );
+          op2.setStatus( Particle::Unfragmented ); op2.setMass( MY_ );
           break;    
         default: { FatalError( "This kT factorisation process is intended for p-on-p collisions! Aborting!" ); } break;
       }
 
       //=================================================================
       //     incoming partons (photons, pomerons, ...)
       //=================================================================
       //FIXME ensure the validity of this approach
       Particle& g1 = event_->getOneByRole( Particle::Parton1 ),
                &g2 = event_->getOneByRole( Particle::Parton2 );
       g1.setMomentum( event_->getOneByRole( Particle::IncomingBeam1 ).momentum()-PX_, true );
       g1.setStatus( Particle::Incoming );
       g2.setMomentum( event_->getOneByRole( Particle::IncomingBeam2 ).momentum()-PY_, true );
       g2.setStatus( Particle::Incoming );
     }
 
     double
     GenericKTProcess::minimalJacobian() const
     {
       double jac = 1.;
       jac *= ( log_qmax_-log_qmin_ )*qt1_; // d(q1t) . q1t
       jac *= ( log_qmax_-log_qmin_ )*qt2_; // d(q2t) . q2t
       jac *= 2.*M_PI; // d(phi1)
       jac *= 2.*M_PI; // d(phi2)
 
       const double mx_range = cuts_.cuts.remnants.at( Cuts::mass ).range();
       switch ( cuts_.mode ) {
         case Kinematics::ElasticElastic: default: break;
         case Kinematics::ElasticInelastic:   jac *= 2.* mx_range * MY_; break;
         case Kinematics::InelasticElastic:   jac *= 2.* mx_range * MX_; break;
         case Kinematics::InelasticInelastic: jac *= 2.* mx_range * MX_;
                                              jac *= 2.* mx_range * MY_; break;
       } // d(mx/y**2)
       return jac;
     }
 
     double
     GenericKTProcess::elasticFlux( double x, double kt2 )
     {
       const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
 
       const double Q2_min = x*x*mp2/( 1.-x ), Q2_ela = Q2_min + kt2/( 1.-x );
       const FormFactors ela = FormFactors::ProtonElastic( Q2_ela );
 
       const double ela1 = ( 1.-x )*( 1.-Q2_min/Q2_ela );
       const double ela2 = ela.FE;
       //const double ela3 = 1.-( Q2_ela-kt2 )/Q2_ela;
       const double f_ela = Constants::alphaEM/M_PI/kt2*( ela1*ela2 + 0.5*x*x*ela.FM );
 
       return f_ela * ( 1.-x ) * kt2 / Q2_ela;
     }
 
 
     double
     GenericKTProcess::inelasticFlux( double x, double kt2, double mx )
     {
       const double mx2 = mx*mx, mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
 
       // F2 structure function
       const double Q2min = 1. / ( 1.-x )*( x*( mx2-mp2 ) + x*x*mp2 ),
                    Q2 = Q2min + kt2 / ( 1.-x );
       float xbj = Q2 / ( Q2 + mx2 - mp2 );
 
       const StructureFunctions sf = StructureFunctionsBuilder::get( cuts_.structure_functions, Q2, xbj );
 
       // Longitudinal/transverse virtual photon cross section R
       // from Sibirtsev and Blunden (Phys Rev C 88,065202 (2013))
       const double R = 0.014 * Q2 * ( exp( -0.07*Q2 ) + 41.*exp( -0.8*Q2 ) );
       const double F1 = sf.F2 * ( 1.+4*xbj*xbj*mp2/Q2 )/( 1.+R )/( 2.*xbj );
 
       const double ine1 = ( 1.-x )*( 1.-Q2min / Q2 );
       const double f_D = sf.F2/( Q2 + mx2 - mp2 ) * ine1, f_C= 2.*F1/Q2;
 
       const double f_ine = Constants::alphaEM/M_PI/kt2*( f_D + 0.5*x*x*f_C );
 
       return f_ine * ( 1.-x ) * kt2 / Q2;
     }
   }
 }
diff --git a/CepGen/Processes/GenericProcess.cpp b/CepGen/Processes/GenericProcess.cpp
index 686f190..b3a988c 100644
--- a/CepGen/Processes/GenericProcess.cpp
+++ b/CepGen/Processes/GenericProcess.cpp
@@ -1,181 +1,187 @@
 #include "GenericProcess.h"
 #include "CepGen/Core/Exception.h"
 
 namespace CepGen
 {
   namespace Process
   {
     GenericProcess::GenericProcess( const std::string& name, const std::string& description, bool has_event ) :
       event_( std::shared_ptr<Event>( new Event ) ),
       is_point_set_( false ), is_incoming_state_set_( false ), is_outgoing_state_set_( false ), is_kinematics_set_( false ),
       name_( name ), description_( description ),
       total_gen_time_( 0. ), num_gen_events_( 0 ), has_event_( has_event )
     {}
 
     GenericProcess::~GenericProcess()
     {}
 
     void
     GenericProcess::setPoint( const unsigned int ndim, double* x )
     {
       if ( ndim != x_.size() ) x_.resize( ndim );
 
       x_ = std::vector<double>( x, x+ndim );
       is_point_set_ = true;
       if ( Logger::get().level>=Logger::DebugInsideLoop ) { dumpPoint(); }
     }
 
     void
     GenericProcess::prepareKinematics()
     {
       if ( !isKinematicsDefined() ) return; // FIXME dump some information...
       const Particle ib1 = event_->getOneByRole( Particle::IncomingBeam1 ),
                      ib2 = event_->getOneByRole( Particle::IncomingBeam2 );
 
       sqs_ = CMEnergy( ib1, ib2 );
       s_ = sqs_*sqs_;
 
       w1_ = ib1.mass2();
       w2_ = ib2.mass2();
 
       Debugging( Form( "Kinematics successfully prepared! sqrt(s) = %.2f", sqs_ ) );
     }
 
     void
     GenericProcess::dumpPoint()
     {
       std::ostringstream os;
       for ( unsigned int i=0; i<x_.size(); i++ ) {
         os << Form( "  x(%2d) = %8.6f\n\t", i, x_[i] );
       }
       Information( Form( "Number of integration parameters: %d\n\t"
                          "%s", x_.size(), os.str().c_str() ) );
     }
 
     void
     GenericProcess::setEventContent( const IncomingState& is, const OutgoingState& os )
     {
       event_->clear();
       //----- add the particles in the event
 
       //--- incoming state
       for ( IncomingState::const_iterator ip=is.begin(); ip!=is.end(); ip++ ) {
-        event_->addParticle( Particle( ip->first, ip->second ) );
+        std::cout << ip->first << "\t" << ip->second << std::endl;
+        Particle& p = event_->addParticle( ip->first );
+        p.setPdgId( ip->second );
+        p.dump();
       }
       //--- central system (if not already there)
       IncomingState::const_iterator central_system = is.find( Particle::CentralSystem );
       if ( central_system == is.end() ) {
-        event_->addParticle( Particle( Particle::Intermediate, Particle::invalidParticle, Particle::Propagator ) );
+        Particle& p = event_->addParticle( Particle::Intermediate );
+        p.setPdgId( Particle::invalidParticle );
+        p.setStatus( Particle::Propagator );
       }
       //--- outgoing state
       for ( OutgoingState::const_iterator op = os.begin(); op != os.end(); ++op ) {
         for ( std::vector<Particle::ParticleCode>::const_iterator it = op->second.begin(); it != op->second.end(); ++it ) {
-          event_->addParticle( Particle( op->first, *it ) );
+          Particle& p = event_->addParticle( op->first );
+          p.setPdgId( *it );
         }
       }
 
       //----- define the particles parentage
 
       const Particles parts = event_->particles();
       for ( Particles::const_iterator p = parts.begin(); p != parts.end(); ++p ) {
         Particle& part = event_->getById( p->id() );
         switch ( part.role() ) {
           case Particle::OutgoingBeam1:
           case Particle::Parton1:
             part.addMother( event_->getOneByRole( Particle::IncomingBeam1 ) );
             break;
           case Particle::OutgoingBeam2:
           case Particle::Parton2:
             part.addMother( event_->getOneByRole( Particle::IncomingBeam2 ) );
             break;
           case Particle::Intermediate:
             part.addMother( event_->getOneByRole( Particle::Parton1 ) );
             part.addMother( event_->getOneByRole( Particle::Parton2 ) );
             break;
           case Particle::CentralSystem:
             part.addMother( event_->getOneByRole( Particle::Intermediate ) );
             break;
           default: break;
         }
       }
 
       //----- freeze the event as it is
 
-      event_->init();
+      event_->freeze();
     }
 
     void
     GenericProcess::setIncomingKinematics( const Particle::Momentum& p1, const Particle::Momentum& p2 )
     {
       event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( p1 );
       event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( p2 );
     }
 
     void
     GenericProcess::formFactors( double q1, double q2, FormFactors& fp1, FormFactors& fp2 ) const
     {
       const double mx2 = MX_*MX_, my2 = MY_*MY_;
 
       switch ( cuts_.mode ) {
         case Kinematics::ElectronElectron: {
           fp1 = FormFactors::Trivial(); // electron (trivial) form factor
           fp2 = FormFactors::Trivial(); // electron (trivial) form factor
         } break;
         case Kinematics::ProtonElectron: {
           fp1 = FormFactors::ProtonElastic( -t1_ ); // proton elastic form factor
           fp2 = FormFactors::Trivial(); // electron (trivial) form factor
         } break;
         case Kinematics::ElectronProton: {
           fp1 = FormFactors::Trivial(); // electron (trivial) form factor
           fp2 = FormFactors::ProtonElastic( -t2_ ); // proton elastic form factor
         } break;
         case Kinematics::ElasticElastic: {
           fp1 = FormFactors::ProtonElastic( -t1_ ); // proton elastic form factor
           fp2 = FormFactors::ProtonElastic( -t2_ ); // proton elastic form factor
         } break;
         case Kinematics::ElasticInelastic: {
           fp1 = FormFactors::ProtonElastic( -t1_ );
           fp2 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t2_, w2_, my2 );
         } break;
         case Kinematics::InelasticElastic: {
           fp1 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t1_, w1_, mx2 );
           fp2 = FormFactors::ProtonElastic( -t2_ );
         } break;
         case Kinematics::InelasticInelastic: {
           fp1 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t1_, w1_, mx2 );
           fp2 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t2_, w2_, my2 );
         } break;
       }
     }
 
     inline bool
     GenericProcess::isKinematicsDefined()
     {
       // check the incoming state
       if ( !particles( Particle::IncomingBeam1 ).empty() && !particles( Particle::IncomingBeam2 ).empty() ) {
         is_incoming_state_set_ = true;
       }
       // check the outgoing state
       if ( !particles( Particle::OutgoingBeam1 ).empty()
         && !particles( Particle::OutgoingBeam2 ).empty()
         && !particles( Particle::CentralSystem ).empty() ) {
         is_outgoing_state_set_ = true;
       }
       // combine both states
       is_kinematics_set_ = is_incoming_state_set_ && is_outgoing_state_set_;
       return is_kinematics_set_;
     }
 
     std::ostream&
     operator<<( std::ostream& os, const GenericProcess& proc )
     {
       return os << proc.name().c_str();
     }
 
     std::ostream&
     operator<<( std::ostream& os, const GenericProcess* proc )
     {
       return os << proc->name().c_str();
     }
   }
 }