Page MenuHomeHEPForge
Files F19250784

No OneTemporary
Actions

Size
37 KB
Referenced Files
None
Subscribers
None

View Options

diff --git a/MatrixElement/MEff2ff.cc b/MatrixElement/MEff2ff.cc
--- a/MatrixElement/MEff2ff.cc
+++ b/MatrixElement/MEff2ff.cc
@@ -1,752 +1,757 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEff2ff class.
//
#include "MEff2ff.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig++/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "Herwig++/Helicity/WaveFunction/TensorWaveFunction.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig++/Helicity/Correlations/HardVertex.h"
#include <numeric>
using namespace Herwig;
using Herwig::Helicity::VectorWaveFunction;
using Herwig::Helicity::ScalarWaveFunction;
using Herwig::Helicity::TensorWaveFunction;
using Herwig::Helicity::incoming;
using Herwig::Helicity::outgoing;
using Herwig::Helicity::HardVertex;
using Herwig::Helicity::HardVertexPtr;
using ThePEG::Helicity::SpinfoPtr;
double MEff2ff::me2() const {
tcPDPtr ina(mePartonData()[0]), inb(mePartonData()[1]),
outa(mePartonData()[2]), outb(mePartonData()[3]);
bool majorana(false);
if( (!outa->CC() && !outb->CC() ) ||
((abs(outa->id()) > 1000000 && abs(outa->id()) < 2000000) &&
(abs(outb->id()) > 1000000 && abs(outb->id()) < 2000000)) )
majorana = true;
double full_me(0.);
vector<SpinorWaveFunction> spA(2), spB(2);
vector<SpinorBarWaveFunction> spbA(2), spbB(2);
if( ina->id() > 0 && inb->id() < 0) {
for(unsigned int ih = 0; ih < 2; ++ih) {
spA[ih] = SpinorWaveFunction(meMomenta()[0], ina, ih,
incoming);
spbA[ih] = SpinorBarWaveFunction(meMomenta()[1], inb, ih,
incoming);
spB[ih] = SpinorWaveFunction(meMomenta()[3], outb, ih, outgoing);
spbB[ih] = SpinorBarWaveFunction(meMomenta()[2], outa, ih, outgoing);
}
if(majorana) {
vector<SpinorWaveFunction> spC(2);
vector<SpinorBarWaveFunction> spbC(2);
for(unsigned int ih = 0; ih < 2; ++ih) {
spC[ih] = SpinorWaveFunction(meMomenta()[2], outa, ih, outgoing);
spbC[ih] = SpinorBarWaveFunction(meMomenta()[3], outb, ih, outgoing);
}
ffb2mfmfHeME(spA, spbA, spbB, spB, spC, spbC, full_me);
SpinorWaveFunction spOut2(meMomenta()[2], outa, outgoing);
SpinorBarWaveFunction spbarOut2(meMomenta()[3], outb, outgoing);
}
else {
ffb2ffbHeME(spA, spbA, spbB, spB, full_me);
}
}
else if( ina->id() > 0 && inb->id() > 0 ) {
SpinorVector spA(2), spB(2);
SpinorBarVector spbA(2), spbB(2);
for(unsigned int ih = 0; ih < 2; ++ih) {
spA[ih] = SpinorWaveFunction(meMomenta()[0], ina, ih,
incoming);
spB[ih] = SpinorWaveFunction(meMomenta()[1], inb, ih,
incoming);
spbA[ih] = SpinorBarWaveFunction(meMomenta()[2], outa, ih,
outgoing);
spbB[ih] = SpinorBarWaveFunction(meMomenta()[3], outb, ih,
outgoing);
}
ff2ffHeME(spA, spB, spbA, spbB, full_me);
}
else if( ina->id() < 0 && inb->id() < 0 ) {
SpinorVector spA(2), spB(2);
SpinorBarVector spbA(2), spbB(2);
for(unsigned int ih = 0; ih < 2; ++ih) {
spbA[ih] = SpinorBarWaveFunction(meMomenta()[0], ina, ih,
incoming);
spbB[ih] = SpinorBarWaveFunction(meMomenta()[1], inb, ih,
incoming);
spA[ih] = SpinorWaveFunction(meMomenta()[2], outa, ih,
outgoing);
spB[ih] = SpinorWaveFunction(meMomenta()[3], outb, ih,
outgoing);
}
fbfb2fbfbHeME(spbA, spbB, spA, spB, full_me);
}
else
throw MEException()
<< "MEff2ff::me2() - Cannot find correct function to deal with process "
<< ina->PDGName() << "," << inb->PDGName() << "->" << outa->PDGName()
<< "," << outb->PDGName() << "\n";
return full_me;
}
ProductionMatrixElement
MEff2ff::ffb2ffbHeME(SpinorVector & fin, SpinorBarVector & fbin,
SpinorBarVector & fbout, SpinorVector & fout,
double & me2) const {
const HPCount ndiags = numberOfDiags();
const size_t ncf(numberOfFlows());
const vector<vector<double> > cfactors(getColourFactors());
const Energy2 m2(scale());
vector<Complex> diag(ndiags, Complex(0.)), flows(ncf, Complex(0.));
vector<double> me(ndiags, 0.);
ScalarWaveFunction interS; VectorWaveFunction interV;
TensorWaveFunction interT;
ProductionMatrixElement prodME(PDT::Spin1Half, PDT::Spin1Half, PDT::Spin1Half,
PDT::Spin1Half);
for(unsigned int ifhel1 = 0; ifhel1 < 2; ++ifhel1) {
for(unsigned int ifhel2 = 0; ifhel2 < 2; ++ifhel2) {
for(unsigned int ofhel1 = 0; ofhel1 < 2; ++ofhel1) {
for(unsigned int ofhel2 = 0; ofhel2 < 2; ++ofhel2) {
flows = vector<Complex>(ncf, Complex(0.));
for(HPCount ix = 0; ix < ndiags; ++ix) {
HPDiagram current = getProcessInfo()[ix];
tcPDPtr offshell = current.intermediate;
if(current.channelType == HPDiagram::tChannel) {
if(offshell->iSpin() == PDT::Spin0) {
interS = theScaV[ix].second->evaluate(m2, 3, offshell,
fout[ofhel2], fbin[ifhel2]);
diag[ix] = theScaV[ix].first->evaluate(m2, fin[ifhel1],
fbout[ofhel1], interS);
}
else if(offshell->iSpin() == PDT::Spin1) {
interV = theVecV[ix].second->evaluate(m2, 3, offshell,
fout[ofhel2], fbin[ifhel2]);
diag[ix] = -theVecV[ix].first->evaluate(m2, fin[ifhel1],
fbout[ofhel1], interV);
}
else if(offshell->iSpin() == PDT::Spin2) {
interT = theTenV[ix].second->evaluate(m2, 3, offshell,
fout[ofhel2], fbin[ifhel2]);
diag[ix] = theTenV[ix].first->evaluate(m2, fin[ifhel1],
fbout[ofhel1], interT);
}
}
else if(current.channelType == HPDiagram::sChannel) {
if(offshell->iSpin() == PDT::Spin0) {
interS = theScaV[ix].second->evaluate(m2, 1, offshell,
fout[ofhel2],fbout[ofhel1]);
diag[ix] = theScaV[ix].first->evaluate(m2, fin[ifhel1],
fbin[ifhel2], interS);
}
else if(offshell->iSpin() == PDT::Spin1) {
interV = theVecV[ix].second->evaluate(m2, 1, offshell,
fout[ofhel2],fbout[ofhel1]);
diag[ix] = theVecV[ix].first->evaluate(m2, fin[ifhel1],
fbin[ifhel2], interV);
}
else if(offshell->iSpin() == PDT::Spin2) {
interT = theTenV[ix].second->evaluate(m2, 1, offshell,
fout[ofhel2],fbout[ofhel1]);
diag[ix] = theTenV[ix].first->evaluate(m2, fin[ifhel1],
fbin[ifhel2], interT);
}
}
else {
throw MEException() << "Incorrect diagram type in matrix element "
<< fullName()
<< Exception::warning;
diag[ix] = 0.;
}
me[ix] += norm(diag[ix]);
//Compute flows
for(size_t iy = 0; iy < current.colourFlow.size(); ++iy)
flows[current.colourFlow[iy].first - 1] +=
current.colourFlow[iy].second * diag[ix];
}//diagram loop
//Now add flows to me2 with appropriate colour factors
for(size_t ii = 0; ii < ncf; ++ii)
for(size_t ij = 0; ij < ncf; ++ij)
me2 += cfactors[ii][ij]*(flows[ii]*conj(flows[ij])).real();
prodME(ifhel1, ifhel2, ofhel1, ofhel2) =
std::accumulate(flows.begin(), flows.end(), Complex(0., 0.));
}//end of first helicity loop
}
}
}
const double identfact = mePartonData()[2]->id() == mePartonData()[3]->id()
? 0.5 : 1.;
const double colfact = mePartonData()[0]->iColour() == PDT::Colour3 ?
1./9. : 1;
DVector save(ndiags);
for(DVector::size_type ix = 0; ix < ndiags; ++ix)
save[ix] = 0.25*identfact*colfact*me[ix];
meInfo(save);
me2 *= 0.25*identfact*colfact;
return prodME;
}
ProductionMatrixElement
MEff2ff:: ff2ffHeME(SpinorVector & fin, SpinorVector & fin2,
SpinorBarVector & fbout, SpinorBarVector & fbout2,
double & me2) const {
const HPCount ndiags = getProcessInfo().size();
const size_t ncf(numberOfFlows());
const vector<vector<double> > cfactors(getColourFactors());
const Energy2 q2(scale());
vector<Complex> diag(ndiags, Complex(0.)), flows(ncf, Complex(0.));
vector<double> me(ndiags, 0.);
ScalarWaveFunction interS; VectorWaveFunction interV;
TensorWaveFunction interT;
ProductionMatrixElement prodME(PDT::Spin1Half, PDT::Spin1Half, PDT::Spin1Half,
PDT::Spin1Half);
for(unsigned int ifhel1 = 0; ifhel1 < 2; ++ifhel1) {
for(unsigned int ifhel2 = 0; ifhel2 < 2; ++ifhel2) {
for(unsigned int ofhel1 = 0; ofhel1 < 2; ++ofhel1) {
for(unsigned int ofhel2 = 0; ofhel2 < 2; ++ofhel2) {
flows = vector<Complex>(ncf, Complex(0.));
for(HPCount ix = 0; ix < ndiags; ++ix) {
HPDiagram current = getProcessInfo()[ix];
tcPDPtr offshell = current.intermediate;
if(current.channelType == HPDiagram::tChannel) {
if(offshell->iSpin() == PDT::Spin0) {
if(current.ordered.second) {
interS = theScaV[ix].second->evaluate(q2, 3, offshell,
fin2[ifhel2],
fbout2[ofhel2]);
diag[ix] = theScaV[ix].first->evaluate(q2, fin[ifhel1],
fbout[ofhel1], interS);
}
else {
interS = theScaV[ix].second->evaluate(q2, 3, offshell,
fin2[ifhel2],
fbout[ofhel1]);
diag[ix] = theScaV[ix].first->evaluate(q2, fin[ifhel1],
fbout2[ofhel2], interS);
}
}
else if(offshell->iSpin() == PDT::Spin1) {
if(current.ordered.second) {
interV = theVecV[ix].second->evaluate(q2, 3, offshell,
fin2[ifhel2],
fbout2[ofhel2]);
diag[ix] = theVecV[ix].first->evaluate(q2, fin[ifhel1],
fbout[ofhel1], interV);
}
else {
interV = theVecV[ix].second->evaluate(q2, 3, offshell,
fin2[ifhel2],
fbout[ofhel1]);
diag[ix] = -theVecV[ix].first->evaluate(q2, fin[ifhel1],
fbout2[ofhel2], interV);
}
}
else if(offshell->iSpin() == PDT::Spin2) {
if(current.ordered.second) {
interT = theTenV[ix].second->evaluate(q2, 3, offshell,
fin2[ifhel2],
fbout2[ofhel2]);
diag[ix] = theTenV[ix].first->evaluate(q2, fin[ifhel1],
fbout[ofhel1], interT);
}
else {
interT = theTenV[ix].second->evaluate(q2, 3, offshell,
fin2[ifhel2],
fbout[ofhel1]);
diag[ix] = theTenV[ix].first->evaluate(q2, fin[ifhel1],
fbout2[ofhel2],
interT);
}
}
}
else {
throw MEException() << "Incorrect diagram type in matrix element "
<< fullName()
<< Exception::warning;
diag[ix] = 0.;
}
me[ix] += norm(diag[ix]);
//Compute flows
for(size_t iy = 0; iy < current.colourFlow.size(); ++iy)
flows[current.colourFlow[iy].first - 1] +=
current.colourFlow[iy].second * diag[ix];
}//diagram loop
//Now add flows to me2 with appropriate colour factors
for(size_t ii = 0; ii < ncf; ++ii)
for(size_t ij = 0; ij < ncf; ++ij)
me2 += cfactors[ii][ij]*(flows[ii]*conj(flows[ij])).real();
prodME(ifhel1, ifhel2, ofhel1, ofhel2) =
std::accumulate(flows.begin(), flows.end(), Complex(0., 0.));
}
}
}
}
const double identfact = mePartonData()[2]->id() == mePartonData()[3]->id()
? 0.5 : 1.;
const double colfact = mePartonData()[0]->iColour() == PDT::Colour3 ? 1./9. : 1;
DVector save(ndiags);
for(DVector::size_type ix = 0; ix < ndiags; ++ix)
save[ix] = 0.25*identfact*colfact*me[ix];
meInfo(save);
me2 = 0.25*identfact*colfact*me2;
return prodME;
}
ProductionMatrixElement
MEff2ff::fbfb2fbfbHeME(SpinorBarVector & fbin, SpinorBarVector & fbin2,
SpinorVector & fout, SpinorVector & fout2,
double & me2) const {
const HPCount ndiags = getProcessInfo().size();
const size_t ncf(numberOfFlows());
const vector<vector<double> > cfactors(getColourFactors());
const Energy2 q2(scale());
vector<Complex> diag(ndiags, Complex(0.)), flows(ncf, Complex(0.));
vector<double> me(ndiags, 0.);
ScalarWaveFunction interS; VectorWaveFunction interV;
TensorWaveFunction interT;
ProductionMatrixElement prodME(PDT::Spin1Half, PDT::Spin1Half, PDT::Spin1Half,
PDT::Spin1Half);
for(unsigned int ifhel1 = 0; ifhel1 < 2; ++ifhel1) {
for(unsigned int ifhel2 = 0; ifhel2 < 2; ++ifhel2) {
for(unsigned int ofhel1 = 0; ofhel1 < 2; ++ofhel1) {
for(unsigned int ofhel2 = 0; ofhel2 < 2; ++ofhel2) {
flows = vector<Complex>(ncf, Complex(0.));
for(HPCount ix = 0; ix < ndiags; ++ix) {
HPDiagram current = getProcessInfo()[ix];
tcPDPtr offshell = current.intermediate;
if(current.channelType == HPDiagram::tChannel) {
if(offshell->iSpin() == PDT::Spin0) {
if(current.ordered.second) {
interS = theScaV[ix].second->evaluate(q2, 3, offshell,
fout2[ofhel2],
fbin2[ifhel2]);
diag[ix] = theScaV[ix].first->evaluate(q2, fout[ofhel1],
fbin[ifhel1], interS);
}
else {
interS = theScaV[ix].second->evaluate(q2, 3, offshell,
fout[ofhel1],
fbin2[ifhel2]);
diag[ix] = -theScaV[ix].first->evaluate(q2, fout2[ofhel2],
fbin[ifhel1], interS);
}
}
else if(offshell->iSpin() == PDT::Spin1) {
if(current.ordered.second) {
interV = theVecV[ix].second->evaluate(q2, 3, offshell,
fout2[ofhel2],
fbin2[ifhel2]);
diag[ix] = theVecV[ix].first->evaluate(q2, fout[ofhel1],
fbin[ifhel1], interV);
}
else {
interV = theVecV[ix].second->evaluate(q2, 3, offshell,
fout[ofhel1],
fbin2[ifhel2]);
diag[ix] = -theVecV[ix].first->evaluate(q2, fout2[ofhel2],
fbin[ifhel1], interV);
}
}
else if(offshell->iSpin() == PDT::Spin2) {
if(current.ordered.second) {
interT = theTenV[ix].second->evaluate(q2, 3, offshell,
fout2[ofhel2],
fbin2[ifhel2]);
diag[ix] = theTenV[ix].first->evaluate(q2, fout[ofhel1],
fbin[ifhel1], interT);
}
else {
interT = theTenV[ix].second->evaluate(q2, 3, offshell,
fout[ofhel1],
fbin2[ifhel2]);
diag[ix] = -theTenV[ix].first->evaluate(q2, fout2[ofhel2],
fbin[ifhel1], interT);
}
}
}
me[ix] += norm(diag[ix]);
//Compute flows
for(size_t iy = 0; iy < current.colourFlow.size(); ++iy)
flows[current.colourFlow[iy].first - 1] +=
current.colourFlow[iy].second * diag[ix];
}//diagram loop
//Now add flows to me2 with appropriate colour factors
for(size_t ii = 0; ii < ncf; ++ii)
for(size_t ij = 0; ij < ncf; ++ij)
me2 += cfactors[ii][ij]*(flows[ii]*conj(flows[ij])).real();
prodME(ifhel1, ifhel2, ofhel1, ofhel2) =
std::accumulate(flows.begin(), flows.end(), Complex(0., 0.));
}
}
}
}
const double identfact = mePartonData()[2]->id() == mePartonData()[3]->id()
? 0.5 : 1.;
const double colfact = (mePartonData()[0]->iColour() == PDT::Colour3bar)
? 1./9. : 1;
DVector save(ndiags);
for(DVector::size_type ix = 0; ix < ndiags; ++ix)
save[ix] = 0.25*identfact*colfact*me[ix];
meInfo(save);
me2 *= 0.25*identfact*colfact;
return prodME;
}
ProductionMatrixElement
MEff2ff::ffb2mfmfHeME(SpinorVector & fin, SpinorBarVector & fbin,
SpinorBarVector & fbout, SpinorVector & fout,
SpinorVector & fout2, SpinorBarVector & fbout2,
double & me2) const {
//useful constants
const HPCount ndiags = numberOfDiags();
const size_t ncf(numberOfFlows());
const vector<vector<double> > cfactors(getColourFactors());
const Energy2 m2(scale());
//store results
vector<Complex> diag(ndiags, Complex(0.)), flows(ncf, Complex(0.));
vector<double> me(ndiags, 0.);
//intermediate wavefunctions
ScalarWaveFunction interS; VectorWaveFunction interV;
TensorWaveFunction interT;
//ProductionMatrixElement object
ProductionMatrixElement prodME(PDT::Spin1Half, PDT::Spin1Half,
PDT::Spin1Half, PDT::Spin1Half);
for(unsigned int ifhel1 = 0; ifhel1 < 2; ++ifhel1) {
for(unsigned int ifhel2 = 0; ifhel2 < 2; ++ifhel2) {
for(unsigned int ofhel1 = 0; ofhel1 < 2; ++ofhel1) {
for(unsigned int ofhel2 = 0; ofhel2 < 2; ++ofhel2) {
flows = vector<Complex>(ncf, Complex(0.));
for(size_t ix = 0; ix < ndiags; ++ix) {
HPDiagram current = getProcessInfo()[ix];
tcPDPtr offshell = current.intermediate;
if(current.channelType == HPDiagram::tChannel) {
if(offshell->iSpin() == PDT::Spin0) {
if(current.ordered.second) {
interS = theScaV[ix].second->evaluate(m2, 3, offshell,
fout[ofhel2],
fbin[ifhel2]);
diag[ix] = theScaV[ix].first->evaluate(m2, fin[ifhel1],
fbout[ofhel1],
interS);
}
else {
interS = theScaV[ix].second->evaluate(m2, 3, offshell,
fout2[ofhel1],
fbin[ifhel2]);
diag[ix] = -theScaV[ix].first->evaluate(m2, fin[ifhel1],
fbout2[ofhel2],
interS);
}
}
else if(offshell->iSpin() == PDT::Spin1) {
if(current.ordered.second) {
interV = theVecV[ix].second->evaluate(m2, 3, offshell,
fout[ofhel2],
fbin[ifhel2]);
diag[ix] = theVecV[ix].first->evaluate(m2, fin[ifhel1],
fbout[ofhel1],
interV);
}
else {
interV = theVecV[ix].second->evaluate(m2, 3, offshell,
fout2[ofhel1],
fbin[ifhel2]);
diag[ix] = theVecV[ix].first->evaluate(m2, fin[ifhel1],
fbout2[ofhel2],
interV);
}
}
else if(offshell->iSpin() == PDT::Spin2) {
if(current.ordered.second) {
interT = theTenV[ix].second->evaluate(m2, 3, offshell,
fout[ofhel2],
fbin[ifhel2]);
diag[ix] = theTenV[ix].first->evaluate(m2, fin[ifhel1],
fbout[ofhel1],
interT);
}
else {
interT = theTenV[ix].second->evaluate(m2, 3, offshell,
fout2[ofhel1],
fbin[ifhel2]);
diag[ix] = theTenV[ix].first->evaluate(m2, fin[ifhel1],
fbout2[ofhel2],
interT);
}
}
}
else if(current.channelType == HPDiagram::sChannel) {
if(offshell->iSpin() == PDT::Spin0) {
interS = theScaV[ix].second->evaluate(m2, 1, offshell,
fout[ofhel2],fbout[ofhel1]);
diag[ix] = theScaV[ix].first->evaluate(m2, fin[ifhel1],
fbin[ifhel2], interS);
}
else if(offshell->iSpin() == PDT::Spin1) {
interV = theVecV[ix].second->evaluate(m2, 1, offshell,
fout[ofhel2],fbout[ofhel1]);
diag[ix] = theVecV[ix].first->evaluate(m2, fin[ifhel1],
fbin[ifhel2], interV);
}
else if(offshell->iSpin() == PDT::Spin2) {
interT = theTenV[ix].second->evaluate(m2, 1, offshell,
fout[ofhel2],fbout[ofhel1]);
diag[ix] = theTenV[ix].first->evaluate(m2, fin[ifhel1],
fbin[ifhel2], interT);
}
}
else {
throw MEException() << "Incorrect diagram type in matrix element "
<< fullName()
<< Exception::warning;
diag[ix] = 0.;
}
me[ix] += norm(diag[ix]);
// Compute flows
for(size_t iy = 0; iy < current.colourFlow.size(); ++iy)
flows[current.colourFlow[iy].first - 1] +=
current.colourFlow[iy].second * diag[ix];
}//diagram loop
//Now add flows to me2 with appropriate colour factors
for(size_t ii = 0; ii < ncf; ++ii)
for(size_t ij = 0; ij < ncf; ++ij)
me2 += cfactors[ii][ij]*(flows[ii]*conj(flows[ij])).real();
prodME(ifhel1, ifhel2, ofhel1, ofhel2) =
std::accumulate(flows.begin(), flows.end(), Complex(0., 0.));
}// ofhel2 loop
}
}
}
const double identfact = mePartonData()[2]->id() == mePartonData()[3]->id()
? 0.5 : 1.;
const double colfact = mePartonData()[0]->coloured() ? 1./9. : 1;
DVector save(ndiags);
for(DVector::size_type ix = 0; ix < ndiags; ++ix)
save[ix] = 0.25*identfact*colfact*me[ix];
meInfo(save);
me2 *= 0.25*identfact*colfact;
return prodME;
}
Selector<const ColourLines *>
MEff2ff::colourGeometries(tcDiagPtr diag) const {
static vector<ColourLines> cf(24);
//33b->11
cf[0] = ColourLines("1 2 -3");
cf[1] = ColourLines("1 -2");
//33b->18
cf[2] = ColourLines("1 2 5, -3 -5");
cf[3] = ColourLines("1 4, -4 2 -3");
//33b->33bar
cf[4] = ColourLines("1 2 -3, 4 -2 -5");
cf[5] = ColourLines("1 3 4, -2 -3 -5");
cf[6] = ColourLines("1 4, -3 -5");
cf[7] = ColourLines("1 -2, 4 -5");
//33b->88
cf[8] = ColourLines("1 4, -4 2 5, -5 -3");
cf[9] = ColourLines("1 5, -5 2 4, -4 -3");
cf[10] = ColourLines("1 3 4, -5 -3 -2, -4 5");
cf[11] = ColourLines("1 3 5, -4 -3 -2, -5 4");
//33->33
cf[12] = ColourLines("1 2 5, 3 -2 4");
cf[13] = ColourLines("1 2 4, 3 -2 5");
cf[14] = ColourLines("1 4, 3 5");
cf[15] = ColourLines("1 5, 3 4");
//3b->3b
cf[16] = ColourLines("-1 -2 -5, -3 2 -4");
cf[17] = ColourLines("-1 -2 -4, -3 2 -5");
cf[18] = ColourLines("-1 -4, -3 -5");
cf[19] = ColourLines("-1 -5, -3 -4");
//33b->81
cf[20]= ColourLines("1 4, -4 2 -3");
cf[21]= ColourLines("-3 -4, 1 2 4");
//11->33bar
cf[22] = ColourLines("4 -5");
//11->11
cf[23] = ColourLines("");
HPDiagram current = getProcessInfo()[abs(diag->id()) - 1];
//select appropriate set of diagrams
PDT::Colour inac(mePartonData()[0]->iColour());
PDT::Colour inbc(mePartonData()[1]->iColour());
PDT::Colour outac(mePartonData()[2]->iColour());
PDT::Colour outbc(mePartonData()[3]->iColour());
vector<ColourLines>::size_type cl(0);
if(inac == PDT::Colour3 && inbc == PDT::Colour3) {
if(current.intermediate->iColour() == PDT::Colour8)
cl = current.ordered.second ? 12 : 13;
else
cl = current.ordered.second ? 14 : 15;
}
else if(inac == PDT::Colour3bar && inbc == PDT::Colour3bar) {
if(current.intermediate->iColour() == PDT::Colour8)
cl = current.ordered.second ? 16 : 17;
else
cl = current.ordered.second ? 18 : 19;
}
else if(inac == PDT::Colour0 && inbc == PDT::Colour0)
cl = (outac == PDT::Colour0) ? 23 : 22;
else {
if(outac == PDT::Colour0 || outbc == PDT::Colour0 ) {
if(current.channelType == HPDiagram::tChannel) {
if(outac == outbc)
cl = current.ordered.second ? 0 : 1;
else if(outac == PDT::Colour0 && outbc == PDT::Colour8)
cl = current.ordered.second ? 2 : 3;
else
cl = current.ordered.second ? 20 : 21;
}
}
else if(outbc == PDT::Colour3bar) {
if(current.channelType == HPDiagram::tChannel)
cl = (current.intermediate->iColour() == PDT::Colour8) ? 4 : 6;
else
cl = (current.intermediate->iColour() == PDT::Colour8) ? 5 : 7;
}
else {
if(current.channelType == HPDiagram::tChannel)
cl = current.ordered.second ? 8 : 9;
else
cl = 10 + rnd2(0.5, 0.5);
}
}
Selector<const ColourLines *> sel;
sel.insert(1., &cf[cl]);
return sel;
}
void MEff2ff::constructVertex(tSubProPtr subp) {
// Hard proces external particles
ParticleVector hardpro(4);
hardpro[0] = subp->incoming().first;
hardpro[1] = subp->incoming().second;
hardpro[2] = subp->outgoing()[0];
hardpro[3] = subp->outgoing()[1];
+
//particle ordering. If q qb initial q = 0, qb = 1
if( hardpro[0]->id() < hardpro[1]->id() )
swap(hardpro[0], hardpro[1]);
if( hardpro[2]->id() < hardpro[3]->id() )
swap(hardpro[2], hardpro[3]);
//pick which process we are doing
if( hardpro[0]->id() > 0) {
- SpinorVector spA, spB;
- SpinorBarVector spbA, spbB;
+ //common spinors
+ SpinorVector spA;
+ SpinorBarVector spbB;
SpinorWaveFunction(spA, hardpro[0], incoming, false, true);
- SpinorBarWaveFunction(spbA, hardpro[1], incoming, false, true);
SpinorBarWaveFunction(spbB, hardpro[2], outgoing,true, true);
- SpinorWaveFunction(spB, hardpro[3], outgoing, true, true);
//majorana
if(!hardpro[2]->dataPtr()->CC() || hardpro[2]->id() == 1000024 ||
hardpro[2]->id() == 1000037) {
- SpinorVector spC;
- SpinorBarVector spbC;
+ SpinorVector spB, spC;
+ SpinorBarVector spbA, spbC;
+ SpinorBarWaveFunction(spbA, hardpro[1], incoming, false, true);
+ SpinorWaveFunction(spB, hardpro[3], outgoing, true, true);
for(unsigned int ix=0;ix<2;++ix) {
spC.push_back(SpinorWaveFunction(-spbB[ix].getMomentum(),
spbB[ix].getParticle(),
spbB[ix].wave().bar().conjugate(),
spbB[ix].direction()));
spbC.push_back(SpinorBarWaveFunction(-spB[ix].getMomentum(),
spB[ix].getParticle(),
spB[ix].wave().bar().conjugate(),
spB[ix].direction()));
}
double me;
ProductionMatrixElement prodME = ffb2mfmfHeME(spA, spbA, spbB, spB, spC,
spbC, me);
HardVertexPtr hardvertex = new_ptr(HardVertex());
hardvertex->ME(prodME);
for(ParticleVector::size_type i = 0; i < 4; ++i) {
dynamic_ptr_cast<SpinfoPtr>(hardpro[i]->spinInfo())->
setProductionVertex(hardvertex);
}
}
//ffbar->ffbar
else if( hardpro[1]->id() < 0 ) {
+ SpinorVector spB;
+ SpinorBarVector spbA;
+ SpinorBarWaveFunction(spbA, hardpro[1], incoming, false, true);
+ SpinorWaveFunction(spB, hardpro[3], outgoing, true, true);
double dummy;
ProductionMatrixElement prodME = ffb2ffbHeME(spA, spbA, spbB, spB,
dummy);
HardVertexPtr hardvertex = new_ptr(HardVertex());
hardvertex->ME(prodME);
for(ParticleVector::size_type i = 0; i < 4; ++i)
dynamic_ptr_cast<SpinfoPtr>(hardpro[i]->spinInfo())->
setProductionVertex(hardvertex);
}
//ff2ff
else {
- SpinorVector spA, spB;
- SpinorBarVector spbA, spbB;
- SpinorWaveFunction(spA,hardpro[0],incoming, false, true);
+ SpinorVector spB;
+ SpinorBarVector spbA;
SpinorWaveFunction(spB,hardpro[1],incoming, false, true);
- SpinorBarWaveFunction(spbA, hardpro[2], outgoing, true, true);
- SpinorBarWaveFunction(spbB, hardpro[3], outgoing, true, true);
+ SpinorBarWaveFunction(spbA, hardpro[3], outgoing, true, true);
+
double dummy;
- ProductionMatrixElement prodME = ff2ffHeME(spA, spB, spbA, spbB,
+ ProductionMatrixElement prodME = ff2ffHeME(spA, spB, spbB, spbA,
dummy);
HardVertexPtr hardvertex = new_ptr(HardVertex());
hardvertex->ME(prodME);
for(ParticleVector::size_type i = 0; i < 4; ++i)
dynamic_ptr_cast<SpinfoPtr>(hardpro[i]->spinInfo())->
setProductionVertex(hardvertex);
}
}
//fbarfbar->fbarfbar
else {
SpinorVector spA, spB;
SpinorBarVector spbA, spbB;
SpinorBarWaveFunction(spbA,hardpro[0],incoming, false, true);
SpinorBarWaveFunction(spbB,hardpro[1],incoming, false, true);
SpinorWaveFunction(spA, hardpro[2], outgoing, true, true);
SpinorWaveFunction(spB, hardpro[3], outgoing, true, true);
double dummy;
ProductionMatrixElement prodME = fbfb2fbfbHeME(spbA, spbB, spA, spB,
dummy);
HardVertexPtr hardvertex = new_ptr(HardVertex());
hardvertex->ME(prodME);
for(ParticleVector::size_type i = 0; i < 4; ++i)
dynamic_ptr_cast<SpinfoPtr>(hardpro[i]->spinInfo())->
setProductionVertex(hardvertex);
}
}
void MEff2ff::persistentOutput(PersistentOStream & os) const {
os << theScaV << theVecV << theTenV;
}
void MEff2ff::persistentInput(PersistentIStream & is, int) {
is >> theScaV >> theVecV >> theTenV;
}
ClassDescription<MEff2ff> MEff2ff::initMEff2ff;
// Definition of the static class description member.
void MEff2ff::Init() {
static ClassDocumentation<MEff2ff> documentation
("This is the implementation of the matrix element for fermion-"
"antifermion -> fermion-antifermion.");
}
diff --git a/MatrixElement/MEfv2vf.cc b/MatrixElement/MEfv2vf.cc
--- a/MatrixElement/MEfv2vf.cc
+++ b/MatrixElement/MEfv2vf.cc
@@ -1,288 +1,288 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEfv2vf class.
//
#include "MEfv2vf.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include <numeric>
using namespace Herwig;
using Herwig::Helicity::FFVVertexPtr;
using Herwig::Helicity::VVVVertexPtr;
using Herwig::Helicity::incoming;
using Herwig::Helicity::outgoing;
using Herwig::Helicity::SpinorWaveFunction;
using Herwig::Helicity::SpinorBarWaveFunction;
using Herwig::Helicity::VectorWaveFunction;
double MEfv2vf::me2() const {
//wavefunctions
SpinorVector sp(2);
VBVector vecIn(2), vecOut(3);
SpinorBarVector spb(2);
double fullme(0.);
if(mePartonData()[0]->id() > 0) {
for(unsigned int i = 0; i < 2; ++i) {
sp[i] = SpinorWaveFunction(meMomenta()[0], mePartonData()[0], i, incoming);
vecIn[i] = VectorWaveFunction(meMomenta()[1], mePartonData()[1], 2*i,
incoming);
vecOut[2*i] = VectorWaveFunction(meMomenta()[2], mePartonData()[2], 2*i,
outgoing);
spb[i] = SpinorBarWaveFunction(meMomenta()[3], mePartonData()[3], i,
outgoing);
}
if(mePartonData()[2]->mass() > 0.0)
vecOut[1] = VectorWaveFunction(meMomenta()[2], mePartonData()[2], 1,
outgoing);
fv2vfHeME(sp, vecIn, vecOut, spb, fullme);
}
else {
for(unsigned int i = 0; i < 2; ++i) {
spb[i] = SpinorBarWaveFunction(meMomenta()[0], mePartonData()[0], i,
incoming);
vecIn[i] = VectorWaveFunction(meMomenta()[1], mePartonData()[1], 2*i,
incoming);
vecOut[2*i] = VectorWaveFunction(meMomenta()[2], mePartonData()[2], 2*i,
outgoing);
sp[i] = SpinorWaveFunction(meMomenta()[3], mePartonData()[3], i, outgoing);
}
if(mePartonData()[2]->mass() > 0.0)
vecOut[1] = VectorWaveFunction(meMomenta()[2], mePartonData()[2], 1,
outgoing);
fbv2vfbHeME(spb, vecIn, vecOut, sp, fullme);
}
return fullme;
}
ProductionMatrixElement
MEfv2vf::fv2vfHeME(SpinorVector & spIn, VBVector & vecIn,
VBVector & vecOut, SpinorBarVector & spbOut,
double & mesq) const {
const HPCount ndiags = numberOfDiags();
const size_t ncf(numberOfFlows());
const vector<vector<double> > cfactors(getColourFactors());
const Energy2 q2(scale());
const bool masslessC = mePartonData()[2]->mass() == 0.0;
ProductionMatrixElement prodME(PDT::Spin1Half, PDT::Spin1, PDT::Spin1,
PDT::Spin1Half);
vector<Complex> diag(ndiags, Complex(0., 0.));
vector<double> me(ndiags, 0.);
//intermediate wave functions
SpinorBarWaveFunction interFB; VectorWaveFunction interV;
SpinorWaveFunction interF;
//loop over helicities
for(unsigned int ifh = 0; ifh < 2; ++ifh) {
for(unsigned int ivh = 0; ivh < 2; ++ivh) {
for(unsigned int ovh = 0; ovh < 3; ++ovh) {
if(masslessC && ovh == 1) ++ovh;
for(unsigned int ofh = 0; ofh < 2; ++ofh) {
//store the flows
vector<Complex> cflows(ncf, Complex(0. ,0.));
for(HPCount ix = 0; ix < ndiags; ++ix) {
HPDiagram diagram = getProcessInfo()[ix];
tcPDPtr offshell = diagram.intermediate;
if(diagram.channelType == HPDiagram::tChannel) {
//t-chan spin-1/2
if(offshell->iSpin() == PDT::Spin1Half) {
interFB = theFerm[ix].second->evaluate(q2, 3, offshell,
spbOut[ofh], vecIn[ivh]);
diag[ix] = theFerm[ix].first->evaluate(q2, spIn[ifh], interFB,
vecOut[ovh]);
}
else if(offshell->iSpin() == PDT::Spin1) {
interV = theVec[ix].second->evaluate(q2, 3, offshell,
vecIn[ivh], vecOut[ovh]);
diag[ix] = theVec[ix].first->evaluate(q2, spIn[ifh],
spbOut[ofh], interV);
}
else
diag[ix] = 0.0;
}
else if(diagram.channelType == HPDiagram::sChannel) {
interFB = theFerm[ix].second->evaluate(q2, 1, offshell,
spbOut[ofh], vecOut[ovh]);
diag[ix] = theFerm[ix].first->evaluate(q2, spIn[ifh], interFB,
vecIn[ivh]);
}
me[ix] += norm(diag[ix]);
//Compute flows
for(size_t iy = 0; iy < diagram.colourFlow.size(); ++iy)
cflows[diagram.colourFlow[iy].first - 1] +=
diagram.colourFlow[iy].second * diag[ix];
}//end of diag loop
//Now add flows to me2 with appropriate colour factors
for(size_t ii = 0; ii < ncf; ++ii)
for(size_t ij = 0; ij < ncf; ++ij)
mesq += cfactors[ii][ij]*(cflows[ii]*conj(cflows[ij])).real();
prodME(ifh, 2*ivh, ovh, ofh) =
std::accumulate(cflows.begin(), cflows.end(), Complex(0., 0.));
}
}
}
}
const double colfact = mePartonData()[0]->iColour() == PDT::Colour3 ?
1./24. : 1;
DVector save(ndiags);
for(DVector::size_type ix = 0; ix < ndiags; ++ix)
save[ix] = 0.25*colfact*me[ix];
meInfo(save);
mesq *= 0.25*colfact;
return prodME;
}
ProductionMatrixElement
MEfv2vf::fbv2vfbHeME(SpinorBarVector & spbIn, VBVector & vecIn,
VBVector & vecOut, SpinorVector & spOut,
double & mesq) const {
const HPCount ndiags = numberOfDiags();
const size_t ncf(numberOfFlows());
const vector<vector<double> > cfactors(getColourFactors());
const Energy2 q2(scale());
const bool masslessC = mePartonData()[2]->mass() == 0.0;
ProductionMatrixElement prodME(PDT::Spin1Half, PDT::Spin1, PDT::Spin1,
PDT::Spin1Half);
vector<Complex> diag(ndiags, Complex(0., 0.));
vector<double> me(ndiags, 0.);
//intermediate wave functions
SpinorBarWaveFunction interFB; VectorWaveFunction interV;
SpinorWaveFunction interF;
//loop over helicities
for(unsigned int ifh = 0; ifh < 2; ++ifh) {
for(unsigned int ivh = 0; ivh < 2; ++ivh) {
for(unsigned int ovh = 0; ovh < 3; ++ovh) {
if(masslessC && ovh == 1) ++ovh;
for(unsigned int ofh = 0; ofh < 2; ++ofh) {
//store the flows
vector<Complex> cflows(ncf, Complex(0. ,0.));
for(HPCount ix = 0; ix < ndiags; ++ix) {
HPDiagram diagram = getProcessInfo()[ix];
tcPDPtr offshell = diagram.intermediate;
if(diagram.channelType == HPDiagram::tChannel) {
if(offshell->iSpin() == PDT::Spin1Half) {
interFB = theFerm[ix].first->evaluate(q2, 3, offshell,
spbIn[ifh], vecOut[ovh]);
diag[ix] = theFerm[ix].second->evaluate(q2, spOut[ofh], interFB,
vecIn[ivh]);
}
else if(offshell->iSpin() == PDT::Spin1) {
interV = theVec[ix].first->evaluate(q2, 3, offshell,
spOut[ofh], spbIn[ifh]);
diag[ix] = theVec[ix].second->evaluate(q2, vecIn[ivh], interV,
vecOut[ovh]);
}
else diag[ix] = 0.0;
}
else if(diagram.channelType == HPDiagram::sChannel) {
if(offshell->iSpin() == PDT::Spin1Half) {
interFB = theFerm[ix].first->evaluate(q2, 1, offshell, spbIn[ifh],
vecIn[ivh]);
diag[ix] = theFerm[ix].second->evaluate(q2, spOut[ofh], interFB,
vecOut[ovh]);
}
}
me[ix] += norm(diag[ix]);
//Compute flows
for(size_t iy = 0; iy < diagram.colourFlow.size(); ++iy)
cflows[diagram.colourFlow[iy].first - 1] +=
diagram.colourFlow[iy].second * diag[ix];
}//end of diag loop
//Now add flows to me2 with appropriate colour factors
for(size_t ii = 0; ii < ncf; ++ii)
for(size_t ij = 0; ij < ncf; ++ij)
mesq += cfactors[ii][ij]*(cflows[ii]*conj(cflows[ij])).real();
prodME(ifh, ivh, ovh, ofh) =
std::accumulate(cflows.begin(), cflows.end(), Complex(0., 0.));
}
}
}
}
const double colfact = mePartonData()[0]->iColour() == PDT::Colour3bar ?
1./24. : 1;
DVector save(ndiags);
for(DVector::size_type ix = 0; ix < ndiags; ++ix)
save[ix] = 0.25*colfact*me[ix];
meInfo(save);
mesq *= 0.25*colfact;
return prodME;
}
Selector<const ColourLines *>
MEfv2vf::colourGeometries(tcDiagPtr diag) const {
static vector<ColourLines> cf(10);
//3 8->8 3
cf[0] = ColourLines("1 4, -4 2 -3, 3 5");
- cf[1] = ColourLines("1 2 -3, 4 -2 5, 3 4");
+ cf[1] = ColourLines("1 2 -3, -4 -2 5, 3 4");
cf[2] = ColourLines("1 -2, 2 3 4, -4 5");
//3b 8 -> 8 3b
- cf[3] = ColourLines("-4 -1, 3 -2 4, -5 -3");
- cf[4] = ColourLines("3 -2 -1, -5 2 4, -3 -4");
+ cf[3] = ColourLines("-4 -1, 3 2 4, -5 -3");
+ cf[4] = ColourLines("3 2 -1, -5 -2 4, -4 -3");
cf[5] = ColourLines("2 -1, -4 -3 -2, -5 4");
//3 8 -> 0 3
cf[6] = ColourLines("1 2 -3, 3 5");
cf[7] = ColourLines("1 -2, 2 3 5");
//3b 8 -> 0 3b
cf[8] = ColourLines("3 2 -1, -3 -5");
cf[9] = ColourLines("2 -1, -5 -3 -2");
HPDiagram current = getProcessInfo()[abs(diag->id()) - 1];
PDT::Colour offcolour = current.intermediate->iColour();
bool octetC = (mePartonData()[2]->iColour() == PDT::Colour8);
bool cc(mePartonData()[0]->id() < 0);
Selector<const ColourLines *> select;
int icf(-1);
if(current.channelType == HPDiagram::sChannel)
if(octetC) icf = cc ? 5 : 2;
else icf = cc ? 9 : 7;
else if(current.channelType == HPDiagram::tChannel) {
if(offcolour == PDT::Colour3 || offcolour == PDT::Colour3bar) {
if(octetC) icf = cc ? 3 : 0;
else icf = cc ? 8 : 6;
}
else if(offcolour == PDT::Colour8)
icf = cc ? 4 : 1;
else
throw MEException() << "MEfv2vf::colourGeometries - There is an incorrect "
<< "coloured object in a t-channel diagram. "
<< "No colour lines set."
<< Exception::warning;
}
else
throw MEException() << "MEfv2vf::colourGeometries - Incorrect diagram type "
<< "encountered. No colour lines set."
<< Exception::warning;
if(icf >= 0) select.insert(1., &cf[icf]);
return select;
}
void MEfv2vf::persistentOutput(PersistentOStream & os) const {
os << theFerm << theVec;
}
void MEfv2vf::persistentInput(PersistentIStream & is, int) {
is >> theFerm >> theVec;
}
ClassDescription<MEfv2vf> MEfv2vf::initMEfv2vf;
// Definition of the static class description member.
void MEfv2vf::Init() {
static ClassDocumentation<MEfv2vf> documentation
("This is the implementation of the matrix element for a fermion-vector boson"
"to a vector-fermion.");
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Sep 30, 5:44 AM (10 h, 10 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6548500
Default Alt Text
(37 KB)

Event Timeline

HEPForge.org · .