HepMC3Interface.cc
No OneTemporary
Actions

Size

6 KB

Referenced Files

None

Subscribers

None

HepMC3Interface.cc
View Options

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

	#include "HEJ/ConfigFlags.hh"
	#include "HEJ/exceptions.hh"

	// include before HepMC3 to override include of "HepMC3/LHEF.h"
	// (theoretically both files should be the same)
	#include "LHEF/LHEF.h"

	#ifdef HEJ_BUILD_WITH_HepMC3

	#include <cassert>
	#include <cmath>
	#include <string>
	#include <utility>
	#include <vector>

	#include "HepMC3/Attribute.h"
	#include "HepMC3/FourVector.h"
	#include "HepMC3/GenCrossSection.h"
	#include "HepMC3/GenCrossSection_fwd.h"
	#include "HepMC3/GenEvent.h"
	#include "HepMC3/GenParticle.h"
	#include "HepMC3/GenParticle_fwd.h"
	#include "HepMC3/GenRunInfo.h"
	#include "HepMC3/GenVertex.h"
	#include "HepMC3/LHEFAttributes.h"
	#include "HepMC3/Units.h"

	#include "HEJ/Event.hh"
	#include "HEJ/Parameters.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/detail/HepMCInterface_common.hh"

	#else // no HepMC3

	#include "HEJ/utility.hh"

	#endif

	#ifdef HEJ_BUILD_WITH_HepMC3

	namespace HEJ {

	namespace detail_HepMC {
	template<>
	struct HepMCVersion<3> {
	using GenEvent = HepMC3::GenEvent;
	using Beam = std::array<HepMC3::GenParticlePtr,2>;
	};

	template<>
	auto make_particle_ptr<3> (
	Particle const & sp, int status
	) {
	return HepMC3::make_shared<HepMC3::GenParticle>(
	to_FourVector<HepMC3::FourVector>(sp),
	static_cast<int> (sp.type),
	status
	);
	}

	template<>
	auto make_vx_ptr<3>() {
	return HepMC3::make_shared<HepMC3::GenVertex>();
	}
	} // namespace detail_HepMC

	namespace {
	void reset_weight_info(LHEF::HEPRUP & heprup){
	heprup.IDWTUP = 2;
	// use placeholders for unknown init block values
	// we can overwrite them after processing all events
	heprup.XSECUP = {0.};
	heprup.XERRUP = {0.};
	heprup.XMAXUP = {0.};
	}

	HepMC3::shared_ptr<HepMC3::GenRunInfo> init_runinfo(LHEF::HEPRUP heprup){
	reset_weight_info(heprup);
	auto runinfo{ HepMC3::make_shared<HepMC3::GenRunInfo>() };

	auto hepr{ HepMC3::make_shared<HepMC3::HEPRUPAttribute>() };
	hepr->heprup = std::move(heprup);
	runinfo->add_attribute(std::string("HEPRUP"), hepr);
	for(auto const & gen: hepr->heprup.generators){
	runinfo->tools().emplace_back(
	HepMC3::GenRunInfo::ToolInfo{gen.name, gen.version, gen.contents} );
	}
	return runinfo;
	}

	std::vector<std::string> get_weight_names(Event const & ev){
	std::vector<std::string> names;
	names.reserve(ev.variations().size()+1); // +1 from central
	names.emplace_back(""); // rivet assumes central band to have no name
	for(auto const & var : ev.variations()){
	if(var.description){
	names.emplace_back( to_simple_string(*var.description) );
	} else {
	names.emplace_back( "" );
	}
	}
	assert(names.size() == ev.variations().size()+1);
	return names;
	}
	} // namespace

	HepMC3Interface::HepMC3Interface(LHEF::HEPRUP heprup):
	beam_particle_{static_cast<ParticleID>(heprup.IDBMUP.first),
	static_cast<ParticleID>(heprup.IDBMUP.second)},
	beam_energy_{heprup.EBMUP.first, heprup.EBMUP.second},
	run_info_{ init_runinfo(std::move(heprup)) },
	event_count_(0.), tot_weight_(0.), tot_weight2_(0.),
	xs_{std::make_shared<HepMC3::GenCrossSection>()}
	{
	}

	HepMC3::GenEvent HepMC3Interface::init_event(Event const & event) const {

	const std::array<HepMC3::GenParticlePtr,2> beam {
	HepMC3::make_shared<HepMC3::GenParticle>(
	HepMC3::FourVector(0,0,-beam_energy_[0],beam_energy_[0]),
	beam_particle_[0], detail_HepMC::Status::beam ),
	HepMC3::make_shared<HepMC3::GenParticle>(
	HepMC3::FourVector(0,0, beam_energy_[1],beam_energy_[1]),
	beam_particle_[1], detail_HepMC::Status::beam )
	};
	auto hepmc_ev{ detail_HepMC::HepMC_init_kinematics<3>(
	event, beam, HepMC3::GenEvent{ HepMC3::Units::GEV, HepMC3::Units::MM }
	) };
	// set up run specific informations
	if( run_info_->weight_names().size() != event.variations().size()+1 ){
	run_info_->set_weight_names( get_weight_names(event) );
	}
	// order matters: weights in hepmc_ev initialised when registering run_info
	hepmc_ev.set_run_info(run_info_);
	assert(hepmc_ev.weights().size() == event.variations().size()+1);
	for(size_t i=0; i<event.variations().size(); ++i){
	hepmc_ev.weights()[i+1] = event.variations()[i].weight;
	//! @TODO set variation specific cross section
	//! the problem is that set_cross_section overwrites everything
	}
	return hepmc_ev;
	}

	void HepMC3Interface::set_central(
	HepMC3::GenEvent & out_ev, Event const & event, int const weight_index
	){
	EventParameters event_param;
	if(weight_index < 0)
	event_param = event.central();
	else if ( static_cast<size_t>(weight_index) < event.variations().size())
	event_param = event.variations(weight_index);
	else
	throw std::invalid_argument{
	"HepMC3Interface tried to access a weight outside of the variation range."
	};
	const double wt = event_param.weight;
	tot_weight_ += wt;
	tot_weight2_ += wt * wt;
	++event_count_;

	// central always on first
	assert(out_ev.weights().size() == event.variations().size()+1);
	out_ev.weights()[0] = wt;

	// out_ev can be setup with a different central scale -> save xs manually
	out_ev.set_cross_section(xs_);
	assert(out_ev.cross_section() && out_ev.cross_section() == xs_);
	// overwrites all previously set xs ...
	xs_->set_cross_section(tot_weight_,std::sqrt(tot_weight2_));

	out_ev.set_event_number(event_count_);
	/// @TODO add number of attempted events
	xs_->set_accepted_events(event_count_);

	/// @TODO add alphaQCD (need function) and alphaQED
	/// @TODO output pdf (currently not avaiable from event alone)
	}

	HepMC3::GenEvent HepMC3Interface::operator()(
	Event const & event, int const weight_index
	){
	HepMC3::GenEvent out_ev(init_event(event));
	set_central(out_ev, event, weight_index);
	return out_ev;
	}

	} // namespace HEJ

	#else // no HepMC3 => empty class

	namespace HepMC3 {
	class GenEvent {};
	class GenCrossSection {};
	class GenRunInfo {};
	}
	namespace HEJ {
	HepMC3Interface::HepMC3Interface(LHEF::HEPRUP /heprup/){
	ignore(beam_particle_,beam_energy_,event_count_,tot_weight_,tot_weight2_);
	throw std::invalid_argument(
	"Failed to create HepMC3Interface: "
	"HEJ 2 was built without HepMC3 support"
	);
	}

	HepMC3::GenEvent HepMC3Interface::operator()(
	Event const & /event/, int /weight_index/
	){return HepMC3::GenEvent();}
	HepMC3::GenEvent HepMC3Interface::init_event(Event const & /event/) const
	{return HepMC3::GenEvent();}
	void HepMC3Interface::set_central(
	HepMC3::GenEvent & /out_ev/, Event const & /event/, int /weight_index/
	){}
	}
	#endif

	namespace HEJ {
	HepMC3Interface::~HepMC3Interface() = default;
	}

File Metadata

Mime Type: text/x-c
Expires: Tue, Sep 30, 5:44 AM (5 h, 52 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6458687
Default Alt Text: HepMC3Interface.cc (6 KB)

HepMC3Interface.ccNo OneTemporaryActions

HepMC3Interface.ccView Options

File Metadata

Event Timeline

HepMC3Interface.cc
No OneTemporary
Actions

HepMC3Interface.cc
View Options