No OneTemporary
Actions

Size

36 KB

Subscribers

None

View Options

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

	#include "SSVDecayer.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 "Herwig/Utilities/Kinematics.h"
	#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"
	using namespace Herwig;
	using namespace ThePEG::Helicity;

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

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

	void SSVDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
	VertexBasePtr vertex,
	map<ShowerInteraction,VertexBasePtr> & inV,
	const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
	map<ShowerInteraction,VertexBasePtr> fourV) {
	decayInfo(incoming,outgoing);
	vertex_ = dynamic_ptr_cast<AbstractVSSVertexPtr>(vertex);
	perturbativeVertex_ = dynamic_ptr_cast<VSSVertexPtr> (vertex);
	vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
	for(auto & inter : itemp) {
	incomingVertex_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(inV.at(inter));
	fourPointVertex_[inter] = dynamic_ptr_cast<AbstractVVSSVertexPtr>(fourV.at(inter));
	outgoingVertexS_[inter] = AbstractVSSVertexPtr();
	outgoingVertexV_[inter] = AbstractVVVVertexPtr();
	if(outV[0].at(inter)) {
	if (outV[0].at(inter)->getName()==VertexType::VSS)
	outgoingVertexS_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(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::VSS)
	outgoingVertexS_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(outV[1].at(inter));
	else
	outgoingVertexV_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[1].at(inter));
	}
	}
	}

	void SSVDecayer::persistentOutput(PersistentOStream & os) const {
	os << vertex_ << perturbativeVertex_
	<< incomingVertex_ << outgoingVertexS_
	<< outgoingVertexV_ << fourPointVertex_;
	}

	void SSVDecayer::persistentInput(PersistentIStream & is, int) {
	is >> vertex_ >> perturbativeVertex_
	>> incomingVertex_ >> outgoingVertexS_
	>> outgoingVertexV_ >> fourPointVertex_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<SSVDecayer,GeneralTwoBodyDecayer>
	describeHerwigSSVDecayer("Herwig::SSVDecayer", "Herwig.so");

	void SSVDecayer::Init() {

	static ClassDocumentation<SSVDecayer> documentation
	("This implements the decay of a scalar to a vector and a scalar");

	}

	double SSVDecayer::me2(const int , const Particle & inpart,
	const ParticleVector & decay,
	MEOption meopt) const {
	unsigned int isc(0),ivec(1);
	if(decay[0]->dataPtr()->iSpin() != PDT::Spin0) swap(isc,ivec);
	if(!ME()) {
	if(ivec==1)
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin1)));
	else
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1,PDT::Spin0)));
	}
	if(meopt==Initialize) {
	ScalarWaveFunction::
	calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&inpart),incoming);
	swave_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),incoming);
	}
	if(meopt==Terminate) {
	ScalarWaveFunction::
	constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),incoming,true);
	ScalarWaveFunction::
	constructSpinInfo(decay[isc],outgoing,true);
	VectorWaveFunction::
	constructSpinInfo(vector_,decay[ivec],outgoing,true,false);
	}
	VectorWaveFunction::
	calculateWaveFunctions(vector_,decay[ivec],outgoing,false);
	ScalarWaveFunction sca(decay[isc]->momentum(),decay[isc]->dataPtr(),outgoing);
	Energy2 scale(sqr(inpart.mass()));
	//make sure decay matrix element is in the correct order
	double output(0.);
	if(ivec == 0) {
	for(unsigned int ix = 0; ix < 3; ++ix)
	(*ME())(0, ix, 0) = vertex_->evaluate(scale,vector_[ix],sca, swave_);
	}
	else {
	for(unsigned int ix = 0; ix < 3; ++ix)
	(*ME())(0, 0, ix) = vertex_->evaluate(scale,vector_[ix],sca,swave_);
	}
	output = (ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
	// colour and identical particle factors
	output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),
	decay[1]->dataPtr());
	// return the answer
	return output;
	}

	Energy SSVDecayer:: partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const {
	if( inpart.second < outa.second + outb.second ) return ZERO;
	if(perturbativeVertex_) {
	double mu1sq(sqr(outa.second/inpart.second)),
	mu2sq(sqr(outb.second/inpart.second));
	tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
	if(outa.first->iSpin() == PDT::Spin0) {
	perturbativeVertex_->setCoupling(sqr(inpart.second), outb.first, outa.first,in);
	}
	else {
	swap(mu1sq,mu2sq);
	perturbativeVertex_->setCoupling(sqr(inpart.second), outa.first, outb.first,in);
	}
	double me2(0.);
	if(mu2sq == 0.)
	me2 = -2.*mu1sq - 2.;
	else
	me2 = ( sqr(mu2sq - mu1sq) - 2.*(mu2sq + mu1sq) + 1. )/mu2sq;
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second, outa.second,
	outb.second);
	Energy output = pcmme2norm(perturbativeVertex_->norm())/8./Constants::pi;
	// colour factor
	output *= colourFactor(inpart.first,outa.first,outb.first);
	// return the answer
	return output;
	}
	else {
	return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
	}
	}


	double SSVDecayer::threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt) {
	int iscal (0), ivect (1), iglu (2);
	// get location of outgoing scalar/vector
	if(decay[1]->dataPtr()->iSpin()==PDT::Spin0) swap(iscal,ivect);

	if(meopt==Initialize) {
	// create scalar wavefunction for decaying particle
	ScalarWaveFunction::calculateWaveFunctions(rho3_,const_ptr_cast<tPPtr>(&inpart),incoming);
	swave3_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),incoming);
	}
	// setup spin information when needed
	if(meopt==Terminate) {
	ScalarWaveFunction::
	constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),incoming,true);
	ScalarWaveFunction::
	constructSpinInfo(decay[iscal],outgoing,true);
	VectorWaveFunction::
	constructSpinInfo(vector3_,decay[ivect],outgoing,true,false);
	VectorWaveFunction::
	constructSpinInfo(gluon_, decay[iglu ],outgoing,true,false);
	return 0.;
	}
	// calculate 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::Spin0, PDT::Spin0,
	PDT::Spin1, PDT::Spin1)));

	// create wavefunctions
	- ScalarWaveFunction scal_(decay[iscal]->momentum(), decay[iscal]->dataPtr(),outgoing);
	+ ScalarWaveFunction scal(decay[iscal]->momentum(), decay[iscal]->dataPtr(),outgoing);
	VectorWaveFunction::calculateWaveFunctions(vector3_,decay[ivect],outgoing,false);
	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] && (outgoingVertexS_[inter] \|\| outgoingVertexV_[inter])) \|\|
	(outgoingVertexS_[inter] && outgoingVertexV_[inter])))
	throw Exception()
	<< "Invalid vertices for radiation in SSV decay in SSVDecayer::threeBodyME"
	<< Exception::runerror;


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

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

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

	assert(swave3_.particle()->id()==scalarInter.particle()->id());
	- Complex diag = vertex_->evaluate(scale,vector3_[iv],scal_,scalarInter);
	+ Complex diag = vertex_->evaluate(scale,vector3_[iv],scal,scalarInter);
	if(!couplingSet) {
	gs = abs(incomingVertex_[inter]->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
	(*ME[colourFlow[0][ix].first])(0, 0, iv, ig) +=
	colourFlow[0][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	// radiation from the outgoing scalar
	if((decay[iscal]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[iscal]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertexS_[inter]);
	// ensure you get correct outgoing 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],scal_,decay[iscal]->mass());
	+ outgoingVertexS_[inter]->evaluate(scale,3,off,gluon_[2*ig],scal,decay[iscal]->mass());

	- assert(scal_.particle()->id()==scalarInter.particle()->id());
	+ assert(scal.particle()->id()==scalarInter.particle()->id());

	if(!couplingSet) {
	gs = abs(outgoingVertexS_[inter]->norm());
	couplingSet = true;
	}
	Complex diag = vertex_->evaluate(scale,vector3_[iv],scalarInter,swave3_);
	for(unsigned int ix=0;ix<colourFlow[S].size();++ix) {
	(*ME[colourFlow[S][ix].first])(0, 0, iv, ig) +=
	colourFlow[S][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 outgoing 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(!couplingSet) {
	gs = abs(outgoingVertexV_[inter]->norm());
	couplingSet = true;
	}
	- Complex diag = vertex_->evaluate(scale,vectorInter,scal_,swave3_);
	+ Complex diag = vertex_->evaluate(scale,vectorInter,scal,swave3_);
	for(unsigned int ix=0;ix<colourFlow[V].size();++ix) {
	(*ME[colourFlow[V][ix].first])(0, 0, iv, ig) +=
	colourFlow[V][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	// radiation from 4 point vertex
	if (fourPointVertex_[inter]) {
	Complex diag = fourPointVertex_[inter]->evaluate(scale, gluon_[2*ig], vector3_[iv],
	- scal_, swave3_);
	+ scal, swave3_);
	for(unsigned int ix=0;ix<colourFlow[3].size();++ix) {
	(*ME[colourFlow[3][ix].first])(0, 0, iv, ig) +=
	colourFlow[3][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/VVSDecayer.cc b/Decay/General/VVSDecayer.cc
	--- a/Decay/General/VVSDecayer.cc
	+++ b/Decay/General/VVSDecayer.cc
	@@ -1,126 +1,302 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the VVSDecayer class.
	//

	#include "VVSDecayer.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"

	using namespace Herwig;
	using namespace ThePEG::Helicity;
	-
	IBPtr VVSDecayer::clone() const {
	return new_ptr(*this);
	}

	IBPtr VVSDecayer::fullclone() const {
	return new_ptr(*this);
	}
	void VVSDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
	VertexBasePtr vertex,
	- map<ShowerInteraction,VertexBasePtr> &,
	- const vector<map<ShowerInteraction,VertexBasePtr> > &,
	+ map<ShowerInteraction,VertexBasePtr> & inV,
	+ const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
	map<ShowerInteraction,VertexBasePtr>) {
	decayInfo(incoming,outgoing);
	vertex_ = dynamic_ptr_cast<AbstractVVSVertexPtr>(vertex);
	perturbativeVertex_ = dynamic_ptr_cast<VVSVertexPtr> (vertex);
	+ vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
	+ for(auto & inter : itemp) {
	+ incomingVertex_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(inV.at(inter));
	+ outgoingVertexS_[inter] = AbstractVSSVertexPtr();
	+ outgoingVertexV_[inter] = AbstractVVVVertexPtr();
	+ if(outV[0].at(inter)) {
	+ if (outV[0].at(inter)->getName()==VertexType::VSS)
	+ outgoingVertexS_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(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::VSS)
	+ outgoingVertexS_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(outV[1].at(inter));
	+ else
	+ outgoingVertexV_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[1].at(inter));
	+ }
	+ }
	}

	void VVSDecayer::persistentOutput(PersistentOStream & os) const {
	- os << vertex_ << perturbativeVertex_;
	+ os << vertex_ << perturbativeVertex_
	+ << incomingVertex_ << outgoingVertexS_ << outgoingVertexV_;
	}

	void VVSDecayer::persistentInput(PersistentIStream & is, int) {
	- is >> vertex_ >> perturbativeVertex_;
	+ is >> vertex_ >> perturbativeVertex_
	+ >> incomingVertex_ >> outgoingVertexS_ >> outgoingVertexV_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<VVSDecayer,GeneralTwoBodyDecayer>
	describeHerwigVVSDecayer("Herwig::VVSDecayer", "Herwig.so");

	void VVSDecayer::Init() {

	static ClassDocumentation<VVSDecayer> documentation
	("The VVSDecayer class implements the decay of a vector"
	" to a vector and a scalar");

	}

	double VVSDecayer::me2(const int , const Particle & inpart,
	const ParticleVector & decay,
	MEOption meopt) const {
	bool massless = ( decay[0]->id()==ParticleID::gamma \|\|
	decay[0]->id()==ParticleID::g );
	if(!ME())
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin1,PDT::Spin0)));
	if(meopt==Initialize) {
	VectorWaveFunction::calculateWaveFunctions(vectors_[0],rho_,
	const_ptr_cast<tPPtr>(&inpart),
	incoming,false);
	}
	if(meopt==Terminate) {
	VectorWaveFunction::constructSpinInfo(vectors_[0],const_ptr_cast<tPPtr>(&inpart),
	incoming,true,false);
	VectorWaveFunction::
	constructSpinInfo(vectors_[1],decay[0],outgoing,true,massless);
	ScalarWaveFunction::
	constructSpinInfo(decay[1],outgoing,true);
	return 0.;
	}
	VectorWaveFunction::
	calculateWaveFunctions(vectors_[1],decay[0],outgoing,massless);
	ScalarWaveFunction sca(decay[1]->momentum(),decay[1]->dataPtr(),outgoing);
	Energy2 scale(sqr(inpart.mass()));
	for(unsigned int in=0;in<3;++in) {
	for(unsigned int out=0;out<3;++out) {
	if(massless&&out==1) ++out;
	(*ME())(in,out,0) =
	vertex_->evaluate(scale,vectors_[0][in],vectors_[1][out],sca);
	}
	}
	double output=(ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
	// colour and identical particle factors
	output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),
	decay[1]->dataPtr());
	// return the answer
	return output;
	}

	Energy VVSDecayer::partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const {
	if( inpart.second < outa.second + outb.second ) return ZERO;
	if(perturbativeVertex_) {
	Energy2 scale(sqr(inpart.second));
	double mu1sq = sqr(outa.second/inpart.second);
	double mu2sq = sqr(outb.second/inpart.second);
	tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
	if( outb.first->iSpin() == PDT::Spin0 )
	perturbativeVertex_->setCoupling(sqr(inpart.second), in,
	outa.first, outb.first);
	else {
	perturbativeVertex_->setCoupling(sqr(inpart.second), in,
	outb.first, outa.first);
	swap(mu1sq, mu2sq);
	}
	double vn = norm(perturbativeVertex_->norm());
	if(vn == ZERO \|\| mu1sq == ZERO) return ZERO;
	double me2 = 2. + 0.25*sqr(1. + mu1sq - mu2sq)/mu1sq;
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,outa.second,
	outb.second);
	Energy output = vnme2pcm/(24.Constants::pi)/scaleUnitRemoval::E2;
	// colour factor
	output *= colourFactor(inpart.first,outa.first,outb.first);
	// return the answer
	return output;
	}
	else {
	return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
	}
	}

	+double VVSDecayer::threeBodyME(const int , const Particle & inpart,
	+ const ParticleVector & decay,
	+ ShowerInteraction inter, MEOption meopt) {
	+ unsigned int ivec(0),isca(1);
	+ if(decay[ivec]->dataPtr()->iSpin()!=PDT::Spin1) swap(ivec,isca);
	+ if(meopt==Initialize) {
	+ // create vector wavefunction for decaying particle
	+ VectorWaveFunction::calculateWaveFunctions(vectors3_[0], rho3_,
	+ const_ptr_cast<tPPtr>(&inpart),
	+ incoming, false);
	+ }
	+ if(meopt==Terminate) {
	+ VectorWaveFunction::
	+ constructSpinInfo(vectors3_[0] ,const_ptr_cast<tPPtr>(&inpart),outgoing,true,false);
	+ VectorWaveFunction::
	+ constructSpinInfo(vectors3_[1],decay[ivec],outgoing,true,false);
	+ ScalarWaveFunction::
	+ constructSpinInfo(decay[isca],outgoing,true);
	+ VectorWaveFunction::
	+ constructSpinInfo(gluon_ ,decay[2],outgoing,true,false);
	+ return 0.;
	+ }
	+ // calculate 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::Spin1, PDT::Spin1,
	+ PDT::Spin0, PDT::Spin1)));
	+ bool massless= decay[ivec]->mass()!=ZERO;
	+ // create wavefunctions
	+ VectorWaveFunction::calculateWaveFunctions(vectors3_[1],decay[0],outgoing,massless);
	+ ScalarWaveFunction scal(decay[isca]->momentum(), decay[isca]->dataPtr(),outgoing);
	+ VectorWaveFunction::calculateWaveFunctions(gluon_ ,decay[2],outgoing,true);
	+
	+ // gauge 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[2]->momentum(),
	+ decay[2]->dataPtr(),10,
	+ outgoing));
	+ }
	+ }
	+#endif
	+
	+ Energy2 scale(sqr(inpart.mass()));
	+
	+ const GeneralTwoBodyDecayer::CFlow & colourFlow
	+ = colourFlows(inpart, decay);
	+ double gs(0.);
	+ bool couplingSet(false);
	+#ifdef GAUGE_CHECK
	+ double total=0.;
	+#endif
	+
	+ for(unsigned int iv0 = 0; iv0 < 3; ++iv0) {
	+ for(unsigned int iv1 = 0; iv1 < 3; ++iv1) {
	+ if(massless && iv1==1) continue;
	+ for(unsigned int ig = 0; ig < 2; ++ig) {
	+ // radiation from the incoming vector
	+ if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	+ (inpart.dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	+ assert(incomingVertex_[inter]);
	+ VectorWaveFunction vectorInter =
	+ incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),vectors3_[0][iv0],
	+ gluon_[2*ig],inpart.mass());
	+
	+ assert(vectors3_[0][iv0].particle()->id()==vectorInter.particle()->id());
	+ Complex diag = vertex_->evaluate(scale,vectorInter,vectors3_[1][iv1],scal);
	+ if(!couplingSet) {
	+ gs = abs(incomingVertex_[inter]->norm());
	+ couplingSet = true;
	+ }
	+ for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
	+ (*ME[colourFlow[0][ix].first])(iv0, iv1, 0, ig) +=
	+ colourFlow[0][ix].second*diag;
	+ }
	+#ifdef GAUGE_CHECK
	+ total+=norm(diag);
	+#endif
	+ }
	+ // radiation from the outgoing vector
	+ if((decay[ivec]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	+ (decay[ivec]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	+ assert(outgoingVertexV_[inter]);
	+ // ensure you get correct outgoing particle from first vertex
	+ tcPDPtr off = decay[ivec]->dataPtr();
	+ if(off->CC()) off = off->CC();
	+ VectorWaveFunction vectorInter =
	+ outgoingVertexV_[inter]->evaluate(scale,3,off,gluon_[2*ig],vectors3_[1][iv1],decay[ivec]->mass());
	+
	+ assert(vectors3_[1][iv1].particle()->id()==vectorInter.particle()->id());
	+
	+ Complex diag =vertex_->evaluate(scale,vectors3_[0][iv0],vectorInter,scal);
	+ if(!couplingSet) {
	+ gs = abs(outgoingVertexV_[inter]->norm());
	+ couplingSet = true;
	+ }
	+ for(unsigned int ix=0;ix<colourFlow[1].size();++ix) {
	+ (*ME[colourFlow[1][ix].first])(iv0, iv1, 0, ig) +=
	+ colourFlow[1][ix].second*diag;
	+ }
	+#ifdef GAUGE_CHECK
	+ total+=norm(diag);
	+#endif
	+ }
	+ // radiation from the outgoing scalar
	+ if((decay[isca]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	+ (decay[isca]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	+ assert(outgoingVertexS_[inter]);
	+ // ensure you get correct outgoing particle from first vertex
	+ tcPDPtr off = decay[isca]->dataPtr();
	+ if(off->CC()) off = off->CC();
	+ ScalarWaveFunction scalarInter =
	+ outgoingVertexS_[inter]->evaluate(scale,3,off,gluon_[2*ig],scal,decay[isca]->mass());
	+
	+ assert(scal.particle()->id()==scalarInter.particle()->id());
	+
	+ Complex diag = vertex_->evaluate(scale,vectors3_[0][iv0],vectors3_[1][iv1],scalarInter);
	+ if(!couplingSet) {
	+ gs = abs(outgoingVertexS_[inter]->norm());
	+ couplingSet = true;
	+ }
	+ for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
	+ (*ME[colourFlow[2][ix].first])(iv0, iv1, 0, ig) +=
	+ colourFlow[2][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/VVSDecayer.h b/Decay/General/VVSDecayer.h
	--- a/Decay/General/VVSDecayer.h
	+++ b/Decay/General/VVSDecayer.h
	@@ -1,142 +1,193 @@
	// -- C++ --
	#ifndef THEPEG_VVSDecayer_H
	#define THEPEG_VVSDecayer_H
	//
	// This is the declaration of the VVSDecayer class.
	//

	#include "GeneralTwoBodyDecayer.h"
	#include "ThePEG/Helicity/Vertex/Scalar/VVSVertex.h"
	#include "ThePEG/Repository/EventGenerator.h"

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

	/** \ingroup Decay
	* The VVSDecayer class implements the decay of a vector to a
	* vector and a scalar in a general model. It holds an VVSVertex pointer
	* that must be typecast from the VertexBase pointer helid in the
	* GeneralTwoBodyDecayer. It implents the virtual functions me2() and
	* partialWidth().
	*
	* @see \ref VVSDecayerInterfaces "The interfaces"
	* defined for VVSDecayer.
	*/
	class VVSDecayer: public GeneralTwoBodyDecayer {

	public:

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

	/** @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 part The decaying Particle.
	* @param decay 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;

	/**
	* 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;

	/**
	* Set the information on the decay
	*/
	virtual void setDecayInfo(PDPtr incoming, PDPair outgoing, VertexBasePtr,
	map<ShowerInteraction,VertexBasePtr> &,
	const vector<map<ShowerInteraction,VertexBasePtr> > &,
	map<ShowerInteraction,VertexBasePtr>);
	+
	+ /**
	+ * Has a POWHEG style correction
	+ */
	+ virtual POWHEGType hasPOWHEGCorrection() {
	+ return (vertex_->orderInGem()+vertex_->orderInGs())==1 ? FSR : No;
	+ }
	+
	+ /**
	+ * Three-body matrix element including additional QCD radiation
	+ */
	+ virtual double threeBodyME(const int , const Particle & inpart,
	+ const ParticleVector & decay,
	+ ShowerInteraction inter, MEOption meopt);
	//@}

	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.
	*/
	VVSDecayer & operator=(const VVSDecayer &);

	private:

	/**
	* Abstract pointer to AbstractVVSVertex
	*/
	AbstractVVSVertexPtr vertex_;

	/**
	* Pointer to the perturbative vertex
	*/
	VVSVertexPtr perturbativeVertex_;

	/**
	+ * Abstract pointer to AbstractVVVVertex for QCD radiation from incoming vector
	+ */
	+ map<ShowerInteraction,AbstractVVVVertexPtr> incomingVertex_;
	+
	+ /**
	+ * Abstract pointer to AbstractVVVVertex for QCD radiation from the first outgoing vector
	+ */
	+ map<ShowerInteraction,AbstractVVVVertexPtr> outgoingVertexV_;
	+
	+ /**
	+ * Abstract pointer to AbstractVVVVertex for QCD radiation from the second outgoing vector
	+ */
	+ map<ShowerInteraction,AbstractVSSVertexPtr> outgoingVertexS_;
	+
	+ /**
	* Spin density matrix
	*/
	mutable RhoDMatrix rho_;

	/**
	* Vector wavefunctions
	*/
	mutable vector<Helicity::VectorWaveFunction> vectors_[2];
	+
	+private:
	+
	+ /**
	+ * Members for the POWHEG correction
	+ */
	+ //@{
	+ /**
	+ * Spin density matrix for 3 body decay
	+ */
	+ mutable RhoDMatrix rho3_;
	+
	+ /**
	+ * Vector wavefunctions
	+ */
	+ mutable vector<Helicity::VectorWaveFunction> vectors3_[2];
	+
	+ /**
	+ * Vector wavefunction for 3 body decay
	+ */
	+ mutable vector<Helicity::VectorWaveFunction> gluon_;
	+ //@}
	};

	}

	#endif /* THEPEG_VVSDecayer_H */
	diff --git a/Decay/General/VVVDecayer.h b/Decay/General/VVVDecayer.h
	--- a/Decay/General/VVVDecayer.h
	+++ b/Decay/General/VVVDecayer.h
	@@ -1,220 +1,220 @@
	// -- C++ --
	//
	// VVVDecayer.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_VVVDecayer_H
	#define HERWIG_VVVDecayer_H
	//
	// This is the declaration of the VVVDecayer class.
	//

	#include "GeneralTwoBodyDecayer.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Helicity/Vertex/Vector/VVVVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractVVVVVertex.h"

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

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

	public:

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

	/** @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 part The decaying Particle.
	* @param decay 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;

	/**
	* 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;

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

	/**
	* Has a POWHEG style correction
	*/
	virtual POWHEGType hasPOWHEGCorrection() {
	return (vertex_->orderInGem()+vertex_->orderInGs())==1 ? FSR : No;
	}

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

	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;
	//@}

	protected:

	/**
	* Find the vertices for the decay
	*/
	void identifyVertices(const Particle & inpart, const ParticleVector & decay,
	AbstractVVVVertexPtr & outgoingVertex1,
	AbstractVVVVertexPtr & outgoingVertex2,
	ShowerInteraction inter);

	private:

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

	private:

	/**
	* Abstract pointer to AbstractVVVVertex
	*/
	AbstractVVVVertexPtr vertex_;

	/**
	* Pointer to the perturbative vertex
	*/
	VVVVertexPtr perturbativeVertex_;

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

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

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

	/**
	* Abstract pointer to AbstractVVVVertex for QCD radiation from the 4-point vertex
	*/
	map<ShowerInteraction,AbstractVVVVVertexPtr> fourPointVertex_;

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

	/**
	* Vector wavefunctions
	*/
	mutable vector<Helicity::VectorWaveFunction> vectors_[3];

	private:

	/**
	- * Member for the POWHEG correction
	+ * Members for the POWHEG correction
	*/
	//@{
	/**
	* Spin density matrix for 3 body decay
	*/
	mutable RhoDMatrix rho3_;

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

	/**
	* Vector wavefunctions
	*/
	mutable vector<Helicity::VectorWaveFunction> vectors3_[2];

	/**
	* Vector wavefunction for 3 body decay
	*/
	mutable vector<Helicity::VectorWaveFunction> gluon_;
	//@}
	};

	}

	#endif /* HERWIG_VVVDecayer_H */

File Metadata

Mime Type: text/x-diff
Expires: Thu, Apr 24, 6:40 AM (1 d, 19 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4880448
Default Alt Text: (36 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions