No OneTemporary
Actions

Size

37 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/General/GeneralHardME.h b/MatrixElement/General/GeneralHardME.h
	--- a/MatrixElement/General/GeneralHardME.h
	+++ b/MatrixElement/General/GeneralHardME.h
	@@ -1,517 +1,522 @@
	// -- C++ --
	//
	// GeneralHardME.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_GeneralHardME_H
	#define HERWIG_GeneralHardME_H
	//
	// This is the declaration of the GeneralHardME class.
	//

	#include "Herwig++/MatrixElement/HwMEBase.h"
	#include "ThePEG/Utilities/Exception.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig++/Models/General/HPDiagram.h"
	#include "Herwig++/MatrixElement/ProductionMatrixElement.h"
	#include "Herwig++/MatrixElement/HardVertex.h"
	#include "ThePEG/EventRecord/SpinInfo.h"
	#include "ThePEG/PDF/PolarizedBeamParticleData.h"
	#include "GeneralHardME.fh"

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

	/**
	* This defines the GeneralHardME class that is designed to serve as a
	* base class for matrix elements of specific spin structures when those
	* structures are created by a a general model, i.e. a SUSY production
	* ME. It stores a vector of diagram structures that contain the required
	* to calculate the matrix element.
	*
	* @see HwMEBase
	*/

	class GeneralHardME: public HwMEBase {

	public:

	/**
	* Convenient typedef for size_type of HPDiagram vector
	*/
	typedef vector<HPDiagram>::size_type HPCount;

	/**
	* Enum for the possible colour structures
	*/
	enum ColourStructure {UNDEFINED,
	Colour11to11,Colour11to33bar,Colour11to88,
	Colour33to33,Colour33barto11,Colour33barto33bar,
	Colour33barto66bar, Colour33barto6bar6,
	Colour33to61, Colour3bar3barto6bar1,
	Colour33to16, Colour3bar3barto16bar,
	Colour38to3bar6, Colour38to63bar,
	Colour33barto18,Colour33barto81,Colour33barto88,
	Colour38to13,Colour38to31,
	Colour38to83,Colour38to38,
	Colour3bar3barto3bar3bar,
	Colour3bar8to13bar,Colour3bar8to3bar1,
	Colour3bar8to83bar,Colour3bar8to3bar8,
	Colour88to11,Colour88to33bar,
	Colour88to66bar,Colour88to88,
	Colour88to18,Colour88to81};

	public:

	/**
	* The default constructor.
	*/
	GeneralHardME();

	public:

	/** @name Virtual functions required by the MEBase class. */
	//@{
	/**
	* Return the order in \f$\alpha_S\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInAlphaS() const;

	/**
	* Return the order in \f$\alpha_{EW}\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInAlphaEW() const;

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

	/**
	* Return the scale associated with the last set phase space point.
	*/
	virtual Energy2 scale() const {
	if(scaleChoice_==0) {
	return scaleFactor_*sHat();
	}
	- else {
	- assert( scaleChoice_== 1 );
	+ else if(scaleChoice_==1) {
	Energy2 mbar = 0.5*(meMomenta()[2].mass2()+meMomenta()[3].mass2());
	Energy2 t = 0.5*(tHat()-mbar);
	Energy2 u = 0.5*(uHat()-mbar);
	Energy2 s = 0.5*sHat();
	return scaleFactor_4.stu/(ss+tt+u*u);
	}
	+ else if(scaleChoice_ ==2) {
	+ Energy2 scale1 = meMomenta()[2].mass2()+meMomenta()[2].perp2();
	+ Energy2 scale2 = meMomenta()[3].mass2()+meMomenta()[3].perp2();
	+ return scaleFactor_sqrt(scale1scale2);
	+ }
	+ else assert(false);
	}

	/**
	* Add all possible diagrams with the add() function.
	*/
	virtual void getDiagrams() const;

	/**
	* Get diagram selector. With the information previously supplied with the
	* setKinematics method, a derived class may optionally
	* override this method to weight the given diagrams with their
	* (although certainly not physical) relative probabilities.
	* @param dv the diagrams to be weighted.
	* @return a Selector relating the given diagrams to their weights.
	*/
	virtual Selector<DiagramIndex>
	diagrams(const DiagramVector & dv) const;

	/**
	* Return a Selector with possible colour geometries for the selected
	* diagram weighted by their relative probabilities.
	* @param diag the diagram chosen.
	* @return the possible colour geometries weighted by their
	* relative probabilities.
	*/
	virtual Selector<const ColourLines *>
	colourGeometries(tcDiagPtr diag) const;
	//@}

	/**
	* Set the diagrams and matrix of colour factors.
	* @param process vector of MEDiagram with information that
	* will allow the diagrams to be created in the specific matrix element
	* @param colour The colour structure for the process
	* @param debug Whether to compare the numerical answer to an analytical
	* formula (This is only stored for certain processes. It is intended
	* for quick checks of the matrix elements).
	* @param scaleOption The option of what scale to use
	* @param scaleFactor The prefactor for the scale
	*/
	void setProcessInfo(const vector<HPDiagram> & process,
	ColourStructure colour, bool debug,
	unsigned int scaleOption,
	double scaleFactor);

	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. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();
	//@}

	protected:

	/**
	* A debugging function to test the value of me2 against an
	* analytic function. This is to be overidden in an inheriting class.
	*/
	virtual void debug(double ) const {}

	protected:

	/**
	* Access the HPDiagrams that store the required information
	* to create the diagrams
	*/
	const vector<HPDiagram> & getProcessInfo() const {
	return diagrams_;
	}

	/**
	* Return the incoming pair
	* @return Pair of particle ids for the incoming particles
	*/
	pair<long, long> getIncoming() const {
	return incoming_;
	}

	/**
	* Return the outgoing pair
	* @return Pair of particle ids for the outgoing particles
	*/
	pair<long, long> getOutgoing() const {
	return outgoing_;
	}

	/**
	* Return the matrix of colour factors
	*/
	const vector<DVector> & getColourFactors() const {
	return colour_;
	}

	/**
	* Get the number of diagrams in this process
	*/
	HPCount numberOfDiags() const {
	return numberOfDiagrams_;
	}

	/**
	* Access number of colour flows
	*/
	size_t numberOfFlows() const {
	return numberOfFlows_;
	}

	/**
	* Whether to print the debug information
	*/
	bool debugME() const {
	return debug_;
	}

	/**
	* Set/Get Info on the selected diagram and colour flow
	*/
	//@{
	/**
	* Colour flow
	*/
	unsigned int colourFlow() const {return flow_;}

	/**
	* Colour flow
	*/
	void colourFlow(unsigned int flow) const {flow_=flow;}

	/**
	* Diagram
	*/
	unsigned int diagram() const {return diagram_;}

	/**
	* Diagram
	*/
	void diagram(unsigned int diag) const {diagram_=diag;}
	//@}

	/**
	* Calculate weight and select colour flow
	*/
	double selectColourFlow(vector<double> & flow,
	vector<double> & me,double average) const;

	/**
	* Access to the colour flow matrix element
	*/
	vector<ProductionMatrixElement> & flowME() const {
	return flowME_;
	}

	/**
	* Access to the diagram matrix element
	*/
	vector<ProductionMatrixElement> & diagramME() const {
	return diagramME_;
	}

	/**
	* Access to the colour structure
	*/
	ColourStructure colour() const {return colourStructure_;}

	/**
	* Extract the paricles from the subprocess
	*/
	ParticleVector hardParticles(tSubProPtr subp) {
	ParticleVector output(4);
	output[0] = subp->incoming().first;
	output[1] = subp->incoming().second;
	output[2] = subp->outgoing()[0];
	output[3] = subp->outgoing()[1];
	//ensure particle ordering is the same as it was when
	//the diagrams were created
	if( output[0]->id() != getIncoming().first )
	swap(output[0], output[1]);
	if( output[2]->id() != getOutgoing().first )
	swap(output[2], output[3]);
	// return answer
	return output;
	}

	/**
	* Set the rescaled momenta
	*/
	void setRescaledMomenta(const ParticleVector & external) {
	cPDVector data(4);
	vector<Lorentz5Momentum> momenta(4);
	for( size_t i = 0; i < 4; ++i ) {
	data[i] = external[i]->dataPtr();
	momenta[i] = external[i]->momentum();
	}
	rescaleMomenta(momenta, data);
	}

	/**
	* Create the vertes
	*/
	void createVertex(ProductionMatrixElement & me,
	ParticleVector & external) {
	HardVertexPtr hardvertex = new_ptr(HardVertex());
	hardvertex->ME(me);
	for(ParticleVector::size_type i = 0; i < 4; ++i) {
	tSpinPtr spin = external[i]->spinInfo();
	if(i<2) {
	tcPolarizedBeamPDPtr beam =
	dynamic_ptr_cast<tcPolarizedBeamPDPtr>(external[i]->dataPtr());
	if(beam) spin->rhoMatrix() = beam->rhoMatrix();
	}
	spin->productionVertex(hardvertex);
	}
	}

	/**
	* Initialize the storage of the helicity matrix elements
	*/
	void initializeMatrixElements(PDT::Spin in1, PDT::Spin in2,
	PDT::Spin out1, PDT::Spin out2) {
	flowME().resize(numberOfFlows(),
	ProductionMatrixElement(in1,in2,out1,out2));
	diagramME().resize(numberOfDiags(),
	ProductionMatrixElement(in1,in2,out1,out2));
	}

	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is an abstract class with persistent data.
	*/
	static AbstractClassDescription<GeneralHardME> initGeneralHardME;

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

	private:

	/**
	* External particles
	*/
	//@{
	/**
	* Store incoming particles
	*/
	pair<long, long> incoming_;

	/**
	* Store the outgoing particles
	*/
	pair<long, long> outgoing_;
	//@}

	/**
	* Diagrams
	*/
	//@{
	/**
	* Store all diagrams as a vector of structures
	*/
	vector<HPDiagram> diagrams_;

	/**
	* Store the number of diagrams for fast retrieval
	*/
	HPCount numberOfDiagrams_;
	//@}

	/**
	* Colour information
	*/
	//@{
	/**
	* The colour structure
	*/
	ColourStructure colourStructure_;

	/**
	* Store colour factors for ME calc.
	*/
	vector<DVector> colour_;

	/**
	* The number of colourflows.
	*/
	unsigned int numberOfFlows_;
	//@}

	/**
	* Whether to test the value of me2 against the analytical function
	*/
	bool debug_;

	/**
	* The scale chocie
	*/
	unsigned int scaleChoice_;

	/**
	* The scale factor
	*/
	double scaleFactor_;

	/**
	* Info on the selected diagram and colour flow
	*/
	//@{
	/**
	* Colour flow
	*/
	mutable unsigned int flow_;

	/**
	* Diagram
	*/
	mutable unsigned int diagram_;
	//@}

	/**
	* Storage of the matrix elements
	*/
	//@{
	/**
	* Matrix elements for the different colour flows
	*/
	mutable vector<ProductionMatrixElement> flowME_;

	/**
	* Matrix elements for the different Feynman diagrams
	*/
	mutable vector<ProductionMatrixElement> diagramME_;
	//@}

	};

	/** Exception class to indicate a problem has occurred with setting
	up to matrix element.*/
	class MEException : public Exception {};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	/** @endcond */

	}

	#endif /* HERWIG_GeneralHardME_H */
	diff --git a/Models/General/TwoToTwoProcessConstructor.cc b/Models/General/TwoToTwoProcessConstructor.cc
	--- a/Models/General/TwoToTwoProcessConstructor.cc
	+++ b/Models/General/TwoToTwoProcessConstructor.cc
	@@ -1,645 +1,650 @@
	// -- C++ --
	//
	// TwoToTwoProcessConstructor.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 TwoToTwoProcessConstructor class.
	//

	#include "TwoToTwoProcessConstructor.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;

	namespace {
	// Helper functor for find_if in duplicate function.
	class SameIncomingAs {
	public:
	SameIncomingAs(tPDPair in) : a(in.first->id()), b(in.second->id()) {}
	bool operator()(tPDPair ppair) const {
	long id1(ppair.first->id()), id2(ppair.second->id());
	return ( id1 == a && id2 == b ) \|\| ( id1 == b && id2 == a );
	}
	private:
	long a, b;
	};

	bool duplicateIncoming(tPDPair ppair, const vector<tPDPair> & incPairs ) {
	vector<tPDPair>::const_iterator it =
	find_if( incPairs.begin(), incPairs.end(), SameIncomingAs(ppair) );
	return it != incPairs.end();
	}
	}

	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&&outgoing_.size()!=2)
	throw InitException()
	<< "Exclusive processes require exactly"
	<< " two outgoing particles but " << outgoing_.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( !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_;
	}

	ClassDescription<TwoToTwoProcessConstructor>
	TwoToTwoProcessConstructor::initTwoToTwoProcessConstructor;
	// Definition of the static class description member.

	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 interfaceIncludeAllDiagramsOff
	(interfaceIncludeAllDiagrams,
	"No",
	"Only include QCD diagrams",
	false);
	static SwitchOption interfaceIncludeAllDiagramsOn
	(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 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,
	+ "GeometicMean",
	+ "Use the geometic mean 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);

	}

	namespace {
	// Helper functor for find_if below.
	class SameProcessAs {
	public:
	SameProcessAs(const HPDiagram & diag) : a(diag) {}
	bool operator()(const HPDiagram & b) const {
	return a.sameProcess(b);
	}
	private:
	HPDiagram a;
	};
	}

	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();
	bool abort=false;
	for( ; dit != dend; ++dit) {
	// check for on-shell s-channel
	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);
	generator()->log() << 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";
	abort = true;
	}
	grouped.insert(*dit);
	}
	if(abort) throw Exception() << "One or more processes with zero width"
	<< " resonant intermediates"
	<< Exception::runerror;
	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) {
	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(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( !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( !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::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: Tue, Nov 19, 7:47 PM (1 d, 9 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3803733
Default Alt Text: (37 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions