HepMCInterface.cc
No OneTemporary
Actions

Size

6 KB

Referenced Files

None

Subscribers

None

HepMCInterface.cc
View Options

	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/HepMCInterface.hh"

	#include "HEJ/exceptions.hh"

	#ifdef HEJ_BUILD_WITH_HepMC_VERSION

	#include <math.h>
	#include <utility>

	#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 {
	HepMC::FourVector to_FourVector(Particle const & sp){
	return {sp.px(), sp.py(), sp.pz(), sp.E()};
	}
	constexpr int status_beam = 4;
	constexpr int status_in = 11;
	constexpr int status_decayed = 2;
	constexpr int status_out = 1;

	template<class HepMCClass, typename... Args>
	auto make_ptr(Args&&... args){
	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	return HepMC::make_shared<HepMCClass>(std::forward<Args>(args)...);
	#else
	return new HepMCClass(std::forward<Args>(args)...);
	#endif
	}
	} // namespace anonymous

	HepMCInterface::HepMCInterface(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 HepMCInterface::cross_section() const {
	HepMC::GenCrossSection xs;
	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	xs.set_cross_section(tot_weight_, sqrt(tot_weight2_), event_count_);
	/// @TODO add number of attempted events
	#else // HepMC 2
	xs.set_cross_section(tot_weight_, sqrt(tot_weight2_));
	#endif
	return xs;
	}

	HepMC::GenEvent HepMCInterface::init_kinematics(Event const & event) {
	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	using GenParticlePtr = HepMC::GenParticlePtr;
	#else
	using GenParticlePtr = HepMC::GenParticle*;
	#endif
	HepMC::GenEvent out_ev{HepMC::Units::GEV, HepMC::Units::MM};

	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	static
	#endif
	const std::array<GenParticlePtr,2> beam {
	make_ptr<HepMC::GenParticle>(
	HepMC::FourVector(0,0,-beam_energy_[0],beam_energy_[0]),
	beam_particle_[0],status_beam ),
	make_ptr<HepMC::GenParticle>(
	HepMC::FourVector(0,0, beam_energy_[1],beam_energy_[1]),
	beam_particle_[1],status_beam )
	};
	out_ev.set_beam_particles(beam[0],beam[1]);

	auto vx = make_ptr<HepMC::GenVertex>();
	for(size_t i=0; i<event.incoming().size(); ++i){
	auto particle{ make_ptr<HepMC::GenParticle>(
	to_FourVector(event.incoming()[i]),
	static_cast<int>(event.incoming()[i].type), status_in
	) };
	auto vx_beam{ make_ptr<HepMC::GenVertex>() };
	vx_beam->add_particle_in(beam[i]);
	vx_beam->add_particle_out(particle);
	out_ev.add_vertex(vx_beam);
	vx->add_particle_in(particle);
	}
	for(size_t i=0; i < event.outgoing().size(); ++i){
	auto const & out = event.outgoing()[i];
	auto particle = make_ptr<HepMC::GenParticle>(
	to_FourVector(out), static_cast<int>(out.type), status_out
	);
	const int status = event.decays().count(i)?status_decayed:status_out;
	particle->set_status(status);
	if( status == status_decayed){
	auto vx_decay = make_ptr<HepMC::GenVertex>();
	vx_decay->add_particle_in(particle);
	for( auto const & out: event.decays().at(i)){
	vx_decay->add_particle_out(
	make_ptr<HepMC::GenParticle>(
	to_FourVector(out), static_cast<int>(out.type), status_out
	)
	);
	}
	out_ev.add_vertex(vx_decay);
	}
	vx->add_particle_out(particle);
	}
	out_ev.add_vertex(vx);

	return out_ev;
	}

	void HepMCInterface::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{
	"HepMCInterface 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;
	}

	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	out_ev.set_cross_section(
	HepMC::make_shared<HepMC::GenCrossSection>(cross_section()) );
	#else // HepMC 2
	out_ev.set_cross_section( cross_section() );
	out_ev.set_signal_process_id(event.type()+1); // "+1": conistent with lhe
	out_ev.set_event_scale(event_param.mur);
	#endif

	++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 HepMCInterface::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 HepMCInterface::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 HepMC => empty class
	namespace HepMC {
	class GenEvent {};
	class GenCrossSection {};
	}
	namespace HEJ{
	HepMCInterface::HepMCInterface(LHEF::HEPRUP const &){
	throw std::invalid_argument(
	"Failed to create HepMCInterface: "
	"HEJ 2 was built without HepMC support"
	);
	}

	HepMC::GenEvent HepMCInterface::operator()(Event const &, ssize_t)
	{return HepMC::GenEvent();}
	HepMC::GenEvent HepMCInterface::init_kinematics(Event const &)
	{return HepMC::GenEvent();}
	void HepMCInterface::add_variation(HepMC::GenEvent &,
	std::vector<EventParameters> const &){}
	void HepMCInterface::set_central(HepMC::GenEvent &, Event const &, ssize_t) {}
	HepMC::GenCrossSection HepMCInterface::cross_section() const
	{return HepMC::GenCrossSection();}
	}
	#endif

File Metadata

Mime Type: text/x-c
Expires: Tue, Sep 30, 6:08 AM (1 d, 10 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6561112
Default Alt Text: HepMCInterface.cc (6 KB)

HepMCInterface.ccNo OneTemporaryActions

HepMCInterface.ccView Options

File Metadata

Event Timeline

HepMCInterface.cc
No OneTemporary
Actions

HepMCInterface.cc
View Options