Page MenuHomeHEPForge
Files F8723943

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/PhaseSpacePoint.hh b/include/HEJ/PhaseSpacePoint.hh
index dbfc131..2ced4da 100644
--- a/include/HEJ/PhaseSpacePoint.hh
+++ b/include/HEJ/PhaseSpacePoint.hh
@@ -1,148 +1,149 @@
/** \file
* \brief Contains the PhaseSpacePoint Class
*/
#pragma once
#include <vector>
#include "HEJ/utility.hh"
#include "HEJ/config.hh"
#include "HEJ/Event.hh"
#include "HEJ/JetSplitter.hh"
#include "HEJ/RNG.hh"
namespace HEJ{
//! A point in resummation phase space
class PhaseSpacePoint{
public:
//! Default PhaseSpacePoint Constructor
PhaseSpacePoint() = default;
//! PhaseSpacePoint Constructor
/**
* @param ev Clustered Jet Event
* @param conf Configuration parameters
* @param ran Random number generator
*/
PhaseSpacePoint(
Event const & ev,
PhaseSpacePointConfig conf,
HEJ::RNG & ran
);
//! Get phase space point weight
double weight() const{
return weight_;
}
//! Access incoming particles
std::array<Particle, 2> const & incoming() const{
return incoming_;
}
//! Access outgoing particles
std::vector<Particle> const & outgoing() const{
return outgoing_;
}
//! Particle decays
/**
* The key in the returned map corresponds to the index in the
* vector returned by outgoing()
*/
std::unordered_map<size_t, std::vector<Particle>> const & decays() const{
return decays_;
}
static constexpr int ng_max = 1000; //< maximum number of extra gluons
private:
std::vector<fastjet::PseudoJet> cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const;
bool pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const;
bool pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const;
int sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets);
int sample_ng_jets(int ng, std::vector<fastjet::PseudoJet> const & Born_jets);
double probability_in_jet(
std::vector<fastjet::PseudoJet> const & Born_jets
) const;
std::vector<fastjet::PseudoJet> gen_non_jet(
int ng_non_jet,
double ptmin, double ptmax
);
void rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
);
std::vector<fastjet::PseudoJet> reshuffle(
std::vector<fastjet::PseudoJet> const & Born_jets,
fastjet::PseudoJet const & q
);
bool jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const;
void reconstruct_incoming(std::array<Particle, 2> const & Born_incoming);
double phase_space_normalisation(
int num_Born_jets,
int num_res_partons
) const;
std::vector<fastjet::PseudoJet> split(
std::vector<fastjet::PseudoJet> const & jets, int ng_jets
);
std::vector<int> distribute_jet_partons(
int ng_jets, std::vector<fastjet::PseudoJet> const & jets
);
std::vector<fastjet::PseudoJet> split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np_in_jet
);
bool split_preserved_jets(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<fastjet::PseudoJet> const & jet_partons
) const;
template<class Particle>
Particle const & most_backward_FKL(
std::vector<Particle> const & partons
) const;
template<class Particle>
Particle const & most_forward_FKL(
std::vector<Particle> const & partons
) const;
template<class Particle>
Particle & most_backward_FKL(std::vector<Particle> & partons) const;
template<class Particle>
Particle & most_forward_FKL(std::vector<Particle> & partons) const;
bool extremal_ok(
std::vector<fastjet::PseudoJet> const & partons
) const;
+ void label_qqx(Event const & event);
void copy_AWZH_boson_from(Event const & event);
bool momentum_conserved() const;
- bool unob_, unof_;
+ bool unob_, unof_, qqxb_, qqxf_, qqxmid_;
double weight_;
PhaseSpacePointConfig param_;
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
//! \internal Particle decays in the format {outgoing index, decay products}
std::unordered_map<size_t, std::vector<Particle>> decays_;
std::reference_wrapper<HEJ::RNG> ran_;
};
}
diff --git a/include/HEJ/event_types.hh b/include/HEJ/event_types.hh
index 7a2353c..1aee605 100644
--- a/include/HEJ/event_types.hh
+++ b/include/HEJ/event_types.hh
@@ -1,69 +1,77 @@
/** \file
* \brief Define different types of events.
*
*/
#pragma once
#include "HEJ/utility.hh"
namespace HEJ{
//! Namespace for event types
namespace event_type{
//! Possible event types
enum EventType: size_t{
FKL, /**< FKL-type event */
unordered_backward, /**< event with unordered backward emission */
unordered_forward, /**< event with unordered forward emission */
extremal_qqxb, /**< event with a backward extremal qqbar */
extremal_qqxf, /**< event with a forward extremal qqbar */
central_qqx, /**< event with a central qqbar */
nonHEJ, /**< event configuration not covered by HEJ */
no_2_jets, /**< event with less than two jets */
bad_final_state, /**< event with an unsupported final state */
unob = unordered_backward,
unof = unordered_forward,
qqxexb = extremal_qqxb,
qqxexf = extremal_qqxf,
qqxmid = central_qqx,
first_type = FKL,
last_type = bad_final_state
};
//! Event type names
/**
* For example, names[FKL] is the string "FKL"
*/
static constexpr auto names = make_array(
"FKL",
"unordered backward",
"unordered forward",
"extremal qqbar backward",
"extremal qqbar forward",
"central qqbar",
"nonHEJ",
"no 2 jets",
"bad final state"
);
inline
bool is_HEJ(EventType type) {
switch(type) {
case FKL:
case unordered_backward:
case unordered_forward:
+ case extremal_qqxb:
+ case extremal_qqxf:
+ case central_qqx:
return true;
default:
return false;
}
}
inline
bool is_uno(EventType type) {
return type == unordered_backward || type == unordered_forward;
}
+ inline
+ bool is_qqx(EventType type) {
+ return type == extremal_qqxb || type == extremal_qqxf || type == central_qqx;
+ }
+
}
}
diff --git a/src/Event.cc b/src/Event.cc
index fe72e65..1caf80f 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,636 +1,650 @@
#include "HEJ/Event.hh"
#include "HEJ/utility.hh"
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
//@{
/**
* \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;
}
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 W emission.
* @param in incoming Particle
* @param out outgoing Particle
*
* Ensures that change type of quark line is possible by a flavour changing
* W emission.
*/
bool is_W_Current(ParticleID in, ParticleID out){
if((in==1 && out==2)||(in==2 && out==1)){
return true;
}
else if((in==-1 && out==-2)||(in==-2 && out==-1)){
return true;
}
else if((in==3 && out==4)||(in==4 && out==3)){
return true;
}
else if((in==-3 && out==-4)||(in==-4 && out==-3)){
return true;
}
else{
return false;
}
}
/**
* \brief checks if particle type remains same from incoming to outgoing
* @param in incoming Particle
* @param out outgoing Particle
*/
bool is_Pure_Current(ParticleID in, ParticleID out){
if(abs(in)<=6 || in==21) return (in==out);
else return false;
}
// @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)
&& is_Pure_Current((end_incoming-1)->type, (end_outgoing-1)->type)){
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
){
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)
&& is_Pure_Current((end_incoming-1)->type, (end_outgoing-1)->type)){
return true;
}
else if(is_Pure_Current(begin_incoming->type, begin_outgoing->type)
&& is_W_Current((end_incoming-1)->type, (end_outgoing-1)->type)){
return true;
}
}
return false;
}
bool is_FKL(
std::array<Particle, 2> const & incoming,
std::vector<Particle> outgoing
){
assert(std::is_sorted(begin(incoming), end(incoming), pz_less{}));
assert(valid_outgoing(begin(outgoing), end(outgoing)));
const auto WEmit = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return abs(s.type) == pid::Wp; }
);
if (abs(WEmit->type) == pid::Wp){
return is_W_FKL(
begin(incoming), end(incoming),
begin(outgoing), end(outgoing)
);
}
else{
return is_FKL(
begin(incoming), end(incoming),
begin(outgoing), end(outgoing)
);
}
}
bool has_2_jets(Event const & event){
return event.jets().size() >= 2;
}
/**
* \brief Checks whether event is unordered backwards
* @param ev Event
* @returns Is Event Unordered Backwards
*
* - 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
*/
bool is_unordered_backward(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() < 3) return false;
if(ev.jets().size() < 3) return false;
if(in.front().type == pid::gluon) return false;
if(out.front().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.
const auto FKL_begin = next(begin(out));
if(is_AWZH_boson(*FKL_begin)) return false;
if(!is_FKL(in, {FKL_begin, end(out)})) return false;
// check that the unordered gluon forms an extra jet
const auto jets = sorted_by_rapidity(ev.jets());
const auto indices = ev.particle_jet_indices({jets.front()});
return indices[0] >= 0 && indices[1] == -1;
}
/**
* \brief Checks for a forward unordered gluon emission
* @param ev Event
* @returns Is the event a forward unordered emission
*
* \see is_unordered_backward
*/
bool is_unordered_forward(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() < 3) return false;
if(ev.jets().size() < 3) return false;
if(in.back().type == pid::gluon) return false;
if(out.back().type != pid::gluon) return false;
// When skipping the unordered emission
// the remainder should be a regular FKL event,
// except that the (new) last outgoing particle must not be a A,W,Z,H.
const auto FKL_end = prev(end(out));
if(is_AWZH_boson(*prev(FKL_end))) return false;
if(!is_FKL(in, {begin(out), FKL_end})) return false;
// check that the unordered gluon forms an extra jet
const auto jets = sorted_by_rapidity(ev.jets());
const auto indices = ev.particle_jet_indices({jets.back()});
return indices.back() >= 0 && indices[indices.size()-2] == -1;
}
/**
* \brief Checks for a forward extremal qqx
* @param ev Event
* @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 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.
*/
bool is_Ex_qqxf(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)));
int fkl_end=2;
if(out.size() < 3) return false;
if(ev.jets().size() < 3) return false;
if(in.back().type != pid::gluon) return false;
if(out.back().type == pid::Higgs || out.front().type == pid::Higgs
|| out.rbegin()[1].type == pid::Higgs) return false;
// if extremal AWZ
if(is_AWZ_boson(out.back())){ // if extremal AWZ
fkl_end++;
if (is_quark(out.rbegin()[1])){ //if second quark
if (!(is_antiquark(out.rbegin()[2]))) return false;// third must be anti-quark
}
else if (is_antiquark(out.rbegin()[1])){ //if second anti-quark
if (!(is_quark(out.rbegin()[2]))) return false;// third must be quark
}
else return false;
}
else if (is_quark(out.rbegin()[0])){ //if extremal quark
if(is_AWZ_boson(out.rbegin()[1])){ // if second AWZ
fkl_end++;
if (!(is_antiquark(out.rbegin()[2]))) return false;// third must be anti-quark
}
else if (!(is_antiquark(out.rbegin()[1]))) return false;// second must be anti-quark
}
else if (is_antiquark(out.rbegin()[0])){ //if extremal anti-quark
if(is_AWZ_boson(out.rbegin()[1])){ // if second AWZ
fkl_end++;
if (!(is_quark(out.rbegin()[2]))) return false;// third must be quark
}
else if (!(is_quark(out.rbegin()[1]))) return false;// second must be quark
}
else return false;
// When skipping the qqbar
// New last outgoing particle must not be a Higgs
if (out.rbegin()[fkl_end].type == pid::Higgs) return false;
+ const auto jets = fastjet::sorted_by_rapidity(ev.jets());
+ const auto indices = ev.particle_jet_indices({jets});
+
+ // Ensure qqbar pair are in separate jets
+ if(indices[indices.size()-2] != indices[indices.size()-1]-1) return false;
+
// Opposite current should be logical to process
if (is_AWZ_boson(out.front().type)){
return (is_Pure_Current(in.front().type, out[1].type)
|| is_W_Current(in.front().type,out[1].type));
}
else
return (is_Pure_Current(in.front().type, out[0].type)
|| is_W_Current(in.front().type,out[0].type));
}
/**
* \brief Checks for a backward extremal qqx
* @param ev Event
* @returns Is the event a backward 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 backwards incoming is gluon
* Checks most extremal particle is not a Higgs (either direction) y
* Checks the second most backwards particle is not Higgs boson y
* Checks the most backwards parton is a either quark or anti-quark. y
* Checks the second most backwards parton is anti-quark or quark. y
+ * Checks the qqbar pair form 2 separate jets.
*/
bool is_Ex_qqxb(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)));
int fkl_start=2;
if(out.size() < 3) return false;
if(ev.jets().size() < 3) return false;
if(in.front().type != pid::gluon) return false;
if(out.back().type == pid::Higgs || out.front().type == pid::Higgs
|| out[1].type == pid::Higgs) return false;
if(is_AWZ_boson(out.front())){ // if extremal AWZ
fkl_start++;
if (is_quark(out[1])){ //if second quark
if (!(is_antiquark(out[2]))) return false;// third must be anti-quark
}
else if (is_antiquark(out[1])){ //if second anti-quark
if (!(is_quark(out[2]))) return false;// third must be quark
}
else return false;
}
else if (is_quark(out[0])){ // if extremal quark
if(is_AWZ_boson(out[1])){ // if second AWZ
fkl_start++;
if (!(is_antiquark(out[2]))) return false;// third must be anti-quark
}
else if (!(is_antiquark(out[1]))) return false;// second must be anti-quark
}
else if (is_antiquark(out[0])){ //if extremal anti-quark
if(is_AWZ_boson(out[1])){ // if second AWZ
fkl_start++;
if (!(is_quark(out[2]))) return false;// third must be quark
}
else if (!(is_quark(out[1]))) return false;// second must be quark
}
else return false;
// When skipping the qqbar
// New last outgoing particle must not be a Higgs.
if (out[fkl_start].type == pid::Higgs) return false;
+ const auto jets = fastjet::sorted_by_rapidity(ev.jets());
+ const auto indices = ev.particle_jet_indices({jets});
+
+ // Ensure qqbar pair form separate jets.
+ if(indices[0] != indices[1]-1) return false;
+
// Other current should be logical to process
if (is_AWZ_boson(out.back())){
return (is_Pure_Current(in.back().type, out.rbegin()[1].type)
|| is_W_Current(in.back().type,out.rbegin()[1].type));
}
else
return (is_Pure_Current(in.back().type, out.rbegin()[0].type)
|| is_W_Current(in.back().type, out.rbegin()[0].type));
}
/**
* \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)
* 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++;
}
if ((is_Pure_Current(in.back().type,out.rbegin()[end_FKL].type)
&& is_Pure_Current(in.front().type,out[start_FKL].type))){
//nothing to do
}
else if (is_W_Current(in.back().type,out.rbegin()[end_FKL].type)
&& is_Pure_Current(in.front().type,out[start_FKL].type)){
//nothing to do
}
else if (!(is_Pure_Current(in.back().type,out.rbegin()[end_FKL].type)
&& is_W_Current(in.front().type,out[start_FKL].type))){
return false;
}
for(size_t i=1+start_FKL; i<out.size()-1-end_FKL;i++){
if (is_quark(out[i])){
if ((is_antiquark(out[i-1]) && i!=1)
|| (is_antiquark(out[i+1]) && i!=out.size()-1-end_FKL))
return true;
else if (is_AWZ_boson(out[i-1]) && (is_antiquark(out[i-2]) && i!=2) )
return true;
else if (is_AWZ_boson(out[i+1]) && (is_antiquark(out[i+2]) && i!=out.size()-2) )
return true;
}
}
return false;
}
using event_type::EventType;
EventType classify(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::nonHEJ;
}
//@}
Particle extract_particle(LHEF::HEPEUP const & hepeup, int i){
return Particle{
static_cast<ParticleID>(hepeup.IDUP[i]),
fastjet::PseudoJet{
hepeup.PUP[i][0], hepeup.PUP[i][1],
hepeup.PUP[i][2], hepeup.PUP[i][3]
}
};
}
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
UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup):
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"
+ "Event has too many incoming particles"
};
}
incoming[in_idx++] = std::move(particle);
}
else outgoing.emplace_back(std::move(particle));
}
std::sort(
begin(incoming), end(incoming),
[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
);
std::sort(begin(outgoing), end(outgoing), rapidity_less{});
// 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 ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
):
ev_{std::move(ev)},
cs_{to_PseudoJet(filter_partons(ev_.outgoing)), jet_def},
min_jet_pt_{min_jet_pt}
{
type_ = classify(*this);
}
std::vector<fastjet::PseudoJet> Event::jets() const{
return cs_.inclusive_jets(min_jet_pt_);
}
/**
* \brief Returns the invarient mass of the event
* @param ev Event
* @returns s hat
*
* Makes use of the FastJet PseudoJet function m2().
* Applies this function to the sum of the incoming partons.
*/
double shat(Event const & ev){
return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
}
namespace{
// colour flow according to Les Houches standard
// TODO: stub
std::vector<std::pair<int, int>> colour_flow(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing
){
std::vector<std::pair<int, int>> result(
incoming.size() + outgoing.size()
);
for(auto & col: result){
col = std::make_pair(-1, -1);
}
return result;
}
}
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()+1; // event ID
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;
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);
}
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);
}
result.ICOLUP = colour_flow(
event.incoming(), filter_partons(event.outgoing())
);
if(result.ICOLUP.size() < num_particles){
const size_t AWZH_boson_idx = std::find_if(
begin(event.outgoing()), end(event.outgoing()),
[](Particle const & s){ return is_AWZH_boson(s); }
) - begin(event.outgoing()) + event.incoming().size();
assert(AWZH_boson_idx <= result.ICOLUP.size());
result.ICOLUP.insert(
begin(result.ICOLUP) + AWZH_boson_idx,
std::make_pair(0,0)
);
}
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 dd254c6..67e4c84 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1655 +1,1650 @@
#include "HEJ/MatrixElement.hh"
#include <unordered_map>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
#include "HEJ/Constants.hh"
#include "HEJ/currents.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/utility.hh"
namespace HEJ{
//cf. last line of eq. (22) in \ref Andersen:2011hs
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j, double lambda
) const {
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,
bool check_momenta
) const {
return tree(event, check_momenta)*virtual_corrections(event);
}
Weights MatrixElement::tree(
Event const & event,
bool check_momenta
) const {
if(! is_HEJ(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
return tree_param(event)*tree_kin(event, check_momenta);
}
Weights MatrixElement::tree_param(
Event const & event
) const {
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_HEJ(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(
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
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 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){
q -= out[1].p;
++first_idx;
}
if(out.back().type == pid::Higgs || event.type() == event_type::unof){
--last_idx;
}
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, CLAMBDA)*(
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
);
return exp(exponent);
}
} // namespace HEJ
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector 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));
// cout << "#Fadin qa : "<<qav<<endl;
// cout << "#Fadin qb : "<<qbv<<endl;
// cout << "#Fadin p1 : "<<p1<<endl;
// cout << "#Fadin p2 : "<<p2<<endl;
// cout << "#Fadin p5 : "<<p5<<endl;
// cout << "#Fadin Gauge Check : "<< CL.dot(p5)<<endl;
// cout << "#Fadin C2L : "<< -CL.dot(CL)<<" "<<-CL.dot(CL)/(qav.m2()*qbv.m2())/(4./p5.perp2())<<endl;
// TODO can this dead test go?
// if (-CL.dot(CL)<0.)
// if (fabs(CL.dot(p5))>fabs(CL.dot(CL))) // not sufficient!
// return 0.;
// else
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
double kperp=(qav-qbv).perp();
if (kperp>HEJ::CLAMBDA)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
return Cls-4./(kperp*kperp);
}
}
//! Lipatov vertex
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip,
CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B
{
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 qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
double kperp=(qav-qbv).perp();
if (kperp>HEJ::CLAMBDA)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
else {
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 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 jM2gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2qg(pn,pb,p1,pa);
else
return jM2qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return jM2qg(p1,pa,pn,pb);
else
return jM2qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2qQ(pn,pb,p1,pa);
else
return jM2qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return jM2qQbar(p1,pa,pn,pb);
else
return jM2qbarQbar(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
* @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
){
// We know it cannot be gg incoming.
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jMWqg(pn,pl,plbar,pb,p1,pa);
else
return jMWqbarg(pn,pl,plbar,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return jMWqg(p1,pl,plbar,pa,pn,pb);
else
return jMWqbarg(p1,pl,plbar,pa,pn,pb);
}
else { // they are both quark
if (wc==true){ // emission off b, (first argument pbout)
if (bptype>0) {
if (aptype>0)
return jMWqQ(pn,pl,plbar,pb,p1,pa);
else
return jMWqQbar(pn,pl,plbar,pb,p1,pa);
}
else {
if (aptype>0)
return jMWqbarQ(pn,pl,plbar,pb,p1,pa);
else
return jMWqbarQbar(pn,pl,plbar,pb,p1,pa);
}
}
else{ // emission off a, (first argument paout)
if (aptype > 0) {
if (bptype > 0)
return jMWqQ(p1,plbar,pl,pa,pn,pb);
else
return jMWqQbar(p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0)
return jMWqbarQ(p1,plbar,pl,pa,pn,pb);
else
return jMWqbarQbar(p1,plbar,pl,pa,pn,pb);
}
}
}
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
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*/
double ME_W_unob_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
){
// we know they are not both gluons
if (bptype == 21 && aptype != 21) { // b gluon => W emission off a
if (aptype > 0)
return jM2Wunogqg(pg,p1,plbar,pl,pa,pn,pb);
else
return jM2Wunogqbarg(pg,p1,plbar,pl,pa,pn,pb);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return junobMWqQg(pn,plbar,pl,pb,p1,pa,pg);
else
return junobMWqQbarg(pn,plbar,pl,pb,p1,pa,pg);
}
else{
if (aptype>0)
return junobMWqbarQg(pn,plbar,pl,pb,p1,pa,pg);
else
return junobMWqbarQbarg(pn,plbar,pl,pb,p1,pa,pg);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return jM2WunogqQ(pg,p1,plbar,pl,pa,pn,pb);
else //qqbar
return jM2WunogqQbar(pg,p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return jM2WunogqbarQ(pg,p1,plbar,pl,pa,pn,pb);
else //qbarqbar
return jM2WunogqbarQbar(pg,p1,plbar,pl,pa,pn,pb);
}
}
}
}
/** Matrix element squared for uno forward 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
* @returns ME Squared for unof Tree-Level Current-Current Scattering
*/
double ME_W_unof_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
){
// we know they are not both gluons
if (aptype==21 && bptype!=21) {//a gluon => W emission off b
if (bptype > 0)
return jM2Wunogqg(pg, pn,plbar, pl, pb, p1, pa);
else
return jM2Wunogqbarg(pg, pn,plbar, pl, pb, p1, pa);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return jM2WunogqQ(pg,pn,plbar,pl,pb,p1,pa);
else
return jM2WunogqQbar(pg,pn,plbar,pl,pb,p1,pa);
}
else{
if (aptype>0)
return jM2WunogqbarQ(pg,pn,plbar,pl,pb,p1,pa);
else
return jM2WunogqbarQbar(pg,pn,plbar,pl,pb,p1,pa);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return junofMWgqQ(pg,pn,pb,p1,plbar,pl,pa);
else //qqbar
return junofMWgqQbar(pg,pn,pb,p1,plbar,pl,pa);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return junofMWgqbarQ(pg,pn,pb,p1,plbar,pl,pa);
else //qbarqbar
return junofMWgqbarQbar(pg,pn,pb,p1,plbar,pl,pa);
}
}
}
}
/** \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
* @returns ME Squared for qqxb Tree-Level Current-Current Scattering
*/
double ME_W_qqxb_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 wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFbackward;
if (pqbar.rapidity() > pq.rapidity()){
swapQuarkAntiquark=true;
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pq.minus())+(pq.minus())/pa.minus())+1./3.)*3./4.;
}
else{
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pqbar.minus())+(pqbar.minus())/pa.minus())+1./3.)*3./4.;
}
// 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 (swapQuarkAntiquark){
return jM2Wggtoqqbarg(pa, pqbar, plbar, pl, pq, pn,pb)*CFbackward;}
else {
return jM2Wggtoqbarqg(pa, pq, plbar, pl, pqbar, pn,pb)*CFbackward;}
}
else if (aptype==21&&bptype!=21 ) {//a gluon => W emission off b leg or qqx
if (wc!=1){ // W Emitted from backwards qqx
if (swapQuarkAntiquark){
return jM2WgQtoqqbarQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
else{
return jM2WgQtoqbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
}
else { // W Must be emitted from forwards leg.
if(bptype > 0){
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pb, pa, pn, pqbar, pq, plbar, pl, false)*CFbackward;}
else{
return jM2WgqtoQQqW(pb, pa, pn, pq, pqbar, plbar, pl, false)*CFbackward;}
} else {
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pb, pa, pn, pqbar, pq, plbar, pl, true)*CFbackward;}
else{
return jM2WgqtoQQqW(pb, pa, pn, pq, pqbar, plbar, pl, true)*CFbackward;}
}
}
}
else{
throw std::logic_error("Incompatible incoming particle types with qqxb");
}
}
/* \brief Matrix element squared for forward 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 p1 Final state 1 Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @returns ME Squared for qqxf Tree-Level Current-Current Scattering
*/
double ME_W_qqxf_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFforward;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pq.plus())+(pq.plus())/pb.plus())+1./3.)*3./4.;
}
else{
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pqbar.plus())+(pqbar.plus())/pb.plus())+1./3.)*3./4.;
}
// 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 (swapQuarkAntiquark){
return jM2Wggtoqqbarg(pb, pqbar, plbar, pl, pq, p1,pa)*CFforward;}
else {
return jM2Wggtoqbarqg(pb, pq, plbar, pl, pqbar, p1,pa)*CFforward;}
}
else if (bptype==21&&aptype!=21) {// b gluon => W emission off a or qqx
if (wc==1){ // W Emitted from forwards qqx
if (swapQuarkAntiquark){
return jM2WgQtoqbarqQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
else {
return jM2WgQtoqqbarQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
}
// W Must be emitted from backwards leg.
if (aptype > 0){
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pa,pb, p1, pqbar, pq, plbar, pl, false)*CFforward;}
else{
return jM2WgqtoQQqW(pa,pb, p1, pq, pqbar, plbar, pl, false)*CFforward;}
} else
{
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pa,pb, p1, pqbar, pq, plbar, pl, true)*CFforward;}
else{
return jM2WgqtoQQqW(pa,pb, p1, pq, pqbar, plbar, pl, true)*CFforward;}
}
}
else{
throw std::logic_error("Incompatible incoming particle types with qqxf");
}
}
/* \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> partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
bool swapQuarkAntiquark=false;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
}
double CFforward = (0.5*(3.-1./3.)*(pb.plus()/(partons[partons.size()-1].plus())+(partons[partons.size()-1].plus())/pb.plus())+1./3.)*3./4.;
double CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(partons[0].minus())+(partons[0].minus())/pa.minus())+1./3.)*3./4.;
double wt=1.;
if (aptype==21) wt*=CFbackward;
if (bptype==21) wt*=CFforward;
if (aptype <=0 && bptype <=0){ // Both External AntiQuark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons,true,true, swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true,true, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true,true, swapQuarkAntiquark, nabove, nbelow, false);
}
} // end both antiquark
else if (aptype<=0){ // a is antiquark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons, false, true, swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons,false,true, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, true, swapQuarkAntiquark, nabove, nbelow, false);
}
} // end a is antiquark
else if (bptype<=0){ // b is antiquark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow, false);
}
} //end b is antiquark
else{ //Both Quark or gluon
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl, plbar, partons, false, false, swapQuarkAntiquark, nabove);}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow, false);
}
}
}
/** \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
){
if (aptype==21&&bptype==21) // gg initial state
return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered forward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param punof Unordered Particle Momentum
* @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 with Higgs and unordered forward emission
*/
double ME_Higgs_current_unof(
int aptype, int bptype,
CLHEP::HepLorentzVector const & punof,
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
){
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (aptype>0)
return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
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 punob 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
*/
double ME_Higgs_current_unob(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & punob,
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
){
if (bptype==21&&aptype!=21) {
if (aptype > 0)
return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (bptype>0)
return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
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
namespace HEJ{
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,
bool check_momenta
) const {
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, check_momenta);
}
switch(AWZH_boson->type){
case pid::Higgs: {
return tree_kin_Higgs(ev, check_momenta);
}
// TODO
case pid::Wp: {
return tree_kin_W(ev, true, check_momenta);
}
case pid::Wm: {
return tree_kin_W(ev, false, check_momenta);
}
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 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)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
Event const & ev, bool check_momenta
) const{
auto out_partons = filter_partons(ev.outgoing());
if(!check_momenta){
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
// TODO: avoid reclustering
fastjet::ClusterSequence cs(to_PseudoJet(out_partons), ev.jet_def());
const auto jets = sorted_by_rapidity(cs.inclusive_jets(ev.min_jet_pt()));
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
- // skip jets caused by unordered emission
- if(ev.type() == event_type::unob){
+ // 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){
+ else if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
assert(jets.size() >= 3);
--most_forward;
}
const auto extremal_jet_indices = cs.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;
}
double MatrixElement::tree_kin_jets(
Event const & ev,
bool check_momenta
) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev, check_momenta);
if(is_uno(ev.type())){
throw not_implemented("unordered emission not implemented for pure jets");
}
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
);
}
namespace{
double tree_kin_W_FKL(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl, bool WPlus
) {
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
bool wc;
if (aptype==partons[0].type) { //leg b emits w
wc = true;
}
else{
wc = false;
q0 -=pl + plbar;
}
double current_factor;
if (WPlus){
current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, pl, plbar, wc
);
}
else{
current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_unob(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl, bool WPlus
) {
auto pg = to_HepLorentzVector(partons[0]);
auto p1 = to_HepLorentzVector(partons[1]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1- pg;
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc;
if (aptype==partons[1].type) { //leg b emits w
wc = true;
}
else{
wc = false;
q0 -=pl + plbar;
}
double current_factor;
if (WPlus){
current_factor = ME_W_unob_current(
aptype, bptype, pn, pb,
p1, pa, pg, pl, plbar, wc
);
}
else{
current_factor = ME_W_unob_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_unof(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl, bool WPlus
) {
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 2]);
auto pg = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc;
if (aptype==partons[0].type) { //leg b emits w
wc = true;}
else{
wc = false;
q0 -=pl + plbar;
}
double current_factor;
if (WPlus){
current_factor = ME_W_unof_current(
aptype, bptype, pn, pb,
p1, pa, pg, pl, plbar, wc
);
}
else{
current_factor = ME_W_unof_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxb(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl, bool WPlus
) {
HLV pq,pqbar;
if(is_quark(partons[0])){
pq = to_HepLorentzVector(partons[0]);
pqbar = to_HepLorentzVector(partons[1]);
}
else{
pq = to_HepLorentzVector(partons[1]);
pqbar = to_HepLorentzVector(partons[0]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - pq - pqbar;
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc;
if (partons[0].type==-partons[1].type) { //leg b emits w
wc = true;}
else{
wc = false;
q0 -=pl + plbar;
}
double current_factor;
if (WPlus){
current_factor = ME_W_qqxb_current(
aptype, bptype, pa, pb,
pq, pqbar, pn, pl, plbar, wc
);
}
else{
current_factor = ME_W_qqxb_current(
aptype, bptype, pa, pb,
pq, pqbar, pn, plbar, pl, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxf(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl, bool WPlus
) {
HLV pq,pqbar;
if(is_quark(partons[partons.size() - 1])){
pq = to_HepLorentzVector(partons[partons.size() - 1]);
pqbar = to_HepLorentzVector(partons[partons.size() - 2]);
}
else{
pq = to_HepLorentzVector(partons[partons.size() - 2]);
pqbar = to_HepLorentzVector(partons[partons.size() - 1]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc;
if (aptype==partons[0].type) { //leg b emits w
wc = true;}
else{
wc = false;
q0 -=pl + plbar;
}
double current_factor;
if (WPlus){
current_factor = ME_W_qqxf_current(
aptype, bptype, pa, pb,
pq, pqbar, p1, pl, plbar, wc
);
}
else{
current_factor = ME_W_qqxf_current(
aptype, bptype, pa, pb,
pq, pqbar, p1, plbar, pl, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn
);
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, bool WPlus
) {
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
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;
}
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
auto first_after_qqx = (backmidquark+2);
for(auto parton_it = begin_ladder; parton_it != first_after_qqx; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
int nabove = std::distance(begin_ladder, backmidquark);
int nbelow = std::distance(first_after_qqx, 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]));
}
double current_factor;
if (WPlus){
current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, pl, plbar, wqq, wc
);
}
else{
current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, backmidquark,
q0, pa, pb, p1, pn
)*FKL_ladder_weight(
first_after_qqx, end_ladder,
qqxt, pa, pb, p1, pn
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
} // namespace anonymous
double MatrixElement::tree_kin_W(
Event const & ev,
bool WPlus,
bool check_momenta
) const {
using namespace event_type;
auto const & incoming(ev.incoming());
auto const & outgoing(ev.outgoing());
auto const & decays(ev.decays());
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));
}
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto the_W = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return abs(s.type) == pid::Wp; }
);
std::vector<Particle> partons(begin(outgoing), the_W);
partons.insert(end(partons), the_W + 1, end(outgoing));
partons = tag_extremal_jet_partons(ev, check_momenta);
- // if(has_unob_gluon(incoming, outgoing)){
if(ev.type() == unordered_backward){
return tree_kin_W_unob(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl, WPlus);
}
- // if(has_unof_gluon(incoming, outgoing)){
if(ev.type() == unordered_forward){
return tree_kin_W_unof(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl, WPlus);
}
- // if(has_Ex_qqxb(incoming, outgoing)){
if(ev.type() == extremal_qqxb){
return tree_kin_W_qqxb(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl, WPlus);
}
- // if(has_Ex_qqxf(incoming, outgoing)){
if(ev.type() == extremal_qqxf){
return tree_kin_W_qqxf(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl, WPlus);
}
- // if(has_mid_qqx(outgoing)){
if(ev.type() == central_qqx){
return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl, WPlus);
}
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl, WPlus);
}
double MatrixElement::tree_kin_Higgs(
Event const & ev,
bool check_momenta
) const {
if(is_uno(ev.type())){
return tree_kin_Higgs_between(ev, check_momenta);
}
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev, check_momenta);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev, check_momenta);
}
return tree_kin_Higgs_between(ev, check_momenta);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
// TODO: code duplication with currents.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;
}
// 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 p1out, CLHEP::HepLorentzVector p1in,
ParticleID type2,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
// 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*MH2gq_outsideH(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(
Event const & ev,
bool check_momenta
) 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, check_momenta);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
ev, check_momenta
);
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
);
}
double MatrixElement::tree_kin_Higgs_last(
Event const & ev,
bool check_momenta
) 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, check_momenta);
}
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
ev, check_momenta
);
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
);
}
double MatrixElement::tree_kin_Higgs_between(
Event const & ev,
bool check_momenta
) 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, check_momenta);
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() == unob){
current_factor = C_A*C_A/2.*ME_Higgs_current_unob( // 1/2 = "K_uno"
incoming[0].type, incoming[1].type,
pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
current_factor = C_A*C_A/2.*ME_Higgs_current_unof( // 1/2 = "K_uno"
incoming[0].type, incoming[1].type,
to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
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
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double MatrixElement::tree_param_partons(
double alpha_s, double mur,
std::vector<Particle> const & partons
) const{
const double gs2 = 4.*M_PI*alpha_s;
double wt = std::pow(gs2, partons.size());
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(partons.size() >= 2);
for(size_t i = 1; i < partons.size()-1; ++i){
wt *= 1 + alpha_s/(2*M_PI)*beta0*log(mur/partons[i].p.perp());
}
}
return wt;
}
double MatrixElement::tree_param(
Event const & ev,
double mur
) const{
assert(is_HEJ(ev.type()));
const auto & out = ev.outgoing();
const double alpha_s = alpha_s_(mur);
auto AWZH_boson = std::find_if(
begin(out), end(out),
[](auto const & p){return is_AWZH_boson(p);}
);
double AWZH_coupling = 1.;
if(AWZH_boson != end(out)){
switch(AWZH_boson->type){
case pid::Higgs: {
AWZH_coupling = alpha_s*alpha_s;
break;
}
// TODO
case pid::Wp:{
AWZH_coupling = alpha_w*alpha_w/2;
break;
}
case pid::Wm:{
AWZH_coupling = alpha_w*alpha_w/2;
break;
}
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
if(ev.type() == event_type::unob){
return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons(
alpha_s, mur, filter_partons({begin(out) + 1, end(out)})
);
}
if(ev.type() == event_type::unof){
return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons(
alpha_s, mur, filter_partons({begin(out), end(out) - 1})
);
}
return AWZH_coupling*tree_param_partons(alpha_s, mur, filter_partons(out));
}
} // namespace HEJ
diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc
index ff32eb3..54f4d52 100644
--- a/src/PhaseSpacePoint.cc
+++ b/src/PhaseSpacePoint.cc
@@ -1,537 +1,595 @@
#include "HEJ/PhaseSpacePoint.hh"
#include <random>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
#include "HEJ/Constants.hh"
#include "HEJ/resummation_jet.hh"
#include "HEJ/utility.hh"
#include "HEJ/kinematics.hh"
namespace HEJ{
namespace {
constexpr int max_jet_user_idx = PhaseSpacePoint::ng_max;
bool is_nonjet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() > max_jet_user_idx;
}
bool is_jet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() <= max_jet_user_idx;
}
// user indices for partons with extremal rapidity
+ constexpr int qqxb_idx = -7;
+ constexpr int qqxf_idx = -6;
constexpr int unob_idx = -5;
constexpr int unof_idx = -4;
constexpr int backward_FKL_idx = -3;
constexpr int forward_FKL_idx = -2;
}
namespace {
double estimate_ng_mean(std::vector<fastjet::PseudoJet> const & Born_jets){
const double delta_y =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
assert(delta_y > 0);
// Formula derived from fit in arXiv:1805.04446 (see Fig. 2)
return 0.975052*delta_y;
}
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
return cs.inclusive_jets(param_.jet_param.min_pt);
}
bool PhaseSpacePoint::pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const{
return cluster_jets(jets).size() == jets.size();
}
int PhaseSpacePoint::sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets){
const double ng_mean = estimate_ng_mean(Born_jets);
std::poisson_distribution<int> dist(ng_mean);
const int ng = dist(ran_.get());
assert(ng >= 0);
assert(ng < ng_max);
weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng);
return ng;
}
void PhaseSpacePoint::copy_AWZH_boson_from(Event const & event){
auto const & from = event.outgoing();
const auto AWZH_boson = std::find_if(
begin(from), end(from),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson == end(from)) return;
auto insertion_point = std::lower_bound(
begin(outgoing_), end(outgoing_), *AWZH_boson, rapidity_less{}
);
outgoing_.insert(insertion_point, *AWZH_boson);
// copy decay products
const int idx = std::distance(begin(from), AWZH_boson);
assert(idx >= 0);
const auto decay_it = event.decays().find(idx);
if(decay_it != end(event.decays())){
const int new_idx = std::distance(begin(outgoing_), insertion_point);
assert(new_idx >= 0);
assert(outgoing_[new_idx].type == AWZH_boson->type);
decays_.emplace(new_idx, decay_it->second);
}
assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
}
+
+ //! \brief relabels qqx-pair with its PDG IDs.
+ //*@param ev Born Event
+ //
+ // This function will label the qqx pair in a qqx event back to
+ // their original types from the input event.
+ // @TODO: Central qqx only works with 4j input and not higher!
+ void PhaseSpacePoint::label_qqx(Event const & event){
+ auto const & bornout = event.outgoing();
+ const auto backquark = std::find_if(
+ begin(bornout) + 1 - ((qqxb_)?1:0), end(bornout) - 1 + ((qqxf_)?1:0) ,
+ [](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
+ );
+ if(backquark == end(bornout)) {
+ weight_=0;
+ return;
+ }
+ auto partons = to_PseudoJet(filter_partons(outgoing_));
+ fastjet::ClusterSequence cs(partons, event.jet_def());
+ const auto jets = fastjet::sorted_by_rapidity(cs.inclusive_jets(event.min_jet_pt()));
+ assert(jets.size() >= (3u + qqxmid_));
+ const auto indices = cs.particle_jet_indices({jets});
+ int pqjet = 0 + qqxmid_ + qqxf_*(jets.size()-2);
+
+ bool jetsadjacent = false;
+ assert(partons.size() == indices.size());
+ for (size_t i = 0; i<indices.size(); i++){
+ if (indices[i]==(pqjet) && indices[i+1] == pqjet+1){
+ outgoing_.at(i).type = backquark->type;
+ outgoing_.at(i+1).type = (backquark+1)->type;
+ jetsadjacent = true;
+ continue;
+ }
+ }
+ if (!jetsadjacent) weight_=0;
+ assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
+ }
+
PhaseSpacePoint::PhaseSpacePoint(
Event const & ev, PhaseSpacePointConfig conf, HEJ::RNG & ran
):
unob_{ev.type() == event_type::unob},
unof_{ev.type() == event_type::unof},
+ qqxb_{ev.type() == event_type::qqxexb},
+ qqxf_{ev.type() == event_type::qqxexf},
+ qqxmid_{ev.type() == event_type::qqxmid},
param_{std::move(conf)},
ran_{ran}
{
weight_ = 1;
const auto Born_jets = sorted_by_rapidity(ev.jets());
const int ng = sample_ng(Born_jets);
weight_ /= std::tgamma(ng + 1);
const int ng_jets = sample_ng_jets(ng, Born_jets);
std::vector<fastjet::PseudoJet> out_partons = gen_non_jet(
ng - ng_jets, CMINPT, param_.jet_param.min_pt
);
- {
+
const auto qperp = std::accumulate(
begin(out_partons), end(out_partons),
fastjet::PseudoJet{}
);
const auto jets = reshuffle(Born_jets, qperp);
if(weight_ == 0.) return;
if(! pass_resummation_cuts(jets)){
weight_ = 0.;
return;
}
std::vector<fastjet::PseudoJet> jet_partons = split(jets, ng_jets);
if(weight_ == 0.) return;
rescale_rapidities(
out_partons,
most_backward_FKL(jet_partons).rapidity(),
most_forward_FKL(jet_partons).rapidity()
);
if(! cluster_jets(out_partons).empty()){
weight_ = 0.;
return;
}
std::sort(begin(out_partons), end(out_partons), rapidity_less{});
assert(
std::is_sorted(begin(jet_partons), end(jet_partons), rapidity_less{})
);
const auto first_jet_parton = out_partons.insert(
end(out_partons), begin(jet_partons), end(jet_partons)
);
std::inplace_merge(
begin(out_partons), first_jet_parton, end(out_partons), rapidity_less{}
);
- }
+
if(! jets_ok(Born_jets, out_partons)){
weight_ = 0.;
return;
}
weight_ *= phase_space_normalisation(Born_jets.size(), out_partons.size());
outgoing_.reserve(out_partons.size() + 1); // one slot for possible A, W, Z, H
for(auto & p: out_partons){
outgoing_.emplace_back(Particle{pid::gluon, std::move(p)});
}
most_backward_FKL(outgoing_).type = ev.incoming().front().type;
most_forward_FKL(outgoing_).type = ev.incoming().back().type;
+ if(qqxmid_||qqxb_||qqxf_){
+ label_qqx(ev);
+ }
copy_AWZH_boson_from(ev);
assert(!outgoing_.empty());
+
reconstruct_incoming(ev.incoming());
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_non_jet(
int count, double ptmin, double ptmax
){
// heuristic parameters for pt sampling
const double ptpar = 1.3 + count/5.;
const double temp1 = atan((ptmax - ptmin)/ptpar);
std::vector<fastjet::PseudoJet> partons(count);
for(size_t i = 0; i < (size_t) count; ++i){
const double r1 = ran_.get().flat();
const double pt = ptmin + ptpar*tan(r1*temp1);
const double temp2 = cos(r1*temp1);
const double phi = 2*M_PI*ran_.get().flat();
weight_ *= 2.0*M_PI*pt*ptpar*temp1/(temp2*temp2);
// we don't know the allowed rapidity span yet,
// set a random value to be rescaled later on
const double y = ran_.get().flat();
partons[i].reset_PtYPhiM(pt, y, phi);
// Set user index higher than any jet-parton index
// in order to assert that these are not inside jets
partons[i].set_user_index(i + 1 + ng_max);
assert(ptmin-1e-5 <= partons[i].pt() && partons[i].pt() <= ptmax+1e-5);
}
assert(std::all_of(partons.cbegin(), partons.cend(), is_nonjet_parton));
return partons;
}
void PhaseSpacePoint::rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
){
constexpr double ep = 1e-7;
for(auto & parton: partons){
assert(0 <= parton.rapidity() && parton.rapidity() <= 1);
const double dy = ymax - ymin - 2*ep;
const double y = ymin + ep + dy*parton.rapidity();
parton.reset_momentum_PtYPhiM(parton.pt(), y, parton.phi());
weight_ *= dy;
assert(ymin <= parton.rapidity() && parton.rapidity() <= ymax);
}
}
namespace {
template<typename T, typename... Rest>
auto min(T const & a, T const & b, Rest&&... r) {
using std::min;
return min(a, min(b, std::forward<Rest>(r)...));
}
}
double PhaseSpacePoint::probability_in_jet(
std::vector<fastjet::PseudoJet> const & Born_jets
) const{
assert(std::is_sorted(begin(Born_jets), end(Born_jets), rapidity_less{}));
assert(Born_jets.size() >= 2);
const double dy =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
const double R = param_.jet_param.def.R();
const int njets = Born_jets.size();
const double p_J_y_large = (njets-1)*R*R/(2.*dy);
const double p_J_y0 = njets*R/M_PI;
return min(p_J_y_large, p_J_y0, 1.);
}
int PhaseSpacePoint::sample_ng_jets(
int ng, std::vector<fastjet::PseudoJet> const & Born_jets
){
const double p_J = probability_in_jet(Born_jets);
std::binomial_distribution<> bin_dist(ng, p_J);
const int ng_J = bin_dist(ran_.get());
weight_ *= std::pow(p_J, -ng_J)*std::pow(1 - p_J, ng_J - ng);
return ng_J;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::reshuffle(
std::vector<fastjet::PseudoJet> const & Born_jets,
fastjet::PseudoJet const & q
){
if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets;
const auto jets = resummation_jet_momenta(Born_jets, q);
if(jets.empty()){
weight_ = 0;
return {};
}
// additional Jacobian to ensure Born integration over delta gives 1
weight_ *= resummation_jet_weight(Born_jets, q);
return jets;
}
std::vector<int> PhaseSpacePoint::distribute_jet_partons(
int ng_jets, std::vector<fastjet::PseudoJet> const & jets
){
size_t first_valid_jet = 0;
size_t num_valid_jets = jets.size();
const double R_eff = 5./3.*param_.jet_param.def.R();
// if there is an unordered jet too far away from the FKL jets
// then extra gluon constituents of the unordered jet would
// violate the FKL rapidity ordering
- if(unob_ && jets[0].delta_R(jets[1]) > R_eff){
+ if((unob_||qqxb_) && jets[0].delta_R(jets[1]) > R_eff){
++first_valid_jet;
--num_valid_jets;
}
- else if(unof_ && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){
+ else if((unof_||qqxf_) && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){
--num_valid_jets;
}
std::vector<int> np(jets.size(), 1);
for(int i = 0; i < ng_jets; ++i){
++np[first_valid_jet + ran_.get().flat() * num_valid_jets];
}
weight_ *= std::pow(num_valid_jets, ng_jets);
return np;
}
#ifndef NDEBUG
namespace{
bool tagged_FKL_backward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return std::find_if(
begin(jet_partons), end(jet_partons),
[](fastjet::PseudoJet const & p){
return p.user_index() == backward_FKL_idx;
}
) != end(jet_partons);
}
bool tagged_FKL_forward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
// the most forward FKL parton is most likely near the end of jet_partons;
// start search from there
return std::find_if(
jet_partons.rbegin(), jet_partons.rend(),
[](fastjet::PseudoJet const & p){
return p.user_index() == forward_FKL_idx;
}
) != jet_partons.rend();
}
bool tagged_FKL_extremal(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return tagged_FKL_backward(jet_partons) && tagged_FKL_forward(jet_partons);
}
} // namespace anonymous
#endif
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
int ng_jets
){
return split(jets, distribute_jet_partons(ng_jets, jets));
}
bool PhaseSpacePoint::pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const{
if(ext_parton.pt() < param_.min_extparton_pt) return false;
return (ext_parton - jet).pt()/jet.pt() < param_.max_ext_soft_pt_fraction;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np
){
assert(! jets.empty());
assert(jets.size() == np.size());
assert(pass_resummation_cuts(jets));
- const size_t most_backward_FKL_idx = 0 + unob_;
- const size_t most_forward_FKL_idx = jets.size() - 1 - unof_;
+ const size_t most_backward_FKL_idx = 0 + unob_ + qqxb_;
+ const size_t most_forward_FKL_idx = jets.size() - 1 - unof_ - qqxf_;
const auto & jet = param_.jet_param;
const JetSplitter jet_splitter{jet.def, jet.min_pt, ran_};
std::vector<fastjet::PseudoJet> jet_partons;
// randomly distribute jet gluons among jets
for(size_t i = 0; i < jets.size(); ++i){
auto split_res = jet_splitter.split(jets[i], np[i]);
weight_ *= split_res.weight;
if(weight_ == 0) return {};
assert(
std::all_of(
begin(split_res.constituents), end(split_res.constituents),
is_jet_parton
)
);
const auto first_new_parton = jet_partons.insert(
end(jet_partons),
begin(split_res.constituents), end(split_res.constituents)
);
// mark uno and extremal FKL emissions here so we can check
// their position once all emissions are generated
auto extremal = end(jet_partons);
- if((unob_ && i == 0) || i == most_backward_FKL_idx){
- // unordered or FKL backward emission
+ if (i == most_backward_FKL_idx){ //FKL backward emission
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
- extremal->set_user_index(
- (i == most_backward_FKL_idx)?backward_FKL_idx:unob_idx
+ extremal->set_user_index(backward_FKL_idx);
+ }
+ else if(((unob_ || qqxb_) && i == 0)){
+ // unordered/qqxb
+ extremal = std::min_element(
+ first_new_parton, end(jet_partons), rapidity_less{}
);
+ extremal->set_user_index((unob_)?unob_idx:qqxb_idx);
}
- else if((unof_ && i == jets.size() - 1) || i == most_forward_FKL_idx){
- // unordered or FKL forward emission
+
+ else if (i == most_forward_FKL_idx){
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
- extremal->set_user_index(
- (i == most_forward_FKL_idx)?forward_FKL_idx:unof_idx
+ extremal->set_user_index(forward_FKL_idx);
+ }
+ else if(((unof_ || qqxf_) && i == jets.size() - 1)){
+ // unordered/qqxf
+ extremal = std::max_element(
+ first_new_parton, end(jet_partons), rapidity_less{}
);
+ extremal->set_user_index((unof_)?unof_idx:qqxf_idx);
}
if(
extremal != end(jet_partons)
&& !pass_extremal_cuts(*extremal, jets[i])
){
weight_ = 0;
return {};
}
}
assert(tagged_FKL_extremal(jet_partons));
std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
if(
!extremal_ok(jet_partons)
|| !split_preserved_jets(jets, jet_partons)
){
weight_ = 0.;
return {};
}
return jet_partons;
}
bool PhaseSpacePoint::extremal_ok(
std::vector<fastjet::PseudoJet> const & partons
) const{
assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
if(unob_ && partons.front().user_index() != unob_idx) return false;
if(unof_ && partons.back().user_index() != unof_idx) return false;
+ if(qqxb_ && partons.front().user_index() != qqxb_idx) return false;
+ if(qqxf_ && partons.back().user_index() != qqxf_idx) return false;
return
most_backward_FKL(partons).user_index() == backward_FKL_idx
&& most_forward_FKL(partons).user_index() == forward_FKL_idx;
}
bool PhaseSpacePoint::split_preserved_jets(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<fastjet::PseudoJet> const & jet_partons
) const{
assert(std::is_sorted(begin(jets), end(jets), rapidity_less{}));
const auto split_jets = sorted_by_rapidity(cluster_jets(jet_partons));
// this can happen if two overlapping jets
// are both split into more than one parton
if(split_jets.size() != jets.size()) return false;
for(size_t i = 0; i < split_jets.size(); ++i){
// this can happen if there are two overlapping jets
// and a parton is assigned to the "wrong" jet
if(!nearby_ep(jets[i].rapidity(), split_jets[i].rapidity(), 1e-2)){
return false;
}
}
return true;
}
template<class Particle>
Particle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> const & partons
) const{
- return partons[0 + unob_];
+ return partons[0 + unob_ + qqxb_];
}
template<class Particle>
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
- const size_t idx = partons.size() - 1 - unof_;
+ const size_t idx = partons.size() - 1 - unof_ - qqxf_;
assert(idx < partons.size());
return partons[idx];
}
template<class Particle>
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
- return partons[0 + unob_];
+ return partons[0 + unob_ + qqxb_];
}
template<class Particle>
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
- const size_t idx = partons.size() - 1 - unof_;
+ const size_t idx = partons.size() - 1 - unof_ - qqxf_;
assert(idx < partons.size());
return partons[idx];
}
namespace {
bool contains_idx(
fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
){
auto const & constituents = jet.constituents();
const int idx = parton.user_index();
return std::find_if(
begin(constituents), end(constituents),
[idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
) != end(constituents);
}
}
/**
* final jet test:
* - number of jets must match Born kinematics
* - no partons designated as nonjet may end up inside jets
* - all other outgoing partons *must* end up inside jets
* - the extremal (in rapidity) partons must be inside the extremal jets
* - rapidities must be the same (by construction)
*/
bool PhaseSpacePoint::jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
if(jets.size() != Born_jets.size()) return false;
int in_jet = 0;
for(size_t i = 0; i < jets.size(); ++i){
assert(jets[i].has_constituents());
for(auto && parton: jets[i].constituents()){
if(is_nonjet_parton(parton)) return false;
}
in_jet += jets[i].constituents().size();
}
const int expect_in_jet = std::count_if(
partons.cbegin(), partons.cend(), is_jet_parton
);
if(in_jet != expect_in_jet) return false;
// note that PseudoJet::contains does not work here
if(! (
contains_idx(most_backward_FKL(jets), most_backward_FKL(partons))
&& contains_idx(most_forward_FKL(jets), most_forward_FKL(partons))
)) return false;
if(unob_ && !contains_idx(jets.front(), partons.front())) return false;
if(unof_ && !contains_idx(jets.back(), partons.back())) return false;
for(size_t i = 0; i < jets.size(); ++i){
assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2));
}
return true;
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<Particle, 2> const & Born_incoming
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
for(size_t i = 0; i < incoming_.size(); ++i){
incoming_[i].type = Born_incoming[i].type;
}
assert(momentum_conserved());
}
double PhaseSpacePoint::phase_space_normalisation(
int num_Born_jets, int num_out_partons
) const{
return pow(16*pow(M_PI,3), num_Born_jets - num_out_partons);
}
bool PhaseSpacePoint::momentum_conserved() const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
const double norm = diff.E();
for(auto const & out: outgoing()) diff -= out.p;
return nearby(diff, fastjet::PseudoJet{}, norm);
}
} //namespace HEJ

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:58 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242308
Default Alt Text
(109 KB)

Event Timeline