Page MenuHomeHEPForge
Files F9501172

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
index 336f49b..1a82773 100644
--- a/include/HEJ/Event.hh
+++ b/include/HEJ/Event.hh
@@ -1,372 +1,372 @@
/** \file
* \brief Declares the Event class and helpers
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <cstddef>
#include <iosfwd>
#include <iterator>
#include <unordered_map>
#include <utility>
#include <vector>
#include "boost/iterator/filter_iterator.hpp"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/PseudoJet.hh"
-#include "HEJ/EWConstants.hh"
#include "HEJ/Parameters.hh"
#include "HEJ/Particle.hh"
#include "HEJ/event_types.hh"
namespace LHEF {
class HEPEUP;
class HEPRUP;
}
namespace fastjet {
class JetDefinition;
}
namespace HEJ {
+ struct EWConstants;
struct RNG;
struct UnclusteredEvent;
/** @brief An event with clustered jets
*
* This is the main HEJ 2 event class.
* It contains kinematic information including jet clustering,
* parameter (e.g. scale) settings and the event weight.
*/
class Event {
public:
class EventData;
//! Iterator over partons
using ConstPartonIterator = boost::filter_iterator<
bool (*)(Particle const &),
std::vector<Particle>::const_iterator
>;
//! Reverse Iterator over partons
using ConstReversePartonIterator = std::reverse_iterator<
ConstPartonIterator>;
//! No default Constructor
Event() = delete;
//! Event Constructor adding jet clustering to an unclustered event
//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0
[[deprecated("UnclusteredEvent will be replaced by EventData")]]
Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
);
//! @name Particle Access
//! @{
//! Incoming particles
std::array<Particle, 2> const & incoming() const{
return incoming_;
}
//! Outgoing particles
std::vector<Particle> const & outgoing() const{
return outgoing_;
}
//! Iterator to the first outgoing parton
ConstPartonIterator begin_partons() const;
//! Iterator to the first outgoing parton
ConstPartonIterator cbegin_partons() const;
//! Iterator to the end of the outgoing partons
ConstPartonIterator end_partons() const;
//! Iterator to the end of the outgoing partons
ConstPartonIterator cend_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator rbegin_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator crbegin_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator rend_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator crend_partons() const;
//! Particle decays
/**
* The key in the returned map corresponds to the index in the
* vector returned by outgoing()
*/
std::unordered_map<std::size_t, std::vector<Particle>> const & decays()
const {
return decays_;
}
//! The jets formed by the outgoing partons, sorted in rapidity
std::vector<fastjet::PseudoJet> const & jets() const{
return jets_;
}
//! @}
//! @name Weight variations
//! @{
//! All chosen parameter, i.e. scale choices (const version)
Parameters<EventParameters> const & parameters() const{
return parameters_;
}
//! All chosen parameter, i.e. scale choices
Parameters<EventParameters> & parameters(){
return parameters_;
}
//! Central parameter choice (const version)
EventParameters const & central() const{
return parameters_.central;
}
//! Central parameter choice
EventParameters & central(){
return parameters_.central;
}
//! Parameter (scale) variations (const version)
std::vector<EventParameters> const & variations() const{
return parameters_.variations;
}
//! Parameter (scale) variations
std::vector<EventParameters> & variations(){
return parameters_.variations;
}
//! Parameter (scale) variation (const version)
/**
* @param i Index of the requested variation
*/
EventParameters const & variations(std::size_t i) const{
return parameters_.variations.at(i);
}
//! Parameter (scale) variation
/**
* @param i Index of the requested variation
*/
EventParameters & variations(std::size_t i){
return parameters_.variations.at(i);
}
//! @}
//! Indices of the jets the outgoing partons belong to
/**
* @param jets Jets to be tested
* @returns A vector containing, for each outgoing parton,
* the index in the vector of jets the considered parton
* belongs to. If the parton is not inside any of the
* passed jets, the corresponding index is set to -1.
*/
std::vector<int> particle_jet_indices(
std::vector<fastjet::PseudoJet> const & jets
) const {
return cs_.particle_jet_indices(jets);
}
//! particle_jet_indices() of the Event jets()
std::vector<int> particle_jet_indices() const {
return particle_jet_indices(jets());
}
//! Jet definition used for clustering
fastjet::JetDefinition const & jet_def() const{
return cs_.jet_def();
}
//! Minimum jet transverse momentum
double min_jet_pt() const{
return min_jet_pt_;
}
//! Event type
event_type::EventType type() const{
return type_;
}
//! Give colours to each particle
/**
* @returns true if new colours are generated, i.e. same as is_resummable()
* @details Colour ordering is done according to leading colour in the MRK
* limit, see \cite Andersen:2011zd. This only affects \ref
* is_resummable() "HEJ" configurations, all other \ref event_type
* "EventTypes" will be ignored.
* @note This overwrites all previously set colours.
*/
bool generate_colours(RNG & /*ran*/);
//! Check that current colours are leading in the high energy limit
/**
* @details Checks that the colour configuration can be split up in
* multiple, rapidity ordered, non-overlapping ladders. Such
* configurations are leading in the MRK limit, see
* \cite Andersen:2011zd
*
* @note This is _not_ to be confused with \ref is_resummable(), however
* for all resummable states it is possible to create a leading colour
* configuration, see generate_colours()
*/
bool is_leading_colour() const;
/**
* @brief Check if given event could have been produced by HEJ
* @details A HEJ state has to fulfil:
* 1. type() has to be \ref is_resummable() "resummable"
* 2. Soft radiation in the tagging jets contributes at most to
* `max_ext_soft_pt_fraction` of the total jet \f$ p_\perp \f$
*
* @note This is true for any resummed stated produced by the
* EventReweighter or any \ref is_resummable() "resummable" Leading
* Order state.
*
* @param max_frac Maximum transverse momentum fraction from soft radiation
* in extremal jets
* @param min_pt Absolute minimal \f$ p_\perp \f$,
* \b deprecated use max_frac instead
* @return True if this state could have been produced by HEJ
*/
bool valid_hej_state(
double max_frac, double min_pt = 0.) const;
private:
//! \internal
//! @brief Construct Event explicitly from input.
/** This is only intended to be called from EventData.
*
* \warning The input is taken _as is_, sorting and classification has to be
* done externally, i.e. by EventData
*/
Event(
std::array<Particle, 2> && incoming,
std::vector<Particle> && outgoing,
std::unordered_map<std::size_t, std::vector<Particle>> && decays,
Parameters<EventParameters> && parameters,
fastjet::JetDefinition const & jet_def,
double min_jet_pt
);
//! Iterator over partons (non-const)
using PartonIterator = boost::filter_iterator<
bool (*)(Particle const &),
std::vector<Particle>::iterator
>;
//! Reverse Iterator over partons (non-const)
using ReversePartonIterator = std::reverse_iterator<PartonIterator>;
//! Iterator to the first outgoing parton (non-const)
PartonIterator begin_partons();
//! Iterator to the end of the outgoing partons (non-const)
PartonIterator end_partons();
//! Reverse Iterator to the first outgoing parton (non-const)
ReversePartonIterator rbegin_partons();
//! Reverse Iterator to the first outgoing parton (non-const)
ReversePartonIterator rend_partons();
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
std::unordered_map<std::size_t, std::vector<Particle>> decays_;
std::vector<fastjet::PseudoJet> jets_;
Parameters<EventParameters> parameters_;
fastjet::ClusterSequence cs_;
double min_jet_pt_;
event_type::EventType type_;
}; // end class Event
//! Class to store general Event setup, i.e. Phase space and weights
class Event::EventData {
public:
//! Default Constructor
EventData() = default;
//! Constructor from LesHouches event information
EventData(LHEF::HEPEUP const & hepeup);
//! Constructor with all values given
EventData(
std::array<Particle, 2> incoming,
std::vector<Particle> outgoing,
std::unordered_map<std::size_t, std::vector<Particle>> decays,
Parameters<EventParameters> parameters
):
incoming(std::move(incoming)), outgoing(std::move(outgoing)),
decays(std::move(decays)), parameters(std::move(parameters))
{}
//! Generate an Event from the stored EventData.
/**
* @details Do jet clustering and classification.
* Use this to generate an Event.
*
* @note Calling this function destroys EventData
*
* @param jet_def Jet definition
* @param min_jet_pt minimal \f$p_T\f$ for each jet
*
* @returns Full clustered and classified event.
*/
Event cluster(
fastjet::JetDefinition const & jet_def, double min_jet_pt);
//! Alias for cluster()
Event operator()(
fastjet::JetDefinition const & jet_def, double const min_jet_pt){
return cluster(jet_def, min_jet_pt);
}
//! Sort particles in rapidity
void sort();
//! Reconstruct intermediate particles from final-state leptons
/**
* Final-state leptons are created from virtual photons, W, or Z bosons.
* This function tries to reconstruct such intermediate bosons if they
* are not part of the event record.
*/
- void reconstruct_intermediate(EWConstants const & ew_parameters);
+ void reconstruct_intermediate(EWConstants const &);
//! Incoming particles
std::array<Particle, 2> incoming;
//! Outcoing particles
std::vector<Particle> outgoing;
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<std::size_t, std::vector<Particle>> decays;
//! Parameters, e.g. scale or inital weight
Parameters<EventParameters> parameters;
}; // end class EventData
//! Print Event
std::ostream& operator<<(std::ostream & os, Event const & ev);
//! Square of the partonic centre-of-mass energy \f$\hat{s}\f$
double shat(Event const & ev);
//! Convert an event to a LHEF::HEPEUP
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * /*heprup*/);
// put deprecated warning at the end, so don't get the warning inside Event.hh,
// additionally doxygen can not identify [[deprecated]] correctly
struct [[deprecated("UnclusteredEvent will be replaced by EventData")]]
UnclusteredEvent;
//! An event before jet clustering
//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0
struct UnclusteredEvent{
//! Default Constructor
UnclusteredEvent() = default;
//! Constructor from LesHouches event information
UnclusteredEvent(LHEF::HEPEUP const & hepeup);
std::array<Particle, 2> incoming; /**< Incoming Particles */
std::vector<Particle> outgoing; /**< Outgoing Particles */
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<std::size_t, std::vector<Particle>> decays;
//! Central parameter (e.g. scale) choice
EventParameters central;
std::vector<EventParameters> variations; /**< For parameter variation */
};
} // namespace HEJ
diff --git a/src/Event.cc b/src/Event.cc
index 44e8f79..57cd345 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,1303 +1,1304 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/Event.hh"
#include <algorithm>
#include <cassert>
#include <cstdlib>
#include <iomanip>
#include <iterator>
#include <memory>
#include <numeric>
#include <ostream>
#include <string>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/PseudoJet.hh"
#include "LHEF/LHEF.h"
#include "HEJ/Constants.hh"
+#include "HEJ/EWConstants.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/RNG.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/optional.hh"
namespace HEJ {
namespace {
using std::size_t;
//! LHE status codes
namespace lhe_status {
enum Status: int {
in = -1,
decay = 2,
out = 1,
};
}
using LHE_Status = lhe_status::Status;
//! true if leptonic W decay
bool valid_W_decay( int const w_type, // sign of W
std::vector<Particle> const & decays
){
if(decays.size() != 2) // no 1->2 decay
return false;
const int pidsum = decays[0].type + decays[1].type;
if( std::abs(pidsum) != 1 || pidsum != w_type ) // correct charge
return false;
// leptonic decay (only check first, second follows from pidsum)
if( w_type == 1 ) // W+
return is_antilepton(decays[0]) || is_neutrino(decays[0]);
// W-
return is_lepton(decays[0]) || is_antineutrino(decays[0]);
}
//! true for Z decay to charged leptons
bool valid_Z_decay(std::vector<Particle> const & decays){
if(decays.size() != 2) // no 1->2 decay
return false;
const int pidsum = decays[0].type + decays[1].type;
if( std::abs(pidsum) != 0 ) // correct charge
return false;
// leptonic decay (only check first, second follows from pidsum)
return is_anylepton(decays[0]) && !is_anyneutrino(decays[0]);
}
//! true if supported decay
bool valid_decay(std::vector<Particle> const & decays){
return valid_W_decay(+1, decays) || // Wp
valid_W_decay(-1, decays) || // Wm
valid_Z_decay( decays) // Z/gamma
;
}
/// @name helper functions to determine event type
//@{
/**
* \brief check if final state valid for HEJ
*
* check final state has the expected number of valid decays for bosons
* and all the rest are quarks or gluons
*/
bool final_state_ok(Event const & ev){
size_t invalid_decays = ev.decays().size();
std::vector<Particle> const & outgoing = ev.outgoing();
for( size_t i=0; i<outgoing.size(); ++i ){
auto const & out{ outgoing[i] };
if(is_AWZH_boson(out.type)){
auto const decay = ev.decays().find(i);
// W decays (required)
if(std::abs(out.type) == ParticleID::Wp){
if( decay != ev.decays().cend() &&
valid_W_decay(out.type>0?+1:-1, decay->second)
){
--invalid_decays;
}
else return false;
}
// Higgs decays (optional)
else if(out.type == ParticleID::h){
if(decay != ev.decays().cend()) --invalid_decays;
}
// Z decays (required)
else if(out.type == ParticleID::Z_photon_mix){
if( decay != ev.decays().cend() &&
valid_Z_decay(decay->second)
){
--invalid_decays;
}
else return false;
}
}
else if(! is_parton(out.type)) return false;
}
// any invalid decays?
if(invalid_decays != 0) return false;
return true;
}
/**
* returns all EventTypes implemented in HEJ
*/
size_t implemented_types(std::vector<Particle> const & bosons){
using namespace event_type;
if(bosons.empty()) return FKL | unob | unof | qqxexb | qqxexf | qqxmid;
if(bosons.size() == 1) {
switch (bosons[0].type) {
case ParticleID::Wp:
case ParticleID::Wm:
return FKL | unob | unof | qqxexb | qqxexf | qqxmid;
case ParticleID::Z_photon_mix:
return FKL | unob | unof;
case ParticleID::h:
return FKL | unob | unof;
default:
return non_resummable;
}
}
if(bosons.size() == 2) {
if(bosons[0].type == ParticleID::Wp && bosons[1].type == ParticleID::Wp) {
return FKL;
}
}
return non_resummable;
}
/**
* \brief function which determines if type change is consistent with Wp emission.
* @param in incoming Particle id
* @param out outgoing Particle id
* @param qqx Current both incoming/both outgoing?
*
* \see is_Wm_Change
*/
bool is_Wp_Change(ParticleID in, ParticleID out, bool qqx){
if(!qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out== (in-1);
if( qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out==-(in+1);
return false;
}
/**
* \brief function which determines if type change is consistent with Wm emission.
* @param in incoming Particle id
* @param out outgoing Particle id
* @param qqx Current both incoming/both outgoing?
*
* Ensures that change type of quark line is possible by a flavour changing
* Wm emission. Allows checking of qqx currents also.
*/
bool is_Wm_Change(ParticleID in, ParticleID out, bool qqx){
if(!qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out== (in+1);
if( qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out==-(in-1);
return false;
}
/**
* \brief checks if particle type remains same from incoming to outgoing
* @param in incoming Particle
* @param out outgoing Particle
* @param qqx Current both incoming/outgoing?
*/
bool no_flavour_change(ParticleID in, ParticleID out, bool qqx){
const int qqxCurrent = qqx?-1:1;
if(std::abs(in)<=pid::top || in==pid::gluon)
return (in==out*qqxCurrent);
return false;
}
bool has_2_jets(Event const & event){
return event.jets().size() >= 2;
}
/**
* \brief check if we have a valid Impact factor
* @param in incoming Particle
* @param out outgoing Particle
* @param qqx Current both incoming/outgoing?
* @param W_change returns +1 if Wp, -1 if Wm, else 0
*/
bool is_valid_impact_factor(
ParticleID in, ParticleID out, bool qqx, int & W_change
){
if( no_flavour_change(in, out, qqx) ){
return true;
}
if( is_Wp_Change(in, out, qqx) ) {
W_change+=1;
return true;
}
if( is_Wm_Change(in, out, qqx) ) {
W_change-=1;
return true;
}
return false;
}
//! Returns all possible classifications from the impact factors
// the beginning points are changed s.t. after the the classification they
// point to the beginning of the (potential) FKL chain
// sets W_change: + if Wp change
// 0 if no change
// - if Wm change
// This function can be used with forward & backwards iterators
template<class OutIterator>
size_t possible_impact_factors(
ParticleID incoming_id, // incoming
OutIterator & begin_out, OutIterator const & end_out, // outgoing
int & W_change, std::vector<Particle> const & boson,
bool const backward // backward?
){
using namespace event_type;
// keep track of all states that we don't test
size_t not_tested = qqxmid;
if(backward)
not_tested |= unof | qqxexf;
else
not_tested |= unob | qqxexb;
// Is this LL current?
if( is_valid_impact_factor(incoming_id, begin_out->type, false, W_change) ){
++begin_out;
return not_tested | FKL;
}
// or NLL current?
// -> needs two partons in two different jets
if( std::distance(begin_out, end_out)>=2
){
auto next = std::next(begin_out);
// Is this unordered emisson?
if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){
if( is_valid_impact_factor(
incoming_id, next->type, false, W_change )
){
// veto Higgs inside uno
assert(next!=end_out);
if( !boson.empty() && boson.front().type == ParticleID::h
){
if( (backward && boson.front().rapidity() < next->rapidity())
||(!backward && boson.front().rapidity() > next->rapidity()))
return non_resummable;
}
begin_out = std::next(next);
return not_tested | (backward?unob:unof);
}
}
// Is this QQbar?
else if( incoming_id==pid::gluon ){
if( is_valid_impact_factor(
begin_out->type, next->type, true, W_change )
){
// veto Higgs inside qqx
assert(next!=end_out);
if( !boson.empty() && boson.front().type == ParticleID::h
){
if( (backward && boson.front().rapidity() < next->rapidity())
||(!backward && boson.front().rapidity() > next->rapidity()))
return non_resummable;
}
begin_out = std::next(next);
return not_tested | (backward?qqxexb:qqxexf);
}
}
}
return non_resummable;
}
//! Returns all possible classifications from central emissions
// the beginning points are changed s.t. after the the classification they
// point to the end of the emission chain
// sets W_change: + if Wp change
// 0 if no change
// - if Wm change
template<class OutIterator>
size_t possible_central(
OutIterator & begin_out, OutIterator const & end_out,
int & W_change, std::vector<Particle> const & boson
){
using namespace event_type;
// if we already passed the central chain,
// then it is not a valid all-order state
if(std::distance(begin_out, end_out) < 0) return non_resummable;
// keep track of all states that we don't test
size_t possible = unob | unof
| qqxexb | qqxexf;
// Find the first non-gluon/non-FKL
while( (begin_out->type==pid::gluon) && (begin_out!=end_out) ){
++begin_out;
}
// end of chain -> FKL
if( begin_out==end_out ){
return possible | FKL;
}
// is this a qqbar-pair?
// needs two partons in two separate jets
auto next = std::next(begin_out);
if( is_valid_impact_factor(
begin_out->type, next->type, true, W_change )
){
// veto Higgs inside qqx
if( !boson.empty() && boson.front().type == ParticleID::h
&& boson.front().rapidity() > begin_out->rapidity()
&& boson.front().rapidity() < next->rapidity()
){
return non_resummable;
}
begin_out = std::next(next);
// remaining chain should be pure gluon/FKL
for(; begin_out!=end_out; ++begin_out){
if(begin_out->type != pid::gluon) return non_resummable;
}
return possible | qqxmid;
}
return non_resummable;
}
/**
* \brief Checks for all event types
* @param ev Event
* @returns Event Type
*
*/
event_type::EventType classify(Event const & ev){
using namespace event_type;
if(! has_2_jets(ev))
return no_2_jets;
// currently we can't handle multiple boson states in the ME. So they are
// considered "bad_final_state" even though the "classify" could work with
// them.
if(! final_state_ok(ev))
return bad_final_state;
// initialise variables
auto const & in = ev.incoming();
// range for current checks
auto begin_out{ev.cbegin_partons()};
auto end_out{ev.crbegin_partons()};
assert(std::distance(begin(in), end(in)) == 2);
assert(std::distance(begin_out, end_out.base()) >= 2);
assert(std::is_sorted(begin_out, end_out.base(), rapidity_less{}));
auto const bosons{ filter_AWZH_bosons(ev.outgoing()) };
// keep track of potential W couplings, at the end the sum should be 0
int remaining_Wp = 0;
int remaining_Wm = 0;
for(auto const & boson : bosons){
if(boson.type == ParticleID::Wp) ++remaining_Wp;
else if(boson.type == ParticleID::Wm) ++remaining_Wm;
}
size_t final_type = ~(no_2_jets | bad_final_state);
// check forward impact factor
int W_change = 0;
final_type &= possible_impact_factors(
in.front().type,
begin_out, end_out.base(),
W_change, bosons, true );
if( final_type == non_resummable )
return non_resummable;
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
// check backward impact factor
W_change = 0;
final_type &= possible_impact_factors(
in.back().type,
end_out, std::make_reverse_iterator(begin_out),
W_change, bosons, false );
if( final_type == non_resummable )
return non_resummable;
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
// check central emissions
W_change = 0;
final_type &= possible_central(
begin_out, end_out.base(), W_change, bosons );
if( final_type == non_resummable )
return non_resummable;
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
// Check whether the right number of Ws are present
if( remaining_Wp != 0 || remaining_Wm != 0 ) return non_resummable;
// result has to be unique
if( (final_type & (final_type-1)) != 0) return non_resummable;
// check that each sub processes is implemented
// (has to be done at the end)
if( (final_type & ~implemented_types(bosons)) != 0 )
return non_resummable;
return static_cast<EventType>(final_type);
}
//@}
Particle extract_particle(LHEF::HEPEUP const & hepeup, size_t i){
auto id = static_cast<ParticleID>(hepeup.IDUP[i]);
auto colour = is_parton(id)?hepeup.ICOLUP[i]:optional<Colour>();
return { id,
{ hepeup.PUP[i][0], hepeup.PUP[i][1],
hepeup.PUP[i][2], hepeup.PUP[i][3] },
colour
};
}
bool is_decay_product(std::pair<int, int> const & mothers){
if(mothers.first == 0) return false;
return mothers.second == 0 || mothers.first == mothers.second;
}
} // namespace
Event::EventData::EventData(LHEF::HEPEUP const & hepeup){
parameters.central = EventParameters{
hepeup.scales.mur, hepeup.scales.muf, hepeup.XWGTUP
};
size_t in_idx = 0;
for (int i = 0; i < hepeup.NUP; ++i) {
// skip decay products
// we will add them later on, but we have to ensure that
// the decayed particle is added before
if(is_decay_product(hepeup.MOTHUP[i])) continue;
auto particle = extract_particle(hepeup, i);
// needed to identify mother particles for decay products
particle.p.set_user_index(i+1);
if(hepeup.ISTUP[i] == LHE_Status::in){
if(in_idx > incoming.size()) {
throw std::invalid_argument{
"Event has too many incoming particles"
};
}
incoming[in_idx++] = std::move(particle);
}
else outgoing.emplace_back(std::move(particle));
}
// add decay products
for (int i = 0; i < hepeup.NUP; ++i) {
if(!is_decay_product(hepeup.MOTHUP[i])) continue;
const int mother_id = hepeup.MOTHUP[i].first;
const auto mother = std::find_if(
begin(outgoing), end(outgoing),
[mother_id](Particle const & particle){
return particle.p.user_index() == mother_id;
}
);
if(mother == end(outgoing)){
throw std::invalid_argument{"invalid decay product parent"};
}
const int mother_idx = std::distance(begin(outgoing), mother);
assert(mother_idx >= 0);
decays[mother_idx].emplace_back(extract_particle(hepeup, i));
}
}
Event::Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double const min_jet_pt
):
Event( Event::EventData{
ev.incoming, ev.outgoing, ev.decays,
Parameters<EventParameters>{ev.central, ev.variations}
}.cluster(jet_def, min_jet_pt) )
{}
//! @TODO remove in HEJ 2.2.0
UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup){
Event::EventData const evData{hepeup};
incoming = evData.incoming;
outgoing = evData.outgoing;
decays = evData.decays;
central = evData.parameters.central;
variations = evData.parameters.variations;
}
void Event::EventData::sort(){
// sort particles
std::sort(
begin(incoming), end(incoming),
[](Particle const & o1, Particle const & o2){return o1.p.pz()<o2.p.pz();}
);
auto old_outgoing = std::move(outgoing);
std::vector<size_t> idx(old_outgoing.size());
std::iota(idx.begin(), idx.end(), 0);
std::sort(idx.begin(), idx.end(), [&old_outgoing](size_t i, size_t j){
return old_outgoing[i].rapidity() < old_outgoing[j].rapidity();
});
outgoing.clear();
outgoing.reserve(old_outgoing.size());
for(size_t i: idx) {
outgoing.emplace_back(std::move(old_outgoing[i]));
}
// find decays again
if(!decays.empty()){
auto old_decays = std::move(decays);
decays.clear();
for(size_t i=0; i<idx.size(); ++i) {
auto decay = old_decays.find(idx[i]);
if(decay != old_decays.end())
decays.emplace(i, std::move(decay->second));
}
assert(old_decays.size() == decays.size());
}
}
namespace {
// use valid_X_decay to determine boson type
ParticleID reconstruct_type(std::vector<Particle> const & progeny) {
if(valid_W_decay(+1, progeny)) { return ParticleID::Wp; }
if(valid_W_decay(-1, progeny)) { return ParticleID::Wm; }
if(valid_Z_decay(progeny)) { return ParticleID::Z_photon_mix; }
throw not_implemented{
"final state with decay X -> "
+ name(progeny[0].type)
+ " + "
+ name(progeny[1].type)
};
}
// reconstruct particle with explicit ParticleID
Particle reconstruct_boson(
std::vector<Particle> const & progeny,
ParticleID const & type
) {
Particle progenitor;
progenitor.p = progeny[0].p + progeny[1].p;
progenitor.type = type;
return progenitor;
}
// reconstruct via call to reconstruct_type
Particle reconstruct_boson(std::vector<Particle> const & progeny) {
Particle progenitor {reconstruct_boson(progeny, reconstruct_type(progeny))};
assert(is_AWZH_boson(progenitor));
return progenitor;
}
typedef std::vector< std::vector<Particle> > GroupedParticles;
typedef std::pair<Particle, std::vector<Particle> > Decay;
typedef std::vector< Decay > Decays;
// return groups of reconstructable progeny
std::vector<GroupedParticles> group_progeny
(
std::vector<Particle> & leptons
) {
/**
Warning: The partition in to charged/neutral leptons is valid ONLY for WW.
**/
assert(leptons.size() == 4);
auto const begin_neutrino = std::partition(
begin(leptons), end(leptons),
[](Particle const & p) {return !is_anyneutrino(p);}
);
std::vector<Particle> neutrinos (begin_neutrino, end(leptons));
leptons.erase(begin_neutrino, end(leptons));
std::sort(begin(leptons), end(leptons), type_less{});
std::sort(begin(neutrinos), end(neutrinos), type_less{});
if(leptons.size() != 2 && neutrinos.size() != 2) { return {}; }
assert(leptons.size() == 2 && neutrinos.size() == 2);
std::vector< GroupedParticles > candidate_grouping {
{{leptons[0], neutrinos[0]}, {leptons[1], neutrinos[1]}},
{{leptons[1], neutrinos[0]}, {leptons[0], neutrinos[1]}}
};
// erase groupings containing invalid decays
candidate_grouping.erase(
std::remove_if(begin(candidate_grouping), end(candidate_grouping),
[] (GroupedParticles & candidate) -> bool {
return ! std::accumulate(
cbegin(candidate), cend(candidate), true,
[] (bool acc_valid, std::vector<Particle> const & decay) -> bool {
return acc_valid && valid_decay(decay);
}
);
}
),
candidate_grouping.end()
);
return candidate_grouping;
}
// 'best' decay ordering measure
double decay_measure(Particle reconstructed, EWConstants const & params) {
ParticleProperties ref = params.prop(reconstructed.type);
return std::abs(reconstructed.p.m() - ref.mass);
}
// decay_measure accumulated over decays
double decay_measure(Decays const & decays, EWConstants const & params) {
return
std::accumulate(
cbegin(decays), cend(decays), 0.,
[&params] (double dm, Decay const & decay) -> double {
return dm + decay_measure(decay.first, params);
}
);
}
// select best combination of decays for the event
Decays select_decays
(
std::vector<Particle> & leptons,
EWConstants const & ew_parameters
) {
std::vector<GroupedParticles> groupings = group_progeny(leptons);
std::vector<Decays> valid_decays;
valid_decays.reserve(groupings.size());
// Reconstruct all groupings
for(GroupedParticles const & group : groupings) {
Decays decays;
for(auto const & progeny : group) {
decays.emplace_back(make_pair(reconstruct_boson(progeny), progeny));
}
valid_decays.emplace_back(decays);
}
if (valid_decays.empty()) {
throw not_implemented{"No supported intermediate reconstruction available"};
}
else if (valid_decays.size() == 1) { return valid_decays[0]; }
else {
// select decay with smallest decay_measure
auto selected = std::min_element(cbegin(valid_decays), cend(valid_decays),
[&ew_parameters] (auto const & d1, auto const & d2) -> bool {
return decay_measure(d1, ew_parameters) < decay_measure(d2, ew_parameters);
}
);
return *selected;
}
}
} // namespace
void Event::EventData::reconstruct_intermediate(EWConstants const & ew_parameters) {
auto const begin_leptons = std::partition(
begin(outgoing), end(outgoing),
[](Particle const & p) {return !is_anylepton(p);}
);
if(begin_leptons == end(outgoing)) return;
assert(is_anylepton(*begin_leptons));
std::vector<Particle> leptons(begin_leptons, end(outgoing));
outgoing.erase(begin_leptons, end(outgoing));
if(leptons.empty()) { return; } // nothing to do
else if(leptons.size() == 2) {
outgoing.emplace_back(reconstruct_boson(leptons));
std::sort(begin(leptons), end(leptons), type_less{});
decays.emplace(outgoing.size()-1, std::move(leptons));
}
else if(leptons.size() == 4) {
// select_decays only supports WpWp
Decays valid_decays = select_decays(leptons, ew_parameters);
for(auto &decay : valid_decays) {
outgoing.emplace_back(decay.first);
std::sort(begin(decay.second), end(decay.second), type_less{});
decays.emplace(outgoing.size()-1, std::move(decay.second));
}
}
else {
throw not_implemented {
std::to_string(leptons.size())
+ " leptons in the final state"
};
}
}
Event Event::EventData::cluster(
fastjet::JetDefinition const & jet_def, double const min_jet_pt
){
sort();
return Event{ std::move(incoming), std::move(outgoing), std::move(decays),
std::move(parameters),
jet_def, min_jet_pt
};
}
Event::Event(
std::array<Particle, 2> && incoming,
std::vector<Particle> && outgoing,
std::unordered_map<size_t, std::vector<Particle>> && decays,
Parameters<EventParameters> && parameters,
fastjet::JetDefinition const & jet_def,
double const min_jet_pt
): incoming_{std::move(incoming)},
outgoing_{std::move(outgoing)},
decays_{std::move(decays)},
parameters_{std::move(parameters)},
cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def },
min_jet_pt_{min_jet_pt}
{
jets_ = sorted_by_rapidity(cs_.inclusive_jets(min_jet_pt_));
assert(std::is_sorted(begin(outgoing_), end(outgoing_),
rapidity_less{}));
type_ = classify(*this);
}
namespace {
//! check that Particles have a reasonable colour
bool correct_colour(Particle const & part){
ParticleID id{ part.type };
if(!is_parton(id))
return !part.colour;
if(!part.colour)
return false;
Colour const & col{ *part.colour };
if(is_quark(id))
return col.first != 0 && col.second == 0;
if(is_antiquark(id))
return col.first == 0 && col.second != 0;
assert(id==ParticleID::gluon);
return col.first != 0 && col.second != 0 && col.first != col.second;
}
//! Connect parton to t-channel colour line & update the line
//! returns false if connection not possible
template<class OutIterator>
bool try_connect_t(OutIterator const & it_part, Colour & line_colour){
if( line_colour.first == it_part->colour->second ){
line_colour.first = it_part->colour->first;
return true;
}
if( line_colour.second == it_part->colour->first ){
line_colour.second = it_part->colour->second;
return true;
}
return false;
}
//! Connect parton to u-channel colour line & update the line
//! returns false if connection not possible
template<class OutIterator>
bool try_connect_u(OutIterator & it_part, Colour & line_colour){
auto it_next = std::next(it_part);
if( try_connect_t(it_next, line_colour)
&& try_connect_t(it_part, line_colour)
){
it_part=it_next;
return true;
}
return false;
}
} // namespace
bool Event::is_leading_colour() const {
if( !correct_colour(incoming()[0]) || !correct_colour(incoming()[1]) )
return false;
Colour line_colour = *incoming()[0].colour;
std::swap(line_colour.first, line_colour.second);
// reasonable colour
if(!std::all_of(outgoing().cbegin(), outgoing().cend(), correct_colour))
return false;
for(auto it_part = cbegin_partons(); it_part!=cend_partons(); ++it_part){
switch (type()) {
case event_type::FKL:
if( !try_connect_t(it_part, line_colour) )
return false;
break;
case event_type::unob:
case event_type::qqxexb: {
if( !try_connect_t(it_part, line_colour)
// u-channel only allowed at impact factor
&& (std::distance(cbegin_partons(), it_part)!=0
|| !try_connect_u(it_part, line_colour)))
return false;
break;
}
case event_type::unof:
case event_type::qqxexf: {
if( !try_connect_t(it_part, line_colour)
// u-channel only allowed at impact factor
&& (std::distance(it_part, cend_partons())!=2
|| !try_connect_u(it_part, line_colour)))
return false;
break;
}
case event_type::qqxmid:{
auto it_next = std::next(it_part);
if( !try_connect_t(it_part, line_colour)
// u-channel only allowed at qqx/qxq pair
&& ( ( !(is_quark(*it_part) && is_antiquark(*it_next))
&& !(is_antiquark(*it_part) && is_quark(*it_next)))
|| !try_connect_u(it_part, line_colour))
)
return false;
break;
}
default:
throw std::logic_error{"unreachable"};
}
// no colour singlet exchange/disconnected diagram
if(line_colour.first == line_colour.second)
return false;
}
return (incoming()[1].colour->first == line_colour.first)
&& (incoming()[1].colour->second == line_colour.second);
}
namespace {
//! connect incoming Particle to colour flow
void connect_incoming(Particle & in, int & colour, int & anti_colour){
in.colour = std::make_pair(anti_colour, colour);
// gluon
if(in.type == pid::gluon)
return;
if(in.type > 0){
// quark
assert(is_quark(in));
in.colour->second = 0;
colour*=-1;
return;
}
// anti-quark
assert(is_antiquark(in));
in.colour->first = 0;
anti_colour*=-1;
}
//! connect outgoing Particle to t-channel colour flow
template<class OutIterator>
void connect_tchannel(
OutIterator & it_part, int & colour, int & anti_colour, RNG & ran
){
assert(colour>0 || anti_colour>0);
if(it_part->type == ParticleID::gluon){
// gluon
if(colour>0 && anti_colour>0){
// on g line => connect to colour OR anti-colour (random)
if(ran.flat() < 0.5){
it_part->colour = std::make_pair(colour+2,colour);
colour+=2;
} else {
it_part->colour = std::make_pair(anti_colour, anti_colour+2);
anti_colour+=2;
}
} else if(colour > 0){
// on q line => connect to available colour
it_part->colour = std::make_pair(colour+2, colour);
colour+=2;
} else {
assert(colour<0 && anti_colour>0);
// on qx line => connect to available anti-colour
it_part->colour = std::make_pair(anti_colour, anti_colour+2);
anti_colour+=2;
}
} else if(is_quark(*it_part)) {
// quark
assert(anti_colour>0);
if(colour>0){
// on g line => connect and remove anti-colour
it_part->colour = std::make_pair(anti_colour, 0);
anti_colour+=2;
anti_colour*=-1;
} else {
// on qx line => new colour
colour*=-1;
it_part->colour = std::make_pair(colour, 0);
}
} else if(is_antiquark(*it_part)) {
// anti-quark
assert(colour>0);
if(anti_colour>0){
// on g line => connect and remove colour
it_part->colour = std::make_pair(0, colour);
colour+=2;
colour*=-1;
} else {
// on q line => new anti-colour
anti_colour*=-1;
it_part->colour = std::make_pair(0, anti_colour);
}
} else { // not a parton
assert(!is_parton(*it_part));
it_part->colour = {};
}
}
//! connect to t- or u-channel colour flow
template<class OutIterator>
void connect_utchannel(
OutIterator & it_part, int & colour, int & anti_colour, RNG & ran
){
OutIterator it_first = it_part++;
if(ran.flat()<.5) {// t-channel
connect_tchannel(it_first, colour, anti_colour, ran);
connect_tchannel(it_part, colour, anti_colour, ran);
}
else { // u-channel
connect_tchannel(it_part, colour, anti_colour, ran);
connect_tchannel(it_first, colour, anti_colour, ran);
}
}
} // namespace
bool Event::generate_colours(RNG & ran){
// generate only for HEJ events
if(!event_type::is_resummable(type()))
return false;
assert(std::is_sorted(
begin(outgoing()), end(outgoing()), rapidity_less{}));
assert(incoming()[0].pz() < incoming()[1].pz());
// positive (anti-)colour -> can connect
// negative (anti-)colour -> not available/used up by (anti-)quark
int colour = COLOUR_OFFSET;
int anti_colour = colour+1;
// initialise first
connect_incoming(incoming_[0], colour, anti_colour);
// reset outgoing colours
std::for_each(outgoing_.begin(), outgoing_.end(),
[](Particle & part){ part.colour = {};});
for(auto it_part = begin_partons(); it_part!=end_partons(); ++it_part){
switch (type()) {
// subleading can connect to t- or u-channel
case event_type::unob:
case event_type::qqxexb: {
if( std::distance(begin_partons(), it_part)==0)
connect_utchannel(it_part, colour, anti_colour, ran);
else
connect_tchannel(it_part, colour, anti_colour, ran);
break;
}
case event_type::unof:
case event_type::qqxexf: {
if( std::distance(it_part, end_partons())==2)
connect_utchannel(it_part, colour, anti_colour, ran);
else
connect_tchannel(it_part, colour, anti_colour, ran);
break;
}
case event_type::qqxmid:{
auto it_next = std::next(it_part);
if( std::distance(begin_partons(), it_part)>0
&& std::distance(it_part, end_partons())>2
&& ( (is_quark(*it_part) && is_antiquark(*it_next))
|| (is_antiquark(*it_part) && is_quark(*it_next)) )
)
connect_utchannel(it_part, colour, anti_colour, ran);
else
connect_tchannel(it_part, colour, anti_colour, ran);
break;
}
default: // rest has to be t-channel
connect_tchannel(it_part, colour, anti_colour, ran);
}
}
// Connect last
connect_incoming(incoming_[1], anti_colour, colour);
assert(is_leading_colour());
return true;
} // generate_colours
namespace {
bool valid_parton(
std::vector<fastjet::PseudoJet> const & jets,
Particle const & parton, int const idx,
double const max_ext_soft_pt_fraction, double const min_extparton_pt
){
// TODO code overlap with PhaseSpacePoint::pass_extremal_cuts
if(min_extparton_pt > parton.pt()) return false;
if(idx<0) return false;
assert(static_cast<int>(jets.size())>=idx);
auto const & jet{ jets[idx] };
return (parton.p - jet).pt()/jet.pt() <= max_ext_soft_pt_fraction;
}
} // namespace
// this should work with multiple types
bool Event::valid_hej_state(double const max_frac,
double const min_pt
) const {
using namespace event_type;
if(!is_resummable(type()))
return false;
auto const & jet_idx{ particle_jet_indices() };
auto idx_begin{ jet_idx.cbegin() };
auto idx_end{ jet_idx.crbegin() };
auto part_begin{ cbegin_partons() };
auto part_end{ crbegin_partons() };
// always seperate extremal jets
if( !valid_parton(jets(), *part_begin, *idx_begin, max_frac, min_pt) )
return false;
++part_begin;
++idx_begin;
if( !valid_parton(jets(), *part_end, *idx_end, max_frac, min_pt) )
return false;
++part_end;
++idx_end;
// unob -> second parton in own jet
if( type() & (unob | qqxexb) ){
if( !valid_parton(jets(), *part_begin, *idx_begin, max_frac, min_pt) )
return false;
++part_begin;
++idx_begin;
}
if( type() & (unof | qqxexf) ){
if( !valid_parton(jets(), *part_end, *idx_end, max_frac, min_pt) )
return false;
++part_end;
// ++idx_end; // last check, we don't need idx_end afterwards
}
if( type() & qqxmid ){
// find qqx pair
auto begin_qqx{ std::find_if( part_begin, part_end.base(),
[](Particle const & part) -> bool {
return part.type != ParticleID::gluon;
}
)};
assert(begin_qqx != part_end.base());
long int qqx_pos{ std::distance(part_begin, begin_qqx) };
assert(qqx_pos >= 0);
idx_begin+=qqx_pos;
if( !( valid_parton(jets(),*begin_qqx, *idx_begin, max_frac,min_pt)
&& valid_parton(jets(),*(++begin_qqx),*(++idx_begin),max_frac,min_pt)
))
return false;
}
return true;
}
Event::ConstPartonIterator Event::begin_partons() const {
return cbegin_partons();
}
Event::ConstPartonIterator Event::cbegin_partons() const {
return boost::make_filter_iterator(
static_cast<bool (*)(Particle const &)>(is_parton),
cbegin(outgoing()),
cend(outgoing())
);
}
Event::ConstPartonIterator Event::end_partons() const {
return cend_partons();
}
Event::ConstPartonIterator Event::cend_partons() const {
return boost::make_filter_iterator(
static_cast<bool (*)(Particle const &)>(is_parton),
cend(outgoing()),
cend(outgoing())
);
}
Event::ConstReversePartonIterator Event::rbegin_partons() const {
return crbegin_partons();
}
Event::ConstReversePartonIterator Event::crbegin_partons() const {
return std::reverse_iterator<ConstPartonIterator>( cend_partons() );
}
Event::ConstReversePartonIterator Event::rend_partons() const {
return crend_partons();
}
Event::ConstReversePartonIterator Event::crend_partons() const {
return std::reverse_iterator<ConstPartonIterator>( cbegin_partons() );
}
Event::PartonIterator Event::begin_partons() {
return boost::make_filter_iterator(
static_cast<bool (*)(Particle const &)>(is_parton),
begin(outgoing_),
end(outgoing_)
);
}
Event::PartonIterator Event::end_partons() {
return boost::make_filter_iterator(
static_cast<bool (*)(Particle const &)>(is_parton),
end(outgoing_),
end(outgoing_)
);
}
Event::ReversePartonIterator Event::rbegin_partons() {
return std::reverse_iterator<PartonIterator>( end_partons() );
}
Event::ReversePartonIterator Event::rend_partons() {
return std::reverse_iterator<PartonIterator>( begin_partons() );
}
namespace {
void print_momentum(std::ostream & os, fastjet::PseudoJet const & part){
constexpr int prec = 6;
const std::streamsize orig_prec = os.precision();
os <<std::scientific<<std::setprecision(prec) << "["
<<std::setw(2*prec+1)<<std::right<< part.px() << ", "
<<std::setw(2*prec+1)<<std::right<< part.py() << ", "
<<std::setw(2*prec+1)<<std::right<< part.pz() << ", "
<<std::setw(2*prec+1)<<std::right<< part.E() << "]"<< std::fixed;
os.precision(orig_prec);
}
void print_colour(std::ostream & os, optional<Colour> const & col){
constexpr int width = 3;
if(!col)
os << "(no color)"; // American spelling for better alignment
else
os << "(" <<std::setw(width)<<std::right<< col->first
<< ", " <<std::setw(width)<<std::right<< col->second << ")";
}
} // namespace
std::ostream& operator<<(std::ostream & os, Event const & ev){
constexpr int prec = 4;
constexpr int wtype = 3; // width for types
const std::streamsize orig_prec = os.precision();
os <<std::setprecision(prec)<<std::fixed;
os << "########## " << name(ev.type()) << " ##########" << std::endl;
os << "Incoming particles:\n";
for(auto const & in: ev.incoming()){
os <<std::setw(wtype)<< in.type << ": ";
print_colour(os, in.colour);
os << " ";
print_momentum(os, in.p);
os << std::endl;
}
os << "\nOutgoing particles: " << ev.outgoing().size() << "\n";
for(auto const & out: ev.outgoing()){
os <<std::setw(wtype)<< out.type << ": ";
print_colour(os, out.colour);
os << " ";
print_momentum(os, out.p);
os << " => rapidity="
<<std::setw(2*prec-1)<<std::right<< out.rapidity() << std::endl;
}
os << "\nForming Jets: " << ev.jets().size() << "\n";
for(auto const & jet: ev.jets()){
print_momentum(os, jet);
os << " => rapidity="
<<std::setw(2*prec-1)<<std::right<< jet.rapidity() << std::endl;
}
if(!ev.decays().empty() ){
os << "\nDecays: " << ev.decays().size() << "\n";
for(auto const & decay: ev.decays()){
os <<std::setw(wtype)<< ev.outgoing()[decay.first].type
<< " (" << decay.first << ") to:\n";
for(auto const & out: decay.second){
os <<" "<<std::setw(wtype)<< out.type << ": ";
print_momentum(os, out.p);
os << " => rapidity="
<<std::setw(2*prec-1)<<std::right<< out.rapidity() << std::endl;
}
}
}
os << std::defaultfloat;
os.precision(orig_prec);
return os;
}
double shat(Event const & ev){
return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
}
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){
LHEF::HEPEUP result;
result.heprup = heprup;
result.weights = {{event.central().weight, nullptr}};
for(auto const & var: event.variations()){
result.weights.emplace_back(var.weight, nullptr);
}
size_t num_particles = event.incoming().size() + event.outgoing().size();
for(auto const & decay: event.decays()) num_particles += decay.second.size();
result.NUP = num_particles;
// the following entries are pretty much meaningless
result.IDPRUP = event.type(); // event type
result.AQEDUP = 1./128.; // alpha_EW
//result.AQCDUP = 0.118 // alpha_QCD
// end meaningless part
result.XWGTUP = event.central().weight;
result.SCALUP = event.central().muf;
result.scales.muf = event.central().muf;
result.scales.mur = event.central().mur;
result.scales.SCALUP = event.central().muf;
result.pdfinfo.p1 = event.incoming().front().type;
result.pdfinfo.p2 = event.incoming().back().type;
result.pdfinfo.scale = event.central().muf;
result.IDUP.reserve(num_particles); // PID
result.ISTUP.reserve(num_particles); // status (in, out, decay)
result.PUP.reserve(num_particles); // momentum
result.MOTHUP.reserve(num_particles); // index mother particle
result.ICOLUP.reserve(num_particles); // colour
// incoming
std::array<Particle, 2> incoming{ event.incoming() };
// First incoming should be positive pz according to LHE standard
// (or at least most (everyone?) do it this way, and Pythia assumes it)
if(incoming[0].pz() < incoming[1].pz())
std::swap(incoming[0], incoming[1]);
for(Particle const & in: incoming){
result.IDUP.emplace_back(in.type);
result.ISTUP.emplace_back(LHE_Status::in);
result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()});
result.MOTHUP.emplace_back(0, 0);
assert(in.colour);
result.ICOLUP.emplace_back(*in.colour);
}
// outgoing
for(size_t i = 0; i < event.outgoing().size(); ++i){
Particle const & out = event.outgoing()[i];
result.IDUP.emplace_back(out.type);
const int status = event.decays().count(i) != 0u
?LHE_Status::decay
:LHE_Status::out;
result.ISTUP.emplace_back(status);
result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
result.MOTHUP.emplace_back(1, 2);
if(out.colour)
result.ICOLUP.emplace_back(*out.colour);
else{
result.ICOLUP.emplace_back(std::make_pair(0,0));
}
}
// decays
for(auto const & decay: event.decays()){
for(auto const & out: decay.second){
result.IDUP.emplace_back(out.type);
result.ISTUP.emplace_back(LHE_Status::out);
result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
const size_t mother_idx = 1 + event.incoming().size() + decay.first;
result.MOTHUP.emplace_back(mother_idx, mother_idx);
result.ICOLUP.emplace_back(0,0);
}
}
assert(result.ICOLUP.size() == num_particles);
static constexpr double unknown_spin = 9.; //per Les Houches accord
result.VTIMUP = std::vector<double>(num_particles, unknown_spin);
result.SPINUP = result.VTIMUP;
return result;
}
} // namespace HEJ
diff --git a/t/cmp_events.cc b/t/cmp_events.cc
index bd4c540..057678f 100644
--- a/t/cmp_events.cc
+++ b/t/cmp_events.cc
@@ -1,67 +1,68 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/Event.hh"
+#include "HEJ/EWConstants.hh"
#include "HEJ/stream.hh"
#include "HEJ/utility.hh"
#include "LHEF/LHEF.h"
namespace {
const HEJ::ParticleProperties Wprop{80.385, 2.085};
const HEJ::ParticleProperties Zprop{91.187, 2.495};
const HEJ::ParticleProperties Hprop{125, 0.004165};
constexpr double vev = 246.2196508;
const HEJ::EWConstants ew_parameters{vev, Wprop, Zprop, Hprop};
}
int main(int argn, char** argv) {
if(argn != 3){
std::cerr << "Usage: cmp_events eventfile eventfile_no_boson ";
return EXIT_FAILURE;
}
HEJ::istream in{argv[1]};
LHEF::Reader reader{in};
HEJ::istream in_no_boson{argv[2]};
LHEF::Reader reader_no_boson{in_no_boson};
while(reader.readEvent()) {
if(! reader_no_boson.readEvent()) {
std::cerr << "wrong number of events in " << argv[2] << '\n';
return EXIT_FAILURE;
}
const auto is_AWZH = [](HEJ::Particle const & p) {
return is_AWZH_boson(p);
};
const HEJ::Event::EventData event_data{reader.hepeup};
const auto boson = std::find_if(
begin(event_data.outgoing), end(event_data.outgoing), is_AWZH
);
if(boson == end(event_data.outgoing)) {
std::cerr << "no boson in event\n";
return EXIT_FAILURE;
}
HEJ::Event::EventData data_no_boson{reader_no_boson.hepeup};
data_no_boson.reconstruct_intermediate(ew_parameters);
const auto new_boson = std::find_if(
begin(event_data.outgoing), end(event_data.outgoing), is_AWZH
);
if(new_boson == end(data_no_boson.outgoing)) {
std::cerr << "failed to find reconstructed boson\n";
return EXIT_FAILURE;
}
if(new_boson->type != boson->type) {
std::cerr << "reconstructed wrong boson type\n";
return EXIT_FAILURE;
}
if(!HEJ::nearby(new_boson->p, boson->p)) {
std::cerr << "reconstructed boson momentum is off\n";
return EXIT_FAILURE;
}
}
return EXIT_SUCCESS;
}
diff --git a/t/hej_test.cc b/t/hej_test.cc
index ecf2dd2..389b7ed 100644
--- a/t/hej_test.cc
+++ b/t/hej_test.cc
@@ -1,624 +1,625 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "hej_test.hh"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <iterator>
#include <memory>
#include <numeric>
#include <random>
#include <utility>
+#include "HEJ/EWConstants.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PDG_codes.hh"
namespace {
const HEJ::ParticleProperties Wprop{80.385, 2.085};
const HEJ::ParticleProperties Zprop{91.187, 2.495};
const HEJ::ParticleProperties Hprop{125, 0.004165};
constexpr double vev = 246.2196508;
const HEJ::EWConstants ew_parameters{vev, Wprop, Zprop, Hprop};
}
HEJ::Event::EventData get_process(int const njet, int const pos_boson){
using namespace HEJ::pid;
HEJ::Event::EventData ev;
if(njet == 0){ // jet idx: -1 -1
ev.outgoing.push_back({gluon, { -24, 12, -57, 63}, {}});
ev.outgoing.push_back({gluon, { 24, -12, 41, 49}, {}});
ev.incoming[0] = {gluon, { 0, 0, -64, 64}, {}};
ev.incoming[1] = {gluon, { 0, 0, 48, 48}, {}};
return ev;
}
else if(njet == 1){ // jet idx: 0 -1 -1
ev.outgoing.push_back({gluon, { 23, 28, -44, 57}, {}});
ev.outgoing.push_back({gluon, { -11, -24, -12, 29}, {}});
ev.outgoing.push_back({gluon, { -12, -4, 39, 41}, {}});
ev.incoming[0] = {gluon, { 0, 0, -72, 72}, {}};
ev.incoming[1] = {gluon, { 0, 0, 55, 55}, {}};
return ev;
}
else if(njet == 2){
switch(pos_boson){
case 0:
ev.outgoing.push_back({higgs, { 198, 33, -170, 291}, {}});
ev.outgoing.push_back({gluon, {-154, 68, 44, 174}, {}});
ev.outgoing.push_back({gluon, { -44, -101, 88, 141}, {}});
ev.incoming[0] = {gluon, { 0, 0, -322, 322}, {}};
ev.incoming[1] = {gluon, { 0, 0, 284, 284}, {}};
return ev;
case 1:
ev.outgoing.push_back({gluon, { -6, 82, -159, 179}, {}});
ev.outgoing.push_back({higgs, { 195, -106, 74, 265}, {}});
ev.outgoing.push_back({gluon, {-189, 24, 108, 219}, {}});
ev.incoming[0] = {gluon, { 0, 0, -320, 320}, {}};
ev.incoming[1] = {gluon, { 0, 0, 343, 343}, {}};
return ev;
case 2:
ev.outgoing.push_back({gluon, { -80, -80, -140, 180}, {}});
ev.outgoing.push_back({gluon, { -60, -32, 0, 68}, {}});
ev.outgoing.push_back({higgs, { 140, 112, 177, 281}, {}});
ev.incoming[0] = {gluon, { 0, 0, -246, 246}, {}};
ev.incoming[1] = {gluon, { 0, 0, 283, 283}, {}};
return ev;
default:
ev.outgoing.push_back({gluon, { -72, 24, 18, 78}, {}});
ev.outgoing.push_back({gluon, { 72, -24, 74, 106}, {}});
ev.incoming[0] = {gluon, { 0, 0, -46, 46}, {}};
ev.incoming[1] = {gluon, { 0, 0, 138, 138}, {}};
return ev;
}
}
if(njet == 3){
switch(pos_boson){
case 0:
ev.outgoing.push_back({higgs, { 152, -117, -88, 245}, {}});
ev.outgoing.push_back({gluon, {-146, 62, -11, 159}, {}});
ev.outgoing.push_back({gluon, { 126, -72, 96, 174}, {}});
ev.outgoing.push_back({gluon, {-132, 127, 144, 233}, {}});
ev.incoming[0] = {gluon, { 0, 0, -335, 335}, {}};
ev.incoming[1] = {gluon, { 0, 0, 476, 476}, {}};
return ev;
case 1:
ev.outgoing.push_back({gluon, {-191, 188, -128, 297}, {}});
ev.outgoing.push_back({higgs, { 199, 72, -76, 257}, {}});
ev.outgoing.push_back({gluon, { 184, -172, -8, 252}, {}});
ev.outgoing.push_back({gluon, {-192, -88, 54, 218}, {}});
ev.incoming[0] = {gluon, { 0, 0, -591, 591}, {}};
ev.incoming[1] = {gluon, { 0, 0, 433, 433}, {}};
return ev;
case 2:
ev.outgoing.push_back({gluon, { -42, 18, -49, 67}, {}});
ev.outgoing.push_back({gluon, { -12, -54, -28, 62}, {}});
ev.outgoing.push_back({higgs, { 99, 32, -16, 163}, {}});
ev.outgoing.push_back({gluon, { -45, 4, 72, 85}, {}});
ev.incoming[0] = {gluon, { 0, 0, -199, 199}, {}};
ev.incoming[1] = {gluon, { 0, 0, 178, 178}, {}};
return ev;
case 3:
ev.outgoing.push_back({gluon, { -65, -32, -76, 105}, {}});
ev.outgoing.push_back({gluon, { -22, 31, -34, 51}, {}});
ev.outgoing.push_back({gluon, { -12, -67, -36, 77}, {}});
ev.outgoing.push_back({higgs, { 99, 68, -4, 173}, {}});
ev.incoming[0] = {gluon, { 0, 0, -278, 278}, {}};
ev.incoming[1] = {gluon, { 0, 0, 128, 128}, {}};
return ev;
default:
ev.outgoing.push_back({gluon, { -90, -135, 30, 165}, {}});
ev.outgoing.push_back({gluon, {-108, 198, 76, 238}, {}});
ev.outgoing.push_back({gluon, { 198, -63, 126, 243}, {}});
ev.incoming[0] = {gluon, { 0, 0, -207, 207}, {}};
ev.incoming[1] = {gluon, { 0, 0, 439, 439}, {}};
return ev;
}
}
if(njet == 4){
switch(pos_boson){
case 0:
ev.outgoing.push_back({higgs, { 199, 72, -76, 257}, {}});
ev.outgoing.push_back({gluon, {-200, -155, -64, 261}, {}});
ev.outgoing.push_back({gluon, { 198, 194, 57, 283}, {}});
ev.outgoing.push_back({gluon, { 1, 32, 8, 33}, {}});
ev.outgoing.push_back({gluon, {-198, -143, 186, 307}, {}});
ev.incoming[0] = {gluon, { 0, 0, -515, 515}, {}};
ev.incoming[1] = {gluon, { 0, 0, 626, 626}, {}};
return ev;
case 1:
ev.outgoing.push_back({gluon, { 198, 61, -162, 263}, {}});
ev.outgoing.push_back({higgs, { 199, 72, -76, 257}, {}});
ev.outgoing.push_back({gluon, {-200, 135, 144, 281}, {}});
ev.outgoing.push_back({gluon, {-198, -186, 171, 321}, {}});
ev.outgoing.push_back({gluon, { 1, -82, 122, 147}, {}});
ev.incoming[0] = {gluon, { 0, 0, -535, 535}, {}};
ev.incoming[1] = {gluon, { 0, 0, 734, 734}, {}};
return ev;
case 2:
ev.outgoing.push_back({gluon, {-180, -27, -164, 245}, {}});
ev.outgoing.push_back({gluon, {-108, 78, -36, 138}, {}});
ev.outgoing.push_back({higgs, { 196, -189, 68, 307}, {}});
ev.outgoing.push_back({gluon, {-107, 136, 76, 189}, {}});
ev.outgoing.push_back({gluon, { 199, 2, 178, 267}, {}});
ev.incoming[0] = {gluon, { 0, 0, -512, 512}, {}};
ev.incoming[1] = {gluon, { 0, 0, 634, 634}, {}};
return ev;
case 3:
ev.outgoing.push_back({gluon, { -12, -30, -84, 90}, {}});
ev.outgoing.push_back({gluon, { -72, 22, -96, 122}, {}});
ev.outgoing.push_back({gluon, { 68, 0, -51, 85}, {}});
ev.outgoing.push_back({higgs, { 64, 72, -81, 177}, {}});
ev.outgoing.push_back({gluon, { -48, -64, 84, 116}, {}});
ev.incoming[0] = {gluon, { 0, 0, -409, 409}, {}};
ev.incoming[1] = {gluon, { 0, 0, 181, 181}, {}};
return ev;
case 4:
ev.outgoing.push_back({gluon, { -72, -49, -72, 113}, {}});
ev.outgoing.push_back({gluon, { -48, 0, -36, 60}, {}});
ev.outgoing.push_back({gluon, { -12, 54, -36, 66}, {}});
ev.outgoing.push_back({gluon, { 68, -77, -56, 117}, {}});
ev.outgoing.push_back({higgs, { 64, 72, -81, 177}, {}});
ev.incoming[0] = {gluon, { 0, 0, -407, 407}, {}};
ev.incoming[1] = {gluon, { 0, 0, 126, 126}, {}};
return ev;
default:
ev.outgoing.push_back({gluon, { 248, -56, -122, 282}, {}});
ev.outgoing.push_back({gluon, { 249, 30, -10, 251}, {}});
ev.outgoing.push_back({gluon, {-249, -18, 26, 251}, {}});
ev.outgoing.push_back({gluon, {-248, 44, 199, 321}, {}});
ev.incoming[0] = {gluon, { 0, 0, -506, 506}, {}};
ev.incoming[1] = {gluon, { 0, 0, 599, 599}, {}};
return ev;
}
}
if(njet == 6){
switch(pos_boson){
case 0:
ev.outgoing.push_back({higgs, { 349, 330, -94, 505}, {}});
ev.outgoing.push_back({gluon, {-315, -300, 0, 435}, {}});
ev.outgoing.push_back({gluon, { 347, 306, 18, 463}, {}});
ev.outgoing.push_back({gluon, {-249, -342, 162, 453}, {}});
ev.outgoing.push_back({gluon, { 345, 312, 284, 545}, {}});
ev.outgoing.push_back({gluon, {-324, -126, 292, 454}, {}});
ev.outgoing.push_back({gluon, {-153, -180, 304, 385}, {}});
ev.incoming[0] = {gluon, { 0, 0, -1137, 1137}, {}};
ev.incoming[1] = {gluon, { 0, 0, 2103, 2103}, {}};
return ev;
case 1:
ev.outgoing.push_back({gluon, { 242, 241, -182, 387}, {}});
ev.outgoing.push_back({higgs, { 243, 238, -190, 409}, {}});
ev.outgoing.push_back({gluon, {-218, -215, -74, 315}, {}});
ev.outgoing.push_back({gluon, {-224, -224, 112, 336}, {}});
ev.outgoing.push_back({gluon, { 241, 182, 154, 339}, {}});
ev.outgoing.push_back({gluon, { -53, -234, 126, 271}, {}});
ev.outgoing.push_back({gluon, {-231, 12, 156, 279}, {}});
ev.incoming[0] = {gluon, { 0, 0, -1117, 1117}, {}};
ev.incoming[1] = {gluon, { 0, 0, 1219, 1219}, {}};
return ev;
case 2:
ev.outgoing.push_back({gluon, { 151, 102, -42, 187}, {}});
ev.outgoing.push_back({gluon, { -86, -46, -17, 99}, {}});
ev.outgoing.push_back({higgs, { 152, 153, 0, 249}, {}});
ev.outgoing.push_back({gluon, { -60, -135, 64, 161}, {}});
ev.outgoing.push_back({gluon, { 150, 123, 110, 223}, {}});
ev.outgoing.push_back({gluon, {-154, -49, 98, 189}, {}});
ev.outgoing.push_back({gluon, {-153, -148, 144, 257}, {}});
ev.incoming[0] = {gluon, { 0, 0, -504, 504}, {}};
ev.incoming[1] = {gluon, { 0, 0, 861, 861}, {}};
return ev;
case 3:
ev.outgoing.push_back({gluon, { 198, 197, -66, 287}, {}});
ev.outgoing.push_back({gluon, {-198, -189, -54, 279}, {}});
ev.outgoing.push_back({gluon, {-200, -64, 2, 210}, {}});
ev.outgoing.push_back({higgs, { 199, 158, 6, 283}, {}});
ev.outgoing.push_back({gluon, {-199, -184, 172, 321}, {}});
ev.outgoing.push_back({gluon, { 196, 168, 177, 313}, {}});
ev.outgoing.push_back({gluon, { 4, -86, 92, 126}, {}});
ev.incoming[0] = {gluon, { 0, 0, -745, 745}, {}};
ev.incoming[1] = {gluon, { 0, 0, 1074, 1074}, {}};
return ev;
case 4:
ev.outgoing.push_back({gluon, { 151, 102, -42, 187}, {}});
ev.outgoing.push_back({gluon, { -86, -133, -14, 159}, {}});
ev.outgoing.push_back({gluon, {-154, -104, -8, 186}, {}});
ev.outgoing.push_back({gluon, { -60, 11, 0, 61}, {}});
ev.outgoing.push_back({higgs, { 152, 153, 0, 249}, {}});
ev.outgoing.push_back({gluon, { 150, 125, 90, 215}, {}});
ev.outgoing.push_back({gluon, {-153, -154, 126, 251}, {}});
ev.incoming[0] = {gluon, { 0, 0, -578, 578}, {}};
ev.incoming[1] = {gluon, { 0, 0, 730, 730}, {}};
return ev;
case 5:
ev.outgoing.push_back({gluon, { -15, -90, -94, 131}, {}});
ev.outgoing.push_back({gluon, { -11, 82, -74, 111}, {}});
ev.outgoing.push_back({gluon, { 23, -80, -64, 105}, {}});
ev.outgoing.push_back({gluon, { -48, -25, -36, 65}, {}});
ev.outgoing.push_back({gluon, { -12, 99, -16, 101}, {}});
ev.outgoing.push_back({higgs, { 68, 92, -18, 170}, {}});
ev.outgoing.push_back({gluon, { -5, -78, 54, 95}, {}});
ev.incoming[0] = {gluon, { 0, 0, -513, 513}, {}};
ev.incoming[1] = {gluon, { 0, 0, 265, 265}, {}};
return ev;
case 6:
ev.outgoing.push_back({gluon, { 198, 197, -66, 287}, {}});
ev.outgoing.push_back({gluon, { 4, -84, -18, 86}, {}});
ev.outgoing.push_back({gluon, {-198, -60, -36, 210}, {}});
ev.outgoing.push_back({gluon, { 196, -78, -36, 214}, {}});
ev.outgoing.push_back({gluon, {-200, 45, 0, 205}, {}});
ev.outgoing.push_back({gluon, {-199, -178, 2, 267}, {}});
ev.outgoing.push_back({higgs, { 199, 158, 6, 283}, {}});
ev.incoming[0] = {gluon, { 0, 0, -850, 850}, {}};
ev.incoming[1] = {gluon, { 0, 0, 702, 702}, {}};
return ev;
default:
ev.outgoing.push_back({gluon, {-350, -112, -280, 462}, {}});
ev.outgoing.push_back({gluon, { 347, 266, -322, 543}, {}});
ev.outgoing.push_back({gluon, {-349, -314, -38, 471}, {}});
ev.outgoing.push_back({gluon, { 349, 348, 12, 493}, {}});
ev.outgoing.push_back({gluon, {-342, -54, 23, 347}, {}});
ev.outgoing.push_back({gluon, { 345, -134, 138, 395}, {}});
ev.incoming[0] = {gluon, { 0, 0, -1589, 1589}, {}};
ev.incoming[1] = {gluon, { 0, 0, 1122, 1122}, {}};
return ev;
}
}
if(njet == 7){
switch(pos_boson){
case -1: // jet idx: -1 0 1 2 3 4 5
ev.outgoing.push_back({gluon, { -15, -18, -54, 59}, {}});
ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}});
ev.outgoing.push_back({gluon, { 23, -100, -64, 121}, {}});
ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}});
ev.outgoing.push_back({gluon, { -5, -92, -12, 93}, {}});
ev.outgoing.push_back({gluon, { -48, -76, -2, 90}, {}});
ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
ev.incoming[0] = {gluon, { 0, 0, -439, 439}, {}};
ev.incoming[1] = {gluon, { 0, 0, 273, 273}, {}};
return ev;
case -2: // jet idx: 0 1 2 3 4 -1 -1
ev.outgoing.push_back({gluon, { -5, -86, -82, 119}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.outgoing.push_back({gluon, { -48, -14, 20, 54}, {}});
ev.outgoing.push_back({gluon, { 23, -50, 26, 61}, {}});
ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
ev.outgoing.push_back({gluon, { -15, -18, 54, 59}, {}});
ev.outgoing.push_back({gluon, { -11, -20, 88, 91}, {}});
ev.incoming[0] = {gluon, { 0, 0, -215, 215}, {}};
ev.incoming[1] = {gluon, { 0, 0, 397, 397}, {}};
return ev;
case -3: // jet idx: 0 0 1 2 2 3 4
// jet pt fraction: 0.6 0.38 1 0.49 0.51 1 1
ev.outgoing.push_back({gluon, { 23, -94, -62, 1.2e+02}, {}});
ev.outgoing.push_back({gluon, { -5, -62, -34, 71}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 1.2e+02}, {}});
ev.outgoing.push_back({gluon, { -12, 95, 56, 1.1e+02}, {}});
ev.outgoing.push_back({gluon, { -11, 98, 70, 1.2e+02}, {}});
ev.outgoing.push_back({gluon, { -48, -1e+02, 82, 1.4e+02}, {}});
ev.outgoing.push_back({gluon, { -15, -30, 78, 85}, {}});
ev.incoming[0] = {gluon, { 0, 0, -2.7e+02, 2.7e+02}, {}};
ev.incoming[1] = {gluon, { 0, 0, 4.8e+02, 4.8e+02}, {}};
return ev;
case -4: // jet idx: 0 1 1 2 3 4 4
// jet pt fraction: 1 0.51 0.49 1 1 0.25 0.75
ev.outgoing.push_back({gluon, { -5, -88, -64, 109}, {}});
ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}});
ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.outgoing.push_back({gluon, { 23, -70, 22, 77}, {}});
ev.outgoing.push_back({gluon, { -15, -32, 16, 39}, {}});
ev.outgoing.push_back({gluon, { -48, -96, 75, 131}, {}});
ev.incoming[0] = {gluon, { 0, 0, -381, 381}, {}};
ev.incoming[1] = {gluon, { 0, 0, 324, 324}, {}};
return ev;
case -5: // jet idx: 0 1 -1 -1 2 3 4
ev.outgoing.push_back({gluon, { -15, -26, -62, 69}, {}});
ev.outgoing.push_back({gluon, { -48, -60, -54, 94}, {}});
ev.outgoing.push_back({gluon, { 23, 10, -14, 29}, {}});
ev.outgoing.push_back({gluon, { -5, -20, 0, 21}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.outgoing.push_back({gluon, { -11, -92, 40, 101}, {}});
ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
ev.incoming[0] = {gluon, { 0, 0, -278, 278}, {}};
ev.incoming[1] = {gluon, { 0, 0, 264, 264}, {}};
return ev;
case -6: // jet idx: 0 1 1 2 -1 2 3
// jet pt fraction: 1 0.63 0.36 0.49 1 0.51 1
ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}});
ev.outgoing.push_back({gluon, { -48, -100, 26, 114}, {}});
ev.outgoing.push_back({gluon, { -15, -62, 26, 69}, {}});
ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
ev.outgoing.push_back({gluon, { -5, -28, 20, 35}, {}});
ev.outgoing.push_back({gluon, { -11, 98, 70, 121}, {}});
ev.outgoing.push_back({gluon, { 23, -96, 92, 135}, {}});
ev.incoming[0] = {gluon, { 0, 0, -216, 216}, {}};
ev.incoming[1] = {gluon, { 0, 0, 486, 486}, {}};
return ev;
case -7: // jet idx: 0 1 2 2 3 3 4
// jet pt fraction: 1 1 0.51 0.49 0.18 0.82 1
ev.outgoing.push_back({gluon, { -15, -80, -100, 129}, {}});
ev.outgoing.push_back({gluon, { 23, -96, -92, 135}, {}});
ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}});
ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
ev.outgoing.push_back({gluon, { -5, -22, -10, 25}, {}});
ev.outgoing.push_back({gluon, { -48, -88, -31, 105}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.incoming[0] = {gluon, { 0, 0, -541, 541}, {}};
ev.incoming[1] = {gluon, { 0, 0, 202, 202}, {}};
return ev;
case -8: // jet idx: 0 1 2 2 2 3 4
// jet pt fraction: 1 1 0.21 0.37 0.41 1 1
ev.outgoing.push_back({gluon, { -48, -44, -62, 90}, {}});
ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
ev.outgoing.push_back({gluon, { -5, -50, -22, 55}, {}});
ev.outgoing.push_back({gluon, { 23, -90, -34, 99}, {}});
ev.outgoing.push_back({gluon, { -15, -100, -28, 105}, {}});
ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}});
ev.outgoing.push_back({gluon, { -11, 96, 76, 123}, {}});
ev.incoming[0] = {gluon, { 0, 0, -423, 423}, {}};
ev.incoming[1] = {gluon, { 0, 0, 277, 277}, {}};
return ev;
case -9: // jet idx: 0 1 2 1 3 0 4
// jet pt fraction: 0.72 0.51 1 0.49 1 0.28 1
ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}});
ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
ev.outgoing.push_back({gluon, { 23, -76, -40, 89}, {}});
ev.outgoing.push_back({gluon, { -11, 92, -40, 101}, {}});
ev.outgoing.push_back({gluon, { -48, -68, -34, 90}, {}});
ev.outgoing.push_back({gluon, { -5, -38, -14, 41}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.incoming[0] = {gluon, { 0, 0, -446, 446}, {}};
ev.incoming[1] = {gluon, { 0, 0, 220, 220}, {}};
return ev;
case -10: // jet idx: 0 1 3 2 4 3 1
// jet pt fraction: 1 0.33 0.51 1 1 0.49 0.67
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.outgoing.push_back({gluon, { -5, -48, 16, 51}, {}});
ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
ev.outgoing.push_back({gluon, { 23, -76, 52, 95}, {}});
ev.outgoing.push_back({gluon, { -48, -60, 54, 94}, {}});
ev.outgoing.push_back({gluon, { -11, 92, 68, 115}, {}});
ev.outgoing.push_back({gluon, { -15, -96, 72, 121}, {}});
ev.incoming[0] = {gluon, { 0, 0, -183, 183}, {}};
ev.incoming[1] = {gluon, { 0, 0, 521, 521}, {}};
return ev;
case -11: // jet idx: 0 1 2 3 4 -1 5
// jet pt fraction: 1 1 1 1 1 1 1
ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}});
ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}});
ev.outgoing.push_back({gluon, { 23, -90, -2, 93}, {}});
ev.outgoing.push_back({gluon, { -48, -76, 2, 90}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.outgoing.push_back({gluon, { -5, -22, 10, 25}, {}});
ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
ev.incoming[0] = {gluon, { 0, 0, -360, 360}, {}};
ev.incoming[1] = {gluon, { 0, 0, 314, 314}, {}};
return ev;
case -12: // jet idx: 0 1 -1 2 3 4 3
// jet pt fraction: 1 1 1 1 0.35 1 0.65
ev.outgoing.push_back({gluon, { 23, -94, -62, 115}, {}});
ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
ev.outgoing.push_back({gluon, { -5, -28, 4, 29}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.outgoing.push_back({gluon, { -15, -58, 34, 69}, {}});
ev.outgoing.push_back({gluon, { -11, 92, 68, 115}, {}});
ev.outgoing.push_back({gluon, { -48, -100, 82, 138}, {}});
ev.incoming[0] = {gluon, { 0, 0, -302, 302}, {}};
ev.incoming[1] = {gluon, { 0, 0, 392, 392}, {}};
return ev;
case -13: // jet idx: 0 1 2 3 3 4 2
// jet pt fraction: 1 1 0.5 0.35 0.65 1 0.5
ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
ev.outgoing.push_back({gluon, { -5, -28, 20, 35}, {}});
ev.outgoing.push_back({gluon, { -48, -96, 75, 131}, {}});
ev.outgoing.push_back({gluon, { 23, -62, 50, 83}, {}});
ev.outgoing.push_back({gluon, { -11, 96, 76, 123}, {}});
ev.incoming[0] = {gluon, { 0, 0, -241, 241}, {}};
ev.incoming[1] = {gluon, { 0, 0, 476, 476}, {}};
return ev;
case -14: // jet idx: 0 1 2 3 3 4 2
// jet pt fraction: 1 1 0.52 0.35 0.65 1 0.48
ev.outgoing.push_back({gluon, { 23, -94, -62, 115}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
ev.outgoing.push_back({gluon, { -15, -96, 72, 121}, {}});
ev.outgoing.push_back({gluon, { -5, -42, 38, 57}, {}});
ev.outgoing.push_back({gluon, { -48, -44, 62, 90}, {}});
ev.outgoing.push_back({gluon, { -11, 88, 88, 125}, {}});
ev.incoming[0] = {gluon, { 0, 0, -231, 231}, {}};
ev.incoming[1] = {gluon, { 0, 0, 505, 505}, {}};
return ev;
case -15: // jet idx: 0 -1 1 0 2 3 4
// jet pt fraction: 0.51 1 1 0.49 1 1 1
ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}});
ev.outgoing.push_back({gluon, { -5, -16, -12, 21}, {}});
ev.outgoing.push_back({gluon, { 23, -94, -62, 115}, {}});
ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
ev.outgoing.push_back({gluon, { -48, -76, 70, 114}, {}});
ev.outgoing.push_back({gluon, { -15, -100, 80, 129}, {}});
ev.incoming[0] = {gluon, { 0, 0, -379, 379}, {}};
ev.incoming[1] = {gluon, { 0, 0, 349, 349}, {}};
return ev;
}
}
throw HEJ::unknown_option{"unknown process"};
}
HEJ::Event::EventData parse_configuration(
std::array<std::string,2> const & in, std::vector<std::string> const & out,
int const overwrite_boson
){
auto boson = std::find_if(out.cbegin(), out.cend(),
[](std::string id){ return !HEJ::is_parton(HEJ::to_ParticleID(id)); });
int const pos_boson = (overwrite_boson!=0)?overwrite_boson:
((boson == out.cend())?-1:std::distance(out.cbegin(), boson));
std::size_t njets = out.size();
if( (overwrite_boson == 0) && boson != out.cend()) --njets;
HEJ::Event::EventData ev{get_process(njets, pos_boson)};
ASSERT((pos_boson<0) || (ev.outgoing[pos_boson].type == HEJ::ParticleID::higgs));
for(std::size_t i=0; i<out.size(); ++i){
ev.outgoing[i].type = HEJ::to_ParticleID(out[i]);
// decay W
if( std::abs(ev.outgoing[i].type) == HEJ::ParticleID::Wp )
ev.decays[i]=decay_W(ev.outgoing[i]);
// decay Z
if( ev.outgoing[i].type == HEJ::ParticleID::Z_photon_mix )
ev.decays[i]=decay_Z(ev.outgoing[i]);
}
for(std::size_t i=0; i<in.size(); ++i){
ev.incoming[i].type = HEJ::to_ParticleID(in[i]);
}
shuffle_particles(ev);
return ev;
}
HEJ::Event::EventData rapidity_order_ps(
std::array<std::string,2> const & in,
std::vector<std::string> const & out,
bool reconstruct /* = false */,
std::unordered_map< size_t, std::vector<std::string> > decays /* = {} */
) {
using namespace HEJ::pid;
HEJ::Event::EventData evd;
const size_t n_out {out.size()};
// Parameters
double pT {40.};
std::array<double,2> cs_phi {1, 0};
std::array<double,2> cs_dphi {cos(2.*M_PI/n_out), sin(2.*M_PI/n_out)};
double y {-4.5}, dy {std::max(-2.*y/(n_out-1), 0.5)};
double sum_pE {0}, sum_pz {0};
// Outgoing
evd.outgoing.reserve(n_out);
for(size_t i = 0; i < n_out; ++i) {
auto pid = HEJ::to_ParticleID(out[i]);
// Kinematics
double mT2 = pT*pT;
if(HEJ::is_AWZH_boson(pid)) {
auto const ref = ew_parameters.prop(pid);
mT2 += ref.mass*ref.mass;
}
double pz {sqrt(mT2)*sinh(y)}, pE {sqrt(mT2)*cosh(y)};
evd.outgoing.push_back({pid, {pT*cs_phi[0], pT*cs_phi[1], pz, pE}, {}});
sum_pE+=pE; sum_pz+=pz;
// Update cos,sin, y
cs_phi = {cs_phi[0] * cs_dphi[0] - cs_phi[1] * cs_dphi[1],
cs_phi[1] * cs_dphi[0] + cs_phi[0] * cs_dphi[1]};
y += dy;
}
// Specified decays
evd.decays.reserve(decays.size());
for(auto const & decay : decays) {
auto const parent = evd.outgoing.at(decay.first);
std::vector<HEJ::Particle> progeny = decay_kinematics(parent);
for(std::size_t i = 0; i < decay.second.size(); ++i) {
progeny[i].type = HEJ::to_ParticleID(decay.second[i]);
}
evd.decays.emplace(decay.first, std::move(progeny));
}
// Reconstruct decays
if(reconstruct) { evd.reconstruct_intermediate(ew_parameters); }
// Incoming
double p0 {0.5*(sum_pE + sum_pz)};
evd.incoming[0] = {HEJ::to_ParticleID(in[0]),{0,0,-p0,p0},{}};
double p1 {0.5*(sum_pE - sum_pz)};
evd.incoming[1] = {HEJ::to_ParticleID(in[1]),{0,0,+p1,p1},{}};
return evd;
}
namespace {
static std::mt19937_64 ran{0};
}
void shuffle_particles(HEJ::Event::EventData & ev) {
// incoming
std::shuffle(ev.incoming.begin(), ev.incoming.end(), ran);
// outgoing (through index)
auto old_outgoing = std::move(ev).outgoing;
std::vector<std::size_t> idx(old_outgoing.size());
std::iota(idx.begin(), idx.end(), 0);
std::shuffle(idx.begin(), idx.end(), ran);
ev.outgoing.clear();
ev.outgoing.reserve(old_outgoing.size());
for(std::size_t i: idx) {
ev.outgoing.emplace_back(std::move(old_outgoing[i]));
}
// find decays again
if(!ev.decays.empty()){
auto old_decays = std::move(ev).decays;
ev.decays.clear();
for(std::size_t i=0; i<idx.size(); ++i) {
auto decay = old_decays.find(idx[i]);
if(decay != old_decays.end())
ev.decays.emplace(i, std::move(decay->second));
}
for(auto & decay: ev.decays){
std::shuffle(decay.second.begin(), decay.second.end(), ran);
}
}
}
std::vector<HEJ::Particle> decay_kinematics( HEJ::Particle const & parent ) {
std::vector<HEJ::Particle> decay_products(2);
const double E = parent.m()/2;
const double theta = 2.*M_PI*ran()/static_cast<double>(ran.max());
const double cos_phi = 2.*ran()/static_cast<double>(ran.max())-1.;
const double sin_phi = std::sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
const double px = E*std::cos(theta)*sin_phi;
const double py = E*std::sin(theta)*sin_phi;
const double pz = E*cos_phi;
decay_products[0].p.reset(px, py, pz, E);
decay_products[1].p.reset(-px, -py, -pz, E);
for(auto & particle: decay_products) particle.p.boost(parent.p);
return decay_products;
}
std::vector<HEJ::Particle> decay_W( HEJ::Particle const & parent ){
if(parent.m() == 0.) // we can't decay massless partons
return {};
std::array<HEJ::ParticleID, 2> decays;
if(parent.type==HEJ::ParticleID::Wp){
// order matters: first particle, second anti
decays[0] = HEJ::ParticleID::nu_e;
decays[1] = HEJ::ParticleID::e_bar;
} else {
// this function is for testing: we don't check that parent==W boson
decays[0] = HEJ::ParticleID::e;
decays[1] = HEJ::ParticleID::nu_e_bar;
}
std::vector<HEJ::Particle> decay_products = decay_kinematics(parent);
for(std::size_t i = 0; i < decay_products.size(); ++i){
decay_products[i].type = decays[i];
}
return decay_products;
}
std::vector<HEJ::Particle> decay_Z( HEJ::Particle const & parent ){
if(parent.m() == 0.) // we can't decay massless partons
return {};
std::array<HEJ::ParticleID, 2> decays;
// order matters: first particle, second anti
decays[0] = HEJ::ParticleID::electron;
decays[1] = HEJ::ParticleID::positron;
std::vector<HEJ::Particle> decay_products = decay_kinematics(parent);
for(std::size_t i = 0; i < decay_products.size(); ++i){
decay_products[i].type = decays[i];
}
return decay_products;
}
diff --git a/t/test_classify_ref.cc b/t/test_classify_ref.cc
index bdb6bf4..d5256a3 100644
--- a/t/test_classify_ref.cc
+++ b/t/test_classify_ref.cc
@@ -1,85 +1,86 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "hej_test.hh"
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>
#include <fastjet/JetDefinition.hh>
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReader.hh"
+#include "HEJ/EWConstants.hh"
namespace {
// this is deliberately chosen bigger than in the generation,
// to cluster multiple partons in one jet
constexpr double min_jet_pt = 40.;
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.6};
const HEJ::ParticleProperties Wprop{80.385, 2.085};
const HEJ::ParticleProperties Zprop{91.187, 2.495};
const HEJ::ParticleProperties Hprop{125, 0.004165};
constexpr double vev = 246.2196508;
const HEJ::EWConstants ew_parameters{vev, Wprop, Zprop, Hprop};
}
int main(int argn, char** argv) {
if(argn != 3 && argn != 4){
std::cerr << "Usage: " << argv[0]
<< " reference_classification input_file.lhe\n";
return EXIT_FAILURE;
}
bool OUTPUT_MODE = false;
if(argn == 4 && std::string("OUTPUT")==std::string(argv[3]))
OUTPUT_MODE = true;
std::fstream ref_file;
if ( OUTPUT_MODE ) {
std::cout << "_______________________USING OUTPUT MODE!_______________________" << std::endl;
ref_file.open(argv[1], std::fstream::out);
} else {
ref_file.open(argv[1], std::fstream::in);
}
auto reader{ HEJ::make_reader(argv[2]) };
std::size_t nevent{0};
while(reader->read_event()){
++nevent;
// We don't need to test forever, the first "few" are enough
if(nevent>4000) break;
HEJ::Event::EventData data{ reader->hepeup() };
shuffle_particles(data);
data.reconstruct_intermediate(ew_parameters);
const HEJ::Event ev{
data.cluster(
jet_def, min_jet_pt
)
};
if ( OUTPUT_MODE ) {
ref_file << ev.type() << std::endl;
} else {
std::string line;
if(!std::getline(ref_file,line)) break;
const auto expected{static_cast<HEJ::event_type::EventType>(std::stoi(line))};
if(ev.type() != expected){
std::cerr << "wrong classification of event " << nevent << ":\n" << ev
<< "classified as " << name(ev.type())
<< ", expected " << name(expected) << "\nJet indices: ";
for(auto const & idx: ev.particle_jet_indices())
std::cerr << idx << " ";
std::cerr << "\n";
return EXIT_FAILURE;
}
}
}
return EXIT_SUCCESS;
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:01 PM (2 h, 18 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486211
Default Alt Text
(92 KB)

Event Timeline