No OneTemporary
Actions

Size

50 KB

Subscribers

None

View Options

	diff --git a/Decay/General/FFSDecayer.cc b/Decay/General/FFSDecayer.cc
	--- a/Decay/General/FFSDecayer.cc
	+++ b/Decay/General/FFSDecayer.cc
	@@ -1,465 +1,472 @@
	// -- C++ --
	//
	// FFSDecayer.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 FFSDecayer class.
	//

	#include "FFSDecayer.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	#include "Herwig/Utilities/Kinematics.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"

	using namespace Herwig;
	using namespace ThePEG::Helicity;

	IBPtr FFSDecayer::clone() const {
	return new_ptr(*this);
	}

	IBPtr FFSDecayer::fullclone() const {
	return new_ptr(*this);
	}

	void FFSDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
	vector<VertexBasePtr> vertex,
	map<ShowerInteraction,VertexBasePtr> & inV,
	const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
	map<ShowerInteraction,VertexBasePtr> ) {
	decayInfo(incoming,outgoing);
	for(auto vert : vertex) {
	vertex_ .push_back(dynamic_ptr_cast<AbstractFFSVertexPtr>(vert));
	perturbativeVertex_ .push_back(dynamic_ptr_cast<FFSVertexPtr> (vert));
	}
	vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
	for(auto & inter : itemp) {
	incomingVertex_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(inV.at(inter));
	outgoingVertexF_[inter] = AbstractFFVVertexPtr();
	outgoingVertexS_[inter] = AbstractVSSVertexPtr();
	if(outV[0].at(inter)) {
	if (outV[0].at(inter)->getName()==VertexType::FFV)
	outgoingVertexF_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[0].at(inter));
	else
	outgoingVertexS_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(outV[0].at(inter));
	}
	if(outV[1].at(inter)) {
	if (outV[1].at(inter)->getName()==VertexType::FFV)
	outgoingVertexF_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[1].at(inter));
	else
	outgoingVertexS_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(outV[1].at(inter));
	}
	}
	}

	void FFSDecayer::persistentOutput(PersistentOStream & os) const {
	os << perturbativeVertex_ << vertex_
	<< incomingVertex_ << outgoingVertexF_
	<< outgoingVertexS_;
	}

	void FFSDecayer::persistentInput(PersistentIStream & is, int) {
	is >> perturbativeVertex_ >> vertex_
	>> incomingVertex_ >> outgoingVertexF_
	>> outgoingVertexS_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	-DescribeClass<FFSDecayer,GeneralTwoBodyDecayer>
	+DescribeClass<FFSDecayer,GeneralTwoBodyDecayer2>
	describeHerwigFFSDecayer("Herwig::FFSDecayer", "Herwig.so");

	void FFSDecayer::Init() {

	static ClassDocumentation<FFSDecayer> documentation
	("The FFSDecayer class implements the decay of a fermion to "
	"a fermion and a scalar.");

	}

	-double FFSDecayer::me2(const int , const Particle & inpart,
	- const ParticleVector & decay,
	+void FFSDecayer::
	+constructSpinInfo(const Particle & part, ParticleVector decay) const {
	+ int itype[2];
	+ if(part.dataPtr()->CC()) itype[0] = part.id() > 0? 0:1;
	+ else itype[0] = 2;
	+ if(decay[0]->dataPtr()->CC()) itype[1] = decay[0]->id() > 0? 0:1;
	+ else itype[1] = 2;
	+ bool ferm(itype[0] == 0 \|\| itype[1] == 0 \|\| (itype[0] == 2 && itype[1] == 2));
	+ // for the decaying particle
	+ if(ferm) {
	+ SpinorWaveFunction::
	+ constructSpinInfo(wave_,const_ptr_cast<tPPtr>(&part),incoming,true);
	+ SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true);
	+ }
	+ else {
	+ SpinorBarWaveFunction::
	+ constructSpinInfo(wavebar_,const_ptr_cast<tPPtr>(&part),incoming,true);
	+ SpinorWaveFunction::constructSpinInfo(wave_,decay[0],outgoing,true);
	+ }
	+ ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true);
	+}
	+
	+double FFSDecayer::me2(const int,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	MEOption meopt) const {
	if(!ME())
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin0)));
	//Need to use different barred or unbarred spinors depending on
	//whether particle is cc or not.
	int itype[2];
	- if(inpart.dataPtr()->CC()) itype[0] = inpart.id() > 0? 0:1;
	- else itype[0] = 2;
	- if(decay[0]->dataPtr()->CC()) itype[1] = decay[0]->id() > 0? 0:1;
	+ if(part.dataPtr()->CC()) itype[0] = part.id() > 0? 0:1;
	+ else itype[0] = 2;
	+ if(outgoing[0]->CC()) itype[1] = outgoing[0]->id() > 0? 0:1;
	else itype[1] = 2;
	bool ferm(itype[0] == 0 \|\| itype[1] == 0 \|\| (itype[0] == 2 && itype[1] == 2));

	if(meopt==Initialize) {
	// spinors and rho
	if(ferm) {
	SpinorWaveFunction ::calculateWaveFunctions(wave_,rho_,
	- const_ptr_cast<tPPtr>(&inpart),
	+ const_ptr_cast<tPPtr>(&part),
	incoming);
	if(wave_[0].wave().Type() != SpinorType::u)
	for(unsigned int ix = 0; ix < 2; ++ix) wave_ [ix].conjugate();
	}
	else {
	SpinorBarWaveFunction::calculateWaveFunctions(wavebar_,rho_,
	- const_ptr_cast<tPPtr>(&inpart),
	+ const_ptr_cast<tPPtr>(&part),
	incoming);
	if(wavebar_[0].wave().Type() != SpinorType::v)
	for(unsigned int ix = 0; ix < 2; ++ix) wavebar_[ix].conjugate();
	}
	// fix rho if no correlations
	fixRho(rho_);
	}
	- // setup spin info when needed
	- if(meopt==Terminate) {
	- // for the decaying particle
	- if(ferm) {
	- SpinorWaveFunction::
	- constructSpinInfo(wave_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
	- SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true);
	- }
	- else {
	- SpinorBarWaveFunction::
	- constructSpinInfo(wavebar_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
	- SpinorWaveFunction::constructSpinInfo(wave_,decay[0],outgoing,true);
	- }
	- ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true);
	- }
	if(ferm)
	SpinorBarWaveFunction::
	- calculateWaveFunctions(wavebar_,decay[0],outgoing);
	+ calculateWaveFunctions(wavebar_,momenta[0],outgoing[0],Helicity::outgoing);
	else
	SpinorWaveFunction::
	- calculateWaveFunctions(wave_ ,decay[0],outgoing);
	- ScalarWaveFunction scal(decay[1]->momentum(),decay[1]->dataPtr(),outgoing);
	- Energy2 scale(sqr(inpart.mass()));
	+ calculateWaveFunctions(wave_ ,momenta[0],outgoing[0],Helicity::outgoing);
	+ ScalarWaveFunction scal(momenta[1],outgoing[1],Helicity::outgoing);
	+ Energy2 scale(sqr(part.mass()));
	for(unsigned int if1 = 0; if1 < 2; ++if1) {
	for(unsigned int if2 = 0; if2 < 2; ++if2) {
	if(ferm) (*ME())(if1, if2, 0) = 0.;
	else (*ME())(if2, if1, 0) = 0.;
	for(auto vert : vertex_) {
	if(ferm) (*ME())(if1, if2, 0) +=
	vert->evaluate(scale,wave_[if1],wavebar_[if2],scal);
	else (*ME())(if2, if1, 0) +=
	vert->evaluate(scale,wave_[if1],wavebar_[if2],scal);
	}
	}
	}
	double output = (ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
	// colour and identical particle factors
	- output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),
	- decay[1]->dataPtr());
	+ output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
	// return the answer
	return output;
	}

	Energy FFSDecayer::partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const {
	if( inpart.second < outa.second + outb.second ) return ZERO;
	if(perturbativeVertex_.size()==1 &&
	perturbativeVertex_[0]) {
	double mu1(0.),mu2(0.);
	tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
	if(outa.first->iSpin() == PDT::Spin1Half) {
	mu1 = outa.second/inpart.second;
	mu2 = outb.second/inpart.second;
	perturbativeVertex_[0]->setCoupling(sqr(inpart.second), in, outa.first, outb.first);
	}
	else {
	mu1 = outb.second/inpart.second;
	mu2 = outa.second/inpart.second;
	perturbativeVertex_[0]->setCoupling(sqr(inpart.second), in, outb.first, outa.first);

	}
	double c2 = norm(perturbativeVertex_[0]->norm());
	Complex cl = perturbativeVertex_[0]->left();
	Complex cr = perturbativeVertex_[0]->right();
	double me2 = c2( (norm(cl) + norm(cr))(1. + sqr(mu1) - sqr(mu2))
	+ 2.mu1(conj(cl)cr + conj(cr)cl).real() );
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second, outa.second,
	outb.second);
	Energy output = me2*pcm/16./Constants::pi;
	// colour factor
	output *= colourFactor(inpart.first,outa.first,outb.first);
	// return the answer
	return output;
	}
	else {
	- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
	+ return GeneralTwoBodyDecayer2::partialWidth(inpart,outa,outb);
	}
	}


	double FFSDecayer::threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt) {
	int iscal (0), iferm (1), iglu (2);
	// get location of outgoing fermion/scalar
	if(decay[1]->dataPtr()->iSpin()==PDT::Spin0) swap(iscal,iferm);
	// work out whether inpart is a fermion or antifermion
	int itype[2];
	if(inpart.dataPtr()->CC()) itype[0] = inpart.id() > 0 ? 0 : 1;
	else itype[0] = 2;
	if(decay[iferm]->dataPtr()->CC()) itype[1] = decay[iferm]->id() > 0 ? 0 : 1;
	else itype[1] = 2;

	bool ferm(false);
	if(itype[0] == itype[1] ) {
	ferm = itype[0]==0 \|\| (itype[0]==2 && decay[iscal]->id() < 0);
	}
	else if(itype[0] == 2) {
	ferm = itype[1]==0;
	}
	else if(itype[1] == 2) {
	ferm = itype[0]==0;
	}
	else if((itype[0] == 1 && itype[1] == 0) \|\|
	(itype[0] == 0 && itype[1] == 1)) {
	if(abs(inpart.id())<=16) {
	ferm = itype[0]==0;
	}
	else if(abs(decay[iferm]->id())<=16) {
	ferm = itype[1]==0;
	}
	else {
	ferm = true;
	}
	}
	else
	assert(false);
	if(meopt==Initialize) {
	// create spinor (bar) for decaying particle
	if(ferm) {
	SpinorWaveFunction::calculateWaveFunctions(wave3_, rho3_, const_ptr_cast<tPPtr>(&inpart),
	incoming);
	if(wave3_[0].wave().Type() != SpinorType::u)
	for(unsigned int ix = 0; ix < 2; ++ix) wave3_[ix].conjugate();
	}
	else {
	SpinorBarWaveFunction::calculateWaveFunctions(wavebar3_,rho3_, const_ptr_cast<tPPtr>(&inpart),
	incoming);
	if(wavebar3_[0].wave().Type() != SpinorType::v)
	for(unsigned int ix = 0; ix < 2; ++ix) wavebar3_[ix].conjugate();
	}
	}
	// setup spin information when needed
	if(meopt==Terminate) {
	if(ferm) {
	SpinorWaveFunction::
	constructSpinInfo(wave3_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
	SpinorBarWaveFunction::constructSpinInfo(wavebar3_,decay[iferm],outgoing,true);
	}
	else {
	SpinorBarWaveFunction::
	constructSpinInfo(wavebar3_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
	SpinorWaveFunction::constructSpinInfo(wave3_,decay[iferm],outgoing,true);
	}
	ScalarWaveFunction::constructSpinInfo( decay[iscal],outgoing,true);
	VectorWaveFunction::constructSpinInfo(gluon_, decay[iglu ],outgoing,true,false);
	return 0.;
	}

	// calulate colour factors and number of colour flows
	unsigned int nflow;
	vector<DVector> cfactors = getColourFactors(inpart, decay, nflow);

	vector<GeneralDecayMEPtr> ME(nflow,new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half, PDT::Spin0,
	PDT::Spin1Half, PDT::Spin1)));
	// create wavefunctions
	if (ferm) SpinorBarWaveFunction::calculateWaveFunctions(wavebar3_, decay[iferm],outgoing);
	else SpinorWaveFunction:: calculateWaveFunctions(wave3_ , decay[iferm],outgoing);

	ScalarWaveFunction swave3_(decay[iscal]->momentum(), decay[iscal]->dataPtr(),outgoing);
	VectorWaveFunction::calculateWaveFunctions(gluon_, decay[iglu ],outgoing,true);

	// gauge invariance test
	#ifdef GAUGE_CHECK
	gluon_.clear();
	for(unsigned int ix=0;ix<3;++ix) {
	if(ix==1) gluon_.push_back(VectorWaveFunction());
	else {
	gluon_.push_back(VectorWaveFunction(decay[iglu ]->momentum(),decay[iglu ]->dataPtr(),10,
	outgoing));
	}
	}
	#endif

	if (! ((incomingVertex_[inter] && (outgoingVertexF_[inter] \|\| outgoingVertexS_[inter])) \|\|
	(outgoingVertexF_[inter] && outgoingVertexS_[inter])))
	throw Exception()
	<< "Invalid vertices for radiation in FFS decay in FFSDecayer::threeBodyME"
	<< Exception::runerror;


	// sort out colour flows
	int F(1), S(2);
	if (decay[iscal]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[iferm]->dataPtr()->iColour()==PDT::Colour8)
	swap(F,S);
	else if (decay[iferm]->dataPtr()->iColour()==PDT::Colour3bar &&
	decay[iscal]->dataPtr()->iColour()==PDT::Colour8)
	swap(F,S);


	Complex diag;
	Energy2 scale(sqr(inpart.mass()));

	- const GeneralTwoBodyDecayer::CFlow & colourFlow
	+ const GeneralTwoBodyDecayer2::CFlow & colourFlow
	= colourFlows(inpart, decay);
	double gs(0.);
	bool couplingSet(false);
	#ifdef GAUGE_CHECK
	double total=0.;
	#endif
	for(unsigned int ifi = 0; ifi < 2; ++ifi) {
	for(unsigned int ifo = 0; ifo < 2; ++ifo) {
	for(unsigned int ig = 0; ig < 2; ++ig) {
	// radiation from the incoming fermion
	if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(inpart.dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(incomingVertex_[inter]);
	if (ferm) {
	SpinorWaveFunction spinorInter =
	incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),wave3_[ifi],
	gluon_[2*ig],inpart.mass());

	assert(wave3_[ifi].particle()->id()==spinorInter.particle()->id());
	diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,spinorInter,wavebar3_[ifo],swave3_);
	}
	else {
	SpinorBarWaveFunction spinorBarInter =
	incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),wavebar3_[ifi],
	gluon_[2*ig],inpart.mass());

	assert(wavebar3_[ifi].particle()->id()==spinorBarInter.particle()->id());
	diag = 0.;
	for(auto vertex :vertex_)
	diag+= vertex->evaluate(scale,wave3_[ifo], spinorBarInter,swave3_);
	}
	if(!couplingSet) {
	gs = abs(incomingVertex_[inter]->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
	(*ME[colourFlow[0][ix].first])(ifi, 0, ifo, ig) +=
	colourFlow[0][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}

	// radiation from outgoing fermion
	if((decay[iferm]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[iferm]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertexF_[inter]);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[iferm]->dataPtr();
	if(off->CC()) off = off->CC();
	if (ferm) {
	SpinorBarWaveFunction spinorBarInter =
	outgoingVertexF_[inter]->evaluate(scale,3,off,wavebar3_[ifo],
	gluon_[2*ig],decay[iferm]->mass());

	assert(wavebar3_[ifo].particle()->id()==spinorBarInter.particle()->id());
	diag = 0.;
	for(auto vertex :vertex_)
	diag+= vertex->evaluate(scale,wave3_[ifi],spinorBarInter,swave3_);
	}
	else {
	SpinorWaveFunction spinorInter =
	outgoingVertexF_[inter]->evaluate(scale,3,off,wave3_[ifo],
	gluon_[2*ig],decay[iferm]->mass());

	assert(wave3_[ifo].particle()->id()==spinorInter.particle()->id());
	diag = 0.;
	for(auto vertex :vertex_)
	diag+= vertex->evaluate(scale,spinorInter,wavebar3_[ifi],swave3_);
	}
	if(!couplingSet) {
	gs = abs(outgoingVertexF_[inter]->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[F].size();++ix) {
	(*ME[colourFlow[F][ix].first])(ifi, 0, ifo, ig) +=
	colourFlow[F][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}

	// radiation from outgoing scalar
	if((decay[iscal]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[iscal]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertexS_[inter]);
	// ensure you get correct ougoing particle from first vertex
	tcPDPtr off = decay[iscal]->dataPtr();
	if(off->CC()) off = off->CC();
	ScalarWaveFunction scalarInter =
	outgoingVertexS_[inter]->evaluate(scale,3,off,gluon_[2*ig],
	swave3_,decay[iscal]->mass());

	assert(swave3_.particle()->id()==scalarInter.particle()->id());
	if (ferm){
	diag = 0.;
	for(auto vertex :vertex_)
	diag += vertex->evaluate(scale,wave3_[ifi],wavebar3_[ifo],scalarInter);
	}
	else {
	diag = 0.;
	for(auto vertex :vertex_)
	diag += vertex->evaluate(scale,wave3_[ifo],wavebar3_[ifi],scalarInter);
	}
	if(!couplingSet) {
	gs = abs(outgoingVertexS_[inter]->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[S].size();++ix) {
	(*ME[colourFlow[S][ix].first])(ifi, 0, ifo, ig) +=
	colourFlow[S][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	}
	}
	}

	// contract matrices
	double output=0.;
	for(unsigned int ix=0; ix<nflow; ++ix){
	for(unsigned int iy=0; iy<nflow; ++iy){
	output+=cfactors[ix][iy](ME[ix]->contract(ME[iy],rho3_)).real();
	}
	}
	// divide by alpha(S,EM)
	output=(4.Constants::pi)/sqr(gs);
	#ifdef GAUGE_CHECK
	double ratio = output/total;
	if(abs(ratio)>1e-20) {
	generator()->log() << "Test of gauge invariance in decay\n" << inpart << "\n";
	for(unsigned int ix=0;ix<decay.size();++ix)
	generator()->log() << *decay[ix] << "\n";
	generator()->log() << "Test of gauge invariance " << ratio << "\n";
	}
	#endif

	// return the answer
	return output;
	}
	diff --git a/Decay/General/FFSDecayer.h b/Decay/General/FFSDecayer.h
	--- a/Decay/General/FFSDecayer.h
	+++ b/Decay/General/FFSDecayer.h
	@@ -1,218 +1,227 @@
	// -- C++ --
	//
	// FFSDecayer.h 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.
	//
	#ifndef HERWIG_FFSDecayer_H
	#define HERWIG_FFSDecayer_H
	//
	// This is the declaration of the FFSDecayer class.
	//

	-#include "GeneralTwoBodyDecayer.h"
	+#include "GeneralTwoBodyDecayer2.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Helicity/Vertex/Scalar/FFSVertex.h"
	#include "ThePEG/Helicity/Vertex/Scalar/VSSVertex.h"
	#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"


	namespace Herwig {
	using namespace ThePEG;
	using Helicity::FFSVertexPtr;

	/** \ingroup Decay
	* The FFSDecayer class implements the decay of a fermion
	* to a fermion and a vector in a general model. It holds an FFVVertex
	* pointer that must be typecast from the VertexBase pointer held in
	- * GeneralTwoBodyDecayer. It implents the virtual functions me2() and
	+ * GeneralTwoBodyDecayer2. It implents the virtual functions me2() and
	* partialWidth().
	*
	- * @see GeneralTwoBodyDecayer
	+ * @see GeneralTwoBodyDecayer2
	*/
	-class FFSDecayer: public GeneralTwoBodyDecayer {
	+class FFSDecayer: public GeneralTwoBodyDecayer2 {

	public:

	/**
	* The default constructor.
	*/
	FFSDecayer() {}

	/** @name Virtual functions required by the Decayer class. */
	//@{
	/**
	* Return the matrix element squared for a given mode and phase-space channel.
	- * @param ichan The channel we are calculating the matrix element for.
	+ * @param ichan The channel we are calculating the matrix element for.
	* @param part The decaying Particle.
	- * @param decay The particles produced in the decay.
	+ * @param outgoing The particles produced in the decay
	+ * @param momenta The momenta of the particles produced in the decay
	* @param meopt Option for the calculation of the matrix element
	* @return The matrix element squared for the phase-space configuration.
	*/
	- virtual double me2(const int ichan, const Particle & part,
	- const ParticleVector & decay, MEOption meopt) const;
	+ double me2(const int ichan,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	+ MEOption meopt) const;
	+
	+ /**
	+ * Construct the SpinInfos for the particles produced in the decay
	+ */
	+ virtual void constructSpinInfo(const Particle & part,
	+ ParticleVector outgoing) const;

	/**
	* Function to return partial Width
	* @param inpart The decaying particle.
	* @param outa One of the decay products.
	* @param outb The other decay product.
	*/
	virtual Energy partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const;
	/**
	* Has a POWHEG style correction
	*/
	virtual POWHEGType hasPOWHEGCorrection() {
	POWHEGType output = FSR;
	for(auto vertex : vertex_) {
	if(vertex->orderInAllCouplings()!=1) {
	output = No;
	break;
	}
	}
	return output;
	}

	/**
	* Three-body matrix element including additional QCD radiation
	*/
	virtual double threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt);

	/**
	* Set the information on the decay
	*/
	virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
	vector<VertexBasePtr>,
	map<ShowerInteraction,VertexBasePtr> &,
	const vector<map<ShowerInteraction,VertexBasePtr> > &,
	map<ShowerInteraction,VertexBasePtr>);
	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	FFSDecayer & operator=(const FFSDecayer &);

	private:

	/**
	* Abstract pointer to AbstractFFSVertex
	*/
	vector<AbstractFFSVertexPtr> vertex_;

	/**
	* Pointer to the perturbative vertex
	*/
	vector<FFSVertexPtr> perturbativeVertex_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from incoming (anti)fermion
	*/
	map<ShowerInteraction,AbstractFFVVertexPtr> incomingVertex_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
	*/
	map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertexF_;

	/**
	* Abstract pointer to AbstractVSSVertex for QCD radiation from outgoing scalar
	*/
	map<ShowerInteraction,AbstractVSSVertexPtr> outgoingVertexS_;

	/**
	* Spin density matrix
	*/
	mutable RhoDMatrix rho_;

	/**
	* Spinor wavefunctions
	*/
	mutable vector<SpinorWaveFunction> wave_ ;

	/**
	* Barred spinor wavefunctions
	*/
	mutable vector<SpinorBarWaveFunction> wavebar_;

	/**
	* Spin density matrix for 3 body decay
	*/
	mutable RhoDMatrix rho3_;

	/**
	* Scalar wavefunction for 3 body decay
	*/
	mutable ScalarWaveFunction swave3_;

	/**
	* Spinor wavefunction for 3 body decay
	*/
	mutable vector<SpinorWaveFunction> wave3_;

	/**
	* Barred spinor wavefunction for 3 body decay
	*/
	mutable vector<SpinorBarWaveFunction> wavebar3_;

	/**
	* Vector wavefunction for 3 body decay
	*/
	mutable vector<VectorWaveFunction> gluon_;

	};

	}

	#endif /* HERWIG_FFSDecayer_H */
	diff --git a/Decay/General/FFVDecayer.cc b/Decay/General/FFVDecayer.cc
	--- a/Decay/General/FFVDecayer.cc
	+++ b/Decay/General/FFVDecayer.cc
	@@ -1,458 +1,469 @@
	// -- C++ --
	//
	// FFVDecayer.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 FFVDecayer class.
	//

	#include "FFVDecayer.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	#include "Herwig/Utilities/Kinematics.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"
	using namespace Herwig;
	using namespace ThePEG::Helicity;

	IBPtr FFVDecayer::clone() const {
	return new_ptr(*this);
	}

	IBPtr FFVDecayer::fullclone() const {
	return new_ptr(*this);
	}

	void FFVDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
	vector<VertexBasePtr> vertex,
	map<ShowerInteraction,VertexBasePtr> & inV,
	const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
	map<ShowerInteraction,VertexBasePtr> ) {
	decayInfo(incoming,outgoing);
	for(auto vert : vertex) {
	vertex_ .push_back(dynamic_ptr_cast<AbstractFFVVertexPtr>(vert));
	perturbativeVertex_ .push_back(dynamic_ptr_cast<FFVVertexPtr> (vert));
	}
	vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
	for(auto & inter : itemp) {
	incomingVertex_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(inV.at(inter));
	if(outV[0].at(inter)) {
	if (outV[0].at(inter)->getName()==VertexType::FFV)
	outgoingVertexF_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[0].at(inter));
	else
	outgoingVertexV_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[0].at(inter));
	}
	if(outV[1].at(inter)) {
	if (outV[1].at(inter)->getName()==VertexType::FFV)
	outgoingVertexF_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[1].at(inter));
	else
	outgoingVertexV_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[1].at(inter));
	}
	}
	}

	void FFVDecayer::persistentOutput(PersistentOStream & os) const {
	os << vertex_ << perturbativeVertex_
	<< incomingVertex_ << outgoingVertexF_
	<< outgoingVertexV_;
	}

	void FFVDecayer::persistentInput(PersistentIStream & is, int) {
	is >> vertex_ >> perturbativeVertex_
	>> incomingVertex_ >> outgoingVertexF_
	>> outgoingVertexV_;
	}

	-double FFVDecayer::me2(const int , const Particle & inpart,
	- const ParticleVector & decay,
	+void FFVDecayer::
	+constructSpinInfo(const Particle & part, ParticleVector decay) const {
	+ // type of process
	+ int itype[2];
	+ if(part.dataPtr()->CC()) itype[0] = part.id() > 0 ? 0 : 1;
	+ else itype[0] = 2;
	+ if(decay[0]->dataPtr()->CC()) itype[1] = decay[0]->id() > 0 ? 0 : 1;
	+ else itype[1] = 2;
	+ //Need to use different barred or unbarred spinors depending on
	+ //whether particle is cc or not.
	+ bool ferm(itype[0] == 0 \|\| itype[1] == 0 \|\| (itype[0] == 2 && itype[1] == 2));
	+ // for the decaying particle
	+ if(ferm) {
	+ SpinorWaveFunction::
	+ constructSpinInfo(wave_,const_ptr_cast<tPPtr>(&part),incoming,true);
	+ SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true);
	+ }
	+ else {
	+ SpinorBarWaveFunction::
	+ constructSpinInfo(wavebar_,const_ptr_cast<tPPtr>(&part),incoming,true);
	+ SpinorWaveFunction::constructSpinInfo(wave_,decay[0],outgoing,true);
	+ }
	+ VectorWaveFunction::
	+ constructSpinInfo(vector_,decay[1],outgoing,true,false);
	+}
	+
	+double FFVDecayer::me2(const int,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	MEOption meopt) const {
	if(!ME())
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin1)));
	// type of process
	int itype[2];
	- if(inpart.dataPtr()->CC()) itype[0] = inpart.id() > 0 ? 0 : 1;
	- else itype[0] = 2;
	- if(decay[0]->dataPtr()->CC()) itype[1] = decay[0]->id() > 0 ? 0 : 1;
	- else itype[1] = 2;
	+ if(part.dataPtr()->CC()) itype[0] = part.id() > 0 ? 0 : 1;
	+ else itype[0] = 2;
	+ if(outgoing[0]->CC()) itype[1] = outgoing[0]->id() > 0 ? 0 : 1;
	+ else itype[1] = 2;
	//Need to use different barred or unbarred spinors depending on
	//whether particle is cc or not.
	bool ferm(itype[0] == 0 \|\| itype[1] == 0 \|\| (itype[0] == 2 && itype[1] == 2));
	if(meopt==Initialize) {
	// spinors and rho
	if(ferm) {
	SpinorWaveFunction ::calculateWaveFunctions(wave_,rho_,
	- const_ptr_cast<tPPtr>(&inpart),
	+ const_ptr_cast<tPPtr>(&part),
	incoming);
	if(wave_[0].wave().Type() != SpinorType::u)
	for(unsigned int ix = 0; ix < 2; ++ix) wave_ [ix].conjugate();
	}
	else {
	SpinorBarWaveFunction::calculateWaveFunctions(wavebar_,rho_,
	- const_ptr_cast<tPPtr>(&inpart),
	+ const_ptr_cast<tPPtr>(&part),
	incoming);
	if(wavebar_[0].wave().Type() != SpinorType::v)
	for(unsigned int ix = 0; ix < 2; ++ix) wavebar_[ix].conjugate();
	}
	// fix rho if no correlations
	fixRho(rho_);
	}
	- // setup spin info when needed
	- if(meopt==Terminate) {
	- // for the decaying particle
	- if(ferm) {
	- SpinorWaveFunction::
	- constructSpinInfo(wave_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
	- SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true);
	- }
	- else {
	- SpinorBarWaveFunction::
	- constructSpinInfo(wavebar_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
	- SpinorWaveFunction::constructSpinInfo(wave_,decay[0],outgoing,true);
	- }
	- VectorWaveFunction::
	- constructSpinInfo(vector_,decay[1],outgoing,true,false);
	- }
	- Energy2 scale(sqr(inpart.mass()));
	+ Energy2 scale(sqr(part.mass()));
	if(ferm)
	SpinorBarWaveFunction::
	- calculateWaveFunctions(wavebar_,decay[0],outgoing);
	+ calculateWaveFunctions(wavebar_,momenta[0],outgoing[0],Helicity::outgoing);
	else
	SpinorWaveFunction::
	- calculateWaveFunctions(wave_ ,decay[0],outgoing);
	- bool massless = decay[1]->dataPtr()->mass()==ZERO;
	+ calculateWaveFunctions(wave_ ,momenta[0],outgoing[0],Helicity::outgoing);
	+ bool massless = outgoing[1]->mass()==ZERO;
	VectorWaveFunction::
	- calculateWaveFunctions(vector_,decay[1],outgoing,massless);
	+ calculateWaveFunctions(vector_,momenta[1],outgoing[1],Helicity::outgoing,massless);
	for(unsigned int if1 = 0; if1 < 2; ++if1) {
	for(unsigned int if2 = 0; if2 < 2; ++if2) {
	for(unsigned int vhel = 0; vhel < 3; ++vhel) {
	if(massless && vhel == 1) ++vhel;
	if(ferm)
	(*ME())(if1, if2,vhel) = 0.;
	else
	(*ME())(if2, if1, vhel) = 0.;
	for(auto vertex : vertex_) {
	if(ferm)
	(*ME())(if1, if2,vhel) +=
	vertex->evaluate(scale,wave_[if1],wavebar_[if2],vector_[vhel]);
	else
	(*ME())(if2, if1, vhel) +=
	vertex->evaluate(scale,wave_[if1],wavebar_[if2],vector_[vhel]);
	}
	}
	}
	}
	double output=(ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
	// colour and identical particle factors
	- output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),decay[1]->dataPtr());
	+ output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
	// return the answer
	return output;
	}

	Energy FFVDecayer::partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const {
	if( inpart.second < outa.second + outb.second ) return ZERO;
	if(perturbativeVertex_.size()==1 &&
	perturbativeVertex_[0]) {
	double mu1(outa.second/inpart.second),mu2(outb.second/inpart.second);
	tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
	if( outa.first->iSpin() == PDT::Spin1Half)
	perturbativeVertex_[0]->setCoupling(sqr(inpart.second), in,
	outa.first, outb.first);
	else {
	swap(mu1,mu2);
	perturbativeVertex_[0]->setCoupling(sqr(inpart.second),in,
	outb.first,outa.first);
	}
	Complex cl(perturbativeVertex_[0]->left()),cr(perturbativeVertex_[0]->right());
	double me2(0.);
	if( mu2 > 0. ) {
	me2 = (norm(cl) + norm(cr))(1. + sqr(mu1mu2) + sqr(mu2)
	- 2.sqr(mu1) - 2.sqr(mu2*mu2)
	+ sqr(mu1*mu1))
	- 6.mu1sqr(mu2)(conj(cl)cr + conj(cr)*cl).real();
	me2 /= sqr(mu2);
	}
	else
	me2 = 2.( (norm(cl) + norm(cr))(sqr(mu1) + 1.)
	- 4.mu1(conj(cl)cr + conj(cr)cl).real() );
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second, outa.second,
	outb.second);
	Energy output = norm(perturbativeVertex_[0]->norm())me2pcm/16./Constants::pi;
	// colour factor
	output *= colourFactor(inpart.first,outa.first,outb.first);
	// return the answer
	return output;
	}
	else {
	- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
	+ return GeneralTwoBodyDecayer2::partialWidth(inpart,outa,outb);
	}
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	-DescribeClass<FFVDecayer,GeneralTwoBodyDecayer>
	+DescribeClass<FFVDecayer,GeneralTwoBodyDecayer2>
	describeHerwigFFVDecayer("Herwig::FFVDecayer", "Herwig.so");

	void FFVDecayer::Init() {

	static ClassDocumentation<FFVDecayer> documentation
	("The FFVDecayer class implements the decay of a fermion to a fermion and a vector boson");

	}

	double FFVDecayer::threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt) {

	int iferm (0), ivect (1), iglu (2);
	// get location of outgoing lepton/vector
	if(decay[1]->dataPtr()->iSpin()==PDT::Spin1Half) swap(iferm,ivect);
	// work out whether inpart is a fermion or antifermion
	int itype[2];
	if(inpart.dataPtr()->CC()) itype[0] = inpart.id() > 0 ? 0 : 1;
	else itype[0] = 2;
	if(decay[iferm]->dataPtr()->CC()) itype[1] = decay[iferm]->id() > 0 ? 0 : 1;
	else itype[1] = 2;

	bool ferm(itype[0] == 0 \|\| itype[1] == 0 \|\|
	(itype[0] == 2 && itype[1] == 2 && decay[ivect]->id() < 0));

	// no emissions from massive vectors
	bool massless = decay[ivect]->dataPtr()->mass()==ZERO;

	if(meopt==Initialize) {
	// create spinor (bar) for decaying particle
	if(ferm) {
	SpinorWaveFunction::calculateWaveFunctions(wave3_, rho3_, const_ptr_cast<tPPtr>(&inpart),
	incoming);
	if(wave3_[0].wave().Type() != SpinorType::u)
	for(unsigned int ix = 0; ix < 2; ++ix) wave3_[ix].conjugate();
	}
	else {
	SpinorBarWaveFunction::calculateWaveFunctions(wavebar3_,rho3_, const_ptr_cast<tPPtr>(&inpart),
	incoming);
	if(wavebar3_[0].wave().Type() != SpinorType::v)
	for(unsigned int ix = 0; ix < 2; ++ix) wavebar3_[ix].conjugate();
	}
	}
	// setup spin information when needed
	if(meopt==Terminate) {
	if(ferm) {
	SpinorWaveFunction::
	constructSpinInfo(wave3_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
	SpinorBarWaveFunction::constructSpinInfo(wavebar3_,decay[iferm],outgoing,true);
	}
	else {
	SpinorBarWaveFunction::
	constructSpinInfo(wavebar3_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
	SpinorWaveFunction::constructSpinInfo(wave3_,decay[iferm],outgoing,true);
	}
	VectorWaveFunction::constructSpinInfo(vector3_, decay[ivect],outgoing,true,massless);
	VectorWaveFunction::constructSpinInfo(gluon_, decay[iglu ],outgoing,true,false);
	return 0.;
	}

	// calulate colour factors and number of colour flows
	unsigned int nflow;
	vector<DVector> cfactors = getColourFactors(inpart, decay, nflow);

	vector<GeneralDecayMEPtr> ME(nflow,new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half, PDT::Spin1Half,
	PDT::Spin1, PDT::Spin1)));

	// create wavefunctions
	if (ferm) SpinorBarWaveFunction::calculateWaveFunctions(wavebar3_, decay[iferm],outgoing);
	else SpinorWaveFunction:: calculateWaveFunctions(wave3_ , decay[iferm],outgoing);

	VectorWaveFunction::calculateWaveFunctions(vector3_, decay[ivect],outgoing,massless);
	VectorWaveFunction::calculateWaveFunctions(gluon_, decay[iglu ],outgoing,true );

	// gauge invariance test
	#ifdef GAUGE_CHECK
	gluon_.clear();
	for(unsigned int ix=0;ix<3;++ix) {
	if(ix==1) gluon_.push_back(VectorWaveFunction());
	else {
	gluon_.push_back(VectorWaveFunction(decay[iglu ]->momentum(),
	decay[iglu ]->dataPtr(),10,
	outgoing));
	}
	}
	#endif

	if (! ((incomingVertex_[inter] && (outgoingVertexF_[inter] \|\| outgoingVertexV_[inter])) \|\|
	(outgoingVertexF_[inter] && outgoingVertexV_[inter])))
	throw Exception()
	<< "Invalid vertices for QCD radiation in FFV decay in FFVDecayer::threeBodyME"
	<< Exception::runerror;


	// sort out colour flows
	int F(1), V(2);
	if (decay[iferm]->dataPtr()->iColour()==PDT::Colour3bar &&
	decay[ivect]->dataPtr()->iColour()==PDT::Colour8)
	swap(F,V);
	else if (decay[ivect]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[iferm]->dataPtr()->iColour()==PDT::Colour8)
	swap(F,V);

	Complex diag;
	Energy2 scale(sqr(inpart.mass()));

	- const GeneralTwoBodyDecayer::CFlow & colourFlow = colourFlows(inpart, decay);
	+ const GeneralTwoBodyDecayer2::CFlow & colourFlow = colourFlows(inpart, decay);
	double gs(0.);
	bool couplingSet(false);
	#ifdef GAUGE_CHECK
	double total=0.;
	#endif
	for(unsigned int ifi = 0; ifi < 2; ++ifi) {
	for(unsigned int ifo = 0; ifo < 2; ++ifo) {
	for(unsigned int iv = 0; iv < 3; ++iv) {
	for(unsigned int ig = 0; ig < 2; ++ig) {
	// radiation from the incoming fermion
	if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(inpart.dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(incomingVertex_[inter]);
	if (ferm){
	SpinorWaveFunction spinorInter =
	incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),wave3_[ifi],
	gluon_[2*ig],inpart.mass());

	assert(wave3_[ifi].particle()->id()==spinorInter.particle()->id());
	diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,spinorInter,wavebar3_[ifo],vector3_[iv]);
	}
	else {
	SpinorBarWaveFunction spinorBarInter =
	incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),wavebar3_[ifi],
	gluon_[2*ig],inpart.mass());

	assert(wavebar3_[ifi].particle()->id()==spinorBarInter.particle()->id());
	diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,wave3_[ifo], spinorBarInter,vector3_[iv]);
	}
	if(!couplingSet) {
	gs = abs(incomingVertex_[inter]->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
	(*ME[colourFlow[0][ix].first])(ifi, ifo, iv, ig) +=
	colourFlow[0][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	// radiation from outgoing fermion
	if((decay[iferm]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[iferm]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertexF_[inter]);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[iferm]->dataPtr();
	if(off->CC()) off = off->CC();
	if (ferm) {
	SpinorBarWaveFunction spinorBarInter =
	outgoingVertexF_[inter]->evaluate(scale,3,off,wavebar3_[ifo],
	gluon_[2*ig],decay[iferm]->mass());

	assert(wavebar3_[ifo].particle()->id()==spinorBarInter.particle()->id());
	diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,wave3_[ifi],spinorBarInter,vector3_[iv]);
	}
	else {
	SpinorWaveFunction spinorInter =
	outgoingVertexF_[inter]->evaluate(scale,3,off,wave3_[ifo],
	gluon_[2*ig],decay[iferm]->mass());

	assert(wave3_[ifo].particle()->id()==spinorInter.particle()->id());

	diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,spinorInter,wavebar3_[ifi],vector3_[iv]);
	}
	if(!couplingSet) {
	gs = abs(outgoingVertexF_[inter]->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[F].size();++ix) {
	(*ME[colourFlow[F][ix].first])(ifi, ifo, iv, ig) +=
	colourFlow[F][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}

	// radiation from outgoing vector
	if((decay[ivect]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[ivect]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertexV_[inter]);
	// ensure you get correct ougoing particle from first vertex
	tcPDPtr off = decay[ivect]->dataPtr();
	if(off->CC()) off = off->CC();
	VectorWaveFunction vectorInter =
	outgoingVertexV_[inter]->evaluate(scale,3,off,gluon_[2*ig],
	vector3_[iv],decay[ivect]->mass());

	assert(vector3_[iv].particle()->id()==vectorInter.particle()->id());
	if (ferm) {
	diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,wave3_[ifi],wavebar3_[ifo],vectorInter);
	}
	else {
	diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,wave3_[ifo],wavebar3_[ifi],vectorInter);
	}
	if(!couplingSet) {
	gs = abs(outgoingVertexV_[inter]->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[V].size();++ix) {
	(*ME[colourFlow[V][ix].first])(ifi, ifo, iv, ig) +=
	colourFlow[V][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	}
	if(massless) ++iv;
	}
	}
	}

	// contract matrices
	double output=0.;
	for(unsigned int ix=0; ix<nflow; ++ix){
	for(unsigned int iy=0; iy<nflow; ++iy){
	output+=cfactors[ix][iy](ME[ix]->contract(ME[iy],rho3_)).real();
	}
	}
	// divide by alpha(S,eM)
	output = (4.Constants::pi)/sqr(gs);
	#ifdef GAUGE_CHECK
	double ratio = output/total;
	if(abs(ratio)>1e-20) {
	generator()->log() << "Test of gauge invariance in decay\n" << inpart << "\n";
	for(unsigned int ix=0;ix<decay.size();++ix)
	generator()->log() << *decay[ix] << "\n";
	generator()->log() << "Test of gauge invariance " << ratio << "\n";
	}
	#endif
	// return the answer
	return output;
	}
	diff --git a/Decay/General/FFVDecayer.h b/Decay/General/FFVDecayer.h
	--- a/Decay/General/FFVDecayer.h
	+++ b/Decay/General/FFVDecayer.h
	@@ -1,223 +1,232 @@
	// -- C++ --
	//
	// FFVDecayer.h 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.
	//
	#ifndef HERWIG_FFVDecayer_H
	#define HERWIG_FFVDecayer_H
	//
	// This is the declaration of the FFVDecayer class.
	//

	-#include "GeneralTwoBodyDecayer.h"
	+#include "GeneralTwoBodyDecayer2.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
	#include "ThePEG/Helicity/Vertex/Vector/VVVVertex.h"

	namespace Herwig {
	using namespace ThePEG;
	using Helicity::FFVVertexPtr;

	/** \ingroup Decay
	* The FFVDecayer class implements the decay of a fermion
	* to a fermion and a vector in a general model. It holds an FFVVertex
	* pointer that must be typecast from the VertexBase pointer held in
	- * GeneralTwoBodyDecayer. It implents the virtual functions me2() and
	+ * GeneralTwoBodyDecayer2. It implents the virtual functions me2() and
	* partialWidth().
	*
	- * @see GeneralTwoBodyDecayer
	+ * @see GeneralTwoBodyDecayer2
	*/
	-class FFVDecayer: public GeneralTwoBodyDecayer {
	+class FFVDecayer: public GeneralTwoBodyDecayer2 {

	public:

	/**
	* The default constructor.
	*/
	FFVDecayer() {}

	public:

	/** @name Virtual functions required by the Decayer class. */
	//@{
	- /**
	- * Return the matrix element squared for a given mode and phase-space channel.
	- * @param ichan The channel we are calculating the matrix element for.
	+ /**
	+ * Return the matrix element squared for a given mode and phase-space channel.
	+ * @param ichan The channel we are calculating the matrix element for.
	* @param part The decaying Particle.
	- * @param decay The particles produced in the decay.
	- * @param meopt Option for the matrix element
	+ * @param outgoing The particles produced in the decay
	+ * @param momenta The momenta of the particles produced in the decay
	+ * @param meopt Option for the calculation of the matrix element
	* @return The matrix element squared for the phase-space configuration.
	*/
	- virtual double me2(const int ichan, const Particle & part,
	- const ParticleVector & decay, MEOption meopt) const;
	+ double me2(const int ichan,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	+ MEOption meopt) const;
	+
	+ /**
	+ * Construct the SpinInfos for the particles produced in the decay
	+ */
	+ virtual void constructSpinInfo(const Particle & part,
	+ ParticleVector outgoing) const;

	/**
	* Function to return partial Width
	* @param inpart The decaying particle.
	* @param outa One of the decay products.
	* @param outb The other decay product.
	*/
	virtual Energy partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const;

	/**
	* Has a POWHEG style correction
	*/
	virtual POWHEGType hasPOWHEGCorrection() {
	POWHEGType output = FSR;
	for(auto vertex : vertex_) {
	if(vertex->orderInAllCouplings()!=1) {
	output = No;
	break;
	}
	}
	return output;
	}

	/**
	* Three-body matrix element including additional QCD radiation
	*/
	virtual double threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt);

	/**
	* Set the information on the decay
	*/
	virtual void setDecayInfo(PDPtr incoming, PDPair outgoing, vector<VertexBasePtr>,
	map<ShowerInteraction,VertexBasePtr> &,
	const vector<map<ShowerInteraction,VertexBasePtr> > &,
	map<ShowerInteraction,VertexBasePtr>);
	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	FFVDecayer & operator=(const FFVDecayer &);

	private:

	/**
	* Abstract pointer to AbstractFFVVertex
	*/
	vector<AbstractFFVVertexPtr> vertex_;

	/**
	* Pointer to the perturbative vertex
	*/
	vector<FFVVertexPtr> perturbativeVertex_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from incoming (anti)fermion
	*/
	map<ShowerInteraction,AbstractFFVVertexPtr> incomingVertex_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
	*/
	map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertexF_;

	/**
	* Abstract pointer to AbstractVVVVertex for QCD radiation from outgoing vector
	*/
	map<ShowerInteraction,AbstractVVVVertexPtr> outgoingVertexV_;

	/**
	* Spin density matrix
	*/
	mutable RhoDMatrix rho_;

	/**
	* Spinor wavefunction
	*/
	mutable vector<SpinorWaveFunction> wave_ ;

	/**
	* Barred spinor wavefunction
	*/
	mutable vector<SpinorBarWaveFunction> wavebar_;

	/**
	* Polarization vectors
	*/
	mutable vector<VectorWaveFunction> vector_;

	/**
	* Spin density matrix for 3 body decay
	*/
	mutable RhoDMatrix rho3_;

	/**
	* vector wavefunction for 3 body decay
	*/
	mutable vector<VectorWaveFunction> vector3_;

	/**
	* Spinor wavefunction for 3 body decay
	*/
	mutable vector<SpinorWaveFunction> wave3_;

	/**
	* Barred spinor wavefunction for 3 body decay
	*/
	mutable vector<SpinorBarWaveFunction> wavebar3_;

	/**
	* Vector wavefunction for gluon in 3 body decay
	*/
	mutable vector<VectorWaveFunction> gluon_;

	};

	}

	#endif /* HERWIG_FFVDecayer_H */

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 1:54 PM (15 h, 35 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4018828
Default Alt Text: (50 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions