diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
index 77f9f53..741832f 100644
--- a/include/HEJ/Event.hh
+++ b/include/HEJ/Event.hh
@@ -1,391 +1,402 @@
/** \file
* \brief Declares the Event class and helpers
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <cstddef>
#include <iosfwd>
#include <iterator>
#include <unordered_map>
#include <utility>
#include <vector>
#include "boost/iterator/filter_iterator.hpp"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/PseudoJet.hh"
#include "HEJ/Constants.hh"
#include "HEJ/Parameters.hh"
#include "HEJ/Particle.hh"
#include "HEJ/event_types.hh"
namespace LHEF {
class HEPEUP;
class HEPRUP;
}
namespace fastjet {
class JetDefinition;
}
namespace HEJ {
struct RNG;
struct UnclusteredEvent;
/** @brief An event with clustered jets
*
* This is the main HEJ 2 event class.
* It contains kinematic information including jet clustering,
* parameter (e.g. scale) settings and the event weight.
*/
class Event {
public:
class EventData;
//! Iterator over partons
using ConstPartonIterator = boost::filter_iterator<
bool (*)(Particle const &),
std::vector<Particle>::const_iterator
>;
//! Reverse Iterator over partons
using ConstReversePartonIterator = std::reverse_iterator<
ConstPartonIterator>;
//! No default Constructor
Event() = delete;
//! Event Constructor adding jet clustering to an unclustered event
//! @deprecated UnclusteredEvent got superseded by EventData
//! and will be removed in HEJ 2.2.0
[[deprecated("UnclusteredEvent got superseded by EventData")]]
Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
);
//! @name Particle Access
//! @{
//! Incoming particles
std::array<Particle, 2> const & incoming() const{
return incoming_;
}
//! Outgoing particles
std::vector<Particle> const & outgoing() const{
return outgoing_;
}
//! Iterator to the first outgoing parton
ConstPartonIterator begin_partons() const;
//! Iterator to the first outgoing parton
ConstPartonIterator cbegin_partons() const;
//! Iterator to the end of the outgoing partons
ConstPartonIterator end_partons() const;
//! Iterator to the end of the outgoing partons
ConstPartonIterator cend_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator rbegin_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator crbegin_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator rend_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator crend_partons() const;
//! Particle decays
/**
* The key in the returned map corresponds to the index in the
* vector returned by outgoing()
*/
std::unordered_map<std::size_t, std::vector<Particle>> const & decays()
const {
return decays_;
}
//! The jets formed by the outgoing partons, sorted in rapidity
std::vector<fastjet::PseudoJet> const & jets() const{
return jets_;
}
//! @}
//! @name Weight variations
//! @{
//! All chosen parameter, i.e. scale choices (const version)
Parameters<EventParameters> const & parameters() const{
return parameters_;
}
//! All chosen parameter, i.e. scale choices
Parameters<EventParameters> & parameters(){
return parameters_;
}
//! Central parameter choice (const version)
EventParameters const & central() const{
return parameters_.central;
}
//! Central parameter choice
EventParameters & central(){
return parameters_.central;
}
//! Parameter (scale) variations (const version)
std::vector<EventParameters> const & variations() const{
return parameters_.variations;
}
//! Parameter (scale) variations
std::vector<EventParameters> & variations(){
return parameters_.variations;
}
//! Parameter (scale) variation (const version)
/**
* @param i Index of the requested variation
*/
EventParameters const & variations(std::size_t i) const{
return parameters_.variations.at(i);
}
//! Parameter (scale) variation
/**
* @param i Index of the requested variation
*/
EventParameters & variations(std::size_t i){
return parameters_.variations.at(i);
}
//! @}
//! Indices of the jets the outgoing partons belong to
/**
* @param jets Jets to be tested
* @returns A vector containing, for each outgoing parton,
* the index in the vector of jets the considered parton
* belongs to. If the parton is not inside any of the
* passed jets, the corresponding index is set to -1.
*/
std::vector<int> particle_jet_indices(
std::vector<fastjet::PseudoJet> const & jets
) const {
return cs_.particle_jet_indices(jets);
}
//! particle_jet_indices() of the Event jets()
std::vector<int> particle_jet_indices() const {
return particle_jet_indices(jets());
}
//! Jet definition used for clustering
fastjet::JetDefinition const & jet_def() const{
return cs_.jet_def();
}
//! Minimum jet transverse momentum
double min_jet_pt() const{
return min_jet_pt_;
}
//! Event type
event_type::EventType type() const{
return type_;
}
//! Give colours to each particle
/**
* @returns true if new colours are generated, i.e. same as is_resummable()
* @details Colour ordering is done according to leading colour in the MRK
* limit, see \cite Andersen:2011zd. This only affects \ref
* is_resummable() "HEJ" configurations, all other \ref event_type
* "EventTypes" will be ignored.
* @note This overwrites all previously set colours.
*/
bool generate_colours(RNG & /*ran*/);
//! Check that current colours are leading in the high energy limit
/**
* @details Checks that the colour configuration can be split up in
* multiple, rapidity ordered, non-overlapping ladders. Such
* configurations are leading in the MRK limit, see
* \cite Andersen:2011zd
*
* @note This is _not_ to be confused with \ref is_resummable(), however
* for all resummable states it is possible to create a leading colour
* configuration, see generate_colours()
*/
bool is_leading_colour() const;
/**
* @brief Check if given event could have been produced by HEJ
* @details A HEJ state has to fulfil:
* 1. type() has to be \ref is_resummable() "resummable"
* 2. Soft radiation in the tagging jets contributes at most to
* `soft_pt_regulator` of the total jet \f$ p_\perp \f$
*
* @note This is true for any resummed stated produced by the
* EventReweighter or any \ref is_resummable() "resummable" Leading
* Order state.
*
* @param soft_pt_regulator Maximum transverse momentum fraction from soft
* radiation in tagging jets
* @param min_pt Absolute minimal \f$ p_\perp \f$,
* \b deprecated use soft_pt_regulator instead
* @return True if this state could have been produced by HEJ
*/
bool valid_hej_state(
double soft_pt_regulator = DEFAULT_SOFT_PT_REGULATOR,
double min_pt = 0.
) const;
+ //! Check that the incoming momenta are valid
+ /**
+ * @details Checks that the incoming parton momenta are on-shell and have
+ * vanishing transverse components.
+ *
+ */
+ bool valid_incoming() const;
+
private:
//! \internal
//! @brief Construct Event explicitly from input.
/** This is only intended to be called from EventData.
*
* \warning The input is taken _as is_, sorting and classification has to be
* done externally, i.e. by EventData
*/
Event(
std::array<Particle, 2> && incoming,
std::vector<Particle> && outgoing,
std::unordered_map<std::size_t, std::vector<Particle>> && decays,
Parameters<EventParameters> && parameters,
fastjet::JetDefinition const & jet_def,
double min_jet_pt
);
//! Iterator over partons (non-const)
using PartonIterator = boost::filter_iterator<
bool (*)(Particle const &),
std::vector<Particle>::iterator
>;
//! Reverse Iterator over partons (non-const)
using ReversePartonIterator = std::reverse_iterator<PartonIterator>;
//! Iterator to the first outgoing parton (non-const)
PartonIterator begin_partons();
//! Iterator to the end of the outgoing partons (non-const)
PartonIterator end_partons();
//! Reverse Iterator to the first outgoing parton (non-const)
ReversePartonIterator rbegin_partons();
//! Reverse Iterator to the first outgoing parton (non-const)
ReversePartonIterator rend_partons();
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
std::unordered_map<std::size_t, std::vector<Particle>> decays_;
std::vector<fastjet::PseudoJet> jets_;
Parameters<EventParameters> parameters_;
fastjet::ClusterSequence cs_;
double min_jet_pt_;
event_type::EventType type_;
}; // end class Event
//! Detect if a backward incoming gluon turns into a backward outgoing Higgs boson
inline
bool is_backward_g_to_h(Event const & ev) {
return ev.outgoing().front().type == pid::higgs
&& ev.incoming().front().type == pid::gluon;
}
//! Detect if a forward incoming gluon turns into a forward outgoing Higgs boson
inline
bool is_forward_g_to_h(Event const & ev) {
return ev.outgoing().back().type == pid::higgs
&& ev.incoming().back().type == pid::gluon;
}
//! Class to store general Event setup, i.e. Phase space and weights
class Event::EventData {
public:
//! Default Constructor
EventData() = default;
//! Constructor from LesHouches event information
EventData(LHEF::HEPEUP const & hepeup);
//! Constructor with all values given
EventData(
std::array<Particle, 2> incoming,
std::vector<Particle> outgoing,
std::unordered_map<std::size_t, std::vector<Particle>> decays,
Parameters<EventParameters> parameters
):
incoming(std::move(incoming)), outgoing(std::move(outgoing)),
decays(std::move(decays)), parameters(std::move(parameters))
{}
//! Generate an Event from the stored EventData.
/**
* @details Do jet clustering and classification.
* Use this to generate an Event.
*
* @note Calling this function destroys EventData
*
* @param jet_def Jet definition
* @param min_jet_pt minimal \f$p_T\f$ for each jet
*
* @returns Full clustered and classified event.
*/
Event cluster(
fastjet::JetDefinition const & jet_def, double min_jet_pt);
//! Alias for cluster()
Event operator()(
fastjet::JetDefinition const & jet_def, double const min_jet_pt){
return cluster(jet_def, min_jet_pt);
}
//! Sort particles in rapidity
void sort();
//! Reconstruct intermediate particles from final-state leptons
/**
* Final-state leptons are created from virtual photons, W, or Z bosons.
* This function tries to reconstruct such intermediate bosons if they
* are not part of the event record.
*/
void reconstruct_intermediate();
//! Incoming particles
std::array<Particle, 2> incoming;
//! Outcoing particles
std::vector<Particle> outgoing;
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<std::size_t, std::vector<Particle>> decays;
//! Parameters, e.g. scale or inital weight
Parameters<EventParameters> parameters;
}; // end class EventData
//! Print Event
std::ostream& operator<<(std::ostream & os, Event const & ev);
//! Square of the partonic centre-of-mass energy \f$\hat{s}\f$
double shat(Event const & ev);
+ //! Tolerance parameter for validity check on incoming momenta
+ static constexpr double TOL = 1e-6;
+
//! Convert an event to a LHEF::HEPEUP
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * /*heprup*/);
// put deprecated warning at the end, so don't get the warning inside Event.hh,
// additionally doxygen can not identify [[deprecated]] correctly
struct [[deprecated("UnclusteredEvent will be replaced by EventData")]]
UnclusteredEvent;
//! An event before jet clustering
//! @deprecated UnclusteredEvent got superseded by EventData
//! and will be removed in HEJ 2.2.0
struct UnclusteredEvent{
//! Default Constructor
UnclusteredEvent() = default;
//! Constructor from LesHouches event information
UnclusteredEvent(LHEF::HEPEUP const & hepeup);
std::array<Particle, 2> incoming; /**< Incoming Particles */
std::vector<Particle> outgoing; /**< Outgoing Particles */
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<std::size_t, std::vector<Particle>> decays;
//! Central parameter (e.g. scale) choice
EventParameters central;
std::vector<EventParameters> variations; /**< For parameter variation */
};
} // namespace HEJ
diff --git a/include/HEJ/MatrixElement.hh b/include/HEJ/MatrixElement.hh
index e03dc19..37c2e8a 100644
--- a/include/HEJ/MatrixElement.hh
+++ b/include/HEJ/MatrixElement.hh
@@ -1,202 +1,206 @@
/** \file
* \brief Contains the MatrixElement Class
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <functional>
#include <vector>
#include "fastjet/PseudoJet.hh"
#include "HEJ/Config.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Parameters.hh"
namespace CLHEP {
class HepLorentzVector;
}
namespace HEJ {
class Event;
struct Particle;
//! Class to calculate the squares of matrix elements
class MatrixElement{
public:
/** \brief MatrixElement Constructor
* @param alpha_s Function taking the renormalisation scale
* and returning the strong coupling constant
* @param conf General matrix element settings
*/
MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
);
/**
* \brief squares of regulated HEJ matrix elements
* @param event The event for which to calculate matrix elements
* @returns The squares of HEJ matrix elements including virtual corrections
*
* This function returns one value for the central parameter choice
* and one additional value for each entry in \ref Event.variations().
* See eq. (22) in \cite Andersen:2011hs for the definition of the squared
* matrix element.
*
+ * Incoming momenta must be lightlike and have vanishing transverse momentum.
+ *
* \internal Relation to standard HEJ Met2: MatrixElement = Met2*shat^2/(pdfta*pdftb)
*/
Weights operator()(Event const & event) const;
//! Squares of HEJ tree-level matrix elements
/**
* @param event The event for which to calculate matrix elements
* @returns The squares of HEJ matrix elements without virtual corrections
*
* cf. eq. (22) in \cite Andersen:2011hs
*/
Weights tree(Event const & event) const;
/**
* \brief Virtual corrections to matrix element squares
* @param event The event for which to calculate matrix elements
* @returns The virtual corrections to the squares of the matrix elements
*
* The all order virtual corrections to LL in the MRK limit is
* given by replacing 1/t in the scattering amplitude according to the
* lipatov ansatz.
*
* If there is more than one entry in the returned vector, each entry
* corresponds to the contribution from the interference of two
* channels. The order of these entries matches the one returned by
* the tree_kin member function, but is otherwise unspecified.
*
* cf. second-to-last line of eq. (22) in \cite Andersen:2011hs
* note that indices are off by one, i.e. out[0].p corresponds to p_1
*/
std::vector<Weights> virtual_corrections(Event const & event) const;
/**
* \brief Scale-dependent part of tree-level matrix element squares
* @param event The event for which to calculate matrix elements
* @returns The scale-dependent part of the squares of the
* tree-level matrix elements
*
* The tree-level matrix elements factorises into a renormalisation-scale
* dependent part, given by the strong coupling to some power, and a
* scale-independent remainder. This function only returns the former parts
* for the central scale choice and all \ref Event.variations().
*
* @see tree, tree_kin
*/
Weights tree_param(Event const & event) const;
/**
* \brief Kinematic part of tree-level matrix element squares
* @param event The event for which to calculate matrix elements
* @returns The kinematic part of the squares of the
* tree-level matrix elements
*
* The tree-level matrix elements factorises into a renormalisation-scale
* dependent part, given by the strong coupling to some power, and a
* scale-independent remainder. This function only returns the latter part.
*
* If there is more than one entry in the returned vector, each entry
* corresponds to the contribution from the interference of two
* channels. The order of these entries matches the one returned by
* the virtual_corrections member function, but is otherwise unspecified.
*
+ * Incoming momenta must be lightlike and have vanishing transverse momentum.
+ *
* @see tree, tree_param
*/
std::vector<double> tree_kin(Event const & event) const;
private:
double tree_param(
Event const & event,
double mur
) const;
double virtual_corrections_W(
Event const & event,
double mur,
Particle const & WBoson
) const;
std::vector <double> virtual_corrections_Z_qq(
Event const & event,
double mur,
Particle const & ZBoson
) const;
double virtual_corrections_Z_qg(
Event const & event,
double mur,
Particle const & ZBoson,
bool is_gq_event
) const;
std::vector<double> virtual_corrections(
Event const & event,
double mur
) const;
//! \internal cf. last line of eq. (22) in \cite Andersen:2011hs
double omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const;
double tree_kin_jets(
Event const & ev
) const;
double tree_kin_W(
Event const & ev
) const;
std::vector <double> tree_kin_Z(
Event const & ev
) const;
double tree_kin_Higgs(
Event const & ev
) const;
double tree_kin_Higgs_first(
Event const & ev
) const;
double tree_kin_Higgs_last(
Event const & ev
) const;
double tree_kin_Higgs_between(
Event const & ev
) const;
double tree_param_partons(
double alpha_s, double mur,
std::vector<Particle> const & partons
) const;
std::vector<int> in_extremal_jet_indices(
std::vector<fastjet::PseudoJet> const & partons
) const;
double MH2_backwardH(
ParticleID type_forward,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pH,
CLHEP::HepLorentzVector const & pn
) const;
double MH2_unob_forwardH(
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pH
) const;
std::function<double (double)> alpha_s_;
MatrixElementConfig param_;
};
} // namespace HEJ
diff --git a/src/Event.cc b/src/Event.cc
index 35d38a0..7296958 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,1218 +1,1228 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/Event.hh"
#include <algorithm>
#include <cassert>
#include <cstdlib>
#include <iomanip>
#include <iterator>
#include <memory>
#include <numeric>
#include <ostream>
#include <string>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/PseudoJet.hh"
#include "LHEF/LHEF.h"
#include "HEJ/Constants.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/RNG.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/optional.hh"
+#include "HEJ/utility.hh"
namespace HEJ {
namespace {
using std::size_t;
//! LHE status codes
namespace lhe_status {
enum Status: int {
in = -1,
decay = 2,
out = 1,
};
}
using LHE_Status = lhe_status::Status;
//! true if leptonic W decay
bool valid_W_decay( int const w_type, // sign of W
std::vector<Particle> const & decays
){
if(decays.size() != 2) // no 1->2 decay
return false;
const int pidsum = decays[0].type + decays[1].type;
if( std::abs(pidsum) != 1 || pidsum != w_type ) // correct charge
return false;
// leptonic decay (only check first, second follows from pidsum)
if( w_type == 1 ) // W+
return is_antilepton(decays[0]) || is_neutrino(decays[0]);
// W-
return is_lepton(decays[0]) || is_antineutrino(decays[0]);
}
//! true for Z decay to charged leptons
bool valid_Z_decay(std::vector<Particle> const & decays){
if(decays.size() != 2) // no 1->2 decay
return false;
const int pidsum = decays[0].type + decays[1].type;
if( std::abs(pidsum) != 0 ) // correct charge
return false;
// leptonic decay (only check first, second follows from pidsum)
return is_anylepton(decays[0]) && !is_anyneutrino(decays[0]);
}
/// @name helper functions to determine event type
//@{
/**
* \brief check if final state valid for HEJ
*
* check if there is at most one photon, W, H, Z in the final state
* and all the rest are quarks or gluons
*/
bool final_state_ok(Event const & ev){
std::vector<Particle> const & outgoing = ev.outgoing();
if(ev.decays().size() > 1) // at most one decay
return false;
bool has_AWZH_boson = false;
for( size_t i=0; i<outgoing.size(); ++i ){
auto const & out{ outgoing[i] };
if(is_AWZH_boson(out.type)){
// at most one boson
if(has_AWZH_boson) return false;
has_AWZH_boson = true;
// valid decay for W
if(std::abs(out.type) == ParticleID::Wp){
// exactly 1 decay of W
if( ev.decays().size() != 1 || ev.decays().cbegin()->first != i )
return false;
if( !valid_W_decay(out.type>0?+1:-1, ev.decays().cbegin()->second) )
return false;
}
// valid decay for Z
if(out.type == ParticleID::Z_photon_mix){
// exactly 1 decay
if( ev.decays().size() != 1 || ev.decays().cbegin()->first != i )
return false;
if( !valid_Z_decay(ev.decays().cbegin()->second) )
return false;
}
}
else if(! is_parton(out.type)) return false;
}
return true;
}
/**
* returns all EventTypes implemented in HEJ
*/
size_t implemented_types(std::vector<Particle> const & bosons){
using namespace event_type;
if(bosons.empty()) return FKL | unob | unof | qqbar_exb | qqbar_exf | qqbar_mid;
if(bosons.size()>1) return non_resummable; // multi boson
switch (bosons[0].type) {
case ParticleID::Wp:
case ParticleID::Wm:
return FKL | unob | unof | qqbar_exb | qqbar_exf | qqbar_mid;
case ParticleID::Z_photon_mix:
return FKL | unob | unof;
case ParticleID::h:
return FKL | unob | unof;
default:
return non_resummable;
}
}
/**
* \brief function which determines if type change is consistent with Wp emission.
* @param in incoming Particle id
* @param out outgoing Particle id
* @param is_qqbar Current both incoming/both outgoing?
*
* \see is_Wm_Change
*/
bool is_Wp_Change(ParticleID in, ParticleID out, bool is_qqbar){
using namespace pid;
if(!is_qqbar && (in==d_bar || in==u || in==s_bar || in==c))
return out == (in-1);
if( is_qqbar && (in==d || in==u_bar || in==s || in==c_bar))
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 is_qqbar Current both incoming/both outgoing?
*
* Ensures that change type of quark line is possible by a flavour changing
* Wm emission. Allows checking of is_qqbar currents also.
*/
bool is_Wm_Change(ParticleID in, ParticleID out, bool is_qqbar){
using namespace pid;
if(!is_qqbar && (in==d || in==u_bar || in==s || in==c_bar))
return out == (in+1);
if( is_qqbar && (in==d_bar || in==u || in==s_bar || in==c))
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 is_qqbar Current both incoming/outgoing?
*/
bool no_flavour_change(ParticleID in, ParticleID out, bool is_qqbar){
const int qqbarCurrent = is_qqbar?-1:1;
if(std::abs(in)<=pid::top || in==pid::gluon)
return (in==out*qqbarCurrent);
return false;
}
bool has_enough_jets(Event const & event){
if(event.jets().size() >= 2) return true;
if(event.jets().empty()) return false;
// check for h+jet
const auto the_higgs = std::find_if(
begin(event.outgoing()), end(event.outgoing()),
[](const auto & particle) { return particle.type == pid::higgs; }
);
return the_higgs != end(event.outgoing());
}
bool is_gluon_to_Higgs(const ParticleID in, const ParticleID out) {
return in == pid::gluon && out == pid::Higgs;
}
/**
* \brief check if we have a valid Impact factor
* @param in incoming Particle
* @param out outgoing Particle
* @param is_qqbar Current both incoming/outgoing?
* @param W_change returns +1 if Wp, -1 if Wm, else 0
*/
bool is_valid_impact_factor(
ParticleID in, ParticleID out, bool is_qqbar, int & W_change
){
if( no_flavour_change(in, out, is_qqbar) || is_gluon_to_Higgs(in, out)) {
return true;
}
if( is_Wp_Change(in, out, is_qqbar) ) {
W_change+=1;
return true;
}
if( is_Wm_Change(in, out, is_qqbar) ) {
W_change-=1;
return true;
}
return false;
}
bool is_extremal_higgs_off_quark(
const ParticleID in,
const ParticleID extremal_out,
const ParticleID out
) {
return in == out && extremal_out == pid::higgs && is_anyquark(in);
}
//! Returns all possible classifications from the impact factors
// the beginning points are changed s.t. after the the classification they
// point to the beginning of the (potential) FKL chain
// sets W_change: + if Wp change
// 0 if no change
// - if Wm change
// This function can be used with forward & backwards iterators
template<class OutIterator>
size_t possible_impact_factors(
ParticleID incoming_id, // incoming
OutIterator & begin_out, OutIterator const & end_out, // outgoing
int & W_change, std::vector<Particle> const & boson,
bool const backward // backward?
){
using namespace event_type;
assert(boson.size() < 2);
if(begin_out == end_out) return non_resummable;
// keep track of all states that we don't test
size_t not_tested = qqbar_mid;
if(backward)
not_tested |= unof | qqbar_exf;
else
not_tested |= unob | qqbar_exb;
// Is this LL current?
if( is_valid_impact_factor(incoming_id, begin_out->type, false, W_change) ){
++begin_out;
return not_tested | FKL;
}
// q -> H q and qbar -> H qbar are technically not LL,
// but we treat them as such anyway
if(is_extremal_higgs_off_quark(incoming_id, begin_out->type, std::next(begin_out)->type)) {
std::advance(begin_out, 2);
return not_tested | FKL;
}
// or NLL current?
// -> needs two partons in two different jets
if( std::distance(begin_out, end_out)>=2
){
auto next = std::next(begin_out);
// Is this unordered emisson?
if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){
if( is_valid_impact_factor(
incoming_id, next->type, false, W_change )
){
// veto Higgs inside uno
assert(next!=end_out);
if( !boson.empty() && boson.front().type == ParticleID::h
){
if( (backward && boson.front().rapidity() < next->rapidity())
||(!backward && boson.front().rapidity() > next->rapidity()))
return non_resummable;
}
begin_out = std::next(next);
return not_tested | (backward?unob:unof);
}
}
// Is this QQbar?
else if( incoming_id==pid::gluon ){
if( is_valid_impact_factor(
begin_out->type, next->type, true, W_change )
){
// veto Higgs inside qqbar
assert(next!=end_out);
if( !boson.empty() && boson.front().type == ParticleID::h
){
if( (backward && boson.front().rapidity() < next->rapidity())
||(!backward && boson.front().rapidity() > next->rapidity()))
return non_resummable;
}
begin_out = std::next(next);
return not_tested | (backward?qqbar_exb:qqbar_exf);
}
}
}
return non_resummable;
}
//! Returns all possible classifications from central emissions
// the beginning points are changed s.t. after the the classification they
// point to the end of the emission chain
// sets W_change: + if Wp change
// 0 if no change
// - if Wm change
template<class OutIterator>
size_t possible_central(
OutIterator & begin_out, OutIterator const & end_out,
int & W_change, std::vector<Particle> const & boson
){
using namespace event_type;
assert(boson.size() < 2);
// keep track of all states that we don't test
size_t possible = unob | unof
| qqbar_exb | qqbar_exf;
// Find the first quark or antiquark emission
begin_out = std::find_if(
begin_out, end_out,
[](Particle const & p) { return is_anyquark(p); }
);
// end of chain -> FKL
if( begin_out==end_out ){
return possible | FKL;
}
// is this a qqbar-pair?
// needs two partons in two separate jets
auto next = std::next(begin_out);
if( is_valid_impact_factor(
begin_out->type, next->type, true, W_change )
){
// veto Higgs inside qqbar
if( !boson.empty() && boson.front().type == ParticleID::h
&& boson.front().rapidity() > begin_out->rapidity()
&& boson.front().rapidity() < next->rapidity()
){
return non_resummable;
}
begin_out = std::next(next);
// remaining chain should be pure FKL (gluon or higgs)
if(std::any_of(
begin_out, end_out,
[](Particle const & p) { return is_anyquark(p); }
)) {
return non_resummable;
}
return possible | qqbar_mid;
}
return non_resummable;
}
namespace {
bool is_parton_or_higgs(Particle const & p) {
return is_parton(p) || p.type == pid::higgs;
}
}
/**
* \brief Checks for all event types
* @param ev Event
* @returns Event Type
*
*/
event_type::EventType classify(Event const & ev){
using namespace event_type;
if(! has_enough_jets(ev))
return not_enough_jets;
// currently we can't handle multiple boson states in the ME. So they are
// considered "bad_final_state" even though the "classify" could work with
// them.
if(! final_state_ok(ev))
return bad_final_state;
// initialise variables
auto const & in = ev.incoming();
// range for current checks
auto begin_out = boost::make_filter_iterator(
is_parton_or_higgs, cbegin(ev.outgoing()), cend(ev.outgoing())
);
auto rbegin_out = std::make_reverse_iterator(
boost::make_filter_iterator(
is_parton_or_higgs, cend(ev.outgoing()), cend(ev.outgoing())
)
);
assert(std::distance(begin(in), end(in)) == 2);
assert(std::distance(begin_out, rbegin_out.base()) >= 2);
assert(std::is_sorted(begin_out, rbegin_out.base(), rapidity_less{}));
auto const boson{ filter_AWZH_bosons(ev.outgoing()) };
// we only allow one boson through final_state_ok
assert(boson.size()<=1);
// keep track of potential W couplings, at the end the sum should be 0
int remaining_Wp = 0;
int remaining_Wm = 0;
if(!boson.empty() && std::abs(boson.front().type) == ParticleID::Wp ){
if(boson.front().type>0) ++remaining_Wp;
else ++remaining_Wm;
}
int W_change = 0;
size_t final_type = ~(not_enough_jets | bad_final_state);
// check forward impact factor
final_type &= possible_impact_factors(
in.front().type,
begin_out, rbegin_out.base(),
W_change, boson, true );
if( final_type == non_resummable )
return non_resummable;
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
W_change = 0;
// check backward impact factor
final_type &= possible_impact_factors(
in.back().type,
rbegin_out, std::make_reverse_iterator(begin_out),
W_change, boson, false );
if( final_type == non_resummable )
return non_resummable;
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
W_change = 0;
// check central emissions
final_type &= possible_central(
begin_out, rbegin_out.base(), W_change, boson );
if( final_type == non_resummable )
return non_resummable;
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
// Check whether the right number of Ws are present
if( remaining_Wp != 0 || remaining_Wm != 0 ) return non_resummable;
// result has to be unique
if( (final_type & (final_type-1)) != 0) return non_resummable;
// check that each sub processes is implemented
// (has to be done at the end)
if( (final_type & ~implemented_types(boson)) != 0 )
return non_resummable;
return static_cast<EventType>(final_type);
}
//@}
Particle extract_particle(LHEF::HEPEUP const & hepeup, size_t i){
auto id = static_cast<ParticleID>(hepeup.IDUP[i]);
auto colour = is_parton(id)?hepeup.ICOLUP[i]:optional<Colour>();
return { id,
{ hepeup.PUP[i][0], hepeup.PUP[i][1],
hepeup.PUP[i][2], hepeup.PUP[i][3] },
colour
};
}
bool is_decay_product(std::pair<int, int> const & mothers){
if(mothers.first == 0) return false;
return mothers.second == 0 || mothers.first == mothers.second;
}
} // namespace
Event::EventData::EventData(LHEF::HEPEUP const & hepeup){
parameters.central = EventParameters{
hepeup.scales.mur, hepeup.scales.muf, hepeup.XWGTUP
};
size_t in_idx = 0;
for (int i = 0; i < hepeup.NUP; ++i) {
// skip decay products
// we will add them later on, but we have to ensure that
// the decayed particle is added before
if(is_decay_product(hepeup.MOTHUP[i])) continue;
auto particle = extract_particle(hepeup, i);
// needed to identify mother particles for decay products
particle.p.set_user_index(i+1);
if(hepeup.ISTUP[i] == LHE_Status::in){
if(in_idx > incoming.size()) {
throw std::invalid_argument{
"Event has too many incoming particles"
};
}
incoming[in_idx++] = std::move(particle);
}
else outgoing.emplace_back(std::move(particle));
}
// add decay products
for (int i = 0; i < hepeup.NUP; ++i) {
if(!is_decay_product(hepeup.MOTHUP[i])) continue;
const int mother_id = hepeup.MOTHUP[i].first;
const auto mother = std::find_if(
begin(outgoing), end(outgoing),
[mother_id](Particle const & particle){
return particle.p.user_index() == mother_id;
}
);
if(mother == end(outgoing)){
throw std::invalid_argument{"invalid decay product parent"};
}
const int mother_idx = std::distance(begin(outgoing), mother);
assert(mother_idx >= 0);
decays[mother_idx].emplace_back(extract_particle(hepeup, i));
}
}
Event::Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double const min_jet_pt
):
Event( Event::EventData{
ev.incoming, ev.outgoing, ev.decays,
Parameters<EventParameters>{ev.central, ev.variations}
}.cluster(jet_def, min_jet_pt) )
{}
//! @TODO remove in HEJ 2.2.0
UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup){
Event::EventData const evData{hepeup};
incoming = evData.incoming;
outgoing = evData.outgoing;
decays = evData.decays;
central = evData.parameters.central;
variations = evData.parameters.variations;
}
void Event::EventData::sort(){
// sort particles
std::sort(
begin(incoming), end(incoming),
[](Particle const & o1, Particle const & o2){return o1.p.pz()<o2.p.pz();}
);
auto old_outgoing = std::move(outgoing);
std::vector<size_t> idx(old_outgoing.size());
std::iota(idx.begin(), idx.end(), 0);
std::sort(idx.begin(), idx.end(), [&old_outgoing](size_t i, size_t j){
return old_outgoing[i].rapidity() < old_outgoing[j].rapidity();
});
outgoing.clear();
outgoing.reserve(old_outgoing.size());
for(size_t i: idx) {
outgoing.emplace_back(std::move(old_outgoing[i]));
}
// find decays again
if(!decays.empty()){
auto old_decays = std::move(decays);
decays.clear();
for(size_t i=0; i<idx.size(); ++i) {
auto decay = old_decays.find(idx[i]);
if(decay != old_decays.end())
decays.emplace(i, std::move(decay->second));
}
assert(old_decays.size() == decays.size());
}
}
namespace {
Particle reconstruct_boson(std::vector<Particle> const & leptons) {
Particle decayed_boson;
decayed_boson.p = leptons[0].p + leptons[1].p;
const int pidsum = leptons[0].type + leptons[1].type;
if(pidsum == +1) {
assert(is_antilepton(leptons[0]));
if(is_antineutrino(leptons[0])) {
throw not_implemented{"lepton-flavour violating final state"};
}
assert(is_neutrino(leptons[1]));
// charged antilepton + neutrino means we had a W+
decayed_boson.type = pid::Wp;
}
else if(pidsum == -1) {
assert(is_antilepton(leptons[0]));
if(is_neutrino(leptons[1])) {
throw not_implemented{"lepton-flavour violating final state"};
}
assert(is_antineutrino(leptons[0]));
// charged lepton + antineutrino means we had a W-
decayed_boson.type = pid::Wm;
}
else if(pidsum == 0) {
assert(is_anylepton(leptons[0]));
if(is_anyneutrino(leptons[0])) {
throw not_implemented{"final state with two neutrinos"};
}
// charged lepton-antilepton pair means we had a Z/photon
decayed_boson.type = pid::Z_photon_mix;
}
else {
throw not_implemented{
"final state with leptons "
+ name(leptons[0].type)
+ " and "
+ name(leptons[1].type)
};
}
return decayed_boson;
}
} // namespace
void Event::EventData::reconstruct_intermediate() {
const auto begin_leptons = std::partition(
begin(outgoing), end(outgoing),
[](Particle const & p) {return !is_anylepton(p);}
);
// We can only reconstruct FS with 2 leptons
if(std::distance(begin_leptons, end(outgoing)) != 2) return;
std::vector<Particle> leptons(begin_leptons, end(outgoing));
std::sort(
begin(leptons), end(leptons),
[](Particle const & p0, Particle const & p1) {
assert(is_anylepton(p0) && is_anylepton(p1));
return p0.type < p1.type;
}
);
// `reconstruct_boson` can throw, it should therefore be called before
// changing `outgoing` to allow the user to recover the original EventData
auto boson = reconstruct_boson(leptons);
outgoing.erase(begin_leptons, end(outgoing));
outgoing.emplace_back(std::move(boson));
decays.emplace(outgoing.size()-1, std::move(leptons));
}
Event Event::EventData::cluster(
fastjet::JetDefinition const & jet_def, double const min_jet_pt
){
sort();
return Event{ std::move(incoming), std::move(outgoing), std::move(decays),
std::move(parameters),
jet_def, min_jet_pt
};
}
Event::Event(
std::array<Particle, 2> && incoming,
std::vector<Particle> && outgoing,
std::unordered_map<size_t, std::vector<Particle>> && decays,
Parameters<EventParameters> && parameters,
fastjet::JetDefinition const & jet_def,
double const min_jet_pt
): incoming_{std::move(incoming)},
outgoing_{std::move(outgoing)},
decays_{std::move(decays)},
parameters_{std::move(parameters)},
cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def },
min_jet_pt_{min_jet_pt}
{
jets_ = sorted_by_rapidity(cs_.inclusive_jets(min_jet_pt_));
assert(std::is_sorted(begin(outgoing_), end(outgoing_),
rapidity_less{}));
type_ = classify(*this);
}
namespace {
//! check that Particles have a reasonable colour
bool correct_colour(Particle const & part){
ParticleID id{ part.type };
if(!is_parton(id))
return !part.colour;
if(!part.colour)
return false;
Colour const & col{ *part.colour };
if(is_quark(id))
return col.first != 0 && col.second == 0;
if(is_antiquark(id))
return col.first == 0 && col.second != 0;
assert(id==ParticleID::gluon);
return col.first != 0 && col.second != 0 && col.first != col.second;
}
//! Connect parton to t-channel colour line & update the line
//! returns false if connection not possible
template<class OutIterator>
bool try_connect_t(OutIterator const & it_part, Colour & line_colour){
if( line_colour.first == it_part->colour->second ){
line_colour.first = it_part->colour->first;
return true;
}
if( line_colour.second == it_part->colour->first ){
line_colour.second = it_part->colour->second;
return true;
}
return false;
}
//! Connect parton to u-channel colour line & update the line
//! returns false if connection not possible
template<class OutIterator>
bool try_connect_u(OutIterator & it_part, Colour & line_colour){
auto it_next = std::next(it_part);
if( try_connect_t(it_next, line_colour)
&& try_connect_t(it_part, line_colour)
){
it_part=it_next;
return true;
}
return false;
}
} // namespace
bool Event::is_leading_colour() const {
if( !correct_colour(incoming()[0]) || !correct_colour(incoming()[1]) )
return false;
Colour line_colour = *incoming()[0].colour;
std::swap(line_colour.first, line_colour.second);
// reasonable colour
if(!std::all_of(outgoing().cbegin(), outgoing().cend(), correct_colour))
return false;
for(auto it_part = cbegin_partons(); it_part!=cend_partons(); ++it_part){
switch (type()) {
case event_type::FKL:
if( !try_connect_t(it_part, line_colour) )
return false;
break;
case event_type::unob:
case event_type::qqbar_exb: {
if( !try_connect_t(it_part, line_colour)
// u-channel only allowed at impact factor
&& (std::distance(cbegin_partons(), it_part)!=0
|| !try_connect_u(it_part, line_colour)))
return false;
break;
}
case event_type::unof:
case event_type::qqbar_exf: {
if( !try_connect_t(it_part, line_colour)
// u-channel only allowed at impact factor
&& (std::distance(it_part, cend_partons())!=2
|| !try_connect_u(it_part, line_colour)))
return false;
break;
}
case event_type::qqbar_mid:{
auto it_next = std::next(it_part);
if( !try_connect_t(it_part, line_colour)
// u-channel only allowed at q-qbar/qbar-q pair
&& ( ( !(is_quark(*it_part) && is_antiquark(*it_next))
&& !(is_antiquark(*it_part) && is_quark(*it_next)))
|| !try_connect_u(it_part, line_colour))
)
return false;
break;
}
default:
throw std::logic_error{"unreachable"};
}
// no colour singlet exchange/disconnected diagram
if(line_colour.first == line_colour.second)
return false;
}
return (incoming()[1].colour->first == line_colour.first)
&& (incoming()[1].colour->second == line_colour.second);
}
namespace {
//! connect incoming Particle to colour flow
void connect_incoming(Particle & in, int & colour, int & anti_colour){
in.colour = std::make_pair(anti_colour, colour);
// gluon
if(in.type == pid::gluon)
return;
if(in.type > 0){
// quark
assert(is_quark(in));
in.colour->second = 0;
colour*=-1;
return;
}
// anti-quark
assert(is_antiquark(in));
in.colour->first = 0;
anti_colour*=-1;
}
//! connect outgoing Particle to t-channel colour flow
template<class OutIterator>
void connect_tchannel(
OutIterator & it_part, int & colour, int & anti_colour, RNG & ran
){
assert(colour>0 || anti_colour>0);
if(it_part->type == ParticleID::gluon){
// gluon
if(colour>0 && anti_colour>0){
// on g line => connect to colour OR anti-colour (random)
if(ran.flat() < 0.5){
it_part->colour = std::make_pair(colour+2,colour);
colour+=2;
} else {
it_part->colour = std::make_pair(anti_colour, anti_colour+2);
anti_colour+=2;
}
} else if(colour > 0){
// on q line => connect to available colour
it_part->colour = std::make_pair(colour+2, colour);
colour+=2;
} else {
assert(colour<0 && anti_colour>0);
// on qbar line => connect to available anti-colour
it_part->colour = std::make_pair(anti_colour, anti_colour+2);
anti_colour+=2;
}
} else if(is_quark(*it_part)) {
// quark
assert(anti_colour>0);
if(colour>0){
// on g line => connect and remove anti-colour
it_part->colour = std::make_pair(anti_colour, 0);
anti_colour+=2;
anti_colour*=-1;
} else {
// on qbar line => new colour
colour*=-1;
it_part->colour = std::make_pair(colour, 0);
}
} else if(is_antiquark(*it_part)) {
// anti-quark
assert(colour>0);
if(anti_colour>0){
// on g line => connect and remove colour
it_part->colour = std::make_pair(0, colour);
colour+=2;
colour*=-1;
} else {
// on q line => new anti-colour
anti_colour*=-1;
it_part->colour = std::make_pair(0, anti_colour);
}
} else { // not a parton
assert(!is_parton(*it_part));
it_part->colour = {};
}
}
//! connect to t- or u-channel colour flow
template<class OutIterator>
void connect_utchannel(
OutIterator & it_part, int & colour, int & anti_colour, RNG & ran
){
OutIterator it_first = it_part++;
if(ran.flat()<.5) {// t-channel
connect_tchannel(it_first, colour, anti_colour, ran);
connect_tchannel(it_part, colour, anti_colour, ran);
}
else { // u-channel
connect_tchannel(it_part, colour, anti_colour, ran);
connect_tchannel(it_first, colour, anti_colour, ran);
}
}
} // namespace
bool Event::generate_colours(RNG & ran){
// generate only for HEJ events
if(!event_type::is_resummable(type()))
return false;
assert(std::is_sorted(
begin(outgoing()), end(outgoing()), rapidity_less{}));
assert(incoming()[0].pz() < incoming()[1].pz());
// positive (anti-)colour -> can connect
// negative (anti-)colour -> not available/used up by (anti-)quark
int colour = COLOUR_OFFSET;
int anti_colour = colour+1;
// initialise first
connect_incoming(incoming_[0], colour, anti_colour);
// reset outgoing colours
std::for_each(outgoing_.begin(), outgoing_.end(),
[](Particle & part){ part.colour = {};});
for(auto it_part = begin_partons(); it_part!=end_partons(); ++it_part){
switch (type()) {
// subleading can connect to t- or u-channel
case event_type::unob:
case event_type::qqbar_exb: {
if( std::distance(begin_partons(), it_part)==0)
connect_utchannel(it_part, colour, anti_colour, ran);
else
connect_tchannel(it_part, colour, anti_colour, ran);
break;
}
case event_type::unof:
case event_type::qqbar_exf: {
if( std::distance(it_part, end_partons())==2)
connect_utchannel(it_part, colour, anti_colour, ran);
else
connect_tchannel(it_part, colour, anti_colour, ran);
break;
}
case event_type::qqbar_mid:{
auto it_next = std::next(it_part);
if( std::distance(begin_partons(), it_part)>0
&& std::distance(it_part, end_partons())>2
&& ( (is_quark(*it_part) && is_antiquark(*it_next))
|| (is_antiquark(*it_part) && is_quark(*it_next)) )
)
connect_utchannel(it_part, colour, anti_colour, ran);
else
connect_tchannel(it_part, colour, anti_colour, ran);
break;
}
default: // rest has to be t-channel
connect_tchannel(it_part, colour, anti_colour, ran);
}
}
// Connect last
connect_incoming(incoming_[1], anti_colour, colour);
assert(is_leading_colour());
return true;
} // generate_colours
namespace {
bool valid_parton(
std::vector<fastjet::PseudoJet> const & jets,
Particle const & parton, int const idx,
double const soft_pt_regulator, double const min_extparton_pt
){
// TODO code overlap with PhaseSpacePoint::pass_extremal_cuts
if(min_extparton_pt > parton.pt()) return false;
if(idx<0) return false;
assert(static_cast<int>(jets.size())>=idx);
auto const & jet{ jets[idx] };
return (parton.p - jet).pt()/jet.pt() <= soft_pt_regulator;
}
} // namespace
// this should work with multiple types
bool Event::valid_hej_state(double const soft_pt_regulator,
double const min_pt
) const {
using namespace event_type;
if(!is_resummable(type()))
return false;
auto const & jet_idx{ particle_jet_indices() };
auto idx_begin{ jet_idx.cbegin() };
auto idx_end{ jet_idx.crbegin() };
auto part_begin{ cbegin_partons() };
auto part_end{ crbegin_partons() };
// always seperate extremal jets
if(!is_backward_g_to_h(*this)) {
if(! valid_parton(jets(), *part_begin, *idx_begin, soft_pt_regulator, min_pt)) {
return false;
}
++part_begin;
++idx_begin;
// unob -> second parton in own jet
if( type() & (unob | qqbar_exb) ){
if( !valid_parton(jets(), *part_begin, *idx_begin,
soft_pt_regulator, min_pt) )
return false;
++part_begin;
++idx_begin;
}
}
if(!is_forward_g_to_h(*this)) {
if(!valid_parton(jets(), *part_end, *idx_end, soft_pt_regulator, min_pt)) {
return false;
}
++part_end;
++idx_end;
if( type() & (unof | qqbar_exf) ){
if( !valid_parton(jets(), *part_end, *idx_end,
soft_pt_regulator, min_pt) )
return false;
++part_end;
// ++idx_end; // last check, we don't need idx_end afterwards
}
}
if( type() & qqbar_mid ){
// find qqbar pair
auto begin_qqbar{ std::find_if( part_begin, part_end.base(),
[](Particle const & part) -> bool {
return part.type != ParticleID::gluon;
}
)};
assert(begin_qqbar != part_end.base());
long int qqbar_pos{ std::distance(part_begin, begin_qqbar) };
assert(qqbar_pos >= 0);
idx_begin+=qqbar_pos;
if( !( valid_parton(jets(), *begin_qqbar, *idx_begin,
soft_pt_regulator, min_pt)
&& valid_parton(jets(), *std::next(begin_qqbar), *std::next(idx_begin),
soft_pt_regulator, min_pt)
))
return false;
}
return true;
}
+ bool Event::valid_incoming() const{
+ for(std::size_t i=0; i < incoming_.size(); ++i){
+ if(!(HEJ::nearby_ep(std::abs(incoming_[i].pz()), incoming_[i].E(), TOL*incoming_[i].E())
+ && (incoming_[i].pt()==0.)))
+ return false;
+ }
+ return true;
+ }
+
Event::ConstPartonIterator Event::begin_partons() const {
return cbegin_partons();
}
Event::ConstPartonIterator Event::cbegin_partons() const {
return {HEJ::is_parton, cbegin(outgoing()), cend(outgoing())};
}
Event::ConstPartonIterator Event::end_partons() const {
return cend_partons();
}
Event::ConstPartonIterator Event::cend_partons() const {
return {HEJ::is_parton, cend(outgoing()), cend(outgoing())};
}
Event::ConstReversePartonIterator Event::rbegin_partons() const {
return crbegin_partons();
}
Event::ConstReversePartonIterator Event::crbegin_partons() const {
return std::reverse_iterator<ConstPartonIterator>( cend_partons() );
}
Event::ConstReversePartonIterator Event::rend_partons() const {
return crend_partons();
}
Event::ConstReversePartonIterator Event::crend_partons() const {
return std::reverse_iterator<ConstPartonIterator>( cbegin_partons() );
}
Event::PartonIterator Event::begin_partons() {
return {HEJ::is_parton, begin(outgoing_), end(outgoing_)};
}
Event::PartonIterator Event::end_partons() {
return {HEJ::is_parton, end(outgoing_), end(outgoing_)};
}
Event::ReversePartonIterator Event::rbegin_partons() {
return std::reverse_iterator<PartonIterator>( end_partons() );
}
Event::ReversePartonIterator Event::rend_partons() {
return std::reverse_iterator<PartonIterator>( begin_partons() );
}
namespace {
void print_momentum(std::ostream & os, fastjet::PseudoJet const & part){
constexpr int prec = 6;
const std::streamsize orig_prec = os.precision();
os <<std::scientific<<std::setprecision(prec) << "["
<<std::setw(2*prec+1)<<std::right<< part.px() << ", "
<<std::setw(2*prec+1)<<std::right<< part.py() << ", "
<<std::setw(2*prec+1)<<std::right<< part.pz() << ", "
<<std::setw(2*prec+1)<<std::right<< part.E() << "]"<< std::fixed;
os.precision(orig_prec);
}
void print_colour(std::ostream & os, optional<Colour> const & col){
constexpr int width = 3;
if(!col)
os << "(no color)"; // American spelling for better alignment
else
os << "(" <<std::setw(width)<<std::right<< col->first
<< ", " <<std::setw(width)<<std::right<< col->second << ")";
}
} // namespace
std::ostream& operator<<(std::ostream & os, Event const & ev){
constexpr int prec = 4;
constexpr int wtype = 3; // width for types
const std::streamsize orig_prec = os.precision();
os <<std::setprecision(prec)<<std::fixed;
os << "########## " << name(ev.type()) << " ##########" << std::endl;
os << "Incoming particles:\n";
for(auto const & in: ev.incoming()){
os <<std::setw(wtype)<< in.type << ": ";
print_colour(os, in.colour);
os << " ";
print_momentum(os, in.p);
os << std::endl;
}
os << "\nOutgoing particles: " << ev.outgoing().size() << "\n";
for(auto const & out: ev.outgoing()){
os <<std::setw(wtype)<< out.type << ": ";
print_colour(os, out.colour);
os << " ";
print_momentum(os, out.p);
os << " => rapidity="
<<std::setw(2*prec-1)<<std::right<< out.rapidity() << std::endl;
}
os << "\nForming Jets: " << ev.jets().size() << "\n";
for(auto const & jet: ev.jets()){
print_momentum(os, jet);
os << " => rapidity="
<<std::setw(2*prec-1)<<std::right<< jet.rapidity() << std::endl;
}
if(!ev.decays().empty() ){
os << "\nDecays: " << ev.decays().size() << "\n";
for(auto const & decay: ev.decays()){
os <<std::setw(wtype)<< ev.outgoing()[decay.first].type
<< " (" << decay.first << ") to:\n";
for(auto const & out: decay.second){
os <<" "<<std::setw(wtype)<< out.type << ": ";
print_momentum(os, out.p);
os << " => rapidity="
<<std::setw(2*prec-1)<<std::right<< out.rapidity() << std::endl;
}
}
}
os << std::defaultfloat;
os.precision(orig_prec);
return os;
}
double shat(Event const & ev){
return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
}
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){
LHEF::HEPEUP result;
result.heprup = heprup;
result.weights = {{event.central().weight, nullptr}};
for(auto const & var: event.variations()){
result.weights.emplace_back(var.weight, nullptr);
}
size_t num_particles = event.incoming().size() + event.outgoing().size();
for(auto const & decay: event.decays()) num_particles += decay.second.size();
result.NUP = num_particles;
// the following entries are pretty much meaningless
result.IDPRUP = event.type(); // event type
result.AQEDUP = 1./128.; // alpha_EW
//result.AQCDUP = 0.118 // alpha_QCD
// end meaningless part
result.XWGTUP = event.central().weight;
result.SCALUP = event.central().muf;
result.scales.muf = event.central().muf;
result.scales.mur = event.central().mur;
result.scales.SCALUP = event.central().muf;
result.pdfinfo.p1 = event.incoming().front().type;
result.pdfinfo.p2 = event.incoming().back().type;
result.pdfinfo.scale = event.central().muf;
result.IDUP.reserve(num_particles); // PID
result.ISTUP.reserve(num_particles); // status (in, out, decay)
result.PUP.reserve(num_particles); // momentum
result.MOTHUP.reserve(num_particles); // index mother particle
result.ICOLUP.reserve(num_particles); // colour
// incoming
std::array<Particle, 2> incoming{ event.incoming() };
// First incoming should be positive pz according to LHE standard
// (or at least most (everyone?) do it this way, and Pythia assumes it)
if(incoming[0].pz() < incoming[1].pz())
std::swap(incoming[0], incoming[1]);
for(Particle const & in: incoming){
result.IDUP.emplace_back(in.type);
result.ISTUP.emplace_back(LHE_Status::in);
result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()});
result.MOTHUP.emplace_back(0, 0);
assert(in.colour);
result.ICOLUP.emplace_back(*in.colour);
}
// outgoing
for(size_t i = 0; i < event.outgoing().size(); ++i){
Particle const & out = event.outgoing()[i];
result.IDUP.emplace_back(out.type);
const int status = event.decays().count(i) != 0u
?LHE_Status::decay
:LHE_Status::out;
result.ISTUP.emplace_back(status);
result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
result.MOTHUP.emplace_back(1, 2);
if(out.colour)
result.ICOLUP.emplace_back(*out.colour);
else{
result.ICOLUP.emplace_back(std::make_pair(0,0));
}
}
// decays
for(auto const & decay: event.decays()){
for(auto const & out: decay.second){
result.IDUP.emplace_back(out.type);
result.ISTUP.emplace_back(LHE_Status::out);
result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
const size_t mother_idx = 1 + event.incoming().size() + decay.first;
result.MOTHUP.emplace_back(mother_idx, mother_idx);
result.ICOLUP.emplace_back(0,0);
}
}
assert(result.ICOLUP.size() == num_particles);
static constexpr double unknown_spin = 9.; //per Les Houches accord
result.VTIMUP = std::vector<double>(num_particles, unknown_spin);
result.SPINUP = result.VTIMUP;
return result;
}
} // namespace HEJ
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 92dfaaf..8d76a0c 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,2164 +1,2171 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <iterator>
#include <limits>
#include <unordered_map>
#include <utility>
#include "fastjet/PseudoJet.hh"
#include "HEJ/ConfigFlags.hh"
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/Event.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/Hjets.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Particle.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Zjets.hh"
#include "HEJ/event_types.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/jets.hh"
#include "HEJ/utility.hh"
namespace HEJ {
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const {
const double lambda = param_.regulator_lambda;
const double result = - alpha_s*N_C/M_PI*std::log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
return (
1. + alpha_s/(4.*M_PI)*BETA0*std::log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(Event const & event) const {
std::vector <double> tree_kin_part=tree_kin(event);
std::vector <Weights> virtual_part=virtual_corrections(event);
if(tree_kin_part.size() != virtual_part.size()) {
throw std::logic_error("tree and virtuals have different sizes");
}
Weights sum = Weights{0., std::vector<double>(event.variations().size(), 0.)};
for(size_t i=0; i<tree_kin_part.size(); ++i) {
sum += tree_kin_part.at(i)*virtual_part.at(i);
}
return tree_param(event)*sum;
}
Weights MatrixElement::tree(Event const & event) const {
std::vector <double> tree_kin_part=tree_kin(event);
double sum = 0.;
for(double i : tree_kin_part) {
sum += i;
}
return tree_param(event)*sum;
}
Weights MatrixElement::tree_param(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
std::vector<Weights> MatrixElement::virtual_corrections(Event const & event) const {
if(! is_resummable(event.type())) {
return {Weights{0., std::vector<double>(event.variations().size(), 0.)}};
}
// only compute once for each renormalisation scale
std::unordered_map<double, std::vector<double> > known_vec;
std::vector<double> central_vec=virtual_corrections(event, event.central().mur);
known_vec.emplace(event.central().mur, central_vec);
for(auto const & var: event.variations()) {
const auto ME_it = known_vec.find(var.mur);
if(ME_it == end(known_vec)) {
known_vec.emplace(var.mur, virtual_corrections(event, var.mur));
}
}
// At this stage known_vec contains one vector of virtual corrections for each mur value
// Now put this into a vector of Weights
std::vector<Weights> result_vec;
for(size_t i=0; i<central_vec.size(); ++i) {
Weights result;
result.central = central_vec.at(i);
for(auto const & var: event.variations()) {
const auto ME_it = known_vec.find(var.mur);
result.variations.emplace_back(ME_it->second.at(i));
}
result_vec.emplace_back(result);
}
return result_vec;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
const double mur,
Particle const & WBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - partons[0].p;
std::size_t first_idx = 0;
std::size_t last_idx = partons.size() - 1;
#ifndef NDEBUG
bool wc = true;
#endif
bool wqq = false;
// With extremal qqbar or unordered gluon outside the extremal
// partons then it is not part of the FKL ladder and does not
// contribute to the virtual corrections. W emitted from the
// most backward leg must be taken into account in t-channel
if (event.type() == event_type::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else if (event.type() == event_type::qqbar_exb) {
q -= partons[1].p;
++first_idx;
if (std::abs(partons[0].type) != std::abs(partons[1].type)){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else {
if(event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf){
--last_idx;
}
if (in[0].type != partons[0].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
std::size_t first_idx_qqbar = last_idx;
std::size_t last_idx_qqbar = last_idx;
//if qqbarMid event, virtual correction do not occur between qqbar pair.
if(event.type() == event_type::qqbar_mid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(std::abs(backquark->type) != std::abs((backquark+1)->type)) {
wqq=true;
#ifndef NDEBUG
wc=false;
#endif
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqbar = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(std::size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -=partons[j+1].p;
} // End Loop one
if (last_idx != first_idx_qqbar) q -= partons[last_idx+1].p;
if (wqq) q -= WBoson.p;
for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -= partons[j+1].p;
}
#ifndef NDEBUG
if (wc) q -= WBoson.p;
assert(
nearby(q, -1*pb, norm)
|| is_AWZH_boson(partons.back().type)
|| event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf
);
#endif
return std::exp(exponent);
}
std::vector <double> MatrixElement::virtual_corrections_Z_qq(
Event const & event,
const double mur,
Particle const & ZBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q_t = pa - partons[0].p - ZBoson.p;
fastjet::PseudoJet q_b = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q_t -= partons[1].p;
q_b -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--last_idx;
}
double sum_top=0.;
double sum_bot=0.;
double sum_mix=0.;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
const double dy = partons[j+1].rapidity() - partons[j].rapidity();
const double tmp_top = omega0(alpha_s, mur, q_t)*dy;
const double tmp_bot = omega0(alpha_s, mur, q_b)*dy;
sum_top += tmp_top;
sum_bot += tmp_bot;
sum_mix += (tmp_top + tmp_bot) / 2.;
q_t -= partons[j+1].p;
q_b -= partons[j+1].p;
}
return {exp(sum_top), exp(sum_bot), exp(sum_mix)};
}
double MatrixElement::virtual_corrections_Z_qg(
Event const & event,
const double mur,
Particle const & ZBoson,
const bool is_gq_event
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
// If this is a gq event, don't subtract the Z momentum from first q
fastjet::PseudoJet q = (is_gq_event ? pa - partons[0].p : pa - partons[0].p - ZBoson.p);
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--last_idx;
}
double sum=0.;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
sum += omega0(alpha_s, mur, q)*(partons[j+1].rapidity()
- partons[j].rapidity());
q -= partons[j+1].p;
}
return exp(sum);
}
std::vector<double> MatrixElement::virtual_corrections(
Event const & event,
const double mur
) const{
auto const & in = event.incoming();
auto const & out = event.outgoing();
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
const auto AWZH_boson = std::find_if(
begin(out), end(out),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson != end(out) && std::abs(AWZH_boson->type) == pid::Wp){
return {virtual_corrections_W(event, mur, *AWZH_boson)};
}
if(AWZH_boson != end(out) && AWZH_boson->type == pid::Z_photon_mix){
if(is_gluon(in.back().type)){
// This is a qg event
return {virtual_corrections_Z_qg(event, mur, *AWZH_boson, false)};
}
if(is_gluon(in.front().type)){
// This is a gq event
return {virtual_corrections_Z_qg(event, mur, *AWZH_boson, true)};
}
// This is a qq event
return virtual_corrections_Z_qq(event, mur, *AWZH_boson);
}
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
std::size_t first_idx = 0;
std::size_t last_idx = out.size() - 1;
// if there is a Higgs boson _not_ emitted off an incoming gluon,
// extremal qqbar 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 && in.front().type != pid::gluon)
|| event.type() == event_type::unob
|| event.type() == event_type::qqbar_exb){
q -= out[1].p;
++first_idx;
}
if((out.back().type == pid::Higgs && in.back().type != pid::gluon)
|| event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf){
--last_idx;
}
std::size_t first_idx_qqbar = last_idx;
std::size_t last_idx_qqbar = last_idx;
//if central qqbar event, virtual correction do not occur between q-qbar.
if(event.type() == event_type::qqbar_mid){
const auto backquark = std::find_if(
begin(out) + 1, end(out) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
);
if(backquark == end(out) || (backquark+1)->type==pid::gluon) return {0.};
last_idx = std::distance(begin(out), backquark);
first_idx_qqbar = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(std::size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
if (last_idx != first_idx_qqbar) q -= out[last_idx+1].p;
for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| (out.back().type == pid::Higgs && in.back().type != pid::gluon)
|| event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf
);
return {std::exp(exponent)};
}
namespace {
//! Lipatov vertex for partons emitted into extremal jets
CLHEP::HepLorentzVector CLipatov(
CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
const CLHEP::HepLorentzVector p5 = qav-qbv;
const CLHEP::HepLorentzVector CL = -(qav+qbv)
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
return CL;
}
double C2Lipatov(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2
){
const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, p1, p2);
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
const double kperp=(qav-qbv).perp();
if (kperp>lambda)
return Cls;
return Cls-4./(kperp*kperp);
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, p1, p2);
const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, p1, p2);
return -CL_t.dot(CL_b);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, p1, p2)
/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
const double kperp = (qav_t - qbv_t).perp();
if (kperp > lambda){
return Cls;
}
return Cls - 4.0 / (kperp * kperp);
}
CLHEP::HepLorentzVector CLipatov(
CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
){
const CLHEP::HepLorentzVector p5 = qav-qbv;
const CLHEP::HepLorentzVector CL = -(qav+qbv)
+ 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;
}
//! Lipatov vertex
double C2Lipatov( // B
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim,
CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom,
CLHEP::HepLorentzVector const & pop
){
const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, pim, pip, pom, pop);
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
const double kperp=(qav-qbv).perp();
if (kperp>lambda)
return Cls;
return Cls-4./(kperp*kperp);
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
) {
const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, pim, pip, pom, pop);
const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, pim, pip, pom, pop);
return -CL_t.dot(CL_b);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, pa, pb, p1, p2)
/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
const double kperp = (qav_t - qbv_t).perp();
if (kperp > lambda) {
return Cls;
}
return Cls - 4.0 / (kperp * kperp);
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pg Unordered gluon momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_uno_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
using namespace currents;
assert(aptype!=pid::gluon); // aptype cannot be gluon
if (bptype==pid::gluon) {
if (is_quark(aptype))
return ME_unob_qg(pg,p1,pa,pn,pb);
return ME_unob_qbarg(pg,p1,pa,pn,pb);
}
if (is_antiquark(bptype)) {
if (is_quark(aptype))
return ME_unob_qQbar(pg,p1,pa,pn,pb);
return ME_unob_qbarQbar(pg,p1,pa,pn,pb);
}
//bptype == quark
if (is_quark(aptype))
return ME_unob_qQ(pg,p1,pa,pn,pb);
return ME_unob_qbarQ(pg,p1,pa,pn,pb);
}
/** Matrix element squared for tree-level current-current scattering
* @param bptype Particle b PDG ID
* @param pgin Incoming gluon momentum
* @param pq Quark from splitting Momentum
* @param pqbar Anti-quark from splitting Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param swap_qqbar Ordering of qqbar pair. False: pqbar extremal.
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The forward qqbar contribution can be calculated by reversing the
* argument ordering.
*/
double ME_qqbar_current(
ParticleID bptype,
CLHEP::HepLorentzVector const & pgin,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
bool const swap_qqbar
){
using namespace currents;
if (bptype==pid::gluon) {
if (swap_qqbar) // pq extremal
return ME_Exqqbar_qqbarg(pgin,pq,pqbar,pn,pb);
// pqbar extremal
return ME_Exqqbar_qbarqg(pgin,pq,pqbar,pn,pb);
}
// b leg quark line
if (swap_qqbar) //extremal pq
return ME_Exqqbar_qqbarQ(pgin,pq,pqbar,pn,pb);
return ME_Exqqbar_qbarqQ(pgin,pq,pqbar,pn,pb);
}
/* \brief Matrix element squared for central qqbar tree-level current-current
* scattering
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqbarpair
* @param nbelow Number of gluons emitted after central qqbarpair
* @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
* @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering
*/
double ME_qqbar_mid_current(
ParticleID aptype, ParticleID bptype, int nabove,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<CLHEP::HepLorentzVector> const & partons
){
using namespace currents;
// CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards)
const bool swap_qqbar=pqbar.rapidity() < pq.rapidity();
double wt=1.;
if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F;
return wt*ME_Cenqqbar_qq(pa, pb, partons, is_antiquark(bptype),
is_antiquark(aptype), swap_qqbar, nabove);
}
/** 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(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
using namespace currents;
if (aptype==pid::gluon && bptype==pid::gluon) {
return ME_gg(pn,pb,p1,pa);
}
if (aptype==pid::gluon && bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_qg(pn,pb,p1,pa);
return ME_qbarg(pn,pb,p1,pa);
}
if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_qg(p1,pa,pn,pb);
return ME_qbarg(p1,pa,pn,pb);
}
// they are both quark
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_qQ(pn,pb,p1,pa);
return ME_qQbar(pn,pb,p1,pa);
}
if (is_quark(aptype))
return ME_qQbar(p1,pa,pn,pb);
return ME_qbarQbar(pn,pb,p1,pa);
}
/** Matrix element squared for tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_W_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc, ParticleProperties const & Wprop
){
using namespace currents;
// We know it cannot be gg incoming.
assert(!(aptype==pid::gluon && bptype==pid::gluon));
if (aptype==pid::gluon && bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_W_qg(pn,plbar,pl,pb,p1,pa,Wprop);
return ME_W_qbarg(pn,plbar,pl,pb,p1,pa,Wprop);
}
if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_W_qg(p1,plbar,pl,pa,pn,pb,Wprop);
return ME_W_qbarg(p1,plbar,pl,pa,pn,pb,Wprop);
}
// they are both quark
if (wc){ // emission off b, (first argument pbout)
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_W_qQ(pn,plbar,pl,pb,p1,pa,Wprop);
return ME_W_qQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
if (is_quark(aptype))
return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa,Wprop);
return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
// emission off a, (first argument paout)
if (is_quark(aptype)) {
if (is_quark(bptype))
return ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop);
return ME_W_qQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
// a is anti-quark
if (is_quark(bptype))
return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb,Wprop);
return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_W_uno_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc, ParticleProperties const & Wprop
){
using namespace currents;
// we know they are not both gluons
assert(bptype != pid::gluon || aptype != pid::gluon);
if (bptype == pid::gluon && aptype != pid::gluon) {
// b gluon => W emission off a
if (is_quark(aptype))
return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
// they are both quark
if (wc) {// emission off b, i.e. b is first current
if (is_quark(bptype)){
if (is_quark(aptype))
return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
if (is_quark(aptype))
return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
// wc == false, emission off a, i.e. a is first current
if (is_quark(aptype)) {
if (is_quark(bptype)) //qq
return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
//qqbar
return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
// a is anti-quark
if (is_quark(bptype)) //qbarq
return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
//qbarqbar
return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
/** \brief Matrix element squared for backward qqbar tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param swap_qqbar Ordering of qqbar pair. False: pqbar extremal.
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqbarb Tree-Level Current-Current Scattering
*
* @note calculate forwards qqbar contribution by reversing argument ordering.
*/
double ME_W_qqbar_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const swap_qqbar, bool const wc,
ParticleProperties const & Wprop
){
using namespace currents;
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==pid::gluon && bptype==pid::gluon) {
//a gluon, b gluon gg->qqbarWg
// This will be a wqqbar emission as there is no other possible W Emission
// Site.
if (swap_qqbar)
return ME_WExqqbar_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb, Wprop);
return ME_WExqqbar_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop);
}
assert(aptype==pid::gluon && bptype!=pid::gluon );
//a gluon => W emission off b leg or qqbar
if (!wc){ // W Emitted from backwards qqbar
if (swap_qqbar)
return ME_WExqqbar_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop);
return ME_WExqqbar_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb, Wprop);
}
// W Must be emitted from forwards leg.
return ME_W_Exqqbar_QQq(pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop);
}
/* \brief Matrix element squared for central qqbar 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 qqbarpair
* @param nbelow Number of gluons emitted after central qqbarpair
* @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 qqbar
* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
* @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering
*/
double ME_W_qqbar_mid_current(
ParticleID aptype, ParticleID bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<CLHEP::HepLorentzVector> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc,
ParticleProperties const & Wprop
){
using namespace currents;
// CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards)
const bool swap_qqbar=pqbar.rapidity() < pq.rapidity();
double wt=1.;
if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F;
if(wqq)
return wt*ME_WCenqqbar_qq(pa, pb, pl, plbar, partons,
is_antiquark(aptype),is_antiquark(bptype),
swap_qqbar, nabove, Wprop);
return wt*ME_W_Cenqqbar_qq(pa, pb, pl, plbar, partons,
is_antiquark(aptype), is_antiquark(bptype),
swap_qqbar, nabove, nbelow, wc, Wprop);
}
/** Matrix element squared for tree-level current-current scattering With Z+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 plbar Final state positron momentum
* @param pl Final state electron momentum
* @param Zprop Z properties
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
std::vector<double> ME_Z_current(
const ParticleID aptype, const ParticleID 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,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
if(is_anyquark(aptype) && is_gluon(bptype)){
// This is a qg event
return { ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw) };
}
if(is_gluon(aptype) && is_anyquark(bptype)){
// This is a gq event
return { ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw) };
}
assert(is_anyquark(aptype) && is_anyquark(bptype));
// This is a qq event
return ME_Z_qQ(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With Z+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param plbar Final state positron momentum
* @param pl Final state electron momentum
* @param Zprop Z properties
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
std::vector<double> ME_Z_uno_current(
const ParticleID aptype, const ParticleID 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,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
if (is_anyquark(aptype) && is_gluon(bptype)) {
// This is a qg event
return { ME_Zuno_qg(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw) };
}
if (is_gluon(aptype) && is_anyquark(bptype)) {
// This is a gq event
return { ME_Zuno_qg(pb,pa,pg,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw) };
}
assert(is_anyquark(aptype) && is_anyquark(bptype));
// This is a qq event
return ME_Zuno_qQ(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
}
/** \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(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb, double vev
){
using namespace currents;
if (aptype==pid::gluon && bptype==pid::gluon)
// gg initial state
return ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev);
if (aptype==pid::gluon&&bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_H_qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
return ME_H_qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
}
if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_H_qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
return ME_H_qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
}
// they are both quark
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
return ME_H_qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
if (is_quark(aptype))
return ME_H_qbarQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
return ME_H_qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param pg Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*
* @note This function assumes unordered gluon backwards from pa-p1 current.
* For unof, reverse call order
*/
double ME_Higgs_current_uno(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb, double vev
){
using namespace currents;
if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_H_unob_gQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
return ME_H_unob_gQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
// they are both quark
if (is_quark(aptype)) {
if (is_quark(bptype))
return ME_H_unob_qQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
return ME_H_unob_qbarQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
if (is_quark(bptype))
return ME_H_unob_qQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
return ME_H_unob_qbarQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
void validate(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
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
std::vector<double> MatrixElement::tree_kin(
Event const & ev
) const {
if(!ev.valid_hej_state(param_.soft_pt_regulator)) return {0.};
+ if(!ev.valid_incoming()){
+ throw std::invalid_argument{
+ "Invalid momentum for one or more incoming particles. "
+ "Incoming momenta must have vanishing mass and transverse momentum."
+ };
+ }
+
auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing()))
return {tree_kin_jets(ev)};
switch(AWZH_boson->type){
case pid::Higgs:
return {tree_kin_Higgs(ev)};
case pid::Wp:
case pid::Wm:
return {tree_kin_W(ev)};
case pid::Z_photon_mix:
return tree_kin_Z(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
namespace {
constexpr int EXTREMAL_JET_IDX = 1;
constexpr int NO_EXTREMAL_JET_IDX = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == EXTREMAL_JET_IDX;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
double lambda
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
template<class InputIterator>
std::vector <double> FKL_ladder_weight_mix(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0_t, CLHEP::HepLorentzVector const & q0_b,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
const double lambda
){
double wt_top = 1;
double wt_bot = 1;
double wt_mix = 1;
auto qi_t = q0_t;
auto qi_b = q0_b;
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_t = qi_t - g;
const auto qip1_b = qi_b - g;
if(treat_as_extremal(*gluon_it)){
wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, lambda)*C_A;
wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, lambda)*C_A;
wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, lambda)*C_A;
} else{
wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, p1, pn, lambda)*C_A;
wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A;
wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A;
}
qi_t = qip1_t;
qi_b = qip1_b;
}
return {wt_top, wt_bot, wt_mix};
}
std::vector<Particle> tag_extremal_jet_partons( Event const & ev ){
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(NO_EXTREMAL_JET_IDX);
}
return out_partons;
}
auto const & jets = ev.jets();
std::vector<fastjet::PseudoJet> extremal_jets;
if(! is_backward_g_to_h(ev)) {
auto most_backward = begin(jets);
// skip jets caused by unordered emission or qqbar
if(ev.type() == event_type::unob || ev.type() == event_type::qqbar_exb){
assert(jets.size() >= 2);
++most_backward;
}
extremal_jets.emplace_back(*most_backward);
}
if(! is_forward_g_to_h(ev)) {
auto most_forward = end(jets) - 1;
if(ev.type() == event_type::unof || ev.type() == event_type::qqbar_exf){
assert(jets.size() >= 2);
--most_forward;
}
extremal_jets.emplace_back(*most_forward);
}
const auto extremal_jet_indices = ev.particle_jet_indices(
extremal_jets
);
assert(extremal_jet_indices.size() == out_partons.size());
for(std::size_t i = 0; i < out_partons.size(); ++i){
assert(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 tree_kin_jets_qqbar_mid(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
double lambda
){
CLHEP::HepLorentzVector pq;
CLHEP::HepLorentzVector pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqbar
auto qqbart = q0;
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqbart -= to_HepLorentzVector(*parton_it);
}
const int nabove = std::distance(begin_ladder, backmidquark);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for (std::size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_qqbar_mid_current(
aptype, bptype, nabove, pa, pb,
pq, pqbar, partonsHLV
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqbart, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_qqbar(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda){
const bool swap_qqbar = is_quark(*BeginPart);
const auto pgin = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pq = to_HepLorentzVector(*(BeginPart+(swap_qqbar?0:1)));
const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_qqbar?1:0)));
const auto p1 = to_HepLorentzVector(*(BeginPart));
const auto pn = to_HepLorentzVector(*(EndPart-1));
assert((BeginIn)->type==pid::gluon); // Incoming a must be gluon.
const double current_factor = ME_qqbar_current(
(EndIn-1)->type, pgin, pq, pqbar, pn, pb, swap_qqbar
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pgin-pq-pqbar, pgin, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const double current_factor = ME_uno_current(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa
);
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pa-p1-pg, pa, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
} // namespace
double MatrixElement::tree_kin_jets(Event const & ev) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
if (ev.type()==event_type::FKL){
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
pa - p1, pa, pb, p1, pn,
param_.regulator_lambda
);
}
if (ev.type()==event_type::unordered_backward){
return tree_kin_jets_uno(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
if (ev.type()==event_type::unordered_forward){
return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
if (ev.type()==event_type::extremal_qqbar_backward){
return tree_kin_jets_qqbar(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
if (ev.type()==event_type::extremal_qqbar_forward){
return tree_kin_jets_qqbar(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
if (ev.type()==event_type::central_qqbar){
return tree_kin_jets_qqbar_mid(incoming[0].type, incoming[1].type,
to_HepLorentzVector(incoming[0]),
to_HepLorentzVector(incoming[1]),
partons, param_.regulator_lambda);
}
throw std::logic_error("Cannot reweight non-resummable processes in Pure Jets");
}
namespace {
double tree_kin_W_FKL(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
double lambda, ParticleProperties const & Wprop
){
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_uno(InIter BeginIn, partIter BeginPart,
partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
double lambda, ParticleProperties const & Wprop
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
bool wc = (BeginIn)->type==(BeginPart+1)->type; //leg b emits w
auto q0 = pa - p1 - pg;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_uno_current(
(BeginIn)->type, (BeginIn+1)->type, pn, pb,
p1, pa, pg, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_qqbar(InIter BeginIn, partIter BeginPart,
partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
double lambda, ParticleProperties const & Wprop
){
const bool swap_qqbar=is_quark(*BeginPart);
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pq = to_HepLorentzVector(*(BeginPart+(swap_qqbar?0:1)));
const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_qqbar?1:0)));
const auto p1 = to_HepLorentzVector(*(BeginPart));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const bool wc = (BeginIn+1)->type!=(EndPart-1)->type; //leg b emits w
auto q0 = pa - pq - pqbar;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_qqbar_current(
(BeginIn)->type, (BeginIn+1)->type, pa, pb,
pq, pqbar, pn, plbar, pl, swap_qqbar, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
double tree_kin_W_qqbar_mid(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
double lambda, ParticleProperties const & Wprop
){
CLHEP::HepLorentzVector pq;
CLHEP::HepLorentzVector pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons.front());
auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
// t-channel momentum after qqbar
auto qqbart = q0;
bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqbar emit W
bool wc = !wqq && (aptype==partons.front().type); //leg b emits w
assert(!wqq || !wc);
if(wqq){ // emission from qqbar
qqbart -= pl + plbar;
} else if(!wc) { // emission from leg a
q0 -= pl + plbar;
qqbart -= pl + plbar;
}
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqbart -= to_HepLorentzVector(*parton_it);
}
const int nabove = std::distance(begin_ladder, backmidquark);
const int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for (std::size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqbar_mid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqbart, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
} // namespace
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
#ifndef NDEBUG
// assert that there is exactly one decay corresponding to the W
assert(ev.decays().size() == 1);
auto const & w_boson{
std::find_if(ev.outgoing().cbegin(), ev.outgoing().cend(),
[] (Particle const & p) -> bool {
return std::abs(p.type) == ParticleID::Wp;
}) };
assert(w_boson != ev.outgoing().cend());
assert( static_cast<long int>(ev.decays().cbegin()->first)
== std::distance(ev.outgoing().cbegin(), w_boson) );
#endif
// find decay products of W
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert( ( is_anylepton(decay.at(0)) && is_anyneutrino(decay.at(1)) )
|| ( is_anylepton(decay.at(1)) && is_anyneutrino(decay.at(0)) ) );
// get lepton & neutrino
CLHEP::HepLorentzVector plbar;
CLHEP::HepLorentzVector pl;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
}
else{
pl = to_HepLorentzVector(decay.at(0));
plbar = to_HepLorentzVector(decay.at(1));
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
if(ev.type() == FKL){
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_backward){
return tree_kin_W_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_forward){
return tree_kin_W_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqbar_backward){
return tree_kin_W_qqbar(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqbar_forward){
return tree_kin_W_qqbar(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
assert(ev.type() == central_qqbar);
return tree_kin_W_qqbar_mid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
namespace {
std::vector <double> tree_kin_Z_FKL(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw
){
const auto p1 = to_HepLorentzVector(partons[0]);
const auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
const std::vector <double> current_factor = ME_Z_current(
aptype, bptype, pn, pb, p1, pa,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
q0, pa, pb, p1, pn, lambda));
} else if(is_gluon(aptype)){
// This is a gq event
const auto q0 = pa-p1;
ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
q0, pa, pb, p1, pn, lambda));
} else {
// This is a qq event
const auto q0 = pa-p1-plbar-pl;
const auto q1 = pa-p1;
ladder_factor=FKL_ladder_weight_mix(begin_ladder, end_ladder,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); ++i){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
template<class InIter, class partIter>
std::vector <double> tree_kin_Z_uno(InIter BeginIn, partIter BeginPart, partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const ParticleID aptype = (BeginIn)->type;
const ParticleID bptype = (BeginIn+1)->type;
const std::vector <double> current_factor = ME_Z_uno_current(
aptype, bptype, pn, pb, p1, pa, pg,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-pg-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else if(is_gluon(aptype)){
// This is a gq event
const auto q0 = pa-pg-p1;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a qq event
const auto q0 = pa-pg-p1-plbar-pl;
const auto q1 = pa-pg-p1;
ladder_factor=FKL_ladder_weight_mix(BeginPart+2, EndPart-1,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); ++i){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
} // namespace
std::vector<double> MatrixElement::tree_kin_Z(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
// find decay products of Z
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert(is_anylepton(decay.at(0)) && !is_anyneutrino(decay.at(0))
&& decay.at(0).type==-decay.at(1).type);
// get leptons
CLHEP::HepLorentzVector plbar;
CLHEP::HepLorentzVector pl;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
}
else{
pl = to_HepLorentzVector(decay.at(0));
plbar = to_HepLorentzVector(decay.at(1));
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
const double stw2 = param_.ew_parameters.sin2_tw();
const double ctw = param_.ew_parameters.cos_tw();
if(ev.type() == FKL){
return tree_kin_Z_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_backward){
return tree_kin_Z_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_forward){
return tree_kin_Z_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
throw std::logic_error("Can only reweight FKL or uno processes in Z+Jets");
}
double MatrixElement::tree_kin_Higgs(Event const & ev) const {
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev);
}
return tree_kin_Higgs_between(ev);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
){
if(type == pid::gluon) return currents::K_g(pout, pin);
return C_F;
}
} // namespace
// kinetic matrix element square for backward Higgs emission
// cf. eq:ME_h_jets_peripheral in developer manual,
// but without factors \alpha_s and the FKL ladder
double MatrixElement::MH2_backwardH(
const ParticleID type_forward,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pH,
CLHEP::HepLorentzVector const & pn
) const {
using namespace currents;
const double vev = param_.ew_parameters.vev();
return K(type_forward, pn, pb)*ME_jgH_j(
pH, pa, pn, pb,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb, vev
)/(4.*(N_C*N_C - 1)*(pa-pH).m2()*(pb-pn).m2());
}
// kinetic matrix element square for backward unordered emission
// and forward g -> Higgs
double MatrixElement::MH2_unob_forwardH(
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pH
) const {
using namespace currents;
const double vev = param_.ew_parameters.vev();
// Colour Acceleration Modifiers for quarks,
// see eq:K_q in developer manual
constexpr double K_f1 = C_F;
constexpr double nhel = 4.;
return K_f1*ME_juno_jgH(
pg, p1, pa, pH, pb,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb, vev
)/(nhel*(N_C*N_C - 1)*(pa - p1 - pg).m2()*(pb - pH).m2());
}
double MatrixElement::tree_kin_Higgs_first(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(is_anyquark(incoming.front())) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming.front());
const auto pb = to_HepLorentzVector(incoming.back());
const auto pH = to_HepLorentzVector(outgoing.front());
const auto end_ladder = end(partons) - ((ev.type() == event_type::unof)?2:1);
const auto pn = to_HepLorentzVector(*end_ladder);
const double ladder = FKL_ladder_weight(
begin(partons), end_ladder,
pa - pH, pa, pb, pH, pn,
param_.regulator_lambda
);
if(ev.type() == event_type::unof) {
const auto pg = to_HepLorentzVector(outgoing.back());
return MH2_unob_forwardH(
pb, pa, pg, pn, pH
)*ladder;
}
return MH2_backwardH(
incoming.back().type,
pa, pb, pH, pn
)*ladder;
}
double MatrixElement::tree_kin_Higgs_last(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.back().type == pid::Higgs);
if(is_anyquark(incoming.back())) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(ev);
}
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming.front());
const auto pb = to_HepLorentzVector(incoming.back());
const auto pH = to_HepLorentzVector(outgoing.back());
auto begin_ladder = begin(partons) + 1;
auto q0 = pa - to_HepLorentzVector(partons.front());
if(ev.type() == event_type::unob) {
q0 -= to_HepLorentzVector(*begin_ladder);
++begin_ladder;
}
const auto p1 = to_HepLorentzVector(*(begin_ladder - 1));
const double ladder = FKL_ladder_weight(
begin_ladder, end(partons),
q0, pa, pb, p1, pH,
param_.regulator_lambda
);
if(ev.type() == event_type::unob) {
const auto pg = to_HepLorentzVector(outgoing.front());
return MH2_unob_forwardH(
pa, pb, pg, p1, pH
)*ladder;
}
return MH2_backwardH(
incoming.front().type,
pb, pa, pH, p1
)*ladder;
}
namespace {
template<class InIter, class partIter>
double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart,
CLHEP::HepLorentzVector const & qH,
CLHEP::HepLorentzVector const & qHp1,
double mt, bool inc_bot, double mb, double vev
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
return ME_Higgs_current_uno(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa,
qH, qHp1, mt, inc_bot, mb, vev
);
}
} // namespace
double MatrixElement::tree_kin_Higgs_between(Event const & ev) const {
using namespace event_type;
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[(ev.type() == unob)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - ((ev.type() == unof)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
(ev.type() == unob)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
(ev.type() == unof)
|| partons.front().type != pid::gluon
))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor = NAN;
if(ev.type() == FKL){
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
}
else if(ev.type() == unob){
current_factor = C_A*C_A/2*tree_kin_Higgs_uno(
begin(incoming), end(incoming), begin(partons),
end(partons), qH, qH-pH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
current_factor = C_A*C_A/2*tree_kin_Higgs_uno(
rbegin(incoming), rend(incoming), rbegin(partons),
rend(partons), qH-pH, qH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
throw std::logic_error("Can only reweight FKL or uno processes in H+Jets");
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn,
param_.regulator_lambda
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn,
param_.regulator_lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
namespace {
double get_AWZH_coupling(Event const & ev, double alpha_s, double alpha_w) {
const auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](auto const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing())) return 1.;
switch(AWZH_boson->type){
case pid::Higgs:
return alpha_s*alpha_s;
case pid::Wp:
case pid::Wm:
case pid::Z_photon_mix:
return alpha_w*alpha_w;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
} // namespace
double MatrixElement::tree_param(Event const & ev, double mur) const {
assert(is_resummable(ev.type()));
const auto begin_partons = ev.begin_partons();
const auto end_partons = ev.end_partons();
const auto num_partons = std::distance(begin_partons, end_partons);
const double alpha_s = alpha_s_(mur);
const double gs2 = 4.*M_PI*alpha_s;
double res = std::pow(gs2, num_partons);
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(num_partons >= 2);
const auto first_emission = std::next(begin_partons);
const auto last_emission = std::prev(end_partons);
for(auto parton = first_emission; parton != last_emission; ++parton){
res *= 1. + alpha_s/(2.*M_PI)*BETA0*std::log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res;
}
} // namespace HEJ