diff --git a/Cards/Config/Hadronisation/pythia6_cff.py b/Cards/Config/Hadronisation/pythia6_cff.py
index bd5e419..1f37834 100644
--- a/Cards/Config/Hadronisation/pythia6_cff.py
+++ b/Cards/Config/Hadronisation/pythia6_cff.py
@@ -1,13 +1,16 @@
from Config.containers_cfi import Module, Parameters
pythia6 = Module('pythia6',
seed = 1000,
maxTrials = 1,
preConfiguration = (
'MSTU(21)=1',
),
+ tuningSettings = (
+ 'MSEL=0', # user-defined process
+ ),
processConfiguration = (
- #'commonSettings', 'tuningSettings',
+ 'tuningSettings',
),
)
diff --git a/CepGen/Cards/PythonTypes.cpp b/CepGen/Cards/PythonTypes.cpp
index c4f07fa..07d204b 100644
--- a/CepGen/Cards/PythonTypes.cpp
+++ b/CepGen/Cards/PythonTypes.cpp
@@ -1,246 +1,231 @@
#include "CepGen/Cards/PythonHandler.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/utils.h"
#if PY_MAJOR_VERSION < 3
# define PYTHON2
#endif
namespace cepgen
{
namespace card
{
//------------------------------------------------------------------
// typed retrieval helpers
//------------------------------------------------------------------
template<> bool
PythonHandler::is<int>( PyObject* obj ) const
{
if ( !obj )
throwPythonError( "Failed to retrieve integer object!" );
#ifdef PYTHON2
return PyInt_Check( obj );
#else
return PyLong_Check( obj );
#endif
}
template<> bool
PythonHandler::is<bool>( PyObject* obj ) const
{
if ( !obj )
throwPythonError( "Failed to retrieve boolean object!" );
return PyBool_Check( obj );
}
template<> bool
PythonHandler::is<long>( PyObject* obj ) const
{
#ifdef PYTHON2
return PyInt_Check( obj ) || PyLong_Check( obj );
#else
return PyLong_Check( obj );
#endif
}
template<> int
PythonHandler::get<int>( PyObject* obj ) const
{
- if ( !obj )
- throwPythonError( "Failed to retrieve integer object!" );
if ( !is<int>( obj ) )
throw CG_ERROR( "PythonHandler:get" ) << "Object has invalid type: integer \"" << obj->ob_type->tp_name << "\".";
#ifdef PYTHON2
return PyInt_AsLong( obj );
#else
return PyLong_AsLong( obj );
#endif
}
template<> unsigned long
PythonHandler::get<unsigned long>( PyObject* obj ) const
{
- if ( !obj )
- throwPythonError( "Failed to retrieve unsigned long object!" );
if ( !is<long>( obj ) )
throw CG_ERROR( "PythonHandler:get" ) << "Object has invalid type: unsigned long \"" << obj->ob_type->tp_name << "\".";
#ifdef PYTHON2
return PyInt_AsUnsignedLongMask( obj );
#else
if ( !PyLong_Check( obj ) )
throwPythonError( Form( "Object \"%s\" has invalid type: unsigned long != %s", key, obj->ob_type->tp_name ) );
return PyLong_AsUnsignedLong( obj );
#endif
}
template<> long long
PythonHandler::get<long long>( PyObject* obj ) const
{
- if ( !obj )
- throwPythonError( "Failed to retrieve long long object!" );
if ( !is<long>( obj ) )
throw CG_ERROR( "PythonHandler:get" ) << "Object has invalid type: long long != \"" << obj->ob_type->tp_name << "\".";
return PyLong_AsLongLong( obj );
}
template<> bool
PythonHandler::is<double>( PyObject* obj ) const
{
if ( !obj )
throwPythonError( "Failed to retrieve float object!" );
return PyFloat_Check( obj );
}
template<> double
PythonHandler::get<double>( PyObject* obj ) const
{
- if ( !obj )
- throwPythonError( "Failed to retrieve float object!" );
+ CG_WARNING("");
if ( !is<double>( obj ) )
throw CG_ERROR( "PythonHandler:get" ) << "Object has invalid type: double != \"" << obj->ob_type->tp_name << "\".";
return PyFloat_AsDouble( obj );
}
template<> bool
PythonHandler::is<std::string>( PyObject* obj ) const
{
if ( !obj )
throwPythonError( "Failed to retrieve string object!" );
#ifdef PYTHON2
return PyString_Check( obj );
#else
return PyUnicode_Check( obj );
#endif
}
template<> std::string
PythonHandler::get<std::string>( PyObject* obj ) const
{
- if ( !obj )
- throwPythonError( "Failed to retrieve string object!" );
if ( !is<std::string>( obj ) )
throw CG_ERROR( "PythonHandler:get" ) << "Object has invalid type: string != \"" << obj->ob_type->tp_name << "\".";
#ifdef PYTHON2
const std::string out = PyString_AsString( obj ); // deprecated in python v3+
#else
PyObject* pstr = PyUnicode_AsEncodedString( obj, "utf-8", "strict" ); // new
if ( !pstr )
throwPythonError( "Failed to decode a Python object!" );
const std::string out = PyBytes_AS_STRING( pstr );
Py_CLEAR( pstr );
#endif
return out;
}
template<> bool
PythonHandler::is<Limits>( PyObject* obj ) const
{
if ( !obj )
throwPythonError( "Failed to retrieve limits object!" );
if ( !isVector<double>( obj ) )
return false;
const size_t size = getVector<double>( obj ).size();
return ( size == 1 || size == 2 );
}
template<> Limits
PythonHandler::get<Limits>( PyObject* obj ) const
{
- if ( !obj )
- throwPythonError( "Failed to retrieve limits object!" );
if ( !is<Limits>( obj ) )
throw CG_ERROR( "PythonHandler:get" ) << "Object has invalid type: limits != \"" << obj->ob_type->tp_name << "\".";
const auto vec = getVector<double>( obj );
if ( vec.size() == 1 )
return Limits{ vec.at( 0 ) };
return Limits{ vec.at( 0 ), vec.at( 1 ) };
}
template<> bool
PythonHandler::is<ParametersList>( PyObject* obj ) const
{
if ( !obj )
throwPythonError( "Failed to retrieve parameters list object!" );
return PyDict_Check( obj );
}
template<> ParametersList
PythonHandler::get<ParametersList>( PyObject* obj ) const
{
- if ( !obj )
- throwPythonError( "Failed to retrieve parameters list object!" );
if ( !is<ParametersList>( obj ) )
throw CG_ERROR( "PythonHandler:get" ) << "Object has invalid type: parameters list != \"" << obj->ob_type->tp_name << "\".";
ParametersList out;
PyObject* pkey = nullptr, *pvalue = nullptr;
Py_ssize_t pos = 0;
while ( PyDict_Next( obj, &pos, &pkey, &pvalue ) ) {
const std::string skey = is<std::string>( pkey )
? get<std::string>( pkey )
: is<int>( pkey )
? std::to_string( get<int>( pkey ) ) // integer-type key
: "invalid";
if ( is<int>( pvalue ) )
out.set<int>( skey, get<int>( pvalue ) );
else if ( is<double>( pvalue ) )
out.set<double>( skey, get<double>( pvalue ) );
else if ( is<std::string>( pvalue ) )
out.set<std::string>( skey, get<std::string>( pvalue ) );
else if ( is<ParametersList>( pvalue ) )
out.set<ParametersList>( skey, get<ParametersList>( pvalue ) );
else if ( is<Limits>( pvalue ) )
out.set<Limits>( skey, get<Limits>( pvalue ) );
else if ( PyTuple_Check( pvalue ) || PyList_Check( pvalue ) ) { // vector
if ( isVector<int>( pvalue ) )
out.set<std::vector<int> >( skey, getVector<int>( pvalue ) );
else if ( isVector<double>( pvalue ) )
out.set<std::vector<double> >( skey, getVector<double>( pvalue ) );
else if ( isVector<std::string>( pvalue ) )
out.set<std::vector<std::string> >( skey, getVector<std::string>( pvalue ) );
else if ( isVector<ParametersList>( pvalue ) )
out.set<std::vector<ParametersList> >( skey, getVector<ParametersList>( pvalue ) );
}
else
throwPythonError( "Invalid object retrieved as parameters list value!" );
}
return out;
}
template<typename T> bool
PythonHandler::isVector( PyObject* obj ) const
{
if ( !obj )
return false;
if ( !PyTuple_Check( obj ) && !PyList_Check( obj ) )
return false;
PyObject* pfirst = PyTuple_Check( obj ) ? PyTuple_GetItem( obj, 0 ) : PyList_GetItem( obj, 0 );
if ( !is<T>( pfirst ) )
return false;
return true;
}
template<typename T> std::vector<T>
PythonHandler::getVector( PyObject* obj ) const
{
- if ( !obj )
- throwPythonError( "Failed to retrieve list/tuple object!" );
if ( !isVector<T>( obj ) )
throw CG_ERROR( "PythonHandler:get" ) << "Object has invalid type: list/tuple != \"" << obj->ob_type->tp_name << "\".";
std::vector<T> vec;
const bool tuple = PyTuple_Check( obj );
const Py_ssize_t num_entries = tuple ? PyTuple_Size( obj ) : PyList_Size( obj );
//--- check every single element inside the list/tuple
for ( Py_ssize_t i = 0; i < num_entries; ++i ) {
PyObject* pit = tuple ? PyTuple_GetItem( obj, i ) : PyList_GetItem( obj, i );
if ( !is<T>( pit ) )
throwPythonError( "Mixed types detected in vector" );
vec.emplace_back( get<T>( pit ) );
}
return vec;
}
}
}
diff --git a/CepGen/Hadronisers/GenericHadroniser.cpp b/CepGen/Hadronisers/GenericHadroniser.cpp
index 9781bca..4ea2614 100644
--- a/CepGen/Hadronisers/GenericHadroniser.cpp
+++ b/CepGen/Hadronisers/GenericHadroniser.cpp
@@ -1,53 +1,53 @@
#include "CepGen/Hadronisers/GenericHadroniser.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/ParametersList.h"
namespace cepgen
{
namespace hadr
{
GenericHadroniser::GenericHadroniser( const ParametersList& plist, const std::string& name ) :
name_( name ),
seed_ ( plist.get<int>( "seed", -1ll ) ),
max_trials_( plist.get<int>( "maxTrials", 1 ) )
{
CG_DEBUG( "Hadroniser:init" )
<< "\"" << name_ << "\"-type hadroniser built with:\n\t"
<< "* seed = " << seed_ << "\n\t"
<< "* maximum trials: " << max_trials_;
}
void
GenericHadroniser::readStrings( const std::vector<std::string>& params ) {
if ( params.empty() )
return;
std::ostringstream os;
for ( const auto& p : params ) {
readString( p );
os << "\n\t '" << p << "'";
}
CG_DEBUG( "Hadroniser:configure" )
<< "Feeding \"" << name_ << "\" hadroniser with:"
<< os.str();
}
std::string
GenericHadroniser::name() const
{
return name_;
}
}
std::ostream&
operator<<( std::ostream& os, const hadr::GenericHadroniser& hadr )
{
- return os << hadr.name().c_str();
+ return os << hadr.name();
}
std::ostream&
operator<<( std::ostream& os, const hadr::GenericHadroniser* hadr )
{
- return os << hadr->name().c_str();
+ return os << hadr->name();
}
}
diff --git a/CepGen/Hadronisers/Pythia6Hadroniser.cpp b/CepGen/Hadronisers/Pythia6Hadroniser.cpp
index 16cfd44..4a6d035 100644
--- a/CepGen/Hadronisers/Pythia6Hadroniser.cpp
+++ b/CepGen/Hadronisers/Pythia6Hadroniser.cpp
@@ -1,341 +1,335 @@
#include "CepGen/Hadronisers/GenericHadroniser.h"
#include "CepGen/Hadronisers/HadronisersHandler.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Event/Particle.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Core/ParametersList.h" //FIXME
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/utils.h"
#include <algorithm>
extern "C"
{
/// Get the particle's mass in GeV from the Pythia6 module
extern double pymass_( int& );
/// Get the resonant particle's width in GeV from the Pythia6 module
//extern double pywidt_( int& );
/// Launch the Pythia6 fragmentation
extern void pyexec_();
/// Set a parameter value to the Pythia6 module
extern void pygive_( const char*, int );
extern void pyckbd_();
/// List all the particles in the event in a human-readable format
extern void pylist_( int& );
/// Join two coloured particles in a colour singlet
extern void pyjoin_( int&, int& );
/// Get a particle's human-readable name from the Pythia6 module
extern void pyname_( int&, char*, int );
/// Get integer-valued event information from the Pythia6 module
extern int pyk_( int&, int& );
/// Get real-valued event information from the Pythia6 module
extern double pyp_( int&, int& );
/// Purely virtual method to call at the end of the run
void pystop_() { CG_INFO( "Pythia6Hadroniser" ) << "End of run"; }
/// Particles content of the event
extern struct
{
/// Number of particles in the event
int n;
int npad;
/// Particles' general information (status, PDG id, mother, daughter 1, daughter 2)
int k[5][4000];
/// Particles' kinematics, in GeV (px, py, pz, E, M)
double p[5][4000];
/// Primary vertex for the particles
double v[5][4000];
} pyjets_;
}
namespace cepgen
{
namespace hadr
{
/**
* Full interface to the Pythia6 @cite Sjostrand:2006za 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 Pythia6 hadronisation algorithm
*/
class Pythia6Hadroniser : public GenericHadroniser
{
public:
Pythia6Hadroniser( const ParametersList& );
void setParameters( const Parameters& ) override {}
- void readString( const char* param ) override;
+ inline void readString( const char* param ) override { pygive( param ); }
void init() override {}
bool run( Event& ev, double& weight, bool full ) override;
void setCrossSection( double xsec, double xsec_err ) override {}
//bool hadronise( const Particle* );
private:
static constexpr unsigned short MAX_PART_STRING = 3;
static constexpr unsigned short MAX_STRING_EVENT = 2;
/// Maximal number of characters to fetch for the particle's name
static constexpr unsigned short NAME_CHR = 16;
inline static double pymass(int pdgid_) { return pymass_(pdgid_); }
//inline static double pywidt(int pdgid_) { return pywidt_(pdgid_); }
inline static void pyckbd() { pyckbd_(); }
inline static void pygive( const std::string& line ) { pygive_( line.c_str(), line.length() ); }
inline static void pylist( int mlist ) { pylist_( mlist ); }
inline static int pyk( int id, int qty ) { return pyk_( id, qty ); }
inline static double pyp( int id, int qty ) { return pyp_( id, qty ); }
inline static std::string pyname( int pdgid ) {
char out[NAME_CHR];
std::string s;
pyname_( pdgid, out, NAME_CHR );
s = std::string( out, NAME_CHR );
s.erase( remove( s.begin(), s.end(), ' ' ), s.end() );
return s;
}
/**
* \brief Connect entries with colour flow information
- * \param[in] njoin_ Number of particles to join in the colour flow
- * \param[in] ijoin_ List of particles unique identifier to join in the colour flow
+ * \param[in] njoin Number of particles to join in the colour flow
+ * \param[in] ijoin List of particles unique identifier to join in the colour flow
*/
inline static void pyjoin( int njoin, int ijoin[2] ) { return pyjoin_( njoin, *ijoin ); }
bool prepareHadronisation( Event& );
std::pair<short,short> pickPartonsContent() const;
struct EventProperties
{
unsigned int str_in_evt = 0;
unsigned int num_part_in_str[MAX_STRING_EVENT] = { 0 };
};
EventProperties fillParticles( const Event& ) const;
};
Pythia6Hadroniser::Pythia6Hadroniser( const ParametersList& plist ) :
GenericHadroniser( plist, "pythia6" )
- {}
-
- void
- Pythia6Hadroniser::readString( const char* param )
- {
- pygive( param );
- }
+ {CG_INFO("");}
bool
Pythia6Hadroniser::run( Event& ev, double& weight, bool full )
{
full = true; //FIXME
weight = 1.;
if ( full )
prepareHadronisation( ev );
if ( utils::Logger::get().level >= utils::Logger::Level::debug ) {
CG_DEBUG( "Pythia6Hadroniser" ) << "Dump of the event before the hadronisation:";
ev.dump();
}
//--- fill Pythia 6 common blocks
auto prop = fillParticles( ev );
CG_DEBUG( "Pythia6Hadroniser" )
<< "Passed the string construction stage.\n\t"
<< " " << prop.str_in_evt << " string objects were identified and constructed.";
unsigned int oldnpart = pyjets_.n;
pyexec_();
int criteria = oldnpart+1;
for ( unsigned int i = 0; i < MAX_STRING_EVENT; ++i )
criteria += prop.num_part_in_str[i];
if ( pyjets_.k[1][criteria] == 2212
&& pyjets_.k[0][criteria] == 1 ) {
CG_WARNING( "Pythia6Hadroniser" ) << "System is non-inelastic.";
return false;
}
// We filter the first particles already present in the event
for ( unsigned int p = oldnpart; p < (unsigned int)pyjets_.n; ++p ) {
const pdgid_t pdg_id = abs( pyjets_.k[1][p] );
const short charge = pyjets_.k[1][p]/(short)pdg_id;
ParticleProperties prop;
if ( full )
try { prop = PDG::get()( pdg_id ); } catch ( const Exception& ) {
prop = ParticleProperties{ pdg_id,
pyname( pdg_id ), pyname( pdg_id ),
(short)pyk( p+1, 12 ), // colour factor
pymass( pdg_id ), -1., //pmas( pdg_id, 2 ),
(short)pyk( p+1, 6 ), // charge
false
};
PDG::get().define( prop );
}
const Particle::Role role = pyjets_.k[2][p] != 0
? ev[pyjets_.k[2][p]-1].role() // child particle inherits its mother's role
: Particle::Role::UnknownRole;
auto& pa = ev.addParticle( role );
pa.setId( p );
pa.setPdgId( pdg_id, charge );
pa.setStatus( (Particle::Status)pyjets_.k[0][p] );
pa.setMomentum( Particle::Momentum( pyjets_.p[0][p], pyjets_.p[1][p], pyjets_.p[2][p], pyjets_.p[3][p] ) );
pa.setMass( pyjets_.p[4][p] );
if ( role != Particle::Role::UnknownRole ) {
auto& moth = ev[pyjets_.k[2][p]-1];
moth.setStatus( role == Particle::Role::CentralSystem
? Particle::Status::Resonance
: Particle::Status::Fragmented );
pa.addMother( moth );
}
}
return true;
}
bool
Pythia6Hadroniser::prepareHadronisation( Event& ev )
{
CG_DEBUG( "Pythia6Hadroniser" ) << "Hadronisation preparation called.";
for ( const auto& part : ev.particles() ) {
if ( part.status() != Particle::Status::Unfragmented )
continue;
//--- only loop over all protons to be fragmented
part.dump();
const auto partons = pickPartonsContent();
const double mx2 = part.mass2();
const double mq = pymass( partons.first ), mq2 = mq*mq;
const double mdq = pymass( partons.second ), mdq2 = mdq*mdq;
//--- choose random direction in MX frame
const double phi = 2.*M_PI*drand(), theta = acos( 2.*drand()-1. ); // theta angle
//--- compute momentum of decay particles from MX
const double px = std::sqrt( std::pow( mx2-mdq2+mq2, 2 ) / ( 4.*mx2 ) - mq2 );
//--- build 4-vectors and boost decay particles
// const auto pq = Particle::Momentum::fromPThetaPhi( px, theta, phi, std::hypot( px, mq ) );
const auto pq = Particle::Momentum::fromPThetaPhi( px, theta, phi, part.mass() );
CG_INFO("")<<pq << "|" << -pq;
//--- singlet
auto singl_mom = part.momentum();
singl_mom.lorentzBoost( pq );
auto& quark = ev.addParticle( part.role() );
quark.addMother( ev[part.id()] );
quark.setPdgId( partons.first );
quark.setStatus( Particle::Status::FinalState );
quark.setMomentum( singl_mom );
//--- doublet
auto doubl_mom = part.momentum();
doubl_mom.lorentzBoost( -pq );
auto& diquark = ev.addParticle( part.role() );
diquark.addMother( ev[part.id()] );
diquark.setPdgId( partons.second );
diquark.setStatus( Particle::Status::FinalState );
diquark.setMomentum( doubl_mom );
std::cout << part.momentum()-(singl_mom+doubl_mom) << "|" << part.numDaughters() << std::endl;
ev[part.id()].setStatus( Particle::Status::Fragmented );
}
ev.dump();
return true;
}
Pythia6Hadroniser::EventProperties
Pythia6Hadroniser::fillParticles( const Event& ev ) const
{
pyjets_.n = 0;
//--- initialising the string fragmentation variables
EventProperties out;
out.str_in_evt = 0;
int jlpsf[MAX_STRING_EVENT][MAX_PART_STRING] = { 0 };
for ( const auto& role : ev.roles() ) { // loop on roles
unsigned int part_in_str = 0;
bool role_has_string = false;
for ( const auto& part : ev[role] ) {
unsigned int np = part.id();
pyjets_.p[0][np] = part.momentum().px();
pyjets_.p[1][np] = part.momentum().py();
pyjets_.p[2][np] = part.momentum().pz();
pyjets_.p[3][np] = part.energy();
pyjets_.p[4][np] = part.mass();
pyjets_.k[0][np] = part.status() <= Particle::Status::Undefined
? 21 // incoming beam
: (int)part.status();
pyjets_.k[1][np] = part.integerPdgId();
pyjets_.k[2][np] = part.mothers().empty()
? 0 // no mother
: *( part.mothers().begin() )+1; // mother
const auto& daug = part.daughters();
if ( daug.empty() )
pyjets_.k[3][np] = pyjets_.k[4][np] = 0; // no daughters
else {
pyjets_.k[3][np] = *daug.begin()+1; // daughter 1
pyjets_.k[4][np] = *daug.rbegin()+1; // daughter 2
}
for ( int i = 0; i < 5; ++i )
pyjets_.v[i][np] = 0.;
if ( part.status() == Particle::Status::DebugResonance ) {
pyjets_.k[0][np] = 1; //FIXME PYTHIA/JETSET workaround
jlpsf[out.str_in_evt][part_in_str++] = part.id()+1;
out.num_part_in_str[out.str_in_evt]++;
role_has_string = true;
}
else if ( part.status() == Particle::Status::Undecayed )
pyjets_.k[0][np] = 2; // intermediate resonance
pyjets_.n++;
}
//--- at most one string per role
if ( role_has_string )
out.str_in_evt++;
}
//--- loop over the strings to bind everything together
for ( unsigned short i = 0; i < out.str_in_evt; ++i ) {
if ( out.num_part_in_str[i] < 2 )
continue;
std::ostringstream dbg;
for ( unsigned short j = 0; j < out.num_part_in_str[i]; ++j )
if ( jlpsf[i][j] != -1 )
dbg << Form( "\n\t * %2d (pdgId=%4d)", jlpsf[i][j], pyjets_.k[1][jlpsf[i][j]-1] );
CG_INFO( "Pythia6Hadroniser" )
<< "Joining " << out.num_part_in_str[i] << " particle" << utils::s( out.num_part_in_str[i] )
<< " with " << ev[jlpsf[i][0]].role() << " role"
<< " in a same string (id=" << i << ")"
<< dbg.str();
pyjoin( out.num_part_in_str[i], jlpsf[i] );
}
return out;
}
std::pair<short,short>
Pythia6Hadroniser::pickPartonsContent() const
{
const double ranudq = drand();
if ( ranudq < 1./9. )
return { PDG::down, 2203 }; // (d,uu1)
if ( ranudq < 5./9. )
return { PDG::up, 2101 }; // (u,ud0)
return { PDG::up, 2103 }; // (u,ud1)
}
}
}
// register hadroniser and define alias
REGISTER_HADRONISER( pythia6, Pythia6Hadroniser )