diff --git a/CepGen/IO/HepMCEventInterface.cpp b/CepGen/IO/HepMCEventInterface.cpp index 0086952..64b2cbd 100644 --- a/CepGen/IO/HepMCEventInterface.cpp +++ b/CepGen/IO/HepMCEventInterface.cpp @@ -1,111 +1,128 @@ #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 ); + GenVertexPtr v1 = make_shared( origin ), v2 = make_shared( origin ); + GenVertexPtr vcm = make_shared( origin ); #else - GenVertex* v1 = new GenVertex( origin ), - *v2 = new GenVertex( origin ), - *vcm = new GenVertex( origin ); + GenVertex* v1 = new GenVertex( origin ), *v2 = new GenVertex( origin ); + GenVertex* 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() ); + const auto& mom_orig = part_orig.momentum(); + FourVector pmom( mom_orig.px(), mom_orig.py(), mom_orig.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 ); + 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: + // skip the two-parton system and propagate the parentage + cm_id = i; + continue; + case cepgen::Particle::CentralSystem: default: { + const auto& moth = part_orig.mothers(); + if ( moth.empty() ) + // skip disconnected lines + continue; + short m1 = *moth.begin(), m2 = moth.size() > 1 ? *moth.rbegin() : -1; + if ( m1 == cm_id || ( m2 >= 0 && ( m1 < cm_id && m2 >= cm_id ) ) ) + // if connected to the central system + 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 index 36018eb..47845f7 100644 --- a/CepGen/IO/HepMCEventInterface.h +++ b/CepGen/IO/HepMCEventInterface.h @@ -1,51 +1,29 @@ #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 cepgen { class Event; } + 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 1badab6..fc446c5 100644 --- a/CepGen/IO/HepMCHandler.cpp +++ b/CepGen/IO/HepMCHandler.cpp @@ -1,148 +1,140 @@ #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 +# define HEPMC_VERSION_CODE HEPMC3_VERSION_CODE # include "HepMC3/WriterAscii.h" # include "HepMC3/WriterAsciiHepMC2.h" # include "HepMC3/WriterHEPEVT.h" # ifdef HEPMC3_ROOTIO # include "HepMC3/WriterRoot.h" # include "HepMC3/WriterRootTree.h" # endif -#elif !defined( HEPMC_VERSION_CODE ) // HepMC v2 -# include "HepMC/IO_GenEvent.h" #else -# include "HepMC/WriterAscii.h" -# include "HepMC/WriterHEPEVT.h" -# define HEPMC3 +# include "HepMC/Version.h" +# if !defined( HEPMC_VERSION_CODE ) // HepMC v2 +# include "HepMC/IO_GenEvent.h" +# else +# include "HepMC/WriterAscii.h" +# include "HepMC/WriterHEPEVT.h" +# define HEPMC3 +# endif #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 {} /// 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_; }; 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 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::setCrossSection( double xsect, double xsect_err ) { xs_->set_cross_section( xsect, xsect_err ); } template void HepMCHandler::fillEvent( const Event& evt ) { 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_++ ); } -#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 ) +REGISTER_IO_MODULE( hepmc, HepMCHandler ) +REGISTER_IO_MODULE( hepmc3, HepMCHandler ) +REGISTER_IO_MODULE( hepevt, HepMCHandler ) +# if HEPMC_VERSION_CODE >= 3001000 +REGISTER_IO_MODULE( hepmc2, HepMCHandler ) +# endif # ifdef HEPMC3_ROOTIO -REGISTER_IO_MODULE( hepmc_root, RootHandler ) -REGISTER_IO_MODULE( hepmc_root_tree, RootTreeHandler ) +REGISTER_IO_MODULE( hepmc_root, HepMCHandler ) +REGISTER_IO_MODULE( hepmc_root_tree, HepMCHandler ) # endif #else -REGISTER_IO_MODULE( hepmc, HepMC2Handler ) -#endif -#if defined( HEPMC3 ) && HEPMC_VERSION_CODE >= 3001000 -REGISTER_IO_MODULE( hepmc2, HepMC2Handler ) +REGISTER_IO_MODULE( hepmc, HepMCHandler ) #endif