Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8310252
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
230 KB
Subscribers
None
View Options
diff --git a/src/Event.cc b/src/Event.cc
index e9ee7cf..0a9544b 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,889 +1,942 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Event.hh"
#include <algorithm>
#include <assert.h>
#include <numeric>
#include <unordered_set>
#include <utility>
#include "LHEF/LHEF.h"
#include "fastjet/JetDefinition.hh"
#include "HEJ/Constants.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/PDG_codes.hh"
namespace HEJ{
namespace {
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(std::vector<Particle> const & outgoing){
+ 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(auto const & out: outgoing){
+ 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; // pure jets
+ 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(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 qqx 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
){
// Is this unordered emisson?
if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){
if( is_valid_impact_factor(
incoming_id, (begin_out+1)->type, false, W_change )
){
// veto Higgs inside uno
assert((begin_out+1)<end_out);
if( !boson.empty() && boson.front().type == ParticleID::h
){
if( (backward && boson.front().rapidity() < (begin_out+1)->rapidity())
||(!backward && boson.front().rapidity() > (begin_out+1)->rapidity()))
return non_resummable;
}
begin_out+=2;
return not_tested | (backward?unob:unof);
}
}
// Is this QQbar?
else if( incoming_id==pid::gluon ){
if( is_valid_impact_factor(
begin_out->type, (begin_out+1)->type, true, W_change )
){
// veto Higgs inside qqx
assert((begin_out+1)<end_out);
if( !boson.empty() && boson.front().type == ParticleID::h
){
if( (backward && boson.front().rapidity() < (begin_out+1)->rapidity())
||(!backward && boson.front().rapidity() > (begin_out+1)->rapidity()))
return non_resummable;
}
begin_out+=2;
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,
OutIterator & qqx_pos
){
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
if( is_valid_impact_factor(
begin_out->type, (begin_out+1)->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() < (begin_out+1)->rapidity()
){
return non_resummable;
}
qqx_pos=begin_out;
begin_out+=2;
// 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;
}
bool invalid_jet(std::unordered_set<int> & other, int const idx){
if(idx<0) return true;
if(other.find(idx) != other.cend()) return true;
other.insert(idx);
return false;
}
bool jets_ok( size_t const final_type,
std::vector<int> const & jet_idx, size_t const qqx_pos
){
using namespace event_type;
std::unordered_set<int> other;
auto idx_begin{jet_idx.cbegin()};
auto idx_end{jet_idx.crbegin()};
// always seperate extremal jets
if(invalid_jet(other, *idx_begin)) return false;
if(invalid_jet(other, *idx_end)) return false;
// unob -> second parton in own jet
if( (final_type & (unob | qqxexb))
&& invalid_jet(other, *(idx_begin+1)) ) return false;
if( (final_type & (unof | qqxexf))
&& invalid_jet(other, *(idx_end+1)) ) return false;
assert( !(final_type & qqxmid) || jet_idx.size()>qqx_pos+1 );
if( (final_type & qqxmid)
&& ( invalid_jet(other, *(idx_begin+qqx_pos))
|| invalid_jet(other, *(idx_begin+qqx_pos+1)) ) ) return false;
return true;
}
/**
* \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.outgoing()))
+ if(! final_state_ok(ev))
return bad_final_state;
// initialise variables
auto const & in = ev.incoming();
auto const & out = filter_partons(ev.outgoing());
assert(std::distance(begin(in), end(in)) == 2);
assert(out.size() >= 2);
assert(std::distance(begin(out), end(out)) >= 2);
assert(std::is_sorted(begin(out), end(out), 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() && abs(boson.front().type) == ParticleID::Wp ){
if(boson.front().type>0) ++remaining_Wp;
else ++remaining_Wm;
}
int W_change = 0;
// range for current checks
auto begin_out{out.cbegin()};
auto end_out{out.crbegin()};
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
auto qqx_pos{out.cend()};
final_type &= possible_central(
begin_out, end_out.base(), W_change, boson, qqx_pos );
if( final_type == non_resummable )
return non_resummable;
assert( !(final_type&qqxmid) || qqx_pos != out.cend() );
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 jet configurations
if(!jets_ok( final_type,
ev.particle_jet_indices( ev.jets() ),
std::distance( out.cbegin(), qqx_pos) ))
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, int 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.weight()
};
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();
Event ev{ std::move(incoming), std::move(outgoing), std::move(decays),
std::move(parameters),
jet_def, min_jet_pt
};
assert(std::is_sorted(begin(ev.outgoing_), end(ev.outgoing_),
rapidity_less{}));
ev.type_ = classify(ev);
return ev;
}
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_));
}
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;
}
}
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);
for(auto const & part: outgoing()){
// reasonable colour
if(!correct_colour(part))
return false;
if(!is_parton(part)) // skip colour neutral particles
continue;
// if possible connect to line
if( line_colour.first == part.colour->second )
line_colour.first = part.colour->first;
else if( line_colour.second == part.colour->first )
line_colour.second = part.colour->second;
else
return false;
// 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 {
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;
}
}
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);
for(auto & part: outgoing_){
assert(colour>0 || anti_colour>0);
if(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){
part.colour = std::make_pair(colour+2,colour);
colour+=2;
} else {
part.colour = std::make_pair(anti_colour, anti_colour+2);
anti_colour+=2;
}
} else if(colour > 0){
// on q line => connect to available colour
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
part.colour = std::make_pair(anti_colour, anti_colour+2);
anti_colour+=2;
}
} else if(is_quark(part)) {
// quark
assert(anti_colour>0);
if(colour>0){
// on g line => connect and remove anti-colour
part.colour = std::make_pair(anti_colour, 0);
anti_colour+=2;
anti_colour*=-1;
} else {
// on qx line => new colour
colour*=-1;
part.colour = std::make_pair(colour, 0);
}
} else if(is_antiquark(part)) {
// anti-quark
assert(colour>0);
if(anti_colour>0){
// on g line => connect and remove colour
part.colour = std::make_pair(0, colour);
colour+=2;
colour*=-1;
} else {
// on q line => new anti-colour
anti_colour*=-1;
part.colour = std::make_pair(0, anti_colour);
}
} else { // not a parton
assert(!is_parton(part));
part.colour = {};
}
}
// Connect last
connect_incoming(incoming_[1], anti_colour, colour);
assert(is_leading_colour());
return true;
} // generate_colours
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())
);
};
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 << "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{
assert(is_AWZH_boson(out));
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/src/MatrixElement.cc b/src/MatrixElement.cc
index 1cb9f03..143d1d4 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1523 +1,1542 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <unordered_map>
#include <utility>
#include "CLHEP/Vector/LorentzVector.h"
#include "HEJ/Constants.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Hjets.hh"
#include "HEJ/jets.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/utility.hh"
namespace HEJ{
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const {
const double lambda = param_.regulator_lambda;
const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
// use alpha_s(sqrt(q_j*lambda)), evolved to mur
return (
1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(Event const & event) const {
return tree(event)*virtual_corrections(event);
}
Weights MatrixElement::tree(Event const & event) const {
return tree_param(event)*tree_kin(event);
}
Weights MatrixElement::tree_param(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
Weights MatrixElement::virtual_corrections(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = virtual_corrections(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = virtual_corrections(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
double mur,
Particle const & WBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
#ifndef NDEBUG
bool wc = true;
#endif
bool wqq = false;
// With extremal qqx or unordered gluon outside the extremal
// partons then it is not part of the FKL ladder and does not
// contribute to the virtual corrections. W emitted from the
// most backward leg must be taken into account in t-channel
if (event.type() == event_type::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else if (event.type() == event_type::qqxexb) {
q -= partons[1].p;
++first_idx;
if (abs(partons[0].type) != abs(partons[1].type)){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else {
if(event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
if (in[0].type != partons[0].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(abs(backquark->type) != abs((backquark+1)->type)) {
wqq=true;
#ifndef NDEBUG
wc=false;
#endif
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -=partons[j+1].p;
} // End Loop one
if (last_idx != first_idx_qqx) q -= partons[last_idx+1].p;
if (wqq) q -= WBoson.p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -= partons[j+1].p;
}
#ifndef NDEBUG
if (wc) q -= WBoson.p;
assert(
nearby(q, -1*pb, norm)
|| is_AWZH_boson(partons.back().type)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
#endif
return exp(exponent);
}
double MatrixElement::virtual_corrections(
Event const & event,
double mur
) const{
auto const & in = event.incoming();
auto const & out = event.outgoing();
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
const auto AWZH_boson = std::find_if(
begin(out), end(out),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson != end(out) && abs(AWZH_boson->type) == pid::Wp){
return virtual_corrections_W(event, mur, *AWZH_boson);
}
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
size_t first_idx = 0;
size_t last_idx = out.size() - 1;
// if there is a Higgs boson, extremal qqx or unordered gluon
// outside the extremal partons then it is not part of the FKL
// ladder and does not contribute to the virtual corrections
if((out.front().type == pid::Higgs)
|| event.type() == event_type::unob
|| event.type() == event_type::qqxexb){
q -= out[1].p;
++first_idx;
}
if((out.back().type == pid::Higgs)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(out) + 1, end(out) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
);
if(backquark == end(out) || (backquark+1)->type==pid::gluon) return 0;
last_idx = std::distance(begin(out), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
if (last_idx != first_idx_qqx) q -= out[last_idx+1].p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| out.back().type == pid::Higgs
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2
){
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
return Cls-4./(kperp*kperp);
}
//! Lipatov vertex
double C2Lipatov( // B
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim,
CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom,
CLHEP::HepLorentzVector const & pop
){
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
double temp=Cls-4./(kperp*kperp);
return temp;
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pg Unordered gluon momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_uno_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
assert(aptype!=21); // aptype cannot be gluon
if (bptype==21) {
if (aptype > 0)
return ME_unob_qg(pg,p1,pa,pn,pb);
else
return ME_unob_qbarg(pg,p1,pa,pn,pb);
}
else if (bptype<0) { // ----- || -----
if (aptype > 0)
return ME_unob_qQbar(pg,p1,pa,pn,pb);
else
return ME_unob_qbarQbar(pg,p1,pa,pn,pb);
}
else { //bptype == quark
if (aptype > 0)
return ME_unob_qQ(pg,p1,pa,pn,pb);
else
return ME_unob_qbarQ(pg,p1,pa,pn,pb);
}
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
if (aptype==21&&bptype==21) {
return ME_gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_qg(pn,pb,p1,pa);
else
return ME_qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return ME_qg(p1,pa,pn,pb);
else
return ME_qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_qQ(pn,pb,p1,pa);
else
return ME_qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return ME_qQbar(p1,pa,pn,pb);
else
return ME_qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_W_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc, ParticleProperties const & Wprop
){
// We know it cannot be gg incoming.
assert(!(aptype==21 && bptype==21));
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_W_qg(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qbarg(pn,plbar,pl,pb,p1,pa,Wprop);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return ME_W_qg(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qbarg(p1,plbar,pl,pa,pn,pb,Wprop);
}
else { // they are both quark
if (wc==true){ // emission off b, (first argument pbout)
if (bptype>0) {
if (aptype>0)
return ME_W_qQ(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
else {
if (aptype>0)
return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
}
else{ // emission off a, (first argument paout)
if (aptype > 0) {
if (bptype > 0)
return ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
else { // a is anti-quark
if (bptype > 0)
return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_W_uno_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc, ParticleProperties const & Wprop
){
// we know they are not both gluons
if (bptype == 21 && aptype != 21) { // b gluon => W emission off a
if (aptype > 0)
return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else { // they are both quark
if (wc) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else{
if (aptype>0)
return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else //qqbar
return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else //qbarqbar
return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Matrix element squared for backward qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxb Tree-Level Current-Current Scattering
*
* @note calculate forwards qqx contribution by reversing argument ordering.
*/
double ME_W_qqx_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const swap_q_qx, bool const wc,
ParticleProperties const & Wprop
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
// const bool swap_q_qx= pqbar.rapidity() > pq.rapidity();
const double CFbackward = K_g( (swap_q_qx)?pq:pqbar ,pa)/HEJ::C_F;
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swap_q_qx)
return ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb, Wprop)*CFbackward;
else
return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)*CFbackward;
}
else if (aptype==21&&bptype!=21 ) {//a gluon => W emission off b leg or qqx
if (!wc){ // W Emitted from backwards qqx
if (swap_q_qx)
return ME_WExqqx_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop)*CFbackward;
else
return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)*CFbackward;
}
else { // W Must be emitted from forwards leg.
if (swap_q_qx)
return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, bptype<0, Wprop)*CFbackward;
else
return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, bptype<0, Wprop)*CFbackward;
}
}
throw std::logic_error("Incompatible incoming particle types with qqxb");
}
/* \brief Matrix element squared for central qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqxpair
* @param nbelow Number of gluons emitted after central qqxpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum\
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wqq Boolean. True siginfies W boson is emitted from Central qqx
* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
*/
double ME_W_qqxmid_current(
int aptype, int bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<HLV> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc,
ParticleProperties const & Wprop
){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
const bool swap_q_qx=pqbar.rapidity() < pq.rapidity();
double wt=1.;
if (aptype==21) wt*=K_g(partons.front(),pa)/HEJ::C_F;
if (bptype==21) wt*=K_g(partons.back(),pb)/HEJ::C_F;
if(wqq)
return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons,(bptype<0),(aptype<0),
swap_q_qx, nabove, Wprop);
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, (bptype<0), (aptype<0),
swap_q_qx, nabove, nbelow, wc, Wprop);
}
/** \brief Matrix element squared for tree-level current-current scattering with Higgs
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
*/
double ME_Higgs_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb, double vev
){
if (aptype==21&&bptype==21) // gg initial state
return ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_H_qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
else
return ME_H_qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return ME_H_qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
else
return ME_H_qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
else
return ME_H_qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return ME_H_qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
else
return ME_H_qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param pg Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*
* @note This function assumes unordered gluon backwards from pa-p1 current.
* For unof, reverse call order
*/
double ME_Higgs_current_uno(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb, double vev
){
if (bptype==21&&aptype!=21) {
if (aptype > 0)
return ME_H_unob_gQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_gQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return ME_H_unob_qQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_qbarQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
else {
if (bptype>0)
return ME_H_unob_qQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_qbarQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
}
throw std::logic_error("unknown particle types");
}
CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
void validate(HEJ::MatrixElementConfig const & config) {
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(!config.Higgs_coupling.use_impact_factors) {
throw std::invalid_argument{
"Invalid Higgs coupling settings.\n"
"HEJ without QCDloop support can only use impact factors.\n"
"Set use_impact_factors to true or recompile HEJ.\n"
};
}
#endif
if(config.Higgs_coupling.use_impact_factors
&& config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) {
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
} // namespace anonymous
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
double MatrixElement::tree_kin(
Event const & ev
) const {
if(! is_resummable(ev.type())) return 0.;
auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing()))
return tree_kin_jets(ev);
switch(AWZH_boson->type){
case pid::Higgs:
return tree_kin_Higgs(ev);
case pid::Wp:
case pid::Wm:
return tree_kin_W(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
namespace{
constexpr int extremal_jet_idx = 1;
constexpr int no_extremal_jet_idx = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == extremal_jet_idx;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
double lambda
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
Event const & ev
) const{
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
const auto & jets = ev.jets();
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission or qqx
if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){
assert(jets.size() >= 3);
++most_backward;
}
else if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
assert(jets.size() >= 3);
--most_forward;
}
const auto extremal_jet_indices = ev.particle_jet_indices(
{*most_backward, *most_forward}
);
assert(extremal_jet_indices.size() == out_partons.size());
for(size_t i = 0; i < out_partons.size(); ++i){
assert(HEJ::is_parton(out_partons[i]));
const int idx = (extremal_jet_indices[i]>=0)?
extremal_jet_idx:
no_extremal_jet_idx;
out_partons[i].p.set_user_index(idx);
}
return out_partons;
}
namespace {
template<class InIter, class partIter>
double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const double current_factor = ME_uno_current(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pa-p1-pg, pa, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
}
double MatrixElement::tree_kin_jets(Event const & ev) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
if (ev.type()==HEJ::event_type::FKL){
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
pa - p1, pa, pb, p1, pn,
param_.regulator_lambda
);
}
else if (ev.type()==HEJ::event_type::unordered_backward){
return tree_kin_jets_uno(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
else if (ev.type()==HEJ::event_type::unordered_forward){
return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
else {
throw std::logic_error("Can only reweight FKL or uno processes in Pure Jets");
}
}
namespace{
double tree_kin_W_FKL(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda, ParticleProperties const & Wprop
){
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_uno(InIter BeginIn, partIter BeginPart,
partIter EndPart, const HLV & plbar, const HLV & pl,
double lambda, ParticleProperties const & Wprop){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
bool wc = (BeginIn)->type==(BeginPart+1)->type; //leg b emits w
auto q0 = pa - p1 - pg;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_uno_current(
(BeginIn)->type, (BeginIn+1)->type, pn, pb,
p1, pa, pg, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_qqx(InIter BeginIn, partIter BeginPart,
partIter EndPart, const HLV & plbar, const HLV & pl,
double lambda, ParticleProperties const & Wprop){
const bool swap_q_qx=is_quark(*BeginPart);
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pq = to_HepLorentzVector(*(BeginPart+(swap_q_qx?0:1)));
const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_q_qx?1:0)));
const auto p1 = to_HepLorentzVector(*(BeginPart));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const bool wc = (BeginIn+1)->type!=(EndPart-1)->type; //leg b emits w
auto q0 = pa - pq - pqbar;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_qqx_current(
(BeginIn)->type, (BeginIn+1)->type, pa, pb,
pq, pqbar, pn, plbar, pl, swap_q_qx, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxmid(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda, ParticleProperties const & Wprop
){
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
bool wc, wqq;
if (backmidquark->type == -(backmidquark+1)->type){ // Central qqx does not emit
wqq=false;
if (aptype==partons[0].type) {
wc = true;
}
else{
wc = false;
q0-=pl+plbar;
}
}
else{
wqq = true;
wc = false;
qqxt-=pl+plbar;
}
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
int nabove = std::distance(begin_ladder, backmidquark);
int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqxt, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
} // namespace anonymous
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
- auto const & decays(ev.decays());
+
+ #ifndef NDEBUG
+ // assert that there is exactly one decay corresponding to the W
+ assert(ev.decays().size() == 1);
+ auto const & w_boson{
+ std::find_if(ev.outgoing().cbegin(), ev.outgoing().cend(),
+ [] (Particle const & p) -> bool {
+ return std::abs(p.type) == ParticleID::Wp;
+ }) };
+ assert(w_boson != ev.outgoing().cend());
+ assert( (long int) ev.decays().cbegin()->first
+ == std::distance(ev.outgoing().cbegin(), w_boson) );
+ #endif
+
+ // find decay products of W
+ auto const & decay{ ev.decays().cbegin()->second };
+ assert(decay.size() == 2);
+ assert( ( is_anylepton(decay.at(0)) && is_anyneutrino(decay.at(1)) )
+ || ( is_anylepton(decay.at(1)) && is_anyneutrino(decay.at(0)) ) );
+
+ // get lepton & neutrino
HLV plbar, pl;
- for (auto& x: decays) {
- if (x.second.at(0).type < 0){
- plbar = to_HepLorentzVector(x.second.at(0));
- pl = to_HepLorentzVector(x.second.at(1));
- }
- else{
- pl = to_HepLorentzVector(x.second.at(0));
- plbar = to_HepLorentzVector(x.second.at(1));
- }
+ if (decay.at(0).type < 0){
+ plbar = to_HepLorentzVector(decay.at(0));
+ pl = to_HepLorentzVector(decay.at(1));
}
+ else{
+ pl = to_HepLorentzVector(decay.at(0));
+ plbar = to_HepLorentzVector(decay.at(1));
+ }
+
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
if(ev.type() == FKL){
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_backward){
return tree_kin_W_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_forward){
return tree_kin_W_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqxb){
return tree_kin_W_qqx(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqxf){
return tree_kin_W_qqx(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
assert(ev.type() == central_qqx);
return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
double MatrixElement::tree_kin_Higgs(Event const & ev) const {
if(is_uno(ev.type())){
return tree_kin_Higgs_between(ev);
}
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev);
}
return tree_kin_Higgs_between(ev);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
#ifdef HEJ_BUILD_WITH_QCDLOOP
// TODO: code duplication with jets.cc
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(type == ParticleID::gluon) return K_g(pout, pin);
return C_F;
}
#endif
// Colour factor in strict MRK limit
double K_MRK(ParticleID type) {
return (type == ParticleID::gluon)?C_A:C_F;
}
}
double MatrixElement::MH2_forwardH(
CLHEP::HepLorentzVector const & p1out,
CLHEP::HepLorentzVector const & p1in,
ParticleID type2,
CLHEP::HepLorentzVector const & p2out,
CLHEP::HepLorentzVector const & p2in,
CLHEP::HepLorentzVector const & pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
const double vev = param_.ew_parameters.vev();
// gluon case
#ifdef HEJ_BUILD_WITH_QCDLOOP
if(!param_.Higgs_coupling.use_impact_factors){
return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*ME_Houtside_gq(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb, vev
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH,vev) + C2gHgm(p1in,p1out,pH,vev)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(outgoing[1].type != pid::gluon) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
const auto q0 = pa - p1 - pH;
const double t1 = q0.m2();
const double t2 = (pn - pb).m2();
return MH2_forwardH(
p1, pa, incoming[1].type, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_last(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.back().type == pid::Higgs);
if(outgoing[outgoing.size()-2].type != pid::gluon) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
const double t1 = q0.m2();
const double t2 = (pn + pH - pb).m2();
return MH2_forwardH(
pn, pb, incoming[0].type, p1, pa, pH,
t2, t1
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
namespace {
template<class InIter, class partIter>
double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, const HLV & qH, const HLV & qHp1,
double mt, bool inc_bot, double mb, double vev){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
return ME_Higgs_current_uno(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa,
qH, qHp1, mt, inc_bot, mb, vev
);
}
}
double MatrixElement::tree_kin_Higgs_between(Event const & ev) const {
using namespace event_type;
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[(ev.type() == unob)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - ((ev.type() == unof)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
(ev.type() == unob)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
(ev.type() == unof)
|| partons.front().type != pid::gluon
))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor;
if(ev.type() == FKL){
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
}
else if(ev.type() == unob){
current_factor = HEJ::C_A*HEJ::C_A/2*tree_kin_Higgs_uno(
begin(incoming), end(incoming), begin(partons),
end(partons), qH, qH-pH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
current_factor = HEJ::C_A*HEJ::C_A/2*tree_kin_Higgs_uno(
rbegin(incoming), rend(incoming), rbegin(partons),
rend(partons), qH-pH, qH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
throw std::logic_error("Can only reweight FKL or uno processes in H+Jets");
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn,
param_.regulator_lambda
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn,
param_.regulator_lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
namespace {
double get_AWZH_coupling(Event const & ev, double alpha_s, double alpha_w) {
const auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](auto const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing())) return 1.;
switch(AWZH_boson->type){
case pid::Higgs:
return alpha_s*alpha_s;
case pid::Wp:
case pid::Wm:
return alpha_w*alpha_w;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
}
double MatrixElement::tree_param(Event const & ev, double mur) const {
assert(is_resummable(ev.type()));
const auto begin_partons = ev.begin_partons();
const auto end_partons = ev.end_partons();
const auto num_partons = std::distance(begin_partons, end_partons);
const double alpha_s = alpha_s_(mur);
const double gs2 = 4.*M_PI*alpha_s;
double res = std::pow(gs2, num_partons);
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(num_partons >= 2);
const auto first_emission = std::next(begin_partons);
const auto last_emission = std::prev(end_partons);
for(auto parton = first_emission; parton != last_emission; ++parton){
res *= 1. + alpha_s/(2.*M_PI)*beta0*log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res;
}
} // namespace HEJ
diff --git a/t/CMakeLists.txt b/t/CMakeLists.txt
index df4ce46..e704bef 100644
--- a/t/CMakeLists.txt
+++ b/t/CMakeLists.txt
@@ -1,277 +1,289 @@
set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}")
set(tst_ME_data_dir "${tst_dir}/ME_data")
# test event classification
add_executable(test_classify ${tst_dir}/test_classify.cc)
target_link_libraries(test_classify HEJ)
add_test(
NAME t_classify
COMMAND test_classify
)
add_executable(test_classify_ref ${tst_dir}/test_classify_ref.cc)
target_link_libraries(test_classify_ref HEJ)
add_test(
NAME t_classify_ref
COMMAND test_classify_ref ${tst_dir}/classify_ref ${tst_dir}/classify.lhe.gz
)
add_test(
NAME t_classify_ref_4j
COMMAND test_classify_ref ${tst_dir}/classify_ref_4j ${tst_dir}/4j.lhe.gz
)
+add_test(
+ NAME t_classify_ref_W4j
+ COMMAND test_classify_ref ${tst_dir}/classify_ref_W4j ${tst_dir}/MGemnubar4j_unweighted.lhe.gz
+ )
+
+# test for valid W decays
+add_executable(test_decay ${tst_dir}/test_decay.cc)
+target_link_libraries(test_decay HEJ)
+add_test(
+ NAME t_valid_decay
+ COMMAND test_decay
+ )
# test phase space point
add_executable(test_psp ${tst_dir}/test_psp.cc)
target_link_libraries(test_psp HEJ)
add_test(
NAME t_psp
COMMAND test_psp ${tst_dir}/psp_gen.lhe.gz
)
# test importing scales
add_library(scales SHARED ${tst_dir}/scales.cc)
target_link_libraries(scales HEJ)
add_executable(test_scale_import ${tst_dir}/test_scale_import)
target_link_libraries(test_scale_import HEJ)
add_test(
NAME t_scale_import
COMMAND test_scale_import ${tst_dir}/jet_config_with_import.yml
)
# test scale arithmetic (e.g. 2*H_T/4)
add_executable(test_scale_arithmetics ${tst_dir}/test_scale_arithmetics)
target_link_libraries(test_scale_arithmetics HEJ)
add_test(
NAME t_scale_arithmetics
COMMAND test_scale_arithmetics ${tst_dir}/jet_config.yml ${tst_dir}/2j.lhe.gz
)
# test "ParameterDescription"
add_executable(test_descriptions ${tst_dir}/test_descriptions)
target_link_libraries(test_descriptions HEJ)
add_test(
NAME t_descriptions
COMMAND test_descriptions
)
# test "EventParameters*Weight"
add_executable(test_parameters ${tst_dir}/test_parameters)
target_link_libraries(test_parameters HEJ)
add_test(
NAME test_parameters
COMMAND test_parameters
)
# test unweighting
add_executable(test_unweighter ${tst_dir}/test_unweighter)
target_link_libraries(test_unweighter HEJ)
add_test(
NAME test_unweighter
COMMAND test_unweighter ${tst_dir}/4j.lhe.gz
)
# test colour generation
add_executable(test_colours ${tst_dir}/test_colours)
target_link_libraries(test_colours HEJ)
add_test(
NAME t_colour_flow
COMMAND test_colours
)
# test matrix elements
add_executable(test_ME_generic ${tst_dir}/test_ME_generic.cc)
target_link_libraries(test_ME_generic HEJ)
add_test(
NAME t_ME_j
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_jets_tree.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz
)
add_test(
NAME t_ME_j_virt
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_jets_tree_virt.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz
)
add_test(
NAME t_ME_h
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_mtinf_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
)
add_test(
NAME t_ME_h_virt
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_mtinf_virt.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
)
if(${QCDloop_FOUND})
add_test(
NAME t_ME_h_mt
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mt.yml ${tst_ME_data_dir}/ME_h_mt_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
)
add_test(
NAME t_ME_h_mtmb
COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtmb.yml ${tst_ME_data_dir}/ME_h_mtmb_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
)
endif()
add_test(
NAME t_ME_w_FKL
COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_w_FKL_tree.dat ${tst_ME_data_dir}/PSP_w_FKL.lhe.gz
)
add_test(
NAME t_ME_w_FKL_virt
COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_w_FKL_virt.dat ${tst_ME_data_dir}/PSP_w_FKL.lhe.gz
)
add_test(
NAME t_ME_Wp
COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wp.dat ${tst_ME_data_dir}/PSP_Wp.lhe.gz
)
add_test(
NAME t_ME_Wp_virt
COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wp_virt.dat ${tst_ME_data_dir}/PSP_Wp.lhe.gz
)
add_test(
NAME t_ME_Wm
COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wm.dat ${tst_ME_data_dir}/PSP_Wm.lhe.gz
)
add_test(
NAME t_ME_Wm_virt
COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wm_virt.dat ${tst_ME_data_dir}/PSP_Wm.lhe.gz
)
# test main executable
file(COPY "${tst_dir}/jet_config.yml" DESTINATION "${CMAKE_BINARY_DIR}")
set(test_config "${CMAKE_BINARY_DIR}/jet_config.yml")
if(${HepMC3_FOUND})
file(READ ${test_config} config)
file(WRITE ${test_config} "${config} - tst.hepmc\n")
endif()
if(${HepMC_FOUND})
file(READ ${test_config} config)
file(WRITE ${test_config} "${config} - tst.hepmc2\n")
if(${rivet_FOUND})
file(READ ${test_config} config)
file(WRITE ${test_config} "${config}\nanalysis:\n rivet: MC_XS\n output: tst\n")
endif()
endif()
set(test_cmd_main "$<TARGET_FILE:HEJ_main>\\\;${test_config}\\\;${tst_dir}/2j.lhe.gz")
# check that HepMC3 output is correct
if(${HepMC3_FOUND})
add_executable(check_hepmc ${tst_dir}/check_hepmc.cc)
target_link_libraries(check_hepmc HEJ ${HEPMC3_LIBRARIES})
target_include_directories(check_hepmc PRIVATE ${HEPMC3_INCLUDE_DIR})
set(test_cmd_hepmc "$<TARGET_FILE:check_hepmc>\\\;tst.hepmc")
else()
set(test_cmd_hepmc "")
endif()
# check that LHEF output is correct
add_executable(check_lhe ${tst_dir}/check_lhe.cc)
target_link_libraries(check_lhe HEJ)
set(test_cmd_lhe "$<TARGET_FILE:check_lhe>\\\;tst.lhe")
# Run dependent tests in one command to ensure correct execution order
# Note: The commands are concatenated with "\;" to escape CMake lists.
# Thus arguments have to be escaped twice "\\\;".
# e.g. "cmd1\;cmd2\\\;arg1\\\;arg2" is executed like "cmd1 && cmd2 arg1 arg2"
add_test(
NAME t_main
COMMAND ${CMAKE_COMMAND}
-DCMDS=${test_cmd_main}\;${test_cmd_hepmc}\;${test_cmd_lhe}
-P ${PROJECT_SOURCE_DIR}/cmake/run_multiple_tests.cmake
)
# check HDF5 reader
if(${HighFive_FOUND})
add_executable(test_hdf5 ${tst_dir}/test_hdf5.cc)
target_link_libraries(test_hdf5 HEJ)
add_test(
NAME t_hdf5
COMMAND test_hdf5 ${tst_dir}/Wm9-g4-repack.hdf5
)
endif()
# check rivet interface
if(${RIVET_FOUND})
add_executable(check_rivet ${tst_dir}/check_rivet.cc)
target_link_libraries(check_rivet HEJ rivet::rivet HepMC::HepMC)
add_test(
NAME t_rivet
COMMAND check_rivet
)
endif()
# test boson reconstruction
add_executable(cmp_events ${tst_dir}/cmp_events.cc)
target_link_libraries(cmp_events HEJ)
add_test(
NAME t_epnu_2j_noW
COMMAND cmp_events ${tst_dir}/epnu2jLOFKL_unweight.lhe.tar.gz ${tst_dir}/epnu2jLOFKL_unweight_noW.lhe.gz
)
# test resummed result
add_executable(check_res ${tst_dir}/check_res.cc)
target_link_libraries(check_res HEJ)
if(${TEST_ALL}) # deactivate long tests by default
add_test(
NAME t_2j
COMMAND check_res ${tst_dir}/2j.lhe.gz 3.49391e+07 419684
)
add_test(
NAME t_3j
COMMAND check_res ${tst_dir}/3j.lhe.gz 2.37902e+06 25746.6
)
add_test(
NAME t_3j_unof
COMMAND check_res ${tst_dir}/3j.lhe.gz 133399 4688.83 unof
)
add_test(
NAME t_3j_unob
COMMAND check_res ${tst_dir}/3j.lhe.gz 105247 3449.45 unob
)
add_test(
NAME t_4j
COMMAND check_res ${tst_dir}/4j.lhe.gz 603713 72822.6
)
add_test(
NAME t_h_3j
COMMAND check_res ${tst_dir}/h_3j.lhe.gz 0.821622 0.0220334
)
add_test(
NAME t_h_3j_unof
COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0127362 0.000271555 unof
)
add_test(
NAME t_h_3j_unob
COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0130615 0.000224793 unob
)
add_test(
NAME t_epnu_2j
COMMAND check_res ${tst_dir}/epnu2jLOFKL_unweight.lhe.tar.gz 262.7 3
)
add_test(
NAME t_MGepnu_3j
COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 38.9512 1
)
add_test(
NAME t_MGemnubar_3j
COMMAND check_res ${tst_dir}/MGemnubar3j_unweighted.lhe.gz 24.1575 1
)
add_test(
NAME t_MGepnu_3j_unof
COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 9.63702 0.128355 unof
)
add_test(
NAME t_MGepnu_3j_unob
COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 9.70119 0.108436 unob
)
add_test(
NAME t_MGepnu_3j_splitf
COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 2.91995 0.0463182 splitf
)
add_test(
NAME t_MGepnu_3j_splitb
COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 3.40708 0.0550975 splitb
)
add_test(
NAME t_MGepnu_4j
COMMAND check_res ${tst_dir}/MGepnu4j_unweighted.lhe.gz 10.2542 0.135106
)
add_test(
NAME t_MGemnubar_4j
COMMAND check_res ${tst_dir}/MGemnubar4j_unweighted.lhe.gz 5.57909 0.0300496
)
add_test(
NAME t_MGepnu_4j_qqxmid
COMMAND check_res ${tst_dir}/MGepnu4j_unweighted.lhe.gz 0.732084 0.005 qqxmid
)
endif()
diff --git a/t/classify_ref_4j b/t/classify_ref_4j
index 701338a..056dd0f 100644
--- a/t/classify_ref_4j
+++ b/t/classify_ref_4j
@@ -1,3255 +1,3255 @@
4
4
0
0
4
0
4
4
4
16
4
0
0
0
4
0
0
4
0
4
0
0
4
0
0
4
0
4
0
4
0
0
0
-64
+0
4
4
2
0
0
2
4
0
0
4
4
4
4
0
4
0
0
4
0
2
4
4
0
0
0
4
0
4
0
4
4
-64
+0
4
4
0
4
4
4
0
4
4
0
0
-32
+0
0
4
4
-64
+0
2
0
0
4
0
0
0
0
4
0
0
4
4
0
0
0
4
0
0
0
0
0
0
0
4
2
4
0
0
4
0
0
4
0
4
4
0
0
0
4
0
4
4
4
4
-64
+0
0
4
4
0
4
4
4
2
0
0
4
4
0
0
0
4
0
4
4
0
0
4
4
2
4
2
4
4
0
0
0
4
4
0
4
0
8
4
-128
+0
4
0
4
4
2
0
-64
+0
2
2
8
4
4
-128
+0
4
0
8
4
4
16
8
0
4
0
0
0
0
4
0
0
4
4
0
0
4
4
0
2
0
0
4
4
0
4
4
0
0
4
4
4
0
2
0
4
4
4
4
0
4
0
0
0
0
4
0
4
4
-64
+0
0
0
4
0
0
4
0
4
0
0
0
4
4
0
4
4
4
0
4
0
0
4
0
0
-64
+0
0
8
4
4
4
4
0
4
8
4
8
4
4
4
4
0
4
4
4
0
0
4
0
0
4
0
4
0
4
16
0
0
0
2
0
0
4
0
0
4
0
-64
+0
8
4
0
0
0
0
-32
+0
4
4
4
16
0
0
4
16
4
4
0
0
4
0
0
2
0
0
-64
+0
0
0
0
0
0
4
0
4
0
0
4
2
4
4
0
0
4
4
4
0
0
4
0
0
4
0
0
4
0
4
0
0
4
4
4
4
0
0
-64
+0
4
4
4
2
0
2
0
0
0
4
0
4
4
0
4
0
4
4
0
4
0
2
0
0
4
0
0
0
-64
+0
4
0
4
0
0
0
4
0
4
4
0
4
0
4
4
16
0
4
4
0
0
4
0
2
4
0
4
4
-32
+0
4
0
4
4
2
0
4
0
8
8
4
0
0
0
4
2
0
4
0
4
0
4
4
4
4
0
4
0
0
2
0
4
0
4
0
4
4
0
0
4
4
4
0
4
4
0
0
4
4
4
4
0
2
0
0
0
4
2
4
4
4
4
4
2
0
4
4
0
4
0
0
4
2
0
4
0
4
4
4
4
0
4
0
0
4
0
0
0
0
0
-64
+0
0
4
0
0
0
4
4
4
4
0
0
4
4
4
4
4
4
16
4
4
4
0
4
4
4
4
0
2
4
4
0
16
0
0
-128
+0
0
2
4
4
4
-64
+0
4
0
0
4
0
0
4
0
0
4
4
0
4
4
4
0
4
4
4
4
0
4
0
-64
+0
4
0
0
0
0
4
4
4
2
4
0
0
4
0
0
0
4
0
0
0
0
4
-32
-128
+0
+0
4
2
4
4
0
4
0
0
4
0
16
0
4
0
4
-128
+0
2
8
4
0
0
0
0
0
4
0
4
4
4
0
0
0
0
4
4
0
0
0
4
0
0
0
2
4
4
0
0
4
2
0
0
4
4
4
-64
+0
4
4
0
0
2
4
4
0
0
0
0
4
4
4
2
0
-32
+0
0
0
4
4
4
4
2
4
4
0
0
0
-64
+0
0
0
0
0
4
0
2
2
4
4
4
4
-32
+0
0
0
4
0
0
0
4
4
0
0
8
4
2
4
0
4
0
4
2
0
4
0
0
0
4
4
0
4
0
4
0
4
0
4
4
-64
+0
4
4
4
4
4
4
2
2
2
0
0
0
4
0
4
0
4
0
0
4
4
4
4
8
0
8
0
0
0
0
4
4
0
2
2
0
2
4
-32
+0
0
4
4
4
0
4
2
0
4
0
4
0
0
4
0
4
0
4
0
4
0
4
4
4
0
0
0
4
4
4
4
0
0
0
0
4
0
4
0
0
2
4
4
0
0
2
4
0
2
8
4
4
0
0
4
0
0
4
0
4
4
4
4
-32
+0
4
4
4
4
4
0
0
4
0
0
4
2
4
0
4
0
2
0
0
0
4
4
4
4
0
4
0
0
4
4
0
4
4
8
0
4
0
4
0
4
4
-128
+0
4
0
0
0
2
0
16
0
0
0
4
4
0
0
4
0
0
4
0
2
0
0
4
0
4
0
4
0
4
-64
+0
4
16
4
4
4
0
2
2
0
0
8
4
0
4
0
0
4
4
4
0
-64
+0
4
-64
-32
+0
+0
0
0
4
-32
+0
0
0
4
0
8
0
0
0
4
4
2
4
0
2
4
0
2
0
0
0
4
0
4
4
4
4
0
0
-64
+0
4
0
4
0
4
4
2
4
4
4
4
0
0
4
4
-128
+0
0
2
0
0
0
4
0
0
0
0
4
4
0
4
4
0
0
4
4
0
0
0
0
0
0
0
4
0
0
-32
+0
0
4
0
0
0
4
0
0
4
-128
+0
0
0
4
0
0
4
0
0
0
0
-64
+0
4
0
4
0
4
2
0
4
4
4
0
0
0
4
4
-32
+0
0
2
0
4
4
16
4
4
4
0
0
0
4
0
0
0
4
0
0
0
-32
+0
4
0
0
4
0
0
4
16
0
4
2
0
0
0
0
4
0
-128
+0
0
0
4
4
4
0
0
0
-32
+0
0
4
2
0
2
0
4
2
4
4
8
16
0
4
4
4
4
0
4
0
16
0
4
4
0
0
0
4
0
0
0
4
0
4
0
4
4
4
4
0
4
4
0
0
0
0
4
4
4
0
0
8
4
2
4
4
4
4
0
4
4
0
0
2
4
16
0
4
2
0
0
4
4
0
0
4
4
4
0
0
0
-32
+0
4
0
4
0
0
0
0
4
0
0
-128
+0
0
4
0
0
0
4
0
0
4
0
2
0
0
0
4
0
0
-64
+0
4
0
4
0
4
4
4
0
0
4
4
0
0
4
4
2
0
4
0
0
4
4
0
4
4
0
0
0
4
0
2
16
0
2
4
4
2
-32
+0
0
0
4
0
2
4
-64
+0
0
0
0
0
4
4
4
2
4
0
0
-32
+0
0
4
0
0
0
0
4
0
0
0
4
0
4
4
0
-64
+0
4
4
4
0
-64
+0
0
4
4
0
4
4
0
0
4
2
4
4
0
0
4
0
0
4
0
0
4
16
0
4
0
4
4
0
2
4
4
0
4
0
4
0
0
4
8
0
0
0
4
0
4
4
4
8
8
-64
+0
0
4
0
4
4
4
16
0
2
-128
+0
4
4
0
0
4
0
0
0
-64
+0
0
4
0
-64
-128
+0
+0
2
8
4
0
2
0
0
0
4
4
4
4
0
0
4
2
0
0
4
4
0
2
0
4
4
4
0
4
4
4
0
0
4
4
0
0
0
4
4
4
2
0
4
4
4
0
4
4
4
4
0
0
0
4
4
4
0
4
0
4
0
0
4
0
0
0
-64
+0
0
4
4
2
0
0
0
4
0
4
0
4
4
4
0
16
4
0
4
8
0
2
4
4
4
4
2
4
0
0
4
4
0
-32
+0
0
0
0
0
4
4
0
4
0
0
4
-64
+0
4
0
4
0
0
0
4
-32
+0
0
0
4
4
8
4
4
0
0
0
0
-64
+0
4
4
0
0
2
-32
+0
2
0
4
4
4
0
-32
+0
2
4
4
4
2
4
4
2
2
0
0
4
0
4
0
4
4
4
0
0
4
4
4
0
2
4
2
4
2
0
4
4
2
0
4
2
0
4
0
0
4
8
0
2
4
0
16
0
4
0
4
4
2
0
2
0
0
4
0
0
0
0
4
4
4
4
2
-64
+0
4
0
4
2
0
0
0
4
4
4
4
4
4
0
0
4
4
4
4
0
-64
+0
4
4
0
0
0
0
0
2
4
0
4
4
4
0
0
4
0
0
4
2
0
4
4
8
2
0
0
4
4
16
4
0
0
0
4
4
0
0
-64
+0
4
4
4
16
4
4
0
4
0
0
2
0
4
4
0
4
4
0
4
0
4
4
0
4
0
4
4
0
0
0
0
0
0
0
4
-128
+0
0
4
0
0
-128
+0
4
-128
+0
0
0
0
0
0
16
16
4
0
0
4
4
4
4
2
4
0
0
4
0
4
0
0
4
-128
+0
4
0
0
4
4
4
4
2
4
-128
0
-64
+0
+0
0
4
16
0
0
0
4
0
0
4
0
4
0
4
8
4
0
4
4
0
4
0
4
0
0
2
4
4
4
0
0
4
0
4
4
0
0
0
4
4
2
0
0
4
0
0
0
4
4
4
0
4
0
4
0
4
16
4
8
4
0
4
0
4
4
2
0
4
4
0
4
0
0
4
0
4
4
0
4
0
4
0
4
4
4
4
4
-64
+0
4
2
4
0
0
2
-32
+0
0
4
-128
+0
4
0
4
4
4
2
-64
+0
0
4
4
0
4
4
0
2
0
0
4
0
4
0
0
0
4
4
4
0
0
0
4
4
0
0
0
4
0
0
0
0
0
0
0
0
0
2
0
4
0
4
0
4
0
4
0
0
0
4
0
0
0
4
0
0
0
0
4
0
0
8
4
0
0
0
4
4
4
4
0
4
4
4
4
4
0
4
4
4
0
0
0
4
4
4
4
8
4
4
4
4
0
4
8
4
0
-64
+0
0
2
0
0
0
0
4
0
0
0
4
4
4
4
4
0
4
4
0
4
4
4
0
0
0
4
4
0
0
0
-32
+0
0
4
4
4
0
0
0
4
0
-128
+0
4
4
4
4
4
2
16
-32
+0
4
4
0
0
0
2
0
-64
+0
4
0
4
0
0
4
4
4
4
2
4
16
4
0
0
4
2
0
4
4
0
4
0
4
4
0
4
0
4
4
4
4
0
4
0
4
0
0
4
4
4
0
4
0
-64
+0
4
4
4
0
4
0
2
4
4
0
4
0
4
4
4
4
0
0
0
0
4
4
4
4
0
2
4
0
8
0
4
0
0
0
0
0
4
4
4
0
0
0
0
0
2
0
4
4
4
16
4
4
4
4
4
0
0
0
0
4
4
4
-64
+0
4
4
0
2
0
4
4
4
0
4
2
0
4
4
4
0
0
2
0
4
0
0
4
0
0
4
4
4
4
4
2
4
0
4
0
4
4
0
-64
+0
4
2
4
4
0
0
4
0
0
0
0
4
4
8
0
0
2
4
4
0
0
0
4
0
4
0
4
0
4
4
2
4
4
-128
+0
4
0
4
4
0
4
2
0
4
4
4
4
2
0
0
4
0
0
0
4
4
0
4
0
4
4
0
0
4
4
2
4
0
4
4
0
0
2
0
0
0
8
0
0
4
4
0
0
0
0
0
0
4
0
2
4
0
0
4
0
0
4
0
4
4
4
4
8
0
0
0
0
4
4
4
0
4
4
0
4
4
0
8
4
4
2
2
0
0
0
4
0
0
4
4
4
4
0
0
4
0
4
4
0
4
4
4
4
8
0
2
4
0
-64
+0
2
4
4
4
0
0
0
0
4
0
0
4
0
0
2
4
4
4
0
0
4
4
0
2
-128
+0
4
0
0
0
4
8
4
4
0
4
0
0
4
0
4
0
4
0
4
4
16
0
0
4
0
0
0
0
16
0
4
0
0
4
4
0
4
4
4
4
0
4
-64
0
-64
+0
+0
16
4
0
4
4
0
4
0
0
0
4
0
0
2
4
0
4
2
0
4
0
4
4
4
0
0
4
4
4
4
4
2
0
16
0
0
4
4
4
4
0
4
0
4
4
2
0
0
4
4
2
-32
+0
0
4
0
-128
+0
0
0
0
2
-32
+0
4
4
0
4
-64
+0
4
0
2
4
0
4
0
0
4
0
0
4
0
4
0
0
4
0
2
4
4
4
0
0
4
0
4
4
4
0
4
2
0
4
4
0
0
0
-32
+0
0
0
4
4
4
0
0
4
0
4
4
0
4
0
0
0
4
0
0
-32
+0
4
4
0
4
4
2
0
4
4
4
-64
+0
16
0
0
4
4
4
0
4
0
0
4
2
0
2
4
0
4
0
4
0
0
4
0
0
0
0
4
0
0
4
0
4
0
4
0
0
0
0
4
4
4
4
4
0
0
0
4
0
4
0
2
0
4
0
0
-64
-64
-128
+0
+0
+0
0
0
4
0
4
0
16
4
0
4
2
0
0
0
0
4
4
0
0
0
0
0
4
4
4
4
4
2
16
4
0
0
0
0
-128
+0
4
4
0
0
4
4
4
4
4
4
4
4
0
4
0
4
4
16
-128
+0
4
0
0
4
4
4
4
4
0
4
4
4
2
4
2
4
2
4
0
4
4
0
4
4
0
-32
+0
4
4
16
0
2
0
4
4
2
4
4
4
0
16
4
0
4
4
4
0
0
0
0
0
4
0
4
0
4
0
4
4
4
4
0
0
4
4
0
16
0
0
0
0
4
0
16
0
0
0
4
0
0
0
4
4
4
0
0
0
0
0
4
4
0
4
0
4
4
0
2
0
2
0
0
0
0
0
4
0
0
4
4
4
0
0
0
4
0
4
4
0
4
4
4
0
4
0
2
0
0
4
0
4
0
0
0
0
4
0
4
4
4
0
0
4
4
4
0
4
2
0
4
4
-128
+0
2
4
4
4
16
0
4
0
4
0
0
4
4
2
4
-128
+0
2
0
0
4
2
0
0
0
4
0
4
4
0
0
0
0
0
0
2
4
0
4
0
4
8
0
0
4
0
0
4
0
0
4
16
4
4
0
0
0
4
0
0
0
0
4
16
0
4
4
0
4
0
0
4
4
4
4
0
0
4
0
4
0
0
4
4
-64
+0
4
4
0
4
8
0
0
0
4
4
4
0
4
0
0
4
0
0
8
16
2
0
4
0
4
0
4
4
0
0
4
4
0
4
0
8
2
4
0
0
4
4
4
4
0
0
0
0
0
0
4
4
0
16
4
0
4
2
0
4
0
4
0
4
0
0
0
4
4
0
4
4
0
4
0
-32
+0
4
0
0
0
0
4
0
0
4
4
0
0
4
4
0
4
2
4
0
0
0
4
-64
+0
4
-128
+0
0
4
4
8
0
4
8
0
-32
+0
4
4
16
2
4
4
0
4
0
0
4
16
4
4
0
0
4
0
4
4
4
8
4
0
4
4
4
2
4
2
4
4
0
4
4
0
4
0
4
4
0
2
0
4
16
4
4
4
0
0
0
4
4
0
4
4
0
4
4
4
4
4
4
0
4
4
4
4
4
2
0
4
4
0
0
4
0
4
4
4
0
2
0
4
0
4
4
0
0
2
0
4
0
4
0
4
4
0
4
4
-128
+0
2
0
4
4
4
4
-128
+0
0
0
8
4
-32
-128
+0
+0
4
16
2
0
4
4
2
4
4
-32
+0
4
2
4
0
0
4
0
0
0
0
4
0
0
0
4
4
4
4
0
0
4
4
0
4
0
4
0
4
4
4
4
0
-128
+0
16
-32
+0
4
8
0
4
0
0
0
0
4
0
-128
+0
4
0
0
0
4
4
4
0
0
-128
+0
4
4
4
0
0
4
0
2
4
0
0
4
4
0
4
0
0
4
0
-64
+0
0
2
0
4
0
0
4
4
4
4
2
4
0
4
-64
+0
4
4
0
0
0
0
4
4
0
4
0
4
0
0
0
4
4
4
4
4
4
4
4
4
4
4
0
4
4
0
4
16
0
4
4
16
4
4
4
4
4
4
4
4
2
4
0
0
0
4
0
0
4
0
4
4
4
4
0
0
0
0
0
4
8
0
4
2
4
4
4
4
0
4
0
2
4
4
4
4
4
0
4
4
0
4
4
4
0
0
0
4
-64
+0
0
4
0
0
2
4
4
4
4
8
4
diff --git a/t/classify_ref_4j b/t/classify_ref_W4j
similarity index 73%
copy from t/classify_ref_4j
copy to t/classify_ref_W4j
index 701338a..0e37a12 100644
--- a/t/classify_ref_4j
+++ b/t/classify_ref_W4j
@@ -1,3255 +1,4000 @@
4
-4
0
0
-4
0
4
4
+8
4
-16
-4
-0
-0
-0
4
0
0
4
0
+8
4
0
+8
+16
0
4
0
0
4
0
4
+32
+4
0
4
0
0
0
-64
4
-4
-2
0
-0
-2
4
0
0
4
-4
+128
4
4
0
4
0
+16
+128
0
-4
0
-2
-4
-4
0
0
0
4
0
-4
0
-4
-4
-64
-4
-4
0
4
-4
-4
0
-4
-4
0
0
+4
32
-0
+128
4
4
-64
-2
0
0
4
-0
-0
-0
+8
0
4
0
+8
+2
+8
0
4
-4
-0
-0
0
4
+64
0
+4
+32
0
0
+4
+4
0
+8
+4
+16
0
+64
0
0
4
-2
4
-0
-0
4
0
-0
-4
+16
0
4
4
0
0
0
4
0
-4
-4
-4
+0
4
64
0
+32
4
4
-0
-4
4
4
-2
0
+32
0
4
-4
0
+128
0
+4
0
4
0
+0
+4
+4
4
4
0
0
4
4
+32
2
4
+0
+0
2
-4
-4
+8
0
0
0
-4
+2
4
0
-4
0
-8
-4
-128
4
0
-4
-4
-2
0
-64
-2
2
-8
-4
4
-128
+0
+0
4
0
-8
4
+8
4
-16
8
0
+0
+0
4
0
0
+16
0
+2
+8
0
-4
0
0
+8
+64
+4
4
4
0
0
-4
+0
4
0
-2
0
0
-4
-4
0
4
4
-0
-0
4
4
4
+16
+0
+0
+0
0
-2
0
-4
-4
-4
4
0
+0
4
0
+16
0
0
+1
0
-4
0
-4
-4
-64
0
0
-4
0
+32
0
4
+4
0
4
+16
0
0
0
+128
+32
4
+8
+0
+0
4
+1
+8
0
4
+0
+0
+0
+32
4
4
+16
+0
+8
+0
0
4
0
0
4
+128
+0
+0
0
0
-64
0
-8
-4
-4
4
4
-0
+1
4
+16
8
+0
4
-8
4
4
+0
4
+0
4
0
+32
+0
4
4
4
+16
0
0
4
0
-0
4
0
-4
0
4
-16
0
0
+16
0
-2
+8
0
0
-4
0
+4
0
4
0
-64
-8
+4
+4
4
0
0
+2
+0
0
0
-32
4
+8
+0
4
4
-16
0
0
-4
-16
+0
4
4
+32
+0
+0
+0
+64
0
0
4
0
0
-2
0
0
-64
+16
+0
0
+128
0
0
0
+16
0
-4
0
-4
+8
+8
0
0
4
2
-4
-4
0
0
4
4
-4
-0
0
-4
0
0
-4
0
0
-4
+8
0
-4
0
0
-4
-4
-4
-4
0
0
-64
-4
-4
-4
-2
0
-2
0
0
0
4
+8
0
-4
-4
+8
0
4
-0
4
4
0
-4
+8
+0
0
-2
0
0
4
0
0
+8
+0
+0
0
-64
-4
0
-4
0
0
0
4
0
+0
4
4
+128
0
4
0
4
4
+4
16
0
-4
-4
0
0
-4
+16
0
-2
-4
0
-4
-4
+1
+0
+0
+0
+8
32
+0
4
0
4
-4
-2
-0
-4
0
-8
-8
4
0
0
0
-4
2
0
-4
0
-4
0
+0
+2
4
4
+0
4
+8
4
+1
0
-4
0
0
-2
0
-4
0
4
-0
4
+0
+16
+8
4
+32
0
+2
0
-4
-4
-4
0
4
-4
0
0
+16
4
4
4
4
0
-2
-0
0
0
-4
+8
2
4
+0
4
+0
+16
+0
4
-4
-4
-2
+0
+0
0
4
4
+128
0
-4
0
0
-4
-2
+16
0
4
0
+0
4
4
+0
4
4
0
+0
4
0
0
-4
0
0
0
0
+16
0
-64
+8
+0
+16
+4
+8
0
4
0
0
0
-4
-4
-4
-4
0
0
4
4
4
4
+32
4
4
-16
-4
+1
4
+8
4
0
-4
-4
-4
-4
0
-2
-4
4
0
16
0
+64
0
-128
+1
0
-2
+0
+128
4
+32
4
+0
4
-64
4
-0
-0
4
0
-0
+4
+4
4
0
+8
+4
+32
+4
+4
+2
+4
0
+8
+32
4
4
0
4
+0
+0
4
4
0
4
+0
4
4
4
+128
+16
+0
+8
+16
0
+8
+4
4
0
-64
4
0
+8
0
+4
+16
0
+8
+4
0
4
4
+0
4
-2
+0
4
0
0
-4
+8
+8
0
+4
0
0
4
+4
+4
0
0
0
0
4
-32
-128
-4
-2
+8
4
4
0
4
-0
-0
4
+4
+8
0
-16
0
4
0
+8
4
-128
-2
+0
8
4
0
0
0
+4
+4
+1
+0
0
0
-4
+8
+32
0
4
4
4
0
0
0
-0
+16
4
4
0
0
-0
4
0
0
0
-2
-4
+8
4
0
0
4
-2
0
-0
-4
4
4
-64
4
+0
+8
+0
4
0
0
-2
4
+0
4
+2
+0
0
0
0
0
+1
+0
+8
4
+0
4
4
-2
0
-32
0
+8
0
4
4
4
-4
-2
-4
-4
0
+8
0
0
64
0
+4
0
0
+16
+0
0
-4
0
-2
-2
-4
-4
-4
-4
-32
0
0
+0
+0
+0
+0
+16
4
0
0
+4
0
4
4
0
0
-8
+128
4
-2
+0
+16
4
+128
+0
+32
0
4
0
+0
+0
4
2
+64
0
+16
4
+4
+0
0
0
0
4
+8
4
0
+0
4
0
4
0
+0
+0
+0
+0
+0
+0
4
0
4
+1
4
-64
4
4
+0
4
+0
4
4
4
-2
-2
-2
-0
-0
-0
+8
+16
4
0
4
-0
4
0
0
4
4
-4
-4
8
0
-8
+128
0
0
0
0
-4
-4
0
-2
-2
0
-2
-4
-32
0
4
+2
+0
4
4
0
4
-2
-0
+8
4
-0
4
0
-0
4
0
-4
0
4
0
4
0
4
4
-4
0
+16
0
+16
0
4
4
+16
4
4
-0
-0
+64
0
0
4
-0
+2
+1
4
0
+4
+4
0
-2
4
4
0
0
-2
-4
0
-2
-8
4
4
-0
-0
+4
4
0
+16
0
4
0
4
-4
-4
-4
32
-4
-4
-4
-4
-4
-0
-0
-4
0
+64
0
4
-2
4
0
4
-0
-2
-0
+16
0
0
4
4
-4
-4
-0
-4
-0
+8
+32
0
-4
-4
+2
0
-4
-4
-8
0
4
+64
0
4
0
4
4
-128
-4
0
0
0
-2
-0
-16
-0
0
0
-4
-4
0
0
4
0
0
-4
0
2
0
0
-4
0
4
0
-4
0
-4
-64
-4
-16
-4
-4
-4
0
-2
+4
2
0
0
-8
-4
-0
4
0
0
-4
-4
-4
-0
-64
-4
-64
32
0
+8
0
-4
-32
0
0
4
0
8
0
+8
0
-0
-4
-4
2
-4
0
-2
4
-0
-2
+4
0
0
0
-4
0
4
-4
-4
-4
+8
0
0
64
-4
-0
-4
0
+128
4
4
-2
4
+0
+0
4
4
4
0
-0
-4
4
-128
0
2
0
-0
-0
4
-0
-0
+8
0
0
4
4
-0
+16
4
4
0
0
-4
-4
0
0
0
+64
+8
+4
0
+2
0
+8
0
+4
0
4
+2
0
+8
+8
+4
0
-32
0
4
0
0
-0
+4
4
0
0
-4
-128
0
0
4
+2
0
0
4
+4
+0
0
0
0
0
+0
+4
+4
64
4
-0
4
0
-4
2
0
4
4
4
0
-0
-0
4
+16
4
-32
0
+0
+8
2
0
4
+2
4
-16
+8
4
4
4
-0
-0
+16
0
4
0
0
0
-4
0
+4
+16
+4
+32
0
0
-32
+16
4
0
+16
0
4
0
+4
+4
+128
+4
0
+8
+8
4
-16
0
4
2
+4
0
-0
+4
+4
+16
+4
+4
0
0
4
-0
-128
+16
0
0
4
4
+0
+128
4
0
+4
0
+4
0
-32
0
+16
+16
+8
4
-2
0
-2
+8
+0
+0
0
4
+0
2
4
-4
-8
-16
0
-4
-4
-4
-4
0
-4
0
16
0
4
4
0
-0
+4
0
4
0
0
0
-4
0
-4
0
4
4
4
-4
+8
0
+16
+0
+1
4
-4
+16
0
0
0
0
4
+0
+32
4
4
0
-0
8
4
+32
+0
+8
+0
2
+0
4
-4
+0
+0
+0
+16
+0
+1
4
4
0
4
+64
4
-0
-0
2
-4
-16
0
4
-2
0
0
4
4
+4
+0
0
0
-4
-4
4
0
0
0
-32
+16
+0
+4
4
0
+0
+64
4
+8
0
+8
+8
0
0
0
4
+4
+64
+2
0
+8
0
-128
+2
+16
+16
+8
0
4
+2
+32
0
0
0
-4
0
0
-4
0
+0
+8
2
0
0
+8
+128
+0
0
4
0
0
-64
-4
0
-4
0
-4
-4
-4
0
0
+16
4
+0
4
0
0
-4
-4
-2
+8
0
-4
+16
0
+128
0
4
-4
0
4
-4
0
+4
+2
+64
0
0
4
0
-2
+16
16
0
-2
-4
-4
+8
+8
+64
+0
2
32
-0
-0
4
0
2
-4
-64
0
+4
+8
0
+8
+8
0
0
4
-4
-4
-2
-4
-0
-0
-32
0
4
0
+4
0
+4
0
0
4
+4
0
0
0
-4
0
4
4
-0
-64
4
+0
4
4
+8
+0
0
-64
0
-4
-4
0
-4
-4
0
0
4
2
+0
4
-4
+2
0
+16
0
-4
+0
+64
+32
0
0
+16
+4
4
0
+64
0
-4
16
0
-4
+2
0
4
-4
0
-2
-4
-4
0
-4
+8
+0
0
4
0
0
4
-8
0
0
0
+32
4
0
4
+0
+0
4
+0
4
+0
+0
8
-8
-64
0
-4
0
4
+32
+8
4
-4
-16
0
-2
-128
4
4
0
0
-4
-0
-0
0
-64
0
-4
0
-64
-128
-2
8
-4
-0
-2
-0
0
-0
-4
4
4
+0
4
0
+4
0
4
-2
0
0
+16
+0
4
4
0
-2
0
4
-4
-4
0
4
4
+16
4
0
-0
-4
4
+32
+8
0
0
+32
0
4
-4
-4
-2
+0
0
4
4
-4
+16
0
+0
+0
+16
4
4
4
+0
4
0
+8
0
+8
0
+128
+0
+8
4
-4
+1
4
0
+8
4
-0
4
-0
-0
4
0
0
+4
+8
0
-64
0
4
4
+0
+0
2
0
0
0
-4
0
-4
0
4
-4
-4
0
+1
16
-4
0
-4
-8
0
-2
-4
+0
+0
+0
+0
+0
4
4
+32
4
2
4
0
-0
4
4
0
-32
0
+16
+4
0
0
+16
+4
+0
+16
0
4
4
0
4
0
+4
0
4
-64
4
+4
+16
+0
+4
+4
+64
0
4
0
0
0
-4
-32
0
0
4
-4
8
4
4
+4
0
0
0
-0
-64
4
4
+8
0
0
-2
-32
-2
0
4
4
-4
+8
+64
0
-32
2
-4
-4
-4
-2
-4
-4
-2
-2
-0
-0
-4
-0
-4
+16
0
+2
4
4
4
+16
0
0
-4
-4
-4
0
-2
-4
-2
4
2
0
4
-4
-2
0
-4
2
+8
+16
+0
0
+16
+4
+1
+128
4
0
0
4
8
+16
0
+128
2
4
+128
0
-16
+8
0
-4
0
4
+64
4
-2
-0
-2
-0
-0
4
0
+4
0
+4
0
0
4
4
4
4
-2
-64
-4
-0
-4
-2
0
0
0
4
+0
+64
4
-4
+0
4
4
4
0
0
-4
-4
-4
-4
0
-64
+0
+16
+0
4
4
0
0
+4
0
0
0
-2
-4
0
-4
-4
-4
+8
0
0
4
0
+64
+4
+8
0
4
-2
0
+8
4
4
-8
-2
+128
+0
0
0
4
4
-16
+4
4
0
+16
+2
0
+1
0
+16
4
4
+16
0
+16
0
-64
-4
4
4
16
-4
-4
0
-4
+2
+64
0
0
-2
0
4
4
0
4
+0
+0
4
0
+0
+0
+0
4
0
4
4
-0
+32
4
-0
4
4
+8
0
0
-0
-0
-0
-0
+32
+4
0
4
-128
0
4
0
0
-128
-4
-128
0
0
0
+2
+0
0
0
-16
-16
4
0
0
4
-4
-4
-4
-2
-4
0
0
4
-0
4
0
0
-4
128
-4
0
+64
0
4
+0
+16
+0
4
-4
-4
-2
-4
-128
0
64
0
-4
-16
0
+8
0
0
+16
+4
4
0
0
-4
0
+16
4
0
-4
8
4
-0
-4
4
0
+0
+0
+16
+0
4
0
4
0
0
-2
+0
+64
+0
+16
4
4
4
0
-0
+1
+8
4
-0
4
4
0
0
0
-4
-4
+0
2
+16
0
0
+128
4
0
-0
-0
-4
4
4
0
+16
4
+16
+128
+8
+8
0
-4
-0
-4
16
4
-8
4
0
4
0
4
4
-2
0
-4
-4
0
-4
0
0
-4
+8
+0
+0
0
4
4
0
+2
+8
4
0
4
0
4
-4
-4
-4
-4
-64
-4
2
+0
4
0
+8
0
-2
-32
0
4
-128
-4
0
-4
-4
-4
-2
-64
0
-4
-4
0
4
-4
+16
+8
0
-2
0
0
-4
0
4
+8
+0
0
0
0
4
4
4
0
0
+32
0
+8
+16
4
-4
-0
0
+64
+4
0
4
0
+4
+2
0
+128
0
0
0
0
0
+4
+4
0
+4
0
-2
0
4
+1
0
-4
0
+0
+0
+8
4
+32
0
4
+4
+32
+4
+64
+0
0
0
0
-4
0
0
0
+32
+4
4
+4
+32
0
+4
0
+16
0
+2
+2
0
4
+32
0
0
-8
-4
0
0
+16
0
4
4
+0
+0
4
4
0
4
4
4
4
-4
+16
0
+32
4
4
4
0
0
0
-4
-4
-4
-4
-8
-4
-4
-4
-4
-0
-4
-8
-4
+1
0
-64
0
-2
0
0
0
0
4
0
0
0
+8
+0
+0
+16
+8
4
+1
+8
+16
+0
4
+128
+64
+0
4
-4
+0
+8
4
0
4
4
+4
+4
+2
0
+16
4
4
4
+4
+0
0
+2
0
0
4
+0
+8
+16
+1
4
0
0
+8
+0
0
-32
0
4
-4
+0
+0
+8
4
0
0
0
-4
0
-128
-4
-4
-4
-4
-4
-2
16
-32
-4
-4
0
0
+8
+0
0
-2
0
64
-4
0
-4
0
0
+0
+0
+1
4
4
-4
-4
+0
+0
+64
+0
+16
2
+0
+32
+4
+0
4
-16
4
0
0
+8
4
+0
+8
2
0
+0
+4
4
4
0
+0
4
0
4
4
+8
0
-4
+8
0
4
-4
-4
-4
0
-4
+2
+0
0
-4
0
0
4
-4
+0
+0
4
0
4
0
-64
+0
4
4
4
0
-4
0
-2
-4
-4
0
4
0
4
4
-4
-4
+1
+0
0
0
0
0
4
+0
+16
4
+0
4
+16
4
-0
-2
4
0
-8
0
-4
+16
0
0
0
0
0
-4
-4
-4
+8
0
0
0
+8
0
0
-2
+4
0
4
+0
+0
4
+8
4
+0
+0
16
4
4
4
-4
-4
0
+4
0
0
0
4
4
-4
-64
-4
-4
0
-2
0
4
-4
-4
0
-4
-2
0
-4
-4
+8
4
0
0
-2
0
-4
0
0
-4
0
0
-4
-4
-4
-4
+16
+0
+0
4
2
-4
0
-4
0
-4
-4
0
-64
+0
4
-2
4
4
0
-0
4
0
0
0
-0
-4
+16
4
-8
-0
0
-2
-4
-4
0
0
0
-4
0
-4
+8
0
-4
0
4
+64
+8
4
+1
+128
2
+0
4
4
-128
-4
-0
4
4
0
4
-2
-0
+128
4
4
4
4
-2
+0
0
0
4
0
0
0
+0
+8
+0
4
4
+4
+4
+0
+64
0
4
0
4
4
+8
+2
+0
0
0
4
+0
+0
+0
+0
+0
4
-2
4
0
4
+16
4
+8
+0
+1
+8
+0
0
0
-2
0
0
0
-8
+0
+0
+0
+16
0
0
4
+16
+0
+0
+0
+128
4
0
0
+16
+8
+16
+0
+0
0
0
0
0
4
+4
+4
0
2
+16
+0
+4
4
0
0
-4
0
0
-4
0
+16
+32
+16
4
+0
+8
4
+0
+0
4
4
-8
0
0
0
0
4
-4
-4
0
4
+128
+128
+4
4
0
+32
4
-4
+2
0
8
+8
4
4
-2
-2
-0
0
0
4
+4
0
0
+0
+8
+0
4
+0
+8
4
-4
+0
4
0
0
+16
4
0
-4
-4
0
4
4
+0
+0
+0
+0
4
4
-8
0
-2
4
-0
-64
-2
4
4
4
-0
+128
+4
0
0
0
4
+2
0
+4
+2
+1
0
4
0
0
+32
2
4
+8
4
4
-0
-0
+8
4
4
-0
-2
-128
+4
+8
4
0
-0
+32
+4
0
4
8
+64
4
4
0
+0
+32
4
0
+8
+2
0
-4
+32
0
-4
0
4
-0
4
4
-16
+4
+8
0
0
4
0
0
+4
+0
0
0
16
0
+16
+0
+0
+4
+0
+1
+8
4
0
0
4
+0
4
0
+0
+8
+0
+128
+0
+8
4
+0
4
4
4
0
4
-64
+4
+0
0
64
-16
4
-0
4
4
0
+8
+0
+0
4
0
0
0
-4
0
0
-2
-4
+0
0
4
-2
0
+0
+16
4
0
4
4
+0
+8
4
0
0
4
+8
+8
4
-4
-4
-4
+0
2
0
-16
0
0
-4
-4
-4
-4
0
4
+8
+1
+0
0
4
-4
-2
0
0
+64
+0
4
4
-2
-32
-0
4
0
+0
+0
+0
+0
128
0
0
0
-2
-32
4
4
+128
0
4
-64
-4
0
-2
-4
0
-4
0
0
4
0
0
-4
0
4
0
-0
4
0
-2
4
+128
4
+0
4
0
0
-4
+8
0
4
4
4
0
+0
+4
4
-2
0
+8
4
4
0
0
+4
+0
0
-32
0
0
-4
4
4
0
-0
4
0
4
4
0
-4
0
0
0
4
-0
-0
-32
4
4
0
4
-4
+64
+0
+0
2
+4
+4
0
4
4
4
-64
-16
+4
+2
0
0
+32
4
+0
4
4
0
-4
+0
+0
+32
0
0
4
-2
+4
0
-2
+4
+4
+4
4
0
4
+8
+4
+8
+4
0
4
+8
+0
0
0
-4
0
0
+8
+4
+4
0
0
+64
4
+4
+8
0
0
+0
+128
+4
4
0
4
+2
+8
+0
0
4
0
+16
+64
0
0
+4
+2
+4
0
+128
4
4
+16
+4
+16
4
4
4
0
0
0
4
0
4
0
+4
+32
2
0
-4
+16
0
+16
0
-64
-64
128
-0
-0
-4
-0
4
+16
0
+8
16
4
0
+8
+0
+8
+0
4
-2
0
0
+2
+2
+2
+0
0
+8
0
4
-4
-0
0
0
0
0
+32
4
-4
-4
-4
-4
+0
2
-16
-4
0
0
+32
+0
+4
+8
+0
0
0
-128
4
4
0
0
4
4
4
4
+0
4
4
-4
+8
4
0
-4
0
4
-4
-16
-128
-4
0
0
-4
-4
-4
-4
-4
0
-4
-4
-4
-2
-4
-2
-4
-2
-4
0
-4
-4
0
4
4
0
32
4
-4
16
+128
0
-2
-0
-4
-4
-2
-4
-4
4
0
+0
16
4
0
4
-4
-4
-0
-0
+8
0
+16
+4
0
0
-4
0
4
0
-4
0
4
-4
-4
-4
-0
0
4
+64
+64
+1
4
0
16
-0
-0
-0
+64
0
4
0
-16
0
+16
0
0
4
0
+16
0
-0
-4
4
4
0
+64
+0
0
0
+4
0
+2
0
+128
4
4
-0
+1
+128
+8
4
0
4
-4
0
-2
0
-2
0
0
+4
0
0
0
-4
+16
+8
0
0
-4
-4
+64
4
0
0
+64
0
4
0
4
-4
0
4
-4
+16
4
0
4
-0
+64
+4
+4
2
0
0
-4
0
4
0
0
0
0
-4
0
+64
4
-4
-4
-0
0
4
4
4
0
4
-2
-0
-4
-4
-128
-2
-4
-4
4
-16
0
4
-0
4
0
0
+0
4
4
-2
-4
-128
-2
-0
0
-4
-2
+32
0
0
0
4
-0
4
+0
4
0
+4
0
0
0
+32
0
0
-2
4
-0
4
0
-4
8
-0
-0
4
-0
-0
4
0
0
4
16
4
-4
-0
-0
0
4
-0
-0
+4
+4
0
0
4
-16
-0
4
4
-0
4
0
-0
-4
-4
-4
4
-0
-0
4
-0
4
-0
-0
4
4
-64
4
4
+16
0
4
8
+4
0
0
0
4
-4
-4
+8
0
-4
+8
0
+128
+16
0
-4
0
0
+4
+4
8
-16
-2
0
-4
0
+32
+0
+4
4
0
+16
+4
4
4
0
0
4
4
+4
+4
0
4
0
8
-2
4
0
-0
-4
4
+0
+8
4
+0
+8
4
0
+2
+4
0
+64
0
0
0
0
4
4
0
-16
-4
0
4
-2
0
-4
0
4
+64
+4
+8
+4
+4
+4
+0
0
4
0
0
0
4
4
-0
+2
4
+1
+16
+0
4
+1
+2
0
+64
4
+64
+1
0
32
4
0
0
0
0
+1
4
-0
-0
4
+0
4
0
0
-4
-4
0
-4
-2
-4
0
+4
0
+8
+2
0
4
-64
4
-128
-0
4
4
8
-0
4
-8
0
-32
4
4
+8
+0
16
-2
-4
-4
0
-4
+0
+0
0
0
4
16
-4
-4
0
0
4
+64
+32
0
4
-4
-4
-8
-4
0
4
-4
-4
-2
+16
+8
4
2
4
-4
0
+16
4
-4
+0
0
4
0
4
4
+32
+32
0
-2
0
4
-16
-4
+0
4
+0
+0
4
0
0
0
4
-4
0
+2
4
-4
+64
0
-4
-4
-4
-4
-4
-4
0
4
+0
+0
+0
4
4
-4
-4
+16
2
0
-4
-4
+64
+32
0
0
-4
0
-4
-4
-4
0
2
0
-4
+8
0
4
-4
0
0
-2
0
-4
0
+16
4
0
4
4
0
+0
+8
+0
4
+0
+0
4
-128
-2
+0
0
4
+128
4
+0
4
4
-128
-0
-0
8
4
-32
-128
4
-16
-2
0
+1
4
4
-2
-4
4
+8
32
4
-2
4
0
0
-4
+8
0
0
0
0
+2
4
0
0
0
4
-4
+0
4
4
0
0
-4
-4
+32
0
-4
0
-4
+16
+0
0
4
+2
4
+2
+1
+0
+0
4
4
+8
0
+4
128
+2
16
-32
4
-8
-0
4
0
-0
-0
-0
4
0
-128
-4
0
0
+128
+8
0
-4
-4
-4
0
0
-128
4
4
+8
4
0
0
4
+8
0
-2
4
0
0
-4
+16
4
0
-4
0
0
-4
0
64
+8
0
-2
+0
+32
0
4
0
0
-4
-4
-4
-4
-2
-4
0
+8
4
-64
+32
+0
+0
4
4
+16
0
+16
0
0
0
-4
-4
+1
0
+16
4
0
-4
0
0
+2
0
4
4
4
4
4
-4
-4
-4
-4
-4
-4
0
4
-4
0
-4
-16
0
-4
-4
-16
-4
-4
-4
-4
-4
-4
-4
-4
2
-4
0
0
0
4
0
-0
4
+8
0
+2
4
4
+8
4
-4
0
0
0
0
+16
+128
+0
0
4
+4
+0
+0
+4
+0
+4
+0
+0
+16
+0
+4
+4
+4
+4
+4
+4
+0
+8
+4
+0
+0
+0
+0
+8
+0
+0
+0
+16
+0
+4
+32
+0
+4
+4
+4
+4
+0
+8
+0
+0
+0
+0
+4
+4
+64
+4
+4
+4
+4
+0
+0
+0
+4
+0
+4
+4
+0
+8
+1
+4
+0
+0
+0
+8
+0
+0
+0
+0
+4
+0
+4
+4
+8
+16
+8
+8
+0
+0
+64
+8
+4
+8
+4
+0
+4
+0
+4
+4
+16
+0
+32
+4
+0
+8
+4
+0
+0
+0
+0
+16
+0
+4
+8
+64
+4
+8
+4
+4
+0
+4
+0
+32
+4
+4
+4
+0
+1
+8
+0
+4
+0
+0
8
+4
+2
0
4
+64
+8
+4
+2
+0
+0
2
4
+0
+0
+0
+0
+4
+0
4
4
4
0
+0
4
+8
+0
+128
+0
+0
0
2
+8
+0
4
4
+0
4
+0
4
+0
4
+8
+0
+0
0
4
+0
+8
4
0
4
+0
+0
4
+0
+8
+0
+128
+8
+0
+64
+0
4
0
0
0
+0
+4
+4
+0
+8
+4
4
+0
+0
+8
+0
64
+2
+4
+8
+4
+16
+4
+2
+0
+8
+0
+0
+0
+4
+4
+4
+4
+0
+16
+8
+16
+0
+0
+8
+4
+0
+4
+128
+0
+0
0
4
0
+4
0
2
+0
+0
+16
+0
+4
+0
+0
+4
+4
+0
+0
+4
+4
+64
+4
+4
+0
+0
+4
+0
+4
+0
+16
+8
+0
+4
+0
+0
+4
+4
+4
+0
+0
+4
+0
+0
+16
+0
+0
+0
4
+0
+0
+4
+4
+4
+32
+0
+0
+0
+8
4
+0
4
+16
+0
+0
4
8
4
+0
+0
+0
+0
+0
+8
+0
+4
+0
+0
+4
+4
+0
+16
+4
+8
+4
+0
+4
+4
+128
+64
+0
+0
+8
+0
+0
+0
+32
+0
+0
+0
+4
+0
+0
+0
+0
+0
+0
+0
+0
+4
+4
+4
+64
+0
+4
+16
+4
+0
+2
+0
+16
+0
+0
+0
+16
+16
+0
+0
+0
+16
+0
+0
+4
+4
+0
+0
+0
+0
+0
+4
+0
+0
+1
+32
+16
+16
+0
+4
+64
+1
+0
+16
+32
+0
+4
+0
+8
+0
+2
+4
+0
+0
+0
+4
+0
+0
+0
+4
+4
+0
+4
+0
+8
+0
+0
+4
+0
+0
+0
+4
+0
+0
+4
+64
+16
+0
+16
+16
+32
+32
+0
+1
+0
+4
+0
+0
+0
+4
+0
+4
+64
+4
+0
+0
+4
+4
+4
+0
+8
+0
+0
+4
+0
+0
+0
+0
+0
+4
+0
+0
+64
+0
+4
+8
+0
+0
+4
+4
+0
+4
+0
+4
+4
+4
+4
+0
+128
+4
+0
+4
+0
+0
+4
+0
+0
+128
+0
+0
+32
+128
+4
+0
+0
+0
+4
+0
+0
+0
+0
+16
+4
+16
+0
+0
+0
+4
+4
+4
+16
+4
+0
+4
+4
+1
+4
+0
+4
+0
+0
+4
+4
+4
+0
+4
+0
+0
+0
+4
+4
+0
+0
+0
+16
+4
+0
+4
+16
+0
+4
+8
+0
+0
+4
+0
+0
+4
+0
+4
+4
+4
+8
+0
+0
+0
+16
+8
+0
+0
+16
+16
+4
+0
+4
+0
+0
+4
+4
+16
+0
+16
+0
+4
+0
+4
+0
+4
+0
+0
+2
+16
+0
+0
+0
+64
+4
+0
+32
+0
+64
+0
+4
+0
+128
+0
+0
+0
+0
+16
+0
+8
+32
+0
+64
+128
+4
+8
+0
+0
+8
+4
+0
+0
+2
+0
+0
+0
+4
+0
+4
+4
+16
+0
+4
+32
+0
+4
+4
+0
+0
+16
+0
+4
+4
+0
+4
+0
+0
+0
+0
+4
+0
+0
+4
+8
+4
+0
+0
+0
+0
+0
+0
+0
+4
+2
+4
+4
+0
+0
+0
+4
+128
+4
+0
+4
+0
+0
+16
+8
+0
+4
+0
+4
+0
+0
+0
+4
+0
+4
+4
+4
+64
+0
+4
+16
+0
+32
+4
+0
+0
+0
+0
+0
+0
+0
+4
+0
+0
+4
+0
+0
+0
+16
+0
+0
+4
+32
+0
+0
+2
+0
+0
+0
+0
+0
+8
+4
+4
+0
+4
+0
+8
+0
+0
+0
+0
+4
+4
+4
+0
+0
+4
+4
+4
+0
+16
+4
+0
+0
+32
+0
+4
+0
+16
+4
+4
+2
+16
+4
+2
+8
+0
+0
+0
+0
+8
+0
+8
+4
+0
+0
+0
+4
+4
+0
+0
+0
+4
+64
+4
+4
+0
+64
+0
+4
+0
+4
+0
+4
+4
+0
+4
+0
+0
+0
+4
+0
+0
+128
+0
+16
+128
+0
+4
+0
+4
+8
+0
+4
+4
+16
+4
+0
+4
+4
+64
+4
+64
diff --git a/t/test_ME_generic.cc b/t/test_ME_generic.cc
index e22b945..9e94c9d 100644
--- a/t/test_ME_generic.cc
+++ b/t/test_ME_generic.cc
@@ -1,199 +1,202 @@
/**
* \brief Generic tester for the ME for a given set of PSP
*
* \note reference weights and PSP (as LHE file) have to be given as
* _individual_ files
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <algorithm>
#include <cmath>
#include <fstream>
#include <random>
#include "HEJ/Event.hh"
#include "HEJ/EventReader.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/stream.hh"
#include "HEJ/YAMLreader.hh"
constexpr double alpha_s = 0.118;
constexpr double ep = 1e-5;
constexpr double ep_mirror = 1e-3;
void shuffle_particles(HEJ::Event::EventData & ev) {
static std::mt19937_64 ran{0};
// incoming
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
// outgoing (through index)
auto old_outgoing = std::move(ev.outgoing);
std::vector<size_t> idx(old_outgoing.size());
std::iota(idx.begin(), idx.end(), 0);
std::shuffle(begin(idx), end(idx), ran);
ev.outgoing.clear();
ev.outgoing.reserve(old_outgoing.size());
for(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(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(begin(decay.second), end(decay.second), ran);
+ }
}
}
enum MEComponent {tree, virt};
MEComponent guess_component(std::string const & data_file) {
if(data_file.find("virt") != data_file.npos) return MEComponent::virt;
return MEComponent::tree;
}
HEJ::Event::EventData mirror_event(HEJ::Event::EventData ev){
for(auto & part: ev.incoming){
auto & p{ part.p };
p.reset(p.px(),p.py(),-p.pz(),p.E());
}
for(auto & part: ev.outgoing){
auto & p{ part.p };
p.reset(p.px(),p.py(),-p.pz(),p.E());
}
for(auto & decay: ev.decays){
for(auto & part: decay.second){
auto & p{ part.p };
p.reset(p.px(),p.py(),-p.pz(),p.E());
}
}
return ev;
}
int main(int argn, char** argv){
if(argn != 4 && argn != 5){
std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml ME_weights input_file.lhe\n\n";
return EXIT_FAILURE;
}
bool OUTPUT_MODE = false;
if(argn == 5 && std::string("OUTPUT")==std::string(argv[4]))
OUTPUT_MODE = true;
const HEJ::Config config = HEJ::load_config(argv[1]);
std::fstream wgt_file;
if ( OUTPUT_MODE ) {
std::cout << "_______________________USING OUTPUT MODE!_______________________" << std::endl;
wgt_file.open(argv[2], std::fstream::out);
wgt_file.precision(10);
} else {
wgt_file.open(argv[2], std::fstream::in);
}
auto reader{ HEJ::make_reader(argv[3])};
const auto component = guess_component(argv[2]);
HEJ::MatrixElement ME{
[](double){ return alpha_s; },
HEJ::to_MatrixElementConfig(config)
};
double max_ratio = 0.;
size_t idx_max_ratio = 0;
HEJ::Event ev_max_ratio(HEJ::Event::EventData{}.cluster(
config.resummation_jets.def,0
)
);
double av_ratio = 0;
size_t i = 0;
while(reader->read_event()){
++i;
HEJ::Event::EventData data{reader->hepeup()};
shuffle_particles(data);
HEJ::Event::EventData data_mirror{mirror_event(data)};
shuffle_particles(data_mirror);
HEJ::Event event{
data.cluster(
config.resummation_jets.def,
config.resummation_jets.min_pt
)
};
HEJ::Event event_mirror{
data_mirror.cluster(
config.resummation_jets.def,
config.resummation_jets.min_pt
)
};
const double our_ME = (component == MEComponent::tree)?
ME.tree(event).central:
ME.virtual_corrections(event).central
;
if(!std::isfinite(our_ME)){
std::cerr << "Found non-finite ME ("<< our_ME <<")\n" << event << std::endl;
return EXIT_FAILURE;
}
const double ME_mirror = (component == MEComponent::tree)?
ME.tree(event_mirror).central:
ME.virtual_corrections(event_mirror).central
;
if(!std::isfinite(ME_mirror)){
std::cerr << "Found non-finite ME ("<< ME_mirror <<")\n" << event_mirror << std::endl;
return EXIT_FAILURE;
}
if(std::abs(our_ME/ME_mirror-1.)>ep_mirror){
size_t precision(std::cout.precision());
std::cerr.precision(16);
std::cerr<< "z-Mirrored ME gives different result " << i << "\n"
<<our_ME << " vs " << ME_mirror << " => difference: "
<< std::abs(our_ME/ME_mirror-1.) << "\n" << event
<< "\nmirrored:\n" << event_mirror << std::endl;
std::cerr.precision(precision);
return EXIT_FAILURE;
}
if ( OUTPUT_MODE ) {
wgt_file << our_ME << std::endl;
} else {
std::string line;
if(!std::getline(wgt_file,line)) break;
const double ref_ME = std::stod(line);
const double diff = std::abs(our_ME/ref_ME-1.);
av_ratio+=diff;
if( diff > max_ratio ) {
max_ratio = diff;
idx_max_ratio = i;
ev_max_ratio = event;
}
if( diff > ep ){
size_t precision(std::cout.precision());
std::cerr.precision(16);
std::cerr<< "Large difference in PSP " << i << "\nis: "<<our_ME
<< " should: " << ref_ME << " => difference: " << diff << "\n"
<< event << std::endl;
std::cerr.precision(precision);
return EXIT_FAILURE;
}
}
}
wgt_file.close();
if ( i<100 )
throw std::invalid_argument{"Not enough PSP tested"};
if ( !OUTPUT_MODE ) {
size_t precision(std::cout.precision());
std::cout.precision(16);
std::cout << "Avg ratio after " << i << " PSP: " << av_ratio/i << std::endl;
std::cout << "maximal ratio at " << idx_max_ratio << ": " << max_ratio << std::endl;
std::cout.precision(precision);
}
return EXIT_SUCCESS;
}
diff --git a/t/test_classify.cc b/t/test_classify.cc
index ede2728..eddd9f8 100644
--- a/t/test_classify.cc
+++ b/t/test_classify.cc
@@ -1,983 +1,1128 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <iostream>
#include <random>
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#define ASSERT(x) if(!(x)) { \
throw std::logic_error("Assertion '" #x "' failed."); \
}
namespace {
const fastjet::JetDefinition jet_def{fastjet::JetAlgorithm::antikt_algorithm, 0.4};
const double min_jet_pt{30.};
- const std::array<std::string, 6> all_quarks{"-4","-1","1","2","3","4"};
- const std::array<std::string, 7> all_partons{"g","-2","-1","1","2","3","4"};
- const std::array<std::string, 3> all_bosons{"h", "Wp", "Wm"};
+ 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};
void shuffle_particles(HEJ::Event::EventData & ev) {
// incoming
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
// outgoing (through index)
auto old_outgoing = std::move(ev.outgoing);
std::vector<size_t> idx(old_outgoing.size());
std::iota(idx.begin(), idx.end(), 0);
std::shuffle(begin(idx), end(idx), ran);
ev.outgoing.clear();
ev.outgoing.reserve(old_outgoing.size());
for(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(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(begin(decay.second), end(decay.second), ran);
+ }
}
}
// if pos_boson = -1 (or not implemented) -> no boson
// njet==7 is special: has less jets, i.e. multiple parton in one jet,
+ // all partons are massive (4 GeV) -> can be boson/decay
// pos_boson < 0 to select process (see list for details)
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, { -11, -18, -42, 47}, {}});
- ev.outgoing.push_back({gluon, { -15, 26, -18, 35}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, -24, 113}, {}});
- ev.outgoing.push_back({gluon, { 23, -54, -6, 59}, {}});
- ev.outgoing.push_back({gluon, { -5, -44, 8, 45}, {}});
- ev.outgoing.push_back({gluon, { -48, -96, 44, 116}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, 88, 133}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -249, 249}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 299, 299}, {}};
+ 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, { 23, -80, -64, 105}, {}});
- ev.outgoing.push_back({gluon, { -5, -84, -12, 85}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, 24, 113}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, 88, 133}, {}});
- ev.outgoing.push_back({gluon, { -48, -24, 62, 82}, {}});
- ev.outgoing.push_back({gluon, { -11, -18, 42, 47}, {}});
- ev.outgoing.push_back({gluon, { -15, 20, 60, 65}, {}});
+ 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, 415, 415}, {}};
+ ev.incoming[1] = {gluon, { 0, 0, 397, 397}, {}};
return ev;
case -3: // jet idx: 0 0 1 2 3 4 5
- ev.outgoing.push_back({gluon, { -5, -86, -70, 111}, {}});
- ev.outgoing.push_back({gluon, { -15, -52, -36, 65}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, -44, 109}, {}});
- ev.outgoing.push_back({gluon, { -48, -60, 5, 77}, {}});
- ev.outgoing.push_back({gluon, { -11, 92, 8, 93}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, 24, 113}, {}});
- ev.outgoing.push_back({gluon, { 23, -80, 64, 105}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -361, 361}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 312, 312}, {}};
+ ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}});
+ ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
+ ev.outgoing.push_back({gluon, { -15, -58, -34, 69}, {}});
+ ev.outgoing.push_back({gluon, { 23, -76, -40, 89}, {}});
+ ev.outgoing.push_back({gluon, { -5, -52, 8, 53}, {}});
+ ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
+ ev.outgoing.push_back({gluon, { -48, -100, 82, 138}, {}});
+ ev.incoming[0] = {gluon, { 0, 0, -394, 394}, {}};
+ ev.incoming[1] = {gluon, { 0, 0, 304, 304}, {}};
return ev;
case -4: // jet idx: 0 1 2 3 4 5 5
- ev.outgoing.push_back({gluon, { -5, -40, -56, 69}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, -88, 133}, {}});
- ev.outgoing.push_back({gluon, { 23, -84, -72, 113}, {}});
- ev.outgoing.push_back({gluon, { -11, 92, 8, 93}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, 24, 113}, {}});
- ev.outgoing.push_back({gluon, { -15, -58, 30, 67}, {}});
- ev.outgoing.push_back({gluon, { -48, -96, 96, 144}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -395, 395}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 337, 337}, {}};
+ ev.outgoing.push_back({gluon, { 23, -52, -76, 95}, {}});
+ ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}});
+ ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}});
+ ev.outgoing.push_back({gluon, { -48, -44, -10, 66}, {}});
+ ev.outgoing.push_back({gluon, { -5, -92, -12, 93}, {}});
+ ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
+ ev.outgoing.push_back({gluon, { -11, 98, 70, 121}, {}});
+ ev.incoming[0] = {gluon, { 0, 0, -387, 387}, {}};
+ ev.incoming[1] = {gluon, { 0, 0, 333, 333}, {}};
return ev;
case -5: // jet idx: 0 1 -1 -1 2 3 4
- ev.outgoing.push_back({gluon, { 23, -64, -80, 105}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, -16, 101}, {}});
- ev.outgoing.push_back({gluon, { -15, 20, 0, 25}, {}});
- ev.outgoing.push_back({gluon, { -11, -10, 2, 15}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, 24, 113}, {}});
- ev.outgoing.push_back({gluon, { -48, -72, 54, 102}, {}});
- ev.outgoing.push_back({gluon, { -5, -60, 48, 77}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -253, 253}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 285, 285}, {}};
+ 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 0 0 1 2 2 3
- ev.outgoing.push_back({gluon, { -5, -60, -48, 77}, {}});
- ev.outgoing.push_back({gluon, { -15, -52, -36, 65}, {}});
- ev.outgoing.push_back({gluon, { -48, -96, -58, 122}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, -24, 113}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, -16, 101}, {}});
- ev.outgoing.push_back({gluon, { -11, 92, 8, 93}, {}});
- ev.outgoing.push_back({gluon, { 23, -70, 14, 75}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -403, 403}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 243, 243}, {}};
+ ev.outgoing.push_back({gluon, { -15, -62, -26, 69}, {}});
+ ev.outgoing.push_back({gluon, { -48, -88, -31, 105}, {}});
+ ev.outgoing.push_back({gluon, { -5, -38, -6, 39}, {}});
+ ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
+ ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
+ ev.outgoing.push_back({gluon, { -11, 98, 70, 121}, {}});
+ ev.outgoing.push_back({gluon, { 23, -98, 74, 125}, {}});
+ ev.incoming[0] = {gluon, { 0, 0, -265, 265}, {}};
+ ev.incoming[1] = {gluon, { 0, 0, 422, 422}, {}};
return ev;
case -7: // jet idx: 0 1 1 2 2 3 4
- ev.outgoing.push_back({gluon, { 23, -46, -46, 69}, {}});
- ev.outgoing.push_back({gluon, { -15, -90, -70, 115}, {}});
- ev.outgoing.push_back({gluon, { -5, -78, -54, 95}, {}});
- ev.outgoing.push_back({gluon, { -11, 88, -28, 93}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, -16, 101}, {}});
- ev.outgoing.push_back({gluon, { -48, -60, 5, 77}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, 24, 113}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -424, 424}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 239, 239}, {}};
+ ev.outgoing.push_back({gluon, { 23, -70, -66, 99}, {}});
+ ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}});
+ ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
+ ev.outgoing.push_back({gluon, { -5, -38, -14, 41}, {}});
+ ev.outgoing.push_back({gluon, { -48, -100, -37, 117}, {}});
+ ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}});
+ ev.outgoing.push_back({gluon, { -15, -78, 30, 85}, {}});
+ ev.incoming[0] = {gluon, { 0, 0, -462, 462}, {}};
+ ev.incoming[1] = {gluon, { 0, 0, 229, 229}, {}};
return ev;
case -8: // jet idx: 0 1 2 2 2 3 4
- ev.outgoing.push_back({gluon, { 23, -84, -84, 121}, {}});
- ev.outgoing.push_back({gluon, { -11, 92, -8, 93}, {}});
- ev.outgoing.push_back({gluon, { -15, -36, 0, 39}, {}});
- ev.outgoing.push_back({gluon, { -5, -62, 10, 63}, {}});
- ev.outgoing.push_back({gluon, { -48, -96, 19, 109}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, 24, 113}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, 88, 133}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -311, 311}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 360, 360}, {}};
+ ev.outgoing.push_back({gluon, { -48, -40, -56, 84}, {}});
+ ev.outgoing.push_back({gluon, { -11, 92, -68, 115}, {}});
+ ev.outgoing.push_back({gluon, { 23, -94, -62, 115}, {}});
+ ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}});
+ ev.outgoing.push_back({gluon, { -5, -48, -16, 51}, {}});
+ ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
+ ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}});
+ ev.incoming[0] = {gluon, { 0, 0, -449, 449}, {}};
+ ev.incoming[1] = {gluon, { 0, 0, 261, 261}, {}};
return ev;
case -9: // jet idx: 0 1 2 1 3 0 4
- ev.outgoing.push_back({gluon, { -12, 99, -44, 109}, {}});
- ev.outgoing.push_back({gluon, { -48, -80, -30, 98}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, -24, 113}, {}});
- ev.outgoing.push_back({gluon, { -15, -90, -18, 93}, {}});
- ev.outgoing.push_back({gluon, { 23, -70, -14, 75}, {}});
- ev.outgoing.push_back({gluon, { -11, 92, -8, 93}, {}});
- ev.outgoing.push_back({gluon, { -5, -38, 50, 63}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -366, 366}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 278, 278}, {}};
+ ev.outgoing.push_back({gluon, { -15, -96, -72, 121}, {}});
+ ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
+ ev.outgoing.push_back({gluon, { -48, -72, -39, 95}, {}});
+ ev.outgoing.push_back({gluon, { -11, 92, -40, 101}, {}});
+ ev.outgoing.push_back({gluon, { 23, -20, -12, 33}, {}});
+ ev.outgoing.push_back({gluon, { -5, -92, -36, 99}, {}});
+ ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}});
+ ev.incoming[0] = {gluon, { 0, 0, -476, 476}, {}};
+ ev.incoming[1] = {gluon, { 0, 0, 201, 201}, {}};
return ev;
case -10: // jet idx: 0 1 3 2 4 3 1
- ev.outgoing.push_back({gluon, { 23, -80, -64, 105}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, -44, 109}, {}});
- ev.outgoing.push_back({gluon, { -15, -60, -20, 65}, {}});
- ev.outgoing.push_back({gluon, { -48, -54, -16, 74}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, -24, 113}, {}});
- ev.outgoing.push_back({gluon, { -5, -84, -12, 85}, {}});
- ev.outgoing.push_back({gluon, { -11, 92, -8, 93}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -416, 416}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 228, 228}, {}};
+ 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, -203, 203}, {}};
+ ev.incoming[1] = {gluon, { 0, 0, 501, 501}, {}};
return ev;
case -11: // jet idx: 0 1 2 3 3 4 2
- ev.outgoing.push_back({gluon, { -48, -56, -48, 88}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, -24, 113}, {}});
- ev.outgoing.push_back({gluon, { -5, -62, -10, 63}, {}});
- ev.outgoing.push_back({gluon, { -11, 92, 8, 93}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, 16, 101}, {}});
- ev.outgoing.push_back({gluon, { 23, -70, 14, 75}, {}});
- ev.outgoing.push_back({gluon, { -15, -90, 18, 93}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -326, 326}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 300, 300}, {}};
+ 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 -12: // jet idx: 0 1 0 2 3 4 5
- ev.outgoing.push_back({gluon, { -15, -58, -30, 67}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, -24, 113}, {}});
- ev.outgoing.push_back({gluon, { -48, -96, -19, 109}, {}});
- ev.outgoing.push_back({gluon, { 23, -54, -6, 59}, {}});
- ev.outgoing.push_back({gluon, { -11, 92, -8, 93}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, 88, 133}, {}});
- ev.outgoing.push_back({gluon, { -5, -70, 86, 111}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -299, 299}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 386, 386}, {}};
+ ev.outgoing.push_back({gluon, { -5, -72, -92, 117}, {}});
+ ev.outgoing.push_back({gluon, { -48, -44, -62, 90}, {}});
+ ev.outgoing.push_back({gluon, { -15, -100, -80, 129}, {}});
+ ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
+ ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}});
+ ev.outgoing.push_back({gluon, { -11, 98, 70, 121}, {}});
+ ev.outgoing.push_back({gluon, { 23, -70, 66, 99}, {}});
+ ev.incoming[0] = {gluon, { 0, 0, -459, 459}, {}};
+ ev.incoming[1] = {gluon, { 0, 0, 325, 325}, {}};
return ev;
- case -13: // jet idx: 0 1 2 2 2 3 0
- ev.outgoing.push_back({gluon, { -11, 88, -28, 93}, {}});
- ev.outgoing.push_back({gluon, { 68, 87, -24, 113}, {}});
- ev.outgoing.push_back({gluon, { -48, -84, -21, 99}, {}});
- ev.outgoing.push_back({gluon, { -15, -90, -18, 93}, {}});
- ev.outgoing.push_back({gluon, { -5, -30, -6, 31}, {}});
- ev.outgoing.push_back({gluon, { 23, -70, -14, 75}, {}});
- ev.outgoing.push_back({gluon, { -12, 99, -16, 101}, {}});
- ev.incoming[0] = {gluon, { 0, 0, -366, 366}, {}};
- ev.incoming[1] = {gluon, { 0, 0, 239, 239}, {}};
+ case -13: // jet idx: 0 1 2 0 1 3 0
+ ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}});
+ ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}});
+ ev.outgoing.push_back({gluon, { 23, -40, -16, 49}, {}});
+ ev.outgoing.push_back({gluon, { -48, -100, -37, 117}, {}});
+ ev.outgoing.push_back({gluon, { -11, 88, -20, 91}, {}});
+ ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}});
+ ev.outgoing.push_back({gluon, { -5, -38, -6, 39}, {}});
+ ev.incoming[0] = {gluon, { 0, 0, -429, 429}, {}};
+ ev.incoming[1] = {gluon, { 0, 0, 212, 212}, {}};
return ev;
}
}
throw HEJ::unknown_option{"unkown process"};
}
bool couple_quark(std::string const & boson, std::string & quark){
if(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 && !(abs(qflav)%2)) return false; // not anti-down
if(W_charge*qflav > 0 && (abs(qflav)%2)) return false; // not up
quark=std::to_string(qflav-W_charge);
}
return true;
}
+ 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){
+ decays[0] = HEJ::ParticleID::e_bar;
+ decays[1] = HEJ::ParticleID::nu_e;
+ } else {
+ decays[0] = HEJ::ParticleID::e;
+ decays[1] = HEJ::ParticleID::nu_e_bar;
+ }
+ std::vector<HEJ::Particle> decay_products(decays.size());
+ for(size_t i = 0; i < decays.size(); ++i){
+ decay_products[i].type = decays[i];
+ }
+ // choose polar and azimuth angle in parent rest frame
+ const double E = parent.m()/2;
+ const double theta = 0.;
+ const double cos_phi = 1.;
+ const double sin_phi = sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
+ const double px = E*cos(theta)*sin_phi;
+ const double py = E*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;
+ }
+
// create event corresponding from given Configuration
// overwrite_boson to force a specific boson position, indepentent from input
// (useful for njet == 7)
HEJ::Event parse_configuration(
std::array<std::string,2> const & in, std::vector<std::string> const & out,
int const overwrite_boson = 0
){
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));
size_t njets = out.size();
- if(boson != out.cend()) --njets;
+ 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(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]);
}
for(size_t i=0; i<in.size(); ++i){
ev.incoming[i].type = HEJ::to_ParticleID(in[i]);
}
shuffle_particles(ev);
return ev.cluster(jet_def, min_jet_pt);
}
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)};
if(ev.type() != expected){
std::cerr << "Expected type " << HEJ::event_type::name(expected)
<< " but found " << HEJ::event_type::name(ev.type()) << "\n" << ev;
auto jet_idx{ ev.particle_jet_indices(ev.jets()) };
std::cout << "Particle Jet indices: ";
for(int const i: jet_idx)
std::cout << i << " ";
std::cout << std::endl;
return false;
}
return true;
}
- bool check_fkl(){
+ //! 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(!match_expectation(event_type::FKL, base_in, base_out))
+ if(implemented && !match_expectation(type, base_in, base_out))
return false;
- for(auto const & boson: all_bosons) // any boson
+ 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(event_type::FKL, in, out))
+ if(!match_expectation(type, in, out))
return false;
}
}
- for(int i=-3; i>-13;--i){ // extremal parton inside jet
- std::array<std::string,2> base_in{first,last};
- std::vector<std::string> base_out(7, "g");
- base_out.front() = first;
- base_out.back() = last;
- if(!match_expectation(event_type::FKL, base_in, base_out, i))
+ if(implemented){
+ for(int i=-3; i>-13;--i){ // extremal parton inside jet
+ std::array<std::string,2> base_in{first,last};
+ std::vector<std::string> base_out(7, "g");
+ base_out.front() = first;
+ base_out.back() = last;
+ if(!match_expectation(type, base_in, base_out, i))
+ return false;
+ }
+ if(! match_expectation(type,{first,last},{"h",first,"g","g","g","g",last}, -13)
+ && match_expectation(type,{first,last},{first,"g","g","g","g",last,"h"}, -13)
+ && match_expectation(type,{"2",last}, {"1","g","g","g","g",last,"Wp"}, -13)
+ && match_expectation(type,{first,"1"}, {"Wp",first,"g","g","g","g","2"}, -13)
+ )
return false;
}
}
return true;
}
- bool check_uno(){
+ //! 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{ event_type::unob };
- auto const f{ event_type::unof };
+ 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( !match_expectation(expectation, base_in, base_out) )
+ if( implemented
+ && !match_expectation(expectation, base_in, base_out) )
return false;
- for(auto const & boson: all_bosons){ // any boson
+ 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;
}
}
}
}
// test allowed jet configurations
- if(!(match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -3)
- && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -4)
- && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -5)
- && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -5)
- && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -6)
- && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -7)
- && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -7)
- && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -8)
- && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -8)
- && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -9)
- && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -10)
- && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -11)
- && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -11)
- && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -12)
- && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -12)
+ if( implemented && !(
+ match_expectation(f,{fkl,uno},{"h",fkl,"g","g","g",uno,"g"}, -1)
+ && match_expectation(b,{uno,"2"},{"Wp","g",uno,"g","g","g","1"}, -1)
+ && match_expectation(b,{"1",fkl},{"Wm","g","2","g","g","g",fkl}, -1)
+ && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -3)
+ && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -4)
+ && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -5)
+ && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -5)
+ && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -6)
+ && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -7)
+ && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -7)
+ && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -8)
+ && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -8)
+ && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -9)
+ && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -10)
+ && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -11)
+ && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -11)
+ && match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -12)
+ && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -12)
+ && match_expectation(f,{fkl,uno},{"h",fkl,"g","g","g",uno,"g"}, -13)
+ && match_expectation(b,{uno,fkl},{"g",uno,"g","g","g",fkl,"h"}, -13)
+ && match_expectation(f,{"2",uno},{"1","g","g","g",uno,"g","Wp"}, -13)
+ && match_expectation(f,{fkl,"1"},{"Wm",fkl,"g","g","g","2","g"}, -13)
+ && match_expectation(b,{uno,"1"},{"Wm","g",uno,"g","g","g","2"}, -13)
))
return false;
}
return true;
}
- bool check_extremal_qqx(){
+ //! 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{ event_type::qqxexb };
- auto const f{ event_type::qqxexf };
+ 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( !match_expectation(expectation, base_in, base_out) )
+ if( !implemented
+ && !match_expectation(expectation, base_in, base_out) )
return false;
- for(auto const & boson: all_bosons){ // any boson
+ 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(!(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( !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;
+ } else 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;
}
- bool check_central_qqx(){
+ //! 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{ event_type::qqxmid };
+ 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( !match_expectation(t, base_in, base_out) )
+ if( !implemented && !match_expectation(t, base_in, base_out) )
return false;
- for(auto const & boson: all_bosons) // any boson
+ 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;
}
}
}
// test allowed jet configurations (non exhaustive collection)
- if(!(match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -3)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -3)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -4)
- && match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -4)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -5)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -6)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -7)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -7)
- && match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -8)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -8)
- && match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -9)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -9)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -10)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -10)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -11)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -12)
- && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -12)
- ))
- return false;
+ if( !implemented ){
+ if( !( match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -3)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -3)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -4)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -4)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -5)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -6)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -7)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -7)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -8)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -8)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -9)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -9)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -10)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -10)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -11)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -12)
+ && match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -12)
+ ))
+ return false;
+ } else if (fkl1 == "1") {
+ if( !( match_expectation(t,{"1",fkl2},{"2","Wm",qqx,qqx2,"g","g",fkl2}, -3)
+ && match_expectation(t,{"1",fkl2},{"2","g","Wm",qqx,qqx2,"g",fkl2}, -3)
+ && match_expectation(t,{"1",fkl2},{"2","g","g",qqx,qqx2,"Wm",fkl2}, -4)
+ && match_expectation(t,{"1",fkl2},{"2",qqx,qqx2,"g","Wm","g",fkl2}, -4)
+ && match_expectation(t,{"1",fkl2},{"2","g","g","Wm",qqx,qqx2,fkl2}, -5)
+ && match_expectation(t,{"1",fkl2},{"2","g","Wm",qqx,qqx2,"g",fkl2}, -6)
+ && match_expectation(t,{"1",fkl2},{"2","Wm",qqx,qqx2,"g","g",fkl2}, -7)
+ && match_expectation(t,{"1",fkl2},{"2","g","Wm","g",qqx,qqx2,fkl2}, -7)
+ && match_expectation(t,{"1",fkl2},{"2",qqx,qqx2,"Wm","g","g",fkl2}, -8)
+ && match_expectation(t,{"1",fkl2},{"2","g","g","Wm",qqx,qqx2,fkl2}, -8)
+ && match_expectation(t,{"1",fkl2},{"2",qqx,qqx2,"g","Wm","g",fkl2}, -9)
+ && match_expectation(t,{"1",fkl2},{"2","g",qqx,qqx2,"g","Wm",fkl2}, -9)
+ && match_expectation(t,{"1",fkl2},{"2","g","g",qqx,qqx2,"Wm",fkl2}, -10)
+ && match_expectation(t,{"1",fkl2},{"2","g","g","Wm",qqx,qqx2,fkl2}, -10)
+ && match_expectation(t,{"1",fkl2},{"2","g","Wm","g",qqx,qqx2,fkl2}, -11)
+ && match_expectation(t,{"1",fkl2},{"2","Wm","g",qqx,qqx2,"g",fkl2}, -12)
+ && match_expectation(t,{"1",fkl2},{"2","g","Wm","g",qqx,qqx2,fkl2}, -12)
+ ))
+ 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
;
}
// Events failing the jet requirements, e.g. qqx inside one jet
- //! TODO qqx currently not allow for pure jet -> need W for qqx
bool check_illegal_jets(){
auto type{ HEJ::event_type::non_resummable};
return true
// uno backward not in jet
&& match_expectation(type, {"1","g"}, {"g","1","g","g","g","g","g"}, -1)
// & also legal uno on other side
&& match_expectation(type, {"1","1"}, {"g","1","g","g","g","1","g"}, -1)
// qqx backward not in jet
&& match_expectation(type, {"g","2"}, {"-1","1","g","g","g","g","2"}, -1)
// uno forward not in jet
&& match_expectation(type, {"3","3"}, {"3","g","g","g","g","3","g"}, -2)
// qqx backward not in jet
&& match_expectation(type, {"g","g"}, {"g","g","g","g","g","-2","2"}, -2)
// uno backward in same jet
&& match_expectation(type, {"1","g"}, {"g","1","g","g","g","g","g"}, -3)
// & also legal uno on other side
&& match_expectation(type, {"1","1"}, {"g","1","g","g","g","1","g"}, -3)
// qqx backward in same jet
&& match_expectation(type, {"g","2"}, {"-4","4","g","g","g","g","2"}, -3)
// uno forward in same jet
&& match_expectation(type, {"3","2"}, {"3","g","g","g","g","2","g"}, -4)
// qqx backward in same jet
&& match_expectation(type, {"g","g"}, {"g","g","g","g","g","-2","2"}, -4)
// central qqx not in jet
&& match_expectation(type, {"1","2"}, {"1","g","-1","1","g","g","2"}, -5)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","-1","1","g","g","g","2"}, -6)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","g","g","g","2","-2","2"}, -6)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","-1","1","g","g","g","2"}, -7)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","g","g","3","-3","g","2"}, -7)
// central qqx in same jet
&& match_expectation(type, {"g","3"}, {"g","g","-2","2","g","g","3"}, -8)
// central qqx in same jet
&& match_expectation(type, {"g","-2"}, {"g","g","g","2","-2","g","-2"}, -8)
// FKL not in own jet
&& match_expectation(type, {"1","1"}, {"1","g","g","g","g","g","1"}, -1)
&& match_expectation(type, {"g","g"}, {"g","g","g","g","g","g","g"}, -2)
&& match_expectation(type, {"1","-3"}, {"1","g","g","g","g","-3","g"}, -1)
&& match_expectation(type, {"2","g"}, {"g","2","g","g","g","g","g"}, -2)
&& match_expectation(type, {"2","g"}, {"2","g","g","g","g","2","-2"}, -1)
&& match_expectation(type, {"g","-2"}, {"4","-4","g","g","g","g","-2"}, -2)
&& match_expectation(type, {"-1","g"}, {"-1","1","-1","g","g","g","g"}, -1)
&& match_expectation(type, {"4","-3"}, {"4","g","g","1","-1","g","-3"}, -2)
// FKL in same jet
&& match_expectation(type, {"1","1"}, {"1","g","g","g","g","g","1"}, -13)
// uno in same jet as FKL
&& match_expectation(type, {"-1","1"}, {"-1","g","g","g","g","1","g"}, -9)
&& match_expectation(type, {"-1","1"}, {"g","-1","g","g","g","g","1"}, -10)
&& match_expectation(type, {"-1","1"}, {"-1","g","g","g","g","1","g"}, -13)
&& match_expectation(type, {"-1","1"}, {"g","-1","g","g","g","g","1"}, -13)
// extremal qqx in same jet as FKL
&& match_expectation(type, {"-1","g"}, {"-1","g","g","g","g","1","-1"}, -9)
&& match_expectation(type, {"g","1"}, {"1","-1","g","g","g","g","1"}, -10)
&& match_expectation(type, {"-1","g"}, {"-1","g","g","g","g","-1","1"}, -13)
&& match_expectation(type, {"g","g"}, {"g","g","g","g","g","-1","1"}, -13)
// central qqx in same jet as FKL
&& match_expectation(type,{"2","2"}, {"2","-2","2","g","g","g","2"}, -3)
&& match_expectation(type,{"-1","g"}, {"-1","g","g","g","2","-2","g"}, -4)
&& match_expectation(type,{"g","4"}, {"g","g","-4","4","g","g","4"}, -6)
&& match_expectation(type,{"g","g"}, {"g","3","-3","g","g","g","g"}, -6)
&& match_expectation(type,{"g","1"}, {"g","g","g","g","-2","2","1"}, -9)
&& match_expectation(type,{"g","g"}, {"g","1","-1","g","g","g","g"}, -10)
&& match_expectation(type,{"g","1"}, {"g","g","-2","2","g","g","1"}, -11)
&& match_expectation(type,{"g","1"}, {"g","-2","2","g","g","g","1"}, -11)
&& match_expectation(type,{"3","2"}, {"3","g","-1","1","g","g","2"}, -12)
&& match_expectation(type,{"3","-2"}, {"3","-1","1","g","g","g","-2"}, -12)
&& match_expectation(type,{"3","-2"}, {"3","-1","1","g","g","g","-2"}, -13)
;
}
// 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_illegal_jets()) 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 a74f634..96aab58 100644
--- a/t/test_classify_ref.cc
+++ b/t/test_classify_ref.cc
@@ -1,100 +1,105 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <fstream>
#include <random>
#include <algorithm>
#include "LHEF/LHEF.h"
#include "HEJ/YAMLreader.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};
void shuffle_particles(HEJ::Event::EventData & ev) {
static std::mt19937_64 ran{0};
// incoming
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
// outgoing (through index)
auto old_outgoing = std::move(ev.outgoing);
std::vector<size_t> idx(old_outgoing.size());
std::iota(idx.begin(), idx.end(), 0);
std::shuffle(begin(idx), end(idx), ran);
ev.outgoing.clear();
ev.outgoing.reserve(old_outgoing.size());
for(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(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(begin(decay.second), end(decay.second), ran);
+ }
}
}
}
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]) };
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(ev.jets()))
std::cerr << idx << " ";
std::cerr << "\n";
return EXIT_FAILURE;
}
}
}
return EXIT_SUCCESS;
}
diff --git a/t/test_colours.cc b/t/test_colours.cc
index a7b96c4..8aa4031 100644
--- a/t/test_colours.cc
+++ b/t/test_colours.cc
@@ -1,368 +1,405 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <random>
#include <stdexcept>
#include <utility>
#include "HEJ/Constants.hh"
#include "HEJ/Event.hh"
#include "HEJ/RNG.hh"
#define ASSERT(x) if(!(x)) { \
throw std::logic_error("Assertion '" #x "' failed."); \
}
/// biased RNG to connect always to colour
class dum_rnd: public HEJ::DefaultRNG {
public:
dum_rnd() = default;
double flat() override {
return 0.;
};
};
void shuffle_particles(HEJ::Event::EventData & ev) {
static std::mt19937_64 ran{0};
// incoming
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
// outgoing (through index)
auto old_outgoing = std::move(ev.outgoing);
std::vector<size_t> idx(old_outgoing.size());
std::iota(idx.begin(), idx.end(), 0);
std::shuffle(begin(idx), end(idx), ran);
ev.outgoing.clear();
ev.outgoing.reserve(old_outgoing.size());
for(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(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(begin(decay.second), end(decay.second), ran);
+ }
+ }
+}
+
+std::vector<HEJ::Particle> decay_Wp( HEJ::Particle const & parent ){
+ const std::array<HEJ::ParticleID, 2> decays{ HEJ::ParticleID::e_bar,
+ HEJ::ParticleID::nu_e };
+ std::vector<HEJ::Particle> decay_products(decays.size());
+ for(size_t i = 0; i < decays.size(); ++i){
+ decay_products[i].type = decays[i];
+ }
+ // choose polar and azimuth angle in parent rest frame
+ const double E = parent.m()/2;
+ const double theta = 0.;
+ const double cos_phi = 1.;
+ const double sin_phi = sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
+ const double px = E*cos(theta)*sin_phi;
+ const double py = E*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;
+}
+
+HEJ::Event::EventData decay_boson( HEJ::Event::EventData ev ){
+ for( size_t i=0; i<ev.outgoing.size(); ++i ){
+ if( ev.outgoing[i].type == HEJ::ParticleID::Wp){
+ ev.decays[i] = decay_Wp(ev.outgoing[i]);
+ }
}
+ return ev;
}
+
void dump_event(HEJ::Event const & ev){
for(auto const & in: ev.incoming()){
std::cerr << "in type=" << in.type
<< ", colour={" << (*in.colour).first
<< ", " << (*in.colour).second << "}\n";
}
for(auto const & out: ev.outgoing()){
std::cerr << "out type=" << out.type << ", colour={";
if(out.colour)
std::cerr << (*out.colour).first << ", " << (*out.colour).second;
else
std::cerr << "non, non";
std::cerr << "}\n";
}
}
/// true if colour is allowed for particle
bool correct_colour(HEJ::Particle const & part){
if(!HEJ::is_parton(part) && !part.colour) return true;
if(!part.colour) return false;
int const colour = part.colour->first;
int const anti_colour = part.colour->second;
if(part.type == HEJ::ParticleID::gluon)
return colour != anti_colour
&& colour >= HEJ::COLOUR_OFFSET
&& anti_colour >= HEJ::COLOUR_OFFSET;
if(HEJ::is_quark(part))
return anti_colour == 0 && colour >= HEJ::COLOUR_OFFSET;
return colour == 0 && anti_colour >= HEJ::COLOUR_OFFSET;
}
bool correct_colour(HEJ::Event const & ev){
if(!ev.is_leading_colour())
return false;
// some of these additional checks are also in ev.is_leading_colour()
for(auto const & part: ev.incoming()){
if(!correct_colour(part))
return false;
}
for(auto const & part: ev.outgoing()){
if(!correct_colour(part))
return false;
}
return true;
}
bool match_expected(
HEJ::Event const & ev,
std::vector<HEJ::Colour> const & expected
){
ASSERT(ev.outgoing().size()+2==expected.size());
for(size_t i=0; i<ev.incoming().size(); ++i){
ASSERT(ev.incoming()[i].colour);
if( *ev.incoming()[i].colour != expected[i])
return false;
}
for(size_t i=2; i<ev.outgoing().size()+2; ++i){
if( ev.outgoing()[i-2].colour ){
if( *ev.outgoing()[i-2].colour != expected[i] )
return false;
} else if( expected[i].first != 0 || expected[i].second != 0)
return false;
}
return true;
}
void check_event(
HEJ::Event::EventData unc_ev,
std::vector<HEJ::Colour> const & expected_colours
){
+ unc_ev = decay_boson(std::move(unc_ev));
shuffle_particles(unc_ev); // make sure incoming order doesn't matter
HEJ::Event ev{unc_ev.cluster(
fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30.)
};
ASSERT(HEJ::event_type::is_resummable(ev.type()));
dum_rnd rng;
ASSERT(!ev.is_leading_colour());
ASSERT(ev.generate_colours(rng));
if(!correct_colour(ev)){
std::cerr << "Found illegal colours for event\n";
dump_event(ev);
throw std::invalid_argument("Illegal colour set");
}
if(!match_expected(ev, expected_colours)){
std::cerr << "Colours didn't match expectation. Found\n";
dump_event(ev);
std::cerr << "but expected\n";
for(auto const & col: expected_colours){
std::cerr << "colour={" << col.first << ", " << col.second << "}\n";
}
throw std::logic_error("Colours did not match expectation");
}
}
HEJ::Event::EventData reset_colour(
HEJ::Event::EventData ev, std::vector<HEJ::Colour> const & goal
){
for(size_t i=0; i<2; ++i){
ev.incoming[i].colour = goal[i];
}
for(size_t i=0; i<ev.outgoing.size(); ++i){
auto const & col_goal{ goal[i+2] };
if(col_goal.first == 0 && col_goal.second == 0)
ev.outgoing[i].colour = HEJ::optional<HEJ::Colour>{};
else
ev.outgoing[i].colour = col_goal;
}
return ev;
}
int main() {
HEJ::Event::EventData ev;
std::vector<HEJ::Colour> expected_colours(7);
- /// pure gluon
- ev.incoming[0] = { HEJ::ParticleID::gluon, { 0, 0,-427, 427}, {}};
- ev.incoming[1] = { HEJ::ParticleID::gluon, { 0, 0, 851, 851}, {}};
- ev.outgoing.push_back({ HEJ::ParticleID::gluon, { 196, 124, -82, 246}, {}});
- ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-167,-184, 16, 249}, {}});
- ev.outgoing.push_back({ HEJ::ParticleID::higgs, { 197, 180, 168, 339}, {}});
- ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-190, -57, 126, 235}, {}});
- ev.outgoing.push_back({ HEJ::ParticleID::gluon, { -36, -63, 196, 209}, {}});
+ /// pure gluon (they all have a mass of 4 GeV to allow decays)
+ ev.incoming[0] = { HEJ::ParticleID::gluon, { 0, 0, -205, 205}, {}};
+ ev.incoming[1] = { HEJ::ParticleID::gluon, { 0, 0, 279, 279}, {}};
+ ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-15, -82, -82, 117}, {}});
+ ev.outgoing.push_back({ HEJ::ParticleID::gluon, { 68, 93, 20, 117}, {}});
+ ev.outgoing.push_back({ HEJ::ParticleID::higgs, {-30, -65, 22, 75}, {}});
+ ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-12, 92, 76, 120}, {}});
+ ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-11, -38, 38, 55}, {}});
expected_colours[0] = {502, 501};
expected_colours[1] = {509, 502};
expected_colours[2] = {503, 501};
expected_colours[3] = {505, 503};
expected_colours[4] = { 0, 0};
expected_colours[5] = {507, 505};
expected_colours[6] = {509, 507};
// default colours is always forbidden!
// default: swap last two (anti-)colour -> crossing
ev=reset_colour(ev, expected_colours);
std::swap(ev.outgoing[4].colour, ev.outgoing[3].colour);
check_event(ev, expected_colours);
/// last g to Qx (=> gQx -> g ... Qx)
ev.incoming[1].type = HEJ::ParticleID::d_bar;
ev.outgoing[4].type = HEJ::ParticleID::d_bar;
// => only end changes
expected_colours[1].first = 0;
expected_colours[6].first = 0;
// default: swap last two anti-colours -> last gluon colour singlet
ev=reset_colour(ev, expected_colours);
std::swap(ev.outgoing[4].colour->second, ev.outgoing[3].colour->second);
check_event(ev, expected_colours);
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno forward (=> gQx -> g ... Qx g)
std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
// => uno quarks eats colour and gluon connects to anti-colour
new_expected[5] = {0, expected_colours[3].first};
new_expected[6] = {expected_colours[0].first, expected_colours[0].first+2};
new_expected[1].second += 2; // one more anti-colour in line
// default: swap last two anti-colours -> crossing
new_ev=reset_colour(new_ev, new_expected);
std::swap(new_ev.outgoing[4].colour->second, new_ev.outgoing[3].colour->second);
check_event(new_ev, new_expected);
}
/// swap Qx <-> Q (=> gQ -> g ... Q)
ev.incoming[1].type = HEJ::ParticleID::d;
ev.outgoing[4].type = HEJ::ParticleID::d;
// => swap: colour<->anti && initial<->final
std::swap(expected_colours[1], expected_colours[6]);
std::swap(expected_colours[1].first, expected_colours[1].second);
std::swap(expected_colours[6].first, expected_colours[6].second);
// default: swap incoming <-> outgoing
ev=reset_colour(ev, expected_colours);
std::swap(ev.incoming[0].colour, ev.outgoing[0].colour);
check_event(ev, expected_colours);
/// first g to qx (=> qxQ -> qx ... Q)
ev.incoming[0].type = HEJ::ParticleID::u_bar;
ev.outgoing[0].type = HEJ::ParticleID::u_bar;
expected_colours[0] = { 0, 501};
// => shift anti-colour index one up
expected_colours[1].first -= 2;
expected_colours[5] = expected_colours[3];
expected_colours[3] = expected_colours[2];
expected_colours[2] = { 0, 502};
// default: closed qx->qx g
ev=reset_colour(ev, expected_colours);
ev.outgoing[1].colour->first = ev.outgoing[0].colour->second;
ev.outgoing[1].colour->second = ev.incoming[0].colour->second;
ev.outgoing[4].colour->first = ev.outgoing[3].colour->second;
check_event(ev, expected_colours);
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno backward (=> qxQ -> g qx ... Q)
std::swap(new_ev.outgoing[0].type, new_ev.outgoing[1].type);
// => uno gluon connects to quark colour
new_expected[3] = expected_colours[2];
new_expected[2] = {expected_colours[0].second+2, expected_colours[0].second};
// default: Colourful Higgs
new_ev=reset_colour(new_ev, new_expected);
new_ev.outgoing[2].colour = std::make_pair(1,1);
check_event(new_ev, new_expected);
/// swap qx <-> q (=> qQ -> g q ... Q)
new_ev.incoming[0].type = HEJ::ParticleID::u;
new_ev.outgoing[1].type = HEJ::ParticleID::u;
// => swap: colour<->anti && inital<->final
std::swap(new_expected[0], new_expected[3]);
std::swap(new_expected[0].first, new_expected[0].second);
std::swap(new_expected[3].first, new_expected[3].second);
// => & connect first gluon with remaining anti-colour
new_expected[2] = {new_expected[0].first, new_expected[0].first+2};
// shift colour line one down
new_expected[1].first-=2;
new_expected[5].first-=2;
new_expected[5].second-=2;
// shift anti-colour line one up
new_expected[6].first+=2;
// default: swap 2 quarks -> disconnected
new_ev=reset_colour(new_ev, new_expected);
std::swap(new_ev.outgoing[1].colour, new_ev.outgoing[4].colour);
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno forward (=> qxQ -> qx ... Q g)
std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
// => uno gluon connects to remaining colour
new_expected[5] = expected_colours[6];
new_expected[6] = {expected_colours[3].first+2, expected_colours[3].first};
// default: no colour on last gluon
new_ev=reset_colour(new_ev, new_expected);
new_ev.incoming[1].colour->first = new_ev.outgoing[4].colour->second;
new_ev.outgoing[4].colour = {};
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
- /// qqx backward (=> gQ -> qx q ... Q)
+ /// qqx backward (=> gQ -> qx q ... Q) with Wp
// => swap: incoming q <-> outgoing gluon
std::swap(new_ev.incoming[0].type, new_ev.outgoing[1].type);
new_ev.outgoing[1].type=static_cast<HEJ::ParticleID>(
- -1*new_ev.outgoing[1].type);
+ -(new_ev.outgoing[1].type+1) );
+ new_ev.outgoing[2].type = HEJ::ParticleID::Wp;
// incoming q -> outgoing q (colour<->anti)
std::swap(new_expected[0], new_expected[3]);
std::swap(new_expected[3].first, new_expected[3].second);
new_expected[3].first+=2;
new_expected[0].first-=1; // skip one index
// couple first in to first out
new_expected[2].second=new_expected[0].second;
// default: swap qqx <-> first g
new_ev=reset_colour(new_ev, new_expected);
std::swap(new_ev.outgoing[0].colour->second, new_ev.outgoing[3].colour->second);
std::swap(new_ev.outgoing[1].colour->first, new_ev.outgoing[3].colour->first);
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
- /// qqx forward (=> qx g -> qx ... Qx Q)
+ /// qqx forward (=> qx g -> qx ... Qx Q) with Wp
// => swap: incoming Q <-> outgoing gluon
std::swap(new_ev.incoming[1].type, new_ev.outgoing[3].type);
new_ev.outgoing[3].type=static_cast<HEJ::ParticleID>(
- -1*new_ev.outgoing[3].type);
+ -(new_ev.outgoing[3].type+1));
+ new_ev.outgoing[2].type = HEJ::ParticleID::Wp;
// incoming q -> outgoing q (colour<->anti)
std::swap(new_expected[1], new_expected[5]);
std::swap(new_expected[5].first, new_expected[5].second);
new_expected[5].second-=2;
new_expected[1].second-=1; // skip one index
// couple last in to last out
new_expected[6].first=new_expected[1].first;
// default: uncoloured quark
new_ev=reset_colour(new_ev, new_expected);
new_ev.outgoing[0].colour = {};
check_event(new_ev, new_expected);
// move Higgs to position 1 (=> qx g -> qx h g Qx Q)
std::swap(new_ev.outgoing[1].type, new_ev.outgoing[2].type);
std::swap(new_expected[3], new_expected[4]); // trivial
// default: incoming qx wrong colour
new_ev=reset_colour(new_ev, new_expected);
new_ev.incoming[0].colour->first = 1;
check_event(new_ev, new_expected);
// central qqx (=> qx g -> qx h Q Qx g)
// => swap: Q <-> g
std::swap(new_ev.outgoing[2].type, new_ev.outgoing[4].type);
std::swap(new_expected[4], new_expected[6]);
// gluon was connected on left side, i.e. doesn't matter for QQx
// => couple Q to out qx
new_expected[4].first = new_expected[2].second;
// Qx next in line
new_expected[5].second = new_expected[4].first+2;
// incoming g shifted by one position in line
new_expected[1].first-=2;
new_expected[1].second+=2;
// default: wrong colour in last incoming
new_ev=reset_colour(new_ev, new_expected);
std::swap(new_ev.incoming[1].colour->first,
new_ev.incoming[1].colour->second);
check_event(new_ev, new_expected);
}
return EXIT_SUCCESS;
}
diff --git a/t/test_decay.cc b/t/test_decay.cc
new file mode 100644
index 0000000..a7be172
--- /dev/null
+++ b/t/test_decay.cc
@@ -0,0 +1,193 @@
+/**
+ * \authors The HEJ collaboration (see AUTHORS for details)
+ * \date 2019
+ * \copyright GPLv2 or later
+ *
+ * \brief Test classification for (invalid) W decays
+ */
+#include <random>
+
+#include "HEJ/Event.hh"
+
+namespace {
+ const fastjet::JetDefinition jet_def{fastjet::JetAlgorithm::antikt_algorithm, 0.4};
+ const double min_jet_pt{30.};
+
+ void shuffle_particles(HEJ::Event::EventData & ev) {
+ static std::mt19937_64 ran{0};
+ // incoming
+ std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
+ // outgoing (through index)
+ auto old_outgoing = std::move(ev.outgoing);
+ std::vector<size_t> idx(old_outgoing.size());
+ std::iota(idx.begin(), idx.end(), 0);
+ std::shuffle(begin(idx), end(idx), ran);
+ ev.outgoing.clear();
+ ev.outgoing.reserve(old_outgoing.size());
+ for(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(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(begin(decay.second), end(decay.second), ran);
+ }
+ }
+ }
+
+ HEJ::Event::EventData new_event() {
+ HEJ::Event::EventData ev;
+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, -96, -76, 123}, {}});
+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -70, -22, 75}, {}});
+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 93, -20, 117}, {}});
+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 95, 56, 111}, {}});
+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { -30, -22, 25, 45}, {}});
+ ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -254, 254}, {}};
+ ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 217, 217}, {}};
+ return ev;
+ }
+
+ 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){
+ decays[0] = HEJ::ParticleID::nu_e;
+ decays[1] = HEJ::ParticleID::e_bar;
+ } else {
+ decays[0] = HEJ::ParticleID::e;
+ decays[1] = HEJ::ParticleID::nu_e_bar;
+ }
+ std::vector<HEJ::Particle> decay_products(decays.size());
+ for(size_t i = 0; i < decays.size(); ++i){
+ decay_products[i].type = decays[i];
+ }
+ // choose polar and azimuth angle in parent rest frame
+ const double E = parent.m()/2;
+ const double theta = 0.;
+ const double cos_phi = 1.;
+ const double sin_phi = sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
+ const double px = E*cos(theta)*sin_phi;
+ const double py = E*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;
+ }
+
+ bool test_event(HEJ::Event::EventData data, bool const valid
+ ){
+ using namespace HEJ::event_type;
+ EventType const expected{ valid?FKL:bad_final_state };
+ shuffle_particles(data);
+ auto const ev = std::move(data).cluster(jet_def, min_jet_pt);
+ if(ev.type() != expected){
+ std::cerr << "Event does not match expectation, expected "
+ << name(expected) << "\n" << ev << std::endl;
+ return false;
+ }
+ return true;
+ }
+} // namespace anonymous
+
+int main() {
+ using namespace HEJ::pid;
+ auto ev = new_event();
+ // basic FKL
+ test_event(ev, true);
+ // W position shouldn't matter
+ for(auto const W_type: {Wp, Wm}){
+ for(size_t w_pos = 1; w_pos<ev.outgoing.size()-1; ++w_pos){
+ ev = new_event();
+ ev.outgoing[w_pos].type = W_type;
+ ev.outgoing.back().type = (W_type==Wp)?d:u;
+ ev.incoming.back().type = (W_type==Wp)?u:d;
+
+ // no decay
+ if(!test_event(ev, false))
+ return EXIT_FAILURE;
+
+ // working decay
+ ev.decays[w_pos] = decay_W(ev.outgoing[w_pos]);
+ if(!test_event(ev, true))
+ return EXIT_FAILURE;
+
+ // swapped W+ <-> W-
+ ev.decays[w_pos].at(0).type = static_cast<ParticleID>(
+ -ev.decays[w_pos].at(0).type );
+ ev.decays[w_pos].at(1).type = static_cast<ParticleID>(
+ -ev.decays[w_pos].at(1).type );
+ if(!test_event(ev, false))
+ return EXIT_FAILURE;
+ ev.decays[w_pos].at(0).type = static_cast<ParticleID>(
+ -ev.decays[w_pos].at(0).type );
+ ev.decays[w_pos].at(1).type = static_cast<ParticleID>(
+ -ev.decays[w_pos].at(1).type );
+
+ // replace e -> mu (normal)
+ ev.decays[w_pos].at(0).type = static_cast<ParticleID>(
+ ev.decays[w_pos].at(0).type+2 );
+ if(!test_event(ev, false))
+ return EXIT_FAILURE;
+ ev.decays[w_pos].at(0).type = static_cast<ParticleID>(
+ ev.decays[w_pos].at(0).type-2 );
+
+ // replace e -> mu (anti)
+ ev.decays[w_pos].at(1).type = static_cast<ParticleID>(
+ ev.decays[w_pos].at(1).type-2 );
+ if(!test_event(ev, false))
+ return EXIT_FAILURE;
+
+ // all mu
+ ev.decays[w_pos].at(0).type = static_cast<ParticleID>(
+ ev.decays[w_pos].at(0).type+2 );
+ if(!test_event(ev, true))
+ return EXIT_FAILURE;
+ ev.decays[w_pos].at(0).type = static_cast<ParticleID>(
+ ev.decays[w_pos].at(0).type-2 );
+ ev.decays[w_pos].at(1).type = static_cast<ParticleID>(
+ ev.decays[w_pos].at(1).type+2 );
+
+ // partonic
+ ev.decays[w_pos].at(0).type = static_cast<ParticleID>(
+ ev.decays[w_pos].at(0).type-10 );
+ ev.decays[w_pos].at(1).type = static_cast<ParticleID>(
+ ev.decays[w_pos].at(1).type+10 );
+ if(!test_event(ev, false))
+ return EXIT_FAILURE;
+ ev.decays[w_pos].at(0).type = static_cast<ParticleID>(
+ ev.decays[w_pos].at(0).type+10 );
+ ev.decays[w_pos].at(1).type = static_cast<ParticleID>(
+ ev.decays[w_pos].at(1).type-10 );
+
+ // double check that we undid all changes
+ if(!test_event(ev, true))
+ return EXIT_FAILURE;
+
+ // 1->3 decay
+ ev.decays[w_pos].emplace_back(
+ HEJ::Particle{photon, fastjet::PseudoJet(0,0,0,0), {}}
+ );
+ if(!test_event(ev, false))
+ return EXIT_FAILURE;
+ ev.decays[w_pos].pop_back();
+
+ // second decay
+ ev.decays[0] = decay_W(ev.outgoing[0]);
+ ev.decays[0].at(0).type = ev.outgoing[0].type;
+ ev.decays[0].at(1).type = gluon;
+ if(!test_event(ev, false))
+ return EXIT_FAILURE;
+ }
+ }
+
+ return EXIT_SUCCESS;
+}
diff --git a/t/test_scale_arithmetics.cc b/t/test_scale_arithmetics.cc
index 933971e..49f8a09 100644
--- a/t/test_scale_arithmetics.cc
+++ b/t/test_scale_arithmetics.cc
@@ -1,120 +1,123 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <fstream>
#include <random>
#include <algorithm>
#include "LHEF/LHEF.h"
#include "HEJ/EventReweighter.hh"
#include "HEJ/make_RNG.hh"
#include "HEJ/Event.hh"
#include "HEJ/YAMLreader.hh"
#include "HEJ/stream.hh"
constexpr double alpha_s = 0.118;
constexpr double ep = 1e-13;
void dump(HEJ::Event const & ev){
{
LHEF::Writer writer{std::cout};
std::cout << std::setprecision(6);
writer.hepeup = to_HEPEUP(std::move(ev), nullptr);
writer.writeEvent();
}
std::cout << "Rapidity ordering:\n";
for(const auto & part: ev.outgoing()){
std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl;
}
}
void shuffle_particles(HEJ::Event::EventData & ev) {
static std::mt19937_64 ran{0};
// incoming
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
// outgoing (through index)
auto old_outgoing = std::move(ev.outgoing);
std::vector<size_t> idx(old_outgoing.size());
std::iota(idx.begin(), idx.end(), 0);
std::shuffle(begin(idx), end(idx), ran);
ev.outgoing.clear();
ev.outgoing.reserve(old_outgoing.size());
for(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(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(begin(decay.second), end(decay.second), ran);
+ }
}
}
int main(int argn, char** argv){
if(argn != 3){
std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml input_file.lhe\n\n";
return EXIT_FAILURE;
}
HEJ::Config config = HEJ::load_config(argv[1]);
config.scales = HEJ::to_ScaleConfig(
YAML::Load("scales: [H_T, 1 * H_T, 2/2 * H_T, 2*H_T/2, H_T/2*2, H_T/2/2*4, H_T*H_T/H_T]")
);
HEJ::istream in{argv[2]};
LHEF::Reader reader{in};
auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);
HEJ::ScaleGenerator scale_gen{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
};
HEJ::EventReweighter resum{
reader.heprup,
std::move(scale_gen),
to_EventReweighterConfig(config),
*ran
};
size_t i = 0;
while(reader.readEvent()){
++i;
HEJ::Event::EventData data{reader.hepeup};
shuffle_particles(data);
HEJ::Event event{
data.cluster(
config.resummation_jets.def,
config.resummation_jets.min_pt
)
};
auto resummed = resum.reweight(event, config.trials);
for(auto && ev: resummed) {
for(auto &&var: ev.variations()) {
if(std::abs(var.muf - ev.central().muf) > ep) {
std::cerr
<< std::setprecision(15)
<< "unequal scales: " << var.muf
<< " != " << ev.central().muf << '\n'
<< "in resummed event:\n";
dump(ev);
std::cerr << "\noriginal event:\n";
dump(event);
return EXIT_FAILURE;
}
}
}
}
return EXIT_SUCCESS;
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:52 PM (9 h, 52 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023707
Default Alt Text
(230 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment