diff --git a/include/HEJ/HDF5Reader.hh b/include/HEJ/HDF5Reader.hh
index 5b00fed..3b8e9e5 100644
--- a/include/HEJ/HDF5Reader.hh
+++ b/include/HEJ/HDF5Reader.hh
@@ -1,47 +1,48 @@
/** \file
* \brief Header file for reading events in the HDF5 event format.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include "HEJ/EventReader.hh"
namespace HEJ{
//! Class for reading events from a file in the HDF5 file format
/**
* @details This format is specified in \cite Hoeche:2019rti.
*/
class HDF5Reader : public EventReader{
public:
+ HDF5Reader() = delete;
+
//! Contruct object reading from the given file
explicit HDF5Reader(std::string const & filename);
//! Read an event
bool read_event() override;
//! Access header text
std::string const & header() const override;
//! Access run information
LHEF::HEPRUP const & heprup() const override;
//! Access last read event
LHEF::HEPEUP const & hepeup() const override;
//! Get number of events
HEJ::optional<size_t> number_events() const override;
+ ~HDF5Reader();
+
private:
struct HDF5ReaderImpl;
- struct HDF5ReaderImplDeleter {
- void operator()(HDF5ReaderImpl* p);
- };
- std::unique_ptr<HDF5ReaderImpl, HDF5ReaderImplDeleter> impl_;
+ std::unique_ptr<HDF5ReaderImpl> impl_;
};
}
diff --git a/include/HEJ/HDF5Writer.hh b/include/HEJ/HDF5Writer.hh
index 01103bf..718ea33 100644
--- a/include/HEJ/HDF5Writer.hh
+++ b/include/HEJ/HDF5Writer.hh
@@ -1,55 +1,54 @@
/** \file
* \brief Contains the EventWriter for HDF5 Output.
*
* The output format is specified in arXiv:1905.05120.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#include "HEJ/EventWriter.hh"
namespace LHEF {
class HEPRUP;
}
namespace HEJ{
class Event;
//! This is an event writer specifically for HDF5 output.
/**
* \internal Implementation note: This uses the pimpl ("pointer to
* implementation") idiom. HDF5 support is optional. Without pimpl,
* we would have to specify whether HDF5 is available via the
* preprocessor whenever this header is included. We don't want to
* burden users of the HEJ library (for example the HEJ fixed-order
* generator) with those details
*/
class HDF5Writer: public EventWriter{
public:
//! Constructor
/**
* @param file name of the output file
* @param heprup general process information
*/
HDF5Writer(std::string const & file, LHEF::HEPRUP heprup);
- ~HDF5Writer() override = default;
+ HDF5Writer() = delete;
//! Write an event to the output file
void write(Event const & ev) override;
+ ~HDF5Writer();
+
private:
struct HDF5WriterImpl;
- struct HDF5WriterImplDeleter {
- void operator()(HDF5WriterImpl* p);
- };
- std::unique_ptr<HDF5WriterImpl, HDF5WriterImplDeleter> impl_;
+ std::unique_ptr<HDF5WriterImpl> impl_;
};
}
diff --git a/include/HEJ/HepMC2Writer.hh b/include/HEJ/HepMC2Writer.hh
index 6dcc2e3..b64bc6b 100644
--- a/include/HEJ/HepMC2Writer.hh
+++ b/include/HEJ/HepMC2Writer.hh
@@ -1,54 +1,53 @@
/** \file
* \brief Contains the EventWriter for HepMC Output.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#include "HEJ/EventWriter.hh"
namespace LHEF {
class HEPRUP;
}
namespace HEJ{
class Event;
//! This is an event writer specifically for HepMC output.
/**
* \internal Implementation note:
* This uses the pimpl ("pointer to implementation") idiom.
* HepMC support is optional and the implementation depends on the
* HepMC version. Without pimpl, we would have to specify the HepMC version
* via the preprocessor whenever this header is included. We don't want to
* burden users of the HEJ library (for example the HEJ fixed-order generator)
* with those details
*/
class HepMC2Writer: public EventWriter{
public:
//! Constructor
/**
* @param file name of the output file
* @param heprup general process information
*/
HepMC2Writer(std::string const & file, LHEF::HEPRUP heprup);
- ~HepMC2Writer() override = default;
+ HepMC2Writer() = delete;
//! Write an event to the output file
void write(Event const & ev) override;
+ ~HepMC2Writer();
+
private:
struct HepMC2WriterImpl;
- struct HepMC2WriterImplDeleter {
- void operator()(HepMC2WriterImpl* p);
- };
- std::unique_ptr<HepMC2WriterImpl, HepMC2WriterImplDeleter> impl_;
+ std::unique_ptr<HepMC2WriterImpl> impl_;
};
}
diff --git a/include/HEJ/HepMC3Writer.hh b/include/HEJ/HepMC3Writer.hh
index fdc658c..05d3ca5 100644
--- a/include/HEJ/HepMC3Writer.hh
+++ b/include/HEJ/HepMC3Writer.hh
@@ -1,54 +1,52 @@
/** \file
* \brief Contains the EventWriter for HepMC3 Output.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#include "HEJ/EventWriter.hh"
namespace LHEF {
class HEPRUP;
}
namespace HEJ{
class Event;
//! This is an event writer specifically for HepMC3 output.
/**
* \internal Implementation note:
* This uses the pimpl ("pointer to implementation") idiom.
* HepMC3 support is optional and the implementation depends on the
* HepMC3 version. Without pimpl, we would have to specify the HepMC3 version
* via the preprocessor whenever this header is included. We don't want to
* burden users of the HEJ library (for example the HEJ fixed-order generator)
* with those details
*/
class HepMC3Writer: public EventWriter{
public:
//! Constructor
/**
* @param file name of the output file
* @param heprup general process information
*/
HepMC3Writer(std::string const & file, LHEF::HEPRUP heprup);
- ~HepMC3Writer() override = default;
+ HepMC3Writer() = delete;
//! Write an event to the output file
void write(Event const & ev) override;
+ ~HepMC3Writer();
private:
struct HepMC3WriterImpl;
- struct HepMC3WriterImplDeleter {
- void operator()(HepMC3WriterImpl* p);
- };
- std::unique_ptr<HepMC3WriterImpl, HepMC3WriterImplDeleter> impl_;
+ std::unique_ptr<HepMC3WriterImpl> impl_;
};
}
diff --git a/src/HDF5Reader.cc b/src/HDF5Reader.cc
index 4689ed2..3ed4b14 100644
--- a/src/HDF5Reader.cc
+++ b/src/HDF5Reader.cc
@@ -1,309 +1,304 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/HDF5Reader.hh"
#ifdef HEJ_BUILD_WITH_HDF5
#include <numeric>
#include <iterator>
#include "highfive/H5File.hpp"
namespace HEJ {
namespace {
// buffer size for reader
// each "reading from disk" reads "chunk_size" many event at once
constexpr std::size_t chunk_size = 10000;
struct ParticleData {
std::vector<int> id;
std::vector<int> status;
std::vector<int> mother1;
std::vector<int> mother2;
std::vector<int> color1;
std::vector<int> color2;
std::vector<double> px;
std::vector<double> py;
std::vector<double> pz;
std::vector<double> e;
std::vector<double> m;
std::vector<double> lifetime;
std::vector<double> spin;
};
struct EventRecords {
std::vector<int> particle_start;
std::vector<int> nparticles;
std::vector<int> pid;
std::vector<double> weight;
std::vector<double> scale;
std::vector<double> fscale;
std::vector<double> rscale;
std::vector<double> aqed;
std::vector<double> aqcd;
double trials;
ParticleData particles;
};
class ConstEventIterator {
public:
// iterator traits
using iterator_category = std::bidirectional_iterator_tag;
using value_type = LHEF::HEPEUP;
using difference_type = std::ptrdiff_t;
using pointer = const LHEF::HEPEUP*;
using reference = LHEF::HEPEUP const &;
using iterator = ConstEventIterator;
friend iterator cbegin(EventRecords const & records) noexcept;
friend iterator cend(EventRecords const & records) noexcept;
iterator& operator++() {
particle_offset_ += records_.get().nparticles[idx_];
++idx_;
return *this;
}
iterator& operator--() {
--idx_;
particle_offset_ -= records_.get().nparticles[idx_];
return *this;
}
iterator operator--(int) {
auto res = *this;
--(*this);
return res;
}
bool operator==(iterator const & other) const {
return idx_ == other.idx_;
}
bool operator!=(iterator other) const {
return !(*this == other);
}
value_type operator*() const {
value_type hepeup{};
auto const & r = records_.get();
hepeup.NUP = r.nparticles[idx_];
hepeup.IDPRUP = r.pid[idx_];
hepeup.XWGTUP = r.weight[idx_]/r.trials;
hepeup.weights.emplace_back(hepeup.XWGTUP, nullptr);
hepeup.SCALUP = r.scale[idx_];
hepeup.scales.muf = r.fscale[idx_];
hepeup.scales.mur = r.rscale[idx_];
hepeup.AQEDUP = r.aqed[idx_];
hepeup.AQCDUP = r.aqcd[idx_];
const size_t start = particle_offset_;
const size_t end = start + hepeup.NUP;
auto const & p = r.particles;
hepeup.IDUP = std::vector<long>( begin(p.id)+start, begin(p.id)+end );
hepeup.ISTUP = std::vector<int>( begin(p.status)+start, begin(p.status)+end );
hepeup.VTIMUP = std::vector<double>( begin(p.lifetime)+start, begin(p.lifetime)+end );
hepeup.SPINUP = std::vector<double>( begin(p.spin)+start, begin(p.spin)+end );
hepeup.MOTHUP.resize(hepeup.NUP);
hepeup.ICOLUP.resize(hepeup.NUP);
hepeup.PUP.resize(hepeup.NUP);
for(size_t i = 0; i < hepeup.MOTHUP.size(); ++i) {
const size_t idx = start + i;
assert(idx < end);
hepeup.MOTHUP[i] = std::make_pair(p.mother1[idx], p.mother2[idx]);
hepeup.ICOLUP[i] = std::make_pair(p.color1[idx], p.color2[idx]);
hepeup.PUP[i] = std::vector<double>{
p.px[idx], p.py[idx], p.pz[idx], p.e[idx], p.m[idx]
};
}
return hepeup;
}
private:
explicit ConstEventIterator(EventRecords const & records):
records_{records} {}
std::reference_wrapper<const EventRecords> records_;
size_t idx_ = 0;
size_t particle_offset_ = 0;
}; // end ConstEventIterator
ConstEventIterator cbegin(EventRecords const & records) noexcept {
return ConstEventIterator{records};
}
ConstEventIterator cend(EventRecords const & records) noexcept {
auto it =ConstEventIterator{records};
it.idx_ = records.aqcd.size(); // or size of any other records member
return it;
}
} // end anonymous namespace
struct HDF5Reader::HDF5ReaderImpl{
HighFive::File file;
std::size_t event_idx;
std::size_t particle_idx;
std::size_t nevents;
EventRecords records;
ConstEventIterator cur_event;
LHEF::HEPRUP heprup;
LHEF::HEPEUP hepeup;
explicit HDF5ReaderImpl(std::string const & filename):
file{filename},
event_idx{0},
particle_idx{0},
nevents{
file.getGroup("event")
.getDataSet("nparticles") // or any other dataset
.getSpace().getDimensions().front()
},
records{},
cur_event{cbegin(records)},
heprup{},
hepeup{}
{
read_heprup();
read_event_records(chunk_size);
}
void read_heprup() {
const auto init = file.getGroup("init");
init.getDataSet( "PDFgroupA" ).read(heprup.PDFGUP.first);
init.getDataSet( "PDFgroupB" ).read(heprup.PDFGUP.second);
init.getDataSet( "PDFsetA" ).read(heprup.PDFSUP.first);
init.getDataSet( "PDFsetB" ).read(heprup.PDFSUP.second);
init.getDataSet( "beamA" ).read(heprup.IDBMUP.first);
init.getDataSet( "beamB" ).read(heprup.IDBMUP.second);
init.getDataSet( "energyA" ).read(heprup.EBMUP.first);
init.getDataSet( "energyB" ).read(heprup.EBMUP.second);
init.getDataSet( "numProcesses" ).read(heprup.NPRUP);
init.getDataSet( "weightingStrategy" ).read(heprup.IDWTUP);
const auto proc_info = file.getGroup("procInfo");
proc_info.getDataSet( "procId" ).read(heprup.LPRUP);
proc_info.getDataSet( "xSection" ).read(heprup.XSECUP);
proc_info.getDataSet( "error" ).read(heprup.XERRUP);
// TODO: is this identification correct?
proc_info.getDataSet( "unitWeight" ).read(heprup.XMAXUP);
std::vector<double> trials;
file.getGroup("event").getDataSet("trials").read(trials);
records.trials = std::accumulate(begin(trials), end(trials), 0.);
}
std::size_t read_event_records(std::size_t count) {
count = std::min(count, nevents-event_idx);
auto events = file.getGroup("event");
events.getDataSet("nparticles").select({event_idx}, {count}).read(records.nparticles);
assert(records.nparticles.size() == count);
events.getDataSet("pid").select( {event_idx}, {count} ).read( records.pid );
events.getDataSet("weight").select( {event_idx}, {count} ).read( records.weight );
events.getDataSet("scale").select( {event_idx}, {count} ).read( records.scale );
events.getDataSet("fscale").select( {event_idx}, {count} ).read( records.fscale );
events.getDataSet("rscale").select( {event_idx}, {count} ).read( records.rscale );
events.getDataSet("aqed").select( {event_idx}, {count} ).read( records.aqed );
events.getDataSet("aqcd").select( {event_idx}, {count} ).read( records.aqcd );
const std::size_t particle_count = std::accumulate(
begin(records.nparticles), end(records.nparticles), 0
);
auto pdata = file.getGroup("particle");
auto & particles = records.particles;
pdata.getDataSet("id").select( {particle_idx}, {particle_count} ).read( particles.id );
pdata.getDataSet("status").select( {particle_idx}, {particle_count} ).read( particles.status );
pdata.getDataSet("mother1").select( {particle_idx}, {particle_count} ).read( particles.mother1 );
pdata.getDataSet("mother2").select( {particle_idx}, {particle_count} ).read( particles.mother2 );
pdata.getDataSet("color1").select( {particle_idx}, {particle_count} ).read( particles.color1 );
pdata.getDataSet("color2").select( {particle_idx}, {particle_count} ).read( particles.color2 );
pdata.getDataSet("px").select( {particle_idx}, {particle_count} ).read( particles.px );
pdata.getDataSet("py").select( {particle_idx}, {particle_count} ).read( particles.py );
pdata.getDataSet("pz").select( {particle_idx}, {particle_count} ).read( particles.pz );
pdata.getDataSet("e").select( {particle_idx}, {particle_count} ).read( particles.e );
pdata.getDataSet("m").select( {particle_idx}, {particle_count} ).read( particles.m );
pdata.getDataSet("lifetime").select( {particle_idx}, {particle_count} ).read( particles.lifetime );
pdata.getDataSet("spin").select( {particle_idx}, {particle_count} ).read( particles.spin );
event_idx += count;
particle_idx += particle_count;
return count;
}
};
HDF5Reader::HDF5Reader(std::string const & filename):
- impl_{
- new HDF5Reader::HDF5ReaderImpl{filename},
- HDF5Reader::HDF5ReaderImplDeleter{}
- }
+ impl_{std::make_unique<HDF5ReaderImpl>(filename)}
{}
bool HDF5Reader::read_event() {
if(impl_->cur_event == cend(impl_->records)) {
// end of active chunk, read new events from file
const auto events_read = impl_->read_event_records(chunk_size);
impl_->cur_event = cbegin(impl_->records);
if(events_read == 0) return false;
}
impl_->hepeup = *impl_->cur_event;
++impl_->cur_event;
return true;
}
namespace {
static const std::string nothing = "";
}
std::string const & HDF5Reader::header() const {
return nothing;
}
LHEF::HEPRUP const & HDF5Reader::heprup() const {
return impl_->heprup;
}
LHEF::HEPEUP const & HDF5Reader::hepeup() const {
return impl_->hepeup;
}
HEJ::optional<size_t> HDF5Reader::number_events() const {
return impl_->nevents;
}
}
#else // no HDF5 support
namespace HEJ {
class HDF5Reader::HDF5ReaderImpl{};
HDF5Reader::HDF5Reader(std::string const &){
throw std::invalid_argument{
"Failed to create HDF5 reader: "
"HEJ 2 was built without HDF5 support"
};
}
bool HDF5Reader::read_event() {
throw std::logic_error{"unreachable"};
}
std::string const & HDF5Reader::header() const {
throw std::logic_error{"unreachable"};
}
LHEF::HEPRUP const & HDF5Reader::heprup() const {
throw std::logic_error{"unreachable"};
}
LHEF::HEPEUP const & HDF5Reader::hepeup() const {
throw std::logic_error{"unreachable"};
}
HEJ::optional<size_t> HDF5Reader::number_events() const {
throw std::logic_error{"unreachable"};
}
}
#endif
namespace HEJ {
- void HDF5Reader::HDF5ReaderImplDeleter::operator()(HDF5ReaderImpl* p) {
- delete p;
- }
+ HDF5Reader::~HDF5Reader() = default;
}
diff --git a/src/HDF5Writer.cc b/src/HDF5Writer.cc
index c61d1e0..032db8d 100644
--- a/src/HDF5Writer.cc
+++ b/src/HDF5Writer.cc
@@ -1,386 +1,381 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/HDF5Writer.hh"
#include <cassert>
#include "LHEF/LHEF.h"
#ifdef HEJ_BUILD_WITH_HDF5
#include <type_traits>
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "highfive/H5File.hpp"
namespace HEJ{
using HighFive::File;
using HighFive::DataSpace;
namespace{
constexpr std::size_t chunk_size = 1000;
constexpr unsigned compression_level = 3;
size_t to_index(event_type::EventType const type){
return type==0?0:floor(log2(type))+1;
}
template<typename T>
void write_dataset(HighFive::Group & group, std::string const & name, T val) {
using data_t = std::decay_t<T>;
group.createDataSet<data_t>(name, DataSpace::From(val)).write(val);
}
template<typename T>
void write_dataset(
HighFive::Group & group, std::string const & name,
std::vector<T> const & val
) {
using data_t = std::decay_t<T>;
group.createDataSet<data_t>(name, DataSpace::From(val)).write(val);
}
struct Cache {
explicit Cache(size_t capacity):
capacity{capacity}
{
nparticles.reserve(capacity);
start.reserve(capacity);
pid.reserve(capacity);
weight.reserve(capacity);
scale.reserve(capacity);
fscale.reserve(capacity);
rscale.reserve(capacity);
aqed.reserve(capacity);
aqcd.reserve(capacity);
trials.reserve(capacity);
}
void fill(HEJ::Event ev) {
const auto hepeup = to_HEPEUP(ev, nullptr);
fill_event_params(hepeup);
fill_event_particles(hepeup);
}
void fill_event_params(LHEF::HEPEUP const & ev) {
nparticles.emplace_back(ev.NUP);
start.emplace_back(particle_pos);
pid.emplace_back(ev.IDPRUP);
weight.emplace_back(ev.XWGTUP);
scale.emplace_back(ev.SCALUP);
fscale.emplace_back(ev.scales.muf);
rscale.emplace_back(ev.scales.mur);
aqed.emplace_back(ev.AQEDUP);
aqcd.emplace_back(ev.AQCDUP);
// set first trial=1 for first event
// -> sum(trials) = 1 -> xs=sum(weights)/sum(trials) as in Sherpa
if(particle_pos == 0){
trials.emplace_back(1.);
} else {
trials.emplace_back(0.);
}
particle_pos += ev.NUP;
}
void fill_event_particles(LHEF::HEPEUP const & ev) {
id.insert(end(id), begin(ev.IDUP), end(ev.IDUP));
status.insert(end(status), begin(ev.ISTUP), end(ev.ISTUP));
lifetime.insert(end(lifetime), begin(ev.VTIMUP), end(ev.VTIMUP));
spin.insert(end(spin), begin(ev.SPINUP), end(ev.SPINUP));
for(int i = 0; i < ev.NUP; ++i) {
mother1.emplace_back(ev.MOTHUP[i].first);
mother2.emplace_back(ev.MOTHUP[i].second);
color1.emplace_back(ev.ICOLUP[i].first);
color2.emplace_back(ev.ICOLUP[i].second);
px.emplace_back(ev.PUP[i][0]);
py.emplace_back(ev.PUP[i][1]);
pz.emplace_back(ev.PUP[i][2]);
e.emplace_back(ev.PUP[i][3]);
m.emplace_back(ev.PUP[i][4]);
}
}
bool is_full() const {
return nparticles.size() >= capacity;
}
void clear() {
nparticles.clear();
start.clear();
pid.clear();
id.clear();
status.clear();
mother1.clear();
mother2.clear();
color1.clear();
color2.clear();
weight.clear();
scale.clear();
fscale.clear();
rscale.clear();
aqed.clear();
aqcd.clear();
trials.clear();
px.clear();
py.clear();
pz.clear();
e.clear();
m.clear();
lifetime.clear();
spin.clear();
}
size_t capacity;
std::vector<int> nparticles, start, pid, id, status,
mother1, mother2, color1, color2;
std::vector<double> weight, scale, fscale, rscale, aqed,
aqcd, trials, px, py, pz, e, m, lifetime, spin;
private:
size_t particle_pos = 0;
};
}
struct HDF5Writer::HDF5WriterImpl{
File file;
LHEF::HEPRUP heprup;
Cache cache{chunk_size};
size_t event_idx = 0;
size_t particle_idx = 0;
HDF5WriterImpl(std::string const & filename, LHEF::HEPRUP && hepr):
file{filename, File::ReadWrite | File::Create | File::Truncate},
heprup{std::move(hepr)}
{
// TODO: code duplication with Les Houches Writer
const int max_number_types = to_index(event_type::last_type)+1;
heprup.NPRUP = max_number_types;
// ids of event types
heprup.LPRUP.clear();
heprup.LPRUP.reserve(max_number_types);
heprup.LPRUP.emplace_back(0);
for(size_t i=event_type::first_type+1; i<=event_type::last_type; i*=2) {
heprup.LPRUP.emplace_back(i);
}
heprup.XSECUP = std::vector<double>(max_number_types, 0.);
heprup.XERRUP = std::vector<double>(max_number_types, 0.);
heprup.XMAXUP = std::vector<double>(max_number_types, 0.);
write_init();
create_event_group();
create_particle_group();
}
void write_init() {
auto init = file.createGroup("init");
write_dataset(init, "PDFgroupA" , heprup.PDFGUP.first);
write_dataset(init, "PDFgroupB" , heprup.PDFGUP.second);
write_dataset(init, "PDFsetA" , heprup.PDFSUP.first);
write_dataset(init, "PDFsetB" , heprup.PDFSUP.second);
write_dataset(init, "beamA" , heprup.IDBMUP.first);
write_dataset(init, "beamB" , heprup.IDBMUP.second);
write_dataset(init, "energyA" , heprup.EBMUP.first);
write_dataset(init, "energyB" , heprup.EBMUP.second);
write_dataset(init, "numProcesses" , heprup.NPRUP);
write_dataset(init, "weightingStrategy", heprup.IDWTUP);
auto proc_info = file.createGroup("procInfo");
write_dataset(proc_info, "procId", heprup.LPRUP);
}
static HighFive::DataSetCreateProps const & hdf5_chunk() {
static const auto props = [](){
HighFive::DataSetCreateProps props;
props.add(HighFive::Chunking({chunk_size}));
props.add(HighFive::Deflate(compression_level));
return props;
}();
return props;
}
void create_event_group() {
static const auto dim = DataSpace({0}, {DataSpace::UNLIMITED});
auto events = file.createGroup("event");
events.createDataSet<int>("nparticles", dim, hdf5_chunk());
events.createDataSet<int>("start", dim, hdf5_chunk());
events.createDataSet<int>("pid", dim, hdf5_chunk());
events.createDataSet<double>("weight", dim, hdf5_chunk());
events.createDataSet<double>("scale", dim, hdf5_chunk());
events.createDataSet<double>("fscale", dim, hdf5_chunk());
events.createDataSet<double>("rscale", dim, hdf5_chunk());
events.createDataSet<double>("aqed", dim, hdf5_chunk());
events.createDataSet<double>("aqcd", dim, hdf5_chunk());
events.createDataSet<double>("trials", dim, hdf5_chunk());
}
void resize_event_group(size_t new_size) {
auto events = file.getGroup("event");
events.getDataSet("nparticles").resize({new_size});
events.getDataSet("start").resize({new_size});
events.getDataSet("pid").resize({new_size});
events.getDataSet("weight").resize({new_size});
events.getDataSet("scale").resize({new_size});
events.getDataSet("fscale").resize({new_size});
events.getDataSet("rscale").resize({new_size});
events.getDataSet("aqed").resize({new_size});
events.getDataSet("aqcd").resize({new_size});
events.getDataSet("trials").resize({new_size});
}
void create_particle_group() {
static const auto dim = DataSpace({0}, {DataSpace::UNLIMITED});
auto particles = file.createGroup("particle");
particles.createDataSet<int>("id", dim, hdf5_chunk());
particles.createDataSet<int>("status", dim, hdf5_chunk());
particles.createDataSet<int>("mother1", dim, hdf5_chunk());
particles.createDataSet<int>("mother2", dim, hdf5_chunk());
particles.createDataSet<int>("color1", dim, hdf5_chunk());
particles.createDataSet<int>("color2", dim, hdf5_chunk());
particles.createDataSet<double>("px", dim, hdf5_chunk());
particles.createDataSet<double>("py", dim, hdf5_chunk());
particles.createDataSet<double>("pz", dim, hdf5_chunk());
particles.createDataSet<double>("e", dim, hdf5_chunk());
particles.createDataSet<double>("m", dim, hdf5_chunk());
particles.createDataSet<double>("lifetime", dim, hdf5_chunk());
particles.createDataSet<double>("spin", dim, hdf5_chunk());
}
void resize_particle_group(size_t new_size) {
auto particles = file.getGroup("particle");
particles.getDataSet("id").resize({new_size});
particles.getDataSet("status").resize({new_size});
particles.getDataSet("mother1").resize({new_size});
particles.getDataSet("mother2").resize({new_size});
particles.getDataSet("color1").resize({new_size});
particles.getDataSet("color2").resize({new_size});
particles.getDataSet("px").resize({new_size});
particles.getDataSet("py").resize({new_size});
particles.getDataSet("pz").resize({new_size});
particles.getDataSet("e").resize({new_size});
particles.getDataSet("m").resize({new_size});
particles.getDataSet("lifetime").resize({new_size});
particles.getDataSet("spin").resize({new_size});
}
void write(Event const & ev){
cache.fill(ev);
if(cache.is_full()) {
dump_cache();
}
const double wt = ev.central().weight;
const size_t idx = to_index(ev.type());
heprup.XSECUP[idx] += wt;
heprup.XERRUP[idx] += wt*wt;
if(wt > heprup.XMAXUP[idx]){
heprup.XMAXUP[idx] = wt;
}
}
void dump_cache() {
write_event_params();
write_event_particles();
cache.clear();
}
void write_event_params() {
auto events = file.getGroup("event");
// choose arbitrary dataset to find size
const auto dataset = events.getDataSet("nparticles");
const size_t size = dataset.getSpace().getDimensions().front();
resize_event_group(size + cache.nparticles.size());
#define WRITE_FROM_CACHE(GROUP, PROPERTY) \
GROUP.getDataSet(#PROPERTY).select({size}, {cache.PROPERTY.size()}).write(cache.PROPERTY)
WRITE_FROM_CACHE(events, nparticles);
WRITE_FROM_CACHE(events, start);
WRITE_FROM_CACHE(events, pid);
WRITE_FROM_CACHE(events, weight);
WRITE_FROM_CACHE(events, scale);
WRITE_FROM_CACHE(events, fscale);
WRITE_FROM_CACHE(events, rscale);
WRITE_FROM_CACHE(events, aqed);
WRITE_FROM_CACHE(events, aqcd);
WRITE_FROM_CACHE(events, trials);
}
void write_event_particles() {
auto particles = file.getGroup("particle");
// choose arbitrary dataset to find size
const auto dataset = particles.getDataSet("id");
const size_t size = dataset.getSpace().getDimensions().front();
resize_particle_group(size + cache.id.size());
WRITE_FROM_CACHE(particles, id);
WRITE_FROM_CACHE(particles, status);
WRITE_FROM_CACHE(particles, lifetime);
WRITE_FROM_CACHE(particles, spin);
WRITE_FROM_CACHE(particles, mother1);
WRITE_FROM_CACHE(particles, mother2);
WRITE_FROM_CACHE(particles, color1);
WRITE_FROM_CACHE(particles, color2);
WRITE_FROM_CACHE(particles, px);
WRITE_FROM_CACHE(particles, py);
WRITE_FROM_CACHE(particles, pz);
WRITE_FROM_CACHE(particles, e);
WRITE_FROM_CACHE(particles, m);
}
#undef WRITE_FROM_CACHE
~HDF5WriterImpl(){
dump_cache();
auto proc_info = file.getGroup("procInfo");
write_dataset(proc_info, "xSection", heprup.XSECUP);
write_dataset(proc_info, "error", heprup.XERRUP);
write_dataset(proc_info, "unitWeight", heprup.XMAXUP);
}
};
HDF5Writer::HDF5Writer(std::string const & filename, LHEF::HEPRUP heprup):
- impl_{
- new HDF5Writer::HDF5WriterImpl{filename, std::move(heprup)},
- HDF5Writer::HDF5WriterImplDeleter{}
- }
+ impl_{new HDF5Writer::HDF5WriterImpl{filename, std::move(heprup)}}
{}
void HDF5Writer::write(Event const & ev){
impl_->write(ev);
}
}
#else // no HDF5 support
namespace HEJ{
class HDF5Writer::HDF5WriterImpl{};
HDF5Writer::HDF5Writer(std::string const &, LHEF::HEPRUP){
throw std::invalid_argument{
"Failed to create HDF5 writer: "
"HEJ 2 was built without HDF5 support"
};
}
void HDF5Writer::write(Event const &){
assert(false);
}
}
#endif
namespace HEJ {
- void HDF5Writer::HDF5WriterImplDeleter::operator()(HDF5WriterImpl* p) {
- delete p;
- }
+ HDF5Writer::~HDF5Writer() = default;
}
diff --git a/src/HepMC2Writer.cc b/src/HepMC2Writer.cc
index f603d5d..32d92f0 100644
--- a/src/HepMC2Writer.cc
+++ b/src/HepMC2Writer.cc
@@ -1,85 +1,80 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/HepMC2Writer.hh"
#include <cassert>
#include "LHEF/LHEF.h"
#ifdef HEJ_BUILD_WITH_HepMC2
#include "HepMC/IO_GenEvent.h"
#include <utility>
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/HepMC2Interface.hh"
namespace HEJ{
struct HepMC2Writer::HepMC2WriterImpl{
HepMC2Interface hepmc_;
HepMC2WriterImpl & operator=(HepMC2WriterImpl const & other) = delete;
HepMC2WriterImpl(HepMC2WriterImpl const & other) = delete;
HepMC2WriterImpl & operator=(HepMC2WriterImpl && other) = delete;
HepMC2WriterImpl(HepMC2WriterImpl && other) = delete;
HepMC::IO_GenEvent writer_;
HepMC2WriterImpl(
std::string const & file, LHEF::HEPRUP && heprup
):
hepmc_(heprup),
writer_{file}
{}
void write(Event const & ev){
auto out_ev = hepmc_(ev);
writer_.write_event(&out_ev);
}
};
- void HepMC2Writer::HepMC2WriterImplDeleter::operator()(HepMC2WriterImpl* p) {
- delete p;
- }
-
HepMC2Writer::HepMC2Writer(std::string const & file, LHEF::HEPRUP heprup):
- impl_{std::unique_ptr<HepMC2WriterImpl, HepMC2WriterImplDeleter>{
- new HepMC2WriterImpl(file, std::move(heprup))
- }}
+ impl_{new HepMC2WriterImpl{file, std::move(heprup)}}
{}
void HepMC2Writer::write(Event const & ev){
impl_->write(ev);
}
} // namespace HEJ
#else // no HepMC2
namespace HEJ{
class HepMC2Writer::HepMC2WriterImpl{};
HepMC2Writer::HepMC2Writer(std::string const &, LHEF::HEPRUP){
throw std::invalid_argument(
"Failed to create HepMC writer: "
"HEJ 2 was built without HepMC2 support"
);
}
void HepMC2Writer::write(Event const &){
assert(false);
}
- void HepMC2Writer::HepMC2WriterImplDeleter::operator()(HepMC2WriterImpl* p) {
- delete p;
- }
}
#endif
+
+namespace HEJ{
+ HepMC2Writer::~HepMC2Writer() = default;
+}
diff --git a/src/HepMC3Writer.cc b/src/HepMC3Writer.cc
index d646484..143bd65 100644
--- a/src/HepMC3Writer.cc
+++ b/src/HepMC3Writer.cc
@@ -1,122 +1,117 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/HepMC3Writer.hh"
#include <cassert>
#include "LHEF/LHEF.h"
#ifdef HEJ_BUILD_WITH_HepMC3
#include "HepMC3/LHEFAttributes.h"
#include "HepMC3/WriterAscii.h"
#include "HEJ/Version.hh"
#include <utility>
#include "HepMC3/GenParticle.h"
#include "HepMC3/GenVertex.h"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/HepMC3Interface.hh"
namespace {
void reset_weight_info(LHEF::HEPRUP & heprup){
heprup.IDWTUP = 2;
// use placeholders for unknown init block values
// we can overwrite them after processing all events
heprup.XSECUP = {0.};
heprup.XERRUP = {0.};
heprup.XMAXUP = {0.};
}
HepMC3::shared_ptr<HepMC3::GenRunInfo> init_runinfo(LHEF::HEPRUP && heprup){
reset_weight_info(heprup);
HepMC3::GenRunInfo runinfo;
auto hepr = HepMC3::make_shared<HepMC3::HEPRUPAttribute>();
hepr->heprup = heprup;
runinfo.add_attribute(std::string("HEPRUP"), hepr);
for (int i = 0, N = hepr->heprup.generators.size(); i < N; ++i ){
HepMC3::GenRunInfo::ToolInfo tool;
tool.name = hepr->heprup.generators[i].name;
tool.version = hepr->heprup.generators[i].version;
tool.description = hepr->heprup.generators[i].contents;
runinfo.tools().push_back(tool);
}
return HepMC3::make_shared<HepMC3::GenRunInfo>(runinfo);
}
} // namespace anonymous
namespace HEJ{
struct HepMC3Writer::HepMC3WriterImpl{
HepMC3Interface HepMC3_;
HepMC3WriterImpl & operator=(HepMC3WriterImpl const & other) = delete;
HepMC3WriterImpl(HepMC3WriterImpl const & other) = delete;
HepMC3WriterImpl & operator=(HepMC3WriterImpl && other) = delete;
HepMC3WriterImpl(HepMC3WriterImpl && other) = delete;
HepMC3::WriterAscii writer_;
HepMC3WriterImpl(
std::string const & file, LHEF::HEPRUP && heprup
):
HepMC3_(heprup),
writer_{file, init_runinfo(std::move(heprup))}
{}
~HepMC3WriterImpl(){
writer_.close();
}
void write(Event const & ev){
auto out_ev = HepMC3_(ev);
writer_.write_event(out_ev);
}
};
- void HepMC3Writer::HepMC3WriterImplDeleter::operator()(HepMC3WriterImpl* p) {
- delete p;
- }
-
HepMC3Writer::HepMC3Writer(std::string const & file, LHEF::HEPRUP heprup):
- impl_{std::unique_ptr<HepMC3WriterImpl, HepMC3WriterImplDeleter>{
- new HepMC3WriterImpl(file, std::move(heprup))
- }}
+ impl_{new HepMC3WriterImpl(file, std::move(heprup))}
{}
void HepMC3Writer::write(Event const & ev){
impl_->write(ev);
}
} // namespace HEJ
#else // no HepMC3
namespace HEJ{
class HepMC3Writer::HepMC3WriterImpl{};
HepMC3Writer::HepMC3Writer(std::string const &, LHEF::HEPRUP){
throw std::invalid_argument(
"Failed to create HepMC3 writer: "
"HEJ 2 was built without HepMC3 support"
);
}
void HepMC3Writer::write(Event const &){
assert(false);
}
- void HepMC3Writer::HepMC3WriterImplDeleter::operator()(HepMC3WriterImpl* p) {
- delete p;
- }
}
#endif
+
+namespace HEJ{
+ HepMC3Writer::~HepMC3Writer() = default;
+}
diff --git a/src/Hjets.cc b/src/Hjets.cc
index e2fa700..61e64f9 100644
--- a/src/Hjets.cc
+++ b/src/Hjets.cc
@@ -1,1084 +1,1084 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include "HEJ/Hjets.hh"
#include <assert.h>
#include <limits>
#include "HEJ/Constants.hh"
#ifdef HEJ_BUILD_WITH_QCDLOOP
#include "qcdloop/qcdloop.h"
#endif
const COM looprwfactor = (COM(0.,1.)*M_PI*M_PI)/pow((2.*M_PI),4);
constexpr double infinity = std::numeric_limits<double>::infinity();
namespace {
// Loop integrals
#ifdef HEJ_BUILD_WITH_QCDLOOP
COM B0DD(HLV q, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_B0 = [](){
ql::Bubble<std::complex<double>,double,double> ql_B0;
ql_B0.setCacheSize(100);
return ql_B0;
}();
static std::vector<double> masses(2);
static std::vector<double> momenta(1);
for(auto & m: masses) m = mq*mq;
momenta.front() = q.m2();
ql_B0.integral(result, 1, masses, momenta);
return result[0];
}
COM C0DD(HLV q1, HLV q2, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_C0 = [](){
ql::Triangle<std::complex<double>,double,double> ql_C0;
ql_C0.setCacheSize(100);
return ql_C0;
}();
static std::vector<double> masses(3);
static std::vector<double> momenta(3);
for(auto & m: masses) m = mq*mq;
momenta[0] = q1.m2();
momenta[1] = q2.m2();
momenta[2] = (q1+q2).m2();
ql_C0.integral(result, 1, masses, momenta);
return result[0];
}
COM D0DD(HLV q1,HLV q2, HLV q3, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_D0 = [](){
ql::Box<std::complex<double>,double,double> ql_D0;
ql_D0.setCacheSize(100);
return ql_D0;
}();
static std::vector<double> masses(4);
static std::vector<double> momenta(6);
for(auto & m: masses) m = mq*mq;
momenta[0] = q1.m2();
momenta[1] = q2.m2();
momenta[2] = q3.m2();
momenta[3] = (q1+q2+q3).m2();
momenta[4] = (q1+q2).m2();
momenta[5] = (q2+q3).m2();
ql_D0.integral(result, 1, masses, momenta);
return result[0];
}
COM A1(HLV q1, HLV q2, double mt)
// As given in Eq. (B.2) of VDD
{
double q12,q22,Q2;
HLV Q;
double Delta3,mt2;
COM ans(COM(0.,0.));
q12=q1.m2();
q22=q2.m2();
Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
Q2=Q.m2();
Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
assert(mt > 0.);
mt2=mt*mt;
ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*( 4.*mt2/Delta3*(Q2-q12-q22)
-1.-4.*q12*q22/Delta3-12.*q12*q22*Q2/Delta3/Delta3*(q12+q22-Q2) )
- looprwfactor*COM(0,-1)*( B0DD(q2,mt)-B0DD(Q,mt) )
* ( 2.*q22/Delta3+12.*q12*q22/Delta3/Delta3*(q22-q12+Q2) )
- looprwfactor*COM(0,-1)*( B0DD(q1,mt)-B0DD(Q,mt) )
* ( 2.*q12/Delta3+12.*q12*q22/Delta3/Delta3*(q12-q22+Q2) )
- 2./Delta3/16/M_PI/M_PI*(q12+q22-Q2);
return ans;
}
COM A2(HLV q1, HLV q2, double mt)
// As given in Eq. (B.2) of VDD, but with high energy limit
// of invariants taken.
{
double q12,q22,Q2;
HLV Q;
double Delta3,mt2;
COM ans(COM(0.,0.));
assert(mt > 0.);
mt2=mt*mt;
q12=q1.m2();
q22=q2.m2();
Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
Q2=Q.m2();
Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*( 2.*mt2+1./2.*(q12+q22-Q2)
+2.*q12*q22*Q2/Delta3 )
+looprwfactor*COM(0,-1)*(B0DD(q2,mt)-B0DD(Q,mt))
*q22*(q22-q12-Q2)/Delta3
+looprwfactor*COM(0,-1)*(B0DD(q1,mt)-B0DD(Q,mt))
*q12*(q12-q22-Q2)/Delta3+1./16/M_PI/M_PI;
return ans;
}
#else // no QCDloop
COM A1(HLV, HLV, double) {
throw std::logic_error{"A1 called without QCDloop support"};
}
COM A2(HLV, HLV, double) {
throw std::logic_error{"A2 called without QCDloop support"};
}
#endif
void to_current(const HLV & q, current & ret){
ret[0]=q.e();
ret[1]=q.x();
ret[2]=q.y();
ret[3]=q.z();
}
/**
* @brief Higgs vertex contracted with current @param C1 and @param C2
*/
COM cHdot(const current & C1, const current & C2, const current & q1,
const current & q2, double mt, bool incBot, double mb, double vev)
{
if (mt == infinity) {
return (cdot(C1,C2)*cdot(q1,q2)-cdot(C1,q2)*cdot(C2,q1))/(3*M_PI*vev);
}
else {
HLV vq1,vq2;
vq1.set(q1[1].real(),q1[2].real(),q1[3].real(),q1[0].real());
vq2.set(q2[1].real(),q2[2].real(),q2[3].real(),q2[0].real());
// first minus sign obtained because of q1-difference to VDD
// Factor is because 4 mt^2 g^2/vev A1 -> 16 pi mt^2/vev alphas,
if(!(incBot))
return 16.*M_PI*mt*mt/vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)
-cdot(C1,C2)*A2(-vq1,vq2,mt));
else
return 16.*M_PI*mt*mt/vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)
-cdot(C1,C2)*A2(-vq1,vq2,mt))
+ 16.*M_PI*mb*mb/vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mb)
-cdot(C1,C2)*A2(-vq1,vq2,mb));
}
}
//@{
/**
* @brief Higgs+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param q1 t-channel momenta into higgs vertex
* @param q2 t-channel momenta out of higgs vertex
* @param mt top mass (inf or value)
* @param incBot Bool, to include bottom mass (true) or not (false)?
* @param mb bottom mass (value)
* @param pg Unordered Gluon momenta
*
* Calculates j^\mu H j_\mu. FKL with higgs vertex somewhere in the FKL chain.
* Handles all possible incoming states.
*/
double j_h_j(HLV const & p1out, HLV const & p1in, HLV const & p2out,
HLV const & p2in, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev
){
current j1p,j1m,j2p,j2m, q1v, q2v;
// Note need to flip helicities in anti-quark case.
joi(p1out, false, p1in, false, j1p);
joi(p1out, true, p1in, true, j1m);
joi(p2out, false, p2in, false, j2p);
joi(p2out, true, p2in, true, j2m);
to_current(q1, q1v);
to_current(q2, q2v);
COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb, vev);
COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb, vev);
COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb, vev);
COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb, vev);
// average over helicities
const double sst=(abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm))/4.;
return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
}
} // namespace anonymous
double ME_H_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev)
* K_g(p2out,p2in)/HEJ::C_A;
}
double ME_H_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev)
* K_g(p2out,p2in)/HEJ::C_A;
}
double ME_H_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev)
* K_g(p2out,p2in)/HEJ::C_A * K_g(p1out,p1in)/HEJ::C_A;
}
//@}
namespace {
//@{
/// @brief Higgs vertex contracted with one current
CCurrent jH(HLV const & pout, bool helout, HLV const & pin,
bool helin, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev)
{
CCurrent j2 = joi(pout,helout,pin,helin);
CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz());
if(mt == infinity)
return ((q1.dot(q2))*j2 - j2.dot(q1)*jq2)/(3*M_PI*vev);
else
{
if(incBot)
return (-16.*M_PI*mb*mb/vev*j2.dot(q1)*jq2*A1(-q1,q2,mb)
-16.*M_PI*mb*mb/vev*j2*A2(-q1,q2,mb))
+ (-16.*M_PI*mt*mt/vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)
-16.*M_PI*mt*mt/vev*j2*A2(-q1,q2,mt));
else
return (-16.*M_PI*mt*mt/vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)
-16.*M_PI*mt*mt/vev*j2*A2(-q1,q2,mt));
}
}
//@}
//@{
/**
* @brief Higgs+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (W emission)
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param q1 t-channel momenta into higgs vertex
* @param q2 t-channel momenta out of higgs vertex
* @param mt top mass (inf or value)
* @param incBot Bool, to include bottom mass (true) or not (false)?
* @param mb bottom mass (value)
*
* Calculates j_{uno}^\mu H j_\mu. Unordered with higgs vertex
* somewhere in the FKL chain. Handles all possible incoming states.
*/
double juno_h_j(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in,
HLV const & qH1, HLV const & qH2,
double mt, bool incBot, double mb, double vev
){
// This construction is taking rapidity order: pg > p1out >> p2out
HLV q1=p1in-p1out; // Top End
HLV q2=-(p2in-p2out); // Bottom End
HLV qg=p1in-p1out-pg; // Extra bit post-gluon
// Note <p1|eps|pa> current split into two by gauge choice.
// See James C's Thesis (p72). <p1|eps|pa> -> <p1|pg><pg|pa>
CCurrent mj1p=joi(p1out,false, p1in,false);
CCurrent mj1m=joi(p1out, true, p1in, true);
CCurrent jgap=joi(pg, false, p1in,false);
CCurrent jgam=joi(pg, true, p1in, true);
// Note for function joo(): <p1+|pg+> = <pg-|p1->.
CCurrent j2gp=joo(p1out, false, pg, false);
CCurrent j2gm=joo(p1out, true, pg, true);
CCurrent mjH2p=jH(p2out, true,p2in, true,qH1,qH2,mt,incBot,mb, vev);
CCurrent mjH2m=jH(p2out,false,p2in,false,qH1,qH2,mt,incBot,mb, vev);
// Dot products of these which occur again and again
COM MHmp=mj1m.dot(mjH2p);
COM MHmm=mj1m.dot(mjH2m);
COM MHpp=mj1p.dot(mjH2p);
COM MHpm=mj1p.dot(mjH2m);
CCurrent p2o(p2out), p2i(p2in), p1o(p1out), p1i(p1in), qsum(q1+qg);
CCurrent Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m + 2.*mj1m.dot(pg)*mjH2m
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHmm/2.) )/q1.m2();
CCurrent Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m + 2.*mj1m.dot(pg)*mjH2p
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHmp/2.) )/q1.m2();
CCurrent Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p + 2.*mj1p.dot(pg)*mjH2m
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHpm/2.) )/q1.m2();
CCurrent Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p + 2.*mj1p.dot(pg)*mjH2p
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHpp/2.) )/q1.m2();
CCurrent U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2();
CCurrent U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2();
CCurrent U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2();
CCurrent U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2();
CCurrent U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2();
CCurrent U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
CCurrent U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
CCurrent U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(p1in-pg).m2();
constexpr double cf=HEJ::C_F;
double amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
double amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
double apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
double app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());
// Now add the t-channels for the Higgs
ampsq/=qH1.m2()*qg.m2();
ampsq/=16.;
// Factor of (Cf/Ca) for each quark to match ME_H_qQ.
ampsq*=HEJ::C_F*HEJ::C_F/HEJ::C_A/HEJ::C_A;
return ampsq;
}
} // namespace anonymous
double ME_H_unob_qQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV qH1,
HLV qH2, double mt, bool incBot, double mb, double vev){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_qbarQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb, double vev){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_qQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb, double vev){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_qbarQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb, double vev){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_gQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb, double vev){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev)
*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_H_unob_gQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb, double vev){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev)
*K_g(p2out,p2in)/HEJ::C_F;
}
//@}
// Begin finite mass stuff
#ifdef HEJ_BUILD_WITH_QCDLOOP
namespace {
// All the stuff needed for the box functions in qg->qgH now...
COM E1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2=-(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*(1 - 8.*mq*mq/s12 + S2/(2.*s12) +
S2*(s12 - 8.*mq*mq)*(s34 + S1)/(2.*s12*Delta) +
2.*(s34 + S1)*(s34 + S1)/Delta +
S2*pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)*C0DD(k2,
k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2,
mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
2.*(s34 + S1)/Delta +
S2*pow((s34 + S1),2)/Delta/Delta) + (C0DD(k1, q2, mq) -
C0DD(k1 + k2, q2, mq))*(1. - 4.*mq*mq/s12) -
C0DD(k1 + k2, q2, mq)*2.*s34/
S1 - (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*s34*(s34 +
S1)/(S1*Delta) + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*(2.*s34*(s34 +
S1)*(S1 - S2)/(Delta*Sigma) +
2.*s34*(s34 + S1)/(S1*Delta)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2,
mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
S1)*(2.*s12*s34 -
S2*(S1 + S2))/(Delta*Sigma));
}
COM F1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-S2*D0DD(k1, k2, q2,
mq)*(0.5 - (s12 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) -
s12*pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2,
mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
S2*pow((s34 + S2),2)/Delta/Delta)
- (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))*(1. - 4.*mq*mq/s12)
- C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2,
mq))*2.*pow((s34 + S2),2)/((s12 + S1)*Delta) - (B0DD(
q2, mq) - B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*2.*s34*(s34 +
S2)*(S2 - S1)/(Delta*Sigma) + (B0DD(
k1 + k2, mq) -
B0DD(k1 + k2 + q2,
mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
S2)*(2.*s12*s34 -
S2*(S1 + S2))/(Delta*Sigma));
}
COM G1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor*(S2*D0DD(k1, q2, k2,
mq)*(Delta/s12/s12 - 4.*mq*mq/s12) -
S2*((s12 + S1)*C0DD(k1, k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./
s12/s12 - (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
S2*((s12 + S2)*C0DD(k1 + q2, k2, mq) -
S2*C0DD(k2, q2, mq))*(1./
s12/s12 + (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
C0DD(k1, q2, mq) - (C0DD(k1, k2 + q2, mq) -
C0DD(k1, q2, mq))*4.*mq*mq/
s12 + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2./
s12 + (B0DD(k1 + q2, mq) -
B0DD(q2, mq))*2.*s34/(s12*S1) + (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/(s12*(s12 + S1)));
}
COM E4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (-s12*D0DD(k2, k1, q2,
mq)*(0.5 - (S1 - 8.*mq*mq)*(s34 + S1)/(2.*Delta) -
s12*pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)*C0DD(k2,
k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
s12*pow((s34 + S1),2)/Delta/Delta) -
C0DD(k1 + k2, q2, mq) + (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(2.*s34/Delta +
2.*s12*(s34 + S1)/((s12 + S2)*Delta)) - (B0DD(
q2, mq) - B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*((2.*s34*(2.*s12*s34 - S2*(S1 + S2) +
s12*(S1 - S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))
*((2.*s12*(2.*s12*s34 - S1*(S1 + S2) + s34*(S2 - S1)))/(Delta*Sigma)));
}
COM F4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (-s12*D0DD(k1, k2, q2,
mq)*(0.5 + (S1 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) +
s12*pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
s12*pow((s34 + S2),2)/Delta/Delta) -
C0DD(k1 + k2, q2, mq) - (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*(s34 +
S2)/Delta + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*2.*s34*(2.*s12*s34 -
S1*(S1 + S2) +
s12*(S2 - S1))/(Delta*Sigma) - (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))
*(2.*s12*(2.*s12*s34 - S2*(S1 + S2) + s34*(S1 - S2))/(Delta*Sigma)));
}
COM G4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor* (-D0DD(k1, q2, k2,
mq)*(Delta/s12 + (s12 + S1)/2. -
4.*mq*mq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./
s12 - (S1 - 4.*mq*mq)/(2.*Delta)) + ((s12 + S2)*C0DD(
k1 + q2, k2, mq) -
S2*C0DD(k2, q2, mq))*(1./
s12 + (S1 - 4.*mq*mq)/(2.*Delta)) + (B0DD(
k1 + k2 + q2, mq) -
B0DD(k1 + q2, mq))*2./(s12 + S2));
}
COM E10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta +
12.*mq*mq*S1*(s34 + S1)/Delta/Delta -
4.*s12*S1*pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2, mq))*(1./Delta +
4.*mq*mq*S1/Delta/Delta -
4.*s12*S1*pow((s34 + S1),2)/Delta/Delta/Delta) +
C0DD(k1 + k2, q2, mq)*(4.*s12*s34*(S1 - S2)/(Delta*Sigma) -
4.*(s12 -
2.*mq*mq)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma)) + (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(4.*(s34 + S1)/((s12 + S2)*Delta) +
8.*S1*(s34 + S1)/Delta/Delta) + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*(12.*s34*(2.*s12 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) -
4.*s34*(4.*s12 + 3.*S1 +
S2)/(Delta*Sigma) +
8.*s12*s34*(s34*(s12 + S2) -
S1*(s34 +
S1))/(Delta*Delta*Sigma)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2,
mq))*(12.*s12*(2.*s34 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) +
8.*s12*S1*(s34*(s12 + S2) -
S1*(s34 +
S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma));
}
COM F10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (s12*D0DD(k1, k2, q2,
mq)*((s34 + S2)/Delta - 4.*mq*mq/Delta +
12.*mq*mq*s34*(s12 + S1)/Delta/Delta -
4.*s12*pow((s34 + S2),2)/Delta/Delta -
4.*s12*S1*pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2, mq))*(1./Delta +
4.*mq*mq*S1/Delta/Delta -
4.*s12*(s34 + S2)/Delta/Delta -
4.*s12*S1*pow((s34 + S2),2)/Delta/Delta/Delta) -
C0DD(k1 + k2, q2, mq)*(4.*s12*s34/(S2*Delta) +
4.*s12*s34*(S2 - S1)/(Delta*Sigma) +
4.*(s12 -
2.*mq*mq)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma)) - (B0DD(
k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(4.*s34/(S2*Delta) +
8.*s34*(s12 + S1)/Delta/Delta) - (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*(-12*s34*(2*s12 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) -
4.*s12*s34*s34/(S2*Delta*Delta) +
4.*s34*S1/(Delta*Sigma) -
4.*s34*(s12*s34*(2.*s12 + S2) -
S1*S1*(2.*s12 +
S1))/(Delta*Delta*Sigma)) - (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(-12.*s12*(2.*s34 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) +
8.*s12*(2.*s34 + S1)/(Delta*Sigma) -
8.*s12*s34*(2.*s12*s34 - S1*(S1 + S2) +
s12*(S2 -
S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma));
}
COM G10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor* (-D0DD(k1, q2, k2, mq)*(1. +
4.*S1*mq*mq/Delta) + ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./Delta +
4.*S1*mq*mq/Delta/Delta) - ((s12 + S2)*C0DD(k1 + q2,
k2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta +
4.*S1*mq*mq/Delta/Delta) + (B0DD(k1 + k2 + q2, mq) -
B0DD(k1 + q2, mq))*4.*(s34 +
S1)/(Delta*(s12 + S2)) + (B0DD(q2, mq) -
B0DD(k2 + q2, mq))*4.*s34/(Delta*S2));
}
COM H1(HLV k1, HLV k2, HLV kh, double mq){
return E1(k1,k2,kh,mq)+F1(k1,k2,kh,mq)+G1(k1,k2,kh,mq);
}
COM H4(HLV k1, HLV k2, HLV kh, double mq){
return E4(k1,k2,kh,mq)+F4(k1,k2,kh,mq)+G4(k1,k2,kh,mq);
}
COM H10(HLV k1, HLV k2, HLV kh, double mq){
return E10(k1,k2,kh,mq)+F10(k1,k2,kh,mq)+G10(k1,k2,kh,mq);
}
COM H2(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H1(k2,k1,kh,mq);
}
COM H5(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H4(k2,k1,kh,mq);
}
COM H12(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H10(k2,k1,kh,mq);
}
// FL and FT functions
COM FL(HLV q1, HLV q2, double mq){
HLV Q = q1 + q2;
double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
return -1./(2.*detQ2)*((2.-
3.*q1.m2()*q2.dot(Q)/detQ2)*(B0DD(q1, mq) -
B0DD(Q, mq)) + (2. -
3.*q2.m2()*q1.dot(Q)/detQ2)*(B0DD(q2, mq) -
B0DD(Q, mq)) - (4.*mq*mq + q1.m2() + q2.m2() +
Q.m2() - 3.*q1.m2()*q2.m2()*Q.m2()/detQ2)*C0DD(
q1, q2, mq) - 2.);
}
COM FT(HLV q1, HLV q2, double mq){
HLV Q = q1 + q2;
double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
return -1./(2.*detQ2)*(Q.m2()*(B0DD(q1, mq) + B0DD(q2, mq) - 2.*B0DD(Q, mq) -
2.*q1.dot(q2)*C0DD(q1, q2, mq)) + (q1.m2() -
q2.m2()) *(B0DD(q1, mq) - B0DD(q2, mq))) -
q1.dot(q2)*FL(q1, q2, mq);
}
HLV ParityFlip(HLV p){
HLV flippedVector;
flippedVector.setE(p.e());
flippedVector.setX(-p.x());
flippedVector.setY(-p.y());
flippedVector.setZ(-p.z());
return flippedVector;
}
/// @brief HC amp for qg->qgH with finite top (i.e. j^{++}_H)
void g_gH_HC(HLV pa, HLV p1,
HLV pH, double mq, double vev, current &retAns)
{
current cura1,pacur,p1cur,pHcur,conjeps1,conjepsH1,epsa,epsHa,epsHapart1,
epsHapart2,conjepsH1part1,conjepsH1part2;
COM ang1a,sqa1;
const double F = 4.*mq*mq/vev;
// Easier to have the whole thing as current object so I can use cdot functionality.
// Means I need to write pa,p1 as current objects
to_current(pa, pacur);
to_current(p1,p1cur);
to_current(pH,pHcur);
bool gluonforward = true;
if(pa.z() < 0)
gluonforward = false;
//HEJ gauge
jio(pa,false,p1,false,cura1);
if(gluonforward){
// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
// sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp)
sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
} else {
ang1a = sqrt(pa.minus()*p1.plus());
sqa1 = sqrt(pa.minus()*p1.plus());
}
const double prop = (pa-p1-pH).m2();
cmult(-1./sqrt(2)/ang1a,cura1,conjeps1);
cmult(1./sqrt(2)/sqa1,cura1,epsa);
const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
const COM h4 = H4(p1,-pa,pH,mq);
const COM h5 = H5(p1,-pa,pH,mq);
const COM h10 = H10(p1,-pa,pH,mq);
const COM h12 = H12(p1,-pa,pH,mq);
cmult(Fta*pa.dot(pH), epsa, epsHapart1);
cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2);
cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
cadd(epsHapart1, epsHapart2, epsHa);
cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
const COM aH1 = cdot(pHcur, cura1);
current T1,T2,T3,T4,T5,T6,T7,T8,T9,T10;
if(gluonforward){
cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1);
cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/ang1a, epsHa, T2);
}
else{
cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())
*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/sqa1, conjepsH1, T1);
cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())
*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/ang1a, epsHa, T2);
}
cmult(sqrt(2.)/ang1a*aH1, epsHa, T3);
cmult(sqrt(2.)/sqa1*aH1, conjepsH1, T4);
cmult(-sqrt(2.)*Fta*pa.dot(p1)*aH1/sqa1, conjeps1, T5);
cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/ang1a, epsa, T6);
cmult(-aH1/sqrt(2.)/sqa1*h4*8.*COM(0.,1.)*M_PI*M_PI, conjeps1, T7);
cmult(aH1/sqrt(2.)/ang1a*h5*8.*COM(0.,1.)*M_PI*M_PI, epsa, T8);
cmult(aH1*aH1/2./ang1a/sqa1*h10*8.*COM(0.,1.)*M_PI*M_PI, pacur, T9);
cmult(-aH1*aH1/2./ang1a/sqa1*h12*8.*COM(0.,1.)*M_PI*M_PI, p1cur, T10);
current ans;
- for(int i=0;i<4;i++)
+ for(int i=0;i<4;++i)
{
ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5[i]+T6[i]+T7[i]+T8[i]+T9[i]+T10[i];
}
retAns[0] = F/prop*ans[0];
retAns[1] = F/prop*ans[1];
retAns[2] = F/prop*ans[2];
retAns[3] = F/prop*ans[3];
}
/// @brief HNC amp for qg->qgH with finite top (i.e. j^{+-}_H)
void g_gH_HNC(HLV pa, HLV p1, HLV pH, double mq, double vev, current &retAns)
{
const double F = 4.*mq*mq/vev;
COM ang1a,sqa1;
current conjepsH1,epsHa,p1cur,pacur,pHcur,conjeps1,epsa,paplusp1cur,
p1minuspacur,cur1a,cura1,epsHapart1,epsHapart2,conjepsH1part1,
conjepsH1part2;
// Find here if pa, meaning the gluon, is forward or backward
bool gluonforward = true;
if(pa.z() < 0)
gluonforward = false;
jio(pa,true,p1,true,cura1);
joi(p1,true,pa,true,cur1a);
to_current(pa,pacur);
to_current(p1,p1cur);
to_current(pH,pHcur);
to_current(pa+p1,paplusp1cur);
to_current(p1-pa,p1minuspacur);
const COM aH1 = cdot(pHcur,cura1);
const COM oneHa = std::conj(aH1); // = cdot(pHcur,cur1a)
if(gluonforward){
// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
// sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp)
sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
}
else {
ang1a = sqrt(pa.minus()*p1.plus());
sqa1 = sqrt(pa.minus()*p1.plus());
}
const double prop = (pa-p1-pH).m2();
cmult(1./sqrt(2)/sqa1, cur1a, epsa);
cmult(-1./sqrt(2)/sqa1, cura1, conjeps1);
const COM phase = cdot(conjeps1, epsa);
const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
const COM Falpha = FT(p1-pa,pa-p1-pH,mq);
const COM Fbeta = FL(p1-pa,pa-p1-pH,mq);
const COM h1 = H1(p1,-pa, pH, mq);
const COM h2 = H2(p1,-pa, pH, mq);
const COM h4 = H4(p1,-pa, pH, mq);
const COM h5 = H5(p1,-pa, pH, mq);
const COM h10 = H10(p1,-pa, pH, mq);
const COM h12 = H12(p1,-pa, pH, mq);
cmult(Fta*pa.dot(pH), epsa, epsHapart1);
cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2);
cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
cadd(epsHapart1, epsHapart2, epsHa);
cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
current T1,T2,T3,T4,T5a,T5b,T6,T7,T8a,T8b,T9,T10,T11a,
T11b,T12a,T12b,T13;
if(gluonforward){
cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1);
cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/sqa1, epsHa, T2);
}
else{
cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())
*prop/sqa1, conjepsH1, T1);
cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())*((p1.x()+COM(0.,1.)*p1.y())/p1.perp())
*prop/sqa1, epsHa, T2);
}
const COM boxdiagFact = 8.*COM(0.,1.)*M_PI*M_PI;
cmult(aH1*sqrt(2.)/sqa1, epsHa, T3);
cmult(oneHa*sqrt(2.)/sqa1, conjepsH1, T4);
cmult(-2.*phase*Fta*pa.dot(pH), p1cur, T5a);
cmult(2.*phase*Ft1*p1.dot(pH), pacur, T5b);
cmult(-sqrt(2.)*Fta*p1.dot(pa)*oneHa/sqa1, conjeps1, T6);
cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/sqa1, epsa, T7);
cmult(-boxdiagFact*phase*h2, pacur, T8a);
cmult(boxdiagFact*phase*h1, p1cur, T8b);
cmult(boxdiagFact*aH1/sqrt(2.)/sqa1*h5, epsa, T9);
cmult(-boxdiagFact*oneHa/sqrt(2.)/sqa1*h4, conjeps1, T10);
cmult(boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h10, pacur, T11a);
cmult(-boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h12, p1cur, T11b);
cmult(-phase/(pa-p1).m2()*Falpha*(p1-pa).dot(pa-p1-pH), paplusp1cur, T12a);
cmult(phase/(pa-p1).m2()*Falpha*(pa+p1).dot(pa-p1-pH), p1minuspacur, T12b);
cmult(-phase*Fbeta*(pa-p1-pH).m2(), paplusp1cur, T13);
current ans;
- for(int i=0;i<4;i++)
+ for(int i=0;i<4;++i)
{
ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5a[i]+T5b[i]+T6[i]+T7[i]+T8a[i]+T8b[i]
+T9[i]+T10[i]+T11a[i]+T11b[i]+T12a[i]+T12b[i]+T13[i];
}
retAns[0] = F/prop*ans[0];
retAns[1] = F/prop*ans[1];
retAns[2] = F/prop*ans[2];
retAns[3] = F/prop*ans[3];
}
} // namespace anonymous
// JDC - new amplitude with Higgs emitted close to gluon with full mt effects.
// Keep usual HEJ-style function call
double ME_Houtside_gq(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pH,
double mq, bool includeBottom, double mq2, double vev
){
current cur2bplus,cur2bminus, cur2bplusFlip, cur2bminusFlip;
current retAns,retAnsb;
joi(p2out,true,p2in,true,cur2bplus);
joi(p2out,false,p2in,false,cur2bminus);
joi(ParityFlip(p2out),true,ParityFlip(p2in),true,cur2bplusFlip);
joi(ParityFlip(p2out),false,ParityFlip(p2in),false,cur2bminusFlip);
COM app1,app2,apm1,apm2;
COM app3, app4, apm3, apm4;
if(!includeBottom)
{
g_gH_HC(p1in,p1out,pH,mq,vev,retAns);
app1=cdot(retAns,cur2bplus);
app2=cdot(retAns,cur2bminus);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
app3=cdot(retAns,cur2bplusFlip);
app4=cdot(retAns,cur2bminusFlip);
// And non-conserving bits
g_gH_HNC(p1in,p1out,pH,mq,vev,retAns);
apm1=cdot(retAns,cur2bplus);
apm2=cdot(retAns,cur2bminus);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
apm3=cdot(retAns,cur2bplusFlip);
apm4=cdot(retAns,cur2bminusFlip);
} else {
g_gH_HC(p1in,p1out,pH,mq,vev,retAns);
g_gH_HC(p1in,p1out,pH,mq2,vev,retAnsb);
app1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
app2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,vev,retAnsb);
app3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
app4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
// And non-conserving bits
g_gH_HNC(p1in,p1out,pH,mq,vev,retAns);
g_gH_HNC(p1in,p1out,pH,mq2,vev,retAnsb);
apm1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
apm2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,vev,retAnsb);
apm3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
apm4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
}
return abs2(app1) + abs2(app2) + abs2(app3) + abs2(app4) + abs2(apm1)
+ abs2(apm2) + abs2(apm3) + abs2(apm4);
}
#endif // HEJ_BUILD_WITH_QCDLOOP
double C2gHgm(HLV p2, HLV p1, HLV pH, double vev)
{
const double A=1./(3.*M_PI*vev);
// Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta
double s12,p1p,p2p;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
p1p=p1.plus();
p2p=p2.plus();
} else { // opposite case
p1p=p1.minus();
p2p=p2.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=COM(0,1)*A/(2.*s12)*(p2p/p1p*conj(p1perp)*p3perp+p1p/p2p*p1perp*conj(p3perp));
temp=temp*conj(temp);
return temp.real();
}
double C2gHgp(HLV p2, HLV p1, HLV pH, double vev)
{
const double A=1./(3.*M_PI*vev);
// Implements Eq. (4.23) in hep-ph/0301013
double s12,php,p1p,phm;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
php=pH.plus();
phm=pH.minus();
p1p=p1.plus();
} else { // opposite case
php=pH.minus();
phm=pH.plus();
p1p=p1.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=-COM(0,1)*A/(2.*s12)*( conj(p1perp*p3perp)*pow(php/p1p,2)/(1.+php/p1p)
+s12*(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1p*phm)
-pow(conj(p3perp)
+(1.+php/p1p)*conj(p1perp),2)/((1.+php/p1p)*(pH.m2()+2.*p1.dot(pH)))) );
temp=temp*conj(temp);
return temp.real();
}
double C2qHqm(HLV p2, HLV p1, HLV pH, double vev)
{
const double A=1./(3.*M_PI*vev);
// Implements Eq. (4.22) in hep-ph/0301013
double s12,p2p,p1p;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
p2p=p2.plus();
p1p=p1.plus();
} else { // opposite case
p2p=p2.minus();
p1p=p1.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=A/(2.*s12)*( sqrt(p2p/p1p)*p3perp*conj(p1perp)
+sqrt(p1p/p2p)*p1perp*conj(p3perp) );
temp=temp*conj(temp);
return temp.real();
}
diff --git a/src/Tensor.cc b/src/Tensor.cc
index 9cf5c0f..69d8d75 100644
--- a/src/Tensor.cc
+++ b/src/Tensor.cc
@@ -1,691 +1,691 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Tensor.hh"
#include <array>
#include <iostream>
#include <valarray>
#include <vector>
#include <CLHEP/Vector/LorentzVector.h>
namespace HEJ {
namespace{
// Tensor Template definitions
short int sigma_index5[1024];
short int sigma_index3[64];
std::valarray<COM> permfactor5;
std::valarray<COM> permfactor3;
short int helfactor5[2][1024];
short int helfactor3[2][64];
// map from a list of tensor lorentz indices onto a single index
// in d dimensional spacetime
// TODO: basically the same as detail::accumulate_idx
template<std::size_t N>
std::size_t tensor_to_list_idx(std::array<int, N> const & indices) {
std::size_t list_idx = 0;
for(std::size_t i = 1, factor = 1; i <= N; ++i) {
list_idx += factor*indices[N-i];
factor *= detail::d;
}
#ifndef NDEBUG
constexpr std::size_t max_possible_idx = detail::power(detail::d, N);
assert(list_idx < max_possible_idx);
#endif
return list_idx;
}
// generate all unique perms of vectors {a,a,a,a,b}, return in perms
// set_permfactor is a bool which encodes the anticommutation relations of the
// sigma matrices namely if we have sigma0, set_permfactor= false because it
// commutes with all others otherwise we need to assign a minus sign to odd
// permutations, set in permfactor
// note, inital perm is always even
void perms41(int same4, int diff, std::vector<std::array<int,5>> & perms){
bool set_permfactor(true);
if(same4==0||diff==0)
set_permfactor=false;
- for(int diffpos=0;diffpos<5;diffpos++){
+ for(int diffpos=0;diffpos<5;++diffpos){
std::array<int,5> tempvec={same4,same4,same4,same4,same4};
tempvec[diffpos]=diff;
perms.push_back(tempvec);
if(set_permfactor){
if(diffpos%2==1)
permfactor5[tensor_to_list_idx(tempvec)]=-1.;
}
}
}
// generate all unique perms of vectors {a,a,a,b,b}, return in perms
// note, inital perm is always even
void perms32(int same3, int diff, std::vector<std::array<int,5>> & perms){
bool set_permfactor(true);
if(same3==0||diff==0)
set_permfactor=false;
- for(int diffpos1=0;diffpos1<5;diffpos1++){
- for(int diffpos2=diffpos1+1;diffpos2<5;diffpos2++){
+ for(int diffpos1=0;diffpos1<5;++diffpos1){
+ for(int diffpos2=diffpos1+1;diffpos2<5;++diffpos2){
std::array<int,5> tempvec={same3,same3,same3,same3,same3};
tempvec[diffpos1]=diff;
tempvec[diffpos2]=diff;
perms.push_back(tempvec);
if(set_permfactor){
if((diffpos2-diffpos1)%2==0)
permfactor5[tensor_to_list_idx(tempvec)]=-1.;
}
}
}
}
// generate all unique perms of vectors {a,a,a,b,c}, return in perms
// we have two bools since there are three distinct type of sigma matrices to
// permute, so need to test if the 3xrepeated = sigma0 or if one of the
// singles is sigma0
// if diff1/diff2 are distinct, flipping them results in a change of perm,
// otherwise it'll be symmetric under flip -> encode this in set_permfactor2
// as long as diff2!=0 can ensure inital perm is even
// if diff2=0 then it is not possible to ensure inital perm even -> encode in
// bool signflip
void perms311(int same3, int diff1, int diff2,
std::vector<std::array<int,5>> & perms
){
bool set_permfactor2(true);
bool same0(false);
bool diff0(false);
bool signflip(false); // if true, inital perm is odd
if(same3==0) // is the repeated matrix sigma0?
same0 = true;
else if(diff2==0){ // is one of the single matrices sigma0
diff0=true;
if((diff1%3-same3)!=-1)
signflip=true;
} else if(diff1==0){
std::cerr<<"Note, only first and last argument may be zero"<<std::endl;
return;
}
// possible outcomes: tt, ft, tf
- for(int diffpos1=0;diffpos1<5;diffpos1++){
- for(int diffpos2=diffpos1+1;diffpos2<5;diffpos2++){
+ for(int diffpos1=0;diffpos1<5;++diffpos1){
+ for(int diffpos2=diffpos1+1;diffpos2<5;++diffpos2){
std::array<int,5> tempvec={same3,same3,same3,same3,same3};
tempvec[diffpos1]=diff1;
tempvec[diffpos2]=diff2;
perms.push_back(tempvec);
if(!same0 && !diff0){
// full antisymmetric under exchange of same3,diff1,diff2
if((diffpos2-diffpos1)%2==0){
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1); // odd perm
// if this perm is odd, swapping diff1<->diff2 automatically even
set_permfactor2=false;
} else {
permfactor5[tensor_to_list_idx(tempvec)]=COM(0,1); // even perm
// if this perm is even, swapping diff1<->diff2 automatically odd
set_permfactor2=true;
}
} else if(diff0){// one of the single matrices is sigma0
if(signflip){ // initial config is odd!
if(diffpos1%2==1){
permfactor5[tensor_to_list_idx(tempvec)]=COM(0,1); // even perm
// initally symmetric under diff1<->diff2 =>
// if this perm is even, automatically even for first diffpos2
set_permfactor2=false;
} else {
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1); // odd perm
// initally symmetric under diff1<->diff2 =>
// if this perm is odd, automatically odd for first diffpos2
set_permfactor2=true;
}
} else {// diff0=true, initial config is even
if(diffpos1%2==1){
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1); // odd perm
// initally symmetric under diff1<->diff2 =>
// if this perm is odd, automatically odd for first diffpos2
set_permfactor2=true;
} else {
permfactor5[tensor_to_list_idx(tempvec)]=COM(0,1); // even perm
// initally symmetric under diff1<->diff2 =>
// if this perm is even, automatically even for first diffpos2
set_permfactor2=false;
}
}
if((diffpos2-diffpos1-1)%2==1)
set_permfactor2=!set_permfactor2; // change to account for diffpos2
} else if(same0){
// the repeated matrix is sigma0
// => only relative positions of diff1, diff2 matter.
// always even before flip because diff1pos<diff2pos
permfactor5[tensor_to_list_idx(tempvec)]=COM(0,1);
// if this perm is odd, swapping diff1<->diff2 automatically odd
set_permfactor2=true;
}
tempvec[diffpos1]=diff2;
tempvec[diffpos2]=diff1;
perms.push_back(tempvec);
if(set_permfactor2)
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1);
else
permfactor5[tensor_to_list_idx(tempvec)]=COM(0,1);
}
}
} // end perms311
// generate all unique perms of vectors {a,a,b,b,c}, return in perms
void perms221(int same2a, int same2b, int diff,
std::vector<std::array<int,5>> & perms
){
bool set_permfactor1(true);
bool set_permfactor2(true);
if(same2b==0){
std::cerr<<"second entry in perms221() shouldn't be zero" <<std::endl;
return;
} else if(same2a==0)
set_permfactor1=false;
else if(diff==0)
set_permfactor2=false;
- for(int samepos=0;samepos<5;samepos++){
+ for(int samepos=0;samepos<5;++samepos){
int permcount = 0;
- for(int samepos2=samepos+1;samepos2<5;samepos2++){
- for(int diffpos=0;diffpos<5;diffpos++){
+ for(int samepos2=samepos+1;samepos2<5;++samepos2){
+ for(int diffpos=0;diffpos<5;++diffpos){
if(diffpos==samepos||diffpos==samepos2) continue;
std::array<int,5> tempvec={same2a,same2a,same2a,same2a,same2a};
tempvec[samepos]=same2b;
tempvec[samepos2]=same2b;
tempvec[diffpos]=diff;
perms.push_back(tempvec);
if(set_permfactor1){
if(set_permfactor2){// full anti-symmetry
if(permcount%2==1)
permfactor5[tensor_to_list_idx(tempvec)]=-1.;
} else { // diff is sigma0
if( ((samepos2-samepos-1)%3>0)
&& (abs(abs(samepos2-diffpos)-abs(samepos-diffpos))%3>0) )
permfactor5[tensor_to_list_idx(tempvec)]=-1.;
}
} else { // same2a is sigma0
if((diffpos>samepos) && (diffpos<samepos2))
permfactor5[tensor_to_list_idx(tempvec)]=-1.;
}
- permcount++;
+ ++permcount;
}
}
}
}
// generate all unique perms of vectors {a,a,b,b,c}, return in perms
// there must be a sigma zero if we have 4 different matrices in string
// bool is true if sigma0 is the repeated matrix
void perms2111(int same2, int diff1,int diff2,int diff3,
std::vector<std::array<int,5>> & perms
){
bool twozero(false);
if(same2==0)
twozero=true;
else if (diff1!=0){
std::cerr<<"One of first or second argurments must be a zero"<<std::endl;
return;
} else if(diff2==0|| diff3==0){
std::cerr<<"Only first and second arguments may be a zero."<<std::endl;
return;
}
int permcount = 0;
- for(int diffpos1=0;diffpos1<5;diffpos1++){
- for(int diffpos2=0;diffpos2<5;diffpos2++){
+ for(int diffpos1=0;diffpos1<5;++diffpos1){
+ for(int diffpos2=0;diffpos2<5;++diffpos2){
if(diffpos2==diffpos1) continue;
- for(int diffpos3=0;diffpos3<5;diffpos3++){
+ for(int diffpos3=0;diffpos3<5;++diffpos3){
if(diffpos3==diffpos2||diffpos3==diffpos1) continue;
std::array<int,5> tempvec={same2,same2,same2,same2,same2};
tempvec[diffpos1]=diff1;
tempvec[diffpos2]=diff2;
tempvec[diffpos3]=diff3;
perms.push_back(tempvec);
if(twozero){// don't care about exact positions of singles, just order
if(diffpos2>diffpos3 && diffpos3>diffpos1)
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1);// odd
else if(diffpos1>diffpos2 && diffpos2>diffpos3)
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1);// odd
else if(diffpos3>diffpos1 && diffpos1>diffpos2)
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1);// odd
else
permfactor5[tensor_to_list_idx(tempvec)]=COM(0,1);// evwn
} else {
if(permcount%2==1)
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1);
else
permfactor5[tensor_to_list_idx(tempvec)]=COM(0,1);
}
- permcount++;
+ ++permcount;
}
}
}
}
void perms21(int same, int diff, std::vector<std::array<int,3>> & perms){
bool set_permfactor(true);
if(same==0||diff==0)
set_permfactor=false;
- for(int diffpos=0; diffpos<3;diffpos++){
+ for(int diffpos=0; diffpos<3;++diffpos){
std::array<int,3> tempvec={same,same,same};
tempvec[diffpos]=diff;
perms.push_back(tempvec);
if(set_permfactor && diffpos==1)
permfactor3[tensor_to_list_idx(tempvec)]=-1.;
}
}
void perms111(int diff1, int diff2, int diff3,
std::vector<std::array<int,3>> & perms
){
bool sig_zero(false);
if(diff1==0)
sig_zero=true;
else if(diff2==0||diff3==0){
std::cerr<<"Only first argument may be a zero."<<std::endl;
return;
}
int permcount=0;
- for(int pos1=0;pos1<3;pos1++){
- for(int pos2=pos1+1;pos2<3;pos2++){
+ for(int pos1=0;pos1<3;++pos1){
+ for(int pos2=pos1+1;pos2<3;++pos2){
std::array<int,3> tempvec={diff1,diff1,diff1};
tempvec[pos1]=diff2;
tempvec[pos2]=diff3;
perms.push_back(tempvec);
if(sig_zero){
permfactor3[tensor_to_list_idx(tempvec)]=1.*COM(0,1); // even
} else if(permcount%2==1){
permfactor3[tensor_to_list_idx(tempvec)]=-1.*COM(0,1); // odd
} else {
permfactor3[tensor_to_list_idx(tempvec)]=1.*COM(0,1); // even
}
tempvec[pos1]=diff3;
tempvec[pos2]=diff2;
perms.push_back(tempvec);
if(sig_zero){
permfactor3[tensor_to_list_idx(tempvec)]=-1.*COM(0,1); // odd
} else if(permcount%2==1){
permfactor3[tensor_to_list_idx(tempvec)]=1.*COM(0,1); // even
} else {
permfactor3[tensor_to_list_idx(tempvec)]=-1.*COM(0,1); // odd
}
- permcount++;
+ ++permcount;
}
}
}
COM perp(CLHEP::HepLorentzVector const & p) {
return COM{p.x(), p.y()};
}
COM perp_hat(CLHEP::HepLorentzVector const & p) {
const COM p_perp = perp(p);
if(p_perp == COM{0., 0.}) return -1.;
return p_perp/std::abs(p_perp);
}
} // anonymous namespace
// "angle" product angle(pi, pj) = <i j>
// see eq. (53) (\eqref{eq:angle_product}) in developer manual
COM angle(CLHEP::HepLorentzVector const & pi, CLHEP::HepLorentzVector const & pj) {
return
+ std::sqrt(COM{pi.minus()*pj.plus()})*perp_hat(pi)
- std::sqrt(COM{pi.plus()*pj.minus()})*perp_hat(pj);
}
// "square" product square(pi, pj) = [i j]
// see eq. (54) (\eqref{eq:square_product}) in developer manual
COM square(CLHEP::HepLorentzVector const & pi, CLHEP::HepLorentzVector const & pj) {
return -std::conj(angle(pi, pj));
}
Tensor<2> metric(){
static const Tensor<2> g = [](){
Tensor<2> g(0.);
g(0,0) = 1.;
g(1,1) = -1.;
g(2,2) = -1.;
g(3,3) = -1.;
return g;
}();
return g;
}
// <1|mu|2>
// see eqs. (48), (49) in developer manual
Tensor<1> current(
CLHEP::HepLorentzVector const & p,
CLHEP::HepLorentzVector const & q,
Helicity h
){
using namespace helicity;
const COM p_perp_hat = perp_hat(p);
const COM q_perp_hat = perp_hat(q);
std::array<std::array<COM, 2>, 2> sqrt_pq;
sqrt_pq[minus][minus] = std::sqrt(COM{p.minus()*q.minus()})*p_perp_hat*conj(q_perp_hat);
sqrt_pq[minus][plus] = std::sqrt(COM{p.minus()*q.plus()})*p_perp_hat;
sqrt_pq[plus][minus] = std::sqrt(COM{p.plus()*q.minus()})*conj(q_perp_hat);
sqrt_pq[plus][plus] = std::sqrt(COM{p.plus()*q.plus()});
Tensor<1> j;
j(0) = sqrt_pq[plus][plus] + sqrt_pq[minus][minus];
j(1) = sqrt_pq[plus][minus] + sqrt_pq[minus][plus];
j(2) = COM{0., 1.}*(sqrt_pq[plus][minus] - sqrt_pq[minus][plus]);
j(3) = sqrt_pq[plus][plus] - sqrt_pq[minus][minus];
return (h == minus)?j:j.complex_conj();
}
// <1|mu nu sigma|2>
// TODO: how does this work?
Tensor<3> rank3_current(
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
Helicity h
){
const CLHEP::HepLorentzVector p1new{ p1.e()<0?-p1:p1 };
const CLHEP::HepLorentzVector p2new{ p2.e()<0?-p2:p2 };
auto j = HEJ::current(p1new, p2new, h);
if(h == helicity::minus) {
j *= -1;
j(0) *= -1;
}
Tensor<3> j3;
for(std::size_t tensor_idx=0; tensor_idx < j3.size(); ++tensor_idx){
const double hfact = helfactor3[h][tensor_idx];
j3.components[tensor_idx] = j(sigma_index3[tensor_idx]) * hfact
* permfactor3[tensor_idx];
}
return j3;
}
// <1|mu nu sigma tau rho|2>
Tensor<5> rank5_current(
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
Helicity h
){
const CLHEP::HepLorentzVector p1new{ p1.e()<0?-p1:p1 };
const CLHEP::HepLorentzVector p2new{ p2.e()<0?-p2:p2 };
auto j = HEJ::current(p1new, p2new, h);
if(h == helicity::minus) {
j *= -1;
j(0) *= -1;
}
Tensor<5> j5;
for(std::size_t tensor_idx=0; tensor_idx < j5.size(); ++tensor_idx){
const double hfact = helfactor5[h][tensor_idx];
j5.components[tensor_idx] = j(sigma_index5[tensor_idx]) * hfact
* permfactor5[tensor_idx];
}
return j5;
}
Tensor<1> to_tensor(CLHEP::HepLorentzVector const & p){
Tensor<1> newT;
newT.components={p.e(),p.x(),p.y(),p.z()};
return newT;
}
// polarisation vector according to eq. (55) in developer manual
Tensor<1> eps(
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pr,
Helicity pol
){
if(pol == helicity::plus) {
return current(pr, pg, pol)/(std::sqrt(2)*square(pg, pr));
}
return current(pr, pg, pol)/(std::sqrt(2)*angle(pg, pr));
}
namespace {
// slow function! - but only need to evaluate once.
bool init_sigma_index_helper(){
// initialize permfactor(3) to list of ones (change to minus one for each odd
// permutation and multiply by i for all permutations in perms2111, perms311,
// perms111)
permfactor5.resize(1024,1.);
permfactor3.resize(64,1.);
// first set sigma_index (5)
std::vector<std::array<int,5>> sigma0indices;
std::vector<std::array<int,5>> sigma1indices;
std::vector<std::array<int,5>> sigma2indices;
std::vector<std::array<int,5>> sigma3indices;
// need to generate all possible permuations of {i,j,k,l,m}
// where each index can be {0,1,2,3,4}
// 1024 possibilities
// perms with 5 same
sigma0indices.push_back({0,0,0,0,0});
sigma1indices.push_back({1,1,1,1,1});
sigma2indices.push_back({2,2,2,2,2});
sigma3indices.push_back({3,3,3,3,3});
// perms with 4 same
perms41(1,0,sigma0indices);
perms41(2,0,sigma0indices);
perms41(3,0,sigma0indices);
perms41(0,1,sigma1indices);
perms41(2,1,sigma1indices);
perms41(3,1,sigma1indices);
perms41(0,2,sigma2indices);
perms41(1,2,sigma2indices);
perms41(3,2,sigma2indices);
perms41(0,3,sigma3indices);
perms41(1,3,sigma3indices);
perms41(2,3,sigma3indices);
// perms with 3 same, 2 same
perms32(0,1,sigma0indices);
perms32(0,2,sigma0indices);
perms32(0,3,sigma0indices);
perms32(1,0,sigma1indices);
perms32(1,2,sigma1indices);
perms32(1,3,sigma1indices);
perms32(2,0,sigma2indices);
perms32(2,1,sigma2indices);
perms32(2,3,sigma2indices);
perms32(3,0,sigma3indices);
perms32(3,1,sigma3indices);
perms32(3,2,sigma3indices);
// perms with 3 same, 2 different
perms311(1,2,3,sigma0indices);
perms311(2,3,1,sigma0indices);
perms311(3,1,2,sigma0indices);
perms311(0,2,3,sigma1indices);
perms311(2,3,0,sigma1indices);
perms311(3,2,0,sigma1indices);
perms311(0,3,1,sigma2indices);
perms311(1,3,0,sigma2indices);
perms311(3,1,0,sigma2indices);
perms311(0,1,2,sigma3indices);
perms311(1,2,0,sigma3indices);
perms311(2,1,0,sigma3indices);
perms221(1,2,0,sigma0indices);
perms221(1,3,0,sigma0indices);
perms221(2,3,0,sigma0indices);
perms221(0,2,1,sigma1indices);
perms221(0,3,1,sigma1indices);
perms221(2,3,1,sigma1indices);
perms221(0,1,2,sigma2indices);
perms221(0,3,2,sigma2indices);
perms221(1,3,2,sigma2indices);
perms221(0,1,3,sigma3indices);
perms221(0,2,3,sigma3indices);
perms221(1,2,3,sigma3indices);
perms2111(0,1,2,3,sigma0indices);
perms2111(1,0,2,3,sigma1indices);
perms2111(2,0,3,1,sigma2indices);
perms2111(3,0,1,2,sigma3indices);
if(sigma0indices.size()!=256){
std::cerr<<"sigma_index not set: ";
std::cerr<<"sigma0indices has "<< sigma0indices.size() << " components" << std::endl;
return false;
} else if(sigma1indices.size()!=256){
std::cerr<<"sigma_index not set: ";
std::cerr<<"sigma1indices has "<< sigma0indices.size() << " components" << std::endl;
return false;
} else if(sigma2indices.size()!=256){
std::cerr<<"sigma_index not set: ";
std::cerr<<"sigma2indices has "<< sigma0indices.size() << " components" << std::endl;
return false;
} else if(sigma3indices.size()!=256){
std::cerr<<"sigma_index not set: ";
std::cerr<<"sigma3indices has "<< sigma0indices.size() << " components" << std::endl;
return false;
}
- for(int i=0;i<256;i++){
+ for(int i=0;i<256;++i){
// map each unique set of tensor indices to its position in a list
int index0 = tensor_to_list_idx(sigma0indices.at(i));
int index1 = tensor_to_list_idx(sigma1indices.at(i));
int index2 = tensor_to_list_idx(sigma2indices.at(i));
int index3 = tensor_to_list_idx(sigma3indices.at(i));
sigma_index5[index0]=0;
sigma_index5[index1]=1;
sigma_index5[index2]=2;
sigma_index5[index3]=3;
short int sign[4]={1,-1,-1,-1};
// plus->true->1
helfactor5[1][index0] = sign[sigma0indices.at(i)[1]]
* sign[sigma0indices.at(i)[3]];
helfactor5[1][index1] = sign[sigma1indices.at(i)[1]]
* sign[sigma1indices.at(i)[3]];
helfactor5[1][index2] = sign[sigma2indices.at(i)[1]]
* sign[sigma2indices.at(i)[3]];
helfactor5[1][index3] = sign[sigma3indices.at(i)[1]]
* sign[sigma3indices.at(i)[3]];
// minus->false->0
helfactor5[0][index0] = sign[sigma0indices.at(i)[0]]
* sign[sigma0indices.at(i)[2]]
* sign[sigma0indices.at(i)[4]];
helfactor5[0][index1] = sign[sigma1indices.at(i)[0]]
* sign[sigma1indices.at(i)[2]]
* sign[sigma1indices.at(i)[4]];
helfactor5[0][index2] = sign[sigma2indices.at(i)[0]]
* sign[sigma2indices.at(i)[2]]
* sign[sigma2indices.at(i)[4]];
helfactor5[0][index3] = sign[sigma3indices.at(i)[0]]
* sign[sigma3indices.at(i)[2]]
* sign[sigma3indices.at(i)[4]];
}
// now set sigma_index3
std::vector<std::array<int,3>> sigma0indices3;
std::vector<std::array<int,3>> sigma1indices3;
std::vector<std::array<int,3>> sigma2indices3;
std::vector<std::array<int,3>> sigma3indices3;
// perms with 3 same
sigma0indices3.push_back({0,0,0});
sigma1indices3.push_back({1,1,1});
sigma2indices3.push_back({2,2,2});
sigma3indices3.push_back({3,3,3});
// 2 same
perms21(1,0,sigma0indices3);
perms21(2,0,sigma0indices3);
perms21(3,0,sigma0indices3);
perms21(0,1,sigma1indices3);
perms21(2,1,sigma1indices3);
perms21(3,1,sigma1indices3);
perms21(0,2,sigma2indices3);
perms21(1,2,sigma2indices3);
perms21(3,2,sigma2indices3);
perms21(0,3,sigma3indices3);
perms21(1,3,sigma3indices3);
perms21(2,3,sigma3indices3);
// none same
perms111(1,2,3,sigma0indices3);
perms111(0,2,3,sigma1indices3);
perms111(0,3,1,sigma2indices3);
perms111(0,1,2,sigma3indices3);
if(sigma0indices3.size()!=16){
std::cerr<<"sigma_index3 not set: ";
std::cerr<<"sigma0indices3 has "<< sigma0indices3.size() << " components" << std::endl;
return false;
} else if(sigma1indices3.size()!=16){
std::cerr<<"sigma_index3 not set: ";
std::cerr<<"sigma1indices3 has "<< sigma0indices3.size() << " components" << std::endl;
return false;
} else if(sigma2indices3.size()!=16){
std::cerr<<"sigma_index3 not set: ";
std::cerr<<"sigma2indices3 has "<< sigma0indices3.size() << " components" << std::endl;
return false;
} else if(sigma3indices3.size()!=16){
std::cerr<<"sigma_index3 not set: ";
std::cerr<<"sigma3indices3 has "<< sigma0indices3.size() << " components" << std::endl;
return false;
}
- for(int i=0;i<16;i++){
+ for(int i=0;i<16;++i){
int index0 = tensor_to_list_idx(sigma0indices3.at(i));
int index1 = tensor_to_list_idx(sigma1indices3.at(i));
int index2 = tensor_to_list_idx(sigma2indices3.at(i));
int index3 = tensor_to_list_idx(sigma3indices3.at(i));
sigma_index3[index0]=0;
sigma_index3[index1]=1;
sigma_index3[index2]=2;
sigma_index3[index3]=3;
short int sign[4]={1,-1,-1,-1};
// plus->true->1
helfactor3[1][index0] = sign[sigma0indices3.at(i)[1]];
helfactor3[1][index1] = sign[sigma1indices3.at(i)[1]];
helfactor3[1][index2] = sign[sigma2indices3.at(i)[1]];
helfactor3[1][index3] = sign[sigma3indices3.at(i)[1]];
// minus->false->0
helfactor3[0][index0] = sign[sigma0indices3.at(i)[0]]
* sign[sigma0indices3.at(i)[2]];
helfactor3[0][index1] = sign[sigma1indices3.at(i)[0]]
* sign[sigma1indices3.at(i)[2]];
helfactor3[0][index2] = sign[sigma2indices3.at(i)[0]]
* sign[sigma2indices3.at(i)[2]];
helfactor3[0][index3] = sign[sigma3indices3.at(i)[0]]
* sign[sigma3indices3.at(i)[2]];
}
return true;
} // end init_sigma_index
}
void init_sigma_index(){
static const bool initialised = init_sigma_index_helper();
(void) initialised; // shut up compiler warnings
}
}
diff --git a/src/Wjets.cc b/src/Wjets.cc
index ba1da6c..1e48e5e 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,1666 +1,1666 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Wjets.hh"
#include <array>
#include <iostream>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/jets.hh"
#include "HEJ/Tensor.hh"
using HEJ::Tensor;
using HEJ::init_sigma_index;
using HEJ::metric;
using HEJ::rank3_current;
using HEJ::rank5_current;
using HEJ::eps;
using HEJ::to_tensor;
using HEJ::Helicity;
using HEJ::angle;
using HEJ::square;
using HEJ::flip;
using HEJ::ParticleProperties;
namespace helicity = HEJ::helicity;
namespace { // Helper Functions
// FKL W Helper Functions
double WProp (const HLV & plbar, const HLV & pl, ParticleProperties const & wprop){
COM propW = COM(0.,-1.)/( (pl+plbar).m2() - wprop.mass*wprop.mass
+ COM(0.,1.)*wprop.mass*wprop.width);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
namespace {
// FKL current including W emission off negative helicities
// See eq. (87) {eq:jW-} in developer manual
// Note that the terms are rearranged
Tensor<1> jW_minus(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl
){
using HEJ::helicity::minus;
const double tWin = (pa-pl-plbar).m2();
const double tWout = (p1+pl+plbar).m2();
// C++ arithmetic operators are evaluated left-to-right,
// so the following first computes complex scalar coefficients,
// which then multiply a current, reducing the number
// of multiplications
return 2.*(
+ angle(p1, pl)*square(p1, plbar)/tWout
+ square(pa, plbar)*angle(pa, pl)/tWin
)*HEJ::current(p1, pa, helicity::minus)
+ 2.*angle(p1, pl)*square(pl, plbar)/tWout
*HEJ::current(pl, pa, helicity::minus)
+ 2.*square(pa, plbar)*angle(pl, plbar)/tWin
*HEJ::current(p1, plbar, helicity::minus);
}
}
// FKL current including W emission
// see eqs. (87), (88) {eq:jW-}, {eq:jW+} in developer manual
Tensor<1> jW(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl,
Helicity h
){
if(h == helicity::minus) {
return jW_minus(pa, p1, plbar, pl);
}
return jW_minus(pa, p1, pl, plbar).complex_conj();
}
/**
* @brief Contraction of the \tilde{U}_1 tensor
*
* This is the contraction of the tensor defined in eq:U1tensor
* in the developer manual with the uno gluon polarisation vector,
* the leptonic and the partonic current (see eq:wunocurrent) in the
* developer manual)
*/
COM U1contr(
HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity h1,
HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
) {
const HLV pW = pl+plbar;
const HLV q1g = pa-pW-p1-pg;
const HLV q1 = pa-p1-pW;
const HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+pg).m2();
if(h1 == helicity::plus) {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity +++
return (4.*sqrt(2.)*(-1.*taW*angle(pa,pb)*square(p1,plbar)*
(s1W*angle(p1,pg)*angle(pg,pl)*square(p1,p2)*square(p2,pg) -
1.*(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar))*
((s1g + s1W)*angle(p1,pg)*square(p1,p2) +
s1g*(angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar)))) +
s1gW*s1W*angle(p1,pg)*angle(pa,pl)*pow(square(p1,p2),2.)*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*square(p2,pg));
}
else {
// helicity ++-
return (sqrt(2.)*(taW*angle(pa,pb)*(4.*s1g*square(p1,plbar)*
(angle(p1,pl)*square(p1,pg)*(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar)) +
angle(pl,plbar)*(angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg) +
angle(pb,pl)*square(p2,pl) + angle(pb,plbar)*square(p2,plbar))*square(pg,plbar)) +
square(p1,pg)*(4.*angle(p1,pb)*square(p1,plbar)*
((s1g + s1W)*angle(p1,pl)*square(p1,p2) - 1.*s1W*angle(pg,pl)*square(p2,pg) +
s1W*angle(pl,plbar)*square(p2,plbar)) -
4.*s1W*angle(pb,pg)*(angle(p1,pl)*square(p1,p2) - 1.*angle(pg,pl)*square(p2,pg) +
angle(pl,plbar)*square(p2,plbar))*square(pg,plbar))) +
4.*s1gW*s1W*angle(pa,pl)*square(p1,pg)*
(angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg))*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*angle(pb,pg));
}
}
else {
if(hg == helicity::plus) {
// helicity +-+
return (sqrt(2.)*(4.*taW*angle(p2,pa)*square(p1,plbar)*
(s1W*angle(p1,pg)*angle(pg,pl)*square(p1,pb)*square(pb,pg) -
1.*(angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar))*
((s1g + s1W)*angle(p1,pg)*square(p1,pb) +
s1g*(angle(pg,pl)*square(pb,pl) + angle(pg,plbar)*square(pb,plbar)))) -
4.*s1gW*s1W*angle(p1,pg)*angle(pa,pl)*pow(square(p1,pb),2.)*
(angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*square(pb,pg));
}
else {
// helicity +--
return (-1.*sqrt(2.)*(taW*angle(p2,pa)*(4.*s1g*square(p1,plbar)*
(angle(p1,pl)*square(p1,pg)*(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
angle(p2,plbar)*square(pb,plbar)) +
angle(pl,plbar)*(angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg) +
angle(p2,pl)*square(pb,pl) + angle(p2,plbar)*square(pb,plbar))*square(pg,plbar)) +
square(p1,pg)*(4.*angle(p1,p2)*square(p1,plbar)*
((s1g + s1W)*angle(p1,pl)*square(p1,pb) - 1.*s1W*angle(pg,pl)*square(pb,pg) +
s1W*angle(pl,plbar)*square(pb,plbar)) -
4.*s1W*angle(p2,pg)*(angle(p1,pl)*square(p1,pb) - 1.*angle(pg,pl)*square(pb,pg) +
angle(pl,plbar)*square(pb,plbar))*square(pg,plbar))) +
4.*s1gW*s1W*angle(pa,pl)*square(p1,pg)*
(angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg))*
(angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*angle(p2,pg));
}
}
}
else {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity -++
return (
sqrt(2.)*(-4.*s1gW*s1W*angle(p1,pg)*(angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg))*
(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*square(pa,plbar) +
taW*angle(pg,pl)*square(p2,pa)*(4.*s1g*angle(p1,pl)*
(angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar))*square(pl,plbar) +
4.*s1W*angle(p1,pg)*square(p2,pg)*
(angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pg)*square(pg,plbar) -
1.*angle(pb,pl)*square(pl,plbar))) -
4.*taW*angle(p1,pg)*angle(p1,pl)*square(p2,pa)*
(square(p1,plbar)*((s1g + s1W)*angle(p1,pb)*square(p1,p2) +
s1g*(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar))) -
1.*s1W*square(p1,p2)*(angle(pb,pg)*square(pg,plbar) + angle(pb,pl)*square(pl,plbar))))
)/(s1g*s1gW*s1W*taW*square(p2,pg));
}
else {
// helicity -+-
return (-1.*sqrt(2.)*(4.*s1W*angle(p1,pb)*square(p1,pg)*
(s1gW*angle(p1,pb)*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
square(pa,plbar) - 1.*taW*angle(p1,pl)*angle(pb,pg)*square(p2,pa)*square(pg,plbar)) +
taW*angle(p1,pl)*square(p2,pa)*((s1g + s1W)*angle(p1,pb)*square(p1,pg) +
s1g*(angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar)))*
(4.*angle(p1,pb)*square(p1,plbar) - 4.*angle(pb,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*angle(pb,pg));
}
}
else {
if(hg == helicity::plus) {
// helicity --+
return (sqrt(2.)*(4.*s1gW*s1W*angle(p1,pg)*square(pa,plbar)*
(angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg))*
(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)) +
taW*square(pa,pb)*(-1.*angle(pg,pl)*
(4.*s1g*angle(p1,pl)*(angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg) +
angle(p2,pl)*square(pb,pl) + angle(p2,plbar)*square(pb,plbar))*square(pl,plbar) +
4.*s1W*angle(p1,pg)*square(pb,pg)*
(angle(p1,p2)*square(p1,plbar) - 1.*angle(p2,pg)*square(pg,plbar) -
1.*angle(p2,pl)*square(pl,plbar))) +
4.*angle(p1,pg)*angle(p1,pl)*(square(p1,plbar)*
((s1g + s1W)*angle(p1,p2)*square(p1,pb) +
s1g*(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
angle(p2,plbar)*square(pb,plbar))) -
1.*s1W*square(p1,pb)*(angle(p2,pg)*square(pg,plbar) + angle(p2,pl)*square(pl,plbar)))))
)/(s1g*s1gW*s1W*taW*square(pb,pg));
}
else {
// helicity ---
return (sqrt(2.)*(s1W*angle(p1,p2)*square(p1,pg)*
(4.*s1gW*angle(p1,p2)*square(pa,plbar)*
(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)) -
4.*taW*angle(p1,pl)*angle(p2,pg)*square(pa,pb)*square(pg,plbar)) +
taW*angle(p1,pl)*square(pa,pb)*((s1g + s1W)*angle(p1,p2)*square(p1,pg) +
s1g*(angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))*
(4.*angle(p1,p2)*square(p1,plbar) - 4.*angle(p2,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*angle(p2,pg));
}
}
}
}
/**
* @brief Contraction of the \tilde{U}_2 tensor
*
* This is the contraction of the tensor defined in eq:U2tensor
* in the developer manual with the uno gluon polarisation vector,
* the leptonic and the partonic current (see eq:wunocurrent) in the
* developer manual)
*/
COM U2contr(
HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity h1,
HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
) {
const HLV pW = pl+plbar;
const HLV q1g=pa-pW-p1-pg;
const HLV q1 = pa-p1-pW;
const HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double s1W = (p1+pW).m2();
const double tag = (pa-pg).m2();
const double taWg = (pa-pW-pg).m2();
if(h1 == helicity::plus) {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity +++
return (-4.*sqrt(2.)*(taW*angle(pa,pg)*(taWg*square(p1,plbar)*
(angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg))*
(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) -
1.*s1W*angle(pg,pl)*square(p1,p2)*square(p2,pg)*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pg)*square(pg,plbar) +
angle(pb,pl)*square(pl,plbar))) +
s1W*angle(pa,pl)*square(p1,p2)*(tag*
(angle(pa,pg)*(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar))*square(pa,plbar) +
angle(pg,pl)*(angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg) +
angle(pb,pl)*square(p2,pl) + angle(pb,plbar)*square(p2,plbar))*square(pl,plbar)) +
angle(pa,pg)*square(p2,pa)*((tag + taW)*angle(pa,pb)*square(pa,plbar) +
taW*(angle(pb,pg)*square(pg,plbar) + angle(pb,pl)*square(pl,plbar))))))/
(s1W*tag*taW*taWg*square(p2,pg));
}
else {
// helicity ++-
return (-4.*sqrt(2.)*(s1W*tag*angle(pa,pl)*square(p1,p2)*
(angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar))*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)) -
1.*angle(pa,pb)*square(pa,pg)*(taW*taWg*angle(pa,pb)*square(p1,plbar)*
(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) +
s1W*angle(pa,pl)*square(p1,p2)*
(taW*angle(pb,pg)*square(pg,plbar) +
(tag + taW)*(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))))))/
(s1W*tag*taW*taWg*angle(pb,pg));
}
}
else {
if(hg == helicity::plus) {
// helicity +-+
return (-4.*sqrt(2.)*(taW*angle(pa,pg)*(taWg*square(p1,plbar)*
(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg))*
(angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar)) -
1.*s1W*square(p1,pb)*(angle(pa,pl)*square(pa,pb) + angle(pg,pl)*square(pb,pg))*
(angle(p2,pg)*square(pg,plbar) + angle(p2,pl)*square(pl,plbar))) +
s1W*square(p1,pb)*(tag*angle(pa,pl)*
(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
angle(p2,plbar)*square(pb,plbar))*
(angle(pa,pg)*square(pa,plbar) + angle(pg,pl)*square(pl,plbar)) +
angle(p2,pa)*(angle(pa,pg)*square(pa,plbar)*
((tag + taW)*angle(pa,pl)*square(pa,pb) + taW*angle(pg,pl)*square(pb,pg)) +
tag*angle(pa,pl)*angle(pg,pl)*square(pa,pb)*square(pl,plbar)))))/
(s1W*tag*taW*taWg*square(pb,pg));
}
else {
// helicity +--
return (sqrt(2.)*(4.*taW*angle(p2,pa)*square(pa,pg)*
(taWg*angle(p2,pa)*square(p1,plbar)*
(angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar)) -
1.*s1W*angle(p2,pg)*angle(pa,pl)*square(p1,pb)*square(pg,plbar)) +
s1W*angle(pa,pl)*square(p1,pb)*((tag + taW)*angle(p2,pa)*square(pa,pg) +
tag*(angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))*
(4.*angle(p2,pa)*square(pa,plbar) - 4.*angle(p2,pl)*square(pl,plbar))))/
(s1W*tag*taW*taWg*angle(p2,pg));
}
}
}
else {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity -++
return (sqrt(2.)*(-4.*s1W*angle(p1,pb)*(taW*angle(pa,pg)*angle(pg,pl)*square(p2,pa)*square(p2,pg) -
1.*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
((tag + taW)*angle(pa,pg)*square(p2,pa) +
tag*(angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))))*square(pa,plbar)
+ 4.*taW*taWg*angle(p1,pl)*angle(pa,pg)*pow(square(p2,pa),2.)*
(angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar))))/
(s1W*tag*taW*taWg*square(p2,pg));
}
else {
// helicity -+-
return (sqrt(2.)*(4.*s1W*angle(p1,pb)*(tag*angle(pl,plbar)*
(angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar))*square(pa,plbar)*square(pg,plbar) +
taW*angle(pg,pl)*square(p2,pg)*square(pa,pg)*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pg)*square(pg,plbar)) -
1.*square(pa,pg)*((tag*angle(pa,pl)*
(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar)) +
angle(pa,pb)*((tag + taW)*angle(pa,pl)*square(p2,pa) +
taW*angle(pl,plbar)*square(p2,plbar)))*square(pa,plbar) +
taW*angle(pb,pg)*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
square(pg,plbar))) - 4.*taW*taWg*angle(p1,pl)*
(angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg))*square(pa,pg)*
(angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar))))/
(s1W*tag*taW*taWg*angle(pb,pg));
}
}
else {
if(hg == helicity::plus) {
// helicity --+
return (4.*sqrt(2.)*(angle(p1,p2)*(taW*taWg*angle(p1,pl)*angle(pa,pg)*pow(square(pa,pb),2.)*
square(p1,plbar) + s1W*square(pa,plbar)*
(tag*(angle(pg,pl)*square(pb,pl) + angle(pg,plbar)*square(pb,plbar))*
(-1.*angle(pa,pl)*square(pa,pb) + angle(pl,plbar)*square(pb,plbar)) +
angle(pa,pg)*square(pa,pb)*((tag + taW)*angle(pa,pl)*square(pa,pb) +
taW*angle(pg,pl)*square(pb,pg) - 1.*(tag + taW)*angle(pl,plbar)*square(pb,plbar))))
- 1.*taW*taWg*angle(p1,pl)*angle(p2,pl)*angle(pa,pg)*pow(square(pa,pb),2.)*square(pl,plbar))
)/(s1W*tag*taW*taWg*square(pb,pg));
}
else {
// helicity ---
return (sqrt(2.)*(s1W*angle(p1,p2)*(-4.*square(pa,pg)*square(pa,plbar)*
(tag*angle(pa,pl)*(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
angle(p2,plbar)*square(pb,plbar)) +
angle(p2,pa)*((tag + taW)*angle(pa,pl)*square(pa,pb) +
taW*angle(pg,pl)*square(pb,pg) - 1.*taW*angle(pl,plbar)*square(pb,plbar))) +
4.*tag*angle(pl,plbar)*square(pa,plbar)*
(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
angle(p2,plbar)*square(pb,plbar))*square(pg,plbar) +
4.*taW*angle(p2,pg)*square(pa,pg)*
(angle(pa,pl)*square(pa,pb) + angle(pg,pl)*square(pb,pg) -
1.*angle(pl,plbar)*square(pb,plbar))*square(pg,plbar)) -
4.*taW*taWg*angle(p1,pl)*square(pa,pg)*
(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg))*
(angle(p1,p2)*square(p1,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
(s1W*tag*taW*taWg*angle(p2,pg));
}
}
}
}
/**
* @brief Contraction of the \tilde{L} tensor
*
* This is the contraction of the tensor defined in eq:Ltensor
* in the developer manual with the uno gluon polarisation vector,
* the leptonic and the partonic current (see eq:wunocurrent) in the
* developer manual)
*/
COM Lcontr(
HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity h1,
HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
) {
const HLV pW = pl+plbar;
const HLV q1g=pa-pW-p1-pg;
const HLV q1 = pa-p1-pW;
const HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double taW1 = (pa-pW-p1).m2();
const double s1W = (p1+pW).m2();
const double s2g = (p2+pg).m2();
const double sbg = (pb+pg).m2();
if(h1 == helicity::plus) {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity +++
return (-2.*sqrt(2.)*(angle(pb,pg)*q1g.m2()*square(p2,pb)*
(taW*angle(pa,pb)*square(p1,plbar)*
(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) +
s1W*angle(pa,pl)*square(p1,p2)*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))) +
2.*sbg*(taW*square(p1,plbar)*(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar))*
(angle(pa,pg)*angle(pb,pg)*square(p2,pg) +
angle(pa,pb)*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))) +
s1W*angle(pa,pl)*square(p1,p2)*
((angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) + angle(pg,pl)*square(p2,pl) +
angle(pg,plbar)*square(p2,plbar))*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)) +
angle(pb,pg)*square(p2,pg)*(angle(pa,pg)*square(pa,plbar) +
angle(pg,pl)*square(pl,plbar))))))/(s1W*sbg*taW*taW1*square(p2,pg));
}
else {
// helicity ++-
return (-1.*sqrt(2.)*(s2g*taW*angle(pa,pb)*square(p1,plbar)*
(4.*(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar))*
(angle(p1,pb)*square(p1,pg) - 1.*angle(pa,pb)*square(pa,pg) +
angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar)) +
4.*angle(pb,pg)*square(p2,pg)*(angle(p1,pl)*square(p1,pg) +
angle(pl,plbar)*square(pg,plbar))) +
s1W*s2g*angle(pa,pl)*(4.*angle(pb,pg)*square(p1,pg)*square(p2,pg) -
4.*angle(pa,pb)*square(p1,p2)*square(pa,pg) +
4.*square(p1,p2)*(angle(p1,pb)*square(p1,pg) + angle(pb,pl)*square(pg,pl) +
angle(pb,plbar)*square(pg,plbar)))*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)) -
2.*angle(p2,pb)*q1g.m2()*square(p2,pg)*
(taW*angle(pa,pb)*square(p1,plbar)*
(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) +
s1W*angle(pa,pl)*square(p1,p2)*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)))))/
(s1W*s2g*taW*taW1*angle(pb,pg));
}
}
else {
if(hg == helicity::plus) {
// helicity +-+
return (-1.*sqrt(2.)*(2.*taW*square(p1,plbar)*(angle(p1,pl)*square(p1,pb) +
angle(pl,plbar)*square(pb,plbar))*
(2.*s2g*angle(pa,pg)*(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg)) +
angle(p2,pa)*(angle(p2,pg)*q1g.m2()*square(p2,pb) -
2.*s2g*(angle(p1,pg)*square(p1,pb) + angle(pg,pl)*square(pb,pl) +
angle(pg,plbar)*square(pb,plbar)))) +
s1W*angle(pa,pl)*square(p1,pb)*(4.*s2g*square(pa,plbar)*
(angle(pa,pg)*(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg)) -
1.*angle(p2,pa)*(angle(p1,pg)*square(p1,pb) + angle(pg,pl)*square(pb,pl) +
angle(pg,plbar)*square(pb,plbar))) +
s2g*(-4.*angle(p2,pl)*angle(pa,pg)*square(pa,pb) +
4.*angle(p2,pg)*angle(pg,pl)*square(pb,pg) +
4.*angle(p2,pl)*(angle(p1,pg)*square(p1,pb) + angle(pg,pl)*square(pb,pl) +
angle(pg,plbar)*square(pb,plbar)))*square(pl,plbar) +
2.*angle(p2,pg)*q1g.m2()*square(p2,pb)*
(angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar)))))/
(s1W*s2g*taW*taW1*square(pb,pg));
}
else {
// helicity +--
return (sqrt(2.)*(2.*taW*angle(p2,pa)*square(p1,plbar)*
(2.*sbg*angle(p2,pg)*square(pb,pg)*
(angle(p1,pl)*square(p1,pg) + angle(pl,plbar)*square(pg,plbar)) +
(angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar))*
(angle(p2,pb)*q1g.m2()*square(pb,pg) +
2.*sbg*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))) +
2.*s1W*angle(pa,pl)*(2.*sbg*angle(p1,p2)*square(p1,pb)*square(p1,pg) +
2.*sbg*angle(p2,pa)*square(p1,pb)*square(pa,pg) +
angle(p2,pb)*q1g.m2()*square(p1,pb)*square(pb,pg) +
2.*sbg*angle(p2,pg)*square(p1,pg)*square(pb,pg) +
2.*sbg*angle(p2,pl)*square(p1,pb)*square(pg,pl) +
2.*sbg*angle(p2,plbar)*square(p1,pb)*square(pg,plbar))*
(angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
(s1W*sbg*taW*taW1*angle(p2,pg));
}
}
}
else {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity -++
return (sqrt(2.)*(2.*s1W*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
(2.*sbg*angle(p1,pg)*angle(pb,pg)*square(p2,pg) +
angle(p1,pb)*(angle(pb,pg)*q1g.m2()*square(p2,pb) +
2.*sbg*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))))*square(pa,plbar) +
taW*angle(p1,pl)*square(p2,pa)*(4.*sbg*square(p1,plbar)*
(angle(p1,pg)*angle(pb,pg)*square(p2,pg) +
angle(p1,pb)*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))) -
4.*sbg*(angle(pb,pg)*angle(pg,pl)*square(p2,pg) +
angle(pb,pl)*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar)))*square(pl,plbar) +
2.*angle(pb,pg)*q1g.m2()*square(p2,pb)*
(angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar)))))/
(s1W*sbg*taW*taW1*square(p2,pg));
}
else {
// helicity -+-
return (sqrt(2.)*((angle(p1,pb)*square(pa,plbar)*
(((angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
(4.*angle(p1,pb)*square(p1,pg) - (2.*angle(p2,pb)*q1g.m2()*square(p2,pg))/s2g -
4.*angle(pa,pb)*square(pa,pg) + 4.*angle(pb,pl)*square(pg,pl) +
4.*angle(pb,plbar)*square(pg,plbar)))/angle(pb,pg) +
square(p2,pg)*(-4.*angle(pa,pl)*square(pa,pg) + 4.*angle(pl,plbar)*square(pg,plbar))))/
taW + (2.*angle(p1,pl)*(2.*s2g*angle(p1,pb)*square(p1,pg)*square(p2,pa) -
1.*angle(p2,pb)*q1g.m2()*square(p2,pa)*square(p2,pg) +
2.*s2g*(-1.*angle(pa,pb)*square(p2,pa)*square(pa,pg) -
1.*angle(pb,pg)*square(p2,pg)*square(pa,pg) +
square(p2,pa)*(angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar))))*
(angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar)))/(s1W*s2g*angle(pb,pg))
))/taW1;
}
}
else {
if(hg == helicity::plus) {
// helicity --+
return (2.*sqrt(2.)*(-1.*angle(p1,p2)*(2.*s2g*angle(p1,pg)*square(p1,pb) -
1.*angle(p2,pg)*q1g.m2()*square(p2,pb) +
2.*s2g*(-1.*angle(pa,pg)*square(pa,pb) + angle(pg,pl)*square(pb,pl) +
angle(pg,plbar)*square(pb,plbar)))*
(taW*angle(p1,pl)*square(p1,plbar)*square(pa,pb) +
s1W*square(pa,plbar)*(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)))
+ taW*angle(p1,pl)*square(pa,pb)*
(angle(p2,pg)*(-1.*angle(p2,pl)*q1g.m2()*square(p2,pb) +
2.*s2g*angle(pg,pl)*square(pb,pg)) +
2.*s2g*angle(p2,pl)*(-1.*angle(pa,pg)*square(pa,pb) + angle(pg,pl)*square(pb,pl) +
angle(pg,plbar)*square(pb,plbar)))*square(pl,plbar) -
2.*s2g*angle(p1,pg)*(s1W*angle(p2,pg)*square(pa,plbar)*square(pb,pg)*
(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)) +
taW*angle(p1,pl)*square(pa,pb)*
(angle(p2,pg)*square(p1,plbar)*square(pb,pg) -
1.*angle(p2,pl)*square(p1,pb)*square(pl,plbar)))))/(s1W*s2g*taW*taW1*square(pb,pg));
}
else {
// helicity ---
return (2.*sqrt(2.)*(-1.*angle(p1,p2)*(angle(p2,pb)*q1g.m2()*square(pb,pg)*
(taW*angle(p1,pl)*square(p1,plbar)*square(pa,pb) +
s1W*square(pa,plbar)*(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar))
) + 2.*sbg*(taW*angle(p1,pl)*square(p1,plbar) + s1W*angle(pa,pl)*square(pa,plbar))*
(angle(p2,pg)*square(pa,pg)*square(pb,pg) +
square(pa,pb)*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar))) -
2.*s1W*sbg*angle(pl,plbar)*square(pa,plbar)*
(angle(p2,pg)*square(pb,pg)*square(pg,plbar) +
square(pb,plbar)*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))) +
taW*angle(p1,pl)*angle(p2,pl)*(2.*sbg*angle(p2,pg)*square(pa,pg)*square(pb,pg) +
square(pa,pb)*(angle(p2,pb)*q1g.m2()*square(pb,pg) +
2.*sbg*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar))))*square(pl,plbar)))/
(s1W*sbg*taW*taW1*angle(p2,pg));
}
}
}
}
/**
* @brief W+Jets Unordered Contribution Helper Functions
* @returns result of equation (4.1.28) in Helen's Thesis (p.100)
*/
double jM2Wuno(HLV pg, HLV p1,HLV plbar, HLV pl, HLV pa, Helicity h1,
HLV p2, HLV pb, Helicity h2, Helicity pol,
ParticleProperties const & wprop
){
//@TODO Simplify the below (less Tensor class?)
init_sigma_index();
HLV pW = pl+plbar;
HLV q1g=pa-pW-p1-pg;
HLV q1 = pa-p1-pW;
HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double taW1 = (pa-pW-p1).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+pg).m2();
const double s2g = (p2+pg).m2();
const double sbg = (pb+pg).m2();
const double tag = (pa-pg).m2();
const double taWg = (pa-pW-pg).m2();
//use p1 as ref vec in pol tensor
Tensor<1> epsg = eps(pg,p2,pol);
Tensor<1> epsW = HEJ::current(pl,plbar,helicity::minus);
Tensor<1> j2b = HEJ::current(p2,pb,h2);
Tensor<1> Tq1q2 = to_tensor(
(pb/sbg + p2/s2g)*(q1 - pg).m2() + 2*q1 - pg
);
Tensor<3> J31a = rank3_current(p1, pa, h1);
Tensor<2> J2_qaW =J31a.contract((pa-pW)/taW/taW1, 2);
Tensor<2> J2_p1W =J31a.contract((p1+pW)/s1W/taW1, 2);
Tensor<3> L1a = outer(Tq1q2, J2_qaW);
Tensor<3> L1b = outer(Tq1q2, J2_p1W);
Tensor<3> L2a = outer(-2*pg, J2_qaW);
Tensor<3> L2b = outer(-2*pg, J2_p1W);
Tensor<3> L3 = outer(metric(), J2_qaW.contract(2*pg-q1,1)+J2_p1W.contract(2*pg-q1,2));
Tensor<3> L(0.);
Tensor<5> J51a = rank5_current(p1, pa, h1);
Tensor<4> J_qaW = J51a.contract((pa-pW),4);
Tensor<4> J_qag = J51a.contract(pa-pg,4);
Tensor<4> J_p1gW = J51a.contract(p1+pg+pW,4);
Tensor<3> U1a = J_qaW.contract(p1+pg,2);
Tensor<3> U1b = J_p1gW.contract(p1+pg,2);
Tensor<3> U1c = J_p1gW.contract(p1+pW,2);
Tensor<3> U1(0.);
Tensor<3> U2a = J_qaW.contract(pa-pg-pW,2);
Tensor<3> U2b = J_qag.contract(pa-pg-pW,2);
Tensor<3> U2c = J_qag.contract(p1+pW,2);
Tensor<3> U2(0.);
- for(int nu=0; nu<4;nu++){
- for(int mu=0;mu<4;mu++){
- for(int rho=0;rho<4;rho++){
+ for(int nu=0; nu<4;++nu){
+ for(int mu=0;mu<4;++mu){
+ for(int rho=0;rho<4;++rho){
L(nu, mu, rho) = L1a(nu,mu,rho) + L1b(nu,rho,mu)
+ L2a(mu,nu,rho) + L2b(mu,rho,nu) + L3(mu,nu,rho);
U1(nu, mu, rho) = U1a(nu, mu, rho) / (s1g*taW)
+ U1b(nu,rho,mu) / (s1g*s1gW) + U1c(rho,nu,mu) / (s1W*s1gW);
U2(nu,mu,rho) = U2a(mu,nu,rho) / (taWg*taW)
+ U2b(mu,rho,nu) / (taWg*tag) + U2c(rho,mu,nu) / (s1W*tag);
}
}
}
COM X = ((((U1-L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
COM Y = ((((U2+L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
double amp = HEJ::C_A*HEJ::C_F*HEJ::C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
double t1 = q1g.m2();
double t2 = q2.m2();
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
//Divide by t-channels
amp/=(t1*t2);
return amp;
}
/**
* @brief New implementation of unordered W+Jets current
*
* See eq:wunocurrent in the developer manual for the definition
* of this current
*
* The aim is to replace jM2Wuno and eventually all uses of the Tensor class.
* Explicit tensor contractions are rather slow and the initialisation code
* in Tensor.cc is not very transparent.
*
* At the moment, this function _only_ works for positive-energy spinors,
* so jM2Wuno has to be kept for qqbar emission.
*/
double jM2Wuno_pos_energy(
HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity const h1,
HLV const & p2, HLV const & pb, Helicity const h2, Helicity const pol,
ParticleProperties const & wprop
){
using HEJ::C_A;
using HEJ::C_F;
const COM U1 = U1contr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
const COM U2 = U2contr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
const COM L = Lcontr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
const COM X = U1 - L;
const COM Y = U2 + L;
double amp = C_A*C_F*C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
amp *= WProp(plbar, pl, wprop);
const HLV q1g = pa-pl-plbar-p1-pg;
const HLV q2 = p2-pb;
const double t1 = q1g.m2();
const double t2 = q2.m2();
amp /= (t1*t2);
return amp;
}
// Relevant Wqqx Helper Functions.
//g->qxqlxl (Calculates gluon to qqx Current. See JV_\mu in WSubleading Notes)
Tensor <1> gtqqxW(HLV pq,HLV pqbar,HLV pl,HLV plbar){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
//blank 3 Gamma Current
Tensor<3> JV23 = rank3_current(pq,pqbar,helicity::minus);
// Components of g->qqW before W Contraction
Tensor<2> JV1 = JV23.contract((pq + pl + plbar),2)/(s2AB);
Tensor<2> JV2 = JV23.contract((pqbar + pl + plbar),2)/(s3AB);
// g->qqW Current. Note Minus between terms due to momentum flow.
// Also note: (-I)^2 from W vert. (I) from Quark prop.
Tensor<1> JVCur = (JV1.contract(ABCur,1) - JV2.contract(ABCur,2))*COM(0.,-1.);
return JVCur;
}
// Helper Functions Calculate the Crossed Contribution
Tensor <2> MCrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar, HLV pl,
HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCross?
// Useful propagator factors
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV q1, q3;
q1=pa;
- for(int i=0; i<nabove+1;i++){
+ for(int i=0; i<nabove+1;++i){
q1=q1-partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double tcro1=(q3+pq).m2();
double tcro2=(q1-pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar,helicity::minus);
//Blank 5 gamma Current
Tensor<5> J523 = rank5_current(pq,pqbar,helicity::minus);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_q3q = J523.contract((q3 + pq),2);
Tensor<4> J_2AB = J523.contract((pq + pl + plbar),2);
// Components of Crossed Vertex Contribution
Tensor<3> Xcro1 = J_q3q.contract((pqbar + pl + plbar),3);
Tensor<3> Xcro2 = J_q3q.contract((q1 - pqbar),3);
Tensor<3> Xcro3 = J_2AB.contract((q1 - pqbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xcro1Cont = Xcro1.contract(ABCur,3)/(tcro1*s3AB);
Tensor<2> Xcro2Cont = Xcro2.contract(ABCur,2)/(tcro1*tcro2);
Tensor<2> Xcro3Cont = Xcro3.contract(ABCur,1)/(s2AB*tcro2);
//Initialise the Crossed Vertex Object
Tensor<2> Xcro(0.);
- for(int mu=0; mu<4;mu++){
- for(int nu=0;nu<4;nu++){
+ for(int mu=0; mu<4;++mu){
+ for(int nu=0;nu<4;++nu){
Xcro(mu,nu) = -(-Xcro1Cont(nu,mu)+Xcro2Cont(nu,mu)+Xcro3Cont(nu,mu));
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar,
HLV pl, HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncross?
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV q1, q3;
q1=pa;
- for(int i=0; i<nabove+1;i++){
+ for(int i=0; i<nabove+1;++i){
q1=q1-partons.at(i);
}
q3 = q1 - pl - plbar - pq - pqbar;
double tunc1 = (q1-pq).m2();
double tunc2 = (q3+pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
//Blank 5 gamma Current
Tensor<5> J523 = rank5_current(pq,pqbar,helicity::minus);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_2AB = J523.contract((pq + pl + plbar),2);
Tensor<4> J_q1q = J523.contract((q1 - pq),2);
// 2 Contractions taken care of.
Tensor<3> Xunc1 = J_2AB.contract((q3 + pqbar),3);
Tensor<3> Xunc2 = J_q1q.contract((q3 + pqbar),3);
Tensor<3> Xunc3 = J_q1q.contract((pqbar + pl + plbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xunc1Cont = Xunc1.contract(ABCur,1)/(s2AB*tunc2);
Tensor<2> Xunc2Cont = Xunc2.contract(ABCur,2)/(tunc1*tunc2);
Tensor<2> Xunc3Cont = Xunc3.contract(ABCur,3)/(tunc1*s3AB);
//Initialise the Uncrossed Vertex Object
Tensor<2> Xunc(0.);
- for(int mu=0; mu<4;mu++){
- for(int nu=0;nu<4;nu++){
+ for(int mu=0; mu<4;++mu){
+ for(int nu=0;nu<4;++nu){
Xunc(mu,nu) = -(- Xunc1Cont(mu,nu)+Xunc2Cont(mu,nu) +Xunc3Cont(mu,nu));
}
}
return Xunc;
}
// Helper Functions Calculate the g->qqxW (Eikonal) Contributions
Tensor <2> MSymW(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
HLV pl,HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MSym?
double sa2=(pa+pq).m2();
double s12=(p1+pq).m2();
double sa3=(pa+pqbar).m2();
double s13=(p1+pqbar).m2();
double saA=(pa+pl).m2();
double s1A=(p1+pl).m2();
double saB=(pa+plbar).m2();
double s1B=(p1+plbar).m2();
double sb2=(pb+pq).m2();
double s42=(p4+pq).m2();
double sb3=(pb+pqbar).m2();
double s43=(p4+pqbar).m2();
double sbA=(pb+pl).m2();
double s4A=(p4+pl).m2();
double sbB=(pb+plbar).m2();
double s4B=(p4+plbar).m2();
double s23AB=(pl+plbar+pq+pqbar).m2();
HLV q1,q3;
q1=pa;
- for(int i=0;i<nabove+1;i++){
+ for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
q3=q1-pq-pqbar-plbar-pl;
double t1 = (q1).m2();
double t3 = (q3).m2();
// g->qqW Current (Factors of sqrt2 dealt with in this function.)
Tensor<1> JV = gtqqxW(pq,pqbar,pl,plbar);
// 1a gluon emisson Contribution
Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13+s1A+s1B))
+ pa*(t1/(sa2+sa3+saA+saB)) );
Tensor<2> X1aCont = X1a.contract(JV,3);
//4b gluon emission Contribution
Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43+s4A+s4B))
+ pb*(t3/(sb2+sb3+sbA+sbB)) );
Tensor<2> X4bCont = X4b.contract(JV,3);
//Set up each term of 3G diagram.
Tensor<3> X3g1 = outer(q1+pq+pqbar+pl+plbar, metric());
Tensor<3> X3g2 = outer(q3-pq-pqbar-pl-plbar, metric());
Tensor<3> X3g3 = outer(q1+q3, metric());
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(JV,3);
Tensor<2> X3g2Cont = X3g2.contract(JV,2);
Tensor<2> X3g3Cont = X3g3.contract(JV,1);
// XSym is an amalgamation of x1a, X4b and X3g.
// Makes sense from a colour factor point of view.
Tensor<2>Xsym(0.);
- for(int mu=0; mu<4;mu++){
- for(int nu=0;nu<4;nu++){
+ for(int mu=0; mu<4;++mu){
+ for(int nu=0;nu<4;++nu){
Xsym(mu,nu) = (X3g1Cont(nu,mu) + X3g2Cont(mu,nu) - X3g3Cont(nu,mu))
+ (X1aCont(mu,nu) - X4bCont(mu,nu));
}
}
return Xsym/s23AB;
}
Tensor <2> MCross(HLV pa, HLV pq, HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCrossW?
HLV q1;
q1=pa;
- for(int i=0;i<nabove+1;i++){
+ for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
double t2=(q1-pqbar).m2();
//Blank 3 gamma Current
Tensor<3> J323 = rank3_current(pq,pqbar,hq);
// 2 gamma current (with 1 contraction already).
Tensor<2> XCroCont = J323.contract((q1-pqbar),2)/(t2);
//Initialise the Crossed Vertex
Tensor<2> Xcro(0.);
- for(int mu=0; mu<4;mu++){
- for(int nu=0;nu<4;nu++){
+ for(int mu=0; mu<4;++mu){
+ for(int nu=0;nu<4;++nu){
Xcro(mu,nu) = XCroCont(nu,mu);
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncross(HLV pa, HLV pq,HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncrossW?
HLV q1;
q1=pa;
- for(int i=0;i<nabove+1;i++){
+ for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
double t2 = (q1-pq).m2();
//Blank 3 gamma Current
Tensor<3> J323 = rank3_current(pq,pqbar,hq);
// 2 gamma currents (with 1 contraction already).
Tensor<2> XUncCont = J323.contract((q1-pq),2)/t2;
//Initialise the Uncrossed Vertex
Tensor<2> Xunc(0.);
- for(int mu=0; mu<4;mu++){
- for(int nu=0;nu<4;nu++){
+ for(int mu=0; mu<4;++mu){
+ for(int nu=0;nu<4;++nu){
Xunc(mu,nu) = -XUncCont(mu,nu);
}
}
return Xunc;
}
// Helper Functions Calculate the Eikonal Contributions
Tensor <2> MSym(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
std::vector<HLV> partons, Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MsymW?
HLV q1, q3;
q1=pa;
- for(int i=0;i<nabove+1;i++){
+ for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
q3 = q1-pq-pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
double s23 = (pq+pqbar).m2();
double sa2 = (pa+pq).m2();
double sa3 = (pa+pqbar).m2();
double s12 = (p1+pq).m2();
double s13 = (p1+pqbar).m2();
double sb2 = (pb+pq).m2();
double sb3 = (pb+pqbar).m2();
double s42 = (p4+pq).m2();
double s43 = (p4+pqbar).m2();
Tensor<1> qqxCur = HEJ::current(pq, pqbar, hq);
// // 1a gluon emisson Contribution
Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13))+ pa*(t1/(sa2+sa3)));
Tensor<2> X1aCont = X1a.contract(qqxCur,3);
// //4b gluon emission Contribution
Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43)) + pb*(t3/(sb2+sb3)));
Tensor<2> X4bCont = X4b.contract(qqxCur,3);
// New Formulation Corresponding to New Analytics
Tensor<3> X3g1 = outer(q1+pq+pqbar, metric());
Tensor<3> X3g2 = outer(q3-pq-pqbar, metric());
Tensor<3> X3g3 = outer(q1+q3, metric());
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(qqxCur,3);
Tensor<2> X3g2Cont = X3g2.contract(qqxCur,2);
Tensor<2> X3g3Cont = X3g3.contract(qqxCur,1);
Tensor<2>Xsym(0.);
- for(int mu=0; mu<4;mu++){
- for(int nu=0;nu<4;nu++){
+ for(int mu=0; mu<4;++mu){
+ for(int nu=0;nu<4;++nu){
Xsym(mu, nu) = COM(0,1) * ( (X3g1Cont(nu,mu) + X3g2Cont(mu,nu)
- X3g3Cont(nu,mu)) + (X1aCont(mu,nu) - X4bCont(mu,nu)) );
}
}
return Xsym/s23;
}
//! W+Jets FKL Contributions
/**
* @brief W+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_\mu.
* Handles all possible incoming states.
*/
double jW_j( HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
bool aqlineb, bool /* aqlinef */,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out);
const double WPropfact = WProp(plbar, pl, wprop);
const auto j_W = COM{0,-1}*jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
double Msqr = 0.;
for(const auto h: {plus, minus}) {
const auto j = HEJ::current(p2out, p2in, h);
Msqr += abs2(j_W.contract(j, 1));
}
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*Msqr/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
} // Anonymous Namespace
double ME_W_qQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop);
}
double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true, wprop);
}
double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop);
}
double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true, wprop);
}
double ME_W_qg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param pg Unordered Gluon momenta
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_{uno}_\mu. Ie, unordered with W emission opposite side.
* Handles all possible incoming states.
*/
double jW_juno(HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out,
HLV p2in, HLV pg, bool aqlineb, bool aqlinef,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out-pg);
const HLV q3=-(p2in-p2out);
const Helicity fhel = aqlinef?plus:minus;
const auto j_W = jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
const auto mj2p = HEJ::current(p2out, p2in, flip(fhel));
const auto mj2m = HEJ::current(p2out, p2in, fhel);
const auto jgbp = HEJ::current(pg, p2in, flip(fhel));
const auto jgbm = HEJ::current(pg, p2in, fhel);
const auto j2gp = HEJ::current(p2out, pg, flip(fhel));
const auto j2gm = HEJ::current(p2out, pg, fhel);
// Dot products of these which occur again and again
COM MWmp=j_W.dot(mj2p); // And now for the Higgs ones
COM MWmm=j_W.dot(mj2m);
const auto qsum = to_tensor(q2+q3);
const auto p1o = to_tensor(p1out);
const auto p1i = to_tensor(p1in);
const auto p2o = to_tensor(p2out);
const auto p2i = to_tensor(p2in);
const auto Lmm=( (-1.)*qsum*(MWmm) + (-2.*COM{j_W.dot(pg)})*mj2m + 2.*COM{mj2m.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmm/2. ) )/q3.m2();
const auto Lmp=( (-1.)*qsum*(MWmp) + (-2.*COM{j_W.dot(pg)})*mj2p + 2.*COM{mj2p.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmp/2. ) )/q3.m2();
const auto U1mm=(COM{jgbm.dot(j_W)}*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
const auto U1mp=(COM{jgbp.dot(j_W)}*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
const auto U2mm=((-1.)*COM{j2gm.dot(j_W)}*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
const auto U2mp=((-1.)*COM{j2gp.dot(j_W)}*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
double amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
ampsq*=WProp(plbar, pl, wprop);
return ampsq/((16)*(q2.m2()*q1.m2()));
}
}
double ME_W_unob_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false, wprop);
}
double ME_W_unob_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true, wprop);
}
double ME_W_unob_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false, wprop);
}
double ME_W_unob_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true, wprop);
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (Quark - W and Uno emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (Quark - W and Uno emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
*
* Calculates j_W_{uno} ^\mu j_\mu. Ie, unordered with W emission same side.
* Handles all possible incoming states. Note this handles both forward and back-
* -ward Wuno emission. For forward, ensure p1out is the uno and W emission parton.
* @TODO: Include separate wrapper functions for forward and backward to clean up
* ME_W_unof_current in `MatrixElement.cc`.
*/
double jWuno_j(HLV pg, HLV p1out, HLV plbar, HLV pl, HLV p1in,
HLV p2out, HLV p2in, bool aqlineb,
ParticleProperties const & wprop
){
//Calculate different Helicity choices
const Helicity h = aqlineb?helicity::plus:helicity::minus;
double ME2mpp = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::plus,helicity::plus, wprop);
double ME2mpm = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::plus,helicity::minus, wprop);
double ME2mmp = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::minus,helicity::plus, wprop);
double ME2mmm = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::minus,helicity::minus, wprop);
//Helicity sum and average over initial states
return (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
}
}
double ME_Wuno_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_Wuno_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
/**
* @brief W+Jets Extremal qqx Contributions, function to handle all incoming types.
* @param pgin Incoming gluon which will split into qqx.
* @param pqout Quark of extremal qqx outgoing (W-Emission).
* @param plbar Outgoing anti-lepton momenta
* @param pl Outgoing lepton momenta
* @param pqbarout Anti-quark of extremal qqx pair. (W-Emission)
* @param pout Outgoing Particle 2 (end of FKL chain)
* @param p2in Incoming Particle 2
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W_{qqx} ^\mu j_\mu. Ie, Ex-QQX with W emission same side.
* Handles all possible incoming states. Calculated via crossing symmetry from jWuno_j.
*/
double jWqqx_j(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in, bool aqlinef,
ParticleProperties const & wprop
){
//Calculate Different Helicity Configurations.
const Helicity h = aqlinef?helicity::plus:helicity::minus;
double ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::plus,helicity::plus, wprop);
double ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::plus,helicity::minus, wprop);
double ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::minus,helicity::plus, wprop);
double ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::minus,helicity::minus, wprop);
//Helicity sum and average over initial states.
double ME2 = (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop);
}
double ME_WExqqx_qqbarQ(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop);
}
double ME_WExqqx_qbarqg(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop)
*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_WExqqx_qqbarg(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop)
*K_g(p2out,p2in)/HEJ::C_F;
}
namespace {
//Function to calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm1(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below. (Less Tensor class use?)
double t1 = (p3-pb)*(p3-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p3-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,2)/t1;
return gqqCur*(-1);
}
//Function to calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm2(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double t1 = (p2-pb)*(p2-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb,refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p2-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,1)/t1;
return gqqCur;
}
//Function to calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm3(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s23 = (p2+p3)*(p2+p3);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank1 = outer(p2+p3, metric())/s23;
Tensor<3> qqCurBlank2 = outer(pb, metric())/s23;
Tensor<1> Cur23 = HEJ::current(p2, p3,hel2);
Tensor<2> qqCur1 = qqCurBlank1.contract(Cur23,3);
Tensor<2> qqCur2 = qqCurBlank2.contract(Cur23,3);
Tensor<2> qqCur3 = qqCurBlank2.contract(Cur23,1);
Tensor<1> gqqCur = (qqCur1.contract(epsg,1)
- qqCur2.contract(epsg,2)
+ qqCur3.contract(epsg,1))*2*COM(0,1);
return gqqCur;
}
}
// no wqq emission
double ME_W_Exqqx_QQq(HLV pa, HLV pb, HLV p1, HLV p2,
HLV p3,HLV plbar, HLV pl, bool aqlinepa,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
// 2 independent helicity choices (complex conjugation related).
Tensor<1> TMmmm1 = qggm1(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm2 = qggm2(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm3 = qggm3(pb,p2,p3,minus,minus, pa);
Tensor<1> TMpmm1 = qggm1(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm2 = qggm2(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm3 = qggm3(pb,p2,p3,minus,plus, pa);
// Build the external quark line W Emmision
Tensor<1> cur1a = jW(pa,p1,plbar,pl, aqlinepa?plus:minus);
//Contract with the qqxCurrent.
COM Mmmm1 = TMmmm1.contract(cur1a,1);
COM Mmmm2 = TMmmm2.contract(cur1a,1);
COM Mmmm3 = TMmmm3.contract(cur1a,1);
COM Mpmm1 = TMpmm1.contract(cur1a,1);
COM Mpmm2 = TMpmm2.contract(cur1a,1);
COM Mpmm3 = TMpmm3.contract(cur1a,1);
//Colour factors:
COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
cm1m1=8./3.;
cm2m2=8./3.;
cm3m3=6.;
cm1m2 =-1./3.;
cm1m3 = -3.*COM(0.,1.);
cm2m3 = 3.*COM(0.,1.);
//Sqaure and sum for each helicity config:
double Mmmm = real( cm1m1*pow(abs(Mmmm1),2) + cm2m2*pow(abs(Mmmm2),2)
+ cm3m3*pow(abs(Mmmm3),2) + 2.*real(cm1m2*Mmmm1*conj(Mmmm2))
+ 2.*real(cm1m3*Mmmm1*conj(Mmmm3))
+ 2.*real(cm2m3*Mmmm2*conj(Mmmm3)) );
double Mpmm = real( cm1m1*pow(abs(Mpmm1),2) + cm2m2*pow(abs(Mpmm2),2)
+ cm3m3*pow(abs(Mpmm3),2) + 2.*real(cm1m2*Mpmm1*conj(Mpmm2))
+ 2.*real(cm1m3*Mpmm1*conj(Mpmm3))
+ 2.*real(cm2m3*Mpmm2*conj(Mpmm3)) );
// Divide by WProp
const double WPropfact = WProp(plbar, pl, wprop);
return (2*WPropfact*(Mmmm+Mpmm)/24./4.)/(pa-p1-pl-plbar).m2()/(p2+p3-pb).m2();
}
// W+Jets qqxCentral
double ME_WCenqqx_qq(HLV pa, HLV pb,HLV pl, HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
ParticleProperties const & wprop
){
init_sigma_index();
HLV pq, pqbar, p1, p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am, T4bm, T1ap, T4bp;
if(!(aqlinepa)){
T1ap = HEJ::current(p1, pa, helicity::plus);
T1am = HEJ::current(p1, pa, helicity::minus);}
else if(aqlinepa){
T1ap = HEJ::current(pa, p1, helicity::plus);
T1am = HEJ::current(pa, p1, helicity::minus);}
if(!(aqlinepb)){
T4bp = HEJ::current(p4, pb, helicity::plus);
T4bm = HEJ::current(p4, pb, helicity::minus);}
else if(aqlinepb){
T4bp = HEJ::current(pb, p4, helicity::plus);
T4bm = HEJ::current(pb, p4, helicity::minus);}
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc = MUncrossW(pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xcro = MCrossW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xsym = MSymW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is
// also the choice in qqbar helicity.
// (- - hel choice)
COM M_mmUnc = (((Xunc).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc).contract(T1am,1)).contract(T4bp,1));
COM M_mpCro = (((Xcro).contract(T1am,1)).contract(T4bp,1));
COM M_mpSym = (((Xsym).contract(T1am,1)).contract(T4bp,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc).contract(T1ap,1)).contract(T4bp,1));
COM M_ppCro = (((Xcro).contract(T1ap,1)).contract(T4bp,1));
COM M_ppSym = (((Xsym).contract(T1ap,1)).contract(T4bp,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
- for(int i=0;i<nabove+1;i++){
+ for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
return amp;
}
// no wqq emission
double ME_W_Cenqqx_qq(HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
int nbelow, bool forwards, ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
if (!forwards){ //If Emission from Leg a instead, flip process.
std::swap(pa, pb);
std::reverse(partons.begin(),partons.end());
std::swap(aqlinepa, aqlinepb);
qqxmarker = !qqxmarker;
std::swap(nabove,nbelow);
}
HLV pq, pqbar, p1,p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am(0.), T1ap(0.);
if(!(aqlinepa)){
T1ap = HEJ::current(p1, pa, plus);
T1am = HEJ::current(p1, pa, minus);}
else if(aqlinepa){
T1ap = HEJ::current(pa, p1, plus);
T1am = HEJ::current(pa, p1, minus);}
Tensor <1> T4bm = jW(pb, p4, plbar, pl, aqlinepb?plus:minus);
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc_m = MUncross(pa, pq, pqbar,partons, minus, nabove);
Tensor<2> Xcro_m = MCross( pa, pq, pqbar,partons, minus, nabove);
Tensor<2> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, minus, nabove);
Tensor<2> Xunc_p = MUncross(pa, pq, pqbar,partons, plus, nabove);
Tensor<2> Xcro_p = MCross( pa, pq, pqbar,partons, plus, nabove);
Tensor<2> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, plus, nabove);
// (- - hel choice)
COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
- for(int i=0;i<nabove+1;i++){
+ for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
return amp;
}