diff --git a/CepGen/CMakeLists.txt b/CepGen/CMakeLists.txt index 3d71600..bcc5061 100644 --- a/CepGen/CMakeLists.txt +++ b/CepGen/CMakeLists.txt @@ -1,126 +1,127 @@ file(GLOB core_sources Core/*.cpp *.cpp) file(GLOB phys_sources Physics/*.cpp) file(GLOB sf_sources StructureFunctions/*.cpp) file(GLOB io_sources IO/GenericExportHandler.cpp IO/TextHandler.cpp) file(GLOB hadr_sources Hadronisers/GenericHadroniser.cpp) include_directories(${PROJECT_SOURCE_DIR}) #----- check the external dependencies for SFs set(GRV_PATH ${PROJECT_SOURCE_DIR}/External) file(GLOB grv_sources ${GRV_PATH}/grv_*.f) if(grv_sources) message(STATUS "GRV PDFset found in ${grv_sources}!") add_definitions(-DGRVPDF) list(APPEND sf_sources ${grv_sources}) else() message(STATUS "GRV PDFset not found. Will proceed without it") endif() #----- check the external dependencies for physics utilities if(alphas_sources) list(APPEND phys_sources ${alphas_sources}) endif() set(addons_libraries "") set(strf_libraries ${GSL_LIB} ${GSL_CBLAS_LIB}) #--- linking with LHAPDF if(LHAPDF) message(STATUS "LHAPDF found in ${LHAPDF}") list(APPEND strf_libraries ${LHAPDF}) include_directories(${LHAPDF_INCLUDE}) else() file(GLOB partonic_sf StructureFunctions/Partonic.cpp) list(REMOVE_ITEM sf_sources ${partonic_sf}) endif() #--- linking with Pythia 6 if(PYTHIA6) message(STATUS "Pythia 6 found in ${PYTHIA6}") list(APPEND hadr_sources "Hadronisers/Pythia6Hadroniser.cpp") list(APPEND addons_libraries ${PYTHIA6}) if(PYTHIA6DUMMY) message(STATUS "Pythia 6 addons found in ${PYTHIA6DUMMY}") list(APPEND addons_libraries ${PYTHIA6DUMMY}) endif() elseif(EXISTS $ENV{PYTHIA6_SRC}) file(GLOB pythia6_src $ENV{PYTHIA6_SRC}) message(STATUS "Pythia 6 source found in ${pythia6_src}") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Wno-tabs -Wno-maybe-uninitialized -Wno-integer-division -Wno-unused-variable -Wno-unused-dummy-argument") add_library(pythia6 SHARED ${pythia6_src}) list(APPEND hadr_sources "Hadronisers/Pythia6Hadroniser.cpp") list(APPEND addons_libraries pythia6) endif() #--- linking with Pythia 8 if(PYTHIA8) message(STATUS "Pythia 8 found in ${PYTHIA8}") message(STATUS "Pythia 8 will be used for LHEF output") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-misleading-indentation") include_directories(${PYTHIA8_INCLUDE}) add_definitions(-DPYTHIA8) list(APPEND hadr_sources "Hadronisers/Pythia8Hadroniser.cpp") list(APPEND io_sources "IO/PythiaEventInterface.cpp") list(APPEND io_sources "IO/LHEFPythiaHandler.cpp") list(APPEND addons_libraries ${PYTHIA8} dl) endif() #--- linking with HepMC if(HEPMC_LIB) message(STATUS "HepMC found in ${HEPMC_LIB}") if(HEPMC_LIB MATCHES ".*HepMC3.?.so") message(STATUS "HepMC version 3 found") if(HEPMC_ROOT_LIB) message(STATUS "HepMC ROOT I/O library found") add_definitions(-DHEPMC3_ROOTIO) list(APPEND addons_libraries ${HEPMC_ROOT_LIB}) endif() add_definitions(-DHEPMC3) endif() list(APPEND io_sources "IO/HepMCHandler.cpp") + list(APPEND io_sources "IO/HepMCEventInterface.cpp") if(NOT PYTHIA8) message(STATUS "HepMC will be used for LHEF output") list(APPEND io_sources "IO/LHEFHepMCHandler.cpp") endif() list(APPEND addons_libraries ${HEPMC_LIB}) include_directories(${HEPMC_INCLUDE}) endif() #--- linking with ROOT/Delphes if(ROOT_FOUND) message(STATUS "ROOT found in ${ROOT_LIBRARY_DIR}") list(APPEND addons_libraries ${ROOT_LIBRARIES}) list(APPEND io_sources "IO/ROOTTreeHandler.cpp") include_directories(${ROOT_INCLUDE_DIRS}) link_directories(${ROOT_LIBRARY_DIR}) if(DELPHES) message(STATUS "Delphes found in ${DELPHES}") list(APPEND addons_libraries ${DELPHES}) list(APPEND io_sources "IO/DelphesHandler.cpp") include_directories(${DELPHES_INCLUDE} ${DELPHES_INCLUDE}/external) endif() endif() #----- build the objects add_library(CepGenCore SHARED ${core_sources} ${phys_sources} ${sf_sources}) target_link_libraries(CepGenCore ${CEPGEN_EXTERNAL_CORE_REQS}) target_link_libraries(CepGenCore ${strf_libraries}) add_library(CepGenAddOns SHARED ${io_sources} ${hadr_sources}) target_link_libraries(CepGenAddOns ${addons_libraries}) target_link_libraries(CepGenAddOns CepGenEvent) #----- installation rules install(TARGETS CepGenCore DESTINATION lib) install(TARGETS CepGenAddOns DESTINATION lib) diff --git a/CepGen/Event/Particle.h b/CepGen/Event/Particle.h index 70d944d..b7188dc 100644 --- a/CepGen/Event/Particle.h +++ b/CepGen/Event/Particle.h @@ -1,355 +1,355 @@ #ifndef CepGen_Event_Particle_h #define CepGen_Event_Particle_h #include "CepGen/Core/Hasher.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Physics/ParticleProperties.h" #include #include #include namespace cepgen { /// A set of integer-type particle identifiers typedef std::set ParticlesIds; /// Kinematic information for one particle class Particle { public: /// Internal status code for a particle enum class Status { PrimordialIncoming = -9, ///< Incoming beam particle DebugResonance = -5, ///< Intermediate resonance (for processes developers) Resonance = -4, ///< Already decayed intermediate resonance Fragmented = -3, ///< Already fragmented outgoing beam Propagator = -2, ///< Generic propagator Incoming = -1, ///< Incoming parton Undefined = 0, ///< Undefined particle FinalState = 1, ///< Stable, final state particle Undecayed = 2, ///< Particle to be decayed externally Unfragmented = 3 ///< Particle to be hadronised externally }; /// Role of the particle in the process enum Role { UnknownRole = -1, ///< Undefined role IncomingBeam1 = 1, ///< \f$z>0\f$ incoming beam particle IncomingBeam2 = 2, ///< \f$z<0\f$ incoming beam particle OutgoingBeam1 = 3, ///< \f$z<0\f$ outgoing beam state/particle OutgoingBeam2 = 5, ///< \f$z>0\f$ outgoing beam state/particle CentralSystem = 6, ///< Central particles system Intermediate = 4, ///< Intermediate two-parton system Parton1 = 41, ///< \f$z>0\f$ beam incoming parton Parton2 = 42 ///< \f$z<0\f$ beam incoming parton }; /** * Container for a particle's 4-momentum, along with useful methods to ease the development of any matrix element level generator * \brief 4-momentum for a particle * \date Dec 2015 * \author Laurent Forthomme */ class Momentum { public: /// Build a 4-momentum at rest with an invalid energy (no mass information known) Momentum(); /// Build a 4-momentum using its 3-momentum coordinates and its energy Momentum( double x, double y, double z, double t = -1. ); /// Build a 4-momentum using its 3-momentum coordinates and its energy Momentum( double* p ); // --- static definitions /// Build a 3-momentum from its three pseudo-cylindric coordinates static Momentum fromPtEtaPhi( double pt, double eta, double phi, double e = -1. ); /// Build a 4-momentum from its scalar momentum, and its polar and azimuthal angles static Momentum fromPThetaPhi( double p, double theta, double phi, double e = -1. ); /// Build a 4-momentum from its four momentum and energy coordinates static Momentum fromPxPyPzE( double px, double py, double pz, double e ); /// Build a 4-momentum from its three momentum coordinates and mass static Momentum fromPxPyPzM( double px, double py, double pz, double m ); /// Build a 4-momentum from its transverse momentum, rapidity and mass static Momentum fromPxPyYM( double px, double py, double rap, double m ); // --- vector and scalar operators /// Scalar product of the 3-momentum with another 3-momentum double threeProduct( const Momentum& ) const; /// Scalar product of the 4-momentum with another 4-momentum double fourProduct( const Momentum& ) const; /// Vector product of the 3-momentum with another 3-momentum double crossProduct( const Momentum& ) const; /// Compute the 4-vector sum of two 4-momenta Momentum operator+( const Momentum& ) const; /// Add a 4-momentum through a 4-vector sum Momentum& operator+=( const Momentum& ); /// Unary inverse operator Momentum operator-() const; /// Compute the inverse per-coordinate 4-vector Momentum operator-( const Momentum& ) const; /// Subtract a 4-momentum through a 4-vector sum Momentum& operator-=( const Momentum& ); /// Scalar product of two 3-momenta double operator*( const Momentum& ) const; /// Scalar product of the 3-momentum with another 3-momentum double operator*=( const Momentum& ); /// Multiply all components of a 4-momentum by a scalar Momentum operator*( double c ) const; /// Multiply all 4-momentum coordinates by a scalar Momentum& operator*=( double c ); /// Left-multiply all 4-momentum coordinates by a scalar friend Momentum operator*( double, const Momentum& ); /// Equality operator bool operator==( const Momentum& ) const; /// Human-readable format for a particle's momentum friend std::ostream& operator<<( std::ostream&, const Momentum& ); Momentum& betaGammaBoost( double gamma, double betagamma ); /// Forward Lorentz boost Momentum& lorentzBoost( const Momentum& p ); // --- setters and getters /// Set all the components of the 4-momentum (in GeV) Momentum& setP( double px, double py, double pz, double e ); /// Set all the components of the 3-momentum (in GeV) Momentum& setP( double px, double py, double pz ); /// Set the energy (in GeV) inline Momentum& setEnergy( double e ) { energy_ = e; return *this; } /// Compute the energy from the mass inline Momentum& setMass( double m ) { return setMass2( m*m ); } /// Compute the energy from the mass Momentum& setMass2( double m2 ); /// Get one component of the 4-momentum (in GeV) double operator[]( const unsigned int i ) const; /// Get one component of the 4-momentum (in GeV) double& operator[]( const unsigned int i ); /// Momentum along the \f$x\f$-axis (in GeV) inline double px() const { return px_; } /// Momentum along the \f$y\f$-axis (in GeV) inline double py() const { return py_; } /// Longitudinal momentum (in GeV) inline double pz() const { return pz_; } /// Transverse momentum (in GeV) double pt() const; /// Squared transverse momentum (in GeV\f$^2\f$) double pt2() const; /// 5-vector of double precision floats (in GeV) const std::vector pVector() const; /// 3-momentum norm (in GeV) inline double p() const { return p_; } /// Squared 3-momentum norm (in GeV\f$^2\f$) inline double p2() const { return p_*p_; } /// Energy (in GeV) inline double energy() const { return energy_; } /// Squared energy (in GeV\f$^2\f$) inline double energy2() const { return energy_*energy_; } /// Squared mass (in GeV\f$^2\f$) as computed from its energy and momentum inline double mass2() const { return energy2()-p2(); } /// Mass (in GeV) as computed from its energy and momentum /// \note Returns \f$-\sqrt{|E^2-\mathbf{p}^2|}<0\f$ if \f$\mathbf{p}^2>E^2\f$ double mass() const; /// Polar angle (angle with respect to the longitudinal direction) double theta() const; /// Azimutal angle (angle in the transverse plane) double phi() const; /// Pseudo-rapidity double eta() const; /// Rapidity double rapidity() const; Momentum& truncate( double tolerance = 1.e-10 ); /// Rotate the transverse components by an angle phi (and reflect the y coordinate) Momentum& rotatePhi( double phi, double sign ); /// Rotate the particle's momentum by a polar/azimuthal angle Momentum& rotateThetaPhi( double theta_, double phi_ ); /// Apply a \f$ z\rightarrow -z\f$ transformation inline Momentum& mirrorZ() { pz_ = -pz_; return *this; } private: /// Compute the 3-momentum's norm Momentum& computeP(); /// Momentum along the \f$x\f$-axis double px_; /// Momentum along the \f$y\f$-axis double py_; /// Momentum along the \f$z\f$-axis double pz_; /// 3-momentum's norm (in GeV/c) double p_; /// Energy (in GeV) double energy_; }; /// Human-readable format for a particle's role in the event friend std::ostream& operator<<( std::ostream& os, const Role& rl ); //----- static getters /// Convert a polar angle to a pseudo-rapidity static double thetaToEta( double theta ); /// Convert a pseudo-rapidity to a polar angle static double etaToTheta( double eta ); /// Convert a pseudo-rapidity to a rapidity static double etaToY( double eta_, double m_, double pt_ ); Particle(); /// Build using the role of the particle in the process and its PDG id /// \param[in] role Role of the particle in the process /// \param[in] id PDG identifier /// \param[in] st Current status Particle( Role role, pdgid_t id, Status st = Status::Undefined ); /// Copy constructor Particle( const Particle& ); inline ~Particle() {} /// Comparison operator (from unique identifier) bool operator<( const Particle& rhs ) const; /// Comparison operator (from their reference's unique identifier) //bool operator<( Particle *rhs ) const { return ( id < rhs->id ); } // --- general particle properties /// Unique identifier (in a Event object context) int id() const { return id_; } //void setId( int id ) { id_ = id; } /// Set the particle unique identifier in an event void setId( int id ) { id_ = id; } /// Electric charge (given as a float number, for the quarks and bound states) float charge() const; /// Set the electric charge sign (+-1 for charged or 0 for neutral particles) void setChargeSign( int sign ) { charge_sign_ = sign; } /// Role in the considered process Role role() const { return role_; } /// Set the particle role in the process void setRole( const Role& role ) { role_ = role; } /** * Codes 1-10 correspond to currently existing partons/particles, and larger codes contain partons/particles which no longer exist, or other kinds of event information * \brief Particle status */ Status status() const { return status_; } /// Set the particle decay/stability status void setStatus( Status status ) { status_ = status; } /// Set the PDG identifier (along with the particle's electric charge) /// \param[in] pdg PDG identifier /// \param[in] ch Electric charge (0, 1, or -1) void setPdgId( pdgid_t pdg, short ch = 0 ); /// Set the PDG identifier (along with the particle's electric charge) /// \param[in] pdg_id PDG identifier (incl. electric charge in e) void setPdgId( short pdg_id ); /// Retrieve the objectified PDG identifier pdgid_t pdgId() const; /// Retrieve the integer value of the PDG identifier int integerPdgId() const; /// Particle's helicity float helicity() const { return helicity_; } /// Set the helicity of the particle void setHelicity( float heli ) { helicity_ = heli; } /// Particle mass in GeV/c\f$^2\f$ /// \return Particle's mass inline double mass() const { return mass_; }; /// Compute the particle mass /// \param[in] off_shell Allow the particle to be produced off-shell? /// \note This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority) void computeMass( bool off_shell = false ); /// Set the particle mass, in GeV/c\f$^2\f$ /// \param m Mass in GeV/c\f$^2\f$ /// \note This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority) void setMass( double m = -1. ); /// Particle squared mass, in GeV\f$^2\f$/c\f$^4\f$ inline double mass2() const { return mass_*mass_; }; /// Retrieve the momentum object associated with this particle inline Momentum& momentum() { return momentum_; } /// Retrieve the momentum object associated with this particle inline Momentum momentum() const { return momentum_; } /// Associate a momentum object to this particle void setMomentum( const Momentum& mom, bool offshell = false ); /** * \brief Set the 3- or 4-momentum associated to the particle * \param[in] px Momentum along the \f$x\f$-axis, in GeV/c * \param[in] py Momentum along the \f$y\f$-axis, in GeV/c * \param[in] pz Momentum along the \f$z\f$-axis, in GeV/c * \param[in] e Energy, in GeV */ void setMomentum( double px, double py, double pz, double e = -1. ); /// Set the 4-momentum associated to the particle /// \param[in] p 4-momentum inline void setMomentum( double p[4] ) { setMomentum( p[0], p[1], p[2], p[3] ); } /// Set the particle's energy /// \param[in] e Energy, in GeV void setEnergy( double e = -1. ); /// Get the particle's energy, in GeV double energy() const; /// Get the particle's squared energy, in GeV\f$^2\f$ inline double energy2() const { return energy()*energy(); }; /// Is this particle a valid particle which can be used for kinematic computations? bool valid(); // --- particle relations /// Is this particle a primary particle? inline bool primary() const { return mothers_.empty(); } /// Set the mother particle /// \param[in] part A Particle object containing all the information on the mother particle void addMother( Particle& part ); /// Get the unique identifier to the mother particle from which this particle arises /// \return An integer representing the unique identifier to the mother of this particle in the event inline ParticlesIds mothers() const { return mothers_; } /** * \brief Add a decay product * \param[in] part The Particle object in which this particle will desintegrate or convert * \return A boolean stating if the particle has been added to the daughters list or if it was already present before */ void addDaughter( Particle& part ); /// Gets the number of daughter particles inline unsigned int numDaughters() const { return daughters_.size(); }; /// Get an identifiers list all daughter particles /// \return An integer vector containing all the daughters' unique identifier in the event inline ParticlesIds daughters() const { return daughters_; } // --- global particle information extraction /// Dump all the information on this particle into the standard output stream void dump() const; - private: + protected: /// Unique identifier in an event int id_; /// Electric charge (+-1 or 0) short charge_sign_; /// Momentum properties handler Momentum momentum_; /// Mass, in GeV/c\f$^2\f$ double mass_; /// Helicity float helicity_; /// Role in the process Role role_; /// Decay/stability status Status status_; /// List of mother particles ParticlesIds mothers_; /// List of daughter particles ParticlesIds daughters_; /// PDG id pdgid_t pdg_id_; /// Collection of standard, bare-level physical properties ParticleProperties phys_prop_; }; /// Compute the centre of mass energy of two particles (incoming or outgoing states) double CMEnergy( const Particle& p1, const Particle& p2 ); /// Compute the centre of mass energy of two particles (incoming or outgoing states) double CMEnergy( const Particle::Momentum& m1, const Particle::Momentum& m2 ); //bool operator<( const Particle& a, const Particle& b ) { return a.id Particles; /// List of particles' roles typedef std::vector ParticleRoles; /// Map between a particle's role and its associated Particle object typedef std::unordered_map > ParticlesMap; } #endif diff --git a/CepGen/IO/HepMCEventInterface.cpp b/CepGen/IO/HepMCEventInterface.cpp new file mode 100644 index 0000000..0086952 --- /dev/null +++ b/CepGen/IO/HepMCEventInterface.cpp @@ -0,0 +1,111 @@ +#include "CepGen/IO/HepMCEventInterface.h" + +#include "CepGen/Parameters.h" +#include "CepGen/Core/Exception.h" + +#include "CepGen/Event/Event.h" +#include "CepGen/Physics/Constants.h" + +#ifdef HEPMC3 +# include "HepMC3/Version.h" +# include "HepMC3/FourVector.h" +# include "HepMC3/GenEvent.h" +# include "HepMC3/GenRunInfo.h" +# include "HepMC3/GenVertex.h" +# include "HepMC3/GenParticle.h" +#else +# include "HepMC/Version.h" +# if !defined( HEPMC_VERSION_CODE ) // HepMC v2 +# include "HepMC/SimpleVector.h" +# include "HepMC/GenEvent.h" +# include "HepMC/GenVertex.h" +# include "HepMC/GenParticle.h" +# else +# include "HepMC/FourVector.h" +# include "HepMC/GenEvent.h" +# include "HepMC/GenRunInfo.h" +# include "HepMC/GenVertex.h" +# include "HepMC/GenParticle.h" +# define HEPMC3 +# endif +#endif + +namespace HepMC +{ + CepGenEvent::CepGenEvent() : + GenEvent( Units::GEV, Units::MM ) + { +#ifdef HEPMC3 + GenEvent::add_attribute( "AlphaQCD", make_shared( cepgen::constants::ALPHA_QCD ) ); + GenEvent::add_attribute( "AlphaEM", make_shared( cepgen::constants::ALPHA_EM ) ); +#else + GenEvent::set_alphaQCD( cepgen::constants::ALPHA_QCD ); + GenEvent::set_alphaQED( cepgen::constants::ALPHA_EM ); +#endif + } + + void + CepGenEvent::feedEvent( const cepgen::Event& evt ) + { + GenEvent::clear(); + GenEvent::weights().push_back( 1. ); // unweighted events + + // filling the particles content + const FourVector origin( 0., 0., 0., 0. ); + auto& part_vec = evt.particles(); + + int cm_id = 0, idx = 1; + +#ifdef HEPMC3 + GenVertexPtr v1 = make_shared( origin ), + v2 = make_shared( origin ), + vcm = make_shared( origin ); +#else + GenVertex* v1 = new GenVertex( origin ), + *v2 = new GenVertex( origin ), + *vcm = new GenVertex( origin ); +#endif + for ( unsigned int i = 0; i < part_vec.size(); ++i ) { + const auto& part_orig = part_vec.at( i ); + FourVector pmom( part_orig.momentum().px(), + part_orig.momentum().py(), + part_orig.momentum().pz(), + part_orig.energy() ); +#ifdef HEPMC3 + auto part = make_shared( pmom, part_orig.integerPdgId(), (int)part_orig.status() ); +#else + auto part = new GenParticle( pmom, part_orig.integerPdgId(), (int)part_orig.status() ); + part->suggest_barcode( idx++ ); +#endif + const auto& moth = part_orig.mothers(); + + switch ( part_orig.role() ) { + case cepgen::Particle::IncomingBeam1: { v1->add_particle_in( part ); } break; + case cepgen::Particle::IncomingBeam2: { v2->add_particle_in( part ); } break; + case cepgen::Particle::OutgoingBeam1: { v1->add_particle_out( part ); } break; + case cepgen::Particle::OutgoingBeam2: { v2->add_particle_out( part ); } break; + case cepgen::Particle::Parton1: { v1->add_particle_out( part ); vcm->add_particle_in( part ); } break; + case cepgen::Particle::Parton2: { v2->add_particle_out( part ); vcm->add_particle_in( part ); } break; + case cepgen::Particle::Intermediate: { cm_id = i; continue; } break; + case cepgen::Particle::CentralSystem: + default: { + if ( moth.empty() ) continue; + if ( *moth.begin() == cm_id ) vcm->add_particle_out( part ); + else + throw CG_FATAL( "HepMCHandler:fillEvent" ) + << "Other particle requested! Not yet implemented!"; + } break; + } + idx++; + } + GenEvent::add_vertex( v1 ); + GenEvent::add_vertex( v2 ); + GenEvent::add_vertex( vcm ); + +#ifndef HEPMC3 + GenEvent::set_beam_particles( *v1->particles_in_const_begin(), *v2->particles_in_const_begin() ); + GenEvent::set_signal_process_vertex( vcm ); +#endif + } +} + diff --git a/CepGen/IO/HepMCEventInterface.h b/CepGen/IO/HepMCEventInterface.h new file mode 100644 index 0000000..36018eb --- /dev/null +++ b/CepGen/IO/HepMCEventInterface.h @@ -0,0 +1,51 @@ +#ifndef CepGen_IO_HepMCEventInterface_h +#define CepGen_IO_HepMCEventInterface_h + +namespace cepgen { + class Parameters; + class Event; + class Particle; +} + +#ifdef HEPMC3 +# include "HepMC3/Version.h" +# include "HepMC3/GenEvent.h" +# define HepMC HepMC3 +#else +# include "HepMC/Version.h" +# include "HepMC/GenEvent.h" +#endif + +namespace HepMC +{ + /// Interfacing between CepGen and HepMC event definitions + /// \author Laurent Forthomme + /// \date Jul 2019 + class CepGenEvent : public GenEvent + { + public: + explicit CepGenEvent(); + /// Initialise this conversion object with CepGen parameters + void initialise( const cepgen::Parameters& ); + /// Feed a new CepGen event to this conversion object + /// \param[in] ev CepGen event to be fed + void feedEvent( const cepgen::Event& ev ); + /// Set the cross section for a given process + /// \param[in] id Process identifier + /// \param[in] xsec Process cross section, in pb + /// \param[in] xsec_err Uncertainty on process cross section, in pb + void setCrossSection( int id, double xsec, double xsec_err ); + /// Specify new process attributes + /// \param[in] id Process identifier + /// \param[in] xsec Process cross section, in pb + /// \param[in] q2_scale Hard event scale \f$Q^2\f$, in GeV\f$^2\f$ + /// \param[in] alpha_qed \f$\alpha_{\rm em}\f$ for this process + /// \param[in] alpha_qcd \f$\alpha_{\rm s}\f$ for this process + void setProcess( int id, double xsec, double q2_scale, double alpha_qed, double alpha_qcd ); + + private: + const cepgen::Parameters* params_; // borrowed + }; +} +#endif + diff --git a/CepGen/IO/HepMCHandler.cpp b/CepGen/IO/HepMCHandler.cpp index 300403b..1badab6 100644 --- a/CepGen/IO/HepMCHandler.cpp +++ b/CepGen/IO/HepMCHandler.cpp @@ -1,239 +1,148 @@ #include "CepGen/IO/ExportHandler.h" +#include "CepGen/IO/HepMCEventInterface.h" #include "CepGen/Parameters.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Event/Event.h" #include "CepGen/Physics/Constants.h" #ifdef HEPMC3 -# include "HepMC3/Version.h" # define HEPMC_VERSION HEPMC3_VERSION # include "HepMC3/WriterAscii.h" # include "HepMC3/WriterAsciiHepMC2.h" # include "HepMC3/WriterHEPEVT.h" # ifdef HEPMC3_ROOTIO # include "HepMC3/WriterRoot.h" # include "HepMC3/WriterRootTree.h" # endif -# include "HepMC3/FourVector.h" -# include "HepMC3/GenEvent.h" -# include "HepMC3/GenRunInfo.h" -# include "HepMC3/GenVertex.h" -# include "HepMC3/GenParticle.h" -using namespace HepMC3; +#elif !defined( HEPMC_VERSION_CODE ) // HepMC v2 +# include "HepMC/IO_GenEvent.h" #else -# include "HepMC/Version.h" -# if !defined( HEPMC_VERSION_CODE ) // HepMC v2 -# include "HepMC/IO_GenEvent.h" -# include "HepMC/SimpleVector.h" -# include "HepMC/GenEvent.h" -# include "HepMC/GenVertex.h" -# include "HepMC/GenParticle.h" -# else -# include "HepMC/WriterAscii.h" -# include "HepMC/WriterHEPEVT.h" -# include "HepMC/FourVector.h" -# include "HepMC/GenEvent.h" -# include "HepMC/GenRunInfo.h" -# include "HepMC/GenVertex.h" -# include "HepMC/GenParticle.h" -# define HEPMC3 -# endif -using namespace HepMC; +# include "HepMC/WriterAscii.h" +# include "HepMC/WriterHEPEVT.h" +# define HEPMC3 #endif #include +using namespace HepMC; + namespace cepgen { namespace io { /// Handler for the HepMC file output /// \tparam T HepMC writer handler (format-dependent) /// \author Laurent Forthomme /// \date Sep 2016 template class HepMCHandler : public GenericExportHandler { public: /// Class constructor explicit HepMCHandler( const ParametersList& ); - void initialise( const Parameters& /*params*/ ) override; + void initialise( const Parameters& /*params*/ ) override {} /// Writer operator void operator<<( const Event& ) override; void setCrossSection( double, double ) override; private: /// Clear the associated HepMC event content void clearEvent(); /// Populate the associated HepMC event with a Event object void fillEvent( const Event& ); /// Writer object std::unique_ptr output_; /// Generator cross section and error std::shared_ptr xs_; #ifdef HEPMC3 /// Auxiliary information on run std::shared_ptr runinfo_; #endif /// Associated HepMC event - std::shared_ptr event_; + std::shared_ptr event_; }; template HepMCHandler::HepMCHandler( const ParametersList& params ) : GenericExportHandler( "hepmc" ), output_( new T( params.get( "filename", "output.hepmc" ) ) ), xs_( new GenCrossSection ), #ifdef HEPMC3 runinfo_( new GenRunInfo ), #endif - event_( new GenEvent( Units::GEV, Units::MM ) ) + event_( new CepGenEvent ) { #ifdef HEPMC3 output_->set_run_info( runinfo_ ); runinfo_->set_weight_names( { "Default" } ); #endif CG_INFO( "HepMC" ) << "Interfacing module initialised " << "for HepMC version " << HEPMC_VERSION << "."; } template void HepMCHandler::operator<<( const Event& evt ) { fillEvent( evt ); #ifdef HEPMC3 output_->write_event( *event_ ); #else output_->write_event( event_.get() ); #endif } template void - HepMCHandler::initialise( const Parameters& ) - { -#ifdef HEPMC3 - event_->add_attribute( "AlphaQCD", make_shared( constants::ALPHA_QCD ) ); - event_->add_attribute( "AlphaEM", make_shared( constants::ALPHA_EM ) ); -#else - event_->set_alphaQCD( constants::ALPHA_QCD ); - event_->set_alphaQED( constants::ALPHA_EM ); -#endif - } - - template void HepMCHandler::setCrossSection( double xsect, double xsect_err ) { xs_->set_cross_section( xsect, xsect_err ); } template void HepMCHandler::fillEvent( const Event& evt ) { - event_->clear(); - + event_->feedEvent( evt ); // general information #ifdef HEPMC3 event_->set_cross_section( xs_ ); event_->set_run_info( runinfo_ ); #else event_->set_cross_section( *xs_ ); #endif - - event_->set_event_number( event_num_ ); - event_->weights().push_back( 1. ); // unweighted events - - // filling the particles content - const FourVector origin( 0., 0., 0., 0. ); - Particles part_vec = evt.particles(); - - int cm_id = 0, idx = 1; - -#ifdef HEPMC3 - GenVertexPtr v1 = make_shared( origin ), - v2 = make_shared( origin ), - vcm = make_shared( origin ); -#else - GenVertex* v1 = new GenVertex( origin ), - *v2 = new GenVertex( origin ), - *vcm = new GenVertex( origin ); -#endif - for ( unsigned int i = 0; i < part_vec.size(); ++i ) { - const Particle part_orig = part_vec.at( i ); - FourVector pmom( part_orig.momentum().px(), - part_orig.momentum().py(), - part_orig.momentum().pz(), - part_orig.energy() ); -#ifdef HEPMC3 - auto part = make_shared( pmom, part_orig.integerPdgId(), (int)part_orig.status() ); -#else - auto part = new GenParticle( pmom, part_orig.integerPdgId(), (int)part_orig.status() ); - part->suggest_barcode( idx++ ); -#endif - const ParticlesIds moth = part_orig.mothers(); - - switch ( part_orig.role() ) { - case Particle::IncomingBeam1: { v1->add_particle_in( part ); } break; - case Particle::IncomingBeam2: { v2->add_particle_in( part ); } break; - case Particle::OutgoingBeam1: { v1->add_particle_out( part ); } break; - case Particle::OutgoingBeam2: { v2->add_particle_out( part ); } break; - case Particle::Parton1: { v1->add_particle_out( part ); vcm->add_particle_in( part ); } break; - case Particle::Parton2: { v2->add_particle_out( part ); vcm->add_particle_in( part ); } break; - case Particle::Intermediate: { cm_id = i; continue; } break; - case Particle::CentralSystem: - default: { - if ( moth.empty() ) continue; - if ( *moth.begin() == cm_id ) vcm->add_particle_out( part ); - else - throw CG_FATAL( "HepMCHandler:fillEvent" ) - << "Other particle requested! Not yet implemented!"; - } break; - } - idx++; - } - event_->add_vertex( v1 ); - event_->add_vertex( v2 ); - event_->add_vertex( vcm ); - -#ifndef HEPMC3 - event_->set_beam_particles( *v1->particles_in_const_begin(), *v2->particles_in_const_begin() ); - event_->set_signal_process_vertex( *v1->vertices_begin() ); - event_->set_beam_particles( *v1->particles_in_const_begin(), *v2->particles_in_const_end() ); -#endif - event_num_++; + event_->set_event_number( event_num_++ ); } #ifdef HEPMC3 # if HEPMC_VERSION_CODE >= 3001000 typedef HepMCHandler HepMC2Handler; # endif typedef HepMCHandler HepMC3Handler; typedef HepMCHandler HEPEVTHandler; # ifdef HEPMC3_ROOTIO typedef HepMCHandler RootHandler; typedef HepMCHandler RootTreeHandler; # endif #else typedef HepMCHandler HepMC2Handler; #endif } } #ifdef HEPMC3 REGISTER_IO_MODULE( hepmc3, HepMC3Handler ) REGISTER_IO_MODULE( hepmc, HepMC3Handler ) REGISTER_IO_MODULE( hepevt, HEPEVTHandler ) # ifdef HEPMC3_ROOTIO REGISTER_IO_MODULE( hepmc_root, RootHandler ) REGISTER_IO_MODULE( hepmc_root_tree, RootTreeHandler ) # endif #else REGISTER_IO_MODULE( hepmc, HepMC2Handler ) #endif #if defined( HEPMC3 ) && HEPMC_VERSION_CODE >= 3001000 REGISTER_IO_MODULE( hepmc2, HepMC2Handler ) #endif