Page MenuHomeHEPForge
Files F11221299

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Gamma/MEGammaP2Jets.cc b/MatrixElement/Gamma/MEGammaP2Jets.cc
--- a/MatrixElement/Gamma/MEGammaP2Jets.cc
+++ b/MatrixElement/Gamma/MEGammaP2Jets.cc
@@ -1,400 +1,400 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEGammaP2Jets class.
//
#include "MEGammaP2Jets.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
using namespace Herwig;
MEGammaP2Jets::MEGammaP2Jets() : _process(0), _minflavour(1) ,_maxflavour(5) {
massOption(vector<unsigned int>(2,0));
}
unsigned int MEGammaP2Jets::orderInAlphaS() const {
return 1;
}
unsigned int MEGammaP2Jets::orderInAlphaEW() const {
return 1;
}
Energy2 MEGammaP2Jets::scale() const {
return 2.*sHat()*tHat()*uHat()/
(sqr(sHat())+sqr(tHat())+sqr(uHat()));
}
void MEGammaP2Jets::persistentOutput(PersistentOStream & os) const {
os << _gluonvertex << _photonvertex
<< _process << _minflavour << _maxflavour;
}
void MEGammaP2Jets::persistentInput(PersistentIStream & is, int) {
is >> _gluonvertex >> _photonvertex
>> _process >> _minflavour >> _maxflavour;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEGammaP2Jets,HwMEBase>
describeHerwigMEGammaP2Jets("Herwig::MEGammaP2Jets", "HwMEGammaHadron.so");
void MEGammaP2Jets::Init() {
static ClassDocumentation<MEGammaP2Jets> documentation
("The MEGammaP2Jets class implements the matrix elements"
" for pointlike photon-hadron to jets.");
static Parameter<MEGammaP2Jets,int> interfaceMinimumFlavour
("MinimumFlavour",
"The minimum flavour of the quarks",
&MEGammaP2Jets::_minflavour, 1, 1, 5,
false, false, Interface::limited);
static Parameter<MEGammaP2Jets,int> interfaceMaximumFlavour
("MaximumFlavour",
"The maximum flavour of the quarks",
&MEGammaP2Jets::_maxflavour, 5, 1, 5,
false, false, Interface::limited);
static Switch<MEGammaP2Jets,unsigned int> interfaceProcess
("Process",
"The allowed partonic subprocesses",
&MEGammaP2Jets::_process, 0, false, false);
static SwitchOption interfaceProcessAll
(interfaceProcess,
"All",
"Include all the subprocesses",
0);
static SwitchOption interfaceProcessGluon
(interfaceProcess,
"Gluon",
"Only include the gamma g -> q qbar processes",
1);
static SwitchOption interfaceProcessQuark
(interfaceProcess,
"Quark",
"Only include the gamma q -> gluon q processes",
2);
static SwitchOption interfaceProcessAntiQuark
(interfaceProcess,
"AntiQuark",
"Only include the gamma qbar -> gluon qbar processes",
3);
}
void MEGammaP2Jets::getDiagrams() const {
tcPDPtr g = getParticleData(ParticleID::g);
tcPDPtr gm = getParticleData(ParticleID::gamma);
for ( int i = _minflavour; i <= _maxflavour; ++i ) {
tcPDPtr q = getParticleData(i);
tcPDPtr qb = q->CC();
// gamma g -> q qbar
if(_process==0||_process==1) {
add(new_ptr((Tree2toNDiagram(3), gm, q , g, 1, q, 2, qb, -1)));
add(new_ptr((Tree2toNDiagram(3), gm, qb, g, 2, q, 1, qb, -2)));
}
// gamma q -> g q
if(_process==0||_process==2) {
add(new_ptr((Tree2toNDiagram(3), gm, q , q, 2, g, 1, q, -3)));
add(new_ptr((Tree2toNDiagram(2), gm, q , 1, q, 3, g, 3, q, -4)));
}
// gamma qbar -> g qbar
if(_process==0||_process==3) {
add(new_ptr((Tree2toNDiagram(3), gm, qb, qb, 2, g, 1, qb, -5)));
add(new_ptr((Tree2toNDiagram(2), gm, qb, 1, qb, 3, g, 3, qb, -6)));
}
}
}
Selector<MEBase::DiagramIndex>
MEGammaP2Jets::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
if ( diags[i]->id() == -1 ) sel.insert(meInfo()[0], i);
else if ( diags[i]->id() == -2 ) sel.insert(meInfo()[1], i);
else if ( diags[i]->id() == -3 ) sel.insert(meInfo()[0], i);
else if ( diags[i]->id() == -4 ) sel.insert(meInfo()[1], i);
else if ( diags[i]->id() == -5 ) sel.insert(meInfo()[0], i);
else if ( diags[i]->id() == -6 ) sel.insert(meInfo()[1], i);
}
return sel;
}
IBPtr MEGammaP2Jets::clone() const {
return new_ptr(*this);
}
IBPtr MEGammaP2Jets::fullclone() const {
return new_ptr(*this);
}
void MEGammaP2Jets::doinit() {
// get the vedrtex pointers from the SM object
tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
// do the initialisation
if(hwsm) {
_gluonvertex = hwsm->vertexFFG();
_photonvertex = hwsm->vertexFFP();
}
else throw InitException() << "Wrong type of StandardModel object in "
<< "MEGammaP2Jets::doinit() the Herwig"
<< " version must be used"
<< Exception::runerror;
// call the base class
HwMEBase::doinit();
}
Selector<const ColourLines *>
MEGammaP2Jets::colourGeometries(tcDiagPtr diag) const {
static ColourLines c1("3 2 4,-3 -5");
static ColourLines c2("3 4,-3 -2 -5");
static ColourLines c3("3 4,-4 2 5");
static ColourLines c4("2 3 4,-4 5");
static ColourLines c5("-3 -4,4 -2 -5");
static ColourLines c6("-2 -3 -4,4 -5");
Selector<const ColourLines *> sel;
if ( diag->id() == -1 ) sel.insert(1.0, &c1);
else if ( diag->id() == -2 ) sel.insert(1.0, &c2);
else if ( diag->id() == -3 ) sel.insert(1.0, &c3);
else if ( diag->id() == -4 ) sel.insert(1.0, &c4);
else if ( diag->id() == -5 ) sel.insert(1.0, &c5);
else if ( diag->id() == -6 ) sel.insert(1.0, &c6);
return sel;
}
double MEGammaP2Jets::me2() const {
// total matrix element and the various components
double me(0.);
// gamma g -> q qbar
if(mePartonData()[1]->id()==ParticleID::g) {
VectorWaveFunction gmin (meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction glin (meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction qout (meMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction qbout(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> v1,v2;
vector<SpinorBarWaveFunction> a3;
vector<SpinorWaveFunction> f4;
for(unsigned int ix=0;ix<2;++ix) {
gmin.reset(2*ix);
v1.push_back(gmin);
glin.reset(2*ix);
v2.push_back(glin);
qout.reset(ix);
a3.push_back(qout);
qbout.reset(ix);
f4.push_back(qbout);
}
// calculate the matrix element
me = gammagluonME(v1,v2,a3,f4,false);
}
// gamma q -> g q
else if(mePartonData()[1]->id()>0) {
VectorWaveFunction gmin (meMomenta()[0],mePartonData()[0],incoming);
SpinorWaveFunction qin (meMomenta()[1],mePartonData()[1],incoming);
VectorWaveFunction glout(meMomenta()[2],mePartonData()[2],outgoing);
SpinorBarWaveFunction qout (meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> v1,v3;
vector<SpinorWaveFunction> f2;
vector<SpinorBarWaveFunction> f4;
for(unsigned int ix=0;ix<2;++ix) {
gmin.reset(2*ix);
v1.push_back(gmin);
qin.reset(ix);
f2.push_back(qin);
glout.reset(2*ix);
v3.push_back(glout);
qout.reset(ix);
f4.push_back(qout);
}
me = gammaquarkME(v1,f2,v3,f4,false);
}
// gamma qbar -> g qbar
else {
VectorWaveFunction gmin (meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction qin (meMomenta()[1],mePartonData()[1],incoming);
VectorWaveFunction glout(meMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction qout (meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> v1,v3;
vector<SpinorBarWaveFunction> a2;
vector<SpinorWaveFunction> a4;
for(unsigned int ix=0;ix<2;++ix) {
gmin.reset(2*ix);
v1.push_back(gmin);
qin.reset(ix);
a2.push_back(qin);
glout.reset(2*ix);
v3.push_back(glout);
qout.reset(ix);
a4.push_back(qout);
}
me = gammaantiquarkME(v1,a2,v3,a4,false);
}
return me;
}
double MEGammaP2Jets::gammagluonME(vector<VectorWaveFunction> & gmin,
vector<VectorWaveFunction> & glin,
vector<SpinorBarWaveFunction> & fout,
vector<SpinorWaveFunction> & aout,
bool calc) const {
// the scale
Energy2 mt(scale());
// storage of the helicity me if needed
ProductionMatrixElement newme(PDT::Spin1,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1Half);
// loop over the helicities
SpinorWaveFunction inter;
vector<double> me(3,0.);
vector<Complex> diag(2);
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
- inter = _gluonvertex->evaluate(mt,5,aout[ohel2].particle(),
+ inter = _gluonvertex->evaluate(mt,5,aout[ohel2].particle()->CC(),
aout[ohel2],glin[ihel2]);
diag[0] = _photonvertex->evaluate(0.*GeV2,inter,fout[ohel1],
gmin[ihel1]);
// second diagram
- inter = _photonvertex->evaluate(0.*GeV2,5,aout[ohel2].particle(),
+ inter = _photonvertex->evaluate(0.*GeV2,5,aout[ohel2].particle()->CC(),
aout[ohel2],gmin[ihel1]);
diag[1] = _gluonvertex->evaluate(mt,inter,fout[ohel1],glin[ihel2]);
for(unsigned int ix=0;ix<2;++ix) me[ix] += norm(diag[ix]);
diag[0]+=diag[1];
me[2] += norm(diag[0]);
// matrix element if needed
if(calc) newme(2*ihel1,2*ihel2,ohel1,ohel2)=diag[0];
}
}
}
}
// save the info on the diagrams
if(!calc) {
DVector save;
save.push_back(me[0]);
save.push_back(me[1]);
meInfo(save);
}
// return the answer with colour and spin factors
if(calc) _me.reset(newme);
return 0.125*me[2];
}
double MEGammaP2Jets::gammaquarkME(vector<VectorWaveFunction> & gmin,
vector<SpinorWaveFunction> & fin,
vector<VectorWaveFunction> & gout,
vector<SpinorBarWaveFunction> & fout,
bool calc) const {
// the scale
Energy2 mt(scale());
// storage of the helicity me if needed
ProductionMatrixElement newme(PDT::Spin1,PDT::Spin1Half,
PDT::Spin1,PDT::Spin1Half);
// loop over the helicities
SpinorWaveFunction inter;
vector<double> me(3,0.);
vector<Complex> diag(2);
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
- inter = _gluonvertex->evaluate(mt,5,fin[ihel2].particle(),
+ inter = _gluonvertex->evaluate(mt,5,fin[ihel2].particle()->CC(),
fin[ihel2],gout[ohel1]);
diag[0] = _photonvertex->evaluate(0.*GeV2,inter,fout[ohel2],gmin[ihel1]);
// second diagram
- inter = _photonvertex->evaluate(0.*GeV2,5,fin[ihel2].particle(),
+ inter = _photonvertex->evaluate(0.*GeV2,5,fin[ihel2].particle()->CC(),
fin[ihel2],gmin[ihel1]);
diag[1] = _gluonvertex->evaluate(mt,inter,fout[ohel2],gout[ohel1]);
for(unsigned int ix=0;ix<2;++ix) me[ix] += norm(diag[ix]);
diag[0]+=diag[1];
me[2] += norm(diag[0]);
// matrix element if needed
if(calc) newme(2*ihel1,ihel2,2*ohel1,ohel2)=diag[0];
}
}
}
}
// save the info on the diagrams
if(!calc) {
DVector save;
save.push_back(me[0]);
save.push_back(me[1]);
meInfo(save);
}
// return the answer with colour and spin factors
if(calc) _me.reset(newme);
// // test vs the analytic result
// double test = -8./3.*norm(_photonvertex->getNorm())*norm(_gluonvertex->getNorm())
// *sqr(int(mePartonData()[1]->iCharge()))/9.
// *(uHat()/sHat()+sHat()/uHat());
// cerr << "testing the matrix element " << test << " " << me[2]/3.
// << " " << test/me[2]*3. << "\n";
return me[2]/3.;
}
double MEGammaP2Jets::gammaantiquarkME(vector<VectorWaveFunction> & gmin,
vector<SpinorBarWaveFunction> & fin,
vector<VectorWaveFunction> & gout,
vector<SpinorWaveFunction> & fout,
bool calc) const {
// the scale
Energy2 mt(scale());
// storage of the helicity me if needed
ProductionMatrixElement newme(PDT::Spin1,PDT::Spin1Half,
PDT::Spin1,PDT::Spin1Half);
// loop over the helicities
SpinorBarWaveFunction inter;
vector<double> me(3,0.);
vector<Complex> diag(2);
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
- inter = _gluonvertex->evaluate(mt,5,fin[ihel2].particle(),
+ inter = _gluonvertex->evaluate(mt,5,fin[ihel2].particle()->CC(),
fin[ihel2],gout[ohel1]);
diag[0] = _photonvertex->evaluate(0.*GeV2,fout[ohel2],inter,gmin[ihel1]);
// second diagram
- inter = _photonvertex->evaluate(0.*GeV2,5,fin[ihel2].particle(),
+ inter = _photonvertex->evaluate(0.*GeV2,5,fin[ihel2].particle()->CC(),
fin[ihel2],gmin[ihel1]);
diag[1] = _gluonvertex->evaluate(mt,fout[ohel2],inter,gout[ohel1]);
for(unsigned int ix=0;ix<2;++ix) me[ix] += norm(diag[ix]);
diag[0]+=diag[1];
me[2] += norm(diag[0]);
// matrix element if needed
if(calc) newme(2*ihel1,ihel2,2*ohel1,ohel2)=diag[0];
}
}
}
}
// save the info on the diagrams
if(!calc) {
DVector save;
save.push_back(me[0]);
save.push_back(me[1]);
meInfo(save);
}
// return the answer with colour and spin factors
if(calc) _me.reset(newme);
// test vs the analytic result
// double test = -8./3.*norm(_photonvertex->getNorm())*norm(_gluonvertex->getNorm())
// *sqr(int(mePartonData()[1]->iCharge()))/9.
// *(uHat()/sHat()+sHat()/uHat());
// cerr << "testing the matrix element " << test << " " << me[2]/3.
// << " " << test/me[2]*3. << "\n";
return me[2]/3.;
}
diff --git a/MatrixElement/Hadron/MEPP2QQ.cc b/MatrixElement/Hadron/MEPP2QQ.cc
--- a/MatrixElement/Hadron/MEPP2QQ.cc
+++ b/MatrixElement/Hadron/MEPP2QQ.cc
@@ -1,1051 +1,1051 @@
// -*- C++ -*-
//
// MEPP2QQ.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEPP2QQ class.
//
#include "MEPP2QQ.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig/MatrixElement/HardVertex.h"
#include "ThePEG/Repository/EventGenerator.h"
using namespace Herwig;
MEPP2QQ::MEPP2QQ() : _quarkflavour(6), _process(0),
_bottomopt(1), _topopt(1), _maxflavour(5) {
}
void MEPP2QQ::rebind(const TranslationMap & trans) {
_gggvertex = trans.translate( _gggvertex);
_qqgvertex = trans.translate( _qqgvertex);
_gluon = trans.translate( _gluon);
for(unsigned int ix=0;ix<_quark.size();++ix)
_quark[ix]=trans.translate(_quark[ix]);
for(unsigned int ix=0;ix<_antiquark.size();++ix)
_antiquark[ix]=trans.translate(_quark[ix]);
HwMEBase::rebind(trans);
}
IVector MEPP2QQ::getReferences() {
IVector ret = HwMEBase::getReferences();
ret.push_back(_gggvertex);
ret.push_back(_qqgvertex);
ret.push_back(_gluon);
for(unsigned int ix=0;ix<_quark.size();++ix)
ret.push_back(_quark[ix]);
for(unsigned int ix=0;ix<_antiquark.size();++ix)
ret.push_back(_antiquark[ix]);
return ret;
}
void MEPP2QQ::doinit() {
HwMEBase::doinit();
// handling of masses
if(_quarkflavour==6) {
massOption(vector<unsigned int>(2,_topopt));
}
else {
massOption(vector<unsigned int>(2,_bottomopt));
}
// get the vertex pointers from the SM object
tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
// do the initialisation
if(hwsm) {
_qqgvertex = hwsm->vertexFFG();
_gggvertex = hwsm->vertexGGG();
}
else throw InitException()
<< "Wrong type of StandardModel object in "
<< "MEPP2QQ::doinit() the Herwig version must be used"
<< Exception::runerror;
// get the particle data objects
_gluon=getParticleData(ParticleID::g);
for(int ix=1;ix<=6;++ix) {
_quark.push_back( getParticleData( ix));
_antiquark.push_back(getParticleData(-ix));
}
}
Energy2 MEPP2QQ::scale() const {
Energy2 mq2 = max(meMomenta()[2].mass2(),meMomenta()[3].mass2());
Energy2 s(0.5*sHat()),u(0.5*(uHat()-mq2)),t(0.5*(tHat()-mq2));
return 4.*s*t*u/(s*s+t*t+u*u);
}
void MEPP2QQ::persistentOutput(PersistentOStream & os) const {
os << _gggvertex << _qqgvertex << _quarkflavour << _bottomopt << _topopt
<< _process << _gluon << _quark << _antiquark << _maxflavour;
}
void MEPP2QQ::persistentInput(PersistentIStream & is, int) {
is >> _gggvertex >> _qqgvertex >> _quarkflavour >> _bottomopt >> _topopt
>> _process >> _gluon >> _quark >> _antiquark >> _maxflavour;
}
unsigned int MEPP2QQ::orderInAlphaS() const {
return 2;
}
unsigned int MEPP2QQ::orderInAlphaEW() const {
return 0;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEPP2QQ,HwMEBase>
describeHerwigMEPP2QQ("Herwig::MEPP2QQ", "HwMEHadron.so");
void MEPP2QQ::Init() {
static ClassDocumentation<MEPP2QQ> documentation
("The MEPP2QQ class implements the matrix element for"
" heavy quark pair production");
static Switch<MEPP2QQ,unsigned int> interfaceQuarkType
("QuarkType",
"The type of quark",
&MEPP2QQ::_quarkflavour, 6, false, false);
static SwitchOption interfaceQuarkTypeCharm
(interfaceQuarkType,
"Charm",
"Produce charm-anticharm",
4);
static SwitchOption interfaceQuarkTypeBottom
(interfaceQuarkType,
"Bottom",
"Produce bottom-antibottom",
5);
static SwitchOption interfaceQuarkTypeTop
(interfaceQuarkType,
"Top",
"Produce top-antitop",
6);
static Switch<MEPP2QQ,unsigned int> interfaceProcess
("Process",
"Which subprocesses to include",
&MEPP2QQ::_process, 0, false, false);
static SwitchOption interfaceProcessAll
(interfaceProcess,
"All",
"Include all subprocesses",
0);
static SwitchOption interfaceProcessPair
(interfaceProcess,
"Pair",
"Only include quark-antiquark pair production subprocesses",
1);
static SwitchOption interfaceProcess1
(interfaceProcess,
"gg2qqbar",
"Include only gg -> q qbar processes",
2);
static SwitchOption interfaceProcessqgqg
(interfaceProcess,
"qg2qg",
"Include only q g -> q g processes",
4);
static SwitchOption interfaceProcessqbargqbarg
(interfaceProcess,
"qbarg2qbarg",
"Include only qbar g -> qbar g processes",
5);
static SwitchOption interfaceProcessqqqq
(interfaceProcess,
"qq2qq",
"Include only q q -> q q processes",
6);
static SwitchOption interfaceProcessqbarqbarqbarqbar
(interfaceProcess,
"qbarqbar2qbarqbar",
"Include only qbar qbar -> qbar qbar processes",
7);
static SwitchOption interfaceProcessqqbarqqbar
(interfaceProcess,
"qqbar2qqbar",
"Include only q qbar -> q qbar processes",
8);
static Switch<MEPP2QQ,unsigned int> interfaceTopMassOption
("TopMassOption",
"Option for the treatment of the top quark mass",
&MEPP2QQ::_topopt, 1, false, false);
static SwitchOption interfaceTopMassOptionOnMassShell
(interfaceTopMassOption,
"OnMassShell",
"The top is produced on its mass shell",
1);
static SwitchOption interfaceTopMassOption2
(interfaceTopMassOption,
"OffShell",
"The top is generated off-shell using the mass and width generator.",
2);
static Switch<MEPP2QQ,unsigned int> interfaceBottomMassOption
("CharmBottomMassOption",
"Option for the treatment of bottom and lighter quark masses",
&MEPP2QQ::_bottomopt, 1, false, false);
static SwitchOption interfaceBottomMassOptionOnMassShell
(interfaceBottomMassOption,
"OnMassShell",
"The bottom is produced on its mass shell, this mass used everywhere",
1);
static SwitchOption interfaceBottomMassOptionZeroMass
(interfaceBottomMassOption,
"ZeroMass",
"The bottom is generated on mass shell, but zero mass is used in"
" the matrix element",
0);
static Parameter<MEPP2QQ,unsigned int> interfaceMaximumFlavour
("MaximumFlavour",
"The maximum flavour of the quarks in the process",
&MEPP2QQ::_maxflavour, 5, 1, 5,
false, false, Interface::limited);
}
Selector<MEBase::DiagramIndex>
MEPP2QQ::diagrams(const DiagramVector & diags) const {
// select the diagram, this is easy for us as we have already done it
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
if(diags[i]->id()==-int(_diagram)) sel.insert(1.0, i);
else sel.insert(0., i);
}
return sel;
}
double MEPP2QQ::gg2qqbarME(vector<VectorWaveFunction> &g1,
vector<VectorWaveFunction> &g2,
vector<SpinorBarWaveFunction> & q,
vector<SpinorWaveFunction> & qbar,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
Energy mass = q[0].mass();
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1Half));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorWaveFunction inters;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
interv=_gggvertex->evaluate(mt,5,_gluon,g1[ihel1],g2[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
//first t-channel diagram
- inters =_qqgvertex->evaluate(mt,3,qbar[ohel2].particle(),
+ inters =_qqgvertex->evaluate(mt,3,qbar[ohel2].particle()->CC(),
qbar[ohel2],g2[ihel2],mass);
diag[0]=_qqgvertex->evaluate(mt,inters,q[ohel1],g1[ihel1]);
//second t-channel diagram
- inters =_qqgvertex->evaluate(mt,3,qbar[ohel2].particle(),
+ inters =_qqgvertex->evaluate(mt,3,qbar[ohel2].particle()->CC(),
qbar[ohel2],g1[ihel1],mass);
diag[1]=_qqgvertex->evaluate(mt,inters,q[ohel1],g2[ihel2]);
// s-channel diagram
diag[2]=_qqgvertex->evaluate(mt,qbar[ohel2],q[ohel1],interv);
// colour flows
flow[0]=diag[0]+diag[2];
flow[1]=diag[1]-diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(2*ihel1,2*ihel2,ohel1,ohel2)=flow[iflow-1];
}
}
}
}
// select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=4+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
// Energy2 mq2=sqr(getParticleData(6)->mass());
// double tau1=(tHat()-mq2)/sHat(),tau2=(uHat()-mq2)/sHat(),rho=4*mq2/sHat();
// double test=(1./6./tau1/tau2-3./8.)*sqr(4.*pi*SM().alphaS(mt))*
// (sqr(tau1)+sqr(tau2)+rho-0.25*sqr(rho)/tau1/tau2);
// cerr << "testing matrix element "
// << output/48./test << "\n";
return output/48.;
}
double MEPP2QQ::qg2qgME(vector<SpinorWaveFunction> & qin,
vector<VectorWaveFunction> &g2,
vector<SpinorBarWaveFunction> & qout,
vector<VectorWaveFunction> &g4,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
Energy mass = qout[0].mass();
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorWaveFunction inters,inters2;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
inters=_qqgvertex->evaluate(mt,3,qin[ihel1].particle()->CC(),
qin[ihel1],g2[ihel2],mass);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// s-channel diagram
diag[0]=_qqgvertex->evaluate(mt,inters,qout[ohel1],g4[ohel2]);
// first t-channel
inters2=_qqgvertex->evaluate(mt,3,qin[ihel1].particle()->CC(),
qin[ihel1],g4[ohel2],mass);
diag[1]=_qqgvertex->evaluate(mt,inters2,qout[ohel1],g2[ihel2]);
// second t-channel
interv=_qqgvertex->evaluate(mt,5,_gluon,qin[ihel1],qout[ohel1]);
diag[2]=_gggvertex->evaluate(mt,g2[ihel2],g4[ohel2],interv);
// colour flows
flow[0]=diag[0]-diag[2];
flow[1]=diag[1]+diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,2*ihel2,ohel1,2*ohel2)=flow[iflow-1];
}
}
}
}
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//cerr << "testing matrix element "
// << 18./output*(-4./9./s/u+1./t/t)*(s*s+u*u)*sqr(alphas) << endl;
//select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=10+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return output/18.;
}
double MEPP2QQ::qbarg2qbargME(vector<SpinorBarWaveFunction> & qin,
vector<VectorWaveFunction> &g2,
vector<SpinorWaveFunction> & qout,
vector<VectorWaveFunction> &g4,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
Energy mass = qout[0].mass();
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorBarWaveFunction inters,inters2;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
inters=_qqgvertex->evaluate(mt,3,qin[ihel1].particle()->CC(),
qin[ihel1],g2[ihel2],mass);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// s-channel diagram
diag[0]=_qqgvertex->evaluate(mt,qout[ohel1],inters,g4[ohel2]);
// first t-channel
inters2=_qqgvertex->evaluate(mt,3,qin[ihel1].particle()->CC(),
qin[ihel1],g4[ohel2],mass);
diag[1]=_qqgvertex->evaluate(mt,qout[ohel1],inters2,g2[ihel2]);
// second t-channel
interv=_qqgvertex->evaluate(mt,5,_gluon,qout[ohel1],qin[ihel1]);
diag[2]=_gggvertex->evaluate(mt,g2[ihel2],g4[ohel2],interv);
// colour flows
flow[0]=diag[0]+diag[2];
flow[1]=diag[1]-diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,2*ihel2,ohel1,2*ohel2)=flow[iflow-1];
}
}
}
}
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//cerr << "testing matrix element "
// << 18./output*(-4./9./s/u+1./t/t)*(s*s+u*u)*sqr(alphas) << endl;
//select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=13+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return output/18.;
}
double MEPP2QQ::qq2qqME(vector<SpinorWaveFunction> & q1,
vector<SpinorWaveFunction> & q2,
vector<SpinorBarWaveFunction> & q3,
vector<SpinorBarWaveFunction> & q4,
unsigned int iflow) const {
// identify special case of identical quarks
bool identical(q1[0].id()==q2[0].id());
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));
// calculate the matrix element
double output(0.),sumdiag[2]={0.,0.};
Complex diag[2];
VectorWaveFunction interv;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q3[ohel1]);
diag[0] = _qqgvertex->evaluate(mt,q2[ihel2],q4[ohel2],interv);
// second diagram if identical
if(identical) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q4[ohel2]);
diag[1]=_qqgvertex->evaluate(mt,q2[ihel2],q3[ohel1],interv);
}
else diag[1]=0.;
// sum of diagrams
for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]);
// total
output +=real(diag[0]*conj(diag[0])+diag[1]*conj(diag[1])
+2./3.*diag[0]*conj(diag[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,ihel2,ohel1,ohel2)=diag[iflow-1];
}
}
}
}
// identical particle symmetry factor if needed
if(identical) output*=0.5;
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//if(identical)
// {cerr << "testing matrix element A "
// << 18./output*0.5*(4./9.*((s*s+u*u)/t/t+(s*s+t*t)/u/u)
// -8./27.*s*s/u/t)*sqr(alphas) << endl;}
//else
// {cerr << "testing matrix element B "
// << 18./output*(4./9.*(s*s+u*u)/t/t)*sqr(alphas) << endl;}
//select a colour flow
_flow=1+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=16+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// final part of colour and spin factors
return output/18.;
}
double MEPP2QQ::qbarqbar2qbarqbarME(vector<SpinorBarWaveFunction> & q1,
vector<SpinorBarWaveFunction> & q2,
vector<SpinorWaveFunction> & q3,
vector<SpinorWaveFunction> & q4,
unsigned int iflow) const {
// identify special case of identical quarks
bool identical(q1[0].id()==q2[0].id());
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0)
{_me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));}
// calculate the matrix element
double output(0.),sumdiag[2]={0.,0.};
Complex diag[2];
VectorWaveFunction interv;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
interv = _qqgvertex->evaluate(mt,5,_gluon,q3[ohel1],q1[ihel1]);
diag[0] = _qqgvertex->evaluate(mt,q4[ohel2],q2[ihel2],interv);
// second diagram if identical
if(identical) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q4[ohel2],q1[ihel1]);
diag[1]=_qqgvertex->evaluate(mt,q3[ohel1],q2[ihel2],interv);
}
else diag[1]=0.;
// sum of diagrams
for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]);
// total
output +=real(diag[0]*conj(diag[0])+diag[1]*conj(diag[1])
+2./3.*diag[0]*conj(diag[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,ihel2,ohel1,ohel2)=diag[iflow-1];
}
}
}
}
// identical particle symmetry factor if needed
if(identical){output*=0.5;}
// test code vs me from ESW
// Energy2 u(uHat()),t(tHat()),s(sHat());
// double alphas(4.*pi*SM().alphaS(mt));
// if(identical)
// {cerr << "testing matrix element A "
// << 18./output*0.5*(4./9.*((s*s+u*u)/t/t+(s*s+t*t)/u/u)
// -8./27.*s*s/u/t)*sqr(alphas) << endl;}
// else
// {cerr << "testing matrix element B "
// << 18./output*(4./9.*(s*s+u*u)/t/t)*sqr(alphas) << endl;}
//select a colour flow
_flow=1+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=18+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// final part of colour and spin factors
return output/18.;
}
double MEPP2QQ::qqbar2qqbarME(vector<SpinorWaveFunction> & q1,
vector<SpinorBarWaveFunction> & q2,
vector<SpinorBarWaveFunction> & q3,
vector<SpinorWaveFunction> & q4,
unsigned int iflow) const {
// type of process
bool diagon[2]={q1[0].id()== -q2[0].id(),
q1[0].id()== -q3[0].id()};
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));
// calculate the matrix element
double output(0.),sumdiag[2]={0.,0.};
Complex diag[2];
VectorWaveFunction interv;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
if(diagon[0]) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q2[ihel2]);
diag[0] = _qqgvertex->evaluate(mt,q4[ohel2],q3[ohel1],interv);
}
else diag[0]=0.;
// second diagram
if(diagon[1]) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q3[ohel1]);
diag[1]=_qqgvertex->evaluate(mt,q4[ohel2],q2[ihel2],interv);
}
else diag[1]=0.;
// sum of diagrams
for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]);
// total
output +=real(diag[0]*conj(diag[0])+diag[1]*conj(diag[1])
+2./3.*diag[0]*conj(diag[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,ihel2,ohel1,ohel2)=diag[iflow-1];
}
}
}
}
//select a colour flow
_flow=1+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=20+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// final part of colour and spin factors
// Energy2 mq2=sqr(getParticleData(6)->mass());
// double tau1=(tHat()-mq2)/sHat(),tau2=(uHat()-mq2)/sHat(),rho=4*mq2/sHat();
// cerr << "testing matrix element "
// << output/18./(4./9.*sqr(4.*pi*SM().alphaS(mt))*(sqr(tau1)+sqr(tau2)+0.5*rho))
// << "\n";
return output/18.;
}
void MEPP2QQ::getDiagrams() const {
// gg -> q qbar subprocesses
if(_process==0||_process==1||_process==2) {
// first t-channel
add(new_ptr((Tree2toNDiagram(3),_gluon,_antiquark[_quarkflavour-1],_gluon,
1,_quark[_quarkflavour-1], 2,_antiquark[_quarkflavour-1],-4)));
// interchange
add(new_ptr((Tree2toNDiagram(3),_gluon,_antiquark[_quarkflavour-1],_gluon,
2,_quark[_quarkflavour-1], 1,_antiquark[_quarkflavour-1],-5)));
// s-channel
add(new_ptr((Tree2toNDiagram(2),_gluon,_gluon, 1, _gluon,
3,_quark[_quarkflavour-1], 3, _antiquark[_quarkflavour-1], -6)));
}
// q g -> q g subprocesses
if(_process==0||_process==4) {
// s-channel
add(new_ptr((Tree2toNDiagram(2),_quark[_quarkflavour-1], _gluon,
1, _quark[_quarkflavour-1],
3, _quark[_quarkflavour-1], 3, _gluon,-10)));
// quark t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[_quarkflavour-1],
_quark[_quarkflavour-1],_gluon,
2,_quark[_quarkflavour-1],1,_gluon,-11)));
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[_quarkflavour-1],_gluon,_gluon,
1,_quark[_quarkflavour-1],2,_gluon,-12)));
}
// qbar g -> qbar g subprocesses
if(_process==0||_process==5) {
// s-channel
add(new_ptr((Tree2toNDiagram(2),_antiquark[_quarkflavour-1], _gluon,
1, _antiquark[_quarkflavour-1],
3, _antiquark[_quarkflavour-1], 3, _gluon,-13)));
// quark t-channel
add(new_ptr((Tree2toNDiagram(3),_antiquark[_quarkflavour-1],
_antiquark[_quarkflavour-1],_gluon,
2,_antiquark[_quarkflavour-1],1,_gluon,-14)));
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_antiquark[_quarkflavour-1],_gluon,_gluon,
1,_antiquark[_quarkflavour-1],2,_gluon,-15)));
}
// processes involving two quark lines
for(unsigned int iy=0;iy<_maxflavour;++iy) {
// q q -> q q subprocesses
if(_process==0||_process==6) {
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[_quarkflavour-1],_gluon,_quark[iy],
1,_quark[_quarkflavour-1],2,_quark[iy],-16)));
// exchange for identical quarks
if(_quarkflavour-1==iy)
add(new_ptr((Tree2toNDiagram(3),_quark[_quarkflavour-1],_gluon,_quark[iy],
2,_quark[_quarkflavour-1],1,_quark[iy],-17)));
}
// qbar qbar -> qbar qbar subprocesses
if(_process==0||_process==7) {
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_antiquark[_quarkflavour-1],_gluon,_antiquark[iy],
1,_antiquark[_quarkflavour-1],2,_antiquark[iy],-18)));
// exchange for identical quarks
if(_quarkflavour-1==iy)
add(new_ptr((Tree2toNDiagram(3),_antiquark[_quarkflavour-1],_gluon,_antiquark[iy],
2,_antiquark[_quarkflavour-1],1,_antiquark[iy],-19)));
}
// q qbar -> q qbar
if(_process==0||_process==1||_process==8) {
// gluon s-channel
add(new_ptr((Tree2toNDiagram(2),_quark[iy], _antiquark[iy],
1, _gluon, 3, _quark[_quarkflavour-1],
3, _antiquark[_quarkflavour-1],-20)));
if(_process!=1) {
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[iy],_gluon,_antiquark[_quarkflavour-1],
1,_quark[iy],2,_antiquark[_quarkflavour-1],-21)));
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[_quarkflavour-1],_gluon,_antiquark[iy],
1,_quark[_quarkflavour-1],2,_antiquark[iy],-21)));
}
}
}
}
Selector<const ColourLines *>
MEPP2QQ::colourGeometries(tcDiagPtr diag) const {
// colour lines for gg to q qbar
static const ColourLines cggqq[4]={ColourLines("1 4, -1 -2 3, -3 -5"),
ColourLines("3 4, -3 -2 1, -1 -5"),
ColourLines("2 -1, 1 3 4, -2 -3 -5"),
ColourLines("1 -2, -1 -3 -5, 2 3 4")};
// colour lines for q g to q g
static const ColourLines cqgqg[4]={ColourLines("1 -2, 2 3 5, 4 -5"),
ColourLines("1 5, 3 4,-3 2 -5 "),
ColourLines("1 2 -3, 3 5, -5 -2 4"),
ColourLines("1 -2 5,3 2 4,-3 -5")};
// colour lines for qbar g -> qbar g
static const ColourLines cqbgqbg[4]={ColourLines("-1 2, -2 -3 -5, -4 5"),
ColourLines("-1 -5, -3 -4, 3 -2 5"),
ColourLines("-1 -2 3, -3 -5, 5 2 -4"),
ColourLines("-1 2 -5,-3 -2 -4, 3 5")};
// colour lines for q q -> q q
static const ColourLines cqqqq[2]={ColourLines("1 2 5,3 -2 4"),
ColourLines("1 2 4,3 -2 5")};
// colour lines for qbar qbar -> qbar qbar
static const ColourLines cqbqbqbqb[2]={ColourLines("-1 -2 -5,-3 2 -4"),
ColourLines("-1 -2 -4,-3 2 -5")};
// colour lines for q qbar -> q qbar
static const ColourLines cqqbqqb[2]={ColourLines("1 3 4,-2 -3 -5"),
ColourLines("1 2 -3,4 -2 -5")};
// select the colour flow (as all ready picked just insert answer)
Selector<const ColourLines *> sel;
switch(abs(diag->id())) {
// gg -> q qbar subprocess
case 4: case 5:
sel.insert(1.0, &cggqq[abs(diag->id())-4]);
break;
case 6:
sel.insert(1.0, &cggqq[1+_flow]);
break;
// q g -> q g subprocess
case 10: case 11:
sel.insert(1.0, &cqgqg[abs(diag->id())-10]);
break;
case 12:
sel.insert(1.0, &cqgqg[1+_flow]);
break;
// q g -> q g subprocess
case 13: case 14:
sel.insert(1.0, &cqbgqbg[abs(diag->id())-13]);
break;
case 15:
sel.insert(1.0, &cqbgqbg[1+_flow]);
break;
// q q -> q q subprocess
case 16: case 17:
sel.insert(1.0, &cqqqq[abs(diag->id())-16]);
break;
// qbar qbar -> qbar qbar subprocess
case 18: case 19:
sel.insert(1.0, &cqbqbqbqb[abs(diag->id())-18]);
break;
// q qbar -> q qbar subprocess
case 20: case 21:
sel.insert(1.0, &cqqbqqb[abs(diag->id())-20]);
break;
}
return sel;
}
double MEPP2QQ::me2() const {
// total matrix element
double me(0.);
// gg initiated processes
if(mePartonData()[0]->id()==ParticleID::g&&mePartonData()[1]->id()==ParticleID::g) {
VectorWaveFunction g1w(rescaledMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(rescaledMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction qw(rescaledMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction qbarw(rescaledMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g1,g2;
vector<SpinorBarWaveFunction> q;
vector<SpinorWaveFunction> qbar;
for(unsigned int ix=0;ix<2;++ix) {
g1w.reset(2*ix);g1.push_back(g1w);
g2w.reset(2*ix);g2.push_back(g2w);
qw.reset(ix);q.push_back(qw);
qbarw.reset(ix);qbar.push_back(qbarw);
}
// calculate the matrix element
me=gg2qqbarME(g1,g2,q,qbar,0);
}
// quark first processes
else if(mePartonData()[0]->id()>0) {
// q g -> q g
if(mePartonData()[1]->id()==ParticleID::g) {
SpinorWaveFunction qinw(rescaledMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(rescaledMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction qoutw(rescaledMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction g4w(rescaledMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g2,g4;
vector<SpinorWaveFunction> qin;
vector<SpinorBarWaveFunction> qout;
for(unsigned int ix=0;ix<2;++ix) {
qinw.reset(ix);qin.push_back(qinw);
g2w.reset(2*ix);g2.push_back(g2w);
qoutw.reset(ix);qout.push_back(qoutw);
g4w.reset(2*ix);g4.push_back(g4w);
}
// calculate the matrix element
me = qg2qgME(qin,g2,qout,g4,0);
}
// q qbar to q qbar
else if(mePartonData()[1]->id()<0) {
SpinorWaveFunction q1w(rescaledMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction q2w(rescaledMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction q3w(rescaledMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction q4w(rescaledMomenta()[3],mePartonData()[3],outgoing);
vector<SpinorWaveFunction> q1,q4;
vector<SpinorBarWaveFunction> q2,q3;
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1.push_back(q1w);
q2w.reset(ix);q2.push_back(q2w);
q3w.reset(ix);q3.push_back(q3w);
q4w.reset(ix);q4.push_back(q4w);
}
// calculate the matrix element
me = qqbar2qqbarME(q1,q2,q3,q4,0);
}
// q q -> q q
else if(mePartonData()[1]->id()>0) {
SpinorWaveFunction q1w(rescaledMomenta()[0],mePartonData()[0],incoming);
SpinorWaveFunction q2w(rescaledMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction q3w(rescaledMomenta()[2],mePartonData()[2],outgoing);
SpinorBarWaveFunction q4w(rescaledMomenta()[3],mePartonData()[3],outgoing);
vector<SpinorWaveFunction> q1,q2;
vector<SpinorBarWaveFunction> q3,q4;
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1.push_back(q1w);
q2w.reset(ix);q2.push_back(q2w);
q3w.reset(ix);q3.push_back(q3w);
q4w.reset(ix);q4.push_back(q4w);
}
// calculate the matrix element
me = qq2qqME(q1,q2,q3,q4,0);
}
}
// antiquark first processes
else if(mePartonData()[0]->id()<0) {
// qbar g -> qbar g
if(mePartonData()[1]->id()==ParticleID::g) {
SpinorBarWaveFunction qinw(rescaledMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(rescaledMomenta()[1],mePartonData()[1],incoming);
SpinorWaveFunction qoutw(rescaledMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction g4w(rescaledMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g2,g4;
vector<SpinorBarWaveFunction> qin;
vector<SpinorWaveFunction> qout;
for(unsigned int ix=0;ix<2;++ix) {
qinw.reset(ix);qin.push_back(qinw);
g2w.reset(2*ix);g2.push_back(g2w);
qoutw.reset(ix);qout.push_back(qoutw);
g4w.reset(2*ix);g4.push_back(g4w);
}
// calculate the matrix element
me = qbarg2qbargME(qin,g2,qout,g4,0);
}
// qbar qbar -> qbar qbar
else if(mePartonData()[1]->id()<0) {
SpinorBarWaveFunction q1w(rescaledMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction q2w(rescaledMomenta()[1],mePartonData()[1],incoming);
SpinorWaveFunction q3w(rescaledMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction q4w(rescaledMomenta()[3],mePartonData()[3],outgoing);
vector<SpinorBarWaveFunction> q1,q2;
vector<SpinorWaveFunction> q3,q4;
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1.push_back(q1w);
q2w.reset(ix);q2.push_back(q2w);
q3w.reset(ix);q3.push_back(q3w);
q4w.reset(ix);q4.push_back(q4w);
}
// calculate the matrix element
me = qbarqbar2qbarqbarME(q1,q2,q3,q4,0);
}
}
else throw Exception() << "Unknown process in MEPP2QQ::me2()"
<< Exception::abortnow;
// return the answer
return me;
}
void MEPP2QQ::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard;
hard.push_back(sub->incoming().first);
hard.push_back(sub->incoming().second);
hard.push_back(sub->outgoing()[0]);
hard.push_back(sub->outgoing()[1]);
// identify the process and calculate the matrix element
if(hard[0]->id()==ParticleID::g&&hard[1]->id()==ParticleID::g) {
if(hard[2]->id()<0) swap(hard[2],hard[3]);
vector<VectorWaveFunction> g1,g2;
vector<SpinorBarWaveFunction> q;
vector<SpinorWaveFunction> qbar;
// off-shell wavefunctions for the spin correlations
VectorWaveFunction( g1,hard[0],incoming,false,true,true);
VectorWaveFunction( g2,hard[1],incoming,false,true,true);
SpinorBarWaveFunction(q ,hard[2],outgoing,true ,true);
SpinorWaveFunction( qbar,hard[3],outgoing,true ,true);
g1[1]=g1[2];g2[1]=g2[2];
// on-shell for matrix element
vector<Lorentz5Momentum> momenta;
cPDVector data;
for(unsigned int ix=0;ix<4;++ix) {
momenta.push_back(hard[ix]->momentum());
data .push_back(hard[ix]->dataPtr());
}
rescaleMomenta(momenta,data);
VectorWaveFunction g1w(rescaledMomenta()[0],data[0],incoming);
VectorWaveFunction g2w(rescaledMomenta()[1],data[1],incoming);
SpinorBarWaveFunction qw(rescaledMomenta()[2],data[2],outgoing);
SpinorWaveFunction qbarw(rescaledMomenta()[3],data[3],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
g1w.reset(2*ix); g1 [ix] = g1w ;
g2w.reset(2*ix); g2 [ix] = g2w ;
qw.reset(ix); q [ix] = qw ;
qbarw.reset(ix); qbar[ix] = qbarw;
}
gg2qqbarME(g1,g2,q,qbar,_flow);
}
else if(hard[0]->id()==ParticleID::g||hard[1]->id()==ParticleID::g) {
if(hard[0]->id()==ParticleID::g) swap(hard[0],hard[1]);
if(hard[2]->id()==ParticleID::g) swap(hard[2],hard[3]);
// on-shell for matrix element
vector<Lorentz5Momentum> momenta;
cPDVector data;
for(unsigned int ix=0;ix<4;++ix) {
momenta.push_back(hard[ix]->momentum());
data .push_back(hard[ix]->dataPtr());
}
rescaleMomenta(momenta,data);
vector<VectorWaveFunction> g2,g4;
VectorWaveFunction(g2,hard[1],incoming,false,true,true);
VectorWaveFunction(g4,hard[3],outgoing,true ,true,true);
VectorWaveFunction gin (rescaledMomenta()[1],data[1],incoming);
VectorWaveFunction gout(rescaledMomenta()[3],data[3],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
gin .reset(2*ix); g2 [ix] = gin ;
gout.reset(2*ix); g4 [ix] = gout;
}
// q g -> q g
if(hard[0]->id()>0) {
vector<SpinorWaveFunction> q1;
vector<SpinorBarWaveFunction> q3;
SpinorWaveFunction( q1,hard[0],incoming,false,true);
SpinorBarWaveFunction(q3,hard[2],outgoing,true ,true);
SpinorWaveFunction qin (rescaledMomenta()[0],data[0],incoming);
SpinorBarWaveFunction qout(rescaledMomenta()[2],data[2],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
qin .reset(ix); q1[ix] = qin ;
qout.reset(ix); q3[ix] = qout;
}
qg2qgME(q1,g2,q3,g4,_flow);
}
// qbar g -> qbar g
else {
vector<SpinorBarWaveFunction> q1;
vector<SpinorWaveFunction> q3;
SpinorBarWaveFunction(q1,hard[0],incoming,false,true);
SpinorWaveFunction( q3,hard[2],outgoing,true ,true);
SpinorBarWaveFunction qin (rescaledMomenta()[0],data[0],incoming);
SpinorWaveFunction qout(rescaledMomenta()[2],data[2],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
qin .reset(ix); q1[ix] = qin ;
qout.reset(ix); q3[ix] = qout;
}
qbarg2qbargME(q1,g2,q3,g4,_flow);
}
}
else if(hard[0]->id()>0||hard[1]->id()>0) {
// q q -> q q
if(hard[0]->id()>0&&hard[1]->id()>0) {
if(hard[2]->id()!=hard[0]->id()) swap(hard[2],hard[3]);
vector<SpinorWaveFunction> q1,q2;
vector<SpinorBarWaveFunction> q3,q4;
SpinorWaveFunction( q1,hard[0],incoming,false,true);
SpinorWaveFunction( q2,hard[1],incoming,false,true);
SpinorBarWaveFunction(q3,hard[2],outgoing,true ,true);
SpinorBarWaveFunction(q4,hard[3],outgoing,true ,true);
// on-shell for matrix element
vector<Lorentz5Momentum> momenta;
cPDVector data;
for(unsigned int ix=0;ix<4;++ix) {
momenta.push_back(hard[ix]->momentum());
data .push_back(hard[ix]->dataPtr());
}
rescaleMomenta(momenta,data);
SpinorWaveFunction q1w(rescaledMomenta()[0],data[0],incoming);
SpinorWaveFunction q2w(rescaledMomenta()[1],data[1],incoming);
SpinorBarWaveFunction q3w(rescaledMomenta()[2],data[2],outgoing);
SpinorBarWaveFunction q4w(rescaledMomenta()[3],data[3],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1[ix] = q1w;
q2w.reset(ix);q2[ix] = q2w;
q3w.reset(ix);q3[ix] = q3w;
q4w.reset(ix);q4[ix] = q4w;
}
qq2qqME(q1,q2,q3,q4,_flow);
}
// q qbar -> q qbar
else {
if(hard[0]->id()<0) swap(hard[0],hard[1]);
if(hard[2]->id()<0) swap(hard[2],hard[3]);
vector<SpinorWaveFunction> q1,q4;
vector<SpinorBarWaveFunction> q2,q3;
// off-shell for spin correlations
SpinorWaveFunction( q1,hard[0],incoming,false,true);
SpinorBarWaveFunction(q2,hard[1],incoming,false,true);
SpinorBarWaveFunction(q3,hard[2],outgoing,true ,true);
SpinorWaveFunction( q4,hard[3],outgoing,true ,true);
// on-shell for matrix element
vector<Lorentz5Momentum> momenta;
cPDVector data;
for(unsigned int ix=0;ix<4;++ix) {
momenta.push_back(hard[ix]->momentum());
data .push_back(hard[ix]->dataPtr());
}
rescaleMomenta(momenta,data);
SpinorWaveFunction q1w(rescaledMomenta()[0],data[0],incoming);
SpinorBarWaveFunction q2w(rescaledMomenta()[1],data[1],incoming);
SpinorBarWaveFunction q3w(rescaledMomenta()[2],data[2],outgoing);
SpinorWaveFunction q4w(rescaledMomenta()[3],data[3],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1[ix] = q1w;
q2w.reset(ix);q2[ix] = q2w;
q3w.reset(ix);q3[ix] = q3w;
q4w.reset(ix);q4[ix] = q4w;
}
qqbar2qqbarME(q1,q2,q3,q4,_flow);
}
}
else if (hard[0]->id()<0&&hard[1]->id()<0) {
if(hard[2]->id()!=hard[0]->id()) swap(hard[2],hard[3]);
vector<SpinorBarWaveFunction> q1,q2;
vector<SpinorWaveFunction> q3,q4;
SpinorBarWaveFunction(q1,hard[0],incoming,false,true);
SpinorBarWaveFunction(q2,hard[1],incoming,false,true);
SpinorWaveFunction( q3,hard[2],outgoing,true ,true);
SpinorWaveFunction( q4,hard[3],outgoing,true ,true);
// on-shell for matrix element
vector<Lorentz5Momentum> momenta;
cPDVector data;
for(unsigned int ix=0;ix<4;++ix) {
momenta.push_back(hard[ix]->momentum());
data .push_back(hard[ix]->dataPtr());
}
rescaleMomenta(momenta,data);
SpinorBarWaveFunction q1w(rescaledMomenta()[0],data[0],incoming);
SpinorBarWaveFunction q2w(rescaledMomenta()[1],data[1],incoming);
SpinorWaveFunction q3w(rescaledMomenta()[2],data[2],outgoing);
SpinorWaveFunction q4w(rescaledMomenta()[3],data[3],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1[ix] = q1w;
q2w.reset(ix);q2[ix] = q2w;
q3w.reset(ix);q3[ix] = q3w;
q4w.reset(ix);q4[ix] = q4w;
}
qbarqbar2qbarqbarME(q1,q2,q3,q4,_flow);
}
else throw Exception() << "Unknown process in MEPP2QQ::constructVertex()"
<< Exception::runerror;
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(_me);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<4;++ix)
hard[ix]->spinInfo()->productionVertex(hardvertex);
}
diff --git a/MatrixElement/Hadron/MEPP2QQHiggs.cc b/MatrixElement/Hadron/MEPP2QQHiggs.cc
--- a/MatrixElement/Hadron/MEPP2QQHiggs.cc
+++ b/MatrixElement/Hadron/MEPP2QQHiggs.cc
@@ -1,636 +1,636 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEPP2QQH class.
//
#include "MEPP2QQHiggs.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/Cuts/Cuts.h"
#include "Herwig/MatrixElement/HardVertex.h"
using namespace Herwig;
MEPP2QQHiggs::MEPP2QQHiggs() : quarkFlavour_(6), process_(0), shapeOpt_(2),
mh_(), wh_(), alpha_(1.1)
{}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEPP2QQHiggs,HwMEBase>
describeHerwigMEPP2QQHiggs("Herwig::MEPP2QQHiggs", "HwMEHadron.so");
void MEPP2QQHiggs::Init() {
static ClassDocumentation<MEPP2QQHiggs> documentation
("The MEPP2QQHiggs class implements the matrix elements for the "
"production of the Higgs boson in association with a heavy quark-antiquark pair");
static Switch<MEPP2QQHiggs,unsigned int> interfaceQuarkType
("QuarkType",
"The type of quark",
&MEPP2QQHiggs::quarkFlavour_, 6, false, false);
static SwitchOption interfaceQuarkTypeBottom
(interfaceQuarkType,
"Bottom",
"Produce bottom-antibottom",
5);
static SwitchOption interfaceQuarkTypeTop
(interfaceQuarkType,
"Top",
"Produce top-antitop",
6);
static Switch<MEPP2QQHiggs,unsigned int> interfaceProcess
("Process",
"Which subprocesses to include",
&MEPP2QQHiggs::process_, 0, false, false);
static SwitchOption interfaceProcessAll
(interfaceProcess,
"All",
"Include all subprocesses",
0);
static SwitchOption interfaceProcess1
(interfaceProcess,
"gg",
"Include only gg -> QQbarH processes",
1);
static SwitchOption interfaceProcessqbarqbarqbarqbar
(interfaceProcess,
"qqbar",
"Include only qbar qbar -> QQbarH processes",
2);
static Switch<MEPP2QQHiggs,unsigned int> interfaceShapeOption
("ShapeScheme",
"Option for the treatment of the Higgs resonance shape",
&MEPP2QQHiggs::shapeOpt_, 2, false, false);
static SwitchOption interfaceStandardShapeFixed
(interfaceShapeOption,
"FixedBreitWigner",
"Breit-Wigner s-channel resonance",
1);
static SwitchOption interfaceStandardShapeRunning
(interfaceShapeOption,
"MassGenerator",
"Use the mass generator to give the shape",
2);
static SwitchOption interfaceStandardShapeOnShell
(interfaceShapeOption,
"OnShell",
"Produce an on-shell Higgs boson",
0);
static Parameter<MEPP2QQHiggs,double> interfaceAlpha
("Alpha",
"Power for the generation of the tranverse mass in the pT mapping",
&MEPP2QQHiggs::alpha_, 1.1, 0.0, 10.0,
false, false, Interface::limited);
}
Energy2 MEPP2QQHiggs::scale() const {
return sHat();
// return sqr(mePartonData()[2]->mass()+mePartonData()[3]->mass()+
// mePartonData()[4]->mass());
}
int MEPP2QQHiggs::nDim() const {
return 4 + int(shapeOpt_>0);
}
unsigned int MEPP2QQHiggs::orderInAlphaS() const {
return 2;
}
unsigned int MEPP2QQHiggs::orderInAlphaEW() const {
return 1;
}
IBPtr MEPP2QQHiggs::clone() const {
return new_ptr(*this);
}
IBPtr MEPP2QQHiggs::fullclone() const {
return new_ptr(*this);
}
void MEPP2QQHiggs::setKinematics() {
HwMEBase::setKinematics();
}
void MEPP2QQHiggs::persistentOutput(PersistentOStream & os) const {
os << quarkFlavour_ << process_ << shapeOpt_
<< ounit(mh_,GeV) << ounit(wh_,GeV) << hmass_
<< GGGVertex_ << QQGVertex_ << QQHVertex_
<< gluon_ << higgs_ << quark_ << antiquark_
<< alpha_;
}
void MEPP2QQHiggs::persistentInput(PersistentIStream & is, int) {
is >> quarkFlavour_ >> process_ >> shapeOpt_
>> iunit(mh_,GeV) >> iunit(wh_,GeV) >> hmass_
>> GGGVertex_ >> QQGVertex_ >> QQHVertex_
>> gluon_ >> higgs_ >> quark_ >> antiquark_
>> alpha_;
}
void MEPP2QQHiggs::doinit() {
HwMEBase::doinit();
// stuff for the higgs mass
higgs_=getParticleData(ParticleID::h0);
mh_ = higgs_->mass();
wh_ = higgs_->width();
if(higgs_->massGenerator()) {
hmass_=dynamic_ptr_cast<GenericMassGeneratorPtr>(higgs_->massGenerator());
}
if(shapeOpt_==2&&!hmass_)
throw InitException()
<< "If using the mass generator for the line shape in MEPP2QQHiggs::doinit()"
<< "the mass generator must be an instance of the GenericMassGenerator class"
<< Exception::runerror;
// get the vertex pointers from the SM object
tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
// do the initialisation
if(!hwsm) throw InitException() << "Wrong type of StandardModel object in "
<< "MEPP2QQHiggs::doinit() the Herwig"
<< " version must be used"
<< Exception::runerror;
GGGVertex_ = hwsm->vertexGGG();
QQGVertex_ = hwsm->vertexFFG();
QQHVertex_ = hwsm->vertexFFH();
// get the particle data objects
gluon_=getParticleData(ParticleID::g);
for(int ix=1;ix<=6;++ix) {
quark_.push_back( getParticleData( ix));
antiquark_.push_back(getParticleData(-ix));
}
}
bool MEPP2QQHiggs::generateKinematics(const double * r) {
jacobian(1.);
// CMS energy
Energy rs = sqrt(sHat());
// quark mass
Energy mq(quark_[quarkFlavour_-1]->mass());
// generate the higgs mass
Energy mh(mh_);
if(shapeOpt_!=0) {
Energy mhmax = min(rs-2.*mq,higgs_->massMax());
Energy mhmin = max(ZERO ,higgs_->massMin());
if(mhmax<=mhmin) return false;
double rhomin = atan2((sqr(mhmin)-sqr(mh_)), mh_*wh_);
double rhomax = atan2((sqr(mhmax)-sqr(mh_)), mh_*wh_);
mh = sqrt(mh_*wh_*tan(rhomin+r[4]*(rhomax-rhomin))+sqr(mh_));
jacobian(jacobian()*(rhomax-rhomin));
}
if(rs<mh+2.*mq) return false;
// limits for virtual quark mass
Energy2 mmin(sqr(mq+mh)),mmax(sqr(rs-mq));
double rhomin,rhomax;
if(alpha_==0.) {
rhomin = mmin/sqr(mq);
rhomax = mmax/sqr(mq);
}
else if(alpha_==1.) {
rhomax = log((mmax-sqr(mq))/sqr(mq));
rhomin = log((mmin-sqr(mq))/sqr(mq));
}
else {
rhomin = pow((mmax-sqr(mq))/sqr(mq),1.-alpha_);
rhomax = pow((mmin-sqr(mq))/sqr(mq),1.-alpha_);
jacobian(jacobian()/(alpha_-1.));
}
// branch for mass smoothing
Energy2 m132,m232;
Energy p1,p2;
// first branch
if(r[1]<=0.5) {
double rtemp = 2.*r[1];
double rho = rhomin+rtemp*(rhomax-rhomin);
if(alpha_==0)
m132 = sqr(mq)*rho;
else if(alpha_==1)
m132 = sqr(mq)*(exp(rho)+1.);
else
m132 = sqr(mq)*(pow(rho,1./(1.-alpha_))+1.);
Energy m13 = sqrt(m132);
try {
p1 = SimplePhaseSpace::getMagnitude(sHat(), m13, mq);
p2 = SimplePhaseSpace::getMagnitude(m132,mq,mh);
} catch ( ImpossibleKinematics ) {
return false;
}
Energy ptmin = lastCuts().minKT(mePartonData()[3]);
double ctmin = -1.0, ctmax = 1.0;
if ( ptmin > ZERO ) {
double ctm = 1.0 - sqr(ptmin/p1);
if ( ctm <= 0.0 ) return false;
ctmin = max(ctmin, -sqrt(ctm));
ctmax = min(ctmax, sqrt(ctm));
}
double cos1 = getCosTheta(ctmin,ctmax,r[0]);
double sin1(sqrt(1.-sqr(cos1)));
double phi1 = Constants::twopi*UseRandom::rnd();
Lorentz5Momentum p13(sin1*p1*cos(phi1),sin1*p1*sin(phi1),cos1*p1,
sqrt(sqr(p1)+m132),m13);
meMomenta()[3].setVect(Momentum3(-sin1*p1*cos(phi1),-sin1*p1*sin(phi1),-cos1*p1));
meMomenta()[3].setMass(mq);
meMomenta()[3].rescaleEnergy();
bool test=Kinematics::twoBodyDecay(p13,mq,mh,-1.+2*r[2],r[3]*Constants::twopi,
meMomenta()[2],meMomenta()[4]);
if(!test) return false;
m232 = (meMomenta()[3]+meMomenta()[4]).m2();
double D = 2./(pow(sqr(mq)/(m132-sqr(mq)),alpha_)+
pow(sqr(mq)/(m232-sqr(mq)),alpha_));
jacobian(0.5*jacobian()*rs/m13*sqr(mq)*D*(rhomax-rhomin)/sHat());
}
// second branch
else {
double rtemp = 2.*(r[1]-0.5);
double rho = rhomin+rtemp*(rhomax-rhomin);
if(alpha_==0)
m232 = sqr(mq)*rho;
else if(alpha_==1)
m232 = sqr(mq)*(exp(rho)+1.);
else
m232 = sqr(mq)*(pow(rho,1./(1.-alpha_))+1.);
Energy m23 = sqrt(m232);
try {
p1 = SimplePhaseSpace::getMagnitude(sHat(), m23, mq);
p2 = SimplePhaseSpace::getMagnitude(m232,mq,mh);
} catch ( ImpossibleKinematics ) {
return false;
}
Energy ptmin = lastCuts().minKT(mePartonData()[2]);
double ctmin = -1.0, ctmax = 1.0;
if ( ptmin > ZERO ) {
double ctm = 1.0 - sqr(ptmin/p1);
if ( ctm <= 0.0 ) return false;
ctmin = max(ctmin, -sqrt(ctm));
ctmax = min(ctmax, sqrt(ctm));
}
double cos1 = getCosTheta(ctmin,ctmax,r[0]);
double sin1(sqrt(1.-sqr(cos1)));
double phi1 = Constants::twopi*UseRandom::rnd();
Lorentz5Momentum p23(-sin1*p1*cos(phi1),-sin1*p1*sin(phi1),-cos1*p1,
sqrt(sqr(p1)+m232),m23);
meMomenta()[2].setVect(Momentum3(sin1*p1*cos(phi1),sin1*p1*sin(phi1),cos1*p1));
meMomenta()[2].setMass(mq);
meMomenta()[2].rescaleEnergy();
bool test=Kinematics::twoBodyDecay(p23,mq,mh,-1.+2*r[2],r[3]*Constants::twopi,
meMomenta()[3],meMomenta()[4]);
if(!test) return false;
m132 = (meMomenta()[2]+meMomenta()[4]).m2();
double D = 2./(pow(sqr(mq)/(m132-sqr(mq)),alpha_)+
pow(sqr(mq)/(m232-sqr(mq)),alpha_));
jacobian(0.5*jacobian()*rs/m23*sqr(mq)*D*(rhomax-rhomin)/sHat());
}
// calculate jacobian
jacobian(0.125*jacobian()*p1*p2/sHat());
// check cuts
vector<LorentzMomentum> out;
tcPDVector tout;
for(unsigned int ix=2;ix<5;++ix) {
out .push_back(meMomenta()[ix]);
tout.push_back(mePartonData()[ix]);
}
return lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]);
}
CrossSection MEPP2QQHiggs::dSigHatDR() const {
using Constants::pi;
// jacobian factor for the higgs
InvEnergy2 bwfact;
Energy moff = meMomenta()[4].mass();
if(shapeOpt_==1) {
bwfact = mePartonData()[4]->generateWidth(moff)*moff/pi/
(sqr(sqr(moff)-sqr(mh_))+sqr(mh_*wh_));
}
else {
bwfact = hmass_->BreitWignerWeight(moff);
}
double jac1 = shapeOpt_==0 ?
1. : double(bwfact*(sqr(sqr(moff)-sqr(mh_))+sqr(mh_*wh_))/(mh_*wh_));
return sqr(hbarc)*me2()*jacobian()*jac1/sHat()/pow(Constants::twopi,3);
}
void MEPP2QQHiggs::getDiagrams() const {
tPDPtr Q = quark_[quarkFlavour_-1];
tPDPtr QB = antiquark_[quarkFlavour_-1];
// gg -> q qbar h0 subprocesses
if(process_==0||process_==1) {
// first t-channel
add(new_ptr((Tree2toNDiagram(3), gluon_, QB, gluon_,
1, Q, 4, Q , 2, QB, 4, higgs_, -1)));
add(new_ptr((Tree2toNDiagram(4), gluon_, QB, QB, gluon_,
1, Q, 3, QB, 2, higgs_, -2)));
add(new_ptr((Tree2toNDiagram(3),gluon_,QB,gluon_,
1, Q, 2, QB, 5, QB, 5, higgs_, -3)));
// interchange
add(new_ptr((Tree2toNDiagram(3),gluon_,Q,gluon_,
2, Q, 4, Q , 1, QB, 4, higgs_, -4)));
add(new_ptr((Tree2toNDiagram(4),gluon_,Q,Q,gluon_,
3, Q, 1, QB, 2, higgs_, -5)));
add(new_ptr((Tree2toNDiagram(3),gluon_,Q,gluon_,
2, Q, 1, QB, 5, QB, 5, higgs_, -6)));
// s-channel
add(new_ptr((Tree2toNDiagram(2),gluon_,gluon_, 1, gluon_,
3, Q, 4, Q, 3, QB, 4, higgs_, -7)));
add(new_ptr((Tree2toNDiagram(2),gluon_,gluon_, 1, gluon_,
3,Q, 3, QB, 5, QB, 5, higgs_, -8)));
}
// q qbar -> q qbar
if(process_==0||process_==2) {
for(unsigned int ix=1;ix<5;++ix) {
// gluon s-channel
add(new_ptr((Tree2toNDiagram(2),quark_[ix-1], antiquark_[ix-1],
1, gluon_, 3, Q, 4, Q , 3, QB, 4, higgs_, -9)));
add(new_ptr((Tree2toNDiagram(2),quark_[ix-1], antiquark_[ix-1],
1, gluon_, 3, Q, 3, QB, 5, QB, 5, higgs_, -10)));
}
}
}
double MEPP2QQHiggs::me2() const {
// total matrix element
double me(0.);
// gg initiated processes
if(mePartonData()[0]->id()==ParticleID::g) {
VectorWaveFunction g1w(meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction qw(meMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction qbarw(meMomenta()[3],mePartonData()[3],outgoing);
ScalarWaveFunction higgs(meMomenta()[4],mePartonData()[4],1.,outgoing);
vector<VectorWaveFunction> g1,g2;
vector<SpinorBarWaveFunction> q;
vector<SpinorWaveFunction> qbar;
for(unsigned int ix=0;ix<2;++ix) {
g1w.reset(2*ix);g1.push_back(g1w);
g2w.reset(2*ix);g2.push_back(g2w);
qw.reset(ix);q.push_back(qw);
qbarw.reset(ix);qbar.push_back(qbarw);
}
// calculate the matrix element
me=ggME(g1,g2,q,qbar,higgs,0);
}
// q qbar initiated
else {
SpinorWaveFunction q1w(meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction q2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction q3w(meMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction q4w(meMomenta()[3],mePartonData()[3],outgoing);
ScalarWaveFunction higgs(meMomenta()[4],mePartonData()[4],1.,outgoing);
vector<SpinorWaveFunction> q1,q4;
vector<SpinorBarWaveFunction> q2,q3;
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1.push_back(q1w);
q2w.reset(ix);q2.push_back(q2w);
q3w.reset(ix);q3.push_back(q3w);
q4w.reset(ix);q4.push_back(q4w);
}
// calculate the matrix element
me = qqME(q1,q2,q3,q4,higgs,0);
}
return me*sHat()*UnitRemoval::InvE2;
}
double MEPP2QQHiggs::ggME(vector<VectorWaveFunction> &g1,
vector<VectorWaveFunction> &g2,
vector<SpinorBarWaveFunction> & q,
vector<SpinorWaveFunction> & qbar,
ScalarWaveFunction & hwave,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
Energy mass = q[0].mass();
// matrix element to be stored
if(iflow!=0) me_.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin0));
// calculate the matrix element
double output(0.),sumflow[2]={0.,0.};
double sumdiag[8]={0.,0.,0.,0.,0.,0.,0.,0.};
Complex diag[8],flow[2];
VectorWaveFunction interv;
SpinorWaveFunction inters,QBoff;
SpinorBarWaveFunction intersb,Qoff;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
interv = GGGVertex_->evaluate(mt,5,gluon_,g1[ihel1],g2[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
- Qoff = QQHVertex_->evaluate(mt,3,q[ohel1].particle(),
+ Qoff = QQHVertex_->evaluate(mt,3,q[ohel1].particle()->CC(),
q[ohel1],hwave,mass);
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
- QBoff = QQHVertex_->evaluate(mt,3,qbar[ohel2].particle(),
+ QBoff = QQHVertex_->evaluate(mt,3,qbar[ohel2].particle()->CC(),
qbar[ohel2],hwave,mass);
// 1st diagram
- inters = QQGVertex_->evaluate(mt,1,qbar[ohel2].particle(),
+ inters = QQGVertex_->evaluate(mt,1,qbar[ohel2].particle()->CC(),
qbar[ohel2],g2[ihel2],mass);
diag[0] = QQGVertex_->evaluate(mt,inters,Qoff,g1[ihel1]);
// 2nd diagram
- intersb = QQGVertex_->evaluate(mt,1,q[ohel1].particle(),
+ intersb = QQGVertex_->evaluate(mt,1,q[ohel1].particle()->CC(),
q[ohel1],g1[ihel1],mass);
diag[1] = QQHVertex_->evaluate(mt,inters,intersb,hwave);
// 3rd diagram
diag[2] = QQGVertex_->evaluate(mt,QBoff,intersb,g2[ihel2]);
// 4th diagram
- inters = QQGVertex_->evaluate(mt,1,qbar[ohel2].particle(),
+ inters = QQGVertex_->evaluate(mt,1,qbar[ohel2].particle()->CC(),
qbar[ohel2],g1[ihel1],mass);
diag[3] = QQGVertex_->evaluate(mt,inters,Qoff,g2[ihel2]);
// 5th diagram
- intersb = QQGVertex_->evaluate(mt,1,q[ohel1].particle(),
+ intersb = QQGVertex_->evaluate(mt,1,q[ohel1].particle()->CC(),
q[ohel1],g2[ihel2],mass);
diag[4] = QQHVertex_->evaluate(mt,inters,intersb,hwave);
// 6th diagram
diag[5] = QQGVertex_->evaluate(mt,QBoff,intersb,g1[ihel1]);
// 7th diagram
diag[6] = QQGVertex_->evaluate(mt,qbar[ohel2],Qoff ,interv);
// 8th diagram
diag[7] = QQGVertex_->evaluate(mt,QBoff ,q[ohel1],interv);
// colour flows
flow[0]=diag[0]+diag[1]+diag[2]+(diag[6]+diag[7]);
flow[1]=diag[3]+diag[4]+diag[5]-(diag[6]+diag[7]);
// sums
for(unsigned int ix=0;ix<8;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) me_(2*ihel1,2*ihel2,ohel1,ohel2,0)=flow[iflow-1];
}
}
}
}
// select a colour flow
flow_ = 1 + UseRandom::rnd2(sumflow[0],sumflow[1]);
if(flow_==1) sumdiag[0]=sumdiag[1]=sumdiag[2]=0.;
else sumdiag[3]=sumdiag[4]=sumdiag[5]=0.;
// select a diagram from that flow
double prob = UseRandom::rnd();
for(unsigned int ix=0;ix<8;++ix) {
if(prob<=sumdiag[ix]) {
diagram_=1+ix;
break;
}
prob -= sumdiag[ix];
}
// final part of colour and spin factors
return output/48.;
}
double MEPP2QQHiggs::qqME(vector<SpinorWaveFunction> & q1,
vector<SpinorBarWaveFunction> & q2,
vector<SpinorBarWaveFunction> & q3,
vector<SpinorWaveFunction> & q4,
ScalarWaveFunction & hwave,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
Energy mass = q3[0].mass();
// matrix element to be stored
if(iflow!=0) me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin0));
// calculate the matrix element
double output(0.),sumdiag[2]={0.,0.};
Complex diag[2];
VectorWaveFunction interv;
SpinorWaveFunction QBoff;
SpinorBarWaveFunction Qoff;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
interv = QQGVertex_->evaluate(mt,5,gluon_,q1[ihel1],q2[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
- Qoff = QQHVertex_->evaluate(mt,3,q3[ohel1].particle(),
+ Qoff = QQHVertex_->evaluate(mt,3,q3[ohel1].particle()->CC(),
q3[ohel1],hwave,mass);
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
- QBoff = QQHVertex_->evaluate(mt,3,q4[ohel2].particle(),
+ QBoff = QQHVertex_->evaluate(mt,3,q4[ohel2].particle()->CC(),
q4[ohel2],hwave,mass);
// 1st diagram
diag[0] = QQGVertex_->evaluate(mt,q4[ohel2],Qoff,interv);
// 2nd diagram
diag[1] = QQGVertex_->evaluate(mt,QBoff,q3[ohel1],interv);
// sum of diagrams
for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]);
diag[0] += diag[1];
output += norm(diag[0]);
if(iflow!=0) me_(ihel1,ihel2,ohel1,ohel2,0) = diag[0];
}
}
}
}
// only 1 colour flow
flow_=1;
// select a diagram
diagram_ = 9+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// final part of colour and spin factors
return output/18.;
}
Selector<const ColourLines *>
MEPP2QQHiggs::colourGeometries(tcDiagPtr diag) const {
// colour lines for gg -> Q Qbar H
static const ColourLines cgg[10]=
{ColourLines("1 4 5, -1 -2 3 , -3 -6 "),
ColourLines("1 5 , -1 -2 -3 4, -4 -6 "),
ColourLines("1 4 , -1 -2 3 , -3 -5 -6"),
ColourLines("3 4 5, 1 2 -3 , -1 -6 "),
ColourLines("4 5 , 1 2 3 -4, -1 -6"),
ColourLines("3 4 , 1 2 -3 , -1 -5 -6"),
ColourLines("1 3 4 5, -1 2, -2 -3 -6"),
ColourLines("2 3 4 5, 1 -2, -1 -3 -6"),
ColourLines("1 3 4, -1 2, -2 -3 -5 -6"),
ColourLines("2 3 4, 1 -2, -1 -3 -5 -6")};
// colour lines for q qbar -> Q Qbar H
static const ColourLines cqq[2]=
{ColourLines("1 3 4 5, -2 -3 -6"),
ColourLines("1 3 4 , -2 -3 -5 -6")};
// select the colour flow (as all ready picked just insert answer)
Selector<const ColourLines *> sel;
switch(abs(diag->id())) {
// gg -> q qbar subprocess
case 1: case 2: case 3: case 4: case 5: case 6:
sel.insert(1.0, &cgg[abs(diag->id())-1]);
break;
case 7:
sel.insert(1.0, &cgg[5 + flow_]);
break;
case 8:
sel.insert(1.0, &cgg[7 + flow_]);
break;
// q qbar -> q qbar subprocess
case 9: case 10:
sel.insert(1.0, &cqq[abs(diag->id())-9]);
break;
}
return sel;
}
Selector<MEBase::DiagramIndex>
MEPP2QQHiggs::diagrams(const DiagramVector & diags) const {
// select the diagram, this is easy for us as we have already done it
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
if(diags[i]->id()==-int(diagram_)) sel.insert(1.0, i);
else sel.insert(0., i);
}
return sel;
}
void MEPP2QQHiggs::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard;
hard.push_back(sub->incoming().first);
hard.push_back(sub->incoming().second);
for(unsigned int ix=0;ix<3;++ix) hard.push_back(sub->outgoing()[ix]);
// identify the process and calculate the matrix element
if(hard[0]->id()<0) swap(hard[0],hard[1]);
if(hard[2]->id()==ParticleID::h0) swap(hard[2],hard[4]);
if(hard[3]->id()==ParticleID::h0) swap(hard[3],hard[4]);
if(hard[2]->id()<0) swap(hard[2],hard[3]);
if(hard[0]->id()==ParticleID::g) {
vector<VectorWaveFunction> g1,g2;
vector<SpinorBarWaveFunction> q;
vector<SpinorWaveFunction> qbar;
// off-shell wavefunctions for the spin correlations
VectorWaveFunction( g1,hard[0],incoming,false,true,true);
VectorWaveFunction( g2,hard[1],incoming,false,true,true);
SpinorBarWaveFunction(q ,hard[2],outgoing,true ,true);
SpinorWaveFunction( qbar,hard[3],outgoing,true ,true);
ScalarWaveFunction hwave( hard[4],outgoing,true);
g1[1]=g1[2];g2[1]=g2[2];
ggME(g1,g2,q,qbar,hwave,flow_);
}
// q qbar -> Q Qbar Higgs
else {
vector<SpinorWaveFunction> q1,q4;
vector<SpinorBarWaveFunction> q2,q3;
// off-shell for spin correlations
SpinorWaveFunction( q1,hard[0],incoming,false,true);
SpinorBarWaveFunction(q2,hard[1],incoming,false,true);
SpinorBarWaveFunction(q3,hard[2],outgoing,true ,true);
SpinorWaveFunction( q4,hard[3],outgoing,true ,true);
ScalarWaveFunction hwave( hard[4],outgoing,true);
qqME(q1,q2,q3,q4,hwave,flow_);
}
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(me_);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<5;++ix)
hard[ix]->spinInfo()->productionVertex(hardvertex);
}
diff --git a/MatrixElement/Hadron/MEPP2VGamma.cc b/MatrixElement/Hadron/MEPP2VGamma.cc
--- a/MatrixElement/Hadron/MEPP2VGamma.cc
+++ b/MatrixElement/Hadron/MEPP2VGamma.cc
@@ -1,382 +1,382 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEPP2VGamma class.
//
#include "MEPP2VGamma.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "Herwig/MatrixElement/HardVertex.h"
using namespace Herwig;
MEPP2VGamma::MEPP2VGamma() : process_(0), maxflavour_(5), massOption_(2)
{}
unsigned int MEPP2VGamma::orderInAlphaS() const {
return 0;
}
unsigned int MEPP2VGamma::orderInAlphaEW() const {
return 2;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEPP2VGamma,HwMEBase>
describeHerwigMEPP2VGamma("Herwig::MEPP2VGamma", "HwMEHadron.so");
void MEPP2VGamma::Init() {
static ClassDocumentation<MEPP2VGamma> documentation
("The MEPP2VGamma class simulates the production of"
" W+/-gamma and Z0gamma in hadron-hadron collisions"
" using the 2->2 matrix elements");
static Switch<MEPP2VGamma,unsigned int> interfaceProcess
("Process",
"Which processes to include",
&MEPP2VGamma::process_, 0, false, false);
static SwitchOption interfaceProcessAll
(interfaceProcess,
"All",
"Include all the processes",
0);
static SwitchOption interfaceProcessWGamma
(interfaceProcess,
"WGamma",
"Only include W+/-gamma",
1);
static SwitchOption interfaceProcessZGamma
(interfaceProcess,
"ZGamma",
"Only include ZGamma",
2);
static Parameter<MEPP2VGamma,int> interfaceMaximumFlavour
("MaximumFlavour",
"The maximum flavour allowed for the incoming quarks",
&MEPP2VGamma::maxflavour_, 5, 2, 5,
false, false, Interface::limited);
static Switch<MEPP2VGamma,unsigned int> interfaceMassOption
("MassOption",
"Option for the treatment of the boson masses",
&MEPP2VGamma::massOption_, 1, false, false);
static SwitchOption interfaceMassOptionOnMassShell
(interfaceMassOption,
"OnMassShell",
"The boson is produced on its mass shell",
1);
static SwitchOption interfaceMassOption2
(interfaceMassOption,
"OffShell",
"The bosons are generated off-shell using the mass and width generator.",
2);
}
void MEPP2VGamma::persistentOutput(PersistentOStream & os) const {
os << FFPvertex_ << FFWvertex_ << FFZvertex_ << WWWvertex_
<< process_ << massOption_;
}
void MEPP2VGamma::persistentInput(PersistentIStream & is, int) {
is >> FFPvertex_ >> FFWvertex_ >> FFZvertex_ >> WWWvertex_
>> process_ >> massOption_;
}
Energy2 MEPP2VGamma::scale() const {
return sHat();
}
IBPtr MEPP2VGamma::clone() const {
return new_ptr(*this);
}
IBPtr MEPP2VGamma::fullclone() const {
return new_ptr(*this);
}
void MEPP2VGamma::doinit() {
HwMEBase::doinit();
// mass option
vector<unsigned int> mopt(2,1);
mopt[0]=massOption_;
massOption(mopt);
rescalingOption(2);
// get the vertices we need
// get a pointer to the standard model object in the run
static const tcHwSMPtr hwsm
= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if (!hwsm) throw InitException() << "hwsm pointer is null in"
<< " MEPP2VGamma::doinit()"
<< Exception::abortnow;
// get pointers to all required Vertex objects
FFZvertex_ = hwsm->vertexFFZ();
FFPvertex_ = hwsm->vertexFFP();
WWWvertex_ = hwsm->vertexWWW();
FFWvertex_ = hwsm->vertexFFW();
}
Selector<const ColourLines *>
MEPP2VGamma::colourGeometries(tcDiagPtr diag) const {
static ColourLines cs("1 -2");
static ColourLines ct("1 2 -3");
Selector<const ColourLines *> sel;
if(diag->id()<-2) sel.insert(1.0, &cs);
else sel.insert(1.0, &ct);
return sel;
}
void MEPP2VGamma::getDiagrams() const {
typedef std::vector<pair<tcPDPtr,tcPDPtr> > Pairvector;
tcPDPtr wPlus = getParticleData(ParticleID::Wplus );
tcPDPtr wMinus = getParticleData(ParticleID::Wminus);
tcPDPtr z0 = getParticleData(ParticleID::Z0 );
tcPDPtr gamma = getParticleData(ParticleID::gamma);
// W+/- gamma
if(process_==0||process_==1) {
// possible parents
Pairvector parentpair;
parentpair.reserve(6);
// don't even think of putting 'break' in here!
switch(maxflavour_) {
case 5:
parentpair.push_back(make_pair(getParticleData(ParticleID::b),
getParticleData(ParticleID::cbar)));
parentpair.push_back(make_pair(getParticleData(ParticleID::b),
getParticleData(ParticleID::ubar)));
case 4:
parentpair.push_back(make_pair(getParticleData(ParticleID::s),
getParticleData(ParticleID::cbar)));
parentpair.push_back(make_pair(getParticleData(ParticleID::d),
getParticleData(ParticleID::cbar)));
case 3:
parentpair.push_back(make_pair(getParticleData(ParticleID::s),
getParticleData(ParticleID::ubar)));
case 2:
parentpair.push_back(make_pair(getParticleData(ParticleID::d),
getParticleData(ParticleID::ubar)));
default:
;
}
// W+ gamma
for(unsigned int ix=0;ix<parentpair.size();++ix) {
add(new_ptr((Tree2toNDiagram(3), parentpair[ix].second->CC(),
parentpair[ix].first, parentpair[ix].first->CC(),
1, wPlus, 2, gamma, -1)));
add(new_ptr((Tree2toNDiagram(3), parentpair[ix].second->CC(),
parentpair[ix].second->CC() , parentpair[ix].first->CC(),
2, wPlus, 1, gamma, -2)));
add(new_ptr((Tree2toNDiagram(2), parentpair[ix].second->CC(),
parentpair[ix].first->CC(), 1, wPlus, 3, wPlus, 3, gamma, -3)));
}
// W- gamma
for(unsigned int ix=0;ix<parentpair.size();++ix) {
add(new_ptr((Tree2toNDiagram(3), parentpair[ix].first,
parentpair[ix].second->CC(),
parentpair[ix].second, 1, wMinus, 2, gamma, -1)));
add(new_ptr((Tree2toNDiagram(3), parentpair[ix].first, parentpair[ix].first ,
parentpair[ix].second, 2, wMinus, 1, gamma, -2)));
add(new_ptr((Tree2toNDiagram(2), parentpair[ix].first,
parentpair[ix].second, 1, wMinus, 3, wMinus, 3, gamma, -3)));
}
}
if(process_==0||process_==2) {
for(int ix=1;ix<=maxflavour_;++ix) {
tcPDPtr qk = getParticleData(ix);
tcPDPtr qb = qk->CC();
add(new_ptr((Tree2toNDiagram(3), qk, qk, qb, 1, z0, 2, gamma, -1)));
add(new_ptr((Tree2toNDiagram(3), qk, qk, qb, 2, z0, 1, gamma, -2)));
}
}
}
Selector<MEBase::DiagramIndex>
MEPP2VGamma::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i )
sel.insert(meInfo()[abs(diags[i]->id()) - 1], i);
return sel;
}
double MEPP2VGamma::me2() const {
// setup momenta and particle data for the external wavefunctions
// incoming
SpinorWaveFunction em_in( meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction ep_in( meMomenta()[1],mePartonData()[1],incoming);
// outgoing
VectorWaveFunction v1_out(meMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction v2_out(meMomenta()[3],mePartonData()[3],outgoing);
vector<SpinorWaveFunction> f1;
vector<SpinorBarWaveFunction> a1;
vector<VectorWaveFunction> v1,v2;
// calculate the wavefunctions
for(unsigned int ix=0;ix<3;++ix) {
if(ix<2) {
em_in.reset(ix);
f1.push_back(em_in);
ep_in.reset(ix);
a1.push_back(ep_in);
}
v1_out.reset(ix);
v1.push_back(v1_out);
if(ix!=1) {
v2_out.reset(ix);
v2.push_back(v2_out);
}
}
if(mePartonData()[2]->id()==ParticleID::Z0) {
return ZGammaME(f1,a1,v1,v2,false);
}
else {
return WGammaME(f1,a1,v1,v2,false);
}
}
double MEPP2VGamma::ZGammaME(vector<SpinorWaveFunction> & f1,
vector<SpinorBarWaveFunction> & a1,
vector<VectorWaveFunction> & v1,
vector<VectorWaveFunction> & v2,
bool calc) const {
double output(0.);
vector<double> me(3,0.0);
if(calc) me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1,PDT::Spin1));
vector<Complex> diag(2,0.0);
SpinorWaveFunction inter;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<3;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
- inter = FFZvertex_->evaluate(scale(),5,f1[ihel1].particle(),
+ inter = FFZvertex_->evaluate(scale(),5,f1[ihel1].particle()->CC(),
f1[ihel1],v1[ohel1]);
diag[0] = FFPvertex_->evaluate(scale(),inter,a1[ihel2],v2[ohel2]);
- inter = FFPvertex_->evaluate(scale(),5,f1[ihel1].particle(),
+ inter = FFPvertex_->evaluate(scale(),5,f1[ihel1].particle()->CC(),
f1[ihel1] ,v2[ohel2]);
diag[1] = FFZvertex_->evaluate(scale(),inter,a1[ihel2],v1[ohel1]);
// individual diagrams
for (size_t ii=0; ii<2; ++ii) me[ii] += std::norm(diag[ii]);
// full matrix element
diag[0] += diag[1];
output += std::norm(diag[0]);
// storage of the matrix element for spin correlations
if(calc) me_(ihel1,ihel2,ohel1,ohel2) = diag[0];
}
}
}
}
DVector save(3);
for (size_t i = 0; i < 3; ++i) {
save[i] = 0.25 * me[i];
}
meInfo(save);
return 0.25*output/3.;
}
double MEPP2VGamma::WGammaME(vector<SpinorWaveFunction> & f1,
vector<SpinorBarWaveFunction> & a1,
vector<VectorWaveFunction> & v1,
vector<VectorWaveFunction> & v2,
bool calc) const {
double output(0.);
vector<double> me(3,0.0);
if(calc) me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1,PDT::Spin1));
vector<Complex> diag(3,0.0);
SpinorWaveFunction inter;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
VectorWaveFunction interW =
FFWvertex_->evaluate(scale(),3,v1[0].particle(),
f1[ihel1],a1[ihel2]);
for(unsigned int ohel1=0;ohel1<3;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// t-channel diagrams
inter = FFWvertex_->evaluate(scale(),5,a1[ihel1].particle(),
f1[ihel1],v1[ohel1]);
diag[0] = FFPvertex_->evaluate(scale(),inter,a1[ihel2],v2[ohel2]);
- inter = FFPvertex_->evaluate(scale(),5,f1[ihel1].particle(),
+ inter = FFPvertex_->evaluate(scale(),5,f1[ihel1].particle()->CC(),
f1[ihel1] ,v2[ohel2]);
diag[1] = FFWvertex_->evaluate(scale(),inter,a1[ihel2],v1[ohel1]);
// s-channel diagram
diag[2] = WWWvertex_->evaluate(scale(),interW,v1[ohel1],v2[ohel2]);
// individual diagrams
for (size_t ii=0; ii<3; ++ii) me[ii] += std::norm(diag[ii]);
// full matrix element
diag[0] += diag[1]+diag[2];
output += std::norm(diag[0]);
// storage of the matrix element for spin correlations
if(calc) me_(ihel1,ihel2,ohel1,ohel2) = diag[0];
}
}
}
}
DVector save(3);
for (size_t i = 0; i < 3; ++i) save[i] = 0.25 * me[i];
meInfo(save);
// spin and colour factors
output *= 0.25/3.;
// testing code
// int iu = abs(mePartonData()[0]->id());
// int id = abs(mePartonData()[1]->id());
// if(iu%2!=0) swap(iu,id);
// iu = (iu-2)/2;
// id = (id-1)/2;
// double ckm = SM().CKM(iu,id);
// InvEnergy4 dsigdt = Constants::pi*sqr(SM().alphaEM(scale()))
// /6./sqr(sHat())/SM().sin2ThetaW()*sqr(1./(1.+tHat()/uHat())-1./3.)*
// (sqr(tHat())+sqr(uHat())+2.*sqr(getParticleData(ParticleID::Wplus)->mass())*sHat())/
// tHat()/uHat();
// double test = 16.*ckm*Constants::pi*sqr(sHat())*dsigdt;
// cerr << "testing W gamma " << test << " " << output << " "
// << (test-output)/(test+output) << "\n";
return output;
}
void MEPP2VGamma::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard;
hard.push_back(sub->incoming().first);
hard.push_back(sub->incoming().second);
hard.push_back(sub->outgoing()[0]);
hard.push_back(sub->outgoing()[1]);
// order of particles
unsigned int order[4]={0,1,2,3};
if(hard[order[0]]->id()<0) swap(order[0],order[1]);
vector<SpinorWaveFunction> q;
vector<SpinorBarWaveFunction> qb;
SpinorWaveFunction (q ,hard[order[0]],incoming,false);
SpinorBarWaveFunction(qb,hard[order[1]],incoming,false);
vector<VectorWaveFunction> w1,w2;
if(hard[order[2]]->id()==ParticleID::gamma)
swap(order[2],order[3]);
VectorWaveFunction (w1,hard[order[2]],outgoing,true ,false);
VectorWaveFunction (w2,hard[order[3]],outgoing,true ,true );
w2[1]=w2[2];
// q qbar -> Z gamma
if(hard[order[2]]->id()==ParticleID::Z0) {
ZGammaME(q,qb,w1,w2,true);
}
// q qbar -> W gamma
else {
WGammaME(q,qb,w1,w2,true);
}
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(me_);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<4;++ix)
hard[order[ix]]->spinInfo()->productionVertex(hardvertex);
}
diff --git a/MatrixElement/Hadron/MEPP2WJet.cc b/MatrixElement/Hadron/MEPP2WJet.cc
--- a/MatrixElement/Hadron/MEPP2WJet.cc
+++ b/MatrixElement/Hadron/MEPP2WJet.cc
@@ -1,838 +1,838 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEPP2WJet class.
//
#include "MEPP2WJet.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Cuts/Cuts.h"
#include "Herwig/MatrixElement/HardVertex.h"
using namespace Herwig;
MEPP2WJet::MEPP2WJet() : _process(0), _maxflavour(5), _plusminus(0),
_wdec(0), _widthopt(1)
{}
void MEPP2WJet::doinit() {
HwMEBase::doinit();
_wplus = getParticleData(ThePEG::ParticleID::Wplus );
_wminus = getParticleData(ThePEG::ParticleID::Wminus);
// cast the SM pointer to the Herwig SM pointer
ThePEG::Ptr<Herwig::StandardModel>::transient_const_pointer
hwsm=ThePEG::dynamic_ptr_cast< ThePEG::Ptr<Herwig::StandardModel>
::transient_const_pointer>(standardModel());
// do the initialisation
if(!hwsm)
throw InitException() << "Must be Herwig::StandardModel in MEPP2WJet::doinit()"
<< Exception::runerror;
// set the vertex pointers
_theFFWVertex = hwsm->vertexFFW();
_theQQGVertex = hwsm->vertexFFG();
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEPP2WJet,HwMEBase>
describeHerwigMEPP2WJet("Herwig::MEPP2WJet", "HwMEHadron.so");
void MEPP2WJet::Init() {
static ClassDocumentation<MEPP2WJet> documentation
("The MEPP2WJet class implements the matrix element for W + jet production");
static Parameter<MEPP2WJet,unsigned int> interfaceMaxFlavour
( "MaxFlavour",
"The heaviest incoming quark flavour this matrix element is allowed to handle "
"(if applicable).",
&MEPP2WJet::_maxflavour, 5, 0, 8, false, false, true);
static Switch<MEPP2WJet,unsigned int> interfaceProcess
("Process",
"Which subprocesses to include",
&MEPP2WJet::_process, 0, false, false);
static SwitchOption interfaceProcessAll
(interfaceProcess,
"All",
"Include all subprocesses",
0);
static SwitchOption interfaceProcessqqbar
(interfaceProcess,
"qqbar",
"Only include q qbar -> W g process",
1);
static SwitchOption interfaceProcessqg
(interfaceProcess,
"qg",
"Only include the q g -> W q process",
2);
static SwitchOption interfaceProcessqbarg
(interfaceProcess,
"qbarg",
"Only include the qbar g -> W qbar process",
3);
static Switch<MEPP2WJet,unsigned int> interfacePlusMinus
("Wcharge",
"Which intermediate W bosons to include",
&MEPP2WJet::_plusminus, 0, false, false);
static SwitchOption interfacePlusMinusAll
(interfacePlusMinus,
"Both",
"Include W+ and W-",
0);
static SwitchOption interfacePlusMinusPlus
(interfacePlusMinus,
"Plus",
"Only include W+",
1);
static SwitchOption interfacePlusMinusMinus
(interfacePlusMinus,
"Minus",
"Only include W-",
2);
static Switch<MEPP2WJet,unsigned int> interfaceWDecay
("WDecay",
"Which processes to include",
&MEPP2WJet::_wdec, 0, false, false);
static SwitchOption interfaceWDecayAll
(interfaceWDecay,
"All",
"Include all SM fermions as outgoing particles",
0);
static SwitchOption interfaceWDecayQuarks
(interfaceWDecay,
"Quarks",
"Only include outgoing quarks",
1);
static SwitchOption interfaceWDecayLeptons
(interfaceWDecay,
"Leptons",
"All include outgoing leptons",
2);
static SwitchOption interfaceWDecayElectron
(interfaceWDecay,
"Electron",
"Only include outgoing e nu_e",
3);
static SwitchOption interfaceWDecayMuon
(interfaceWDecay,
"Muon",
"Only include outgoing mu nu_mu",
4);
static SwitchOption interfaceWDecayTau
(interfaceWDecay,
"Tau",
"Only include outgoing tauu nu_tau",
5);
static SwitchOption interfaceWDecayUpDown
(interfaceWDecay,
"UpDown",
"Only include outgoing u dbar/ d ubar",
6);
static SwitchOption interfaceWDecayUpStrange
(interfaceWDecay,
"UpStrange",
"Only include outgoing u sbar/ s ubar",
7);
static SwitchOption interfaceWDecayUpBottom
(interfaceWDecay,
"UpBottom",
"Only include outgoing u bbar/ b ubar",
8);
static SwitchOption interfaceWDecayCharmDown
(interfaceWDecay,
"CharmDown",
"Only include outgoing c dbar/ d cbar",
9);
static SwitchOption interfaceWDecayCharmStrange
(interfaceWDecay,
"CharmStrange",
"Only include outgoing c sbar/ s cbar",
10);
static SwitchOption interfaceWDecayCharmBottom
(interfaceWDecay,
"CharmBottom",
"Only include outgoing c bbar/ b cbar",
11);
static Switch<MEPP2WJet,unsigned int> interfaceWidthOption
("WidthOption",
"The option for handling the width of the off-shell W boson",
&MEPP2WJet::_widthopt, 1, false, false);
static SwitchOption interfaceWidthOptionFixedDenominator
(interfaceWidthOption,
"FixedDenominator",
"Use a fixed with in the W propagator but the full matrix element"
" in the numerator",
1);
static SwitchOption interfaceWidthOptionAllRunning
(interfaceWidthOption,
"AllRunning",
"Use a running width in the W propagator and the full matrix "
"element in the numerator",
2);
}
void MEPP2WJet::getDiagrams() const {
// which intgermediates to include
bool wplus = _plusminus==0 || _plusminus==1;
bool wminus = _plusminus==0 || _plusminus==2;
// possible incoming and outgoing particles
typedef std::vector<pair<long,long> > Pairvector;
// possible parents
Pairvector parentpair;
parentpair.reserve(6);
// don't even think of putting 'break' in here!
switch(_maxflavour) {
case 5:
parentpair.push_back(make_pair(ParticleID::b, ParticleID::cbar));
parentpair.push_back(make_pair(ParticleID::b, ParticleID::ubar));
case 4:
parentpair.push_back(make_pair(ParticleID::s, ParticleID::cbar));
parentpair.push_back(make_pair(ParticleID::d, ParticleID::cbar));
case 3:
parentpair.push_back(make_pair(ParticleID::s, ParticleID::ubar));
case 2:
parentpair.push_back(make_pair(ParticleID::d, ParticleID::ubar));
default:
;
}
// possible children
Pairvector childpair;
childpair.reserve(9);
childpair.push_back(make_pair(ParticleID::eminus, ParticleID::nu_ebar));
childpair.push_back(make_pair(ParticleID::muminus, ParticleID::nu_mubar));
childpair.push_back(make_pair(ParticleID::tauminus, ParticleID::nu_taubar));
childpair.push_back(make_pair(ParticleID::d, ParticleID::ubar));
childpair.push_back(make_pair(ParticleID::s, ParticleID::ubar));
childpair.push_back(make_pair(ParticleID::b, ParticleID::ubar));
childpair.push_back(make_pair(ParticleID::d, ParticleID::cbar));
childpair.push_back(make_pair(ParticleID::s, ParticleID::cbar));
childpair.push_back(make_pair(ParticleID::b, ParticleID::cbar));
// gluon for diagrams
tcPDPtr g = getParticleData(ParticleID::g);
// loop over the children
bool lepton,quark;
Pairvector::const_iterator child = childpair.begin();
for (; child != childpair.end(); ++child) {
// allowed leptonic decay
lepton=child->first>10&&
(_wdec==0||_wdec==2||
(abs(child->first)-5)/2==int(_wdec));
// allowed quark decay
quark =abs(child->second)<10&&
(_wdec==0||_wdec==1||
(abs(child->second)==2&&(abs(child->first)+11)/2==int(_wdec))||
(abs(child->second)==4&&(abs(child->first)+17)/2==int(_wdec)));
// if decay not allowed skip
if(!(quark||lepton)) continue;
// decay products
tcPDPtr lNeg1 = getParticleData(child->first);
tcPDPtr lNeg2 = getParticleData(child->second);
tcPDPtr lPos1 = lNeg2->CC();
tcPDPtr lPos2 = lNeg1->CC();
Pairvector::const_iterator parent = parentpair.begin();
for (; parent != parentpair.end(); ++parent) {
// parents
tcPDPtr qNeg1 = getParticleData(parent->first);
tcPDPtr qNeg2 = getParticleData(parent->second);
tcPDPtr qPos1 = qNeg2->CC();
tcPDPtr qPos2 = qNeg1->CC();
// diagrams
// q qbar annhilation processes
if(_process==0||_process==1) {
// q qbar -> W- g
if(wminus) {
add(new_ptr((Tree2toNDiagram(3), qNeg1, qNeg2, qNeg2, 1, _wminus,
2, g, 4, lNeg1, 4, lNeg2, -1)));
add(new_ptr((Tree2toNDiagram(3), qNeg1, qNeg1, qNeg2, 2, _wminus,
1, g, 4, lNeg1, 4, lNeg2, -2)));
}
// q qbar -> W+ g
if(wplus) {
add(new_ptr((Tree2toNDiagram(3), qPos1, qPos2, qPos2, 1, _wplus,
2, g, 4, lPos1, 4, lPos2, -3)));
add(new_ptr((Tree2toNDiagram(3), qPos1, qPos1, qPos2, 2, _wplus,
1, g, 4, lPos1, 4, lPos2, -4)));
}
}
// q g compton
if(_process==0||_process==2) {
if(wminus) {
add(new_ptr((Tree2toNDiagram(3), qNeg1, qPos1, g , 1, _wminus,
2, qPos1, 4, lNeg1, 4, lNeg2, -5)));
add(new_ptr((Tree2toNDiagram(2), qNeg1, g, 1, qNeg1, 3, _wminus,
3, qPos1, 4, lNeg1, 4, lNeg2, -6)));
}
if(wplus) {
add(new_ptr((Tree2toNDiagram(3), qPos1, qNeg1, g, 1, _wplus,
2, qNeg1, 4, lPos1, 4, lPos2, -7)));
add(new_ptr((Tree2toNDiagram(2), qPos1, g, 1, qNeg1, 3, _wplus,
3, qNeg1, 4, lPos1, 4, lPos2, -8)));
}
}
// qbar g compton
if(_process==0||_process==3) {
if(wminus) {
add(new_ptr((Tree2toNDiagram(3), qNeg2, qPos2, g, 1, _wminus,
2, qPos2, 4, lNeg1, 4, lNeg2, -9 )));
add(new_ptr((Tree2toNDiagram(2), qNeg2, g, 1, qNeg2, 3, _wminus,
3, qPos2, 4, lNeg1, 4, lNeg2, -10)));
}
if(wplus) {
add(new_ptr((Tree2toNDiagram(3), qPos2, qNeg2, g, 1, _wplus,
2, qNeg2, 4, lPos1, 4, lPos2, -11)));
add(new_ptr((Tree2toNDiagram(2), qPos2, g, 1, qPos2, 3, _wplus,
3, qNeg2, 4, lPos1, 4, lPos2, -12)));
}
}
}
}
}
unsigned int MEPP2WJet::orderInAlphaS() const {
return 1;
}
unsigned int MEPP2WJet::orderInAlphaEW() const {
return 2;
}
void MEPP2WJet::persistentOutput(PersistentOStream & os) const {
os << _theFFWVertex << _theQQGVertex << _wplus << _widthopt
<< _wminus << _process << _maxflavour << _plusminus << _wdec;
}
void MEPP2WJet::persistentInput(PersistentIStream & is, int) {
is >> _theFFWVertex >> _theQQGVertex >> _wplus >> _widthopt
>> _wminus >> _process >> _maxflavour >> _plusminus >> _wdec;
}
int MEPP2WJet::nDim() const {
return 5;
}
Selector<const ColourLines *>
MEPP2WJet::colourGeometries(tcDiagPtr diag) const {
// colour lines for q qbar -> W g
static const ColourLines cqqbar[4]={ColourLines("1 -2 5,-3 -5"),
ColourLines("1 5, -5 2 -3"),
ColourLines("1 -2 5,-3 -5,6 -7"),
ColourLines("1 5, -5 2 -3,6 -7")};
// colour lines for q g -> W q
static const ColourLines cqg [4]={ColourLines("1 2 -3,3 5"),
ColourLines("1 -2,2 3 5"),
ColourLines("1 2 -3,3 5,6 -7"),
ColourLines("1 -2,2 3 5,6 -7")};
// colour lines for qbar q -> W qbar
static const ColourLines cqbarg[4]={ColourLines("-1 -2 3,-3 -5"),
ColourLines("-1 2,-2 -3 -5"),
ColourLines("-1 -2 3,-3 -5,6 -7"),
ColourLines("-1 2,-2 -3 -5,6 -7")};
// select the correct line
unsigned int icol = mePartonData()[3]->coloured() ? 2 : 0;
Selector<const ColourLines *> sel;
switch(abs(diag->id())) {
case 1 : case 3:
sel.insert(1.0, &cqqbar[icol]);
break;
case 2 : case 4:
sel.insert(1.0, &cqqbar[icol+1]);
break;
case 5 : case 7:
sel.insert(1.0, &cqg[icol]);
break;
case 6 : case 8:
sel.insert(1.0, &cqg[icol+1]);
break;
case 9 : case 11:
sel.insert(1.0, &cqbarg[icol]);
break;
case 10 : case 12:
sel.insert(1.0, &cqbarg[icol+1]);
break;
}
return sel;
}
Selector<MEBase::DiagramIndex>
MEPP2WJet::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
int id=abs(diags[i]->id());
if (id <= 2 ) sel.insert(meInfo()[id- 1],i);
else if(id <= 4 ) sel.insert(meInfo()[id- 3],i);
else if(id <= 6 ) sel.insert(meInfo()[id- 5],i);
else if(id <= 8 ) sel.insert(meInfo()[id- 7],i);
else if(id <= 10) sel.insert(meInfo()[id- 9],i);
else if(id <= 12) sel.insert(meInfo()[id-11],i);
}
return sel;
}
Energy2 MEPP2WJet::scale() const {
return _scale;
}
CrossSection MEPP2WJet::dSigHatDR() const {
return me2()*jacobian()/(16.0*sqr(Constants::pi)*sHat())*sqr(hbarc);
}
bool MEPP2WJet::generateKinematics(const double * r) {
// initialize jacobian
jacobian(1.);
// cms energy
Energy ecm=sqrt(sHat());
// find the right W pointer
tcPDPtr wdata = mePartonData()[3]->iCharge()+mePartonData()[4]->iCharge() > 0 ?
_wplus :_wminus;
// first generate the mass of the off-shell gauge boson
// minimum mass of the
tcPDVector ptemp;
ptemp.push_back(mePartonData()[3]);
ptemp.push_back(mePartonData()[4]);
Energy2 minMass2 = max(lastCuts().minSij(mePartonData()[3],mePartonData()[4]),
lastCuts().minS(ptemp));
// minimum pt of the jet
Energy ptmin = max(lastCuts().minKT(mePartonData()[2]),
lastCuts().minKT(wdata));
// maximum mass of the gauge boson so pt is possible
Energy2 maxMass2 = min(ecm*(ecm-2.*ptmin),lastCuts().maxS(ptemp));
if(maxMass2<=ZERO||minMass2<ZERO) return false;
// also impose the limits from the ParticleData object
minMass2 = max(minMass2,sqr(wdata->massMin()));
maxMass2 = min(maxMass2,sqr(wdata->massMax()));
// return if not kinematically possible
if(minMass2>maxMass2) return false;
// generation of the mass
Energy M(wdata->mass()),Gamma(wdata->width());
Energy2 M2(sqr(M)),MG(M*Gamma);
double rhomin = atan2((minMass2-M2),MG);
double rhomax = atan2((maxMass2-M2),MG);
_mw2=M2+MG*tan(rhomin+r[1]*(rhomax-rhomin));
Energy mw=sqrt(_mw2);
// jacobian
jacobian(jacobian()*(sqr(_mw2-M2)+sqr(MG))/MG*(rhomax-rhomin)/sHat());
// set the masses of the outgoing particles in the 2-2 scattering
meMomenta()[2].setMass(ZERO);
Lorentz5Momentum pw(mw);
// generate the polar angle of the hard scattering
double ctmin(-1.0), ctmax(1.0);
Energy q(ZERO);
try {
q = SimplePhaseSpace::getMagnitude(sHat(), meMomenta()[2].mass(),mw);
}
catch ( ImpossibleKinematics ) {
return false;
}
Energy2 pq = sqrt(sHat())*q;
if ( ptmin > ZERO ) {
double ctm = 1.0 - sqr(ptmin/q);
if ( ctm <= 0.0 ) return false;
ctmin = max(ctmin, -sqrt(ctm));
ctmax = min(ctmax, sqrt(ctm));
}
if ( ctmin >= ctmax ) return false;
double cth = getCosTheta(ctmin, ctmax, r[0]);
// momenta of particle in hard scattering
Energy pt = q*sqrt(1.0-sqr(cth));
double phi=2.0*Constants::pi*r[2];
meMomenta()[2].setVect(Momentum3( pt*sin(phi), pt*cos(phi), q*cth));
pw.setVect( Momentum3(-pt*sin(phi),-pt*cos(phi),-q*cth));
meMomenta()[2].rescaleEnergy();
pw.rescaleEnergy();
// set the scale
_scale = _mw2+sqr(pt);
// generate the momenta of the W decay products
meMomenta()[3].setMass(mePartonData()[3]->mass());
meMomenta()[4].setMass(mePartonData()[4]->mass());
Energy q2 = ZERO;
try {
q2 = SimplePhaseSpace::getMagnitude(_mw2, meMomenta()[3].mass(),
meMomenta()[4].mass());
} catch ( ImpossibleKinematics ) {
return false;
}
double cth2 =-1.+2.*r[3];
double phi2=Constants::twopi*r[4];
Energy pt2 =q2*sqrt(1.-sqr(cth2));
Lorentz5Momentum pl[2]={Lorentz5Momentum( pt2*cos(phi2), pt2*sin(phi2), q2*cth2,ZERO,
meMomenta()[3].mass()),
Lorentz5Momentum(-pt2*cos(phi2),-pt2*sin(phi2),-q2*cth2,ZERO,
meMomenta()[4].mass())};
pl[0].rescaleEnergy();
pl[1].rescaleEnergy();
Boost boostv(pw.boostVector());
pl[0].boost(boostv);
pl[1].boost(boostv);
meMomenta()[3] = pl[0];
meMomenta()[4] = pl[1];
// check passes all the cuts
vector<LorentzMomentum> out(3);
out[0] = meMomenta()[2];
out[1] = meMomenta()[3];
out[2] = meMomenta()[4];
tcPDVector tout(3);
tout[0] = mePartonData()[2];
tout[1] = mePartonData()[3];
tout[2] = mePartonData()[4];
if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) )
return false;
// jacobian
jacobian((pq/sHat())*Constants::pi*jacobian()/8./sqr(Constants::pi)*q2/mw);
return true;
}
double MEPP2WJet::me2() const {
InvEnergy2 output(ZERO);
// construct spinors for the leptons (always the same)
vector<SpinorBarWaveFunction> lm;
vector<SpinorWaveFunction> lp;
SpinorBarWaveFunction lmout(meMomenta()[3],mePartonData()[3],outgoing);
SpinorWaveFunction lpout(meMomenta()[4],mePartonData()[4],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
lmout.reset(ix);lm.push_back(lmout);
lpout.reset(ix);lp.push_back(lpout);
}
// q g to q W
if(mePartonData()[0]->id()<=6&&mePartonData()[0]->id()>0&&
mePartonData()[1]->id()==ParticleID::g) {
// polarization states for the particles
vector<SpinorWaveFunction> fin;
vector<VectorWaveFunction> gin;
vector<SpinorBarWaveFunction> fout;
SpinorWaveFunction qin (meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction glin(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction qout(meMomenta()[2],mePartonData()[2],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
qin.reset(ix) ; fin.push_back(qin);
glin.reset(2*ix); gin.push_back(glin);
qout.reset(ix);fout.push_back(qout);
}
output=qgME(fin,gin,fout,lm,lp);
}
// qbar g to qbar W
else if(mePartonData()[0]->id()>=-6&&mePartonData()[0]->id()<0&&
mePartonData()[1]->id()==ParticleID::g) {
vector<SpinorBarWaveFunction> ain;
vector<VectorWaveFunction> gin;
vector<SpinorWaveFunction> aout;
SpinorBarWaveFunction qbin (meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction glin (meMomenta()[1],mePartonData()[1],incoming);
SpinorWaveFunction qbout(meMomenta()[2],mePartonData()[2],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
qbin.reset(ix) ; ain.push_back(qbin);
glin.reset(2*ix) ; gin.push_back(glin);
qbout.reset(ix);aout.push_back(qbout);
}
output=qbargME(ain,gin,aout,lm,lp);
}
// q qbar to g W
else {
vector<SpinorWaveFunction> fin;
vector<SpinorBarWaveFunction> ain;
vector<VectorWaveFunction> gout;
SpinorWaveFunction qin (meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction qbin(meMomenta()[1],mePartonData()[1],incoming);
VectorWaveFunction glout(meMomenta()[2],mePartonData()[2],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
qin.reset(ix) ; fin.push_back(qin);
qbin.reset(ix) ; ain.push_back(qbin);
glout.reset(2*ix); gout.push_back(glout);
}
output=qqbarME(fin,ain,gout,lm,lp);
}
return output*sHat();
}
InvEnergy2 MEPP2WJet::qqbarME(vector<SpinorWaveFunction> & fin,
vector<SpinorBarWaveFunction> & ain,
vector<VectorWaveFunction> & gout,
vector<SpinorBarWaveFunction> & lm,
vector<SpinorWaveFunction> & lp,
bool calc) const {
// if calculation spin corrections construct the me
if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1,PDT::Spin1Half,
PDT::Spin1Half));
// some integers
unsigned int ihel1,ihel2,ohel1,ohel2,ohel3;
// find the right W pointer
tcPDPtr wdata = mePartonData()[3]->iCharge()+mePartonData()[4]->iCharge() > 0
? _wplus :_wminus;
// compute the W current for speed
VectorWaveFunction bcurr[2][2];
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
bcurr[ohel2][ohel3] = _theFFWVertex->evaluate(_mw2,_widthopt,wdata,
lp[ohel3],lm[ohel2]);
}
}
double me[3]={0.,0.,0.};
Complex diag[2];
SpinorWaveFunction inters;
SpinorBarWaveFunction interb;
for(ihel1=0;ihel1<2;++ihel1) {
for(ihel2=0;ihel2<2;++ihel2) {
for(ohel1=0;ohel1<2;++ohel1) {
// intermediates for the diagrams
inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[0],
fin[ihel1],gout[ohel1]);
interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[1],
ain[ihel2],gout[ohel1]);
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
diag[0] = _theFFWVertex->evaluate(_mw2,fin[ihel1],interb,
bcurr[ohel2][ohel3]);
diag[1] = _theFFWVertex->evaluate(_mw2,inters,ain[ihel2],
bcurr[ohel2][ohel3]);
// diagram contributions
me[1] += norm(diag[0]);
me[2] += norm(diag[1]);
// total
diag[0] += diag[1];
me[0] += norm(diag[0]);
if(calc) _me(ihel1,ihel2,2*ohel1,ohel2,ohel3) = diag[0];
}
}
}
}
}
// results
// initial state spin and colour average
double colspin=1./9./4.;
// and C_F N_c from matrix element
colspin *= 4.;
// colour factor for the W decay
if(mePartonData()[3]->coloured()) colspin*=3.;
DVector save;
for(unsigned int ix=0;ix<3;++ix) {
me[ix]*=colspin;
if(ix>0) save.push_back(me[ix]);
}
meInfo(save);
return me[0] * UnitRemoval::InvE2;
}
InvEnergy2 MEPP2WJet::qgME(vector<SpinorWaveFunction> & fin,
vector<VectorWaveFunction> & gin,
vector<SpinorBarWaveFunction> & fout,
vector<SpinorBarWaveFunction> & lm,
vector<SpinorWaveFunction> & lp,
bool calc) const {
// if calculation spin corrections construct the me
if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half));
// find the right W pointer
tcPDPtr wdata = mePartonData()[3]->iCharge()+mePartonData()[4]->iCharge() > 0 ?
_wplus :_wminus;
// some integers
unsigned int ihel1,ihel2,ohel1,ohel2,ohel3;
// compute the leptonic W current for speed
VectorWaveFunction bcurr[2][2];
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
bcurr[ohel2][ohel3] = _theFFWVertex->evaluate(_mw2,_widthopt,wdata,
lp[ohel3],lm[ohel2]);
}
}
// compute the matrix elements
double me[3]={0.,0.,0.};
Complex diag[2];
SpinorWaveFunction inters;
SpinorBarWaveFunction interb;
for(ihel1=0;ihel1<2;++ihel1) {
for(ihel2=0;ihel2<2;++ihel2) {
for(ohel1=0;ohel1<2;++ohel1) {
// intermediates for the diagrams
interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[2],
fout[ohel1],gin[ihel2]);
inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[0],
fin[ihel1],gin[ihel2]);
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
diag[0]=_theFFWVertex->evaluate(_mw2,fin[ihel1],interb,
bcurr[ohel2][ohel3]);
diag[1]=_theFFWVertex->evaluate(_mw2,inters,fout[ohel1],
bcurr[ohel2][ohel3]);
// diagram contributions
me[1] += norm(diag[0]);
me[2] += norm(diag[1]);
// total
diag[0] += diag[1];
me[0] += norm(diag[0]);
if(calc) _me(ihel1,2*ihel2,ohel1,ohel2,ohel3) = diag[0];
}
}
}
}
}
// results
// initial state spin and colour average
double colspin=1./24./4.;
// and C_F N_c from matrix element
colspin *=4.;
// colour factor for the W decay
if(mePartonData()[3]->coloured()) colspin*=3.;
DVector save;
for(unsigned int ix=0;ix<3;++ix) {
me[ix]*=colspin;
if(ix>0) save.push_back(me[ix]);
}
meInfo(save);
return me[0] * UnitRemoval::InvE2;
}
InvEnergy2 MEPP2WJet::qbargME(vector<SpinorBarWaveFunction> & fin,
vector<VectorWaveFunction> & gin,
vector<SpinorWaveFunction> & fout,
vector<SpinorBarWaveFunction> & lm,
vector<SpinorWaveFunction> & lp,
bool calc) const {
// if calculation spin corrections construct the me
if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half));
// find the right W pointer
tcPDPtr wdata = mePartonData()[3]->iCharge()+mePartonData()[4]->iCharge() > 0 ?
_wplus :_wminus;
// some integers
unsigned int ihel1,ihel2,ohel1,ohel2,ohel3;
// compute the leptonic W current for speed
VectorWaveFunction bcurr[2][2];
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
bcurr[ohel2][ohel3] = _theFFWVertex->evaluate(_mw2,_widthopt,wdata,
lp[ohel3],lm[ohel2]);
}
}
// compute the matrix elements
double me[3]={0.,0.,0.};
Complex diag[2];
SpinorWaveFunction inters;
SpinorBarWaveFunction interb;
for(ihel1=0;ihel1<2;++ihel1) {
for(ihel2=0;ihel2<2;++ihel2) {
for(ohel1=0;ohel1<2;++ohel1) {
// intermediates for the diagrams
- inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[2],
+ inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[2]->CC(),
fout[ohel1],gin[ihel2]);
interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[0],
fin[ihel1],gin[ihel2]);
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
diag[0]= _theFFWVertex->evaluate(_mw2,inters,fin[ihel1],
bcurr[ohel2][ohel3]);
diag[1]= _theFFWVertex->evaluate(_mw2,fout[ohel1],interb,
bcurr[ohel2][ohel3]);
// diagram contributions
me[1] += norm(diag[0]);
me[2] += norm(diag[1]);
// total
diag[0] += diag[1];
me[0] += norm(diag[0]);
if(calc) _me(ihel1,2*ihel2,ohel1,ohel2,ohel3) = diag[0];
}
}
}
}
}
// results
// initial state spin and colour average
double colspin=1./24./4.;
// and C_F N_c from matrix element
colspin *= 4.;
// colour factor for the W decay
if(mePartonData()[3]->coloured()) colspin*=3.;
DVector save;
for(unsigned int ix=0;ix<3;++ix) {
me[ix]*=colspin;
if(ix>0) save.push_back(me[ix]);
}
meInfo(save);
return me[0] * UnitRemoval::InvE2;
}
void MEPP2WJet::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard(5);
// incoming
hard[0]=sub->incoming().first;
hard[1]=sub->incoming().second;
if((hard[0]->id()<0&&hard[1]->id()<=6)||
hard[0]->id()==ParticleID::g) swap(hard[0],hard[1]);
// outgoing
for(unsigned int ix=0;ix<3;++ix) {
unsigned int iloc;
PPtr mother=sub->outgoing()[ix]->parents()[0];
if(mother&&abs(mother->id())==ParticleID::Wplus) {
if(sub->outgoing()[ix]->id()>0) iloc=3;
else iloc=4;
}
else iloc=2;
hard[iloc]=sub->outgoing()[ix];
}
// wavefunctions for the W decay products
vector<SpinorBarWaveFunction> lm;
vector<SpinorWaveFunction> lp;
SpinorBarWaveFunction(lm,hard[3],outgoing,true,true);
SpinorWaveFunction (lp,hard[4],outgoing,true,true);
// identify hard process and calculate matrix element
// q g to q W
if(hard[0]->id()<=6&&hard[0]->id()>0&&hard[1]->id()==ParticleID::g) {
vector<SpinorWaveFunction> fin;
vector<VectorWaveFunction> gin;
vector<SpinorBarWaveFunction> fout;
SpinorWaveFunction (fin ,hard[0],incoming,false,true);
VectorWaveFunction (gin ,hard[1],incoming,false,true,true);
SpinorBarWaveFunction (fout,hard[2],outgoing,true ,true);
gin[1]=gin[2];
qgME(fin,gin,fout,lm,lp,true);
}
// qbar g to qbar W
else if(hard[0]->id()>=-6&&hard[0]->id()<0&&hard[1]->id()==ParticleID::g) {
vector<SpinorBarWaveFunction> ain;
vector<VectorWaveFunction> gin;
vector<SpinorWaveFunction> aout;
SpinorBarWaveFunction(ain ,hard[0],incoming,false,true);
VectorWaveFunction (gin ,hard[1],incoming,false,true,true);
SpinorWaveFunction (aout,hard[2],outgoing,true ,true);
gin[1]=gin[2];
qbargME(ain,gin,aout,lm,lp,true);
}
// q qbar to g W
else {
vector<SpinorWaveFunction> fin;
vector<SpinorBarWaveFunction> ain;
vector<VectorWaveFunction> gout;
SpinorWaveFunction (fin ,hard[0],incoming,false,true);
SpinorBarWaveFunction(ain ,hard[1],incoming,false,true);
VectorWaveFunction (gout,hard[2],outgoing,true ,true,true);
gout[1]=gout[2];
qqbarME(fin,ain,gout,lm,lp,true);
}
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(_me);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<5;++ix)
(hard[ix]->spinInfo())->productionVertex(hardvertex);
}
diff --git a/MatrixElement/Hadron/MEPP2ZJet.cc b/MatrixElement/Hadron/MEPP2ZJet.cc
--- a/MatrixElement/Hadron/MEPP2ZJet.cc
+++ b/MatrixElement/Hadron/MEPP2ZJet.cc
@@ -1,861 +1,861 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEPP2ZJet class.
//
#include "MEPP2ZJet.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Cuts/Cuts.h"
#include "Herwig/MatrixElement/HardVertex.h"
using namespace Herwig;
MEPP2ZJet::MEPP2ZJet() : _process(0), _maxflavour(5), _zdec(0),
_gammaZ(0), _widthopt(1), _pprob(0.5)
{}
void MEPP2ZJet::doinit() {
HwMEBase::doinit();
_z0 = getParticleData(ThePEG::ParticleID::Z0 );
_gamma = getParticleData(ThePEG::ParticleID::gamma);
// cast the SM pointer to the Herwig SM pointer
ThePEG::Ptr<Herwig::StandardModel>::transient_const_pointer
hwsm=ThePEG::dynamic_ptr_cast< ThePEG::Ptr<Herwig::StandardModel>
::transient_const_pointer>(standardModel());
// do the initialisation
if(!hwsm)
throw InitException() << "Must be Herwig::StandardModel in MEPP2ZJet::doinit()"
<< Exception::runerror;
// set the vertex pointers
_theFFZVertex = hwsm->vertexFFZ();
_theFFPVertex = hwsm->vertexFFP();
_theQQGVertex = hwsm->vertexFFG();
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEPP2ZJet,HwMEBase>
describeHerwigMEPP2ZJet("Herwig::MEPP2ZJet", "HwMEHadron.so");
void MEPP2ZJet::Init() {
static ClassDocumentation<MEPP2ZJet> documentation
("The MEPP2ZJet class implements the matrix element for Z/gamma+ jet production");
static Parameter<MEPP2ZJet,int> interfaceMaxFlavour
( "MaxFlavour",
"The heaviest incoming quark flavour this matrix element is allowed to handle "
"(if applicable).",
&MEPP2ZJet::_maxflavour, 5, 0, 8, false, false, true);
static Switch<MEPP2ZJet,int> interfaceZDecay
("ZDecay",
"Which process to included",
&MEPP2ZJet::_zdec, 0, false, false);
static SwitchOption interfaceZDecayAll
(interfaceZDecay,
"All",
"Include all SM fermions as outgoing particles",
0);
static SwitchOption interfaceZDecayQuarks
(interfaceZDecay,
"Quarks",
"All include the quarks as outgoing particles",
1);
static SwitchOption interfaceZDecayLeptons
(interfaceZDecay,
"Leptons",
"Only include the leptons as outgoing particles",
2);
static SwitchOption interfaceZDecayChargedLeptons
(interfaceZDecay,
"ChargedLeptons",
"Only include the charged leptons as outgoing particles",
3);
static SwitchOption interfaceZDecayNeutrinos
(interfaceZDecay,
"Neutrinos",
"Only include the neutrinos as outgoing particles",
4);
static SwitchOption interfaceZDecayElectron
(interfaceZDecay,
"Electron",
"Only include e+e- as outgoing particles",
5);
static SwitchOption interfaceZDecayMuon
(interfaceZDecay,
"Muon",
"Only include mu+mu- as outgoing particles",
6);
static SwitchOption interfaceZDecayTau
(interfaceZDecay,
"Tau",
"Only include tau+tau- as outgoing particles",
7);
static SwitchOption interfaceZDecayNu_e
(interfaceZDecay,
"Nu_e",
"Only include nu_e ne_ebar as outgoing particles",
8);
static SwitchOption interfaceZDecaynu_mu
(interfaceZDecay,
"Nu_mu",
"Only include nu_mu nu_mubar as outgoing particles",
9);
static SwitchOption interfaceZDecaynu_tau
(interfaceZDecay,
"Nu_tau",
"Only include nu_tau nu_taubar as outgoing particles",
10);
static SwitchOption interfaceZDecayDown
(interfaceZDecay,
"Down",
"Only include d dbar as outgoing particles",
11);
static SwitchOption interfaceZDecayUp
(interfaceZDecay,
"Up",
"Only include u ubar as outgoing particles",
12);
static SwitchOption interfaceZDecayStrange
(interfaceZDecay,
"Strange",
"Only include s sbar as outgoing particles",
13);
static SwitchOption interfaceZDecayCharm
(interfaceZDecay,
"Charm",
"Only include c cbar as outgoing particles",
14);
static SwitchOption interfaceZDecayBottom
(interfaceZDecay,
"Bottom",
"Only include b bbar as outgoing particles",
15);
static SwitchOption interfaceZDecayTop
(interfaceZDecay,
"Top",
"Only include t tbar as outgoing particles",
16);
static Switch<MEPP2ZJet,unsigned int> interfaceProcess
("Process",
"Which subprocesses to include",
&MEPP2ZJet::_process, 0, false, false);
static SwitchOption interfaceProcessAll
(interfaceProcess,
"All",
"Include all subprocesses",
0);
static SwitchOption interfaceProcessqqbar
(interfaceProcess,
"qqbar",
"Only include q qbar -> Z/gamma g process",
1);
static SwitchOption interfaceProcessqg
(interfaceProcess,
"qg",
"Only include the q g -> Z/gamma q process",
2);
static SwitchOption interfaceProcessqbarg
(interfaceProcess,
"qbarg",
"Only include the qbar g -> Z/gamma qbar process",
3);
static Parameter<MEPP2ZJet,double> interfacePhotonProbablity
("PhotonProbablity",
"Probability for using the \\f$1/s^2\\f$ piece for the"
" generation of the gauge boson mass",
&MEPP2ZJet::_pprob, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Switch<MEPP2ZJet,unsigned int> interfaceGammaZ
("GammaZ",
"Which terms to include",
&MEPP2ZJet::_gammaZ, 0, false, false);
static SwitchOption interfaceGammaZAll
(interfaceGammaZ,
"All",
"Include both gamma and Z terms",
0);
static SwitchOption interfaceGammaZGamma
(interfaceGammaZ,
"Gamma",
"Only include the photon",
1);
static SwitchOption interfaceGammaZZ
(interfaceGammaZ,
"Z",
"Only include the Z",
2);
static Switch<MEPP2ZJet,unsigned int> interfaceWidthOption
("WidthOption",
"The option for handling the width of the off-shell W boson",
&MEPP2ZJet::_widthopt, 1, false, false);
static SwitchOption interfaceWidthOptionFixedDenominator
(interfaceWidthOption,
"FixedDenominator",
"Use a fxied with in the W propagator but the full matrix element"
" in the numerator",
1);
static SwitchOption interfaceWidthOptionAllRunning
(interfaceWidthOption,
"AllRunning",
"Use a running width in the W propagator and the full matrix "
"element in the numerator",
2);
}
void MEPP2ZJet::getDiagrams() const {
// which intermediates to include
bool gamma = _gammaZ==0 || _gammaZ==1;
bool Z0 = _gammaZ==0 || _gammaZ==2;
// pointer for gluon
tcPDPtr g = getParticleData(ParticleID::g);
bool quark,lepton;
for ( int ix=1; ix<17; ++ix ) {
// increment counter to switch between quarks and leptons
if(ix==7) ix+=4;
// is it a valid quark process
quark=ix<=6&&(_zdec==0||_zdec==1||_zdec-10==ix);
// is it a valid lepton process
lepton=ix>=11&&ix<=16&&
(_zdec==0||_zdec==2||(_zdec==3&&ix%2==1)||
(_zdec==4&&ix%2==0)||(ix%2==0&&(ix-10)/2==_zdec-7)||
(ix%2==1&&(ix-9)/2==_zdec-4));
// if not a validf process continue
if(!(quark||lepton)) continue;
// pointer for Z decay products
tcPDPtr lm = getParticleData(ix);
tcPDPtr lp = lm->CC();
for (int i = 1; i <= _maxflavour; ++i ) {
tcPDPtr q = getParticleData(i);
tcPDPtr qb = q->CC();
// q qbar -> Z g -> l+l- g
if(_process==0||_process==1) {
if(gamma) add(new_ptr((Tree2toNDiagram(3), q, q, qb, 1, _gamma,
2, g, 4, lm, 4, lp, -1)));
if(Z0) add(new_ptr((Tree2toNDiagram(3), q, q, qb, 1, _z0,
2, g, 4, lm, 4, lp, -2)));
if(gamma) add(new_ptr((Tree2toNDiagram(3), q, q, qb, 2, _gamma,
1, g, 4, lm, 4, lp, -3)));
if(Z0) add(new_ptr((Tree2toNDiagram(3), q, q, qb, 2, _z0,
1, g, 4, lm, 4, lp, -4)));
}
// q g -> Z q -> l+l- qbar
if(_process==0||_process==2) {
if(gamma) add(new_ptr((Tree2toNDiagram(3), q, q, g, 1, _gamma,
2, q, 4, lm, 4, lp, -5)));
if(Z0) add(new_ptr((Tree2toNDiagram(3), q, q, g, 1, _z0,
2, q, 4, lm, 4, lp, -6)));
if(gamma) add(new_ptr((Tree2toNDiagram(2), q, g, 1, q, 3, _gamma,
3, q, 4, lm, 4, lp, -7)));
if(Z0) add(new_ptr((Tree2toNDiagram(2), q, g, 1, q, 3, _z0,
3, q, 4, lm, 4, lp, -8)));
}
// qbar g -> Z qbar -> l+l- qbar
if(_process==0||_process==3) {
if(gamma) add(new_ptr((Tree2toNDiagram(3), qb, qb, g, 1, _gamma,
2, qb, 4, lm, 4, lp, -9 )));
if(Z0) add(new_ptr((Tree2toNDiagram(3), qb, qb, g, 1, _z0,
2, qb, 4, lm, 4, lp, -10)));
if(gamma) add(new_ptr((Tree2toNDiagram(2), qb, g, 1, qb, 3, _gamma,
3, qb, 4, lm, 4, lp, -11)));
if(Z0) add(new_ptr((Tree2toNDiagram(2), qb, g, 1, qb, 3, _z0,
3, qb, 4, lm, 4, lp, -12)));
}
}
}
}
unsigned int MEPP2ZJet::orderInAlphaS() const {
return 1;
}
unsigned int MEPP2ZJet::orderInAlphaEW() const {
return 2;
}
void MEPP2ZJet::persistentOutput(PersistentOStream & os) const {
os << _theFFZVertex << _theFFPVertex << _theQQGVertex << _z0 << _widthopt
<< _gamma << _process << _maxflavour << _zdec << _pprob << _gammaZ;
}
void MEPP2ZJet::persistentInput(PersistentIStream & is, int) {
is >> _theFFZVertex >> _theFFPVertex >> _theQQGVertex >> _z0 >> _widthopt
>> _gamma >> _process >> _maxflavour >> _zdec >> _pprob >> _gammaZ;
}
int MEPP2ZJet::nDim() const {
return 5;
}
Selector<const ColourLines *>
MEPP2ZJet::colourGeometries(tcDiagPtr diag) const {
// colour lines for q qbar -> Z g
static const ColourLines cqqbar[4]={ColourLines("1 2 5,-3 -5"),
ColourLines("1 5,-5 2 -3"),
ColourLines("1 2 5,-3 -5, 6 -7"),
ColourLines("1 5,-5 2 -3, 6 -7")};
// colour lines for q g -> Z q
static const ColourLines cqg [4]={ColourLines("1 2 -3,3 5"),
ColourLines("1 -2,2 3 5"),
ColourLines("1 2 -3,3 5, 6 -7"),
ColourLines("1 -2,2 3 5, 6 -7")};
// colour lines for qbar q -> Z qbar
static const ColourLines cqbarg[4]={ColourLines("-1 -2 3,-3 -5"),
ColourLines("-1 2,-2 -3 -5"),
ColourLines("-1 -2 3,-3 -5, 6 -7"),
ColourLines("-1 2,-2 -3 -5, 6 -7")};
// select the correct line
unsigned int icol = mePartonData()[3]->coloured() ? 2 : 0;
Selector<const ColourLines *> sel;
switch(abs(diag->id())) {
case 1 : case 2:
sel.insert(1.0, &cqqbar[icol]);
break;
case 3 : case 4:
sel.insert(1.0, &cqqbar[icol+1]);
break;
case 5 : case 6:
sel.insert(1.0, &cqg[icol]);
break;
case 7 : case 8:
sel.insert(1.0, &cqg[icol+1]);
break;
case 9 : case 10:
sel.insert(1.0, &cqbarg[icol]);
break;
case 11 : case 12:
sel.insert(1.0, &cqbarg[icol+1]);
break;
}
return sel;
}
Selector<MEBase::DiagramIndex>
MEPP2ZJet::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
int id=abs(diags[i]->id());
if (id <= 4 ) sel.insert(meInfo()[id-1],i);
else if(id <= 8 ) sel.insert(meInfo()[id-5],i);
else if(id <= 12) sel.insert(meInfo()[id-9],i);
}
return sel;
}
Energy2 MEPP2ZJet::scale() const {
return _scale;
}
CrossSection MEPP2ZJet::dSigHatDR() const {
return me2()*jacobian()/(16.0*sqr(Constants::pi)*sHat())*sqr(hbarc);
}
bool MEPP2ZJet::generateKinematics(const double * r) {
// initialize jacobian
jacobian(1.);
// cms energy
Energy ecm=sqrt(sHat());
// first generate the mass of the off-shell gauge boson
// minimum mass of the
tcPDVector ptemp;
ptemp.push_back(mePartonData()[3]);
ptemp.push_back(mePartonData()[4]);
Energy2 minMass2 = max(lastCuts().minSij(mePartonData()[3],mePartonData()[4]),
lastCuts().minS(ptemp));
// minimum pt of the jet
Energy ptmin = max(lastCuts().minKT(mePartonData()[2]),
lastCuts().minKT(_z0));
// maximum mass of the gauge boson so pt is possible
Energy2 maxMass2 = min(ecm*(ecm-2.*ptmin),lastCuts().maxS(ptemp));
if(maxMass2<=ZERO||minMass2<ZERO) return false;
// also impose the limits from the ParticleData object
minMass2 = max(minMass2,sqr(_z0->massMin()));
maxMass2 = min(maxMass2,sqr(_z0->massMax()));
// also impose the limits from the ParticleData object
if(maxMass2<minMass2) return false;
// generation of the mass
Energy M(_z0->mass()),Gamma(_z0->width());
Energy2 M2(sqr(M)),MG(M*Gamma);
double rhomin = atan2((minMass2-M2),MG);
double rhomax = atan2((maxMass2-M2),MG);
if(r[1]<_pprob) {
double rand=r[1]/_pprob;
_mz2=minMass2*maxMass2/(minMass2+rand*(maxMass2-minMass2));
}
else {
double rand=(r[1]-_pprob)/(1.-_pprob);
_mz2=M2+MG*tan(rhomin+rand*(rhomax-rhomin));
}
Energy mz=sqrt(_mz2);
InvEnergy2 emjac1 = _pprob*minMass2*maxMass2/(maxMass2-minMass2)/sqr(_mz2);
InvEnergy2 emjac2 = (1.-_pprob)*MG/(rhomax-rhomin)/(sqr(_mz2-M2)+sqr(MG));
// jacobian
jacobian(jacobian()/sHat()/(emjac1+emjac2));
// set the masses of the outgoing particles to 2-2 scattering
meMomenta()[2].setMass(ZERO);
Lorentz5Momentum pz(mz);
// generate the polar angle of the hard scattering
double ctmin(-1.0), ctmax(1.0);
Energy q(ZERO);
try {
q = SimplePhaseSpace::getMagnitude(sHat(), meMomenta()[2].mass(),mz);
}
catch ( ImpossibleKinematics ) {
return false;
}
Energy2 pq = sqrt(sHat())*q;
if ( ptmin > ZERO ) {
double ctm = 1.0 - sqr(ptmin/q);
if ( ctm <= 0.0 ) return false;
ctmin = max(ctmin, -sqrt(ctm));
ctmax = min(ctmax, sqrt(ctm));
}
if ( ctmin >= ctmax ) return false;
double cth = getCosTheta(ctmin, ctmax, r[0]);
Energy pt = q*sqrt(1.0-sqr(cth));
double phi = 2.0*Constants::pi*r[2];
meMomenta()[2].setVect(Momentum3( pt*sin(phi), pt*cos(phi), q*cth));
pz.setVect( Momentum3(-pt*sin(phi),-pt*cos(phi),-q*cth));
meMomenta()[2].rescaleEnergy();
pz.rescaleEnergy();
// set the scale
_scale = _mz2+sqr(pt);
// generate the momenta of the Z decay products
meMomenta()[3].setMass(mePartonData()[3]->mass());
meMomenta()[4].setMass(mePartonData()[4]->mass());
Energy q2 = ZERO;
try {
q2 = SimplePhaseSpace::getMagnitude(_mz2, meMomenta()[3].mass(),
meMomenta()[4].mass());
} catch ( ImpossibleKinematics ) {
return false;
}
double cth2 =-1.+2.*r[3];
double phi2=Constants::twopi*r[4];
Energy pt2 =q2*sqrt(1.-sqr(cth2));
Lorentz5Momentum pl[2]={Lorentz5Momentum( pt2*cos(phi2), pt2*sin(phi2), q2*cth2,ZERO,
meMomenta()[3].mass()),
Lorentz5Momentum(-pt2*cos(phi2),-pt2*sin(phi2),-q2*cth2,ZERO,
meMomenta()[4].mass())};
pl[0].rescaleEnergy();
pl[1].rescaleEnergy();
Boost boostv(pz.boostVector());
pl[0].boost(boostv);
pl[1].boost(boostv);
meMomenta()[3] = pl[0];
meMomenta()[4] = pl[1];
// check passes all the cuts
vector<LorentzMomentum> out(3);
tcPDVector tout(3);
for(unsigned int ix=0;ix<3;++ix) {
out[ ix] = meMomenta()[ix+2];
tout[ix] = mePartonData()[ix+2];
}
if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) )
return false;
// jacobian
jacobian((pq/sHat())*Constants::pi*jacobian()/8./sqr(Constants::pi)*q2/mz);
return true;
}
double MEPP2ZJet::me2() const {
InvEnergy2 output(ZERO);
// construct spinors for the leptons (always the same)
vector<SpinorBarWaveFunction> lm;
vector<SpinorWaveFunction> lp;
SpinorBarWaveFunction lmout(meMomenta()[3],mePartonData()[3],outgoing);
SpinorWaveFunction lpout(meMomenta()[4],mePartonData()[4],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
lmout.reset(ix);lm.push_back(lmout);
lpout.reset(ix);lp.push_back(lpout);
}
// q g to q Z
if(mePartonData()[0]->id()<=6&&mePartonData()[0]->id()>0&&
mePartonData()[1]->id()==ParticleID::g) {
// polarization states for the particles
vector<SpinorWaveFunction> fin;
vector<VectorWaveFunction> gin;
vector<SpinorBarWaveFunction> fout;
SpinorWaveFunction qin (meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction glin(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction qout(meMomenta()[2],mePartonData()[2],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
qin.reset(ix) ; fin.push_back(qin);
glin.reset(2*ix); gin.push_back(glin);
qout.reset(ix);fout.push_back(qout);
}
output=qgME(fin,gin,fout,lm,lp);
}
// qbar g to qbar Z
else if(mePartonData()[0]->id()>=-6&&mePartonData()[0]->id()<0&&
mePartonData()[1]->id()==ParticleID::g) {
vector<SpinorBarWaveFunction> ain;
vector<VectorWaveFunction> gin;
vector<SpinorWaveFunction> aout;
SpinorBarWaveFunction qbin (meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction glin (meMomenta()[1],mePartonData()[1],incoming);
SpinorWaveFunction qbout(meMomenta()[2],mePartonData()[2],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
qbin .reset(ix ); ain .push_back(qbin );
glin .reset(2*ix); gin .push_back(glin );
qbout.reset(ix ); aout.push_back(qbout);
}
output=qbargME(ain,gin,aout,lm,lp);
}
// q qbar to g Z
else {
vector<SpinorWaveFunction> fin;
vector<SpinorBarWaveFunction> ain;
vector<VectorWaveFunction> gout;
SpinorWaveFunction qin (meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction qbin(meMomenta()[1],mePartonData()[1],incoming);
VectorWaveFunction glout(meMomenta()[2],mePartonData()[2],outgoing);
for(unsigned int ix=0;ix<2;++ix) {
qin.reset(ix) ; fin.push_back(qin);
qbin.reset(ix) ; ain.push_back(qbin);
glout.reset(2*ix); gout.push_back(glout);
}
output=qqbarME(fin,ain,gout,lm,lp);
}
return output*sHat();
}
InvEnergy2 MEPP2ZJet::qqbarME(vector<SpinorWaveFunction> & fin,
vector<SpinorBarWaveFunction> & ain,
vector<VectorWaveFunction> & gout,
vector<SpinorBarWaveFunction> & lm,
vector<SpinorWaveFunction> & lp,
bool calc) const {
// if calculation spin corrections construct the me
if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1,PDT::Spin1Half,
PDT::Spin1Half));
// diagrams to include
bool gamma = _gammaZ==0 || _gammaZ==1;
bool Z0 = _gammaZ==0 || _gammaZ==2;
// some integers
unsigned int ihel1,ihel2,ohel1,ohel2,ohel3;
// compute the leptonic photon and Z currents for speed
VectorWaveFunction bcurr[2][2][2];
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
// photon current
if(gamma) bcurr[0][ohel2][ohel3]=
_theFFPVertex->evaluate(_mz2,1,_gamma,lp[ohel3],lm[ohel2]);
// Z current
if(Z0) bcurr[1][ohel2][ohel3]=
_theFFZVertex->evaluate(_mz2,_widthopt,_z0,lp[ohel3],lm[ohel2]);
}
}
// compute the matrix elements
double me[5]={0.,0.,0.,0.,0.};
Complex diag[4];
SpinorWaveFunction inters;
SpinorBarWaveFunction interb;
for(ihel1=0;ihel1<2;++ihel1) {
for(ihel2=0;ihel2<2;++ihel2) {
for(ohel1=0;ohel1<2;++ohel1) {
// intermediates for the diagrams
inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[0],
fin[ihel1],gout[ohel1]);
interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[1],
ain[ihel2],gout[ohel1]);
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
diag[0] = gamma ?
_theFFPVertex->evaluate(_mz2,fin[ihel1],interb,
bcurr[0][ohel2][ohel3]) : 0.;
diag[1]= Z0 ?
_theFFZVertex->evaluate(_mz2,fin[ihel1],interb,
bcurr[1][ohel2][ohel3]) : 0.;
diag[2]= gamma ?
_theFFPVertex->evaluate(_mz2,inters,ain[ihel2],
bcurr[0][ohel2][ohel3]) : 0.;
diag[3]= Z0 ?
_theFFZVertex->evaluate(_mz2,inters,ain[ihel2],
bcurr[1][ohel2][ohel3]) : 0.;
// diagram contributions
me[1] += norm(diag[0]);
me[2] += norm(diag[1]);
me[3] += norm(diag[2]);
me[4] += norm(diag[3]);
// total
diag[0] += diag[1] + diag[2] + diag[3];
me[0] += norm(diag[0]);
if(calc) _me(ihel1,ihel2,2*ohel1,ohel2,ohel3) = diag[0];
}
}
}
}
}
// results
// initial state spin and colour average
double colspin = 1./9./4.;
// and C_F N_c from matrix element
colspin *= 4.;
// and for Z decay products
if(mePartonData()[3]->coloured()) colspin *= 3.;
DVector save;
for(unsigned int ix=0;ix<5;++ix) {
me[ix] *= colspin;
if(ix>0) save.push_back(me[ix]);
}
meInfo(save);
return me[0]*UnitRemoval::InvE2;
}
InvEnergy2 MEPP2ZJet::qgME(vector<SpinorWaveFunction> & fin,
vector<VectorWaveFunction> & gin,
vector<SpinorBarWaveFunction> & fout,
vector<SpinorBarWaveFunction> & lm,
vector<SpinorWaveFunction> & lp,
bool calc) const {
// diagrams to include
bool gamma = _gammaZ==0 || _gammaZ==1;
bool Z0 = _gammaZ==0 || _gammaZ==2;
// if calculation spin corrections construct the me
if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half));
// some integers
unsigned int ihel1,ihel2,ohel1,ohel2,ohel3;
// compute the leptonic photon and Z currents for speed
VectorWaveFunction bcurr[2][2][2];
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
// photon current
if(gamma) bcurr[0][ohel2][ohel3]=
_theFFPVertex->evaluate(_mz2,1,_gamma,lp[ohel3],lm[ohel2]);
// Z current
if(Z0) bcurr[1][ohel2][ohel3]=
_theFFZVertex->evaluate(_mz2,_widthopt,_z0,lp[ohel3],lm[ohel2]);
}
}
// compute the matrix elements
double me[5]={0.,0.,0.,0.,0.};
Complex diag[4];
SpinorWaveFunction inters;
SpinorBarWaveFunction interb;
Energy2 _scale=scale();
for(ihel1=0;ihel1<2;++ihel1) {
for(ihel2=0;ihel2<2;++ihel2) {
for(ohel1=0;ohel1<2;++ohel1) {
// intermediates for the diagrams
- interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[2],
+ interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[2]->CC(),
fout[ohel1],gin[ihel2]);
inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[0],
fin[ihel1],gin[ihel2]);
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
diag[0]=gamma ?
_theFFPVertex->evaluate(_mz2,fin[ihel1],interb,
bcurr[0][ohel2][ohel3]) : 0.;
diag[1]=Z0 ?
_theFFZVertex->evaluate(_mz2,fin[ihel1],interb,
bcurr[1][ohel2][ohel3]) : 0.;
diag[2]=gamma ?
_theFFPVertex->evaluate(_mz2,inters,fout[ohel1],
bcurr[0][ohel2][ohel3]) : 0.;
diag[3]=Z0 ?
_theFFZVertex->evaluate(_mz2,inters,fout[ohel1],
bcurr[1][ohel2][ohel3]) : 0.;
// diagram contributions
me[1] += norm(diag[0]);
me[2] += norm(diag[1]);
me[3] += norm(diag[2]);
me[4] += norm(diag[3]);
// total
diag[0] += diag[1] + diag[2] + diag[3];
me[0] += norm(diag[0]);
if(calc) _me(ihel1,2*ihel2,ohel1,ohel2,ohel3) = diag[0];
}
}
}
}
}
// results
// initial state spin and colour average
double colspin = 1./24./4.;
// and C_F N_c from matrix element
colspin *= 4.;
// and for Z decay products
if(mePartonData()[3]->coloured()) colspin *= 3.;
DVector save;
for(unsigned int ix=0;ix<5;++ix) {
me[ix] *= colspin;
if(ix>0) save.push_back(me[ix]);
}
meInfo(save);
return me[0]*UnitRemoval::InvE2;
}
InvEnergy2 MEPP2ZJet::qbargME(vector<SpinorBarWaveFunction> & fin,
vector<VectorWaveFunction> & gin,
vector<SpinorWaveFunction> & fout,
vector<SpinorBarWaveFunction> & lm,
vector<SpinorWaveFunction> & lp,
bool calc) const {
// diagrams to include
bool gamma = _gammaZ==0 || _gammaZ==1;
bool Z0 = _gammaZ==0 || _gammaZ==2;
// if calculation spin corrections construct the me
if(calc) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half));
// some integers
unsigned int ihel1,ihel2,ohel1,ohel2,ohel3;
// compute the leptonic photon and Z currents for speed
VectorWaveFunction bcurr[2][2][2];
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
// photon current
if(gamma) bcurr[0][ohel2][ohel3]=
_theFFPVertex->evaluate(_mz2,1,_gamma,lp[ohel3],lm[ohel2]);
// Z current
if(Z0) bcurr[1][ohel2][ohel3]=
_theFFZVertex->evaluate(_mz2,_widthopt,_z0,lp[ohel3],lm[ohel2]);
}
}
// compute the matrix elements
double me[5]={0.,0.,0.,0.,0.};
Complex diag[4];
SpinorWaveFunction inters;
SpinorBarWaveFunction interb;
Energy2 _scale=scale();
for(ihel1=0;ihel1<2;++ihel1) {
for(ihel2=0;ihel2<2;++ihel2) {
for(ohel1=0;ohel1<2;++ohel1) {
// intermediates for the diagrams
- inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[2],
+ inters=_theQQGVertex->evaluate(_scale,5,mePartonData()[2]->CC(),
fout[ohel1],gin[ihel2]);
interb=_theQQGVertex->evaluate(_scale,5,mePartonData()[0],
fin[ihel1],gin[ihel2]);
for(ohel2=0;ohel2<2;++ohel2) {
for(ohel3=0;ohel3<2;++ohel3) {
diag[0]= gamma ?
_theFFPVertex->evaluate(_mz2,inters,fin[ihel1],
bcurr[0][ohel2][ohel3]) : 0.;
diag[1]= Z0 ?
_theFFZVertex->evaluate(_mz2,inters,fin[ihel1],
bcurr[1][ohel2][ohel3]) : 0.;
diag[2]= gamma ?
_theFFPVertex->evaluate(_mz2,fout[ohel1],interb,
bcurr[0][ohel2][ohel3]) : 0.;
diag[3]= Z0 ?
_theFFZVertex->evaluate(_mz2,fout[ohel1],interb,
bcurr[1][ohel2][ohel3]) : 0.;
// diagram contributions
me[1] += norm(diag[0]);
me[2] += norm(diag[1]);
me[3] += norm(diag[2]);
me[4] += norm(diag[3]);
// total
diag[0] += diag[1] + diag[2] + diag[3];
me[0] += norm(diag[0]);
if(calc) _me(ihel1,2*ihel2,ohel1,ohel2,ohel3) = diag[0];
}
}
}
}
}
// results
// initial state spin and colour average
double colspin = 1./24./4.;
// and C_F N_c from matrix element
colspin *= 4.;
// and for Z decay products
if(mePartonData()[3]->coloured()) colspin*=3.;
DVector save;
for(unsigned int ix=0;ix<5;++ix) {
me[ix] *= colspin;
if(ix>0) save.push_back(me[ix]);
}
meInfo(save);
return me[0]*UnitRemoval::InvE2;
}
void MEPP2ZJet::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard(5);
// incoming
hard[0]=sub->incoming().first;
hard[1]=sub->incoming().second;
if((hard[0]->id()<0&&hard[1]->id()<=6)||
hard[0]->id()==ParticleID::g) swap(hard[0],hard[1]);
// outgoing
for(unsigned int ix=0;ix<3;++ix) {
unsigned int iloc;
PPtr mother=sub->outgoing()[ix]->parents()[0];
if(mother&&(mother->id()==ParticleID::gamma||mother->id()==ParticleID::Z0)) {
if(sub->outgoing()[ix]->id()>0) iloc=3;
else iloc=4;
}
else iloc=2;
hard[iloc]=sub->outgoing()[ix];
}
// wavefunctions for the Z decay products
vector<SpinorBarWaveFunction> lm;
vector<SpinorWaveFunction> lp;
SpinorBarWaveFunction(lm,hard[3],outgoing,true,true);
SpinorWaveFunction (lp,hard[4],outgoing,true,true);
// identify hard process and calculate matrix element
// q g to q Z
if(hard[0]->id()<=6&&hard[0]->id()>0&&hard[1]->id()==ParticleID::g) {
vector<SpinorWaveFunction> fin;
vector<VectorWaveFunction> gin;
vector<SpinorBarWaveFunction> fout;
SpinorWaveFunction (fin ,hard[0],incoming,false,true);
VectorWaveFunction (gin ,hard[1],incoming,false,true,true);
SpinorBarWaveFunction (fout,hard[2],outgoing,true ,true);
gin[1]=gin[2];
qgME(fin,gin,fout,lm,lp,true);
}
// qbar g to qbar Z
else if(hard[0]->id()>=-6&&hard[0]->id()<0&&hard[1]->id()==ParticleID::g) {
vector<SpinorBarWaveFunction> ain;
vector<VectorWaveFunction> gin;
vector<SpinorWaveFunction> aout;
SpinorBarWaveFunction(ain ,hard[0],incoming,false,true);
VectorWaveFunction (gin ,hard[1],incoming,false,true,true);
SpinorWaveFunction (aout,hard[2],outgoing,true ,true);
gin[1]=gin[2];
qbargME(ain,gin,aout,lm,lp,true);
}
// q qbar to g Z
else {
vector<SpinorWaveFunction> fin;
vector<SpinorBarWaveFunction> ain;
vector<VectorWaveFunction> gout;
SpinorWaveFunction (fin ,hard[0],incoming,false,true);
SpinorBarWaveFunction(ain ,hard[1],incoming,false,true);
VectorWaveFunction (gout,hard[2],outgoing,true ,true,true);
gout[1]=gout[2];
qqbarME(fin,ain,gout,lm,lp,true);
}
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(_me);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<5;++ix)
(hard[ix]->spinInfo())->productionVertex(hardvertex);
}
diff --git a/MatrixElement/Hadron/MEQCD2to2.cc b/MatrixElement/Hadron/MEQCD2to2.cc
--- a/MatrixElement/Hadron/MEQCD2to2.cc
+++ b/MatrixElement/Hadron/MEQCD2to2.cc
@@ -1,1182 +1,1182 @@
// -*- C++ -*-
//
// MEQCD2to2.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEQCD2to2 class.
//
#include "MEQCD2to2.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Cuts/Cuts.h"
#include "Herwig/MatrixElement/HardVertex.h"
using namespace Herwig;MEQCD2to2::MEQCD2to2():_maxflavour(5),_process(0) {
massOption(vector<unsigned int>(2,0));
}
void MEQCD2to2::rebind(const TranslationMap & trans)
{
_ggggvertex = trans.translate(_ggggvertex);
_gggvertex = trans.translate( _gggvertex);
_qqgvertex = trans.translate( _qqgvertex);
_gluon = trans.translate( _gluon);
for(unsigned int ix=0;ix<_quark.size();++ix)
{_quark[ix]=trans.translate(_quark[ix]);}
for(unsigned int ix=0;ix<_antiquark.size();++ix)
{_antiquark[ix]=trans.translate(_quark[ix]);}
HwMEBase::rebind(trans);
}
IVector MEQCD2to2::getReferences() {
IVector ret = HwMEBase::getReferences();
ret.push_back(_ggggvertex);
ret.push_back(_gggvertex);
ret.push_back(_qqgvertex);
ret.push_back(_gluon);
for(unsigned int ix=0;ix<_quark.size();++ix)
{ret.push_back(_quark[ix]);}
for(unsigned int ix=0;ix<_antiquark.size();++ix)
{ret.push_back(_antiquark[ix]);}
return ret;
}
void MEQCD2to2::doinit() {
// call the base class
HwMEBase::doinit();
// get the vedrtex pointers from the SM object
tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
// do the initialisation
if(hwsm) {
_qqgvertex = hwsm->vertexFFG();
_gggvertex = hwsm->vertexGGG();
_ggggvertex = hwsm->vertexGGGG();
}
else throw InitException() << "Wrong type of StandardModel object in "
<< "MEQCD2to2::doinit() the Herwig version must be used"
<< Exception::runerror;
// get the particle data objects
_gluon=getParticleData(ParticleID::g);
for(int ix=1;ix<=int(_maxflavour);++ix) {
_quark.push_back( getParticleData( ix));
_antiquark.push_back(getParticleData(-ix));
}
}
Energy2 MEQCD2to2::scale() const {
Energy2 s(sHat()),u(uHat()),t(tHat());
return 2.*s*t*u/(s*s+t*t+u*u);
}
void MEQCD2to2::persistentOutput(PersistentOStream & os) const {
os << _ggggvertex << _gggvertex << _qqgvertex << _maxflavour
<< _process << _gluon << _quark << _antiquark;
}
void MEQCD2to2::persistentInput(PersistentIStream & is, int) {
is >> _ggggvertex >> _gggvertex >> _qqgvertex >> _maxflavour
>> _process >> _gluon >> _quark >> _antiquark;
}
unsigned int MEQCD2to2::orderInAlphaS() const {
return 2;
}
unsigned int MEQCD2to2::orderInAlphaEW() const {
return 0;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEQCD2to2,HwMEBase>
describeHerwigMEQCD2to2("Herwig::MEQCD2to2", "HwMEHadron.so");
void MEQCD2to2::Init() {
static ClassDocumentation<MEQCD2to2> documentation
("The MEQCD2to2 class implements the QCD 2->2 processes in hadron-hadron"
" collisions");
static Parameter<MEQCD2to2,unsigned int> interfaceMaximumFlavour
("MaximumFlavour",
"The maximum flavour of the quarks in the process",
&MEQCD2to2::_maxflavour, 5, 1, 5,
false, false, Interface::limited);
static Switch<MEQCD2to2,unsigned int> interfaceProcess
("Process",
"Which subprocesses to include",
&MEQCD2to2::_process, 0, false, false);
static SwitchOption interfaceProcessAll
(interfaceProcess,
"All",
"Include all subprocesses",
0);
static SwitchOption interfaceProcess1
(interfaceProcess,
"gg2gg",
"Include only gg->gg subprocesses",
1);
static SwitchOption interfaceProcess2
(interfaceProcess,
"gg2qqbar",
"Include only gg -> q qbar processes",
2);
static SwitchOption interfaceProcessqqbargg
(interfaceProcess,
"qqbar2gg",
"Include only q qbar -> gg processes",
3);
static SwitchOption interfaceProcessqgqg
(interfaceProcess,
"qg2qg",
"Include only q g -> q g processes",
4);
static SwitchOption interfaceProcessqbargqbarg
(interfaceProcess,
"qbarg2qbarg",
"Include only qbar g -> qbar g processes",
5);
static SwitchOption interfaceProcessqqqq
(interfaceProcess,
"qq2qq",
"Include only q q -> q q processes",
6);
static SwitchOption interfaceProcessqbarqbarqbarqbar
(interfaceProcess,
"qbarqbar2qbarqbar",
"Include only qbar qbar -> qbar qbar processes",
7);
static SwitchOption interfaceProcessqqbarqqbar
(interfaceProcess,
"qqbar2qqbar",
"Include only q qbar -> q qbar processes",
8);
}
Selector<MEBase::DiagramIndex>
MEQCD2to2::diagrams(const DiagramVector & diags) const {
// select the diagram, this is easy for us as we have already done it
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
if(diags[i]->id()==-int(_diagram)) sel.insert(1.0, i);
else sel.insert(0., i);
}
return sel;
}
double MEQCD2to2::gg2qqbarME(vector<VectorWaveFunction> &g1,
vector<VectorWaveFunction> &g2,
vector<SpinorBarWaveFunction> & q,
vector<SpinorWaveFunction> & qbar,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1Half));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorWaveFunction inters;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
interv=_gggvertex->evaluate(mt,5,_gluon,g1[ihel1],g2[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
//first t-channel diagram
- inters =_qqgvertex->evaluate(mt,5,qbar[ohel2].particle(),
+ inters =_qqgvertex->evaluate(mt,5,qbar[ohel2].particle()->CC(),
qbar[ohel2],g2[ihel2]);
diag[0]=_qqgvertex->evaluate(mt,inters,q[ohel1],g1[ihel1]);
//second t-channel diagram
- inters =_qqgvertex->evaluate(mt,5,qbar[ohel2].particle(),
+ inters =_qqgvertex->evaluate(mt,5,qbar[ohel2].particle()->CC(),
qbar[ohel2],g1[ihel1]);
diag[1]=_qqgvertex->evaluate(mt,inters,q[ohel1],g2[ihel2]);
// s-channel diagram
diag[2]=_qqgvertex->evaluate(mt,qbar[ohel2],q[ohel1],interv);
// colour flows
flow[0]=diag[0]+diag[2];
flow[1]=diag[1]-diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(2*ihel1,2*ihel2,ohel1,ohel2)=flow[iflow-1];
}
}
}
}
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//cerr << "testing matrix element "
// << 48.*(1./6./u/t-3./8./s/s)*(t*t+u*u)*sqr(alphas)/output << endl;
// select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=4+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return output/48.;
}
double MEQCD2to2::qqbar2ggME(vector<SpinorWaveFunction> & q,
vector<SpinorBarWaveFunction> & qbar,
vector<VectorWaveFunction> &g1,
vector<VectorWaveFunction> &g2,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1,PDT::Spin1));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorWaveFunction inters;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
interv=_qqgvertex->evaluate(mt,5,_gluon,q[ihel1],qbar[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first t-channel diagram
inters=_qqgvertex->evaluate(mt,5,q[ihel1].particle()->CC(),
q[ihel1],g1[ohel1]);
diag[0]=_qqgvertex->evaluate(mt,inters,qbar[ihel2],g2[ohel2]);
// second t-channel diagram
inters=_qqgvertex->evaluate(mt,5,q[ihel1].particle()->CC(),
q[ihel1],g2[ohel2]);
diag[1]=_qqgvertex->evaluate(mt,inters,qbar[ihel2],g1[ohel1]);
// s-channel diagram
diag[2]=_gggvertex->evaluate(mt,g1[ohel1],g2[ohel2],interv);
// colour flows
flow[0]=diag[0]-diag[2];
flow[1]=diag[1]+diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,ihel2,2*ohel1,2*ohel2)=flow[iflow-1];
}
}
}
}
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//cerr << "testing matrix element "
// << 27./2.*0.5*(32./27./u/t-8./3./s/s)*(t*t+u*u)*sqr(alphas)/output << endl;
//select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=7+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return 2.*output/27.;
}
double MEQCD2to2::qg2qgME(vector<SpinorWaveFunction> & qin,
vector<VectorWaveFunction> &g2,
vector<SpinorBarWaveFunction> & qout,
vector<VectorWaveFunction> &g4,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorWaveFunction inters,inters2;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
inters=_qqgvertex->evaluate(mt,5,qin[ihel1].particle()->CC(),
qin[ihel1],g2[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// s-channel diagram
diag[0]=_qqgvertex->evaluate(mt,inters,qout[ohel1],g4[ohel2]);
// first t-channel
inters2=_qqgvertex->evaluate(mt,5,qin[ihel1].particle()->CC(),
qin[ihel1],g4[ohel2]);
diag[1]=_qqgvertex->evaluate(mt,inters2,qout[ohel1],g2[ihel2]);
// second t-channel
interv=_qqgvertex->evaluate(mt,5,_gluon,qin[ihel1],qout[ohel1]);
diag[2]=_gggvertex->evaluate(mt,g2[ihel2],g4[ohel2],interv);
// colour flows
flow[0]=diag[0]-diag[2];
flow[1]=diag[1]+diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,2*ihel2,ohel1,2*ohel2)=flow[iflow-1];
}
}
}
}
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//cerr << "testing matrix element "
// << 18./output*(-4./9./s/u+1./t/t)*(s*s+u*u)*sqr(alphas) << endl;
//select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=10+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return output/18.;
}
double MEQCD2to2::gg2ggME(vector<VectorWaveFunction> &g1,vector<VectorWaveFunction> &g2,
vector<VectorWaveFunction> &g3,vector<VectorWaveFunction> &g4,
unsigned int iflow) const {
// colour factors for different flows
static const double c1 = 4.*( sqr(9.)/8.-3.*9./8.+1.-0.75/9.);
static const double c2 = 4.*(-0.25*9. +1.-0.75/9.);
// scale
Energy2 mt(scale());
// // matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1,
PDT::Spin1,PDT::Spin1));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[3]={0.,0.,0.};
Complex diag[3],flow[3];
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// s-channel diagram
diag[0]=_ggggvertex->evaluate(mt,1,g3[ohel1],g1[ihel1],
g4[ohel2],g2[ihel2]);
// t-channel
diag[1]=_ggggvertex->evaluate(mt,1,g1[ihel1],g2[ihel2],
g3[ohel1],g4[ohel2]);
// u-channel
diag[2]=_ggggvertex->evaluate(mt,1,g2[ihel2],g1[ihel1],
g3[ohel1],g4[ohel2]);
// colour flows
flow[0] = diag[0]-diag[2];
flow[1] = -diag[0]-diag[1];
flow[2] = diag[1]+diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) {
sumdiag[ix] += norm(diag[ix]);
sumflow[ix] += norm(flow[ix]);
}
// total
output += c1*(norm(flow[0])+norm(flow[1])+norm(flow[2]))
+2.*c2*real(flow[0]*conj(flow[1])+flow[0]*conj(flow[2])+
flow[1]*conj(flow[2]));
// store the me if needed
if(iflow!=0) _me(2*ihel1,2*ihel2,2*ohel1,2*ohel2)=flow[iflow-1];
}
}
}
}
// spin, colour and identical particle factorsxs
output /= 4.*64.*2.;
// test code vs me from ESW
// Energy2 u(uHat()),t(tHat()),s(sHat());
// using Constants::pi;
// double alphas(4.*pi*SM().alphaS(mt));
// cerr << "testing matrix element "
// << 1./output*9./4.*(3.-t*u/s/s-s*u/t/t-s*t/u/u)*sqr(alphas) << endl;
// select a colour flow
_flow=1+UseRandom::rnd3(sumflow[0],sumflow[1],sumflow[2]);
// and diagram
if(_flow==1) _diagram = 1+2*UseRandom::rnd2(sumdiag[0],sumdiag[2]);
else if(_flow==2) _diagram = 1+ UseRandom::rnd2(sumdiag[0],sumdiag[1]);
else if(_flow==3) _diagram = 2+ UseRandom::rnd2(sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return output;
}
double MEQCD2to2::qbarg2qbargME(vector<SpinorBarWaveFunction> & qin,
vector<VectorWaveFunction> &g2,
vector<SpinorWaveFunction> & qout,
vector<VectorWaveFunction> &g4,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorBarWaveFunction inters,inters2;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
inters=_qqgvertex->evaluate(mt,5,qin[ihel1].particle()->CC(),
qin[ihel1],g2[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// s-channel diagram
diag[0]=_qqgvertex->evaluate(mt,qout[ohel1],inters,g4[ohel2]);
// first t-channel
inters2=_qqgvertex->evaluate(mt,5,qin[ihel1].particle()->CC(),
qin[ihel1],g4[ohel2]);
diag[1]=_qqgvertex->evaluate(mt,qout[ohel1],inters2,g2[ihel2]);
// second t-channel
interv=_qqgvertex->evaluate(mt,5,_gluon,qout[ohel1],qin[ihel1]);
diag[2]=_gggvertex->evaluate(mt,g2[ihel2],g4[ohel2],interv);
// colour flows
flow[0]=diag[0]+diag[2];
flow[1]=diag[1]-diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,2*ihel2,ohel1,2*ohel2)=flow[iflow-1];
}
}
}
}
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//cerr << "testing matrix element "
// << 18./output*(-4./9./s/u+1./t/t)*(s*s+u*u)*sqr(alphas) << endl;
//select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=13+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return output/18.;
}
double MEQCD2to2::qq2qqME(vector<SpinorWaveFunction> & q1,
vector<SpinorWaveFunction> & q2,
vector<SpinorBarWaveFunction> & q3,
vector<SpinorBarWaveFunction> & q4,
unsigned int iflow) const {
// identify special case of identical quarks
bool identical(q1[0].id()==q2[0].id());
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));
// calculate the matrix element
double output(0.),sumdiag[2]={0.,0.};
Complex diag[2];
VectorWaveFunction interv;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q3[ohel1]);
diag[0] = _qqgvertex->evaluate(mt,q2[ihel2],q4[ohel2],interv);
// second diagram if identical
if(identical) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q4[ohel2]);
diag[1]=_qqgvertex->evaluate(mt,q2[ihel2],q3[ohel1],interv);
}
else diag[1]=0.;
// sum of diagrams
for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]);
// total
output +=real(diag[0]*conj(diag[0])+diag[1]*conj(diag[1])
+2./3.*diag[0]*conj(diag[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,ihel2,ohel1,ohel2)=diag[iflow-1];
}
}
}
}
// identical particle symmetry factor if needed
if(identical) output*=0.5;
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//if(identical)
// {cerr << "testing matrix element A "
// << 18./output*0.5*(4./9.*((s*s+u*u)/t/t+(s*s+t*t)/u/u)
// -8./27.*s*s/u/t)*sqr(alphas) << endl;}
//else
// {cerr << "testing matrix element B "
// << 18./output*(4./9.*(s*s+u*u)/t/t)*sqr(alphas) << endl;}
//select a colour flow
_flow=1+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=16+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// final part of colour and spin factors
return output/18.;
}
double MEQCD2to2::qbarqbar2qbarqbarME(vector<SpinorBarWaveFunction> & q1,
vector<SpinorBarWaveFunction> & q2,
vector<SpinorWaveFunction> & q3,
vector<SpinorWaveFunction> & q4,
unsigned int iflow) const {
// identify special case of identical quarks
bool identical(q1[0].id()==q2[0].id());
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0)
{_me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));}
// calculate the matrix element
double output(0.),sumdiag[2]={0.,0.};
Complex diag[2];
VectorWaveFunction interv;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
interv = _qqgvertex->evaluate(mt,5,_gluon,q3[ohel1],q1[ihel1]);
diag[0] = _qqgvertex->evaluate(mt,q4[ohel2],q2[ihel2],interv);
// second diagram if identical
if(identical) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q4[ohel2],q1[ihel1]);
diag[1]=_qqgvertex->evaluate(mt,q3[ohel1],q2[ihel2],interv);
}
else diag[1]=0.;
// sum of diagrams
for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]);
// total
output +=real(diag[0]*conj(diag[0])+diag[1]*conj(diag[1])
+2./3.*diag[0]*conj(diag[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,ihel2,ohel1,ohel2)=diag[iflow-1];
}
}
}
}
// identical particle symmetry factor if needed
if(identical){output*=0.5;}
// test code vs me from ESW
// Energy2 u(uHat()),t(tHat()),s(sHat());
// double alphas(4.*pi*SM().alphaS(mt));
// if(identical)
// {cerr << "testing matrix element A "
// << 18./output*0.5*(4./9.*((s*s+u*u)/t/t+(s*s+t*t)/u/u)
// -8./27.*s*s/u/t)*sqr(alphas) << endl;}
// else
// {cerr << "testing matrix element B "
// << 18./output*(4./9.*(s*s+u*u)/t/t)*sqr(alphas) << endl;}
//select a colour flow
_flow=1+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=18+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// final part of colour and spin factors
return output/18.;
}
double MEQCD2to2::qqbar2qqbarME(vector<SpinorWaveFunction> & q1,
vector<SpinorBarWaveFunction> & q2,
vector<SpinorBarWaveFunction> & q3,
vector<SpinorWaveFunction> & q4,
unsigned int iflow) const {
// type of process
bool diagon[2]={q1[0].id()== -q2[0].id(),
q1[0].id()== -q3[0].id()};
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));
// calculate the matrix element
double output(0.),sumdiag[2]={0.,0.};
Complex diag[2];
VectorWaveFunction interv;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
if(diagon[0]) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q2[ihel2]);
diag[0] = _qqgvertex->evaluate(mt,q4[ohel2],q3[ohel1],interv);
}
else diag[0]=0.;
// second diagram
if(diagon[1]) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q3[ohel1]);
diag[1]=_qqgvertex->evaluate(mt,q4[ohel2],q2[ihel2],interv);
}
else diag[1]=0.;
// sum of diagrams
for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]);
// total
output +=real(diag[0]*conj(diag[0])+diag[1]*conj(diag[1])
+2./3.*diag[0]*conj(diag[1]));
// store the me if needed
if(iflow!=0){_me(ihel1,ihel2,ohel1,ohel2)=diag[iflow-1];}
}
}
}
}
// test code vs me from ESW
// Energy2 u(uHat()),t(tHat()),s(sHat());
// double alphas(4.*pi*SM().alphaS(mt));
// if(diagon[0]&&diagon[1]) {
// cerr << "testing matrix element A "
// << q1[0].id() << " " << q2[0].id() << " -> "
// << q3[0].id() << " " << q4[0].id() << " "
// << 18./output*0.5*(4./9.*((s*s+u*u)/t/t+(u*u+t*t)/s/s)
// -8./27.*u*u/s/t)*sqr(alphas) << endl;
// }
// else if(diagon[0]) {
// cerr << "testing matrix element B "
// << q1[0].id() << " " << q2[0].id() << " -> "
// << q3[0].id() << " " << q4[0].id() << " "
// << 18./output*(4./9.*(t*t+u*u)/s/s)*sqr(alphas) << endl;
// }
// else if(diagon[1]) {
// cerr << "testing matrix element C "
// << q1[0].id() << " " << q2[0].id() << " -> "
// << q3[0].id() << " " << q4[0].id() << " "
// << 18./output*(4./9.*(s*s+u*u)/t/t)*sqr(alphas) << endl;
// }
//select a colour flow
_flow=1+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=20+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// final part of colour and spin factors
return output/18.;
}
void MEQCD2to2::getDiagrams() const {
// gg-> gg subprocess
if(_process==0||_process==1) {
// s-channel
add(new_ptr((Tree2toNDiagram(2),_gluon,_gluon, 1, _gluon,
3,_gluon, 3, _gluon, -1)));
// first t-channel
add(new_ptr((Tree2toNDiagram(3),_gluon,_gluon,_gluon,
1,_gluon, 2,_gluon,-2)));
// second t-channel
add(new_ptr((Tree2toNDiagram(3),_gluon,_gluon,_gluon,
2,_gluon, 1,_gluon,-3)));
}
// processes involving one quark line
for(unsigned int ix=0;ix<_maxflavour;++ix) {
// gg -> q qbar subprocesses
if(_process==0||_process==2) {
// first t-channel
add(new_ptr((Tree2toNDiagram(3),_gluon,_antiquark[ix],_gluon,
1,_quark[ix], 2,_antiquark[ix],-4)));
// interchange
add(new_ptr((Tree2toNDiagram(3),_gluon,_antiquark[ix],_gluon,
2,_quark[ix], 1,_antiquark[ix],-5)));
// s-channel
add(new_ptr((Tree2toNDiagram(2),_gluon,_gluon, 1, _gluon,
3,_quark[ix], 3, _antiquark[ix], -6)));
}
// q qbar -> g g subprocesses
if(_process==0||_process==3) {
// first t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_antiquark[ix],_antiquark[ix],
1,_gluon, 2,_gluon,-7)));
// second t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_antiquark[ix],_antiquark[ix],
2,_gluon, 1,_gluon,-8)));
// s-channel
add(new_ptr((Tree2toNDiagram(2),_quark[ix], _antiquark[ix],
1, _gluon, 3, _gluon, 3, _gluon,-9)));
}
// q g -> q g subprocesses
if(_process==0||_process==4) {
// s-channel
add(new_ptr((Tree2toNDiagram(2),_quark[ix], _gluon,
1, _quark[ix], 3, _quark[ix], 3, _gluon,-10)));
// quark t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_quark[ix],_gluon,
2,_quark[ix],1,_gluon,-11)));
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_gluon,_gluon,
1,_quark[ix],2,_gluon,-12)));
}
// qbar g -> qbar g subprocesses
if(_process==0||_process==5) {
// s-channel
add(new_ptr((Tree2toNDiagram(2),_antiquark[ix], _gluon,
1, _antiquark[ix], 3, _antiquark[ix], 3, _gluon,-13)));
// quark t-channel
add(new_ptr((Tree2toNDiagram(3),_antiquark[ix],_antiquark[ix],_gluon,
2,_antiquark[ix],1,_gluon,-14)));
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_antiquark[ix],_gluon,_gluon,
1,_antiquark[ix],2,_gluon,-15)));
}
// processes involving two quark lines
for(unsigned int iy=ix;iy<_maxflavour;++iy) {
// q q -> q q subprocesses
if(_process==0||_process==6) {
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_gluon,_quark[iy],
1,_quark[ix],2,_quark[iy],-16)));
// exchange for identical quarks
if(ix==iy)
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_gluon,_quark[iy],
2,_quark[ix],1,_quark[iy],-17)));
}
// qbar qbar -> qbar qbar subprocesses
if(_process==0||_process==7) {
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_antiquark[ix],_gluon,_antiquark[iy],
1,_antiquark[ix],2,_antiquark[iy],-18)));
// exchange for identical quarks
if(ix==iy)
add(new_ptr((Tree2toNDiagram(3),_antiquark[ix],_gluon,_antiquark[iy],
2,_antiquark[ix],1,_antiquark[iy],-19)));
}
}
for(unsigned int iy=0;iy<_maxflavour;++iy) {
// q qbar -> q qbar
if(_process==0||_process==8) {
// gluon s-channel
add(new_ptr((Tree2toNDiagram(2),_quark[ix], _antiquark[ix],
1, _gluon, 3, _quark[iy], 3, _antiquark[iy],-20)));
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_gluon,_antiquark[iy],
1,_quark[ix],2,_antiquark[iy],-21)));
}
}
}
}
Selector<const ColourLines *>
MEQCD2to2::colourGeometries(tcDiagPtr diag) const {
// colour lines for gg to gg
static const ColourLines cgggg[12]={ColourLines("1 -2, -1 -3 -5, 5 -4, 2 3 4"),// A_2 s
ColourLines("-1 2, 1 3 5, -5 4, -2 -3 -4"),// A_1 s
ColourLines("1 5, -1 -2 3, -3 -4, -5 2 4"),// A_1 u
ColourLines("-1 -5, 1 2 -3, 3 4, 5 -2 -4"),// A_2 u
ColourLines("1 -2, -1 -3 -4, 4 -5, 2 3 5"),// B_2 s
ColourLines("-1 2, 1 3 4, -4 5, -2 -3 -5"),// B_1 s
ColourLines("1 4, -1 -2 3, -3 -5, -4 2 5"),// B_1 t
ColourLines("-1 -4, 1 2 -3, 3 5, 4 -2 -5"),// B_2 t
ColourLines("1 4, -1 -2 -5, 3 5, -3 2 -4"),// C_1 t
ColourLines("-1 -4, 1 2 5, -3 -5, 3 -2 4"),// C_2 t
ColourLines("1 5, -1 -2 -4, 3 4, -3 2 -5"),// C_1 u
ColourLines("-1 -5, 1 2 4, -3 -4, 3 -2 5") // C_2 u
};
// colour lines for gg to q qbar
static const ColourLines cggqq[4]={ColourLines("1 4, -1 -2 3, -3 -5"),
ColourLines("3 4, -3 -2 1, -1 -5"),
ColourLines("2 -1, 1 3 4, -2 -3 -5"),
ColourLines("1 -2, -1 -3 -5, 2 3 4")};
// colour lines for q qbar to gg
static const ColourLines cqqgg[4]={ColourLines("1 4, -4 -2 5, -3 -5"),
ColourLines("1 5, -3 -4, 4 -2 -5"),
ColourLines("1 3 4, -4 5, -2 -3 -5"),
ColourLines("1 3 5, -5 4, -2 -3 -4")};
// colour lines for q g to q g
static const ColourLines cqgqg[4]={ColourLines("1 -2, 2 3 5, 4 -5"),
ColourLines("1 5, 3 4,-3 2 -5 "),
ColourLines("1 2 -3, 3 5, -5 -2 4"),
ColourLines("1 -2 5,3 2 4,-3 -5")};
// colour lines for qbar g -> qbar g
static const ColourLines cqbgqbg[4]={ColourLines("-1 2, -2 -3 -5, -4 5"),
ColourLines("-1 -5, -3 -4, 3 -2 5"),
ColourLines("-1 -2 3, -3 -5, 5 2 -4"),
ColourLines("-1 2 -5,-3 -2 -4, 3 5")};
// colour lines for q q -> q q
static const ColourLines cqqqq[2]={ColourLines("1 2 5,3 -2 4"),
ColourLines("1 2 4,3 -2 5")};
// colour lines for qbar qbar -> qbar qbar
static const ColourLines cqbqbqbqb[2]={ColourLines("-1 -2 -5,-3 2 -4"),
ColourLines("-1 -2 -4,-3 2 -5")};
// colour lines for q qbar -> q qbar
static const ColourLines cqqbqqb[2]={ColourLines("1 3 4,-2 -3 -5"),
ColourLines("1 2 -3,4 -2 -5")};
// select the colour flow (as all ready picked just insert answer)
Selector<const ColourLines *> sel;
switch(abs(diag->id())) {
// gg -> gg subprocess
case 1:
if(_flow==1) {
sel.insert(0.5, &cgggg[0]);
sel.insert(0.5, &cgggg[1]);
}
else {
sel.insert(0.5, &cgggg[4]);
sel.insert(0.5, &cgggg[5]);
}
break;
case 2:
if(_flow==2) {
sel.insert(0.5, &cgggg[6]);
sel.insert(0.5, &cgggg[7]);
}
else {
sel.insert(0.5, &cgggg[8]);
sel.insert(0.5, &cgggg[9]);
}
break;
case 3:
if(_flow==1) {
sel.insert(0.5, &cgggg[2]);
sel.insert(0.5, &cgggg[3]);
}
else {
sel.insert(0.5, &cgggg[10]);
sel.insert(0.5, &cgggg[11]);
}
break;
// gg -> q qbar subprocess
case 4: case 5:
sel.insert(1.0, &cggqq[abs(diag->id())-4]);
break;
case 6:
sel.insert(1.0, &cggqq[1+_flow]);
break;
// q qbar -> gg subprocess
case 7: case 8:
sel.insert(1.0, &cqqgg[abs(diag->id())-7]);
break;
case 9:
sel.insert(1.0, &cqqgg[1+_flow]);
break;
// q g -> q g subprocess
case 10: case 11:
sel.insert(1.0, &cqgqg[abs(diag->id())-10]);
break;
case 12:
sel.insert(1.0, &cqgqg[1+_flow]);
break;
// q g -> q g subprocess
case 13: case 14:
sel.insert(1.0, &cqbgqbg[abs(diag->id())-13]);
break;
case 15:
sel.insert(1.0, &cqbgqbg[1+_flow]);
break;
// q q -> q q subprocess
case 16: case 17:
sel.insert(1.0, &cqqqq[abs(diag->id())-16]);
break;
// qbar qbar -> qbar qbar subprocess
case 18: case 19:
sel.insert(1.0, &cqbqbqbqb[abs(diag->id())-18]);
break;
// q qbar -> q qbar subprocess
case 20: case 21:
sel.insert(1.0, &cqqbqqb[abs(diag->id())-20]);
break;
}
return sel;
}
double MEQCD2to2::me2() const {
// total matrix element
double me(0.);
// gg initiated processes
if(mePartonData()[0]->id()==ParticleID::g&&mePartonData()[1]->id()==ParticleID::g) {
// gg -> gg
if(mePartonData()[2]->id()==ParticleID::g) {
VectorWaveFunction g1w(meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming);
VectorWaveFunction g3w(meMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction g4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g1,g2,g3,g4;
for(unsigned int ix=0;ix<2;++ix) {
g1w.reset(2*ix);g1.push_back(g1w);
g2w.reset(2*ix);g2.push_back(g2w);
g3w.reset(2*ix);g3.push_back(g3w);
g4w.reset(2*ix);g4.push_back(g4w);
}
// calculate the matrix element
me = gg2ggME(g1,g2,g3,g4,0);
}
// gg -> q qbar
else {
VectorWaveFunction g1w(meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction qw(meMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction qbarw(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g1,g2;
vector<SpinorBarWaveFunction> q;
vector<SpinorWaveFunction> qbar;
for(unsigned int ix=0;ix<2;++ix) {
g1w.reset(2*ix);g1.push_back(g1w);
g2w.reset(2*ix);g2.push_back(g2w);
qw.reset(ix);q.push_back(qw);
qbarw.reset(ix);qbar.push_back(qbarw);
}
// calculate the matrix element
me=gg2qqbarME(g1,g2,q,qbar,0);
}
}
// quark first processes
else if(mePartonData()[0]->id()>0) {
// q g -> q g
if(mePartonData()[1]->id()==ParticleID::g) {
SpinorWaveFunction qinw(meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction qoutw(meMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction g4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g2,g4;
vector<SpinorWaveFunction> qin;
vector<SpinorBarWaveFunction> qout;
for(unsigned int ix=0;ix<2;++ix) {
qinw.reset(ix);qin.push_back(qinw);
g2w.reset(2*ix);g2.push_back(g2w);
qoutw.reset(ix);qout.push_back(qoutw);
g4w.reset(2*ix);g4.push_back(g4w);
}
// calculate the matrix element
me = qg2qgME(qin,g2,qout,g4,0);
}
else if(mePartonData()[1]->id()<0) {
// q qbar initiated processes( q qbar -> gg)
if(mePartonData()[2]->id()==ParticleID::g) {
SpinorWaveFunction qw(meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction qbarw(meMomenta()[1],mePartonData()[1],incoming);
VectorWaveFunction g1w(meMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction g2w(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g1,g2;
vector<SpinorWaveFunction> q;
vector<SpinorBarWaveFunction> qbar;
for(unsigned int ix=0;ix<2;++ix) {
qw.reset(ix);q.push_back(qw);
qbarw.reset(ix);qbar.push_back(qbarw);
g1w.reset(2*ix);g1.push_back(g1w);
g2w.reset(2*ix);g2.push_back(g2w);
}
// calculate the matrix element
me = qqbar2ggME(q,qbar,g1,g2,0);
}
// q qbar to q qbar
else {
SpinorWaveFunction q1w(meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction q2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction q3w(meMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction q4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<SpinorWaveFunction> q1,q4;
vector<SpinorBarWaveFunction> q2,q3;
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1.push_back(q1w);
q2w.reset(ix);q2.push_back(q2w);
q3w.reset(ix);q3.push_back(q3w);
q4w.reset(ix);q4.push_back(q4w);
}
// calculate the matrix element
me = qqbar2qqbarME(q1,q2,q3,q4,0);
}
}
// q q -> q q
else if(mePartonData()[1]->id()>0) {
SpinorWaveFunction q1w(meMomenta()[0],mePartonData()[0],incoming);
SpinorWaveFunction q2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction q3w(meMomenta()[2],mePartonData()[2],outgoing);
SpinorBarWaveFunction q4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<SpinorWaveFunction> q1,q2;
vector<SpinorBarWaveFunction> q3,q4;
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1.push_back(q1w);
q2w.reset(ix);q2.push_back(q2w);
q3w.reset(ix);q3.push_back(q3w);
q4w.reset(ix);q4.push_back(q4w);
}
// calculate the matrix element
me = qq2qqME(q1,q2,q3,q4,0);
}
}
// antiquark first processes
else if(mePartonData()[0]->id()<0) {
// qbar g -> qbar g
if(mePartonData()[1]->id()==ParticleID::g) {
SpinorBarWaveFunction qinw(meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorWaveFunction qoutw(meMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction g4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g2,g4;
vector<SpinorBarWaveFunction> qin;
vector<SpinorWaveFunction> qout;
for(unsigned int ix=0;ix<2;++ix) {
qinw.reset(ix);qin.push_back(qinw);
g2w.reset(2*ix);g2.push_back(g2w);
qoutw.reset(ix);qout.push_back(qoutw);
g4w.reset(2*ix);g4.push_back(g4w);
}
// calculate the matrix element
me = qbarg2qbargME(qin,g2,qout,g4,0);
}
// qbar qbar -> qbar qbar
else if(mePartonData()[1]->id()<0) {
SpinorBarWaveFunction q1w(meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction q2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorWaveFunction q3w(meMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction q4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<SpinorBarWaveFunction> q1,q2;
vector<SpinorWaveFunction> q3,q4;
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1.push_back(q1w);
q2w.reset(ix);q2.push_back(q2w);
q3w.reset(ix);q3.push_back(q3w);
q4w.reset(ix);q4.push_back(q4w);
}
// calculate the matrix element
me = qbarqbar2qbarqbarME(q1,q2,q3,q4,0);
}
}
else throw Exception() << "Unknown process in MEQCD2to2::me2()"
<< Exception::abortnow;
// return the answer
return me;
}
void MEQCD2to2::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard;
hard.push_back(sub->incoming().first);hard.push_back(sub->incoming().second);
hard.push_back(sub->outgoing()[0]);hard.push_back(sub->outgoing()[1]);
// order of particles
unsigned int order[4]={0,1,2,3};
// identify the process and calculate the matrix element
if(hard[0]->id()==ParticleID::g&&hard[1]->id()==ParticleID::g) {
// gg -> gg
if(hard[2]->id()==ParticleID::g) {
vector<VectorWaveFunction> g1,g2,g3,g4;
VectorWaveFunction(g1,hard[0],incoming,false,true,true);
VectorWaveFunction(g2,hard[1],incoming,false,true,true);
VectorWaveFunction(g3,hard[2],outgoing,true ,true,true);
VectorWaveFunction(g4,hard[3],outgoing,true ,true,true);
g1[1]=g1[2];g2[1]=g2[2];g3[1]=g3[2];g4[1]=g4[2];
gg2ggME(g1,g2,g3,g4,_flow);
}
// gg -> q qbar
else {
if(hard[2]->id()<0) swap(order[2],order[3]);
vector<VectorWaveFunction> g1,g2;
vector<SpinorBarWaveFunction> q;
vector<SpinorWaveFunction> qbar;
VectorWaveFunction( g1,hard[ 0 ],incoming,false,true,true);
VectorWaveFunction( g2,hard[ 1 ],incoming,false,true,true);
SpinorBarWaveFunction(q ,hard[order[2]],outgoing,true ,true);
SpinorWaveFunction( qbar,hard[order[3]],outgoing,true ,true);
g1[1]=g1[2];g2[1]=g2[2];
gg2qqbarME(g1,g2,q,qbar,_flow);
}
}
else if(hard[0]->id()==ParticleID::g||hard[1]->id()==ParticleID::g) {
if(hard[0]->id()==ParticleID::g) swap(order[0],order[1]);
if(hard[2]->id()==ParticleID::g) swap(order[2],order[3]);
// q g -> q g
if(hard[order[0]]->id()>0) {
vector<VectorWaveFunction> g2,g4;
vector<SpinorWaveFunction> qin;
vector<SpinorBarWaveFunction> qout;
SpinorWaveFunction( qin,hard[order[0]],incoming,false,true);
VectorWaveFunction( g2,hard[order[1]],incoming,false,true,true);
SpinorBarWaveFunction(qout,hard[order[2]],outgoing,true ,true);
VectorWaveFunction( g4,hard[order[3]],outgoing,true ,true,true);
g2[1]=g2[2];g4[1]=g4[2];
qg2qgME(qin,g2,qout,g4,_flow);
}
// qbar g -> qbar g
else {
vector<VectorWaveFunction> g2,g4;
vector<SpinorBarWaveFunction> qin;
vector<SpinorWaveFunction> qout;
SpinorBarWaveFunction( qin,hard[order[0]],incoming,false,true);
VectorWaveFunction( g2,hard[order[1]],incoming,false,true,true);
SpinorWaveFunction( qout,hard[order[2]],outgoing,true ,true);
VectorWaveFunction( g4,hard[order[3]],outgoing,true ,true,true);
g2[1]=g2[2];g4[1]=g4[2];
qbarg2qbargME(qin,g2,qout,g4,_flow);
}
}
else if(hard[0]->id()>0||hard[1]->id()>0) {
if(hard[2]->id()==ParticleID::g) {
if(hard[0]->id()<0) swap(order[0],order[1]);
vector<SpinorBarWaveFunction> qbar;
vector<SpinorWaveFunction> q;
vector<VectorWaveFunction> g3,g4;
SpinorWaveFunction( q ,hard[order[0]],incoming,false,true);
SpinorBarWaveFunction(qbar,hard[order[1]],incoming,false,true);
VectorWaveFunction( g3,hard[ 2 ],outgoing,true ,true,true);
VectorWaveFunction( g4,hard[ 3 ],outgoing,true ,true,true);
g3[1]=g3[2];g4[1]=g4[2];
qqbar2ggME(q,qbar,g3,g4,_flow);
}
// q q -> q q
else if(hard[0]->id()>0&&hard[1]->id()>0) {
if(hard[2]->id()!=hard[0]->id()) swap(order[2],order[3]);
vector<SpinorWaveFunction> q1,q2;
vector<SpinorBarWaveFunction> q3,q4;
SpinorWaveFunction( q1,hard[order[0]],incoming,false,true);
SpinorWaveFunction( q2,hard[order[1]],incoming,false,true);
SpinorBarWaveFunction(q3,hard[order[2]],outgoing,true ,true);
SpinorBarWaveFunction(q4,hard[order[3]],outgoing,true ,true);
qq2qqME(q1,q2,q3,q4,_flow);
}
// q qbar -> q qbar
else {
if(hard[0]->id()<0) swap(order[0],order[1]);
if(hard[2]->id()<0) swap(order[2],order[3]);
vector<SpinorWaveFunction> q1,q4;
vector<SpinorBarWaveFunction> q2,q3;
SpinorWaveFunction( q1,hard[order[0]],incoming,false,true);
SpinorBarWaveFunction(q2,hard[order[1]],incoming,false,true);
SpinorBarWaveFunction(q3,hard[order[2]],outgoing,true ,true);
SpinorWaveFunction( q4,hard[order[3]],outgoing,true ,true);
qqbar2qqbarME(q1,q2,q3,q4,_flow);
}
}
else if (hard[0]->id()<0&&hard[1]->id()<0) {
if(hard[2]->id()!=hard[0]->id()) swap(order[2],order[3]);
vector<SpinorBarWaveFunction> q1,q2;
vector<SpinorWaveFunction> q3,q4;
SpinorBarWaveFunction(q1,hard[order[0]],incoming,false,true);
SpinorBarWaveFunction(q2,hard[order[1]],incoming,false,true);
SpinorWaveFunction( q3,hard[order[2]],outgoing,true ,true);
SpinorWaveFunction( q4,hard[order[3]],outgoing,true ,true);
qbarqbar2qbarqbarME(q1,q2,q3,q4,_flow);
}
else throw Exception() << "Unknown process in MEQCD2to2::constructVertex()"
<< Exception::runerror;
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(_me);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<4;++ix)
hard[order[ix]]->spinInfo()->productionVertex(hardvertex);
}

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 10:13 AM (1 d, 16 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111068
Default Alt Text
(193 KB)

Event Timeline