Commit ce88143d authored by Enrico Bothmann's avatar Enrico Bothmann Committed by Sebastian Liebschner

Remove obsolete EXTAMP/Main

All sources in this directory have been superseded by sources that are
located directly under EXTAMP/
parent fdaa2a79
This diff is collapsed.
#ifndef EXTAMP_BVI_Process_H
#define EXTAMP_BVI_Process_H
#include "EXTAMP/Main/Process.H"
namespace PHASIC {
class Virtual_ME2_Base;
class Color_Correlated_ME2;
class Process_Info;
class KP_Terms;
namespace EXTAMP {
class BVI_Process : public Process {
public :
BVI_Process(const PHASIC::Process_Info&);
double Partonic(const ATOOLS::Vec4D_Vector &p,
const int mode);
void SetNLOMC(PDF::NLOMC_Base *const mc);
private :
/* Subtraction scheme: O: plain CS
1: Dire
2: modfied CS */
int m_subtype;
PHASIC::Virtual_ME2_Base* p_loop_me;
PHASIC::Color_Correlated_ME2* p_corr_me;
PHASIC::KP_Terms* p_kpterms;
double Calc_I(const ATOOLS::Vec4D_Vector& p,
const double& mur) const;
/* Can't be const because we need to store momentum fractions
m_x0,m_x1,m_eta0,m_eta1 of KP_Terms for reweighting */
double Calc_KP(const ATOOLS::Vec4D_Vector& p);
double m_x0, m_x1, m_eta0, m_eta1;
/* Used by PHASIC::Single_Process for reweighting op KP terms */
double KPTerms(int mode, double scalefac2=1.0);
/* Calc all terms explicitly dependent on renormalization scale
for on-the-fly reweighting. First item: first derivative of all
terms with respect to scale log, second item: second derivative
of all terms with respect to scale log. */
Calc_ScaleDependenceTerms(const ATOOLS::Vec4D_Vector& p,
const double& mur) const;
/* epsilon^0, epsilon^{-1}, epsilon^{-2} coefficients of eq.
(5.32) through (5.34) of arXiv:hep-ph/9605323 */
static double Vi_eps0(const ATOOLS::Flavour& fl, int subtype);
static double Vi_eps1(const ATOOLS::Flavour& fl);
static double Vi_eps2(const ATOOLS::Flavour& fl);
constexpr static double m_CF = 4./3.;
constexpr static double m_CA = 3.0;
constexpr static double m_TR = 1./2.;
static double m_NF;
/* beta0 = 11/3 CA - 4/3 TR nf*/
double m_beta0;
static double Ti2(const ATOOLS::Flavour& fl);
#include "EXTAMP/Main/External_ME_Interface.H"
#include "EXTAMP/Main/Born_Process.H"
#include "PHASIC++/Process/Tree_ME2_Base.H"
namespace EXTAMP {
Born_Process::Born_Process(const PHASIC::Process_Info& pi) : Process(pi)
p_born_me = External_ME_Interface::GetExternalBornME(pi);
double Born_Process::Partonic(const ATOOLS::Vec4D_Vector &p,
const int mode)
/* Maybe move to PHASIC::Single_Process */
double dxs = p_born_me->Calc(p)/NormFac();
/* Single_Process derivatives are responsible for storing the
return value in m_lastdxs and for filling the m_mewgtinfo
inherited from PHASIC::Process_Base */
m_mewgtinfo.m_K = 1.0;
m_mewgtinfo.m_B = dxs;
m_lastxs = dxs;
return dxs;
#ifndef EXTAMP_Born_Process_H
#define EXTAMP_Born_Process_H
#include "EXTAMP/Main/Process.H"
namespace PHASIC {
class Tree_ME2_Base;
class Process_Info;
namespace EXTAMP {
class Born_Process : public Process {
public :
Born_Process(const PHASIC::Process_Info& pi);
~Born_Process() {};
double Partonic(const ATOOLS::Vec4D_Vector &p,
const int mode);
PHASIC::Tree_ME2_Base* p_born_me;
#include "EXTAMP/Main/Dipole_Wrapper_Process.H"
#include "EXTAMP/Main/External_ME_Interface.H"
#include "PHASIC++/Selectors/Combined_Selector.H"
#include "PHASIC++/Main/Process_Integrator.H"
#include "PDF/Main/NLOMC_Base.H"
#include "ATOOLS/Org/MyStrStream.H"
#include <assert.h>
using namespace EXTAMP;
Dipole_Wrapper_Process::Dipole_Wrapper_Process(const RS_Process& rsproc,
PHASIC::CS_Dipole* dipole,
BEAM::Beam_Spectra_Handler* beam,
PDF::ISR_Handler* isr) : p_dipole(dipole)
/* Follow sherpa convetions here and set i<j */
m_norm = rsproc.NormFac();
m_born_procinfo = ConstructBornProcessInfo(rsproc.Info(),I(),J(),
m_born_flavs = m_born_procinfo.ExtractFlavours();
/* This wrapper has to call it's Init method itself, since it's not
initialized through the framework like other procs. This
re-orders flavours according to Sherpa conventions (reordering is
also applied to the Process_Info instance). */
PHASIC::Process_Base::Init(rsproc.Info(), beam, isr, 0);
/* Not done in any of the PHASIC base classes */
m_mincpl = rsproc.Info().m_mincpl;
m_maxcpl = rsproc.Info().m_maxcpl;
/* Construct an index maps such that
Flavours[i] == Dipole()->Flavours[m_indexmap[i]]
Flavours[m_inversemap] == Dipole()->Flavours[i] */
m_indexmap = ConstructIndexMapping(Dipole()->Flavours(),
m_inversemap = InvertIndexMapping(m_indexmap);
/* Cross-check index mapping */
for(size_t i(0); i<BornFlavours().size(); i++)
assert(BornFlavours()[i] == dipole->Flavours()[m_indexmap[i]]);
assert(BornFlavours()[m_inversemap[i]] == dipole->Flavours()[i]);
/* Set name in accordance with Sherpa conventions so that it can be
identified by MC@NLO process. Save name of born config as well
for NLO_subevts. */
m_born_name = PHASIC::Process_Base::GenerateName(m_born_procinfo.m_ii,
m_name = rsproc.Name()+"_RS"+
/* Set couplings: copy pointers from dipole */
/* Fill Comninable Info */
= External_ME_Interface::ConstructCombinableMap(Flavours(),Info(),NIn());
m_id_vector = ConstructIDVector();
PHASIC::Process_Info Dipole_Wrapper_Process::ConstructBornProcessInfo
(const PHASIC::Process_Info& rsinfo,
size_t i, size_t j, const ATOOLS::Flavour& flav_ij)
PHASIC::Process_Info ret(rsinfo);
return ret;
std::vector<size_t> Dipole_Wrapper_Process::ConstructIndexMapping
(const ATOOLS::Flavour_Vector& dipole_flavs,
const ATOOLS::Flavour_Vector& process_flavs,
size_t nin)
/* Construct Process_Info with same ordering as in dipole_flavs */
PHASIC::Subprocess_Info in;
PHASIC::Subprocess_Info fi;
for(size_t i(0); i<nin; i++)
for(size_t i(nin); i<dipole_flavs.size(); i++)
PHASIC::Process_Info pi(in,fi);
assert(pi.ExtractFlavours() == dipole_flavs);
/* Tag all flavours */
for(size_t i(0); i<nin; i++)
int cnt(nin); pi.m_fi.SetTags(cnt);
/* Apply Sherpa's re-ordering conventions */
assert(pi.ExtractFlavours() == process_flavs);
/* Extract map from tags and re-ordered flavours */
std::vector<size_t> ret;
for (size_t i(0);i<nin;++i)
std::vector<int> fi_tags;
size_t nout = dipole_flavs.size()-nin;
for (size_t i(0);i<nout;++i) ret[nin+i]=fi_tags[i];
return ret;
std::vector<size_t> Dipole_Wrapper_Process::InvertIndexMapping
(const std::vector<size_t>& map)
std::vector<size_t> ret(map.size());
for(size_t i(0); i<map.size(); i++)
ret[map[i]] = i;
return ret;
void Dipole_Wrapper_Process::SetSubEventProperties(ATOOLS::NLO_subevt& sub)
sub.p_fl = &BornFlavours()[0];
sub.p_mom = &Momenta()[0];
sub.m_n = BornFlavours().size();
sub.m_i = I();
sub.m_j = J();
sub.m_k = K();
sub.m_kt = BornK();
sub.m_ijt = BornIJ();
sub.p_id = &IDVector()[0];
sub.m_me = 0.0;
sub.m_result = 0.0;
sub.m_trig = false;
sub.p_proc = this;
sub.m_pname = m_born_name;
sub.m_pname = sub.m_pname.substr(0,sub.m_pname.rfind("__"));
void Dipole_Wrapper_Process::AssignSubEvent(ATOOLS::NLO_subevt* sub)
p_subevent = sub;
std::vector<size_t> Dipole_Wrapper_Process::ConstructIDVector() const
std::vector<size_t> ret(m_indexmap);
for(size_t i(0); i<ret.size(); i++) ret[i] = 1<<ret[i];
ret[BornIJ()] = (1<<I())|(1<<J());
return ret;
void Dipole_Wrapper_Process::FillProcessMap(PHASIC::NLOTypeStringProcessMap_Map *apmap)
/* Copied from AMEGIC::Single_DipoleTerm */
if (p_apmap->find(ATOOLS::nlo_type::rsub)==p_apmap->end())
(*p_apmap)[ATOOLS::nlo_type::rsub] = new PHASIC::StringProcess_Map();
std::string fname(m_name);
size_t len(m_pinfo.m_addname.length());
if (len) fname=fname.erase(fname.rfind(m_pinfo.m_addname),len);
bool Dipole_Wrapper_Process::Combinable(const size_t &idi,const size_t &idj)
bool ret = m_cluster_flav_map.find(idi | idj)!=m_cluster_flav_map.end();
return ret;
const ATOOLS::Flavour_Vector &Dipole_Wrapper_Process::CombinedFlavour(const size_t &idij)
std::map<size_t, ATOOLS::Flavour_Vector>::const_iterator it = m_cluster_flav_map.find(idij);
THROW(fatal_error, "Internal error");
return it->second;
void Dipole_Wrapper_Process::SetNLOMC(PDF::NLOMC_Base *const nlomc)
void Dipole_Wrapper_Process::CalcKinematics(const ATOOLS::Vec4D_Vector& p)
if(!Dipole()) THROW(fatal_error, "Invalid dipole pointer");
/* Dipole and dipole wrapper share the same flavour and momentum
ordering in the real emission configuration, pass momentum on
as-is */
/* Apply re-mapping of (Born-) momenta, invert incoming momenta for
NLO_subevts */
for(size_t i(0); i<NIn(); i++)
m_moms[i] = -(Dipole()->Momenta()[m_indexmap[i]]);
for(size_t i(NIn()); i<m_moms.size(); i++)
m_moms[i] = (Dipole()->Momenta()[m_indexmap[i]]);
double Dipole_Wrapper_Process::Calc(ATOOLS::NLO_subevt* const evt)
double sign = MCModeSign(evt);
double dxs = (sign==0.0) ? 0.0 : sign*Dipole()->Calc()/NormFac();
/* NLO_subevts of RS proc are added up by PHASIC::Single_Process,
hence add a minus sign here */
evt->m_me = -dxs;
evt->m_mewgt = -dxs;
evt->m_result = -dxs;
return dxs;
void Dipole_Wrapper_Process::CalculateScale(const ATOOLS::Vec4D_Vector& real_p,
const ATOOLS::Vec4D_Vector& born_p,
ATOOLS::NLO_subevt* const evt)
/* The scale setter knows about the RS and all other dipoles via the
subevtlist of the RS. Scale setting might be based on the real
configuration so we have to set the momenta of the RS here before
calling the scale setter. */
RS_Process* rsproc = ((RS_Process*)(evt->p_real->p_proc));
/* Locally set flavours to born configuration in order not to
confuse the scale setter */
ATOOLS::Flavour_Vector tmp = m_flavs;
m_flavs = m_born_flavs;
m_flavs = tmp;
double Dipole_Wrapper_Process::Partonic(const ATOOLS::Vec4D_Vector &p, const int mode)
/* This method is called in MC@NLO matching from outside an
RS_Process. Hence have to calc scales ourselves */
CalculateScale(p, m_moms, p_subevent);
double dxs = Calc(p_subevent);
return m_lastxs = m_mewgtinfo.m_B = dxs;
int Dipole_Wrapper_Process::MCModeSign(ATOOLS::NLO_subevt* const evt) const
/* m_mcmode is nonzero if this process is used to calcualte the
DA-DS term or the DA term as defined in arXiv:1111.1200. Since
DA=DS apart from theta functions implementing the PS starting
scale in DA, we just have to implement a sign \in {-1,0,1}
here. */
double kt2 = GetKT2ofSplitting(*Dipole()->LastKinematics());
double kt2max = GetMaxKT2ForDA();
evt->m_kt2 = kt2;
int DS(1),DA(kt2<=kt2max*(1.0+1e-6));
/* For the calculation of DA-DS, another minus sign is added in
PHASIC::Single_Process. Need to compensate for that here by yet
another minus sign. */
if (m_mcmode==0) return + DS;
if (m_mcmode==1) return -(DA-DS);
if (m_mcmode==2) return + DA;
THROW(fatal_error, "Unknown MC-mode "+ATOOLS::ToString(m_mcmode));
double Dipole_Wrapper_Process::GetKT2ofSplitting(const PHASIC::Dipole_Kinematics& kin) const
/* Get hardness of splitting as defined by shower. Only needed for
MC@NLO. */
if(!p_nlomc) return 0.0;
double x = kin.ShowerX() ;
double y = kin.ShowerY() ;
double Q2 = kin.ShowerQ2();
return p_nlomc->KT2(*p_subevent,x,y,Q2);
double Dipole_Wrapper_Process::GetMaxKT2ForDA() const
/* Get the maximum hardness for a splitting: PS starting scale */
return p_subevent->m_mu2[(ATOOLS::stp::id(ATOOLS::stp::size+ATOOLS::stp::res))];
return p_subevent->m_mu2[(ATOOLS::stp::res)];
#ifndef EXTAMP_Main_Dipole_Wrapper_Process_H
#define EXTAMP_Main_Dipole_Wrapper_Process_H
#include "PHASIC++/Process/CS_Dipole.H"
#include "PHASIC++/Process/Single_Process.H"
#include "EXTAMP/Main/RS_Process.H"
namespace EXTAMP {
A wrapper around a PHASIC::CS_Dipole that basically adds zero
functionality to it, apart from making it look like a
PHASIC::Single_Process with a 'Partonic' method that returns the
value of the dipole. This is necessary in the MC@NLO
implementation due to AMEGIC's legacy structure of representing
dipole terms by processes.
class Dipole_Wrapper_Process : public PHASIC::Single_Process {
Dipole_Wrapper_Process(const RS_Process& rsproc,
PHASIC::CS_Dipole* dipole,
BEAM::Beam_Spectra_Handler* beam,
PDF::ISR_Handler* isr);
/* Calc dipole kinematics without evaluating dipole. Also perform
trigger and set corresponding flag of NLO_subevt */
void CalcKinematics(const ATOOLS::Vec4D_Vector& p);
/* After having called CalcKinematics(), evaluate the dipole. */
double Calc(ATOOLS::NLO_subevt* sub);
void CalculateScale(const ATOOLS::Vec4D_Vector& real_p,
const ATOOLS::Vec4D_Vector& born_p,
ATOOLS::NLO_subevt* const evt);
/* Inherited from PHASIC::Single_Process. Calls CalcKinematics and
Calc, then returns the result. This is used in by Sherpa's
MC@NLO implementation for calculating DA and DA-DS */
double Partonic(const ATOOLS::Vec4D_Vector &p,
const int mode);
/* A sign implementing DA and DA-DS as defined in arXiv:1111.1200 */
int MCModeSign(ATOOLS::NLO_subevt* const evt) const;
double GetMaxKT2ForDA() const;
double GetKT2ofSplitting(const PHASIC::Dipole_Kinematics& kin) const;
bool Combinable(const size_t &idi,
const size_t &idj);
void FillProcessMap(PHASIC::NLOTypeStringProcessMap_Map *apmap);
const ATOOLS::Flavour_Vector &CombinedFlavour(const size_t &idij);
/* Set properties of NLO_subevt according to this dipole */
void SetSubEventProperties(ATOOLS::NLO_subevt& sub);
/* Set a member pointer to the subevent passed */
void AssignSubEvent(ATOOLS::NLO_subevt* sub);
void SetScaleSetter(PHASIC::Scale_Setter_Base* scl) { p_scale = scl; }
void SetNLOMC(PDF::NLOMC_Base *const nlomc);
/* Set p_scale to NULL in order to avoid deleting it in destructor
of PHASIC::Process_Base. It gets deleted in the RS_Process,
which owns it. */
~Dipole_Wrapper_Process() { p_scale = NULL; }
/* Flavours of this process have to be those of real emission
config. NLO_subevts are assigned the born configuration,
however. Store corresponding information in these members: */
std::string m_born_name;
PHASIC::Process_Info m_born_procinfo;
ATOOLS::Flavour_Vector m_born_flavs;
const ATOOLS::Flavour_Vector& BornFlavours() const { return m_born_flavs; }
static PHASIC::Process_Info ConstructBornProcessInfo(const PHASIC::Process_Info& rsinfo,
size_t i, size_t j,
const ATOOLS::Flavour& flav_ij) ;
/* This pointer (passed in constructor) is assumed to be valid
throughout, associated memory is not managed. */
PHASIC::CS_Dipole* p_dipole;
PHASIC::CS_Dipole* Dipole() const { return p_dipole; }
/* RS_Process is supposed to initialize and manages the
subevent. */
ATOOLS::NLO_subevt* p_subevent;
/* Real emission momentum configuration with momenta of initial
state particles reversed. This convention is used in
NLO_subevts. */
ATOOLS::Vec4D_Vector m_moms;
const ATOOLS::Vec4D_Vector& Momenta() const { return m_moms; }
/* Emitter, emitted, spectator index in real emission flavour
vector. Follow Sherpa conventions: i<j */
size_t I() const { return std::min(Dipole()->I(),Dipole()->J()); }
size_t J() const { return std::max(Dipole()->I(),Dipole()->J()); }
size_t K() const { return Dipole()->K(); }
/* Position of combined flavour ij and spectator k in born
configuration */
size_t BornIJ() const { return m_inversemap[Dipole()->BornIJ()]; }
size_t BornK() const { return m_inversemap[Dipole()->BornK() ]; }
/* Same as in EXTAMP::Process */
double m_norm;
const double& NormFac() const { return m_norm; }
/* Mapping from ordering in p_dipole to ordering in this class
(has to comply with Sherpa ordering conventions):
Flavours[i] == Dipole()->Flavours[m_indexmap[i]]
Flavours[m_inversemap] == Dipole()->Flavours[i] */
std::vector<size_t> m_indexmap;
std::vector<size_t> m_inversemap;
static std::vector<size_t> ConstructIndexMapping(const ATOOLS::Flavour_Vector& dipole_flavs,
const ATOOLS::Flavour_Vector& process_flavs,
size_t nin);
static std::vector<size_t> InvertIndexMapping(const std::vector<size_t>& map);
/* ID vector encoding the id's of particles as they are ordered in
the born flavour vector. This is needed for
ATOOLS::NLO_Subevents, which own a pointer to such a vector.
id's are in binary encoding, meaning the i'th particle has id
1<<i and the combined particles i and j have inded
(1<<i|1<<j). */
std::vector<size_t> m_id_vector;
const std::vector<size_t>& IDVector() const {return m_id_vector; }
std::vector<size_t> ConstructIDVector() const;
std::map<size_t, ATOOLS::Flavour_Vector> m_cluster_flav_map;
#include "PDF/Main/Cluster_Definitions_Base.H"
#include "PHASIC++/Process/Tree_ME2_Base.H"
#include "PHASIC++/Process/ME_Generator_Base.H"
#include "PHASIC++/Process/Process_Group.H"
#include "PHASIC++/Process/Subprocess_Info.H"
#include "PHASIC++/Process/CS_Dipole.H"
#include "EXTAMP/Main/External_ME_Interface.H"
#include "EXTAMP/Main/Process_Group.H"
#include "EXTAMP/Main/Process.H"
#include "EXTAMP/Main/Born_Process.H"
#include "EXTAMP/Main/RS_Process.H"
#include "EXTAMP/Main/BVI_Process.H"
#include <assert.h>
namespace EXTAMP{
External_ME_Interface::External_ME_Interface() :
bool External_ME_Interface::Initialize(MODEL::Model_Base *const model,
BEAM::Beam_Spectra_Handler *const beam,
PDF::ISR_Handler *const isr)
/* Pure virtual in PHASIC::ME_Generator_Base. Store beam, isr, and
model since they need to be passed down to any
PHASIC::Process_Base when calling its init method
PHASIC::Process_Base::Init(...) */
p_beam = beam;
p_isr = isr;
return true;
External_ME_Interface::~External_ME_Interface() {}
PHASIC::Process_Base* External_ME_Interface::InitializeProcess(
const PHASIC::Process_Info &pi, bool add)
/* This is pure virtual in ME_Generator_Base and supposed to
instantiate a high-level group-type process deriving from
PHASIC::Process_Group. This inherits from PHASIC::Process_Base,
which requires PHASIC::Process_Base::Init(...) to be called.
Then, PHASIC::Process_Group::ConstructProcesses() is used in
Sherpa to (recursively) initialize the PARTONIC sub-processes
of the group. Within
PHASIC::Process_Group::ConstructProcesses(), several pure
virtuals of PHASIC::Process_Group are used for that purpose:
- Process_Base* PHASIC::Process_Group::GetProcess(const PHASIC::Process_Info &pi)
- bool PHASIC::Process_Group::Initialize(Process_Base *const proc)
The former must be implemented by EXTAMP::Process_Group in such
a way as to instantiate the PARTONIC (i.e. non-group type)
process corresponding to the Process_Info, and return a valid
pointer to it. Does not make any sense structurally, so just
implemented it as a wrapper around
External_ME_Interface::InstantiatePartonicProcess(PHASIC::Process_Info). */
EXTAMP::Process_Group *newproc = new EXTAMP::Process_Group();
/* In case no valid partonic channels are found, return NULL so the
Matrix_Element_Hanlder knows and throws a 'No hard process found' */
return (newproc->Size()>0) ? newproc : NULL;
/* Determine if a PARTONIC born Process as specified by pi exists */
bool External_ME_Interface::PartonicProcessExists(PHASIC::Process_Info &pi)
if (PHASIC::Tree_ME2_Base::GetME2(pi)!=NULL) return true;
/* Could be a real correction, increase order and try again. */
pi.m_mincpl[0]+= 1;
return PHASIC::Tree_ME2_Base::GetME2(pi)!=NULL;
PHASIC::Tree_ME2_Base* External_ME_Interface::GetExternalBornME(const PHASIC::Process_Info& pi)
/* Use the getter to find an implementation of born matrix element */
PHASIC::Process_Info cpi(pi);
PHASIC::Tree_ME2_Base* me2 = PHASIC::Tree_ME2_Base::GetME2(cpi);
if(me2) return me2;
/* Could be a real correction! */
cpi.m_mincpl[0]+= 1;
me2 = PHASIC::Tree_ME2_Base::GetME2(cpi);
if(me2) return me2;
std::stringstream msg; msg << "No ME found for this process:\n" << pi;
THROW(fatal_error, msg.str());
PHASIC::CS_Dipole* External_ME_Interface::GetExternalDipole(const PHASIC::Dipole_Info& di)
/* Use the getter to find implementation of correla