diff --git a/FixedOrderGen/t/W_nj_classify.cc b/FixedOrderGen/t/W_nj_classify.cc
index 906d6f7..db1609e 100644
--- a/FixedOrderGen/t/W_nj_classify.cc
+++ b/FixedOrderGen/t/W_nj_classify.cc
@@ -1,210 +1,207 @@
/**
* \brief check that the PSP generates the all W+jet subleading processes
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <string>
#include <unordered_map>
#include <vector>
#include "HEJ/Event.hh"
#include "HEJ/event_types.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/PDF.hh"
#include "HEJ/PDG_codes.hh"
#include "fastjet/JetDefinition.hh"
#include "Decay.hh"
#include "JetParameters.hh"
#include "PhaseSpacePoint.hh"
#include "Process.hh"
#include "Status.hh"
#include "Subleading.hh"
namespace {
using namespace HEJFOG;
using namespace HEJ;
void print_psp(PhaseSpacePoint const & psp){
std::cerr << "Process:\n"
<< psp.incoming()[0].type << " + "<< psp.incoming()[1].type << " -> ";
for(auto const & out: psp.outgoing()){
std::cerr << out.type << " ";
}
std::cerr << "\n";
}
void bail_out(PhaseSpacePoint const & psp, std::string msg){
print_psp(psp);
throw std::logic_error{msg};
}
}
int main(){
constexpr std::size_t n_psp_base = 10375;
const JetParameters jet_para{
fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30, 5, 30};
PDF pdf(11000, pid::proton, pid::proton);
constexpr double E_cms = 13000.;
constexpr double subl_change = 0.8;
const ParticlesDecayMap boson_decays{
{pid::Wp, {Decay{ {pid::e_bar, pid::nu_e}, 1.} }},
{pid::Wm, {Decay{ {pid::e, pid::nu_e_bar}, 1.} }}
};
const EWConstants ew_constants{246.2196508,
ParticleProperties{80.385, 2.085},
ParticleProperties{91.187, 2.495},
ParticleProperties{125, 0.004165}
};
HEJ::Mixmax ran{};
auto subl_channels = Subleading::all;
std::vector<event_type::EventType> allowed_types{event_type::FKL,
event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf};
std::cout << "Wp3j" << std::endl;
// Wp3j
Process proc {{pid::proton,pid::proton}, 3, pid::Wp, {}};
std::size_t n_psp = n_psp_base;
#if !defined(__clang__) && defined(__GNUC__) && (__GNUC__ < 6)
// gcc version < 6 explicitly needs hash function for enum
// see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=60970
std::unordered_map<event_type::EventType, std::size_t, std::hash<std::size_t>> type_counter;
#else
std::unordered_map<event_type::EventType, std::size_t> type_counter;
#endif
for( std::size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_decays, ew_constants, ran};
if(psp.status()==good){
const Event ev{ to_EventData(psp).cluster(jet_para.def, jet_para.min_pt) };
++type_counter[ev.type()];
if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
== allowed_types.cend()) {
- bail_out(psp, "Found not allowed event of type "
- +std::string(event_type::name(ev.type())));
+ bail_out(psp, "Found not allowed event of type " + name(ev.type()) );
}
} else { // bad process -> try again
++n_psp;
}
}
std::cout << "Wp+3j: Took " << n_psp << " to generate "
<< n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
std::cout << "States by classification:\n";
for(auto const & entry: type_counter){
const double fraction = static_cast<double>(entry.second)/n_psp_base;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(25)
- << (event_type::name(entry.first) + std::string(":"))
+ << (name(entry.first) + std::string(":"))
<< entry.second << " (" << percent << "%)\n";
}
for(auto const & t: allowed_types){
if(type_counter[t] < 0.05 * n_psp_base){
- std::cerr << "Less than 5% of the events are of type " << event_type::name(t) << std::endl;
+ std::cerr << "Less than 5% of the events are of type " << name(t) << std::endl;
return EXIT_FAILURE;
}
}
// Wm3j - only uno
proc = Process{{pid::proton,pid::proton}, 3, pid::Wm, {}};
n_psp = n_psp_base;
subl_channels = Subleading::uno;
allowed_types = {event_type::FKL, event_type::unob, event_type::unof};
type_counter.clear();
for( std::size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_decays, ew_constants, ran};
if(psp.status()==good){
const Event ev{ to_EventData(psp).cluster(jet_para.def, jet_para.min_pt) };
++type_counter[ev.type()];
if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
== allowed_types.cend()) {
- bail_out(psp, "Found not allowed event of type "
- +std::string(event_type::name(ev.type())));
+ bail_out(psp, "Found not allowed event of type " + name(ev.type()) );
}
} else { // bad process -> try again
++n_psp;
}
}
std::cout << "Wm+3j (only uno): Took " << n_psp << " to generate "
<< n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
std::cout << "States by classification:\n";
for(auto const & entry: type_counter){
const double fraction = static_cast<double>(entry.second)/n_psp_base;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(25)
- << (event_type::name(entry.first) + std::string(":"))
+ << (name(entry.first) + std::string(":"))
<< entry.second << " (" << percent << "%)\n";
}
for(auto const & t: allowed_types){
if(type_counter[t] < 0.05 * n_psp_base){
- std::cerr << "Less than 5% of the events are of type " << event_type::name(t) << std::endl;
+ std::cerr << "Less than 5% of the events are of type " << name(t) << std::endl;
return EXIT_FAILURE;
}
}
// Wm4j
proc = Process{{pid::proton,pid::proton}, 4, pid::Wm, {}};
n_psp = n_psp_base;
subl_channels = Subleading::all;
allowed_types = {event_type::FKL,
event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf,
event_type::qqxmid};
type_counter.clear();
for( std::size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_decays, ew_constants, ran};
if(psp.status()==good){
const Event ev{ to_EventData(psp).cluster(jet_para.def, jet_para.min_pt)};
++type_counter[ev.type()];
if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
== allowed_types.cend()) {
- bail_out(psp, "Found not allowed event of type "
- +std::string(event_type::name(ev.type())));
+ bail_out(psp, "Found not allowed event of type " + name(ev.type()) );
}
} else { // bad process -> try again
++n_psp;
}
}
std::cout << "Wm+4j: Took " << n_psp << " to generate "
<< n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
std::cout << "States by classification:\n";
for(auto const & entry: type_counter){
const double fraction = static_cast<double>(entry.second)/n_psp_base;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(25)
- << (event_type::name(entry.first) + std::string(":"))
+ << (name(entry.first) + std::string(":"))
<< entry.second << " (" << percent << "%)\n";
}
for(auto const & t: allowed_types){
if(type_counter[t] < 0.03 * n_psp_base){
- std::cerr << "Less than 3% of the events are of type " << event_type::name(t) << std::endl;
+ std::cerr << "Less than 3% of the events are of type " << name(t) << std::endl;
return EXIT_FAILURE;
}
}
std::cout << "All processes passed." << std::endl;
return EXIT_SUCCESS;
}
diff --git a/src/CrossSectionAccumulator.cc b/src/CrossSectionAccumulator.cc
index 4855f4b..cb3811d 100644
--- a/src/CrossSectionAccumulator.cc
+++ b/src/CrossSectionAccumulator.cc
@@ -1,67 +1,67 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/CrossSectionAccumulator.hh"
#include <iomanip>
#include <cmath>
#include <string>
#include <utility>
#include "HEJ/Event.hh"
#include "HEJ/Parameters.hh"
namespace HEJ {
void CrossSectionAccumulator::fill( double const wt, double const err,
event_type::EventType const type
){
total_.value += wt;
total_.error += err;
auto & entry = xs_[type];
entry.value += wt;
entry.error += err;
}
void CrossSectionAccumulator::fill(
double const wt, event_type::EventType const type
){
fill(wt, wt*wt, type);
}
void CrossSectionAccumulator::fill(Event const & ev){
fill(ev.central().weight, ev.type());
}
void CrossSectionAccumulator::fill_correlated(
double const sum, double const sum2, event_type::EventType const type
){
fill(sum, sum*sum+sum2, type);
}
void CrossSectionAccumulator::fill_correlated(std::vector<Event> const & evts){
if(evts.empty()) return;
fill_correlated(evts.cbegin(), evts.cend());
}
std::ostream& operator<<(std::ostream& os, const CrossSectionAccumulator& xs){
const std::streamsize orig_prec = os.precision();
os << std::scientific << std::setprecision(3)
<< " " << std::left << std::setw(25)
<< "Cross section: " << xs.total().value
<< " +- " << std::sqrt(xs.total().error) << " (pb)\n";
for(auto const & xs_type: xs) {
os << " " << std::left << std::scientific <<std::setw(25)
- << (event_type::name(xs_type.first) + std::string(": "));
+ << (name(xs_type.first) + std::string(": "));
os << xs_type.second.value << " +- "
<< std::sqrt(xs_type.second.error) << " (pb) "
<< std::fixed << std::setprecision(3)
<< "[" <<std::setw(6) <<std::right
<< ((xs_type.second.value)/xs.total().value)*100 << "%]"
<< std::endl;
}
os << std::defaultfloat;
os.precision(orig_prec);
return os;
}
}
diff --git a/src/Event.cc b/src/Event.cc
index 9770986..f9ef5f8 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,1138 +1,1138 @@
/**
* \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/exceptions.hh"
#include "HEJ/optional.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/RNG.hh"
namespace HEJ{
namespace {
using std::size_t;
constexpr int status_in = -1;
constexpr int status_decayed = 2;
constexpr int status_out = 1;
//! 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]);
}
/// @name helper functions to determine event type
//@{
/**
* \brief check if final state valid for HEJ
*
* check if there is at most one photon, W, H, Z in the final state
* and all the rest are quarks or gluons
*/
bool final_state_ok(Event const & ev){
std::vector<Particle> const & outgoing = ev.outgoing();
if(ev.decays().size() > 1) // at most one decay
return false;
bool has_AWZH_boson = false;
for( size_t i=0; i<outgoing.size(); ++i ){
auto const & out{ outgoing[i] };
if(is_AWZH_boson(out.type)){
// at most one boson
if(has_AWZH_boson) return false;
has_AWZH_boson = true;
// valid decay for W
if(std::abs(out.type) == ParticleID::Wp){
// exactly 1 decay of W
if( ev.decays().size() != 1 || ev.decays().cbegin()->first != i )
return false;
if( !valid_W_decay(out.type>0?+1:-1, ev.decays().cbegin()->second) )
return false;
}
}
else if(! is_parton(out.type)) 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) return non_resummable; // multi boson
switch (bosons[0].type) {
case ParticleID::Wp:
case ParticleID::Wm:
return FKL | unob | unof | qqxexb | qqxexf | qqxmid;
case ParticleID::h:
return FKL | unob | unof;
default:
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)<=6 || in==pid::gluon) return (in==out*qqxCurrent);
else 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;
assert(boson.size() < 2);
// 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;
assert(boson.size() < 2);
// 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 boson{ filter_AWZH_bosons(ev.outgoing()) };
// we only allow one boson through final_state_ok
assert(boson.size()<=1);
// keep track of potential W couplings, at the end the sum should be 0
int remaining_Wp = 0;
int remaining_Wm = 0;
if(!boson.empty() && std::abs(boson.front().type) == ParticleID::Wp ){
if(boson.front().type>0) ++remaining_Wp;
else ++remaining_Wm;
}
int W_change = 0;
size_t final_type = ~(no_2_jets | bad_final_state);
// check forward impact factor
final_type &= possible_impact_factors(
in.front().type,
begin_out, end_out.base(),
W_change, boson, 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;
W_change = 0;
// check backward impact factor
final_type &= possible_impact_factors(
in.back().type,
end_out, std::make_reverse_iterator(begin_out),
W_change, boson, 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;
W_change = 0;
// check central emissions
final_type &= possible_central(
begin_out, end_out.base(), W_change, boson );
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(boson)) != 0 )
return non_resummable;
return static_cast<EventType>(final_type);
}
//@}
Particle extract_particle(LHEF::HEPEUP const & hepeup, size_t i){
const ParticleID id = static_cast<ParticleID>(hepeup.IDUP[i]);
const fastjet::PseudoJet momentum{
hepeup.PUP[i][0], hepeup.PUP[i][1],
hepeup.PUP[i][2], hepeup.PUP[i][3]
};
if(is_parton(id))
return Particle{ id, std::move(momentum), hepeup.ICOLUP[i] };
return Particle{ id, std::move(momentum), {} };
}
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 anonymous
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] == 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 o1, Particle 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 {
Particle reconstruct_boson(std::vector<Particle> const & leptons) {
Particle decayed_boson;
decayed_boson.p = leptons[0].p + leptons[1].p;
const int pidsum = leptons[0].type + leptons[1].type;
if(pidsum == +1) {
assert(is_antilepton(leptons[0]));
if(is_antineutrino(leptons[0])) {
throw not_implemented{"lepton-flavour violating final state"};
}
assert(is_neutrino(leptons[1]));
// charged antilepton + neutrino means we had a W+
decayed_boson.type = pid::Wp;
}
else if(pidsum == -1) {
assert(is_antilepton(leptons[0]));
if(is_neutrino(leptons[1])) {
throw not_implemented{"lepton-flavour violating final state"};
}
assert(is_antineutrino(leptons[0]));
// charged lepton + antineutrino means we had a W-
decayed_boson.type = pid::Wm;
}
else {
throw not_implemented{
"final state with leptons "
+ name(leptons[0].type)
+ " and "
+ name(leptons[1].type)
};
}
return decayed_boson;
}
}
void Event::EventData::reconstruct_intermediate() {
const auto 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.size() != 2) {
throw not_implemented{"Final states with one or more than two leptons"};
}
std::sort(
begin(leptons), end(leptons),
[](Particle const & p0, Particle const & p1) {
return p0.type < p1.type;
}
);
outgoing.emplace_back(reconstruct_boson(leptons));
decays.emplace(outgoing.size()-1, std::move(leptons));
}
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 anonymous
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;
return;
}
//! 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 anonymous
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] };
if( (parton.p - jet).pt()/jet.pt() > max_ext_soft_pt_fraction)
return false;
return true;
}
}
// 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;
}
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){
const std::streamsize orig_prec = os.precision();
os <<std::scientific<<std::setprecision(6) << "["
<<std::setw(13)<<std::right<< part.px() << ", "
<<std::setw(13)<<std::right<< part.py() << ", "
<<std::setw(13)<<std::right<< part.pz() << ", "
<<std::setw(13)<<std::right<< part.E() << "]"<< std::fixed;
os.precision(orig_prec);
}
void print_colour(std::ostream & os, optional<Colour> const & col){
if(!col)
os << "(no color)"; // American spelling for better alignment
else
os << "(" <<std::setw(3)<<std::right<< col->first
<< ", " <<std::setw(3)<<std::right<< col->second << ")";
}
}
std::ostream& operator<<(std::ostream & os, Event const & ev){
const std::streamsize orig_prec = os.precision();
os <<std::setprecision(4)<<std::fixed;
- os << "########## " << event_type::name(ev.type()) << " ##########" << std::endl;
+ os << "########## " << name(ev.type()) << " ##########" << std::endl;
os << "Incoming particles:\n";
for(auto const & in: ev.incoming()){
os <<std::setw(3)<< 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(3)<< out.type << ": ";
print_colour(os, out.colour);
os << " ";
print_momentum(os, out.p);
os << " => rapidity="
<<std::setw(7)<<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(7)<<std::right<< jet.rapidity() << std::endl;
}
if(ev.decays().size() > 0 ){
os << "\nDecays: " << ev.decays().size() << "\n";
for(auto const & decay: ev.decays()){
os <<std::setw(3)<< ev.outgoing()[decay.first].type
<< " (" << decay.first << ") to:\n";
for(auto const & out: decay.second){
os <<" "<<std::setw(3)<< out.type << ": ";
print_momentum(os, out.p);
os << " => rapidity="
<<std::setw(7)<<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(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)?status_decayed: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(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;
}
}
diff --git a/t/test_classify.cc b/t/test_classify.cc
index 2098e21..d295be8 100644
--- a/t/test_classify.cc
+++ b/t/test_classify.cc
@@ -1,499 +1,499 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "hej_test.hh"
#include <array>
#include <cstdlib>
#include <iostream>
#include <random>
#include <string>
#include <vector>
#include "fastjet/JetDefinition.hh"
#include "HEJ/Event.hh"
#include "HEJ/event_types.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/PDG_codes.hh"
namespace {
const fastjet::JetDefinition jet_def{fastjet::JetAlgorithm::antikt_algorithm, 0.4};
const double min_jet_pt{30.};
const std::vector<std::string> all_quarks{"-4","-1","1","2","3","4"};
const std::vector<std::string> all_partons{"g","-2","-1","1","2","3","4"};
const std::vector<std::string> all_bosons{"h", "Wp", "Wm"};
const std::vector<std::string> all_gZ{"photon", "Z"};
const std::vector<std::string> all_w{"W+", "W-"};
static std::mt19937_64 ran{0};
bool couple_quark(std::string const & boson, std::string & quark){
if(std::abs(HEJ::to_ParticleID(boson)) == HEJ::ParticleID::Wp){
auto qflav{ HEJ::to_ParticleID(quark) };
if(!HEJ::is_anyquark(qflav)) return false;
const int W_charge = HEJ::to_ParticleID(boson)>0?1:-1;
if(W_charge*qflav < 0 && !(std::abs(qflav)%2)) return false; // not anti-down
if(W_charge*qflav > 0 && (std::abs(qflav)%2)) return false; // not up
quark=std::to_string(qflav-W_charge);
}
return true;
}
bool match_expectation( HEJ::event_type::EventType expected,
std::array<std::string,2> const & in, std::vector<std::string> const & out,
int const overwrite_boson = 0
){
HEJ::Event ev{ parse_configuration(
in,out,overwrite_boson ).cluster(jet_def, min_jet_pt)};
if(ev.type() != expected){
- std::cerr << "Expected type " << HEJ::event_type::name(expected)
- << " but found " << HEJ::event_type::name(ev.type()) << "\n" << ev;
+ std::cerr << "Expected type " << name(expected)
+ << " but found " << name(ev.type()) << "\n" << ev;
auto jet_idx{ ev.particle_jet_indices() };
std::cout << "Particle Jet indices: ";
for(int const i: jet_idx)
std::cout << i << " ";
std::cout << std::endl;
return false;
}
return true;
}
//! test FKL configurations
//! if implemented==false : check processes that are not in HEJ yet
bool check_fkl( bool const implemented=true ){
using namespace HEJ;
auto const type{ implemented?event_type::FKL:event_type::non_resummable };
std::vector<std::string> bosons;
if(implemented)
bosons = all_bosons;
else {
bosons = all_gZ;
}
for(std::string const & first: all_partons) // all quark flavours
for(std::string const & last: all_partons){
for(int njet=2; njet<=6; ++njet){ // all multiplicities
if(njet==5) continue;
std::array<std::string,2> base_in{first,last};
std::vector<std::string> base_out(njet, "g");
base_out.front() = first;
base_out.back() = last;
if(implemented && !match_expectation(type, base_in, base_out))
return false;
for(auto const & boson: bosons) // any boson
for(int pos=0; pos<=njet; ++pos){ // at any position
auto in{base_in};
auto out{base_out};
// change quark flavours for W
const bool couple_idx = std::uniform_int_distribution<int>{0,1}(ran);
if(!couple_quark(boson, couple_idx?out.back():out.front()))
continue;
out.insert(out.begin()+pos, boson);
if(!match_expectation(type, in, out))
return false;
}
}
}
return true;
}
//! test unordered configurations
//! if implemented==false : check processes that are not in HEJ yet
bool check_uno( bool const implemented=true ){
using namespace HEJ;
auto const b{ implemented?event_type::unob:event_type::non_resummable };
auto const f{ implemented?event_type::unof:event_type::non_resummable };
std::vector<std::string> bosons;
if(implemented)
bosons = all_bosons;
else {
bosons = all_gZ;
}
for(std::string const & uno: all_quarks) // all quark flavours
for(std::string const & fkl: all_partons){
for(int njet=3; njet<=6; ++njet){ // all multiplicities >2
if(njet==5) continue;
for(int i=0; i<2; ++i){ // forward & backwards
std::array<std::string,2> base_in;
std::vector<std::string> base_out(njet, "g");
const int uno_pos = i?1:(njet-2);
const int fkl_pos = i?(njet-1):0;
base_in[i?0:1] = uno;
base_in[i?1:0] = fkl;
base_out[uno_pos] = uno;
base_out[fkl_pos] = fkl;
auto expectation{ i?b:f };
if( implemented
&& !match_expectation(expectation, base_in, base_out) )
return false;
for(auto const & boson: bosons){ // any boson
// at any position (higgs only inside FKL chain)
int start = 0;
int end = njet;
if(to_ParticleID(boson) == pid::higgs){
start = i?(uno_pos+1):fkl_pos;
end = i?(fkl_pos+1):uno_pos;
}
for(int pos=start; pos<=end; ++pos){
auto in{base_in};
auto out{base_out};
// change quark flavours for W
const bool couple_idx = std::uniform_int_distribution<int>{0,1}(ran);
if(!couple_quark(boson, couple_idx?out[fkl_pos]:out[uno_pos]))
continue;
out.insert(out.begin()+pos, boson);
if(!match_expectation(expectation, in, out))
return false;
}
}
}
}
}
return true;
}
//! test extremal qqx configurations
//! if implemented==false : check processes that are not in HEJ yet
bool check_extremal_qqx( bool const implemented=true ){
using namespace HEJ;
auto const b{ implemented?event_type::qqxexb:event_type::non_resummable };
auto const f{ implemented?event_type::qqxexf:event_type::non_resummable };
std::vector<std::string> bosons;
if(implemented)
bosons = all_w;
else {
bosons = all_gZ;
bosons.push_back("h");
}
for(std::string const & qqx: all_quarks) // all quark flavours
for(std::string const & fkl: all_partons){
std::string const qqx2{ std::to_string(HEJ::to_ParticleID(qqx)*-1) };
for(int njet=3; njet<=6; ++njet){ // all multiplicities >2
if(njet==5) continue;
for(int i=0; i<2; ++i){ // forward & backwards
std::array<std::string,2> base_in;
std::vector<std::string> base_out(njet, "g");
const int qqx_pos = i?0:(njet-2);
const int fkl_pos = i?(njet-1):0;
base_in[i?0:1] = "g";
base_in[i?1:0] = fkl;
base_out[fkl_pos] = fkl;
base_out[qqx_pos] = qqx;
base_out[qqx_pos+1] = qqx2;
auto expectation{ i?b:f };
if( implemented
&& !match_expectation(expectation, base_in, base_out) )
return false;
for(auto const & boson: bosons){ // all bosons
// at any position (higgs only inside FKL chain)
int start = 0;
int end = njet;
if(to_ParticleID(boson) == pid::higgs){
start = i?(qqx_pos+2):fkl_pos;
end = i?(fkl_pos+1):qqx_pos;
}
for(int pos=start; pos<=end; ++pos){
auto in{base_in};
auto out{base_out};
// change quark flavours for W
const bool couple_idx = std::uniform_int_distribution<int>{0,1}(ran);
if(couple_idx || !couple_quark(boson, out[fkl_pos]) ){
// (randomly) try couple to FKL, else fall-back to qqx
if(!couple_quark(boson, out[qqx_pos]))
couple_quark(boson, out[qqx_pos+1]);
}
out.insert(out.begin()+pos, boson);
if(!match_expectation(expectation, in, out))
return false;
}
}
}
}
// test allowed jet configurations
if( implemented){
if( !( match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -3)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -4)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -5)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -5)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -6)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -7)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -7)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -8)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -8)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -9)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -10)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -11)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -11)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -12)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -12)
))
return false;
if (fkl == "2") {
if( !( match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqx,qqx2}, -3)
&& match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -4)
&& match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqx,qqx2}, -5)
&& match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -5)
&& match_expectation(f,{"2","g"},{"1","g","Wp","g","g",qqx,qqx2}, -6)
&& match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -7)
&& match_expectation(b,{"g","2"},{qqx,qqx2,"g","g","g","Wp","1"}, -7)
&& match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqx,qqx2}, -8)
&& match_expectation(b,{"g","2"},{qqx,qqx2,"Wp","g","g","g","1"}, -8)
&& match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -9)
&& match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -10)
&& match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -11)
&& match_expectation(b,{"g","2"},{qqx,qqx2,"g","g","g","Wp","1"}, -11)
&& match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -12)
&& match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -12)
))
return false;
}
}
}
return true;
}
//! test central qqx configurations
//! if implemented==false : check processes that are not in HEJ yet
bool check_central_qqx(bool const implemented=true){
using namespace HEJ;
auto const t{ implemented?event_type::qqxmid:event_type::non_resummable };
std::vector<std::string> bosons;
if(implemented)
bosons = all_w;
else {
bosons = all_gZ;
bosons.push_back("h");
}
for(std::string const & qqx: all_quarks) // all quark flavours
for(std::string const & fkl1: all_partons)
for(std::string const & fkl2: all_partons){
std::string const qqx2{ std::to_string(HEJ::to_ParticleID(qqx)*-1) };
for(int njet=4; njet<=6; ++njet){ // all multiplicities >3
if(njet==5) continue;
for(int qqx_pos=1; qqx_pos<njet-2; ++qqx_pos){ // any qqx position
std::array<std::string,2> base_in;
std::vector<std::string> base_out(njet, "g");
base_in[0] = fkl1;
base_in[1] = fkl2;
base_out.front() = fkl1;
base_out.back() = fkl2;
base_out[qqx_pos] = qqx;
base_out[qqx_pos+1] = qqx2;
if( implemented && !match_expectation(t, base_in, base_out) )
return false;
for(auto const & boson: bosons) // any boson
for(int pos=0; pos<=njet; ++pos){ // at any position
if( to_ParticleID(boson) == pid::higgs
&& (pos==qqx_pos || pos==qqx_pos+1) )
continue;
auto in{base_in};
auto out{base_out};
// change quark flavours for W
const int couple_idx{ std::uniform_int_distribution<int>{0,2}(ran) };
// (randomly) try couple to FKL, else fall-back to qqx
if( couple_idx == 0 && couple_quark(boson, out.front()) ){}
else if( couple_idx == 1 && couple_quark(boson, out.back()) ){}
else {
if(!couple_quark(boson, out[qqx_pos]))
couple_quark(boson, out[qqx_pos+1]);
}
out.insert(out.begin()+pos, boson);
if(!match_expectation(t, in, out))
return false;
}
}
}
}
return true;
}
// this checks a (non excessive) list of non-resummable states
bool check_non_resummable(){
auto type{ HEJ::event_type::non_resummable};
return
// 2j - crossing lines
match_expectation(type, {"g","2"}, {"2","g"})
&& match_expectation(type, {"-1","g"}, {"g","-1"})
&& match_expectation(type, {"1","-1"}, {"-1","1"})
&& match_expectation(type, {"g","2"}, {"2","g","h"})
&& match_expectation(type, {"1","2"}, {"2","h","1"})
&& match_expectation(type, {"1","-1"}, {"h","-1","1"})
&& match_expectation(type, {"g","2"}, {"Wp","1","g"})
&& match_expectation(type, {"1","-1"}, {"-2","Wp","1"})
&& match_expectation(type, {"4","g"}, {"g","3","Wp"})
&& match_expectation(type, {"1","-2"}, {"-1","Wm","1"})
&& match_expectation(type, {"g","3"}, {"4","g","Wm"})
&& match_expectation(type, {"1","3"}, {"Wm","4","1"})
// 2j - qqx
&& match_expectation(type, {"g","g"}, {"1","-1"})
&& match_expectation(type, {"g","g"}, {"-2","2","h"})
&& match_expectation(type, {"g","g"}, {"-4","Wp","3"})
&& match_expectation(type, {"g","g"}, {"Wm","-1","2"})
// 3j - crossing lines
&& match_expectation(type, {"g","4"}, {"4","g","g"})
&& match_expectation(type, {"-1","g"}, {"g","g","-1"})
&& match_expectation(type, {"1","3"}, {"3","g","1"})
&& match_expectation(type, {"-2","2"}, {"2","g","-2","h"})
&& match_expectation(type, {"-3","g"}, {"g","g","Wp","-4"})
&& match_expectation(type, {"1","-2"}, {"Wm","-1","g","1"})
&& match_expectation(type, {"-1","g"}, {"1","-1","-1"})
// higgs inside uno
&& match_expectation(type, {"-1","g"}, {"g","h","-1","g"})
&& match_expectation(type, {"-1","1"}, {"g","h","-1","1"})
&& match_expectation(type, {"g","2"}, {"g","2","h","g"})
&& match_expectation(type, {"-1","1"}, {"-1","1","h","g"})
// higgs outside uno
&& match_expectation(type, {"-1","g"}, {"h","g","-1","g"})
&& match_expectation(type, {"-1","1"}, {"-1","1","g","h"})
// higgs inside qqx
&& match_expectation(type, {"g","g"}, {"-1","h","1","g","g"})
&& match_expectation(type, {"g","g"}, {"g","-1","h","1","g"})
&& match_expectation(type, {"g","g"}, {"g","g","2","h","-2"})
// higgs outside qqx
&& match_expectation(type, {"g","g"}, {"h","-1","1","g","g"})
&& match_expectation(type, {"g","g"}, {"g","g","2","-2","h"})
// 4j - two uno
&& match_expectation(type, {"-2","2"}, {"g","-2","2","g"})
&& match_expectation(type, {"1","3"}, {"g","1","h","3","g"})
&& match_expectation(type, {"1","2"}, {"g","1","3","Wp","g"})
&& match_expectation(type, {"1","-2"}, {"g","Wm","1","-1","g"})
// 4j - two gluon outside
&& match_expectation(type, {"g","4"}, {"g","4","g","g"})
&& match_expectation(type, {"1","3"}, {"1","3","h","g","g"})
&& match_expectation(type, {"1","2"}, {"1","3","g","Wp","g"})
&& match_expectation(type, {"1","-2"}, {"1","Wm","-1","g","g"})
&& match_expectation(type, {"-1","g"}, {"g","g","-1","g"})
&& match_expectation(type, {"1","3"}, {"g","g","1","3","h"})
&& match_expectation(type, {"1","2"}, {"g","g","1","Wp","3"})
&& match_expectation(type, {"1","-2"}, {"Wm","g","g","1","-1"})
// 4j - ggx+uno
&& match_expectation(type, {"g","4"}, {"1","-1","4","g"})
&& match_expectation(type, {"2","g"}, {"g","2","-3","3"})
&& match_expectation(type, {"g","4"}, {"1","-1","h","4","g"})
&& match_expectation(type, {"2","g"}, {"g","2","-3","3","h"})
&& match_expectation(type, {"g","4"}, {"Wp","1","-1","3","g"})
&& match_expectation(type, {"2","g"}, {"g","2","-4","Wp","3"})
&& match_expectation(type, {"g","4"}, {"2","Wm","-1","4","g"})
&& match_expectation(type, {"2","g"}, {"g","2","Wp","-3","4"})
// 3j - crossing+uno
&& match_expectation(type, {"1","4"}, {"g","4","1"})
&& match_expectation(type, {"1","4"}, {"4","1","g"})
&& match_expectation(type, {"1","4"}, {"g","h","4","1"})
&& match_expectation(type, {"-1","-3"},{"Wm","g","-4","-1"})
&& match_expectation(type, {"1","4"}, {"3","1","Wp","g"})
&& match_expectation(type, {"1","4"}, {"3","1","g","h"})
// 3j - crossing+qqx
&& match_expectation(type, {"1","g"}, {"-1","1","g","1"})
&& match_expectation(type, {"1","g"}, {"-1","1","1","g"})
&& match_expectation(type, {"g","1"}, {"1","g","1","-1"})
&& match_expectation(type, {"g","1"}, {"g","1","1","-1"})
&& match_expectation(type, {"1","g"}, {"2","-2","g","1"})
&& match_expectation(type, {"1","g"}, {"2","-2","1","g"})
&& match_expectation(type, {"g","1"}, {"1","g","-2","2"})
&& match_expectation(type, {"g","1"}, {"g","1","-2","2"})
&& match_expectation(type, {"1","g"}, {"-1","1","h","g","1"})
&& match_expectation(type, {"1","g"}, {"-1","h","1","1","g"})
&& match_expectation(type, {"g","1"}, {"1","g","1","h","-1"})
&& match_expectation(type, {"g","1"}, {"h","g","1","1","-1"})
&& match_expectation(type, {"1","g"}, {"2","-2","1","g","h"})
&& match_expectation(type, {"g","1"}, {"g","h","1","-2","2"})
&& match_expectation(type, {"1","g"}, {"Wp","3","-4","g","1"})
&& match_expectation(type, {"3","g"}, {"-2","Wm","1","3","g"})
&& match_expectation(type, {"g","1"}, {"1","g","Wm","-3","4"})
&& match_expectation(type, {"g","-3"}, {"g","-3","-1","Wp","2"})
// 4j- gluon in qqx
&& match_expectation(type, {"g","1"}, {"1","g","-1","1"})
&& match_expectation(type, {"1","g"}, {"1","-1","g","1"})
&& match_expectation(type, {"g","1"}, {"1","g","Wm","-2","1"})
&& match_expectation(type, {"2","g"}, {"2","-2","g","Wp","1"})
&& match_expectation(type, {"g","g"}, {"Wp","3","g","-4","g"})
&& match_expectation(type, {"1","g"}, {"1","h","-1","g","1"})
// 6j - two qqx
&& match_expectation(type, {"g","g"}, {"1","-1","g","g","1","-1"})
&& match_expectation(type, {"g","g"}, {"1","-1","g","1","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","-1","g","1","-1"})
&& match_expectation(type, {"g","g"}, {"g","1","-1","1","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","1","-1","-1","g"})
&& match_expectation(type, {"g","g"}, {"h","1","-1","g","g","1","-1"})
&& match_expectation(type, {"g","g"}, {"1","Wp","-2","g","1","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","Wp","-1","g","1","-2"})
&& match_expectation(type, {"g","g"}, {"g","1","-1","Wm","2","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","2","-1","Wm","-1","g"})
// random stuff (can be non-physical)
&& match_expectation(type, {"g","g"}, {"1","-2","2","-1"}) // != 2 qqx
&& match_expectation(type, {"g","g"}, {"1","-2","2","g"}) // could be qqx
&& match_expectation(type, {"e+","e-"},{"1","-1"}) // bad initial state
&& match_expectation(type, {"1","e-"}, {"g","1","Wm"}) // bad initial state
&& match_expectation(type, {"h","g"}, {"g","g"}) // bad initial state
&& match_expectation(type, {"-1","g"}, {"-1","1","1"}) // bad qqx
&& match_expectation(type, {"-1","g"}, {"1","1","-1"}) // crossing in bad qqx
&& match_expectation(type, {"-1","g"}, {"-2","1","1","Wp"}) // bad qqx
&& match_expectation(type, {"1","2"}, {"1","-1","g","g","g","2"}) // bad qqx
&& match_expectation(type, {"1","2"}, {"1","-1","-2","g","g","2"}) // gluon in bad qqx
&& match_expectation(type, {"g","g"}, {"-1","2","g","g"}) // wrong back qqx
&& match_expectation(type, {"g","g"}, {"g","g","2","1"}) // wrong forward qqx
&& match_expectation(type, {"g","g"}, {"g","-2","1","g"}) // wrong central qqx
&& match_expectation(type, {"1","g"}, {"1","-2","g","g","Wp"}) // extra quark
&& match_expectation(type, {"g","1"}, {"g","g","-2","1","Wp"}) // extra quark
&& match_expectation(type, {"g","1"}, {"g","g","Wp","-2","1"}) // extra quark
&& match_expectation(type, {"g","1"}, {"g","-2","1","g","Wp"}) // extra quark
&& match_expectation(type, {"g","g"}, {"g","g","g","-2","1","-1","Wp"}) // extra quark
&& match_expectation(type, {"1","g"}, {"g","Wp","1","-2","g"}) // extra quark
&& match_expectation(type, {"g","g"}, {"1","-1","-2","g","g","g","Wp"}) // extra quark
;
}
// Two boson states, that are currently not implemented
bool check_bad_FS(){
auto type{ HEJ::event_type::bad_final_state};
return
match_expectation(type, {"g","g"}, {"g","h","h","g"})
&& match_expectation(type, {"g","g"}, {"h","g","h","g"})
&& match_expectation(type, {"g","-1"}, {"g","h","Wp","-2"})
&& match_expectation(type, {"-3","-1"},{"-4","g","Wp","Wp","-2"})
&& match_expectation(type, {"-4","-1"},{"-3","Wp","g","Wm","-2"})
&& match_expectation(type, {"-4","-1"},{"g","-3","Wp","Wm","-2"})
&& match_expectation(type, {"-4","-1"},{"-4","g","11","-11","-2"})
&& match_expectation(type, {"-4","-1"},{"-4","g","-13","g","-2"})
&& match_expectation(type, {"3","-2"}, {"Wp","3","Wm","g","g","g","-2"}, -13)
;
}
// not 2 jets
bool check_not_2_jets(){
auto type{ HEJ::event_type::no_2_jets};
return
match_expectation(type, {"g","g"}, {})
&& match_expectation(type, {"1","-1"}, {})
&& match_expectation(type, {"g","-1"}, {"-1"})
&& match_expectation(type, {"g","g"}, {"g"})
;
}
// not implemented processes
bool check_not_implemented(){
return check_fkl(false)
&& check_uno(false)
&& check_extremal_qqx(false)
&& check_central_qqx(false);
}
}
int main() {
// tests for "no false negatives"
// i.e. all HEJ-configurations get classified correctly
if(!check_fkl()) return EXIT_FAILURE;
if(!check_uno()) return EXIT_FAILURE;
if(!check_extremal_qqx()) return EXIT_FAILURE;
if(!check_central_qqx()) return EXIT_FAILURE;
// test for "no false positive"
// i.e. non-resummable gives non-resummable
if(!check_non_resummable()) return EXIT_FAILURE;
if(!check_bad_FS()) return EXIT_FAILURE;
if(!check_not_2_jets()) return EXIT_FAILURE;
if(!check_not_implemented()) return EXIT_FAILURE;
return EXIT_SUCCESS;
}
diff --git a/t/test_classify_ref.cc b/t/test_classify_ref.cc
index 6e08432..fd11a46 100644
--- a/t/test_classify_ref.cc
+++ b/t/test_classify_ref.cc
@@ -1,78 +1,77 @@
/**
* \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"
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};
}
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);
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){
- using HEJ::event_type::name;
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;
}