diff --git a/src/HepMC2Interface.cc b/src/HepMC2Interface.cc index ac0cfc0..dbc675b 100644 --- a/src/HepMC2Interface.cc +++ b/src/HepMC2Interface.cc @@ -1,160 +1,160 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include "HEJ/HepMC2Interface.hh" #include "HEJ/exceptions.hh" #ifdef HEJ_BUILD_WITH_HepMC2 #include <math.h> #include <utility> #include "HEJ/detail/HepMCInterface_common.hh" #include "HEJ/Event.hh" #include "HEJ/Particle.hh" #include "LHEF/LHEF.h" #include "HepMC/GenCrossSection.h" #include "HepMC/GenEvent.h" #include "HepMC/GenParticle.h" #include "HepMC/GenVertex.h" namespace HEJ{ namespace detail_HepMC { template<> struct HepMCVersion<2> { using GenEvent = HepMC::GenEvent; using Beam = std::array<HepMC::GenParticle*,2>; }; template<> auto make_particle_ptr<2> ( Particle const & sp, int status ) { return new HepMC::GenParticle( to_FourVector<HepMC::FourVector>(sp), static_cast<int> (sp.type), status ); } template<> auto make_vx_ptr<2>() { return new HepMC::GenVertex(); } } HepMC2Interface::HepMC2Interface(LHEF::HEPRUP const & heprup): event_count_(0.), tot_weight_(0.), tot_weight2_(0.) { beam_particle_[0] = static_cast<ParticleID>(heprup.IDBMUP.first); beam_particle_[1] = static_cast<ParticleID>(heprup.IDBMUP.second); beam_energy_[0] = heprup.EBMUP.first; beam_energy_[1] = heprup.EBMUP.second; } HepMC::GenCrossSection HepMC2Interface::cross_section() const { HepMC::GenCrossSection xs; xs.set_cross_section(tot_weight_, sqrt(tot_weight2_)); return xs; } HepMC::GenEvent HepMC2Interface::init_kinematics(Event const & event){ const std::array<HepMC::GenParticle*,2> beam { new HepMC::GenParticle( HepMC::FourVector(0,0,-beam_energy_[0],beam_energy_[0]), beam_particle_[0], detail_HepMC::status_beam ), new HepMC::GenParticle( HepMC::FourVector(0,0, beam_energy_[1],beam_energy_[1]), beam_particle_[1], detail_HepMC::status_beam ) }; return detail_HepMC::HepMC_init_kinematics<2>( event, beam, HepMC::GenEvent{ HepMC::Units::GEV, HepMC::Units::MM } ); } void HepMC2Interface::set_central(HepMC::GenEvent & out_ev, Event const & event, ssize_t const weight_index ){ EventParameters event_param; if(weight_index < 0) event_param = event.central(); else if ( (size_t) weight_index < event.variations().size()) event_param = event.variations(weight_index); else throw std::invalid_argument{ "HepMC2Interface tried to access a weight outside of the variation range." }; const double wt = event_param.weight; tot_weight_ += wt; tot_weight2_ += wt * wt; if(out_ev.weights().size() == 0){ out_ev.weights().push_back(wt); } else { // central always on first out_ev.weights()[0] = wt; } out_ev.set_cross_section( cross_section() ); - out_ev.set_signal_process_id(event.type()+1); // "+1": conistent with lhe + out_ev.set_signal_process_id(event.type()); out_ev.set_event_scale(event_param.mur); ++event_count_; out_ev.set_event_number(event_count_); /// @TODO add alphaQCD (need function) and alphaQED /// @TODO output pdf (currently not avaiable from event alone) } void HepMC2Interface::add_variation(HepMC::GenEvent & out_ev, std::vector<EventParameters> const & varis ){ for(auto const & var: varis){ out_ev.weights().push_back(var.weight); } /// @TODO add name list for weights } HepMC::GenEvent HepMC2Interface::operator()(Event const & event, ssize_t const weight_index ){ HepMC::GenEvent out_ev(init_kinematics(event)); set_central(out_ev, event, weight_index); add_variation(out_ev, event.variations()); return out_ev; } } #else // no HepMC2 => empty class namespace HepMC { class GenEvent {}; class GenCrossSection {}; } namespace HEJ{ HepMC2Interface::HepMC2Interface(LHEF::HEPRUP const &){ throw std::invalid_argument( "Failed to create HepMC2Interface: " "HEJ 2 was built without HepMC2 support" ); } HepMC::GenEvent HepMC2Interface::operator()(Event const &, ssize_t) {return HepMC::GenEvent();} HepMC::GenEvent HepMC2Interface::init_kinematics(Event const &) {return HepMC::GenEvent();} void HepMC2Interface::add_variation(HepMC::GenEvent &, std::vector<EventParameters> const &){} void HepMC2Interface::set_central(HepMC::GenEvent &, Event const &, ssize_t){} HepMC::GenCrossSection HepMC2Interface::cross_section() const {return HepMC::GenCrossSection();} } #endif