diff --git a/src/Event.cc b/src/Event.cc
index 4e222f7..ef21ce0 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,1257 +1,1259 @@
/**
* \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 <utility>
#include "LHEF/LHEF.h"
#include "fastjet/JetDefinition.hh"
#include "HEJ/Constants.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/PDG_codes.hh"
#define oldcode 0
#define NPRINT 0
namespace HEJ{
namespace{
constexpr int status_in = -1;
constexpr int status_decayed = 2;
constexpr int status_out = 1;
/// @name helper functions to determine event type
//@{
#if oldcode
/**
* \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 has_AWZH_boson = false;
for(auto const & out: outgoing){
if(is_AWZH_boson(out.type)){
if(has_AWZH_boson) return false;
has_AWZH_boson = true;
}
else if(! is_parton(out.type)) return false;
}
return true;
}
#endif
template<class Iterator>
Iterator remove_AWZH(Iterator begin, Iterator end){
return std::remove_if(
begin, end, [](Particle const & p){return is_AWZH_boson(p);}
);
}
template<class Iterator>
bool valid_outgoing(Iterator begin, Iterator end){
return std::distance(begin, end) >= 2
&& std::is_sorted(begin, end, rapidity_less{})
&& std::count_if(
begin, end, [](Particle const & s){return is_AWZH_boson(s);}
) < 2;
}
/**
* \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_Current
*/
bool is_Wp_Current(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_Current(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 is_Pure_Current(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;
+ }
+
+
#if oldcode
/**
* \brief function which determines if type change is consistent with W emission.
* @param in incoming Particle id
* @param out outgoing Particle id
* @param W W boson ID.
* @param qqx Current both incoming/both outgoing?
*
* Ensures that change type of quark line is possible by a flavour changing
* W emission. Allows checking of qqx currents also.
*/
bool is_W_Current(ParticleID in, ParticleID out, ParticleID W, bool qqx){
return ((W==pid::Wp && is_Wp_Current(in,out,qqx))
|| (W==pid::Wm && is_Wm_Current(in,out,qqx)));
}
- /**
- * \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 is_Pure_Current(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;
- }
#endif
// @note that this changes the outgoing range!
template<class ConstIterator, class Iterator>
bool is_FKL(
ConstIterator begin_incoming, ConstIterator end_incoming,
Iterator begin_outgoing, Iterator end_outgoing
){
assert(std::distance(begin_incoming, end_incoming) == 2);
assert(std::distance(begin_outgoing, end_outgoing) >= 2);
// One photon, W, H, Z in the final state is allowed.
// Remove it for remaining tests,
end_outgoing = remove_AWZH(begin_outgoing, end_outgoing);
if(std::all_of(
begin_outgoing + 1, end_outgoing - 1,
[](Particle const & p){ return p.type == pid::gluon; })
){
// Test if this is a standard FKL configuration.
if ( is_Pure_Current(begin_incoming->type, begin_outgoing->type, false)
&& is_Pure_Current((end_incoming-1)->type, (end_outgoing-1)->type, false)){
return true;
}
}
return false;
}
template<class ConstIterator, class Iterator>
bool is_W_FKL(
ConstIterator begin_incoming, ConstIterator end_incoming,
Iterator begin_outgoing, Iterator end_outgoing, ParticleID W
){
assert(std::distance(begin_incoming, end_incoming) == 2);
assert(std::distance(begin_outgoing, end_outgoing) >= 2);
// One photon, W, H, Z in the final state is allowed.
// Remove it for remaining tests,
end_outgoing = remove_AWZH(begin_outgoing, end_outgoing);
if(std::all_of(
begin_outgoing + 1, end_outgoing - 1,
[](Particle const & p){ return p.type == pid::gluon; })
){
// Test if this is a standard FKL configuration.
if( is_W_Current(begin_incoming->type, begin_outgoing->type, W, false)
&& is_Pure_Current((end_incoming-1)->type, (end_outgoing-1)->type, false)){
return true;
}
else if( is_Pure_Current(begin_incoming->type, begin_outgoing->type, false)
&& is_W_Current((end_incoming-1)->type, (end_outgoing-1)->type, W, false)){
return true;
}
}
return false;
}
#if oldcode
bool is_FKL(
std::array<Particle, 2> const & incoming,
std::vector<Particle> outgoing
){
const auto WEmit = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return abs(s.type) == pid::Wp; }
);
if (WEmit != end(outgoing)){
return is_W_FKL(
begin(incoming), end(incoming),
begin(outgoing), end(outgoing), WEmit->type
);
}
else{
return is_FKL(
begin(incoming), end(incoming),
begin(outgoing), end(outgoing)
);
}
}
#endif
bool has_2_jets(Event const & event){
return event.jets().size() >= 2;
}
#if oldcode
/**
* \brief Checks whether event is unordered.
* @returns Is Event Unordered (default back, reverse input for forwards)
*
* - Checks there is more than 3 constuents in the final state
* - Checks there is more than 3 jets
* - Checks the most backwards parton is a gluon
* - Checks the most forwards jet is not a gluon
* - Checks the rest of the event is FKL
* - Checks the second most backwards is not a different boson
* - Checks the unordered gluon actually forms a jet
* - Checks currents at both ends of FKL chain logical to process.
*/
template<class InIterator, class OutIterator, class JetIterator, class IndIterator>
bool is_unordered(
InIterator begin_in, InIterator end_in,
OutIterator begin_out, OutIterator end_out,
JetIterator begin_jets, JetIterator end_jets,
IndIterator begin_indices
){
if(std::distance(begin_out, (end_out)) < 3) return false;
if(std::distance(begin_jets, (end_jets)) < 3) return false;
if((begin_in[0]).type == pid::gluon) return false;
if(((begin_out)->type)==pid::Higgs) return false;
int bosonBegin = (is_AWZ_boson(begin_out->type))?1:0;
int bosonSecond = (is_AWZ_boson((begin_out+bosonBegin+1)->type))?1+bosonBegin:bosonBegin;
if((begin_out+bosonBegin)->type != pid::gluon) return false;
// When skipping the unordered emission
// the remainder should be a regular FKL event,
// except that the (new) first outgoing particle must not be a A,W,Z,H.
if(((begin_out+1)->type)==pid::Higgs) return false;
if(!is_FKL({begin_in[0], (end_in-1)[0]}, {(begin_out+bosonSecond+1), end_out})) return false;
// check that the unordered gluon forms an extra jet
if(!(begin_indices[0] >= 0 && (begin_indices+1)[0] == -1)) return false;
const auto W = std::find_if(
begin_out, end_out,
[](Particle const & s){ return abs(s.type) == pid::Wp; }
);
int BosonEnd = (HEJ::is_AWZH_boson((end_out-1)->type))?1:0;
if (W != end_out){
// Must have a valid W Current
return ((is_Pure_Current(begin_in->type, (begin_out+bosonSecond+1)->type, false)
&& is_W_Current((end_in-1)->type, (end_out-1-BosonEnd)->type, W->type, false))
||( is_W_Current(begin_in->type, (begin_out+bosonSecond+1)->type, W->type, false)
&& is_Pure_Current((end_in-1)->type, (end_out-1-BosonEnd)->type, false)));
}
return ( is_Pure_Current(begin_in->type, (begin_out+bosonSecond+1)->type, false)
&& is_Pure_Current((end_in-1)->type, (end_out-1-BosonEnd)->type, false));
}
/**
* \brief Checks whether event is unordered backwards
* @param ev Event
* @returns Is Event a backwards unordered emission
*
* \see is_unordered()
*/
bool is_unordered_backward(Event const & ev){
assert(std::is_sorted(std::begin(ev.incoming()), std::end(ev.incoming()), pz_less{}));
assert(valid_outgoing(std::begin(ev.outgoing()), std::end(ev.outgoing())));
return is_unordered(cbegin(ev.incoming()), cend(ev.incoming()),
begin(ev.outgoing()), end(ev.outgoing()),
cbegin(ev.jets()), cend(ev.jets()),
cbegin(ev.particle_jet_indices({ev.jets().front()})));
}
#endif
#if oldcode
/**
* \brief Checks for a forward unordered gluon emission
* @param ev Event
* @returns Is the event a forward unordered emission
*
* \see is_unordered()
*/
bool is_unordered_forward(Event const & ev){
assert(std::is_sorted(std::begin(ev.incoming()), std::end(ev.incoming()), pz_less{}));
assert(valid_outgoing(std::begin(ev.outgoing()), std::end(ev.outgoing())));
return is_unordered(crbegin(ev.incoming()), crend(ev.incoming()),
rbegin(ev.outgoing()), rend(ev.outgoing()),
crbegin(ev.jets()), crend(ev.jets()),
crbegin(ev.particle_jet_indices({ev.jets().back()})));
}
#endif
#if oldcode
/**
* \brief Checks for an extremal qqx (forwards, reverse input for back)
* @returns Is the event a forward extremal qqx event
*
* Checks there is 3 or more than 3 constituents in the final state
* Checks there is 3 or more than 3 jets
* Checks most forwards incoming is gluon
* Checks most extremal particle is not a Higgs (either direction)
* Checks the second most forwards particle is not Higgs boson
* Checks remaining chain is pure gluon
* Checks the most forwards parton is a either quark or anti-quark.
* Checks the second most forwards parton is anti-quark or quark.
* Checks the qqbar pair form 2 separate jets.
* Checks other current logically consistent (with/out W emission.)
*/
template<class InIterator, class OutIterator, class JetIterator, class IndIterator>
bool is_Ex_qqx(
InIterator begin_in, InIterator end_in,
OutIterator begin_out, OutIterator end_out,
JetIterator begin_jets, JetIterator end_jets,
IndIterator back_qqx_indice, IndIterator forward_qqx_indice
){
int fkl_start=is_AWZ_boson(begin_out->type);
int fkl_end=0;
int qBosonqx=0;
if(std::distance(begin_out, (end_out)) < 3) return false;
if(std::distance(begin_jets, (end_jets)) < 3) return false;
if((end_in-1)->type != pid::gluon) return false;
if((end_out-1)->type == pid::Higgs || begin_out->type == pid::Higgs
|| (end_out-2)->type == pid::Higgs) return false;
if(is_AWZ_boson((end_out-1)->type)){ // if extremal AWZ
++fkl_end;
if (is_quark((end_out-2)->type)){ //if second quark
if (!(is_antiquark((end_out-3)->type))) return false;// third must be anti-quark
}
else if (is_antiquark((end_out-2)->type)){ //if second anti-quark
if (!(is_quark((end_out-3)->type))) return false;// third must be quark
}
else return false;
}
else if (is_quark((end_out-1)->type)){ //if extremal quark
if(is_AWZ_boson((end_out-2)->type)){ // if second AWZ
++qBosonqx;
if (!(is_antiquark((end_out-3)->type))) return false;// third must be anti-quark
}
else if (!(is_antiquark((end_out-2)->type))) return false;// second must be anti-quark
}
else if (is_antiquark((end_out-1)->type)){ //if extremal anti-quark
if(is_AWZ_boson((end_out-2)->type)){ // if second AWZ
++qBosonqx;
if (!(is_quark((end_out-3)->type))) return false;// third must be quark
}
else if (!(is_quark((end_out-2)->type))) return false;// second must be quark
}
else return false;
// no other quark emission (first doesn't matter)
if(std::any_of(
begin_out+1+fkl_start, end_out-3-fkl_end-qBosonqx,
[](Particle const & p){ return is_anyquark(p); })
) return false;
// Ensure qqbar pair are in jets
if((back_qqx_indice)[0] != 0 || (forward_qqx_indice)[1] != 0){
return false;
}
const auto W = std::find_if(
begin_out, end_out,
[](Particle const & s){ return abs(s.type) == pid::Wp; }
);
if (W != end_out){
// Must have a valid W Current
return ((is_Pure_Current(begin_in->type, (begin_out+fkl_start)->type, false)
&& is_W_Current((end_out-1-fkl_end)->type, (end_out-2-fkl_end-qBosonqx)->type, W->type, true))
||( is_W_Current(begin_in->type, (begin_out+fkl_start)->type, W->type, false)
&& is_Pure_Current((end_out-1-fkl_end)->type, (end_out-2-fkl_end-qBosonqx)->type, true)));
}
return ( is_Pure_Current(begin_in->type, (begin_out+fkl_start)->type, false)
&& is_Pure_Current((end_out-1-fkl_end)->type, (end_out-2-fkl_end-qBosonqx)->type, true));
}
/**
* \brief Checks for a forward extremal qqx
* @param ev Event
* @returns Is the event a forward extremal qqx event
*
* \see is_Ex_qqx()
*/
bool is_Ex_qqxf(Event const & ev){
assert(std::is_sorted(std::begin(ev.incoming()), std::end(ev.incoming()), pz_less{}));
assert(valid_outgoing(std::begin(ev.outgoing()), std::end(ev.outgoing())));
return is_Ex_qqx(cbegin(ev.incoming()), cend(ev.incoming()),
cbegin(ev.outgoing()), cend(ev.outgoing()),
cbegin(ev.jets()), cend(ev.jets()),
crbegin(ev.particle_jet_indices({ev.jets()[ev.jets().size()-1]})),
crbegin(ev.particle_jet_indices({ev.jets()[ev.jets().size()-2]})));
}
#endif
#if oldcode
/**
* \brief Checks for a backward extremal qqx
* @param ev Event
* @returns Is the event a backward extremal qqx event
*
* \see is_Ex_qqx()
*/
bool is_Ex_qqxb(Event const & ev){
assert(std::is_sorted(std::begin(ev.incoming()), std::end(ev.incoming()), pz_less{}));
assert(valid_outgoing(std::begin(ev.outgoing()), std::end(ev.outgoing())));
return is_Ex_qqx(crbegin(ev.incoming()), crend(ev.incoming()),
crbegin(ev.outgoing()), crend(ev.outgoing()),
crbegin(ev.jets()), crend(ev.jets()),
cbegin(ev.particle_jet_indices({ev.jets()[0]})),
cbegin(ev.particle_jet_indices({ev.jets()[1]})));
}
#endif
#if oldcode
/**
* \brief Checks for a central qqx
* @param ev Event
* @returns Is the event a central extremal qqx event
*
* Checks there is 4 or more than 4 constuents in the final state
* Checks there is 4 or more than 4 jets
* Checks most extremal particle is not a Higgs (either direction) y
* Checks for a central quark in the outgoing states
* Checks for adjacent anti-quark parton. (allowing for AWZ boson emission between)
* Check that other partons are only gluons in chain
* Checks external currents are logically sound.
*/
bool is_Mid_qqx(Event const & ev){
auto const & in = ev.incoming();
auto const & out = ev.outgoing();
assert(std::is_sorted(begin(in), end(in), pz_less{}));
assert(valid_outgoing(begin(out), end(out)));
if(out.size() < 4) return false;
if(ev.jets().size() < 4) return false;
if(out.back().type == pid::Higgs || out.front().type == pid::Higgs)
return false;
size_t start_FKL=0;
size_t end_FKL=0;
if (is_AWZ_boson(out.back())) ++end_FKL;
if (is_AWZ_boson(out.front())) ++start_FKL;
const auto WEmit = std::find_if(
begin(out), end(out),
[](Particle const & s){ return abs(s.type) == pid::Wp; }
);
const auto & jets = ev.jets();
const auto indices = ev.particle_jet_indices({jets});
auto const out_partons = filter_partons(out);
HEJ::ParticleID backQuark = HEJ::pid::gluon;
HEJ::ParticleID forwardQuark = HEJ::pid::gluon;
// search for qqx pair
for (size_t i = 1; i<out_partons.size()-2; ++i){
if ((is_quark(out_partons[i]) && (is_antiquark(out_partons[i+1])))
|| (is_antiquark(out_partons[i]) && (is_quark(out_partons[i+1])))
){
// no quarks before qqx (beside first)
if(std::any_of(
out_partons.cbegin()+1, out_partons.cbegin()+i,
[](Particle const & p){ return is_anyquark(p); })
) return false;
// no quarks after qqx (beside last)
if(std::any_of(
out_partons.cbegin()+i+2, out_partons.cend()-1,
[](Particle const & p){ return is_anyquark(p); })
) return false;
// should be in seperate jets
if (indices[i+1] == indices[i]+1 && indices[i] != -1){
backQuark = out_partons[i].type;
forwardQuark = out_partons[i+1].type;
break;}
}
}
if (backQuark == forwardQuark) return false;
if (WEmit != end(out)){
//One W Current Must Exist.
if ( is_W_Current(in.back().type,out.rbegin()[end_FKL].type, WEmit->type, false)
&& is_Pure_Current(backQuark, forwardQuark, true)
&& is_Pure_Current(in.front().type,out[start_FKL].type, false)){
//Backwards Current is W, Nothing to do
}
else if ( is_Pure_Current(in.back().type,out.rbegin()[end_FKL].type, false)
&& is_Pure_Current(backQuark, forwardQuark, true)
&& is_W_Current(in.front().type,out[start_FKL].type, WEmit->type, false)){
//Forwards Current is W, Nothing to do
}
else if ( !( is_Pure_Current(in.back().type,out.rbegin()[end_FKL].type, false)
&& is_W_Current(backQuark, forwardQuark, WEmit->type, true)
&& is_Pure_Current(in.front().type,out[start_FKL].type, false))){
//No W-coupled fermion line, with W present.
return false;
}
}
else{
//Can only be Pure Currents.
return ( is_Pure_Current(in.back().type,out.rbegin()[end_FKL].type, false)
&& is_Pure_Current(backQuark, forwardQuark, true)
&& is_Pure_Current(in.front().type,out[start_FKL].type, false));
}
return true;
}
#endif
using event_type::EventType;
/**
* \brief Checks for all event types
* @param ev Event
* @returns Event Type
*
*/
EventType classify(Event const & ev){
bool IFfuno(false),IFfgqq(false), IFfLL(false),
IFbuno(false), IFbgqq(false), IFbLL(false);
bool IFfunoWp(false),IFfgqqWp(false), IFfLLWp(false),
IFbunoWp(false), IFbgqqWp(false), IFbLLWp(false);
bool IFfunoWm(false),IFfgqqWm(false), IFfLLWm(false),
IFbunoWm(false), IFbgqqWm(false), IFbLLWm(false);
bool Midqqbar(false), MidqqbarWp(false), MidqqbarWm(false);
#if NPRINT
std::cerr <<"got started\n"<<std::endl;
#endif
// if(! final_state_ok(ev.outgoing()))
// return EventType::bad_final_state;
if(! has_2_jets(ev))
return EventType::no_2_jets;
// now check event parton by parton
// check first the impact factors at either end, and loop over the middle partons
auto const & in = ev.incoming();
auto const & out = filter_partons(ev.outgoing());
auto outN=out.size();
assert(std::distance(begin(in), end(in)) == 2);
assert(outN >= 2);
assert(std::distance(begin(out), end(out)) >= 2);
assert(std::is_sorted(begin(out), end(out), rapidity_less{}));
long unsigned int thispartonN=0;
long unsigned int nextpartonN=1;
- if (in.front().type == out.at(thispartonN).type) { // standard jet LL vertex
+ if (is_Pure_Current(in.front().type, out.at(thispartonN).type, false)) { // standard jet LL vertex
IFbLL=true;
++thispartonN;
++nextpartonN;
} else if (is_Wp_Current(in.front().type, out.at(thispartonN).type, false)) {
// if in SU(2)-up and out SU(2)-down then LLWp-type
IFbLLWp=true;
++thispartonN;
++nextpartonN;
} else if (is_Wm_Current(in.front().type, out.at(thispartonN).type, false)) {
// if in SU(2)-down and out SU(2)-up then LLWm-type
IFbLLWm=true;
++thispartonN;
++nextpartonN;
} else if (outN>=3) {
// Possibility of unordered gluon emission or gluon splitting
if ((in.front().type!=pid::gluon)&&out.at(thispartonN).type==pid::gluon) {
// Possibility of unordered gluon emission
- if (in.front().type == out.at(nextpartonN).type) { // unordered jet vertex
+ if (is_Pure_Current(in.front().type, out.at(nextpartonN).type, false)) { // unordered jet vertex
IFbuno=true;
thispartonN+=2;
nextpartonN+=2;
}
else if (is_Wp_Current(in.front().type,out.at(nextpartonN).type, false)) {
// if in SU(2)-up and out SU(2)-down then unoWp-type
IFbunoWp=true;
thispartonN+=2;
nextpartonN+=2;
}
else if (is_Wm_Current(in.front().type, out.at(nextpartonN).type, false)) {
// if in SU(2)-down and out SU(2)-up then unoWm-type
IFbunoWm=true;
thispartonN+=2;
nextpartonN+=2;
}
} else if (in.front().type==pid::gluon) {
// possibility of backward gluon splitting
- if (out.at(thispartonN).type==-out.at(nextpartonN).type) {
+ if (is_Pure_Current(out.at(thispartonN).type, out.at(nextpartonN).type, true)) {
IFbgqq=true;
thispartonN+=2;
nextpartonN+=2;
} else if (is_Wp_Current(out.at(thispartonN).type, out.at(nextpartonN).type, true)) {
// if in SU(2)-up and out SU(2)-down then gqqWp-type
IFbgqqWp=true;
thispartonN+=2;
nextpartonN+=2;
} else if (is_Wm_Current(out.at(thispartonN).type,out.at(nextpartonN).type, true)) {
// if in SU(2)-down and out SU(2)-up then gqqWm-type
IFbgqqWm=true;
thispartonN+=2;
nextpartonN+=2;
}
}
}
// if we reached this far and thisparton==begin(out)
// then it is not a valid all-order state
if (thispartonN==0) return EventType::FixedOrder;
#if NPRINT
std::cerr <<"got this far\n"<<std::endl;
std::cerr << "IFbLLWp : "<<IFbLLWp<<std::endl;
std::cerr << "IFbunoWp : "<<IFbunoWp<<std::endl;
std::cerr << "IFbLL : "<<IFbLL<<std::endl;
#endif
// else check forward particles.
// Need to ensure the same particles don't count for both sub-leading
// impact factors
long unsigned int thispartonfN=outN-1;
long unsigned int prevpartonfN=outN-2;
- if (in.back().type == out.at(thispartonfN).type) {
+ if (is_Pure_Current(in.back().type, out.at(thispartonfN).type, false)) {
// standard forward jet LL vertex
IFfLL=true;
--thispartonfN;
--prevpartonfN;
} else if (is_Wp_Current(in.back().type, out.at(thispartonfN).type, false)) {
// if in SU(2)-up and out SU(2)-down then LLWp-type
IFfLLWp=true;
--thispartonfN;
--prevpartonfN;
} else if (is_Wm_Current(in.back().type, out.at(thispartonfN).type, false)) {
// if in SU(2)-down and out SU(2)-up then LLWm-type
IFfLLWm=true;
--thispartonfN;
--prevpartonfN;
} else if (outN>=3&&(prevpartonfN!=(thispartonN-1))) { // if the prevparton is not
// the last parton included in backward IF
// Possibility of unordered gluon emission or gluon splitting
// unless we already used up all the partons
if ((in.back().type!=pid::gluon)&&out.at(thispartonfN).type==pid::gluon) {
#if NPRINT
std::cerr <<"got this far2\n"<<std::endl;
std::cerr << "thispartonfN, prevpartonfN : "<<thispartonfN <<", "<<prevpartonfN<<std::endl;
std::cerr <<"out.at(prevpartonfN).type : "<<out.at(prevpartonfN).type<<std::endl;
std::cerr <<"out.at(thispartonfN).type : "<<out.at(thispartonfN).type<<std::endl;
#endif
// possibility of unordered emissions
- if (in.back().type == out.at(prevpartonfN).type) { // unordered jet vertex
+ if (is_Pure_Current(in.back().type, out.at(prevpartonfN).type, false)) { // unordered jet vertex
IFfuno=true;
thispartonfN-=2;
prevpartonfN-=2;
}
else if (is_Wp_Current(in.back().type, out.at(prevpartonfN).type, false)) {
// if in SU(2)-up and out SU(2)-down then unoWp-type
IFfunoWp=true;
thispartonfN-=2;
prevpartonfN-=2;
}
else if (is_Wm_Current(in.back().type, out.at(prevpartonfN).type, false)) {
// if in SU(2)-down and out SU(2)-up then unoWm-type
IFfunoWm=true;
thispartonfN-=2;
prevpartonfN-=2;
}
} else if (in.back().type==pid::gluon) {
// possibility of forward gluon splitting
- if (out.at(thispartonfN).type==-out.at(prevpartonfN).type) {
+ if (is_Pure_Current(out.at(thispartonfN).type, out.at(prevpartonfN).type, true)) {
IFfgqq=true;
thispartonfN-=2;
prevpartonfN-=2;
} else if (is_Wp_Current(out.at(thispartonfN).type,out.at(prevpartonfN).type, true)) {
// if in SU(2)-up and out SU(2)-down then gqqWp-type
IFfgqqWp=true;
thispartonfN-=2;
prevpartonfN-=2;
} else if (is_Wm_Current(out.at(thispartonfN).type,out.at(prevpartonfN).type, true)) {
// if in SU(2)-down and out SU(2)-up then gqqWm-type
IFfgqqWm=true;
thispartonfN-=2;
prevpartonfN-=2;
#if NPRINT
std::cerr << "thispartonfN : "<<thispartonfN<<std::endl;
#endif
}
}
} //impact factors done
// if we reached this far and thispartonf==end(out),
// then it is not a valid all-order state
#if NPRINT
std::cerr << "IFfLL : "<<IFfLL<<std::endl;
std::cerr << "IFfLLWp : "<<IFfLLWp<<std::endl;
std::cerr << "IFfLLWm : "<<IFfLLWm<<std::endl;
std::cerr << "IFfuno : "<<IFfuno<<std::endl;
std::cerr << "IFfgqqWm : "<<IFfgqqWm<<std::endl;
#endif
#if NPRINT
std::cerr << "thispartonfN : "<<thispartonfN<<std::endl;
#endif
if (thispartonfN==outN-1) return EventType::FixedOrder;
// Now check the remaining middle partons
// the first (in rapidity) parton in the forward impact factor
long unsigned int firstforwardIFpartonN = thispartonfN+1;
// allow only one qqbar(+Wp/Wm) insertion
while ((out.at(thispartonN).type==pid::gluon)&&(thispartonN!=firstforwardIFpartonN)) {
++thispartonN;
++nextpartonN;
#if NPRINT
std::cerr<<"one gluon found\n";
#endif
}
#if NPRINT
std::cerr<<"thisparton : "<<out.at(thispartonN).type<<std::endl<<"firstforwardIFpartonN : "<<out.at(firstforwardIFpartonN).type<<std::endl;
std::cerr<< "(thispartonN!=firstforwardIFpartonN) : "<<(thispartonN!=firstforwardIFpartonN) <<std::endl;
// std::cerr<< "(thisparton==firstforwardIFparton) : "<<((*thisparton)==(*firstforwardIFparton)) <<std::endl;
#endif
if (thispartonN!=firstforwardIFpartonN) {
// is this a qqbar-pair?
if (out.at(thispartonN).type==-out.at(nextpartonN).type) {
Midqqbar=true;
thispartonN+=2;
nextpartonN+=2;
} else if (is_Wp_Current(out.at(thispartonN).type,out.at(nextpartonN).type, true)) {
// if in SU(2)-up and out SU(2)-down then gqqWp-type
MidqqbarWp=true;
thispartonN+=2;
nextpartonN+=2;
} else if (is_Wm_Current(out.at(thispartonN).type,out.at(nextpartonN).type, true)) {
// if in SU(2)-down and out SU(2)-up then gqqWm-type
MidqqbarWm=true;
thispartonN+=2;
nextpartonN+=2;
}
//continue checking
while ((out.at(thispartonN).type==pid::gluon)&&(thispartonN!=firstforwardIFpartonN)) {
++thispartonN;
++nextpartonN;
}
}
// are we at the end of the partons? if not, this is not resummable
if (thispartonN!=firstforwardIFpartonN) return EventType::FixedOrder;
#if NPRINT
std::cerr<<!(Midqqbar||MidqqbarWm||MidqqbarWp)<<std::endl;
#endif
// Checks using the non-partonic momenta
auto const boson=filter_AWZH_bosons(ev.outgoing());
int WpN(0),WmN(0),HN(0),ZN(0);
if(boson.size()>0){ // if some bosons were found
// count number of W+, W-, H, Z/gamma
for (auto bitr=boson.begin(); bitr!=boson.end(); ++bitr) {
if (bitr->type==24) ++WpN;
else if (bitr->type==-24) ++WmN;
else if (bitr->type==25) ++HN;
else if (bitr->type==23) ++ZN;
}
}
#if NPRINT
std::cerr<<WpN<<" "<<WmN<<" "<<HN<<" "<<ZN<<std::endl;
#endif
// veto Higgs+uno with Higgs mixing in uno particles
// veto Higgs between central qqbar pair
if (HN==1) {
if (IFbuno) { // if the Higgs is in-between the partons of the backward uno impact factor
for (auto bitr=boson.begin(); bitr!=boson.end(); ++bitr) {
if (bitr->type==25) {
if ((bitr->p.rapidity()>out.at(0).p.rapidity())&&(bitr->p.rapidity()<out.at(1).p.rapidity()))
return EventType::FixedOrder;
}
}
}
if (IFfuno) { // if the Higgs is in-between the partons of the backward uno impact factor
for (auto bitr=boson.begin(); bitr!=boson.end(); ++bitr) {
if (bitr->type==25) {
if ((bitr->p.rapidity()<out.at(outN-1).p.rapidity())&&(bitr->p.rapidity()>out.at(outN-2).p.rapidity()))
return EventType::FixedOrder;
}
}
}
}
// Check whether the right number of Ws are present
// count the numbers of Wp Wm required by the identification of partons in the events
int NWpR=IFbLLWp+IFfLLWp+MidqqbarWp+IFbunoWp+IFfunoWp+IFbgqqWp+IFfgqqWp;
int NWmR=IFbLLWm+IFfLLWm+MidqqbarWm+IFbunoWm+IFfunoWm+IFbgqqWm+IFfgqqWm;
#if NPRINT
std::cerr<<NWpR<<" "<<NWmR<<std::endl;
#endif
if (NWpR!=WpN) return EventType::FixedOrder;
if (NWmR!=WmN) return EventType::FixedOrder;
if (WmN>1||WpN>1||HN>1||ZN>1) return EventType::FixedOrder;
// return the event type
if (!(Midqqbar||MidqqbarWm||MidqqbarWp)) { // if there are no mid qqbar
if ((IFbLL||IFbLLWm||IFbLLWp)&&IFfLL) return EventType::FKL;
if ((IFfLLWm||IFfLLWp)&&IFbLL) return EventType::FKL;
if ((IFbuno||IFbunoWm||IFbunoWp)&&IFfLL) return EventType::unordered_backward;
if ((IFbuno)&&(IFfLLWp||IFfLLWm)) return EventType::unordered_backward;
if ((IFfuno||IFfunoWm||IFfunoWp)&&IFbLL) return EventType::unordered_forward;
if ((IFfuno)&&(IFbLLWp||IFbLLWm)) return EventType::unordered_forward;
if ((IFbgqq||IFbgqqWp||IFbgqqWm)&&IFfLL) return EventType::extremal_qqxb;
if ((IFbgqq)&&(IFfLLWp||IFfLLWm)) return EventType::extremal_qqxb;
if ((IFfgqq||IFfgqqWp||IFfgqqWm)&&IFbLL) return EventType::extremal_qqxf;
if ((IFfgqq)&&(IFbLLWp||IFbLLWm)) return EventType::extremal_qqxf;
} else { // there is a middle qqbar pair
if (IFbLL)
if (IFfLL||IFfLLWp||IFfLLWm) return EventType::central_qqx;
if (IFfLL)
if (IFbLL||IFbLLWp||IFbLLWm) return EventType::central_qqx;
}
return EventType::FixedOrder;
}
#if oldcode
EventType classify_old(Event const & ev){
if(! final_state_ok(ev.outgoing()))
return EventType::bad_final_state;
if(! has_2_jets(ev))
return EventType::no_2_jets;
if(is_FKL(ev.incoming(), ev.outgoing()))
return EventType::FKL;
if(is_unordered_backward(ev))
return EventType::unordered_backward;
if(is_unordered_forward(ev))
return EventType::unordered_forward;
if(is_Ex_qqxb(ev))
return EventType::extremal_qqxb;
if(is_Ex_qqxf(ev))
return EventType::extremal_qqxf;
if(is_Mid_qqx(ev))
return EventType::central_qqx;
return EventType::FixedOrder;
}
//@}
#endif
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) {
HEJ::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 HEJ::not_implemented{"lepton-flavour violating final state"};
}
assert(is_neutrino(leptons[1]));
// charged antilepton + neutrino means we had a W+
decayed_boson.type = HEJ::pid::Wp;
}
else if(pidsum == -1) {
assert(is_antilepton(leptons[0]));
if(is_neutrino(leptons[1])) {
throw HEJ::not_implemented{"lepton-flavour violating final state"};
}
assert(is_antineutrino(leptons[0]));
// charged lepton + antineutrino means we had a W-
decayed_boson.type = HEJ::pid::Wm;
}
else {
throw HEJ::not_implemented{
"final state with leptons "
+ HEJ::name(leptons[0].type)
+ " and "
+ HEJ::name(leptons[1].type)
};
}
return decayed_boson;
}
}
void HEJ::Event::EventData::reconstruct_intermediate() {
const auto begin_leptons = std::partition(
begin(outgoing), end(outgoing),
[](HEJ::Particle const & p) {return !HEJ::is_anylepton(p);}
);
if(begin_leptons == end(outgoing)) return;
assert(is_anylepton(*begin_leptons));
std::vector<HEJ::Particle> leptons(begin_leptons, end(outgoing));
outgoing.erase(begin_leptons, end(outgoing));
if(leptons.size() != 2) {
throw HEJ::not_implemented{"Final states with one or more than two leptons"};
}
std::sort(
begin(leptons), end(leptons),
[](HEJ::Particle const & p0, HEJ::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);
// if (classify_new(ev)!=ev.type_) {
// std::cout<< classify_new(ev) << " "<<ev.type_<<std::endl;
// std::cout<<ev;
// }
// assert (classify_new(ev)==ev.type_);
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 {
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_HEJ(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
}
// Connect last
connect_incoming(incoming_[1], anti_colour, colour);
return true;
} // generate_colours
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);
}
}
std::ostream& operator<<(std::ostream & os, Event const & ev){
const std::streamsize orig_prec = os.precision();
os <<std::setprecision(4)<<std::fixed;
std::cout << "########## " << event_type::name(ev.type()) << " ##########" << std::endl;
std::cout << "Incoming particles:\n";
for(auto const & in: ev.incoming()){
std::cout <<std::setw(3)<< in.type << ": ";
print_momentum(os, in.p);
std::cout << std::endl;
}
std::cout << "\nOutgoing particles: " << ev.outgoing().size() << "\n";
for(auto const & out: ev.outgoing()){
std::cout <<std::setw(3)<< out.type << ": ";
print_momentum(os, out.p);
std::cout << " => rapidity="
<<std::setw(7)<<std::right<< out.rapidity() << std::endl;
}
std::cout << "\nForming Jets: " << ev.jets().size() << "\n";
for(auto const & jet: ev.jets()){
print_momentum(os, jet);
std::cout << " => rapidity="
<<std::setw(7)<<std::right<< jet.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;
result.IDPRUP = event.type(); // event type
// the following entries are pretty much meaningless
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
for(Particle const & in: event.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;
}
}