No OneTemporary
Actions

Size

39 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/General/MEvv2ss.cc b/MatrixElement/General/MEvv2ss.cc
	--- a/MatrixElement/General/MEvv2ss.cc
	+++ b/MatrixElement/General/MEvv2ss.cc
	@@ -1,283 +1,309 @@
	// -- C++ --
	//
	// MEvv2ss.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig is licenced under version 2 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 MEvv2ss class.
	//

	#include "MEvv2ss.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace Herwig;
	using ThePEG::Helicity::TensorWaveFunction;
	using ThePEG::Helicity::incoming;
	using ThePEG::Helicity::outgoing;


	void MEvv2ss::doinit() {
	GeneralHardME::doinit();
	scalar1_.resize(numberOfDiags());
	scalar2_.resize(numberOfDiags());
	+ scalar3_.resize(numberOfDiags());
	vector_ .resize(numberOfDiags());
	tensor_ .resize(numberOfDiags());
	initializeMatrixElements(PDT::Spin1, PDT::Spin1,
	PDT::Spin0, PDT::Spin0);
	for(size_t i = 0; i < numberOfDiags(); ++i ) {
	HPDiagram dg = getProcessInfo()[i];
	if( !dg.intermediate ) {
	contact_ = dynamic_ptr_cast<AbstractVVSSVertexPtr>(dg.vertices.first);
	}
	else if(dg.channelType == HPDiagram::tChannel) {
	- AbstractVSSVertexPtr vss1 =
	- dynamic_ptr_cast<AbstractVSSVertexPtr>(dg.vertices.first);
	- AbstractVSSVertexPtr vss2 =
	- dynamic_ptr_cast<AbstractVSSVertexPtr>(dg.vertices.second);
	- scalar2_[i] = make_pair(vss1, vss2);
	+ if (dg.intermediate->iSpin() == PDT::Spin0 ) {
	+ AbstractVSSVertexPtr vss1 =
	+ dynamic_ptr_cast<AbstractVSSVertexPtr>(dg.vertices.first);
	+ AbstractVSSVertexPtr vss2 =
	+ dynamic_ptr_cast<AbstractVSSVertexPtr>(dg.vertices.second);
	+ scalar2_[i] = make_pair(vss1, vss2);
	+ }
	+ else if( dg.intermediate->iSpin() == PDT::Spin1 ) {
	+ AbstractVVSVertexPtr vvs1 =
	+ dynamic_ptr_cast<AbstractVVSVertexPtr>(dg.vertices.first);
	+ AbstractVVSVertexPtr vvs2 =
	+ dynamic_ptr_cast<AbstractVVSVertexPtr>(dg.vertices.second);
	+ scalar3_[i] = make_pair(vvs1, vvs2);
	+ }
	+ else
	+ assert(false);
	}
	else {
	if( dg.intermediate->iSpin() == PDT::Spin0 ) {
	AbstractVVSVertexPtr vvs =
	dynamic_ptr_cast<AbstractVVSVertexPtr>(dg.vertices.first);
	AbstractSSSVertexPtr sss =
	dynamic_ptr_cast<AbstractSSSVertexPtr>(dg.vertices.second);
	scalar1_[i] = make_pair(vvs, sss);
	}
	else if( dg.intermediate->iSpin() == PDT::Spin1 ) {
	AbstractVVVVertexPtr vvv =
	dynamic_ptr_cast<AbstractVVVVertexPtr>(dg.vertices.first);
	AbstractVSSVertexPtr vss =
	dynamic_ptr_cast<AbstractVSSVertexPtr>(dg.vertices.second);
	vector_[i] = make_pair(vvv, vss);
	}
	else if( dg.intermediate->iSpin() == PDT::Spin2 ) {
	AbstractVVTVertexPtr vvt =
	dynamic_ptr_cast<AbstractVVTVertexPtr>(dg.vertices.first);
	AbstractSSTVertexPtr sst =
	dynamic_ptr_cast<AbstractSSTVertexPtr>(dg.vertices.second);
	tensor_[i] = make_pair(vvt, sst);
	}
	}
	}
	}

	double MEvv2ss::me2() const {
	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);
	}
	ScalarWaveFunction sca1(rescaledMomenta()[2],mePartonData()[2],
	Complex(1.,0.),outgoing);
	ScalarWaveFunction sca2(rescaledMomenta()[3],mePartonData()[3],
	Complex(1.,0.),outgoing);
	double full_me(0.);
	vv2ssME(v1, v2, sca1, sca2, full_me , true);

	#ifndef NDEBUG
	if( debugME() ) debug(full_me);
	#endif

	return full_me;
	}

	ProductionMatrixElement
	MEvv2ss::vv2ssME(const VBVector & v1, const VBVector & v2,
	const ScalarWaveFunction & sca1,
	const ScalarWaveFunction & sca2,
	double & me2, bool first) const {
	const Energy2 m2(scale());
	const Energy masst = sca1.mass(), massu = sca2.mass();
	// 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.);
	//loop over vector helicities
	for(unsigned int iv1 = 0; iv1 < 2; ++iv1) {
	for(unsigned int iv2 = 0; iv2 < 2; ++iv2) {
	vector<Complex> flows(numberOfFlows(),0.);
	// loop over diagrams
	for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
	Complex diag(0.);
	const HPDiagram & current = getProcessInfo()[ix];
	// do four-point diag first
	if(current.channelType == HPDiagram::fourPoint) {
	diag = contact_->evaluate(m2, v1[iv1], v2[iv2], sca1, sca2);
	}
	else {
	tcPDPtr offshell = current.intermediate;
	if(current.channelType == HPDiagram::tChannel) {
	- if(current.ordered.second) {
	- ScalarWaveFunction interS = scalar2_[ix].first->
	- evaluate(m2, 3, offshell, v1[iv1], sca1, masst);
	- diag = scalar2_[ix].second->evaluate(m2, v2[iv2], interS, sca2);
	+ if(offshell->iSpin() == PDT::Spin0) {
	+ if(current.ordered.second) {
	+ ScalarWaveFunction interS = scalar2_[ix].first->
	+ evaluate(m2, 3, offshell, v1[iv1], sca1, masst);
	+ diag = scalar2_[ix].second->evaluate(m2, v2[iv2], interS, sca2);
	+ }
	+ else {
	+ ScalarWaveFunction interS = scalar2_[ix].first->
	+ evaluate(m2, 3, offshell, v1[iv1], sca2, massu);
	+ diag = scalar2_[ix].second->evaluate(m2, v2[iv2], interS, sca1);
	+ }
	}
	else {
	- ScalarWaveFunction interS = scalar2_[ix].first->
	- evaluate(m2, 3, offshell, v1[iv1], sca2, massu);
	- diag = scalar2_[ix].second->evaluate(m2, v2[iv2], interS, sca1);
	+ if(current.ordered.second) {
	+ VectorWaveFunction interV = scalar3_[ix].first->
	+ evaluate(m2, 3, offshell, v1[iv1], sca1);
	+ diag = scalar3_[ix].second->evaluate(m2, v2[iv2], interV, sca2);
	+ }
	+ else {
	+ VectorWaveFunction interV = scalar3_[ix].first->
	+ evaluate(m2, 3, offshell, v1[iv1], sca2);
	+ diag = scalar3_[ix].second->evaluate(m2, v2[iv2], interV, sca1);
	+ }
	}
	}
	else if(current.channelType == HPDiagram::sChannel) {
	if(offshell->iSpin() == PDT::Spin0) {
	ScalarWaveFunction interS = scalar1_[ix].first->
	evaluate(m2, 1, offshell, v1[iv1], v2[iv2]);
	diag = scalar1_[ix].second->evaluate(m2, interS, sca1, sca2);
	}
	else if(offshell->iSpin() == PDT::Spin1) {
	VectorWaveFunction interV = vector_[ix].first->
	evaluate(m2, 1, offshell, v1[iv1], v2[iv2]);
	diag = vector_[ix].second->evaluate(m2, interV, sca1, sca2);
	}
	else if(offshell->iSpin() == PDT::Spin2) {
	TensorWaveFunction interT = tensor_[ix].first->
	evaluate(m2, 1, offshell, v1[iv1], v2[iv2]);
	diag = tensor_[ix].second->evaluate(m2, sca1, sca2, interT);
	}
	}
	else
	diag = 0.;
	}
	me[ix] += norm(diag);
	diagramME()[ix](2iv1, 2iv2, 0, 0) = 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, 0, 0) = 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 MEvv2ss::persistentOutput(PersistentOStream & os) const {
	- os << scalar1_ << scalar2_ << vector_ << tensor_ << contact_;
	+ os << scalar1_ << scalar2_ << scalar3_ << vector_ << tensor_ << contact_;
	}

	void MEvv2ss::persistentInput(PersistentIStream & is, int) {
	- is >> scalar1_ >> scalar2_ >> vector_ >> tensor_ >> contact_;
	+ is >> scalar1_ >> scalar2_ >> scalar3_ >> vector_ >> tensor_ >> contact_;
	initializeMatrixElements(PDT::Spin1, PDT::Spin1,
	PDT::Spin0, PDT::Spin0);
	}

	ClassDescription<MEvv2ss> MEvv2ss::initMEvv2ss;
	// Definition of the static class description member.

	void MEvv2ss::Init() {

	static ClassDocumentation<MEvv2ss> documentation
	("This class implements the ME for the vector-vector to scalar-scalar "
	"hard-process");

	}

	void MEvv2ss::constructVertex(tSubProPtr sub) {
	ParticleVector ext = hardParticles(sub);
	VBVector v1, v2;
	// set up the wavefunctions with real momenta
	VectorWaveFunction(v1, ext[0], incoming, false, true);
	VectorWaveFunction(v2, ext[1], incoming, false, true);
	ScalarWaveFunction sca1(ext[2], outgoing, true);
	ScalarWaveFunction sca2(ext[3], outgoing, true);
	// calculate rescaled moment
	setRescaledMomenta(ext);
	// wavefunctions with rescaled momenta
	VectorWaveFunction v1r (rescaledMomenta()[0],
	ext[0]->dataPtr(), incoming);
	VectorWaveFunction v2r (rescaledMomenta()[1],
	ext[1]->dataPtr(), incoming);
	sca1 = ScalarWaveFunction(rescaledMomenta()[2],
	ext[2]->dataPtr(), outgoing);
	sca2 = ScalarWaveFunction(rescaledMomenta()[3],
	ext[3]->dataPtr(), outgoing);
	for( unsigned int ihel = 0; ihel < 2; ++ihel ) {
	v1r.reset(2*ihel);
	v1[ihel] = v1r;
	v2r.reset(2*ihel);
	v2[ihel] = v2r;
	}
	double dummy(0.);
	ProductionMatrixElement pme = vv2ssME(v1, v2, sca1, sca2, dummy , false);

	#ifndef NDEBUG
	if( debugME() ) debug(dummy);
	#endif

	createVertex(pme,ext);
	}

	void MEvv2ss::debug(double me2) const {
	if( !generator()->logfile().is_open() ) return;
	if( mePartonData()[0]->id() != 21 \|\| mePartonData()[1]->id() != 21) return;
	long id3 = abs(mePartonData()[2]->id());
	long id4 = abs(mePartonData()[3]->id());
	int type = -1;
	//SUSY gg>~q~q
	if( ((id3 >= 1000001 && id3 <= 1000006 ) && (id4 >= 1000001 && id4 <= 1000006 ) ) \|\|
	((id3 >= 2000001 && id3 <= 2000006 ) && (id4 >= 2000001 && id4 <= 2000006 ) ) ) {
	type = 0;
	}
	// Sextet production
	else if(mePartonData()[2]->iColour() == PDT::Colour6 &&
	mePartonData()[3]->iColour() == PDT::Colour6bar ) {
	type = 1;
	}
	else {
	return;
	}
	double gs4 = sqr( 4.Constants::piSM().alphaS(scale()));
	int Nc = SM().Nc();
	Energy4 s2 = sqr(sHat());
	Energy2 m3s = meMomenta()[2].m2();
	Energy2 m4s = meMomenta()[3].m2();
	Energy4 spt2 = uHat()tHat() - m3sm4s;
	Energy2 t3 = tHat()-m3s, u4 = uHat()-m4s;
	Energy4 t3s = sqr(t3) , u4s = sqr(u4);
	Energy8 pre = gs4(sqr(spt2) + s2m3s*m4s);
	// matrix element
	double analytic(0.);
	// triplet scalars
	if(type==0) {
	analytic = preNc
	( u4s + t3s - s2/sqr(Nc) )/2./(sqr(Nc) - 1.)/s2/t3s/u4s;
	}
	// sextet scalars
	else if(type==1) {
	analytic = pre(Nc+2.)/(sqr(Nc)-1.)/Nc
	((Nc+2.)*(Nc-1.)/t3s/u4s - sqr(Nc)/t3/u4/s2);
	}
	double diff = abs(analytic - me2)/(analytic+me2);
	if( diff > 1e-10 ) {
	generator()->log()
	<< mePartonData()[0]->PDGName() << ","
	<< mePartonData()[1]->PDGName() << "->"
	<< mePartonData()[2]->PDGName() << ","
	<< mePartonData()[3]->PDGName() << " difference: "
	<< setprecision(10) << diff << " ratio: " << analytic/me2 << '\n';
	}
	}
	diff --git a/MatrixElement/General/MEvv2ss.h b/MatrixElement/General/MEvv2ss.h
	--- a/MatrixElement/General/MEvv2ss.h
	+++ b/MatrixElement/General/MEvv2ss.h
	@@ -1,227 +1,232 @@
	// -- C++ --
	//
	// MEvv2ss.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_MEvv2ss_H
	#define HERWIG_MEvv2ss_H
	//
	// This is the declaration of the MEvv2ss class.
	//

	#include "GeneralHardME.h"
	#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractVVTVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractSSTVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractVVSSVertex.h"
	#include "ThePEG/Helicity/Vertex/AbstractSSSVertex.h"
	#include "Herwig/MatrixElement/ProductionMatrixElement.h"

	namespace Herwig {
	using namespace ThePEG;
	using ThePEG::Helicity::VectorWaveFunction;
	using ThePEG::Helicity::ScalarWaveFunction;

	/**
	* This is the implementation of the matrix element for the process
	* vector-vector to scalar-scalar. It inherits from GeneralHardME and
	* implements the required virtual functions.
	*
	* @see \ref MEff2ffInterfaces "The Interfaces"
	* defined for MEff2ff.
	* @see GeneralHardME
	*/
	class MEvv2ss: public GeneralHardME {

	public:

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

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

	/**
	* Set the Hardvertex for the spin correlations
	* @param sub
	*/
	virtual void constructVertex(tSubProPtr sub);

	private:

	/**
	* Calculate the matrix element.
	* @param v1 A vector of VectorWaveFunction objects for the first boson
	* @param v2 A vector of VectorWaveFunction objects for the second boson
	* @param sca1 A ScalarWaveFunction for the first outgoing
	* @param sca2 A ScalarWaveFunction for the second outgoing
	* @param me2 The value of the spin-summed matrix element squared
	* (to be calculated)
	* @param first Whether or not first call to decide if colour decomposition etc
	* should be calculated
	*/
	ProductionMatrixElement vv2ssME(const VBVector & v1, const VBVector & v2,
	const ScalarWaveFunction & sca1,
	const ScalarWaveFunction & sca2,
	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 static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<MEvv2ss> initMEvv2ss;

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

	private:

	/** @name The dynamically casted vertices. */
	//@{
	/**
	* Intermediate s-channel scalar
	*/
	vector<pair<AbstractVVSVertexPtr, AbstractSSSVertexPtr> > scalar1_;

	/**
	* Intermediate t-channel scalar
	*/
	vector<pair<AbstractVSSVertexPtr, AbstractVSSVertexPtr> > scalar2_;

	/**
	+ * Intermediate t-channel scalar
	+ */
	+ vector<pair<AbstractVVSVertexPtr, AbstractVVSVertexPtr> > scalar3_;
	+
	+ /**
	* Intermediate s-channel vector
	*/
	vector<pair<AbstractVVVVertexPtr, AbstractVSSVertexPtr> > vector_;

	/**
	* Intermediate s-channel tensor
	*/
	vector<pair<AbstractVVTVertexPtr, AbstractSSTVertexPtr> > tensor_;

	/**
	* The contact vertex
	*/
	AbstractVVSSVertexPtr contact_;
	//@}


	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of MEvv2ss. */
	template <>
	struct BaseClassTrait<Herwig::MEvv2ss,1> {
	/** Typedef of the first base class of MEvv2ss. */
	typedef Herwig::GeneralHardME NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the MEvv2ss class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::MEvv2ss>
	: public ClassTraitsBase<Herwig::MEvv2ss> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::MEvv2ss"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_MEvv2ss_H */
	diff --git a/Models/General/HardProcessConstructor.cc b/Models/General/HardProcessConstructor.cc
	--- a/Models/General/HardProcessConstructor.cc
	+++ b/Models/General/HardProcessConstructor.cc
	@@ -1,647 +1,651 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the HardProcessConstructor class.
	//

	#include "HardProcessConstructor.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace Herwig;

	void HardProcessConstructor::persistentOutput(PersistentOStream & os) const {
	os << debug_ << subProcess_ << model_;
	}

	void HardProcessConstructor::persistentInput(PersistentIStream & is, int) {
	is >> debug_ >> subProcess_ >> model_;
	}

	AbstractClassDescription<HardProcessConstructor>
	HardProcessConstructor::initHardProcessConstructor;
	// Definition of the static class description member.

	void HardProcessConstructor::Init() {

	static ClassDocumentation<HardProcessConstructor> documentation
	("Base class for implementation of the automatic generation of hard processes");

	static Switch<HardProcessConstructor,bool> interfaceDebugME
	("DebugME",
	"Print comparison with analytical ME",
	&HardProcessConstructor::debug_, false, false, false);
	static SwitchOption interfaceDebugMEYes
	(interfaceDebugME,
	"Yes",
	"Print the debug information",
	true);
	static SwitchOption interfaceDebugMENo
	(interfaceDebugME,
	"No",
	"Do not print the debug information",
	false);

	}

	void HardProcessConstructor::doinit() {
	Interfaced::doinit();
	EGPtr eg = generator();
	model_ = dynamic_ptr_cast<HwSMPtr>(eg->standardModel());
	if(!model_)
	throw InitException() << "HardProcessConstructor:: doinit() - "
	<< "The model pointer is null!"
	<< Exception::abortnow;
	if(!eg->eventHandler()) {
	throw
	InitException() << "HardProcessConstructor:: doinit() - "
	<< "The eventHandler pointer was null therefore "
	<< "could not get SubProcessHandler pointer "
	<< Exception::abortnow;
	}
	string subProcessName =
	eg->preinitInterface(eg->eventHandler(), "SubProcessHandlers", "get","");
	subProcess_ = eg->getObject<SubProcessHandler>(subProcessName);
	if(!subProcess_) {
	ostringstream s;
	s << "HardProcessConstructor:: doinit() - "
	<< "There was an error getting the SubProcessHandler "
	<< "from the current event handler. ";
	generator()->logWarning( Exception(s.str(), Exception::warning) );
	}
	}

	GeneralHardME::ColourStructure HardProcessConstructor::
	colourFlow(const tcPDVector & extpart) const {
	PDT::Colour ina = extpart[0]->iColour();
	PDT::Colour inb = extpart[1]->iColour();
	PDT::Colour outa = extpart[2]->iColour();
	PDT::Colour outb = extpart[3]->iColour();

	// incoming colour neutral
	if(ina == PDT::Colour0 && inb == PDT::Colour0) {
	if( outa == PDT::Colour0 && outb == PDT::Colour0 ) {
	return GeneralHardME::Colour11to11;
	}
	else if( outa == PDT::Colour3 && outb == PDT::Colour3bar ) {
	return GeneralHardME::Colour11to33bar;
	}
	else if( outa == PDT::Colour8 && outb == PDT::Colour8 ) {
	return GeneralHardME::Colour11to88;
	}
	else
	assert(false);
	}
	// incoming 3 3
	else if(ina == PDT::Colour3 && inb == PDT::Colour3 ) {
	if( outa == PDT::Colour3 && outb == PDT::Colour3 ) {
	return GeneralHardME::Colour33to33;
	}
	else if( outa == PDT::Colour6 && outb == PDT::Colour0 ) {
	return GeneralHardME::Colour33to61;
	}
	else if( outa == PDT::Colour0 && outb == PDT::Colour6 ) {
	return GeneralHardME::Colour33to16;
	}
	else
	assert(false);
	}
	// incoming 3bar 3bar
	else if(ina == PDT::Colour3bar && inb == PDT::Colour3bar ) {
	if( outa == PDT::Colour3bar && outb == PDT::Colour3bar ) {
	return GeneralHardME::Colour3bar3barto3bar3bar;
	}
	else if( outa == PDT::Colour6bar && outb == PDT::Colour0) {
	return GeneralHardME::Colour3bar3barto6bar1;
	}
	else if ( outa == PDT::Colour0 && outb == PDT::Colour6bar ) {
	return GeneralHardME::Colour3bar3barto16bar;
	}
	else
	assert(false);
	}
	// incoming 3 3bar
	else if(ina == PDT::Colour3 && inb == PDT::Colour3bar ) {
	if( outa == PDT::Colour0 && outb == PDT::Colour0 ) {
	return GeneralHardME::Colour33barto11;
	}
	else if( outa == PDT::Colour3 && outb == PDT::Colour3bar ) {
	return GeneralHardME::Colour33barto33bar;
	}
	else if( outa == PDT::Colour8 && outb == PDT::Colour8 ) {
	return GeneralHardME::Colour33barto88;
	}
	else if( outa == PDT::Colour8 && outb == PDT::Colour0 ) {
	return GeneralHardME::Colour33barto81;
	}
	else if( outa == PDT::Colour0 && outb == PDT::Colour8 ) {
	return GeneralHardME::Colour33barto18;
	}
	else if( outa == PDT::Colour6 && outb == PDT::Colour6bar) {
	return GeneralHardME::Colour33barto66bar;
	}
	else if( outa == PDT::Colour6bar && outb == PDT::Colour6) {
	return GeneralHardME::Colour33barto6bar6;
	}
	else
	assert(false);
	}
	// incoming 88
	else if(ina == PDT::Colour8 && inb == PDT::Colour8 ) {
	if( outa == PDT::Colour0 && outb == PDT::Colour0 ) {
	return GeneralHardME::Colour88to11;
	}
	else if( outa == PDT::Colour3 && outb == PDT::Colour3bar ) {
	return GeneralHardME::Colour88to33bar;
	}
	else if( outa == PDT::Colour8 && outb == PDT::Colour8 ) {
	return GeneralHardME::Colour88to88;
	}
	else if( outa == PDT::Colour8 && outb == PDT::Colour0 ) {
	return GeneralHardME::Colour88to81;
	}
	else if( outa == PDT::Colour0 && outb == PDT::Colour8 ) {
	return GeneralHardME::Colour88to18;
	}
	else if( outa == PDT::Colour6 && outb == PDT::Colour6bar ) {
	return GeneralHardME::Colour88to66bar;
	}
	else
	assert(false);
	}
	// incoming 38
	else if(ina == PDT::Colour3 && inb == PDT::Colour8 ) {
	if(outa == PDT::Colour3 && outb == PDT::Colour0) {
	return GeneralHardME::Colour38to31;
	}
	else if(outa == PDT::Colour0 && outb == PDT::Colour3) {
	return GeneralHardME::Colour38to13;
	}
	else if(outa == PDT::Colour3 && outb == PDT::Colour8) {
	return GeneralHardME::Colour38to38;
	}
	else if(outa == PDT::Colour8 && outb == PDT::Colour3) {
	return GeneralHardME::Colour38to83;
	}
	else if(outa == PDT::Colour3bar && outb == PDT::Colour6){
	return GeneralHardME::Colour38to3bar6;
	}
	else if(outa == PDT::Colour6 && outb == PDT::Colour3bar){
	return GeneralHardME::Colour38to63bar;
	}
	else
	assert(false);
	}
	// incoming 3bar8
	else if(ina == PDT::Colour3bar && inb == PDT::Colour8 ) {
	if(outa == PDT::Colour3bar && outb == PDT::Colour0 ) {
	return GeneralHardME::Colour3bar8to3bar1;
	}
	else if(outa == PDT::Colour0 && outb == PDT::Colour3bar) {
	return GeneralHardME::Colour3bar8to13bar;
	}
	else if(outa == PDT::Colour3bar && outb == PDT::Colour8 ) {
	return GeneralHardME::Colour3bar8to3bar8;
	}
	else if(outa == PDT::Colour8 && outb == PDT::Colour3bar) {
	return GeneralHardME::Colour3bar8to83bar;
	}
	else
	assert(false);
	}
	// unknown colour flow
	else
	assert(false);
	return GeneralHardME::UNDEFINED;
	}


	void HardProcessConstructor::fixFSOrder(HPDiagram & diag) {
	tcPDPtr psa = getParticleData(diag.incoming.first);
	tcPDPtr psb = getParticleData(diag.incoming.second);
	tcPDPtr psc = getParticleData(diag.outgoing.first);
	tcPDPtr psd = getParticleData(diag.outgoing.second);

	//fix a spin order
	if( psc->iSpin() < psd->iSpin() ) {
	swap(diag.outgoing.first, diag.outgoing.second);
	if(diag.channelType == HPDiagram::tChannel) {
	diag.ordered.second = !diag.ordered.second;
	}
	return;
	}

	if( psc->iSpin() == psd->iSpin() &&
	psc->id() < 0 && psd->id() > 0 ) {
	swap(diag.outgoing.first, diag.outgoing.second);
	if(diag.channelType == HPDiagram::tChannel) {
	diag.ordered.second = !diag.ordered.second;
	}
	return;
	}
	}

	void HardProcessConstructor::assignToCF(HPDiagram & diag) {
	if(diag.channelType == HPDiagram::tChannel) {
	if(diag.ordered.second) tChannelCF(diag);
	else uChannelCF(diag);
	}
	else if(diag.channelType == HPDiagram::sChannel) {
	sChannelCF(diag);
	}
	else if (diag.channelType == HPDiagram::fourPoint) {
	fourPointCF(diag);
	}
	else
	assert(false);
	}

	void HardProcessConstructor::tChannelCF(HPDiagram & diag) {
	PDT::Colour ina = getParticleData(diag.incoming.first )->iColour();
	PDT::Colour inb = getParticleData(diag.incoming.second)->iColour();
	PDT::Colour outa = getParticleData(diag.outgoing.first )->iColour();
	PDT::Colour outb = getParticleData(diag.outgoing.second)->iColour();
	vector<CFPair> cfv(1, make_pair(0, 1.));
	if(diag.intermediate->iColour() == PDT::Colour0) {
	if(ina==PDT::Colour0) {
	cfv[0] = make_pair(0, 1);
	}
	else if(ina==PDT::Colour3 \|\| ina==PDT::Colour3bar) {
	if( inb == PDT::Colour0 ) {
	cfv[0] = make_pair(0, 1);
	}
	else if(inb==PDT::Colour3 \|\| outb==PDT::Colour3bar) {
	cfv[0] = make_pair(2, 1);
	}
	else if(inb==PDT::Colour8) {
	cfv[0] = make_pair(2, 1);
	}
	}
	else if(ina==PDT::Colour8) {
	if( inb == PDT::Colour0 ) {
	cfv[0] = make_pair(0, 1);
	}
	else if(inb==PDT::Colour3 \|\| outb==PDT::Colour3bar) {
	cfv[0] = make_pair(2, 1);
	}
	else if(inb==PDT::Colour8) {
	cfv[0] = make_pair(5, 1);
	}
	}
	}
	else if(diag.intermediate->iColour() == PDT::Colour8) {
	if(ina==PDT::Colour8&&outa==PDT::Colour8&&
	inb==PDT::Colour8&&outb==PDT::Colour8) {
	cfv.push_back(make_pair(1, -1.));
	}
	}
	else if(diag.intermediate->iColour() == PDT::Colour3 \|\|
	diag.intermediate->iColour() == PDT::Colour3bar) {
	if(outa == PDT::Colour0 \|\| outb == PDT::Colour0) {
	if( outa != PDT::Colour6 && outb != PDT::Colour6 &&
	outa != PDT::Colour6bar && outb != PDT::Colour6bar) {
	cfv[0] = make_pair(0,1.);
	}
	else {
	cfv[0] = make_pair(0,0.5);
	cfv.push_back(make_pair(1,0.5));
	}
	}
	else if(outa==PDT::Colour6 && outb==PDT::Colour3bar) {
	cfv[0] = make_pair(4,1.);
	cfv.push_back(make_pair(5,1.));
	}
	else if(outa==PDT::Colour6 && outb==PDT::Colour6bar) {
	cfv[0] = make_pair(4, 1.);
	for(unsigned int ix=5;ix<8;++ix)
	cfv.push_back(make_pair(ix,1.));
	}
	else if(outa==PDT::Colour6 \|\| outa ==PDT::Colour6bar \|\|
	outb==PDT::Colour6 \|\| outb ==PDT::Colour6bar ) {
	assert(false);
	}
	}
	else if(diag.intermediate->iColour() == PDT::Colour6 \|\|
	diag.intermediate->iColour() == PDT::Colour6bar) {
	if(ina==PDT::Colour8 && inb==PDT::Colour8) {
	cfv[0] = make_pair(0, 1.);
	for(unsigned int ix=1;ix<4;++ix)
	cfv.push_back(make_pair(ix,1.));
	for(unsigned int ix=4;ix<8;++ix)
	cfv.push_back(make_pair(ix,1.));
	}
	else if(outa==PDT::Colour3bar && outb==PDT::Colour6) {
	cfv[0] = make_pair(0,1.);
	for(unsigned int ix=1;ix<4;++ix)
	cfv.push_back(make_pair(ix,1.));
	}
	else if(outa==PDT::Colour6 && outb==PDT::Colour3bar) {
	cfv[0] = make_pair(4,1.);
	cfv.push_back(make_pair(5,1.));
	}
	}
	diag.colourFlow = cfv;
	}

	void HardProcessConstructor::uChannelCF(HPDiagram & diag) {
	PDT::Colour offshell = diag.intermediate->iColour();
	PDT::Colour ina = getParticleData(diag.incoming.first )->iColour();
	PDT::Colour inb = getParticleData(diag.incoming.second)->iColour();
	PDT::Colour outa = getParticleData(diag.outgoing.first )->iColour();
	PDT::Colour outb = getParticleData(diag.outgoing.second)->iColour();
	vector<CFPair> cfv(1, make_pair(1, 1.));
	if(offshell == PDT::Colour8) {
	- if( outa != outb ) {
	+ if(outa == PDT::Colour0 &&
	+ outb == PDT::Colour0) {
	+ cfv[0].first = 0;
	+ }
	+ else if( outa != outb ) {
	if(outa == PDT::Colour0 \|\|
	outb == PDT::Colour0) {
	cfv[0].first = 0;
	}
	else {
	cfv[0].first = 0;
	cfv.push_back(make_pair(1, -1.));
	}
	}
	else if(outa==PDT::Colour8&&ina==PDT::Colour8) {
	cfv[0]=make_pair(1, -1.);
	cfv.push_back(make_pair(2, 1.));
	}
	}
	else if(offshell == PDT::Colour3 \|\| offshell == PDT::Colour3bar) {
	if( outa == PDT::Colour0 \|\| outb == PDT::Colour0 ) {
	if( outa != PDT::Colour6 && outb != PDT::Colour6 &&
	outa != PDT::Colour6bar && outb != PDT::Colour6bar) {
	cfv[0] = make_pair(0,1.);
	}
	else {
	cfv[0] = make_pair(0,0.5);
	cfv.push_back(make_pair(1,0.5));
	}
	}
	else if(outa==PDT::Colour3bar && outb==PDT::Colour6) {
	cfv[0] = make_pair(4,1.);
	cfv.push_back(make_pair(5,1.));
	}
	else if(outa==PDT::Colour6 && outb==PDT::Colour3bar) {
	cfv[0] = make_pair(0,1.);
	for(int ix=0; ix<4;++ix)
	cfv.push_back(make_pair(ix,1.));
	}
	else if(outa==PDT::Colour6bar && outb==PDT::Colour6) {
	cfv[0] = make_pair(4,1.);
	for(int ix=5; ix<8;++ix)
	cfv.push_back(make_pair(ix,1.));
	}
	}
	else if( offshell == PDT::Colour0 ) {
	if(ina==PDT::Colour0) {
	cfv[0] = make_pair(0, 1);
	}
	else if(ina==PDT::Colour3 \|\| ina==PDT::Colour3bar) {
	if( inb == PDT::Colour0 ) {
	cfv[0] = make_pair(0, 1);
	}
	else if(inb==PDT::Colour3 \|\| outb==PDT::Colour3bar) {
	cfv[0] = make_pair(3, 1);
	}
	else if(inb==PDT::Colour8) {
	cfv[0] = make_pair(2, 1);
	}
	}
	else if(ina==PDT::Colour8) {
	if( inb == PDT::Colour0 ) {
	cfv[0] = make_pair(0, 1);
	}
	else if(inb==PDT::Colour3 \|\| outb==PDT::Colour3bar) {
	cfv[0] = make_pair(2, 1);
	}
	else if(inb==PDT::Colour8) {
	cfv[0] = make_pair(4, 1);
	}
	}
	}
	else if(diag.intermediate->iColour() == PDT::Colour6 \|\|
	diag.intermediate->iColour() == PDT::Colour6bar) {
	if(ina==PDT::Colour8 && inb==PDT::Colour8) {
	cfv[0] = make_pair(0, 1.);
	for(unsigned int ix=1;ix<4;++ix)
	cfv.push_back(make_pair(ix,1.));
	for(unsigned int ix=8;ix<12;++ix)
	cfv.push_back(make_pair(ix,1.));
	}
	else if(outa==PDT::Colour3bar && outa==PDT::Colour6) {
	cfv[0] = make_pair(4, 1.);
	cfv.push_back(make_pair(5,1.));
	}
	else if(outa==PDT::Colour6 && outa==PDT::Colour3bar) {
	cfv[0] = make_pair(0, 1.);
	for(unsigned int ix=1;ix<4;++ix)
	cfv.push_back(make_pair(ix,1.));
	}
	}
	diag.colourFlow = cfv;
	}

	void HardProcessConstructor::sChannelCF(HPDiagram & diag) {
	tcPDPtr pa = getParticleData(diag.incoming.first);
	tcPDPtr pb = getParticleData(diag.incoming.second);
	PDT::Colour ina = pa->iColour();
	PDT::Colour inb = pb->iColour();
	PDT::Colour offshell = diag.intermediate->iColour();
	tcPDPtr pc = getParticleData(diag.outgoing.first);
	tcPDPtr pd = getParticleData(diag.outgoing.second);
	PDT::Colour outa = pc->iColour();
	PDT::Colour outb = pd->iColour();
	vector<CFPair> cfv(1);
	if(offshell == PDT::Colour8) {
	if(ina == PDT::Colour0 \|\| inb == PDT::Colour0 \|\|
	outa == PDT::Colour0 \|\| outb == PDT::Colour0) {
	cfv[0] = make_pair(0, 1);
	}
	else {
	bool incol = ina == PDT::Colour8 && inb == PDT::Colour8;
	bool outcol = outa == PDT::Colour8 && outb == PDT::Colour8;
	bool intrip = ina == PDT::Colour3 && inb == PDT::Colour3bar;
	bool outtrip = outa == PDT::Colour3 && outb == PDT::Colour3bar;
	bool outsex = outa == PDT::Colour6 && outb == PDT::Colour6bar;
	bool outsexb = outa == PDT::Colour6bar && outb == PDT::Colour6;
	if(incol \|\| outcol) {
	// Require an additional minus sign for a scalar/fermion
	// 33bar final state due to the way the vertex rules are defined.
	int prefact(1);
	if( ((pc->iSpin() == PDT::Spin1Half && pd->iSpin() == PDT::Spin1Half) \|\|
	(pc->iSpin() == PDT::Spin0 && pd->iSpin() == PDT::Spin0 )) &&
	(outa == PDT::Colour3 && outb == PDT::Colour3bar) )
	prefact = -1;
	if(incol && outcol) {
	cfv[0].first = 0;
	cfv[0].second = -prefact;
	cfv.push_back(make_pair(2, prefact));
	}
	else if(incol && outsex) {
	cfv[0].first = 4;
	cfv[0].second = prefact;
	for(unsigned int ix=1;ix<4;++ix)
	cfv.push_back(make_pair(4+ix, prefact));
	for(unsigned int ix=0;ix<4;++ix)
	cfv.push_back(make_pair(8+ix,-prefact));
	}
	else {
	cfv[0].first = 0;
	cfv[0].second = -prefact;
	cfv.push_back(make_pair(1, prefact));
	}
	}
	else if( ( intrip && !outtrip ) \|\|
	( !intrip && outtrip ) ) {
	if(!outsex)
	cfv[0] = make_pair(0, 1);
	else {
	cfv[0] = make_pair(0, 1.);
	for(unsigned int ix=0;ix<3;++ix)
	cfv.push_back(make_pair(ix+1, 1.));
	}
	}
	else if((intrip && outsex) \|\| (intrip && outsexb)) {
	cfv[0] = make_pair(0,1.);
	for(int ix=1; ix<4; ++ix)
	cfv.push_back(make_pair(ix,1.));
	}
	else
	cfv[0] = make_pair(1, 1);
	}
	}
	else if(offshell == PDT::Colour0) {
	if( ina == PDT::Colour0 ) {
	cfv[0] = make_pair(0, 1);
	}
	else if(ina==PDT::Colour3 \|\| ina==PDT::Colour3bar) {
	if( outa == PDT::Colour0 ) {
	cfv[0] = make_pair(0, 1);
	}
	else if(outa==PDT::Colour3 \|\| outa==PDT::Colour3bar) {
	cfv[0] = make_pair(3, 1);
	}
	else if(outa==PDT::Colour8) {
	cfv[0] = make_pair(2, 1);
	}
	else
	assert(false);
	}
	else if(ina==PDT::Colour8) {
	if( outa == PDT::Colour0 ) {
	cfv[0] = make_pair(0, 1);
	}
	else if(outa==PDT::Colour3 \|\| outb==PDT::Colour3bar) {
	cfv[0] = make_pair(2, 1);
	}
	else if(outa==PDT::Colour8) {
	cfv[0] = make_pair(3, 1);
	}
	}
	}
	else if(offshell == PDT::Colour3 \|\| offshell == PDT::Colour3bar) {
	if(outa == PDT::Colour6 \|\| outa == PDT::Colour6bar \|\|
	outb == PDT::Colour6bar \|\| outb == PDT::Colour6) {
	cfv[0] = make_pair(6, 1.);
	cfv.push_back(make_pair(7,1.));
	}
	else {
	if(outa == PDT::Colour0 \|\| outb == PDT::Colour0)
	cfv[0] = make_pair(0, 1);
	else
	cfv[0] = make_pair(1, 1);
	}
	}
	else if( offshell == PDT::Colour6 \|\| offshell == PDT::Colour6bar) {
	if((ina == PDT::Colour3 && inb == PDT::Colour3 &&
	outa == PDT::Colour3 && outb == PDT::Colour3 ) \|\|
	(ina == PDT::Colour3bar && inb == PDT::Colour3bar &&
	outa == PDT::Colour3bar && outb == PDT::Colour3bar)) {
	cfv[0] = make_pair(2,0.5);
	cfv.push_back(make_pair(3,0.5));
	}
	else if((ina == PDT::Colour3 && inb == PDT::Colour3 &&
	((outa == PDT::Colour6 && outb == PDT::Colour0)\|\|
	(outb == PDT::Colour6 && outa == PDT::Colour0))) \|\|
	(ina == PDT::Colour3bar && inb == PDT::Colour3bar &&
	((outa == PDT::Colour6bar && outb == PDT::Colour0)\|\|
	(outb == PDT::Colour6bar && outa == PDT::Colour0)))) {
	cfv[0] = make_pair(0,0.5);
	cfv.push_back(make_pair(1,0.5));
	}
	else
	assert(false);
	}
	else {
	if(outa == PDT::Colour0 \|\| outb == PDT::Colour0)
	cfv[0] = make_pair(0, 1);
	else
	cfv[0] = make_pair(1, 1);
	}
	diag.colourFlow = cfv;
	}

	void HardProcessConstructor::fourPointCF(HPDiagram & diag) {
	// count the colours
	unsigned int noct(0),ntri(0),nsng(0),nsex(0);
	for(unsigned int ix=0;ix<4;++ix) {
	PDT::Colour col = getParticleData(diag.ids[ix])->iColour();
	if(col==PDT::Colour0) ++nsng;
	else if(col==PDT::Colour3\|\|col==PDT::Colour3bar) ++ntri;
	else if(col==PDT::Colour8) ++noct;
	else if(col==PDT::Colour6\|\|col==PDT::Colour6bar) ++nsex;
	}
	if(nsng==4 \|\| (ntri==2&&nsng==2) \|\|
	(noct==3 && nsng==1) \|\|
	(ntri==2 && noct==1 && nsng==1) ) {
	vector<CFPair> cfv(1,make_pair(0,1));
	diag.colourFlow = cfv;
	}
	else if(noct==4) {
	// flows for SSVV, VVVV is handled in me class
	vector<CFPair> cfv(3);
	cfv[0] = make_pair(0, 1.);
	cfv[1] = make_pair(1,-2.);
	cfv[2] = make_pair(2, 1.);
	diag.colourFlow = cfv;
	}
	else if(ntri==2&&noct==2) {
	vector<CFPair> cfv(2);
	cfv[0] = make_pair(0, 1);
	cfv[1] = make_pair(1, 1);
	diag.colourFlow = cfv;
	}
	else if(nsex==2&&noct==2) {
	vector<CFPair> cfv;
	for(unsigned int ix=0;ix<4;++ix)
	cfv.push_back(make_pair(ix ,2.));
	for(unsigned int ix=0;ix<8;++ix)
	cfv.push_back(make_pair(4+ix,1.));
	diag.colourFlow = cfv;
	}
	else
	assert(false);
	}

	namespace {
	// Helper functor for find_if in duplicate function.
	class SameDiagramAs {
	public:
	SameDiagramAs(const HPDiagram & diag) : a(diag) {}
	bool operator()(const HPDiagram & b) const {
	return a == b;
	}
	private:
	HPDiagram a;
	};
	}

	bool HardProcessConstructor::duplicate(const HPDiagram & diag,
	const HPDVector & group) const {
	//find if a duplicate diagram exists
	HPDVector::const_iterator it =
	find_if(group.begin(), group.end(), SameDiagramAs(diag));
	return it != group.end();
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:53 PM (1 d, 12 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3797502
Default Alt Text: (39 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions