No OneTemporary
Actions

Size

49 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/General/MEvv2ff.h b/MatrixElement/General/MEvv2ff.h
	--- a/MatrixElement/General/MEvv2ff.h
	+++ b/MatrixElement/General/MEvv2ff.h
	@@ -1,193 +1,191 @@
	// -- C++ --
	//
	// MEvv2ff.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_MEvv2ff_H
	#define HERWIG_MEvv2ff_H
	//
	// This is the declaration of the MEvv2ff class.
	//

	#include "GeneralHardME.h"
	#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractVVTVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractFFVVVertex.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	#include "Herwig/MatrixElement/ProductionMatrixElement.h"

	namespace Herwig {
	using namespace ThePEG;
	using ThePEG::Helicity::SpinorWaveFunction;
	using ThePEG::Helicity::SpinorBarWaveFunction;
	using ThePEG::Helicity::VectorWaveFunction;

	/**
	* This class is designed to implement the matrix element for the
	* \f$2 \rightarrow 2\f$ process vector-vector to fermion-antifermion pair. It
	* inherits from GeneralHardME and implements the me2() virtual function.
	*
	- * @see \ref MEvv2ffInterfaces "The Interfaces"
	- * defined for MEvv2ff.
	* @see GeneralHardME
	*
	*/
	class MEvv2ff: public GeneralHardME {

	public:

	/** A Vector of VectorWaveFunction objects. */
	typedef vector<VectorWaveFunction> VBVector;

	/** A vector of SpinorBarWaveFunction objects. */
	typedef vector<SpinorWaveFunction> SpinorVector;

	/** A vector of SpinorBarWaveFunction objects. */
	typedef vector<SpinorBarWaveFunction> SpinorBarVector;

	public:

	/** @name Virtual functions required by the MEBase class. */
	//@{
	/**
	* The matrix element for the kinematical configuration
	* previously provided by the last call to setKinematics(), suitably
	* scaled by sHat() to give a dimension-less number.
	* @return the matrix element scaled with sHat() to give a
	* dimensionless number.
	*/
	virtual double me2() const;
	//@}

	/**
	* Construct the vertex information for the spin correlations
	* @param sub Pointer to the relevent SubProcess
	*/
	virtual void constructVertex(tSubProPtr sub);

	private:

	/**
	* Calculate the value of the matrix element
	*/
	ProductionMatrixElement vv2ffME(const VBVector & v1, const VBVector & v2,
	const SpinorBarVector & sbar,
	const SpinorVector & sp,
	double & me2, bool first) const;

	protected:

	/**
	* A debugging function to test the value of me2 against an
	* analytic function.
	* @param me2 The value of the \f$ \|\bar{\mathcal{M}}\|^2 \f$
	*/
	virtual void debug(double me2) const;

	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 Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();
	//@}


	protected:

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

	/** 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 {return new_ptr(*this);}
	//@}

	private:

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

	private:

	/** @name Dynamically casted vertices. */
	//@{
	/**
	* Intermediate scalar
	*/
	vector<pair<AbstractVVSVertexPtr, AbstractFFSVertexPtr > > scalar_;
	/**
	* Intermediate fermion
	*/
	vector<pair<AbstractFFVVertexPtr, AbstractFFVVertexPtr> > fermion_;

	/**
	* Intermediate vector
	*/
	vector<pair<AbstractVVVVertexPtr, AbstractFFVVertexPtr> > vector_;

	/**
	* Intermediate tensor
	*/
	vector<pair<AbstractVVTVertexPtr, AbstractFFTVertexPtr> > tensor_;

	/**
	* Four point vertices
	*/
	vector<AbstractFFVVVertexPtr> four_;
	//@}
	};

	}

	#endif /* HERWIG_MEvv2ff_H */
	diff --git a/MatrixElement/General/MEvv2rf.cc b/MatrixElement/General/MEvv2rf.cc
	new file mode 100644
	--- /dev/null
	+++ b/MatrixElement/General/MEvv2rf.cc
	@@ -0,0 +1,377 @@
	+// -- C++ --
	+//
	+// This is the implementation of the non-inlined, non-templated member
	+// functions of the MEvv2rf class.
	+//
	+
	+#include "MEvv2rf.h"
	+#include "ThePEG/Interface/ClassDocumentation.h"
	+#include "ThePEG/Utilities/DescribeClass.h"
	+#include "ThePEG/Persistency/PersistentOStream.h"
	+#include "ThePEG/Persistency/PersistentIStream.h"
	+
	+using namespace Herwig;
	+using ThePEG::Helicity::incoming;
	+using ThePEG::Helicity::outgoing;
	+
	+IBPtr MEvv2rf::clone() const {
	+ return new_ptr(*this);
	+}
	+
	+IBPtr MEvv2rf::fullclone() const {
	+ return new_ptr(*this);
	+}
	+
	+void MEvv2rf::doinit() {
	+ GeneralHardME::doinit();
	+ scalar_ .resize(numberOfDiags());
	+ fermion_.resize(numberOfDiags());
	+ vector_ .resize(numberOfDiags());
	+ four_ .resize(numberOfDiags());
	+ initializeMatrixElements(PDT::Spin1 , PDT::Spin1,
	+ PDT::Spin3Half, PDT::Spin1Half);
	+
	+ for( size_t i = 0; i < numberOfDiags(); ++i ) {
	+ HPDiagram dg = getProcessInfo()[i];
	+ if( dg.channelType == HPDiagram::tChannel ) {
	+ AbstractRFVVertexPtr rfv;
	+ AbstractFFVVertexPtr ffv;
	+ if(dg.ordered.second) {
	+ rfv = dynamic_ptr_cast<AbstractRFVVertexPtr>(dg.vertices.first);
	+ ffv = dynamic_ptr_cast<AbstractFFVVertexPtr>(dg.vertices.second);
	+ }
	+ else {
	+ rfv = dynamic_ptr_cast<AbstractRFVVertexPtr>(dg.vertices.second);
	+ ffv = dynamic_ptr_cast<AbstractFFVVertexPtr>(dg.vertices.first );
	+ }
	+ fermion_[i] = make_pair(rfv, ffv);
	+ }
	+ else if( dg.channelType == HPDiagram::sChannel ) {
	+ if( dg.intermediate->iSpin() == PDT::Spin0 ) {
	+ AbstractVVSVertexPtr vvs =
	+ dynamic_ptr_cast<AbstractVVSVertexPtr>(dg.vertices.first );
	+ AbstractRFSVertexPtr rfs =
	+ dynamic_ptr_cast<AbstractRFSVertexPtr>(dg.vertices.second);
	+ scalar_[i] = make_pair(vvs,rfs);
	+ }
	+ else if( dg.intermediate->iSpin() == PDT::Spin1) {
	+ AbstractVVVVertexPtr vvv =
	+ dynamic_ptr_cast<AbstractVVVVertexPtr>(dg.vertices.first);
	+ AbstractRFVVertexPtr rfv =
	+ dynamic_ptr_cast<AbstractRFVVertexPtr>(dg.vertices.second);
	+ vector_[i] = make_pair(vvv,rfv);
	+ }
	+ }
	+ else if ( dg.channelType == HPDiagram::fourPoint) {
	+ four_[i] = dynamic_ptr_cast<AbstractRFVVVertexPtr>(dg.vertices.first);
	+ }
	+ }
	+}
	+
	+void MEvv2rf::persistentOutput(PersistentOStream & os) const {
	+ os << scalar_ << fermion_ << vector_ << four_;
	+}
	+
	+void MEvv2rf::persistentInput(PersistentIStream & is, int) {
	+ is >> scalar_ >> fermion_ >> vector_ >> four_;
	+ initializeMatrixElements(PDT::Spin1 , PDT::Spin1,
	+ PDT::Spin3Half, PDT::Spin1Half);
	+}
	+
	+//The following static variable is needed for the type
	+// description system in ThePEG.
	+DescribeClass<MEvv2rf,GeneralHardME>
	+describeHerwigMEvv2rf("Herwig::MEvv2rf", "Herwig.so");
	+
	+void MEvv2rf::Init() {
	+
	+ static ClassDocumentation<MEvv2rf> documentation
	+ ("The MEvv2rf class handes the ME calculation for vv -> rf");
	+
	+}
	+
	+double MEvv2rf::me2() const {
	+ // set up the vector wavefunctions
	+ VBVector v1(2), v2(2);
	+ for( size_t i = 0; i < 2; ++i ) {
	+ v1[i] = VectorWaveFunction(rescaledMomenta()[0],mePartonData()[0], 2*i,
	+ incoming);
	+ v2[i] = VectorWaveFunction(rescaledMomenta()[1],mePartonData()[1], 2*i,
	+ incoming);
	+ }
	+ // setup spinor wavefunctions and decide which case to use
	+ bool massless = mePartonData()[2]->mass()==ZERO;
	+ double full_me(0.);
	+ if(mePartonData()[2]->id()<0) {
	+ RSSpinorVector sp(4); SpinorBarVector sbar(2);
	+ for( size_t i = 0; i < 4; ++i ) {
	+ if(massless && (i==2\|\|i==3)) continue;
	+ sp[i] = RSSpinorWaveFunction(rescaledMomenta()[2], mePartonData()[2], i,
	+ outgoing);
	+ }
	+ for( size_t i = 0; i < 2; ++i ) {
	+ sbar[i] = SpinorBarWaveFunction(rescaledMomenta()[3], mePartonData()[3], i,
	+ outgoing);
	+ }
	+ vv2frME(v1, v2, sbar, sp, full_me,true);
	+ }
	+ else {
	+ SpinorVector sp(2); RSSpinorBarVector sbar(4);
	+ for( size_t i = 0; i < 4; ++i ) {
	+ if(massless && (i==2\|\|i==3)) continue;
	+ sbar[i] = RSSpinorBarWaveFunction(rescaledMomenta()[2], mePartonData()[2], i,
	+ outgoing);
	+ }
	+ for( size_t i = 0; i < 2; ++i ) {
	+ sp[i] = SpinorWaveFunction(rescaledMomenta()[3], mePartonData()[3], i,
	+ outgoing);
	+ }
	+ vv2rfME(v1, v2, sbar, sp, full_me,true);
	+ }
	+ return full_me;
	+}
	+
	+ProductionMatrixElement
	+MEvv2rf::vv2rfME(const VBVector & v1, const VBVector & v2,
	+ const RSSpinorBarVector & sbar,const SpinorVector & sp,
	+ double & me2, bool first) const {
	+ // scale
	+ const Energy2 q2 = scale();
	+ // whether or not rs fermion is massless
	+ bool massless = mePartonData()[2]->mass()==ZERO;
	+ // weights for the selection of the diagram
	+ vector<double> me(numberOfDiags(), 0.);
	+ // weights for the selection of the colour flow
	+ vector<double> flow(numberOfFlows(),0.);
	+ //sum over vector helicities
	+ for(unsigned int iv1 = 0; iv1 < 2; ++iv1) {
	+ for(unsigned int iv2 = 0; iv2 < 2; ++iv2) {
	+ //sum over fermion helicities
	+ for(unsigned int of1 = 0; of1 < 4; ++of1) {
	+ if(massless && (of1==1 \|\| of1==2) ) continue;
	+ for(unsigned int of2 = 0; of2 < 2; ++of2) {
	+ vector<Complex> flows(numberOfFlows(),0.);
	+ for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
	+ Complex diag(0.);
	+ const HPDiagram & current = getProcessInfo()[ix];
	+ tPDPtr offshell = current.intermediate;
	+ if(current.channelType == HPDiagram::tChannel &&
	+ offshell->iSpin() == PDT::Spin1Half) {
	+ if(current.ordered.second) {
	+ SpinorBarWaveFunction interF = fermion_[ix].first->
	+ evaluate(q2, 3, offshell, sbar[of1], v1[iv1]);
	+ diag = fermion_[ix].second->
	+ evaluate(q2, sp[of2], interF, v2[iv2]);
	+ }
	+ else {
	+ SpinorWaveFunction interF = fermion_[ix].second->
	+ evaluate(q2, 3, offshell, sp[of2], v1[iv1]);
	+ diag = fermion_[ix].first->
	+ evaluate(q2, interF, sbar[of1], v2[iv2]);
	+ }
	+ }
	+ else if(current.channelType == HPDiagram::sChannel) {
	+ if(offshell->iSpin() == PDT::Spin0) {
	+ ScalarWaveFunction interS = scalar_[ix].first->
	+ evaluate(q2, 1, offshell, v1[iv1], v2[iv2]);
	+ diag = scalar_[ix].second->
	+ evaluate(q2, sp[of2], sbar[of1], interS);
	+ }
	+ else if(offshell->iSpin() == PDT::Spin1) {
	+ VectorWaveFunction interV = vector_[ix].first->
	+ evaluate(q2, 1, offshell, v1[iv1], v2[iv2]);
	+ diag = vector_[ix].second->
	+ evaluate(q2, sp[of2], sbar[of1], interV);
	+ }
	+ }
	+ else if(current.channelType == HPDiagram::fourPoint) {
	+ diag = four_[ix]->evaluate(q2, sp[of2], sbar[of1], v1[iv1], v2[iv2]);
	+ }
	+ else
	+ assert(false);
	+ me[ix] += norm(diag);
	+ diagramME()[ix](2iv1, 2iv2, of1, of2) = diag;
	+ //Compute flows
	+ for(size_t iy = 0; iy < current.colourFlow.size(); ++iy) {
	+ assert(current.colourFlow[iy].first<flows.size());
	+ flows[current.colourFlow[iy].first] +=
	+ current.colourFlow[iy].second * diag;
	+ }
	+ }
	+ // MEs for the different colour flows
	+ for(unsigned int iy = 0; iy < numberOfFlows(); ++iy)
	+ flowME()[iy](2iv1, 2iv2, of1, of2) = flows[iy];
	+ //Now add flows to me2 with appropriate colour factors
	+ for(size_t ii = 0; ii < numberOfFlows(); ++ii)
	+ for(size_t ij = 0; ij < numberOfFlows(); ++ij)
	+ me2 += getColourFactors()[ii][ij](flows[ii]conj(flows[ij])).real();
	+ // contribution to the colour flow
	+ for(unsigned int ii = 0; ii < numberOfFlows(); ++ii) {
	+ flow[ii] += getColourFactors()[ii][ii]*norm(flows[ii]);
	+ }
	+ }
	+ }
	+ }
	+ }
	+ // if not computing the cross section return the selected colour flow
	+ if(!first) return flowME()[colourFlow()];
	+ me2 = selectColourFlow(flow,me,me2);
	+ return flowME()[colourFlow()];
	+}
	+
	+ProductionMatrixElement
	+MEvv2rf::vv2frME(const VBVector & v1, const VBVector & v2,
	+ const SpinorBarVector & sbar,const RSSpinorVector & sp,
	+ double & me2, bool first) const {
	+ // scale
	+ const Energy2 q2 = scale();
	+ // whether or not rs fermion is massless
	+ bool massless = mePartonData()[2]->mass()==ZERO;
	+ // weights for the selection of the diagram
	+ vector<double> me(numberOfDiags(), 0.);
	+ // weights for the selection of the colour flow
	+ vector<double> flow(numberOfFlows(),0.);
	+ //sum over vector helicities
	+ for(unsigned int iv1 = 0; iv1 < 2; ++iv1) {
	+ for(unsigned int iv2 = 0; iv2 < 2; ++iv2) {
	+ //sum over fermion helicities
	+ for(unsigned int of1 = 0; of1 < 4; ++of1) {
	+ for(unsigned int of2 = 0; of2 < 2; ++of2) {
	+ if(massless && (of2==1 \|\| of2==2) ) continue;
	+ vector<Complex> flows(numberOfFlows(),0.);
	+ for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
	+ Complex diag(0.);
	+ const HPDiagram & current = getProcessInfo()[ix];
	+ tPDPtr offshell = current.intermediate;
	+ if(current.channelType == HPDiagram::tChannel &&
	+ offshell->iSpin() == PDT::Spin1Half) {
	+ if(current.ordered.second) {
	+ SpinorWaveFunction interF = fermion_[ix].first->
	+ evaluate(q2, 3, offshell, sp[of2], v1[iv1]);
	+ diag = fermion_[ix].second->
	+ evaluate(q2, interF, sbar[of1],v2[iv2]);
	+ }
	+ else {
	+ SpinorBarWaveFunction interF = fermion_[ix].second->
	+ evaluate(q2, 3, offshell, sbar[of1], v1[iv1]);
	+ diag = fermion_[ix].first->
	+ evaluate(q2, sp[of2], interF, v2[iv2]);
	+ }
	+ }
	+ else if(current.channelType == HPDiagram::sChannel) {
	+ if(offshell->iSpin() == PDT::Spin0) {
	+ ScalarWaveFunction interS = scalar_[ix].first->
	+ evaluate(q2, 1, offshell, v1[iv1], v2[iv2]);
	+ diag = scalar_[ix].second->
	+ evaluate(q2, sp[of2], sbar[of1], interS);
	+ }
	+ else if(offshell->iSpin() == PDT::Spin1) {
	+ VectorWaveFunction interV = vector_[ix].first->
	+ evaluate(q2, 1, offshell, v1[iv1], v2[iv2]);
	+ diag = vector_[ix].second->
	+ evaluate(q2, sp[of2], sbar[of1], interV);
	+ }
	+ }
	+ else if(current.channelType == HPDiagram::fourPoint) {
	+ diag = four_[ix]->evaluate(q2, sp[of2], sbar[of1], v1[iv1], v2[iv2]);
	+ }
	+ else
	+ assert(false);
	+ me[ix] += norm(diag);
	+ diagramME()[ix](2iv1, 2iv2, of2, of1) = diag;
	+ //Compute flows
	+ for(size_t iy = 0; iy < current.colourFlow.size(); ++iy) {
	+ assert(current.colourFlow[iy].first<flows.size());
	+ flows[current.colourFlow[iy].first] +=
	+ current.colourFlow[iy].second * diag;
	+ }
	+ }
	+ // MEs for the different colour flows
	+ for(unsigned int iy = 0; iy < numberOfFlows(); ++iy)
	+ flowME()[iy](2iv1, 2iv2, of1, of2) = flows[iy];
	+ //Now add flows to me2 with appropriate colour factors
	+ for(size_t ii = 0; ii < numberOfFlows(); ++ii)
	+ for(size_t ij = 0; ij < numberOfFlows(); ++ij)
	+ me2 += getColourFactors()[ii][ij](flows[ii]conj(flows[ij])).real();
	+ // contribution to the colour flow
	+ for(unsigned int ii = 0; ii < numberOfFlows(); ++ii) {
	+ flow[ii] += getColourFactors()[ii][ii]*norm(flows[ii]);
	+ }
	+ }
	+ }
	+ }
	+ }
	+ // if not computing the cross section return the selected colour flow
	+ if(!first) return flowME()[colourFlow()];
	+ me2 = selectColourFlow(flow,me,me2);
	+ return flowME()[colourFlow()];
	+}
	+
	+void MEvv2rf::constructVertex(tSubProPtr sub) {
	+ ParticleVector ext = hardParticles(sub);
	+ // vector wavefunctions are common do them first
	+ // wavefunction with real momenta
	+ VBVector v1, v2;
	+ VectorWaveFunction(v1, ext[0], incoming, false, true);
	+ VectorWaveFunction(v2, ext[1], incoming, false, true);
	+ // rescale momenta
	+ setRescaledMomenta(ext);
	+ // wavefunctions with rescaled momenta
	+ VectorWaveFunction v1r(rescaledMomenta()[0],
	+ ext[0]->dataPtr(), incoming);
	+ VectorWaveFunction v2r(rescaledMomenta()[1],
	+ ext[1]->dataPtr(), incoming);
	+ ProductionMatrixElement pme;
	+ for( unsigned int ihel = 0; ihel < 2; ++ihel ) {
	+ v1r.reset(2*ihel);
	+ v1[ihel] = v1r;
	+ v2r.reset(2*ihel);
	+ v2[ihel] = v2r;
	+ }
	+ // setup spinor wavefunctions and decide which case to use
	+ bool massless = mePartonData()[2]->mass()==ZERO;
	+ double dummy(0.);
	+ if(ext[2]->id()<0) {
	+ RSSpinorVector sp;
	+ RSSpinorWaveFunction(sp, ext[2], outgoing, true);
	+ SpinorBarVector sbar;
	+ SpinorBarWaveFunction(sbar, ext[3], outgoing, true);
	+ // wavefuncions with rescaled momenta
	+ SpinorBarWaveFunction sbr(rescaledMomenta()[3],
	+ ext[3]->dataPtr(), outgoing);
	+ RSSpinorWaveFunction spr(rescaledMomenta()[2],
	+ ext[2]->dataPtr(), outgoing);
	+ for( unsigned int ihel = 0; ihel < 2; ++ihel ) {
	+ sbr.reset(ihel);
	+ sbar[ihel] = sbr;
	+ }
	+ for( unsigned int ihel = 0; ihel < 4; ++ihel ) {
	+ if(massless && (ihel==1 \|\| ihel==2)) continue;
	+ spr.reset(ihel);
	+ sp[ihel] = spr;
	+ }
	+ pme = vv2frME(v1, v2, sbar, sp, dummy,false);
	+ }
	+ else {
	+ SpinorVector sp;
	+ SpinorWaveFunction(sp, ext[3], outgoing, true);
	+ RSSpinorBarVector sbar;
	+ RSSpinorBarWaveFunction(sbar, ext[2], outgoing, true);
	+ // wavefuncions with rescaled momenta
	+ RSSpinorBarWaveFunction sbr(rescaledMomenta()[2],
	+ ext[2]->dataPtr(), outgoing);
	+ SpinorWaveFunction spr(rescaledMomenta()[3],
	+ ext[3]->dataPtr(), outgoing);
	+ for( unsigned int ihel = 0; ihel < 4; ++ihel ) {
	+ if(massless && (ihel==1 \|\| ihel==2)) continue;
	+ sbr.reset(ihel);
	+ sbar[ihel] = sbr;
	+ }
	+ for( unsigned int ihel = 0; ihel < 2; ++ihel ) {
	+ spr.reset(ihel);
	+ sp[ihel] = spr;
	+ }
	+ pme = vv2rfME(v1, v2, sbar, sp, dummy,false);
	+ }
	+ createVertex(pme,ext);
	+}
	diff --git a/MatrixElement/General/MEvv2rf.h b/MatrixElement/General/MEvv2rf.h
	new file mode 100644
	--- /dev/null
	+++ b/MatrixElement/General/MEvv2rf.h
	@@ -0,0 +1,182 @@
	+// -- C++ --
	+#ifndef Herwig_MEvv2rf_H
	+#define Herwig_MEvv2rf_H
	+//
	+// This is the declaration of the MEvv2rf class.
	+//
	+
	+#include "GeneralHardME.h"
	+#include "ThePEG/Helicity/Vertex/AbstractRFSVertex.h"
	+#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
	+#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
	+#include "ThePEG/Helicity/Vertex/AbstractRFVVertex.h"
	+#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
	+#include "ThePEG/Helicity/Vertex/AbstractRFVVVertex.h"
	+#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	+#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	+#include "ThePEG/Helicity/WaveFunction/RSSpinorWaveFunction.h"
	+#include "ThePEG/Helicity/WaveFunction/RSSpinorBarWaveFunction.h"
	+#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	+#include "Herwig/MatrixElement/ProductionMatrixElement.h"
	+
	+namespace Herwig {
	+
	+using namespace ThePEG;
	+
	+/**
	+ * This class is designed to implement the matrix element for the
	+ * \f$2 \rightarrow 2\f$ process vector-vector to fermion- RS fermion. It
	+ * inherits from GeneralHardME and implements the me2() virtual function.
	+ *
	+ * @see GeneralHardME
	+ *
	+ */
	+class MEvv2rf: public GeneralHardME {
	+
	+public:
	+
	+ /** A Vector of VectorWaveFunction objects. */
	+ typedef vector<VectorWaveFunction> VBVector;
	+
	+ /** A vector of SpinorBarWaveFunction objects. */
	+ typedef vector<SpinorWaveFunction> SpinorVector;
	+
	+ /** A vector of SpinorBarWaveFunction objects. */
	+ typedef vector<SpinorBarWaveFunction> SpinorBarVector;
	+
	+ /** A vector of SpinorBarWaveFunction objects. */
	+ typedef vector<RSSpinorWaveFunction> RSSpinorVector;
	+
	+ /** A vector of SpinorBarWaveFunction objects. */
	+ typedef vector<RSSpinorBarWaveFunction> RSSpinorBarVector;
	+
	+public:
	+
	+ /** @name Virtual functions required by the MEBase class. */
	+ //@{
	+ /**
	+ * The matrix element for the kinematical configuration
	+ * previously provided by the last call to setKinematics(), suitably
	+ * scaled by sHat() to give a dimension-less number.
	+ * @return the matrix element scaled with sHat() to give a
	+ * dimensionless number.
	+ */
	+ virtual double me2() const;
	+ //@}
	+
	+ /**
	+ * Construct the vertex information for the spin correlations
	+ * @param sub Pointer to the relevent SubProcess
	+ */
	+ virtual void constructVertex(tSubProPtr sub);
	+
	+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:
	+
	+ /** @name Standard Interfaced functions. */
	+ //@{
	+ /**
	+ * Initialize this object after the setup phase before saving an
	+ * EventGenerator to disk.
	+ * @throws InitException if object could not be initialized properly.
	+ */
	+ virtual void doinit();
	+ //@}
	+
	+private:
	+
	+ /**
	+ * Calculate the value of the matrix element
	+ */
	+ ProductionMatrixElement vv2rfME(const VBVector & v1, const VBVector & v2,
	+ const RSSpinorBarVector & sbar,
	+ const SpinorVector & sp,
	+ double & me2, bool first) const;
	+
	+ /**
	+ * Calculate the value of the matrix element
	+ */
	+ ProductionMatrixElement vv2frME(const VBVector & v1, const VBVector & v2,
	+ const SpinorBarVector & sbar,
	+ const RSSpinorVector & sp,
	+ double & me2, bool first) const;
	+
	+private:
	+
	+ /**
	+ * The assignment operator is private and must never be called.
	+ * In fact, it should not even be implemented.
	+ */
	+ MEvv2rf & operator=(const MEvv2rf &);
	+private:
	+
	+ /** @name Dynamically casted vertices. */
	+ //@{
	+ /**
	+ * Intermediate scalar
	+ */
	+ vector<pair<AbstractVVSVertexPtr, AbstractRFSVertexPtr > > scalar_;
	+
	+ /**
	+ * Intermediate fermion
	+ */
	+ vector<pair<AbstractRFVVertexPtr, AbstractFFVVertexPtr> > fermion_;
	+
	+ /**
	+ * Intermediate vector
	+ */
	+ vector<pair<AbstractVVVVertexPtr, AbstractRFVVertexPtr> > vector_;
	+
	+ /**
	+ * Four point vertices
	+ */
	+ vector<AbstractRFVVVertexPtr> four_;
	+ //@}
	+
	+};
	+
	+}
	+
	+#endif /* Herwig_MEvv2rf_H */
	diff --git a/MatrixElement/General/Makefile.am b/MatrixElement/General/Makefile.am
	--- a/MatrixElement/General/Makefile.am
	+++ b/MatrixElement/General/Makefile.am
	@@ -1,20 +1,21 @@
	noinst_LTLIBRARIES = libHwGeneralME.la
	libHwGeneralME_la_SOURCES = \
	GeneralHardME.cc GeneralHardME.h GeneralHardME.fh \
	MEvv2ff.cc MEvv2ff.h \
	+MEvv2rf.cc MEvv2rf.h \
	MEvv2ss.cc MEvv2ss.h \
	MEfv2fs.cc MEfv2fs.h \
	MEff2ss.cc MEff2ss.h \
	MEff2ff.cc MEff2ff.h \
	MEff2vv.cc MEff2vv.h \
	MEfv2vf.cc MEfv2vf.h \
	MEff2vs.cc MEff2vs.h \
	MEff2sv.cc MEff2sv.h \
	MEvv2vv.cc MEvv2vv.h \
	MEvv2vs.cc MEvv2vs.h \
	MEff2tv.cc MEff2tv.h \
	MEvv2tv.cc MEvv2tv.h \
	MEfv2tf.cc MEfv2tf.h \
	GeneralfftoVH.cc GeneralfftoVH.h GeneralfftoVH.fh\
	GeneralfftoffH.cc GeneralfftoffH.h GeneralfftoffH.fh\
	GeneralQQHiggs.cc GeneralQQHiggs.h GeneralQQHiggs.fh
	diff --git a/Models/General/TwoToTwoProcessConstructor.cc b/Models/General/TwoToTwoProcessConstructor.cc
	--- a/Models/General/TwoToTwoProcessConstructor.cc
	+++ b/Models/General/TwoToTwoProcessConstructor.cc
	@@ -1,643 +1,644 @@
	// -- C++ --
	//
	// TwoToTwoProcessConstructor.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 TwoToTwoProcessConstructor class.
	//

	#include "TwoToTwoProcessConstructor.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include <sstream>

	using std::stringstream;

	using namespace Herwig;

	TwoToTwoProcessConstructor::TwoToTwoProcessConstructor() :
	Nout_(0), nv_(0), allDiagrams_(true),
	processOption_(0), scaleChoice_(0), scaleFactor_(1.)
	{}

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

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

	void TwoToTwoProcessConstructor::doinit() {
	HardProcessConstructor::doinit();
	if((processOption_==2 \|\| processOption_==4) &&
	outgoing_.size()!=2)
	throw InitException()
	<< "Exclusive processes require exactly"
	<< " two outgoing particles but " << outgoing_.size()
	<< "have been inserted in TwoToTwoProcessConstructor::doinit()."
	<< Exception::runerror;
	if(processOption_==4 && incoming_.size()!=2)
	throw InitException()
	<< "Exclusive processes require exactly"
	<< " two incoming particles but " << incoming_.size()
	<< "have been inserted in TwoToTwoProcessConstructor::doinit()."
	<< Exception::runerror;
	Nout_ = outgoing_.size();
	PDVector::size_type ninc = incoming_.size();
	// exit if nothing to do
	if(Nout_==0\|\|ninc==0) return;
	//create vector of initial-state pairs
	for(PDVector::size_type i = 0; i < ninc; ++i) {
	for(PDVector::size_type j = 0; j < ninc; ++j) {
	tPDPair inc = make_pair(incoming_[i], incoming_[j]);

	if( (inc.first->iSpin() > inc.second->iSpin()) \|\|
	(inc.first->iSpin() == inc.second->iSpin() &&
	inc.first->id() < inc.second->id()) )
	swap(inc.first, inc.second);

	if( !HPC_helper::duplicateIncoming(inc,incPairs_) ) {
	incPairs_.push_back(inc);
	}
	}
	}
	// excluded vertices
	excludedVertexSet_ =
	set<VertexBasePtr>(excludedVertexVector_.begin(),
	excludedVertexVector_.end());
	}


	void TwoToTwoProcessConstructor::persistentOutput(PersistentOStream & os) const {
	os << vertices_ << incoming_ << outgoing_
	<< allDiagrams_ << processOption_
	<< scaleChoice_ << scaleFactor_ << excluded_ << excludedExternal_
	<< excludedVertexVector_ << excludedVertexSet_;
	}

	void TwoToTwoProcessConstructor::persistentInput(PersistentIStream & is, int) {
	is >> vertices_ >> incoming_ >> outgoing_
	>> allDiagrams_ >> processOption_
	>> scaleChoice_ >> scaleFactor_ >> excluded_ >> excludedExternal_
	>> excludedVertexVector_ >> excludedVertexSet_;
	}

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

	void TwoToTwoProcessConstructor::Init() {

	static ClassDocumentation<TwoToTwoProcessConstructor> documentation
	("TwoToTwoProcessConstructor constructs the possible diagrams for "
	"a process given the external particles");

	static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceIn
	("Incoming",
	"Pointers to incoming particles",
	&TwoToTwoProcessConstructor::incoming_, -1, false, false, true, false);

	static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceOut
	("Outgoing",
	"Pointers to incoming particles",
	&TwoToTwoProcessConstructor::outgoing_, -1, false, false, true, false);

	static Switch<TwoToTwoProcessConstructor,bool> interfaceIncludeAllDiagrams
	("IncludeEW",
	"Switch to decide which diagrams to include in ME calc.",
	&TwoToTwoProcessConstructor::allDiagrams_, true, false, false);
	static SwitchOption interfaceIncludeAllDiagramsNo
	(interfaceIncludeAllDiagrams,
	"No",
	"Only include QCD diagrams",
	false);
	static SwitchOption interfaceIncludeAllDiagramsYes
	(interfaceIncludeAllDiagrams,
	"Yes",
	"Include EW+QCD.",
	true);

	static Switch<TwoToTwoProcessConstructor,unsigned int> interfaceProcesses
	("Processes",
	"Whether to generate inclusive or exclusive processes",
	&TwoToTwoProcessConstructor::processOption_, 0, false, false);
	static SwitchOption interfaceProcessesSingleParticleInclusive
	(interfaceProcesses,
	"SingleParticleInclusive",
	"Require at least one particle from the list of outgoing particles"
	" in the hard process",
	0);
	static SwitchOption interfaceProcessesTwoParticleInclusive
	(interfaceProcesses,
	"TwoParticleInclusive",
	"Require that both the particles in the hard processes are in the"
	" list of outgoing particles",
	1);
	static SwitchOption interfaceProcessesExclusive
	(interfaceProcesses,
	"Exclusive",
	"Require that both the particles in the hard processes are in the"
	" list of outgoing particles in every hard process",
	2);
	static SwitchOption interfaceProcessesVeryExclusive
	(interfaceProcesses,
	"VeryExclusive",
	"Require that both the incoming and outgoing particles in the hard processes are in the"
	" list of outgoing particles in every hard process",
	4);

	static Switch<TwoToTwoProcessConstructor,unsigned int> interfaceScaleChoice
	("ScaleChoice",
	"&TwoToTwoProcessConstructor::scaleChoice_",
	&TwoToTwoProcessConstructor::scaleChoice_, 0, false, false);
	static SwitchOption interfaceScaleChoiceDefault
	(interfaceScaleChoice,
	"Default",
	"Use if sHat if intermediates all colour neutral, otherwise the transverse mass",
	0);
	static SwitchOption interfaceScaleChoicesHat
	(interfaceScaleChoice,
	"sHat",
	"Always use sHat",
	1);
	static SwitchOption interfaceScaleChoiceTransverseMass
	(interfaceScaleChoice,
	"TransverseMass",
	"Always use the transverse mass",
	2);
	static SwitchOption interfaceScaleChoiceGeometicMean
	(interfaceScaleChoice,
	"MaxMT",
	"Use the maximum of m^2+p_T^2 for the two particles",
	3);

	static Parameter<TwoToTwoProcessConstructor,double> interfaceScaleFactor
	("ScaleFactor",
	"The prefactor used in the scale calculation. The scale used is"
	" that defined by scaleChoice multiplied by this prefactor",
	&TwoToTwoProcessConstructor::scaleFactor_, 1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceExcluded
	("Excluded",
	"Particles which are not allowed as intermediates",
	&TwoToTwoProcessConstructor::excluded_, -1, false, false, true, false, false);

	static RefVector<TwoToTwoProcessConstructor,ParticleData> interfaceExcludedExternal
	("ExcludedExternal",
	"Particles which are not allowed as outgoing particles",
	&TwoToTwoProcessConstructor::excludedExternal_, -1,
	false, false, true, false, false);

	static RefVector<TwoToTwoProcessConstructor,VertexBase> interfaceExcludedVertices
	("ExcludedVertices",
	"Vertices which are not included in the 2 -> 2 scatterings",
	&TwoToTwoProcessConstructor::excludedVertexVector_, -1, false, false, true, true, false);

	}

	void TwoToTwoProcessConstructor::constructDiagrams() {
	if(incPairs_.empty() \|\| outgoing_.empty() \|\| !subProcess() ) return;
	nv_ = model()->numberOfVertices();
	//make sure vertices are initialised
	for(unsigned int ix = 0; ix < nv_; ++ix ) {
	VertexBasePtr vertex = model()->vertex(ix);
	if(excludedVertexSet_.find(vertex) !=
	excludedVertexSet_.end()) continue;
	vertices_.push_back(vertex);
	}
	nv_ = vertices_.size();
	//Create necessary diagrams
	vector<tcPDPair>::size_type is;
	PDVector::size_type os;
	for(is = 0; is < incPairs_.size(); ++is) {
	tPDPair ppi = incPairs_[is];
	for(os = 0; os < Nout_; ++os) {
	long fs = outgoing_[os]->id();
	for(size_t iv = 0; iv < nv_; ++iv) {
	tVertexBasePtr vertexA = vertices_[iv];
	//This skips an effective vertex and the EW ones if
	// we only want the strong diagrams
	if( !allDiagrams_ && vertexA->orderInGs() == 0 )
	continue;

	if(vertexA->getNpoint() == 3) {
	//scattering diagrams
	createTChannels(ppi, fs, vertexA);

	//resonance diagrams
	if( vertexA->isIncoming(ppi.first) &&
	vertexA->isIncoming(ppi.second) )
	createSChannels(ppi, fs, vertexA);
	}
	else {
	makeFourPointDiagrams(ppi.first->id(), ppi.second->id(),
	fs, vertexA);
	}
	}
	}
	}
	// need to find all of the diagrams that relate to the same process
	// first insert them into a map which uses the '<' operator
	// to sort the diagrams
	multiset<HPDiagram> grouped;
	HPDVector::iterator dit = processes_.begin();
	HPDVector::iterator dend = processes_.end();
	for( ; dit != dend; ++dit) {
	grouped.insert(*dit);
	}
	assert( processes_.size() == grouped.size() );
	processes_.clear();
	typedef multiset<HPDiagram>::const_iterator set_iter;
	set_iter it = grouped.begin(), iend = grouped.end();
	while( it != iend ) {
	pair<set_iter,set_iter> range = grouped.equal_range(*it);
	set_iter itb = range.first;
	HPDVector process;
	for( ; itb != range.second; ++itb ) {
	process.push_back(*itb);
	}
	// if inclusive enforce the exclusivity
	if(processOption_==2 \|\| processOption_==4) {
	if(!((process[0].outgoing. first==outgoing_[0]->id()&&
	process[0].outgoing.second==outgoing_[1]->id())\|\|
	(process[0].outgoing. first==outgoing_[1]->id()&&
	process[0].outgoing.second==outgoing_[0]->id()))) {
	process.clear();
	it = range.second;
	continue;
	}
	if(processOption_==4) {
	if(!((process[0].incoming. first==incoming_[0]->id()&&
	process[0].incoming.second==incoming_[1]->id())\|\|
	(process[0].incoming. first==incoming_[1]->id()&&
	process[0].incoming.second==incoming_[0]->id()))) {
	process.clear();
	it = range.second;
	continue;
	}
	}
	}
	// check no zero width s-channel intermediates
	for( dit=process.begin(); dit != process.end(); ++dit) {
	tPDPtr out1 = getParticleData(dit->outgoing.first );
	tPDPtr out2 = getParticleData(dit->outgoing.second);
	if(dit->channelType == HPDiagram::sChannel &&
	dit->intermediate->width()==ZERO &&
	dit->intermediate->mass() > out1->mass()+ out2->mass()) {
	tPDPtr in1 = getParticleData(dit->incoming.first );
	tPDPtr in2 = getParticleData(dit->incoming.second);
	throw Exception() << "Process with zero width resonant intermediates\n"
	<< dit->intermediate->PDGName()
	<< " can be on-shell in the process "
	<< in1 ->PDGName() << " " << in2->PDGName() << " -> "
	<< out1->PDGName() << " " << out2->PDGName()
	<< " but has zero width.\nEither set the width, enable "
	<< "calculation of its decays, and hence the width,\n"
	<< "or disable it as a potential intermediate using\n"
	<< "insert " << fullName() << ":Excluded 0 "
	<< dit->intermediate->fullName() << "\n---\n"
	<< Exception::runerror;
	}
	}
	if(find(excludedExternal_.begin(),excludedExternal_.end(),
	getParticleData(process[0].outgoing. first))!=excludedExternal_.end()) {
	process.clear();
	it = range.second;
	continue;
	}
	if(find(excludedExternal_.begin(),excludedExternal_.end(),
	getParticleData(process[0].outgoing.second))!=excludedExternal_.end()) {
	process.clear();
	it = range.second;
	continue;
	}
	// finally if the process is allow assign the colour flows
	for(unsigned int ix=0;ix<process.size();++ix) assignToCF(process[ix]);
	// create the matrix element
	createMatrixElement(process);
	process.clear();
	it = range.second;
	}
	}

	void TwoToTwoProcessConstructor::
	createSChannels(tcPDPair inpp, long fs, tVertexBasePtr vertex) {
	//Have 2 incoming particle and a vertex, find the possible offshell
	//particles
	pair<long,long> inc = make_pair(inpp.first->id(), inpp.second->id());
	tPDSet offshells = search(vertex, inpp.first->id(), incoming,
	inpp.second->id(), incoming, outgoing);
	tPDSet::const_iterator it;
	for(it = offshells.begin(); it != offshells.end(); ++it) {
	if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
	for(size_t iv = 0; iv < nv_; ++iv) {
	tVertexBasePtr vertexB = vertices_[iv];
	if( vertexB->getNpoint() != 3) continue;
	if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;

	tPDSet final;
	if( vertexB->isOutgoing(getParticleData(fs)) &&
	vertexB->isIncoming(*it) )
	final = search(vertexB, (*it)->id(), incoming, fs,
	outgoing, outgoing);
	//Now make diagrams
	if(!final.empty())
	makeDiagrams(inc, fs, final, *it, HPDiagram::sChannel,
	make_pair(vertex, vertexB), make_pair(true,true));
	}
	}

	}

	void TwoToTwoProcessConstructor::
	createTChannels(tPDPair inpp, long fs, tVertexBasePtr vertex) {
	pair<long,long> inc = make_pair(inpp.first->id(), inpp.second->id());
	//first try a with c
	tPDSet offshells = search(vertex, inpp.first->id(), incoming, fs,
	outgoing, outgoing);
	tPDSet::const_iterator it;
	for(it = offshells.begin(); it != offshells.end(); ++it) {
	if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
	for(size_t iv = 0; iv < nv_; ++iv) {
	tVertexBasePtr vertexB = vertices_[iv];
	if( vertexB->getNpoint() != 3 ) continue;
	if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;
	tPDSet final;
	if( vertexB->isIncoming(inpp.second) )
	final = search(vertexB, inc.second, incoming, (*it)->id(),
	incoming, outgoing);
	if( !final.empty() )
	makeDiagrams(inc, fs, final, *it, HPDiagram::tChannel,
	make_pair(vertex, vertexB), make_pair(true,true));
	}
	}
	//now try b with c
	offshells = search(vertex, inpp.second->id(), incoming, fs,
	outgoing, incoming);
	for(it = offshells.begin(); it != offshells.end(); ++it) {
	if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
	for(size_t iv = 0; iv < nv_; ++iv) {
	tVertexBasePtr vertexB = vertices_[iv];
	if( vertexB->getNpoint() != 3 ) continue;
	if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;

	tPDSet final;
	if( vertexB->isIncoming(inpp.first) )
	final = search(vertexB, inc.first, incoming, (*it)->id(),
	outgoing, outgoing);
	if( !final.empty() )
	makeDiagrams(inc, fs, final, *it, HPDiagram::tChannel,
	make_pair(vertexB, vertex), make_pair(true, false));
	}
	}

	}

	void TwoToTwoProcessConstructor::makeFourPointDiagrams(long parta, long partb,
	long partc, VBPtr vert) {
	if(processOption_>=1) {
	PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
	getParticleData(partc));
	if(loc==outgoing_.end()) return;
	}
	tPDSet ext = search(vert, parta, incoming, partb,incoming, partc, outgoing);
	if( ext.empty() ) return;
	IDPair in(parta, partb);
	for(tPDSet::const_iterator iter=ext.begin(); iter!=ext.end();
	++iter) {
	if(processOption_>=1) {
	PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
	*iter);
	if(loc==outgoing_.end()) continue;
	}
	HPDiagram nhp(in,make_pair(partc, (*iter)->id()));
	nhp.vertices = make_pair(vert, vert);
	nhp.channelType = HPDiagram::fourPoint;
	fixFSOrder(nhp);
	if(!checkOrder(nhp)) continue;
	if( !duplicate(nhp, processes_) ) processes_.push_back(nhp);
	}
	}

	void
	TwoToTwoProcessConstructor::makeDiagrams(IDPair in, long out1, const tPDSet & out2,
	PDPtr inter, HPDiagram::Channel chan,
	VBPair vertexpair, BPair cross) {
	if(processOption_>=1) {
	PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
	getParticleData(out1));
	if(loc==outgoing_.end()) return;
	}
	for(tPDSet::const_iterator it = out2.begin(); it != out2.end(); ++it) {
	if(processOption_>=1) {
	PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
	*it);
	if(loc==outgoing_.end()) continue;
	}
	HPDiagram nhp( in,make_pair(out1, (*it)->id()) );
	nhp.intermediate = inter;
	nhp.vertices = vertexpair;
	nhp.channelType = chan;
	nhp.ordered = cross;
	fixFSOrder(nhp);
	if(!checkOrder(nhp)) continue;
	if( !duplicate(nhp, processes_) ) processes_.push_back(nhp);
	}
	}

	set<tPDPtr>
	TwoToTwoProcessConstructor::search(VBPtr vertex, long part1, direction d1,
	long part2, direction d2, direction d3) {
	if(vertex->getNpoint() != 3) return tPDSet();
	if(d1 == incoming && getParticleData(part1)->CC()) part1 = -part1;
	if(d2 == incoming && getParticleData(part2)->CC()) part2 = -part2;
	vector<long> ext;
	tPDSet third;
	for(unsigned int ix = 0;ix < 3; ++ix) {
	vector<long> pdlist = vertex->search(ix, part1);
	ext.insert(ext.end(), pdlist.begin(), pdlist.end());
	}
	for(unsigned int ix = 0; ix < ext.size(); ix += 3) {
	long id0 = ext.at(ix);
	long id1 = ext.at(ix+1);
	long id2 = ext.at(ix+2);
	int pos;
	if((id0 == part1 && id1 == part2) \|\|
	(id0 == part2 && id1 == part1))
	pos = ix + 2;
	else if((id0 == part1 && id2 == part2) \|\|
	(id0 == part2 && id2 == part1))
	pos = ix + 1;
	else if((id1 == part1 && id2 == part2) \|\|
	(id1 == part2 && id2 == part1))
	pos = ix;
	else
	pos = -1;
	if(pos >= 0) {
	tPDPtr p = getParticleData(ext[pos]);
	if(d3 == incoming && p->CC()) p = p->CC();
	third.insert(p);
	}
	}

	return third;
	}

	set<tPDPtr>
	TwoToTwoProcessConstructor::search(VBPtr vertex,
	long part1, direction d1,
	long part2, direction d2,
	long part3, direction d3,
	direction d4) {
	if(vertex->getNpoint() != 4) return tPDSet();
	if(d1 == incoming && getParticleData(part1)->CC()) part1 = -part1;
	if(d2 == incoming && getParticleData(part2)->CC()) part2 = -part2;
	if(d3 == incoming && getParticleData(part3)->CC()) part3 = -part3;
	vector<long> ext;
	tPDSet fourth;
	for(unsigned int ix = 0;ix < 4; ++ix) {
	vector<long> pdlist = vertex->search(ix, part1);
	ext.insert(ext.end(), pdlist.begin(), pdlist.end());
	}
	for(unsigned int ix = 0;ix < ext.size(); ix += 4) {
	long id0 = ext.at(ix); long id1 = ext.at(ix + 1);
	long id2 = ext.at(ix + 2); long id3 = ext.at(ix + 3);
	int pos;
	if((id0 == part1 && id1 == part2 && id2 == part3) \|\|
	(id0 == part1 && id1 == part3 && id2 == part2) \|\|
	(id0 == part2 && id1 == part1 && id2 == part3) \|\|
	(id0 == part2 && id1 == part3 && id2 == part1) \|\|
	(id0 == part3 && id1 == part1 && id2 == part2) \|\|
	(id0 == part3 && id1 == part2 && id2 == part1))
	pos = ix + 3;
	else if((id0 == part1 && id1 == part2 && id3 == part3) \|\|
	(id0 == part1 && id1 == part3 && id3 == part2) \|\|
	(id0 == part2 && id1 == part1 && id3 == part3) \|\|
	(id0 == part2 && id1 == part3 && id3 == part1) \|\|
	(id0 == part3 && id1 == part1 && id3 == part2) \|\|
	(id0 == part3 && id1 == part2 && id3 == part1))
	pos = ix + 2;
	else if((id0 == part1 && id2 == part2 && id3 == part3) \|\|
	(id0 == part1 && id2 == part3 && id3 == part2) \|\|
	(id0 == part2 && id2 == part1 && id3 == part3) \|\|
	(id0 == part2 && id2 == part3 && id3 == part1) \|\|
	(id0 == part3 && id2 == part1 && id3 == part2) \|\|
	(id0 == part3 && id2 == part2 && id3 == part1))
	pos = ix + 1;
	else if((id1 == part1 && id2 == part2 && id3 == part3) \|\|
	(id1 == part1 && id2 == part3 && id3 == part2) \|\|
	(id1 == part2 && id2 == part1 && id3 == part3) \|\|
	(id1 == part2 && id2 == part3 && id3 == part1) \|\|
	(id1 == part3 && id2 == part1 && id3 == part2) \|\|
	(id1 == part3 && id2 == part2 && id3 == part1))
	pos = ix;
	else
	pos = -1;

	if(pos >= 0) {
	tPDPtr p = getParticleData(ext[pos]);
	if(d4 == incoming && p->CC())
	p = p->CC();
	fourth.insert(p);
	}
	}
	return fourth;
	}

	void
	TwoToTwoProcessConstructor::createMatrixElement(const HPDVector & process) const {
	if ( process.empty() ) return;
	// external particles
	tcPDVector extpart(4);
	extpart[0] = getParticleData(process[0].incoming.first);
	extpart[1] = getParticleData(process[0].incoming.second);
	extpart[2] = getParticleData(process[0].outgoing.first);
	extpart[3] = getParticleData(process[0].outgoing.second);
	// create the object
	string objectname ("/Herwig/MatrixElements/");
	string classname = MEClassname(extpart, objectname);
	GeneralHardMEPtr matrixElement = dynamic_ptr_cast<GeneralHardMEPtr>
	(generator()->preinitCreate(classname, objectname));
	if( !matrixElement ) {
	std::stringstream message;
	message << "TwoToTwoProcessConstructor::createMatrixElement "
	<< "- No matrix element object could be created for "
	<< "the process "
	<< extpart[0]->PDGName() << " " << extpart[0]->iSpin() << ","
	<< extpart[1]->PDGName() << " " << extpart[1]->iSpin() << "->"
	<< extpart[2]->PDGName() << " " << extpart[2]->iSpin() << ","
	<< extpart[3]->PDGName() << " " << extpart[3]->iSpin()
	<< ". Constructed class name: \"" << classname << "\"";
	generator()->logWarning(TwoToTwoProcessConstructorError(message.str(),Exception::warning));
	return;
	}
	// choice for the scale
	unsigned int scale;
	if(scaleChoice_==0) {
	// check coloured initial and final state
	bool inColour = ( extpart[0]->coloured() \|\|
	extpart[1]->coloured());
	bool outColour = ( extpart[2]->coloured() \|\|
	extpart[3]->coloured());
	if(inColour&&outColour) {
	bool coloured = false;
	for(unsigned int ix=0;ix<process.size();++ix) {
	if(process[ix].intermediate&&
	process[ix].intermediate->coloured()) {
	coloured = true;
	break;
	}
	}
	scale = coloured ? 1 : 0;
	}
	else {
	scale = 0;
	}
	}
	else {
	scale = scaleChoice_-1;
	}
	// set the information
	matrixElement->setProcessInfo(process, colourFlow(extpart),
	debug(), scale, scaleFactor_);
	// insert it
	generator()->preinitInterface(subProcess(), "MatrixElements",
	subProcess()->MEs().size(),
	"insert", matrixElement->fullName());
	}

	string TwoToTwoProcessConstructor::MEClassname(const vector<tcPDPtr> & extpart,
	string & objname) const {
	string classname("Herwig::ME");
	for(vector<tcPDPtr>::size_type ix = 0; ix < extpart.size(); ++ix) {
	if(ix == 2) classname += "2";
	if(extpart[ix]->iSpin() == PDT::Spin0) classname += "s";
	else if(extpart[ix]->iSpin() == PDT::Spin1) classname += "v";
	else if(extpart[ix]->iSpin() == PDT::Spin1Half) classname += "f";
	+ else if(extpart[ix]->iSpin() == PDT::Spin3Half) classname += "r";
	else if(extpart[ix]->iSpin() == PDT::Spin2) classname += "t";
	else {
	std::stringstream message;
	message << "MEClassname() : Encountered an unknown spin for "
	<< extpart[ix]->PDGName() << " while constructing MatrixElement "
	<< "classname " << extpart[ix]->iSpin();
	generator()->logWarning(TwoToTwoProcessConstructorError(message.str(),Exception::warning));
	}
	}
	objname += "ME" + extpart[0]->PDGName() + extpart[1]->PDGName() + "2"
	+ extpart[2]->PDGName() + extpart[3]->PDGName();
	return classname;
	}

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:36 PM (22 h, 32 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486658
Default Alt Text: (49 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions