No OneTemporary
Actions

Size

79 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
	--- a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
	+++ b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
	@@ -1,542 +1,554 @@
	// -- C++ --
	//
	// MatchboxAmplitude.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 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_MatchboxAmplitude_H
	#define HERWIG_MatchboxAmplitude_H
	//
	// This is the declaration of the MatchboxAmplitude class.
	//

	#include "ThePEG/MatrixElement/Amplitude.h"
	#include "ThePEG/Handlers/LastXCombInfo.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/ColourBasis.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxXComb.h"
	#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
	#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.fh"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief Process information with coupling order
	*/
	struct Process {

	PDVector legs;
	unsigned int orderInAlphaS;
	unsigned int orderInAlphaEW;

	Process()
	: orderInAlphaS(0), orderInAlphaEW(0) {}

	Process(const PDVector& p,
	unsigned int oas,
	unsigned int oae)
	: legs(p), orderInAlphaS(oas), orderInAlphaEW(oae) {}

	bool operator==(const Process& other) const {
	return
	legs == other.legs &&
	orderInAlphaS == other.orderInAlphaS &&
	orderInAlphaEW == other.orderInAlphaEW;
	}

	bool operator<(const Process& other) const {
	if ( orderInAlphaS != other.orderInAlphaS )
	return orderInAlphaS < other.orderInAlphaS;
	if ( orderInAlphaEW != other.orderInAlphaEW )
	return orderInAlphaEW < other.orderInAlphaEW;
	return legs < other.legs;
	}

	void persistentOutput(PersistentOStream & os) const {
	os << legs << orderInAlphaS << orderInAlphaEW;
	}

	void persistentInput(PersistentIStream & is) {
	is >> legs >> orderInAlphaS >> orderInAlphaEW;
	}

	};

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief Enumerate the type of calculation required
	*/
	namespace ProcessType {

	enum Types {

	treeME2 = 0,
	colourCorrelatedME2,
	spinColourCorrelatedME2,
	oneLoopInterference

	};

	}


	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief MatchboxAmplitude is the base class for amplitude
	* implementations inside Matchbox.
	*
	* @see \ref MatchboxAmplitudeInterfaces "The interfaces"
	* defined for MatchboxAmplitude.
	*/
	class MatchboxAmplitude:
	public Amplitude,
	public LastXCombInfo<StandardXComb>,
	public LastMatchboxXCombInfo {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	MatchboxAmplitude();

	/**
	* The destructor.
	*/
	virtual ~MatchboxAmplitude();
	//@}

	public:

	/**
	* Return the amplitude. Needs to be implemented from
	* ThePEG::Amplitude but is actually ill-defined, as colours of the
	* external particles are not specified. To this extent, this
	* implementation just asserts.
	*/
	virtual Complex value(const tcPDVector & particles,
	const vector<Lorentz5Momentum> & momenta,
	const vector<int> & helicities);

	/** @name Subprocess information */
	//@{

	/**
	* Return true, if this amplitude can handle the given process.
	*/
	virtual bool canHandle(const PDVector& p,
	Ptr<MatchboxFactory>::tptr) const { return canHandle(p); }

	/**
	* Return true, if this amplitude can handle the given process.
	*/
	virtual bool canHandle(const PDVector&) const { return false; }

	/**
	* Return the number of random numbers required to evaluate this
	* amplitude at a fixed phase space point.
	*/
	virtual int nDimAdditional() const { return 0; }

	/**
	* Return a ME instance appropriate for this amplitude and the given
	* subprocesses
	*/
	virtual Ptr<MatchboxMEBase>::ptr makeME(const PDVector&) const;

	/**
	* Set the (tree-level) order in \f$g_S\f$ in which this matrix
	* element should be evaluated.
	*/
	virtual void orderInGs(unsigned int) {}

	/**
	* Return the (tree-level) order in \f$g_S\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInGs() const = 0;

	/**
	* Set the (tree-level) order in \f$g_{EM}\f$ in which this matrix
	* element should be evaluated.
	*/
	virtual void orderInGem(unsigned int) {}

	/**
	* Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInGem() const = 0;

	/**
	* Return the Herwig++ StandardModel object
	*/
	Ptr<StandardModel>::tcptr standardModel() {
	if ( !hwStandardModel() )
	hwStandardModel(dynamic_ptr_cast<Ptr<StandardModel>::tcptr>(HandlerBase::standardModel()));
	return hwStandardModel();
	}

	/**
	* Tell whether the outgoing partons should be sorted when determining
	* allowed subprocesses. Otherwise, all permutations are counted as
	* separate subprocesses.
	*/
	virtual bool sortOutgoing() { return true; }

	/**
	+ * Return true, if this amplitude already includes averaging over
	+ * incoming parton's quantum numbers.
	+ */
	+ virtual bool hasInitialAverage() const { return false; }
	+
	+ /**
	+ * Return true, if this amplitude already includes symmetry factors
	+ * for identical outgoing particles.
	+ */
	+ virtual bool hasFinalStateSymmetry() const { return false; }
	+
	+ /**
	* Return true, if this amplitude is handled by a BLHA one-loop provider
	*/
	virtual bool isOLPTree() const { return false; }

	/**
	* Return true, if this amplitude is handled by a BLHA one-loop provider
	*/
	virtual bool isOLPLoop() const { return false; }

	/**
	* Write the order file header
	*/
	virtual void olpOrderFileHeader(ostream&) const;

	/**
	* Write the order file process list
	*/
	virtual void olpOrderFileProcessGroup(ostream&,
	const string&,
	const set<Process>&) const;

	/**
	* Write the order file process list
	*/
	virtual void olpOrderFileProcesses(ostream&,
	const map<pair<Process,int>,int>& procs) const;

	/**
	* Start the one loop provider, if appropriate, giving order and
	* contract files
	*/
	virtual void signOLP(const string&, const string&) { }

	/**
	* Start the one loop provider, if appropriate
	*/
	virtual void startOLP(const string&, int& status) { status = -1; }

	/**
	* Start the one loop provider, if appropriate. This default
	* implementation writes an BLHA 2.0 order file and starts the OLP
	*/
	virtual bool startOLP(const map<pair<Process,int>,int>& procs);


	//@}

	/** @name Colour basis. */
	//@{

	/**
	* Return the colour basis.
	*/
	Ptr<ColourBasis>::tptr colourBasis() const { return theColourBasis; }

	/**
	* Return true, if this amplitude will not require colour correlations.
	*/
	virtual bool noCorrelations() const { return !haveOneLoop(); }

	/**
	* Return true, if the colour basis is capable of assigning colour
	* flows.
	*/
	virtual bool haveColourFlows() const {
	return colourBasis() ? colourBasis()->haveColourFlows() : false;
	}

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

	/**
	* Return an ordering identifier for the current subprocess and
	* colour absis tensor index.
	*/
	const string& colourOrderingString(size_t id) const;

	/**
	* Return an ordering identifier for the current subprocess and
	* colour absis tensor index.
	*/
	const vector<vector<size_t> >& colourOrdering(size_t id) const;

	//@}

	/** @name Phasespace point, crossing and helicities */
	//@{

	/**
	* Set the xcomb object.
	*/
	virtual void setXComb(tStdXCombPtr xc);

	/**
	* Return the momentum as crossed appropriate for this amplitude.
	*/
	Lorentz5Momentum amplitudeMomentum(int) const;

	/**
	* Perform a normal ordering of external legs and fill the
	* crossing information as. This default implementation sorts
	* lexicographically in (abs(colour)/spin/abs(charge)), putting pairs
	* of particles/anti-particles where possible.
	*/
	virtual void fillCrossingMap(size_t shift = 0);

	/**
	* Generate the helicity combinations.
	*/
	virtual set<vector<int> > generateHelicities() const;

	//@}

	/** @name Tree-level amplitudes */
	//@{

	/**
	* Calculate the tree level amplitudes for the phasespace point
	* stored in lastXComb.
	*/
	virtual void prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr);

	/**
	* Return the matrix element squared.
	*/
	virtual double me2() const;

	/**
	* Return the colour correlated matrix element.
	*/
	virtual double colourCorrelatedME2(pair<int,int> ij) const;

	/**
	* Return the large-N colour correlated matrix element.
	*/
	virtual double largeNColourCorrelatedME2(pair<int,int> ij,
	Ptr<ColourBasis>::tptr largeNBasis) const;

	/**
	* Return a positive helicity polarization vector for a gluon of
	* momentum p (with reference vector n) to be used when evaluating
	* spin correlations.
	*/
	virtual LorentzVector<Complex> plusPolarization(const Lorentz5Momentum& p,
	const Lorentz5Momentum& n) const;

	/**
	* Return the colour and spin correlated matrix element.
	*/
	virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
	const SpinCorrelationTensor& c) const;


	/**
	* Return true, if tree-level contributions will be evaluated at amplitude level.
	*/
	virtual bool treeAmplitudes() const { return true; }

	/**
	* Evaluate the amplitude for the given colour tensor id and
	* helicity assignment
	*/
	virtual Complex evaluate(size_t, const vector<int>&, Complex&) { return 0.; }

	//@}

	/** @name One-loop amplitudes */
	//@{

	/**
	* Return true, if this amplitude is capable of calculating one-loop
	* (QCD) corrections.
	*/
	virtual bool haveOneLoop() const { return false; }

	/**
	* Return true, if this amplitude only provides
	* one-loop (QCD) corrections.
	*/
	virtual bool onlyOneLoop() const { return false; }

	/**
	* Return true, if one-loop contributions will be evaluated at amplitude level.
	*/
	virtual bool oneLoopAmplitudes() const { return true; }

	/**
	* Return true, if one loop corrections have been calculated in
	* dimensional reduction. Otherwise conventional dimensional
	* regularization is assumed. Note that renormalization is always
	* assumed to be MSbar.
	*/
	virtual bool isDR() const { return false; }

	/**
	* Return true, if one loop corrections are given in the conventions
	* of the integrated dipoles.
	*/
	virtual bool isCS() const { return false; }

	/**
	* Return true, if one loop corrections are given in the conventions
	* of BDK.
	*/
	virtual bool isBDK() const { return false; }

	/**
	* Return true, if one loop corrections are given in the conventions
	* of everything expanded.
	*/
	virtual bool isExpanded() const { return false; }

	/**
	* Return the value of the dimensional regularization
	* parameter. Note that renormalization scale dependence is fully
	* restored in DipoleIOperator.
	*/
	virtual Energy2 mu2() const { return 0.*GeV2; }

	/**
	* If defined, return the coefficient of the pole in epsilon^2
	*/
	virtual double oneLoopDoublePole() const { return 0.; }

	/**
	* If defined, return the coefficient of the pole in epsilon
	*/
	virtual double oneLoopSinglePole() const { return 0.; }

	/**
	* Calculate the one-loop amplitudes for the phasespace point
	* stored in lastXComb, if provided.
	*/
	virtual void prepareOneLoopAmplitudes(Ptr<MatchboxMEBase>::tcptr);

	/**
	* Return the one-loop/tree interference.
	*/
	virtual double oneLoopInterference() const;

	/**
	* Evaluate the amplitude for the given colour tensor id and
	* helicity assignment
	*/
	virtual Complex evaluateOneLoop(size_t, const vector<int>&) { return 0.; }

	//@}

	/** @name Caching and helpers to setup amplitude objects. */
	//@{

	/**
	* Flush all cashes.
	*/
	virtual void flushCaches() {}

	/**
	* Clone this amplitude.
	*/
	Ptr<MatchboxAmplitude>::ptr cloneMe() const {
	return dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(clone());
	}

	/**
	* Clone the dependencies, using a given prefix.
	*/
	virtual void cloneDependencies(const std::string& prefix = "");

	//@}

	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();

	// If needed, insert declarations of virtual function defined in the
	// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).

	private:

	/**
	* Recursively generate helicities
	*/
	void doGenerateHelicities(set<vector<int> >& res,
	vector<int>& current,
	size_t pos) const;

	/**
	* The colour basis implementation to be used.
	*/
	Ptr<ColourBasis>::ptr theColourBasis;

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

	};

	inline PersistentOStream& operator<<(PersistentOStream& os,
	const Process& h) {
	h.persistentOutput(os);
	return os;
	}

	inline PersistentIStream& operator>>(PersistentIStream& is,
	Process& h) {
	h.persistentInput(is);
	return is;
	}

	}

	#endif /* HERWIG_MatchboxAmplitude_H */
	diff --git a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
	--- a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
	+++ b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
	@@ -1,1305 +1,1310 @@
	// -- C++ --
	//
	// MatchboxMEBase.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 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 MatchboxMEBase class.
	//

	#include "MatchboxMEBase.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDF/PDF.h"
	#include "ThePEG/PDT/PDT.h"
	#include "ThePEG/StandardModel/StandardModelBase.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/Handlers/StdXCombGroup.h"
	#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
	#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"

	#include <iterator>
	using std::ostream_iterator;

	using namespace Herwig;

	MatchboxMEBase::MatchboxMEBase()
	: MEBase(),
	theOneLoop(false),
	theOneLoopNoBorn(false) {}

	MatchboxMEBase::~MatchboxMEBase() {}

	Ptr<MatchboxFactory>::tcptr MatchboxMEBase::factory() const { return theFactory; }

	void MatchboxMEBase::factory(Ptr<MatchboxFactory>::tcptr f) { theFactory = f; }

	Ptr<Tree2toNGenerator>::tptr MatchboxMEBase::diagramGenerator() const { return factory()->diagramGenerator(); }

	Ptr<ProcessData>::tptr MatchboxMEBase::processData() const { return factory()->processData(); }

	unsigned int MatchboxMEBase::getNLight() const { return factory()->nLight(); }

	double MatchboxMEBase::factorizationScaleFactor() const { return factory()->factorizationScaleFactor(); }

	double MatchboxMEBase::renormalizationScaleFactor() const { return factory()->renormalizationScaleFactor(); }

	bool MatchboxMEBase::fixedCouplings() const { return factory()->fixedCouplings(); }

	bool MatchboxMEBase::fixedQEDCouplings() const { return factory()->fixedQEDCouplings(); }

	bool MatchboxMEBase::checkPoles() const { return factory()->checkPoles(); }

	bool MatchboxMEBase::verbose() const { return factory()->verbose(); }

	bool MatchboxMEBase::initVerbose() const { return factory()->initVerbose(); }

	void MatchboxMEBase::getDiagrams() const {

	if ( diagramGenerator() && processData() ) {

	vector<Ptr<Tree2toNDiagram>::ptr> diags;

	vector<Ptr<Tree2toNDiagram>::ptr>& res =
	processData()->diagramMap()[subProcess().legs];
	if ( res.empty() ) {
	res = diagramGenerator()->generate(subProcess().legs,orderInAlphaS(),orderInAlphaEW());
	}
	copy(res.begin(),res.end(),back_inserter(diags));
	processData()->fillMassGenerators(subProcess().legs);

	if ( diags.empty() )
	return;

	for ( vector<Ptr<Tree2toNDiagram>::ptr>::iterator d = diags.begin();
	d != diags.end(); ++d ) {
	add(*d);
	}

	return;

	}

	throw Exception()
	<< "MatchboxMEBase::getDiagrams() expects a Tree2toNGenerator and ProcessData object.\n"
	<< "Please check your setup." << Exception::abortnow;

	}

	Selector<MEBase::DiagramIndex>
	MatchboxMEBase::diagrams(const DiagramVector & diags) const {

	if ( phasespace() ) {
	return phasespace()->selectDiagrams(diags);
	}

	throw Exception()
	<< "MatchboxMEBase::diagrams() expects a MatchboxPhasespace object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return Selector<MEBase::DiagramIndex>();

	}

	Selector<const ColourLines *>
	MatchboxMEBase::colourGeometries(tcDiagPtr diag) const {

	if ( matchboxAmplitude() ) {
	if ( matchboxAmplitude()->haveColourFlows() ) {
	if ( matchboxAmplitude()->treeAmplitudes() )
	matchboxAmplitude()->prepareAmplitudes(this);
	return matchboxAmplitude()->colourGeometries(diag);
	}
	}

	Ptr<Tree2toNDiagram>::tcptr tdiag =
	dynamic_ptr_cast<Ptr<Tree2toNDiagram>::tcptr>(diag);
	assert(diag && processData());

	vector<ColourLines*>& flows = processData()->colourFlowMap()[tdiag];

	if ( flows.empty() ) {

	list<list<list<pair<int,bool> > > > cflows =
	ColourBasis::colourFlows(tdiag);

	for ( list<list<list<pair<int,bool> > > >::const_iterator fit =
	cflows.begin(); fit != cflows.end(); ++fit ) {
	flows.push_back(new ColourLines(ColourBasis::cfstring(*fit)));
	}

	}

	Selector<const ColourLines *> res;
	for ( vector<ColourLines*>::const_iterator f = flows.begin();
	f != flows.end(); ++f )
	res.insert(1.0,*f);

	return res;

	}

	unsigned int MatchboxMEBase::orderInAlphaS() const {
	return subProcess().orderInAlphaS;
	}

	unsigned int MatchboxMEBase::orderInAlphaEW() const {
	return subProcess().orderInAlphaEW;
	}

	void MatchboxMEBase::setXComb(tStdXCombPtr xc) {

	MEBase::setXComb(xc);
	lastMatchboxXComb(xc);
	if ( phasespace() )
	phasespace()->setXComb(xc);
	if ( scaleChoice() )
	scaleChoice()->setXComb(xc);
	if ( matchboxAmplitude() )
	matchboxAmplitude()->setXComb(xc);

	}

	double MatchboxMEBase::generateIncomingPartons(const double* r1, const double* r2) {

	// shamelessly stolen from PartonExtractor.cc

	Energy2 shmax = lastCuts().sHatMax();
	Energy2 shmin = lastCuts().sHatMin();
	Energy2 sh = shminpow(shmax/shmin, r1);
	double ymax = lastCuts().yHatMax();
	double ymin = lastCuts().yHatMin();
	double km = log(shmax/shmin);
	ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/sh)));
	ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/sh)));

	double y = ymin + (r2)(ymax - ymin);
	double x1 = exp(-0.5*log(lastS()/sh) + y);
	double x2 = exp(-0.5*log(lastS()/sh) - y);

	Lorentz5Momentum P1 = lastParticles().first->momentum();
	LorentzMomentum p1 = lightCone((P1.rho() + P1.e())*x1, Energy());
	p1.rotateY(P1.theta());
	p1.rotateZ(P1.phi());
	meMomenta()[0] = p1;

	Lorentz5Momentum P2 = lastParticles().second->momentum();
	LorentzMomentum p2 = lightCone((P2.rho() + P2.e())*x2, Energy());
	p2.rotateY(P2.theta());
	p2.rotateZ(P2.phi());
	meMomenta()[1] = p2;

	lastXCombPtr()->lastX1X2(make_pair(x1,x2));
	lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());

	return km*(ymax - ymin);

	}

	bool MatchboxMEBase::generateKinematics(const double * r) {

	if ( phasespace() ) {

	jacobian(phasespace()->generateKinematics(r,meMomenta()));
	if ( jacobian() == 0.0 )
	return false;

	setScale();
	logGenerateKinematics(r);

	assert(lastMatchboxXComb());

	if ( nDimAmplitude() > 0 ) {
	amplitudeRandomNumbers().resize(nDimAmplitude());
	copy(r + nDimPhasespace(),
	r + nDimPhasespace() + nDimAmplitude(),
	amplitudeRandomNumbers().begin());
	}

	if ( nDimInsertions() > 0 ) {
	insertionRandomNumbers().resize(nDimInsertions());
	copy(r + nDimPhasespace() + nDimAmplitude(),
	r + nDimPhasespace() + nDimAmplitude() + nDimInsertions(),
	insertionRandomNumbers().begin());
	}

	return true;

	}

	throw Exception()
	<< "MatchboxMEBase::generateKinematics() expects a MatchboxPhasespace object.\n"
	<< "Please check your setup." << Exception::abortnow;

	return false;

	}

	int MatchboxMEBase::nDim() const {

	if ( lastMatchboxXComb() )
	return nDimPhasespace() + nDimAmplitude() + nDimInsertions();

	int ampAdd = 0;
	if ( matchboxAmplitude() ) {
	ampAdd = matchboxAmplitude()->nDimAdditional();
	}

	int insertionAdd = 0;
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	insertionAdd = max(insertionAdd,(**v).nDimAdditional());
	}

	return nDimBorn() + ampAdd + insertionAdd;

	}

	int MatchboxMEBase::nDimBorn() const {

	if ( lastMatchboxXComb() )
	return nDimPhasespace();

	if ( phasespace() ) {
	size_t nout = diagrams().front()->partons().size()-2;
	int n = phasespace()->nDim(nout);
	if ( phasespace()->useMassGenerators() ) {
	for ( cPDVector::const_iterator pd =
	diagrams().front()->partons().begin();
	pd != diagrams().front()->partons().end(); ++pd ) {
	if ( processData()->massGenerator(*pd) \|\|
	(**pd).width() != ZERO ) {
	++n;
	}
	}
	}
	return n;
	}

	throw Exception()
	<< "MatchboxMEBase::nDim() expects a MatchboxPhasespace object.\n"
	<< "Please check your setup." << Exception::abortnow;

	return 0;

	}

	void MatchboxMEBase::setScale() const {
	if ( haveX1X2() ) {
	lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
	}
	Energy2 fscale = factorizationScale()*sqr(factorizationScaleFactor());
	Energy2 rscale = renormalizationScale()*sqr(renormalizationScaleFactor());
	Energy2 ewrscale = renormalizationScaleQED();
	lastXCombPtr()->lastScale(fscale);
	if ( !fixedCouplings() ) {
	if ( rscale > lastCuts().scaleMin() )
	lastXCombPtr()->lastAlphaS(SM().alphaS(rscale));
	else
	lastXCombPtr()->lastAlphaS(SM().alphaS(lastCuts().scaleMin()));
	} else {
	lastXCombPtr()->lastAlphaS(SM().alphaS());
	}
	if ( !fixedQEDCouplings() ) {
	lastXCombPtr()->lastAlphaEM(SM().alphaEM(ewrscale));
	} else {
	lastXCombPtr()->lastAlphaEM(SM().alphaEMMZ());
	}
	logSetScale();
	}

	Energy2 MatchboxMEBase::factorizationScale() const {
	if ( scaleChoice() ) {
	return scaleChoice()->factorizationScale();
	}

	throw Exception()
	<< "MatchboxMEBase::factorizationScale() expects a MatchboxScaleChoice object.\n"
	<< "Please check your setup." << Exception::abortnow;

	return ZERO;

	}

	Energy2 MatchboxMEBase::renormalizationScale() const {
	if ( scaleChoice() ) {
	return scaleChoice()->renormalizationScale();
	}

	throw Exception()
	<< "MatchboxMEBase::renormalizationScale() expects a MatchboxScaleChoice object.\n"
	<< "Please check your setup." << Exception::abortnow;

	return ZERO;

	}

	Energy2 MatchboxMEBase::renormalizationScaleQED() const {
	if ( scaleChoice() ) {
	return scaleChoice()->renormalizationScaleQED();
	}
	return renormalizationScale();
	}

	void MatchboxMEBase::setVetoScales(tSubProPtr) const {}

	void MatchboxMEBase::getPDFWeight(Energy2 factorizationScale) const {

	if ( !mePartonData()[0]->coloured() &&
	!mePartonData()[1]->coloured() ) {
	lastMEPDFWeight(1.0);
	logPDFWeight();
	return;
	}

	double w = 1.;

	if ( mePartonData()[0]->coloured() && havePDFWeight1() )
	w *= pdf1(factorizationScale);

	if ( mePartonData()[1]->coloured() && havePDFWeight2() )
	w *= pdf2(factorizationScale);

	lastMEPDFWeight(w);

	logPDFWeight();

	}

	double MatchboxMEBase::pdf1(Energy2 fscale, double xEx) const {

	assert(lastXCombPtr()->partonBins().first->pdf());

	if ( xEx < 1. && lastX1() >= xEx ) {
	return
	( ( 1. - lastX1() ) / ( 1. - xEx ) ) *
	lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
	lastPartons().first->dataPtr(),
	fscale == ZERO ? lastScale() : fscale,
	xEx)/xEx;
	}

	return lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
	lastPartons().first->dataPtr(),
	fscale == ZERO ? lastScale() : fscale,
	lastX1())/lastX1();
	}

	double MatchboxMEBase::pdf2(Energy2 fscale, double xEx) const {

	assert(lastXCombPtr()->partonBins().second->pdf());

	if ( xEx < 1. && lastX2() >= xEx ) {
	return
	( ( 1. - lastX2() ) / ( 1. - xEx ) ) *
	lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
	lastPartons().second->dataPtr(),
	fscale == ZERO ? lastScale() : fscale,
	xEx)/xEx;
	}

	return lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
	lastPartons().second->dataPtr(),
	fscale == ZERO ? lastScale() : fscale,
	lastX2())/lastX2();

	}

	double MatchboxMEBase::me2() const {

	if ( matchboxAmplitude() ) {

	if ( matchboxAmplitude()->treeAmplitudes() )
	matchboxAmplitude()->prepareAmplitudes(this);

	lastME2(matchboxAmplitude()->me2()*
	crossingSign()*
	me2Norm());

	logME2();

	return lastME2();

	}

	throw Exception()
	<< "MatchboxMEBase::me2() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0.;

	}

	double MatchboxMEBase::finalStateSymmetry() const {

	if ( symmetryFactor() > 0.0 )
	return symmetryFactor();

	double sFactor = 1.;

	map<long,int> counts;

	cPDVector checkData;
	copy(mePartonData().begin()+2,mePartonData().end(),back_inserter(checkData));

	cPDVector::iterator p = checkData.begin();
	while ( !checkData.empty() ) {
	if ( counts.find((**p).id()) != counts.end() ) {
	counts[(**p).id()] += 1;
	} else {
	counts[(**p).id()] = 1;
	}
	checkData.erase(p);
	p = checkData.begin();
	continue;
	}

	for ( map<long,int>::const_iterator c = counts.begin();
	c != counts.end(); ++c ) {
	if ( c->second == 1 )
	continue;
	if ( c->second == 2 )
	sFactor /= 2.;
	else if ( c->second == 3 )
	sFactor /= 6.;
	else if ( c->second == 4 )
	sFactor /= 24.;
	}

	symmetryFactor(sFactor);

	return symmetryFactor();

	}

	double MatchboxMEBase::me2Norm(unsigned int addAlphaS) const {

	// assume that we always have incoming
	// spin-1/2 or massless spin-1 particles
	double fac = 1./4.;

	+ if ( hasInitialAverage() )
	+ fac = 1.;
	+
	if ( orderInAlphaS() > 0 \|\| addAlphaS != 0 )
	fac *= pow(lastAlphaS()/SM().alphaS(),double(orderInAlphaS()+addAlphaS));
	if ( orderInAlphaEW() > 0 )
	fac *= pow(lastAlphaEM()/SM().alphaEM(),double(orderInAlphaEW()));

	- if ( mePartonData()[0]->iColour() == PDT::Colour3 \|\|
	- mePartonData()[0]->iColour() == PDT::Colour3bar )
	- fac /= SM().Nc();
	- else if ( mePartonData()[0]->iColour() == PDT::Colour8 )
	- fac /= (SM().Nc()*SM().Nc()-1.);
	+ if ( !hasInitialAverage() ) {
	+ if ( mePartonData()[0]->iColour() == PDT::Colour3 \|\|
	+ mePartonData()[0]->iColour() == PDT::Colour3bar )
	+ fac /= SM().Nc();
	+ else if ( mePartonData()[0]->iColour() == PDT::Colour8 )
	+ fac /= (SM().Nc()*SM().Nc()-1.);

	- if ( mePartonData()[1]->iColour() == PDT::Colour3 \|\|
	- mePartonData()[1]->iColour() == PDT::Colour3bar )
	- fac /= SM().Nc();
	- else if ( mePartonData()[1]->iColour() == PDT::Colour8 )
	- fac /= (SM().Nc()*SM().Nc()-1.);
	+ if ( mePartonData()[1]->iColour() == PDT::Colour3 \|\|
	+ mePartonData()[1]->iColour() == PDT::Colour3bar )
	+ fac /= SM().Nc();
	+ else if ( mePartonData()[1]->iColour() == PDT::Colour8 )
	+ fac /= (SM().Nc()*SM().Nc()-1.);
	+ }

	- return finalStateSymmetry()*fac;
	+ return !hasFinalStateSymmetry() ? finalStateSymmetry()*fac : fac;

	}

	CrossSection MatchboxMEBase::dSigHatDR() const {

	getPDFWeight();

	if ( !lastXCombPtr()->willPassCuts() ) {
	lastME2(0.0);
	lastMECrossSection(ZERO);
	return lastMECrossSection();
	}

	double xme2 = me2();
	lastME2(xme2);

	if ( xme2 == 0. && !oneLoopNoBorn() ) {
	lastME2(0.0);
	lastMECrossSection(ZERO);
	return lastMECrossSection();
	}

	double vme2 = 0.;
	if ( oneLoop() )
	vme2 = oneLoopInterference();

	CrossSection res = ZERO;

	if ( !oneLoopNoBorn() )
	res +=
	(sqr(hbarc)/(2.lastSHat()))
	jacobian()* lastMEPDFWeight() * xme2;

	if ( oneLoop() )
	res +=
	(sqr(hbarc)/(2.lastSHat()))
	jacobian()* lastMEPDFWeight() * vme2;

	if ( !onlyOneLoop() ) {
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	(**v).setXComb(lastXCombPtr());
	res += (**v).dSigHatDR();
	}
	if ( checkPoles() )
	logPoles();
	}

	double weight = 0.0;
	bool applied = false;
	for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
	theReweights.begin(); rw != theReweights.end(); ++rw ) {
	(**rw).setXComb(lastXCombPtr());
	if ( !(**rw).apply() )
	continue;
	weight += (**rw).evaluate();
	applied = true;
	}
	if ( applied )
	res *= weight;
	lastMECrossSection(res);

	return lastMECrossSection();

	}

	double MatchboxMEBase::oneLoopInterference() const {

	if ( matchboxAmplitude() ) {

	if ( matchboxAmplitude()->oneLoopAmplitudes() )
	matchboxAmplitude()->prepareOneLoopAmplitudes(this);
	lastME2(matchboxAmplitude()->oneLoopInterference()*
	crossingSign()*
	me2Norm(1));

	logME2();

	return lastME2();

	}

	throw Exception()
	<< "MatchboxMEBase::oneLoopInterference() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0.;

	}

	MatchboxMEBase::AccuracyHistogram::AccuracyHistogram(double low,
	double up,
	unsigned int nbins)
	: lower(low), upper(up),
	sameSign(0), oppositeSign(0), nans(0),
	overflow(0), underflow(0) {

	double step = (up-low)/nbins;

	for ( unsigned int k = 1; k <= nbins; ++k )
	bins[lower + k*step] = 0.0;

	}

	void MatchboxMEBase::AccuracyHistogram::book(double a, double b) {
	if ( isnan(a) \|\| isnan(b) \|\|
	isinf(a) \|\| isinf(b) ) {
	++nans;
	return;
	}
	if ( a*b >= 0. )
	++sameSign;
	if ( a*b < 0. )
	++oppositeSign;
	double r = 1.;
	if ( abs(a) != 0.0 )
	r = abs(1.-abs(b/a));
	else if ( abs(b) != 0.0 )
	r = abs(b);
	if ( log(r) < lower \|\| r == 0.0 ) {
	++underflow;
	return;
	}
	if ( log(r) > upper ) {
	++overflow;
	return;
	}
	map<double,double>::iterator bin =
	bins.upper_bound(log(r));
	if ( bin == bins.end() )
	return;
	bin->second += 1.;
	}

	void MatchboxMEBase::AccuracyHistogram::dump(const std::string& prefix,
	const cPDVector& proc) const {
	ostringstream fname("");
	for ( cPDVector::const_iterator p = proc.begin();
	p != proc.end(); ++p )
	fname << (**p).PDGName();
	ofstream out((prefix+fname.str()+".dat").c_str());
	out << "# same sign : " << sameSign << " opposite sign : "
	<< oppositeSign << " nans : " << nans
	<< " overflow : " << overflow
	<< " underflow : " << underflow << "\n";
	for ( map<double,double>::const_iterator b = bins.begin();
	b != bins.end(); ++b ) {
	map<double,double>::const_iterator bp = b; --bp;
	if ( b->second != 0. ) {
	if ( b != bins.begin() )
	out << bp->first;
	else
	out << lower;
	out << " " << b->first
	<< " " << b->second
	<< "\n" << flush;
	}
	}
	}

	void MatchboxMEBase::AccuracyHistogram::persistentOutput(PersistentOStream& os) const {
	os << lower << upper << bins
	<< sameSign << oppositeSign << nans
	<< overflow << underflow;
	}

	void MatchboxMEBase::AccuracyHistogram::persistentInput(PersistentIStream& is) {
	is >> lower >> upper >> bins
	>> sameSign >> oppositeSign >> nans
	>> overflow >> underflow;
	}

	void MatchboxMEBase::logPoles() const {
	double res2me = oneLoopDoublePole();
	double res1me = oneLoopSinglePole();
	double res2i = 0.;
	double res1i = 0.;
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	res2i += (**v).oneLoopDoublePole();
	res1i += (**v).oneLoopSinglePole();
	}
	epsilonSquarePoleHistograms[mePartonData()].book(res2me,res2i);
	epsilonPoleHistograms[mePartonData()].book(res1me,res1i);
	}

	bool MatchboxMEBase::haveOneLoop() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->haveOneLoop();
	return false;
	}

	bool MatchboxMEBase::onlyOneLoop() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->onlyOneLoop();
	return false;
	}

	bool MatchboxMEBase::isDR() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->isDR();
	return false;
	}

	bool MatchboxMEBase::isCS() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->isCS();
	return false;
	}

	bool MatchboxMEBase::isBDK() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->isBDK();
	return false;
	}

	bool MatchboxMEBase::isExpanded() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->isExpanded();
	return false;
	}

	Energy2 MatchboxMEBase::mu2() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->mu2();
	return 0*GeV2;
	}

	double MatchboxMEBase::oneLoopDoublePole() const {

	if ( matchboxAmplitude() ) {

	return
	matchboxAmplitude()->oneLoopDoublePole()*
	crossingSign()*
	me2Norm(1);

	}

	return 0.;

	}

	double MatchboxMEBase::oneLoopSinglePole() const {

	if ( matchboxAmplitude() ) {

	return
	matchboxAmplitude()->oneLoopSinglePole()*
	crossingSign()*
	me2Norm(1);

	}

	return 0.;

	}

	vector<Ptr<SubtractionDipole>::ptr>
	MatchboxMEBase::getDipoles(const vector<Ptr<SubtractionDipole>::ptr>& dipoles,
	const vector<Ptr<MatchboxMEBase>::ptr>& borns) const {

	vector<Ptr<SubtractionDipole>::ptr> res;

	// keep track of the dipoles we already did set up
	set<pair<pair<pair<int,int>,int>,pair<Ptr<MatchboxMEBase>::tptr,Ptr<SubtractionDipole>::tptr> > > done;

	cPDVector rep = diagrams().front()->partons();
	int nreal = rep.size();

	// now loop over configs
	for ( int emitter = 0; emitter < nreal; ++emitter ) {
	for ( int spectator = 0; spectator < nreal; ++spectator ) {
	if ( emitter == spectator )
	continue;
	for ( int emission = 2; emission < nreal; ++emission ) {
	if ( emission == emitter \|\| emission == spectator )
	continue;
	for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator b =
	borns.begin(); b != borns.end(); ++b ) {
	if ( (**b).onlyOneLoop() )
	continue;
	for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
	dipoles.begin(); d != dipoles.end(); ++d ) {
	if ( !rep[emitter]->coloured() \|\|
	!rep[emission]->coloured() \|\|
	!rep[spectator]->coloured() ) {
	continue;
	}
	if ( noDipole(emitter,emission,spectator) ) {
	continue;
	}
	if ( done.find(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(b,d)))
	!= done.end() ) {
	continue;
	}
	if ( !(**d).canHandle(rep,emitter,emission,spectator) ) {
	continue;
	}
	// now get to work
	(**d).clearBookkeeping();
	(**d).realEmitter(emitter);
	(**d).realEmission(emission);
	(**d).realSpectator(spectator);
	(*d).realEmissionME(const_cast<MatchboxMEBase>(this));
	(*d).underlyingBornME(b);
	(**d).setupBookkeeping();
	if ( !((**d).empty()) ) {
	Ptr<SubtractionDipole>::ptr nDipole = (**d).cloneMe();
	res.push_back(nDipole);
	done.insert(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(b,d)));
	if ( nDipole->isSymmetric() )
	done.insert(make_pair(make_pair(make_pair(emission,emitter),spectator),make_pair(b,d)));
	ostringstream dname;
	dname << fullName() << "." << (**b).name() << "."
	<< (**d).name() << ".[("
	<< emitter << "," << emission << ")," << spectator << "]";
	if ( ! (generator()->preinitRegister(nDipole,dname.str()) ) )
	throw InitException() << "Dipole " << dname.str() << " already existing.";
	nDipole->cloneDependencies(dname.str());
	}
	}
	}
	}
	}
	}

	for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d = res.begin();
	d != res.end(); ++d )
	(**d).partnerDipoles(res);

	return res;

	}

	double MatchboxMEBase::colourCorrelatedME2(pair<int,int> ij) const {

	if ( matchboxAmplitude() ) {

	if ( matchboxAmplitude()->treeAmplitudes() )
	matchboxAmplitude()->prepareAmplitudes(this);
	lastME2(matchboxAmplitude()->colourCorrelatedME2(ij)*
	crossingSign()*
	me2Norm());

	logME2();

	return lastME2();

	}

	throw Exception()
	<< "MatchboxMEBase::colourCorrelatedME2() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0.;

	}

	double MatchboxMEBase::largeNColourCorrelatedME2(pair<int,int> ij,
	Ptr<ColourBasis>::tptr largeNBasis) const {

	if ( matchboxAmplitude() ) {

	if ( matchboxAmplitude()->treeAmplitudes() )
	matchboxAmplitude()->prepareAmplitudes(this);
	largeNBasis->prepare(mePartonData(),false);
	lastME2(matchboxAmplitude()->largeNColourCorrelatedME2(ij,largeNBasis)*
	crossingSign()*
	me2Norm());

	logME2();

	return lastME2();

	}

	throw Exception()
	<< "MatchboxMEBase::largeNColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0.;

	}

	double MatchboxMEBase::spinColourCorrelatedME2(pair<int,int> ij,
	const SpinCorrelationTensor& c) const {

	if ( matchboxAmplitude() ) {

	if ( matchboxAmplitude()->treeAmplitudes() )
	matchboxAmplitude()->prepareAmplitudes(this);
	lastME2(matchboxAmplitude()->spinColourCorrelatedME2(ij,c)*
	crossingSign()*
	me2Norm());

	logME2();

	return lastME2();

	}

	throw Exception()
	<< "MatchboxMEBase::spinColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0.;

	}

	void MatchboxMEBase::flushCaches() {
	MEBase::flushCaches();
	if ( matchboxAmplitude() )
	matchboxAmplitude()->flushCaches();
	for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
	reweights().begin(); r != reweights().end(); ++r ) {
	(**r).flushCaches();
	}
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	(**v).flushCaches();
	}
	}

	void MatchboxMEBase::print(ostream& os) const {

	os << "--- MatchboxMEBase setup -------------------------------------------------------\n";

	os << " '" << name() << "' for subprocess:\n";

	os << " ";
	for ( PDVector::const_iterator pp = subProcess().legs.begin();
	pp != subProcess().legs.end(); ++pp ) {
	os << (**pp).PDGName() << " ";
	if ( pp == subProcess().legs.begin() + 1 )
	os << "-> ";
	}
	os << "\n";

	os << " including " << (oneLoop() ? "" : "no ") << "virtual corrections";
	if ( oneLoopNoBorn() )
	os << " without Born contributions";
	os << "\n";

	if ( oneLoop() && !onlyOneLoop() ) {
	os << " using insertion operators\n";
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	os << " '" << (**v).name() << "' with "
	<< ((**v).isDR() ? "" : "C") << "DR/";
	if ( (**v).isCS() )
	os << "CS";
	if ( (**v).isBDK() )
	os << "BDK";
	if ( (**v).isExpanded() )
	os << "expanded";
	os << " conventions\n";
	}
	}

	os << "--------------------------------------------------------------------------------\n";

	os << flush;

	}

	void MatchboxMEBase::printLastEvent(ostream& os) const {

	os << "--- MatchboxMEBase last event information --------------------------------------\n";

	os << " for matrix element '" << name() << "'\n";

	os << " process considered:\n ";

	int in = 0;
	for ( cPDVector::const_iterator p = mePartonData().begin();
	p != mePartonData().end(); ++p ) {
	os << (**p).PDGName() << " ";
	if ( ++in == 2 )
	os << " -> ";
	}

	os << " kinematic environment as set by the XComb " << lastXCombPtr() << ":\n"
	<< " sqrt(shat)/GeV = " << sqrt(lastSHat()/GeV2)
	<< " x1 = " << lastX1() << " x2 = " << lastX2()
	<< " alphaS = " << lastAlphaS() << "\n";

	os << " momenta/GeV generated from random numbers\n ";
	copy(lastXComb().lastRandomNumbers().begin(),
	lastXComb().lastRandomNumbers().end(),ostream_iterator<double>(os," "));
	os << ":\n ";

	for ( vector<Lorentz5Momentum>::const_iterator p = meMomenta().begin();
	p != meMomenta().end(); ++p ) {
	os << (*p/GeV) << "\n ";
	}

	os << "last cross section/nb calculated was:\n "
	<< (lastMECrossSection()/nanobarn) << " (pdf weight " << lastMEPDFWeight() << ")\n";

	os << "--------------------------------------------------------------------------------\n";

	os << flush;

	}

	void MatchboxMEBase::logGenerateKinematics(const double * r) const {

	if ( !verbose() )
	return;

	generator()->log() << "'" << name() << "' generated kinematics\nfrom "
	<< nDim() << " random numbers:\n";
	copy(r,r+nDim(),ostream_iterator<double>(generator()->log()," "));
	generator()->log() << "\n";

	generator()->log() << "storing phase space information in XComb "
	<< lastXCombPtr() << "\n";

	generator()->log() << "generated phase space point (in GeV):\n";

	vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
	cPDVector::const_iterator dit = mePartonData().begin();

	for ( ; pit != meMomenta().end() ; ++pit, ++dit )
	generator()->log() << (**dit).PDGName() << " : "
	<< (*pit/GeV) << "\n";

	generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
	<< "and Jacobian = " << jacobian() << " sHat/GeV2 = "
	<< (lastSHat()/GeV2) << "\n" << flush;

	}

	void MatchboxMEBase::logSetScale() const {

	if ( !verbose() )
	return;

	generator()->log() << "'" << name() << "' set scales using XComb " << lastXCombPtr() << ":\n"
	<< "scale/GeV2 = " << (scale()/GeV2) << " xi_R = "
	<< renormalizationScaleFactor() << " xi_F = "
	<< factorizationScaleFactor() << "\n"
	<< "alpha_s = " << lastAlphaS() << "\n" << flush;

	}

	void MatchboxMEBase::logPDFWeight() const {

	if ( !verbose() )
	return;

	generator()->log() << "'" << name() << "' calculated pdf weight = "
	<< lastMEPDFWeight() << " from XComb "
	<< lastXCombPtr() << "\n"
	<< "x1 = " << lastX1() << " (" << (mePartonData()[0]->coloured() ? "" : "not ") << "used) "
	<< "x2 = " << lastX2() << " (" << (mePartonData()[1]->coloured() ? "" : "not ") << "used)\n"
	<< flush;

	}

	void MatchboxMEBase::logME2() const {

	if ( !verbose() )
	return;

	generator()->log() << "'" << name() << "' evaluated me2 using XComb "
	<< lastXCombPtr() << "\n"
	<< "and phase space point (in GeV):\n";

	vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
	cPDVector::const_iterator dit = mePartonData().begin();

	for ( ; pit != meMomenta().end() ; ++pit, ++dit )
	generator()->log() << (**dit).PDGName() << " : "
	<< (*pit/GeV) << "\n";

	generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
	<< "sHat/GeV2 = " << (lastSHat()/GeV2)
	<< " me2 = " << lastME2() << "\n" << flush;

	}

	void MatchboxMEBase::logDSigHatDR() const {

	if ( !verbose() )
	return;

	generator()->log() << "'" << name() << "' evaluated cross section using XComb "
	<< lastXCombPtr() << "\n"
	<< "Jacobian = " << jacobian() << " sHat/GeV2 = "
	<< (lastSHat()/GeV2) << " dsig/nb = "
	<< (lastMECrossSection()/nanobarn) << "\n" << flush;

	}

	void MatchboxMEBase::cloneDependencies(const std::string& prefix) {

	if ( phasespace() ) {
	Ptr<MatchboxPhasespace>::ptr myPhasespace = phasespace()->cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << myPhasespace->name();
	if ( ! (generator()->preinitRegister(myPhasespace,pname.str()) ) )
	throw InitException() << "Phasespace generator " << pname.str() << " already existing.";
	myPhasespace->cloneDependencies(pname.str());
	phasespace(myPhasespace);
	}

	theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());

	if ( matchboxAmplitude() ) {
	Ptr<MatchboxAmplitude>::ptr myAmplitude = matchboxAmplitude()->cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << myAmplitude->name();
	if ( ! (generator()->preinitRegister(myAmplitude,pname.str()) ) )
	throw InitException() << "Amplitude " << pname.str() << " already existing.";
	myAmplitude->cloneDependencies(pname.str());
	matchboxAmplitude(myAmplitude);
	amplitude(myAmplitude);
	matchboxAmplitude()->orderInGs(orderInAlphaS());
	matchboxAmplitude()->orderInGem(orderInAlphaEW());
	}

	if ( scaleChoice() ) {
	Ptr<MatchboxScaleChoice>::ptr myScaleChoice = scaleChoice()->cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << myScaleChoice->name();
	if ( ! (generator()->preinitRegister(myScaleChoice,pname.str()) ) )
	throw InitException() << "Scale choice " << pname.str() << " already existing.";
	scaleChoice(myScaleChoice);
	}

	for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
	theReweights.begin(); rw != theReweights.end(); ++rw ) {
	Ptr<MatchboxReweightBase>::ptr myReweight = (**rw).cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << (**rw).name();
	if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) )
	throw InitException() << "Reweight " << pname.str() << " already existing.";
	myReweight->cloneDependencies(pname.str());
	*rw = myReweight;
	}

	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	Ptr<MatchboxInsertionOperator>::ptr myIOP = (**v).cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << (**v).name();
	if ( ! (generator()->preinitRegister(myIOP,pname.str()) ) )
	throw InitException() << "Insertion operator " << pname.str() << " already existing.";
	*v = myIOP;
	}

	}

	void MatchboxMEBase::prepareXComb(MatchboxXCombData& xc) const {

	if ( phasespace() ) {
	size_t nout = diagrams().front()->partons().size()-2;
	xc.nDimPhasespace(phasespace()->nDim(nout));
	}

	if ( matchboxAmplitude() ) {
	xc.nDimAmplitude(matchboxAmplitude()->nDimAdditional());
	if ( matchboxAmplitude()->colourBasis() ) {
	size_t cdim =
	matchboxAmplitude()->colourBasis()->prepare(diagrams(),matchboxAmplitude()->noCorrelations());
	xc.colourBasisDim(cdim);
	}
	}

	int insertionAdd = 0;
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	insertionAdd = max(insertionAdd,(**v).nDimAdditional());
	}
	xc.nDimInsertions(insertionAdd);

	xc.nLight(getNLight());

	xc.olpId(olpProcess());

	if ( initVerbose() ) {
	string fname = name() + ".diagrams";
	ifstream test(fname.c_str());
	if ( !test ) {
	test.close();
	ofstream out(fname.c_str());
	for ( vector<Ptr<DiagramBase>::ptr>::const_iterator d = diagrams().begin();
	d != diagrams().end(); ++d ) {
	DiagramDrawer::drawDiag(out,dynamic_cast<const Tree2toNDiagram&>(**d));
	out << "\n";
	}
	}
	}

	}

	StdXCombPtr MatchboxMEBase::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
	tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
	tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
	const PBPair & newPartonBins, tCutsPtr newCuts,
	const DiagramVector & newDiagrams, bool mir,
	const PartonPairVec&,
	tStdXCombPtr newHead,
	tMEPtr newME) {

	if ( !newME )
	newME = this;

	Ptr<MatchboxXComb>::ptr xc =
	new_ptr(MatchboxXComb(newMaxEnergy, inc,
	newEventHandler, newSubProcessHandler,
	newExtractor, newCKKW,
	newPartonBins, newCuts, newME,
	newDiagrams, mir,
	newHead));

	prepareXComb(*xc);

	return xc;

	}

	StdXCombPtr MatchboxMEBase::makeXComb(tStdXCombPtr newHead,
	const PBPair & newPartonBins,
	const DiagramVector & newDiagrams,
	tMEPtr newME) {
	if ( !newME )
	newME = this;

	Ptr<MatchboxXComb>::ptr xc =
	new_ptr(MatchboxXComb(newHead, newPartonBins, newME, newDiagrams));

	prepareXComb(*xc);

	return xc;
	}

	void MatchboxMEBase::persistentOutput(PersistentOStream & os) const {
	os << theLastXComb << theFactory << thePhasespace
	<< theAmplitude << theScaleChoice << theVirtuals
	<< theReweights << theSubprocess << theOneLoop
	<< theOneLoopNoBorn
	<< epsilonSquarePoleHistograms << epsilonPoleHistograms
	<< theOLPProcess;
	}

	void MatchboxMEBase::persistentInput(PersistentIStream & is, int) {
	is >> theLastXComb >> theFactory >> thePhasespace
	>> theAmplitude >> theScaleChoice >> theVirtuals
	>> theReweights >> theSubprocess >> theOneLoop
	>> theOneLoopNoBorn
	>> epsilonSquarePoleHistograms >> epsilonPoleHistograms
	>> theOLPProcess;
	lastMatchboxXComb(theLastXComb);
	}

	void MatchboxMEBase::Init() {

	static ClassDocumentation<MatchboxMEBase> documentation
	("MatchboxMEBase is the base class for matrix elements "
	"in the context of the matchbox NLO interface.");

	}

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

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

	void MatchboxMEBase::doinit() {
	MEBase::doinit();
	if ( !theAmplitude )
	theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());
	}

	void MatchboxMEBase::dofinish() {
	MEBase::dofinish();
	for ( map<cPDVector,AccuracyHistogram>::const_iterator
	b = epsilonSquarePoleHistograms.begin();
	b != epsilonSquarePoleHistograms.end(); ++b ) {
	b->second.dump(factory()->poleData() + "epsilonSquarePoles-",b->first);
	}
	for ( map<cPDVector,AccuracyHistogram>::const_iterator
	b = epsilonPoleHistograms.begin();
	b != epsilonPoleHistograms.end(); ++b ) {
	b->second.dump(factory()->poleData() + "epsilonPoles-",b->first);
	}
	}

	// * Attention * The following static variable is needed for the type
	// description system in ThePEG. Please check that the template arguments
	// are correct (the class and its base class), and that the constructor
	// arguments are correct (the class name and the name of the dynamically
	// loadable library where the class implementation can be found).
	DescribeClass<MatchboxMEBase,MEBase>
	describeHerwigMatchboxMEBase("Herwig::MatchboxMEBase", "HwMatchbox.so");

	diff --git a/MatrixElement/Matchbox/Base/MatchboxMEBase.h b/MatrixElement/Matchbox/Base/MatchboxMEBase.h
	--- a/MatrixElement/Matchbox/Base/MatchboxMEBase.h
	+++ b/MatrixElement/Matchbox/Base/MatchboxMEBase.h
	@@ -1,957 +1,974 @@
	// -- C++ --
	//
	// MatchboxMEBase.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 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_MatchboxMEBase_H
	#define HERWIG_MatchboxMEBase_H
	//
	// This is the declaration of the MatchboxMEBase class.
	//

	#include "ThePEG/MatrixElement/MEBase.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/ProcessData.h"
	#include "Herwig++/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
	#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
	#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxReweightBase.h"
	#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
	#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
	#include "Herwig++/MatrixElement/Matchbox/InsertionOperators/MatchboxInsertionOperator.h"
	#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.fh"
	#include "Herwig++/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxXComb.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief MatchboxMEBase is the base class for matrix elements
	* in the context of the matchbox NLO interface.
	*
	* @see \ref MatchboxMEBaseInterfaces "The interfaces"
	* defined for MatchboxMEBase.
	*/
	class MatchboxMEBase:
	public MEBase, public LastMatchboxXCombInfo {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	MatchboxMEBase();

	/**
	* The destructor.
	*/
	virtual ~MatchboxMEBase();
	//@}

	public:

	/**
	* Return the factory which produced this matrix element
	*/
	Ptr<MatchboxFactory>::tcptr factory() const;

	/**
	* Set the factory which produced this matrix element
	*/
	void factory(Ptr<MatchboxFactory>::tcptr f);

	/** @name Subprocess and diagram information. */
	//@{

	/**
	* Return the subprocess.
	*/
	const Process& subProcess() const { return theSubprocess; }

	/**
	* Access the subprocess.
	*/
	Process& subProcess() { return theSubprocess; }

	/**
	* Return the diagram generator.
	*/
	Ptr<Tree2toNGenerator>::tptr diagramGenerator() const;

	/**
	* Return the process data.
	*/
	Ptr<ProcessData>::tptr processData() const;

	/**
	* Return true, if this matrix element does not want to
	* make use of mirroring processes; in this case all
	* possible partonic subprocesses with a fixed assignment
	* of incoming particles need to be provided through the diagrams
	* added with the add(...) method.
	*/
	virtual bool noMirror () const { return true; }

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

	/**
	* 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.
	*/
	virtual Selector<DiagramIndex> diagrams(const DiagramVector &) const;
	using MEBase::diagrams;

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

	/**
	* Return the order in \f$\alpha_S\f$ in which this matrix element
	* is given.
	*/
	virtual unsigned int orderInAlphaS() const;
	using MEBase::orderInAlphaS;

	/**
	* Return the order in \f$\alpha_{EM}\f$ in which this matrix
	* element is given. Returns 0.
	*/
	virtual unsigned int orderInAlphaEW() const;
	using MEBase::orderInAlphaEW;

	/**
	+ * Return true, if this amplitude already includes averaging over
	+ * incoming parton's quantum numbers.
	+ */
	+ virtual bool hasInitialAverage() const {
	+ return matchboxAmplitude() ? matchboxAmplitude()->hasInitialAverage() : false;
	+ }
	+
	+ /**
	+ * Return true, if this amplitude already includes symmetry factors
	+ * for identical outgoing particles.
	+ */
	+ virtual bool hasFinalStateSymmetry() const {
	+ return matchboxAmplitude() ? matchboxAmplitude()->hasFinalStateSymmetry() : false;
	+ }
	+
	+
	+ /**
	* Return the number of light flavours, this matrix
	* element is calculated for.
	*/
	virtual unsigned int getNLight() const;

	/**
	* Return true, if this matrix element is handled by a BLHA one-loop provider
	*/
	virtual bool isOLPTree() const {
	return matchboxAmplitude() ? matchboxAmplitude()->isOLPTree() : false;
	}

	/**
	* Return true, if this matrix element is handled by a BLHA one-loop provider
	*/
	virtual bool isOLPLoop() const {
	return matchboxAmplitude() ? matchboxAmplitude()->isOLPLoop() : false;
	}

	/**
	* Return the process index, if this is an OLP handled matrix element
	*/
	const vector<int>& olpProcess() const { return theOLPProcess; }

	/**
	* Set the process index, if this is an OLP handled matrix element
	*/
	void olpProcess(int pType, int id) {
	if ( theOLPProcess.empty() )
	theOLPProcess.resize(4,0);
	theOLPProcess[pType] = id;
	}

	//@}

	/** @name Phasespace generation */
	//@{

	/**
	* Return the phase space generator to be used.
	*/
	Ptr<MatchboxPhasespace>::tptr phasespace() const { return thePhasespace; }

	/**
	* Set the phase space generator to be used.
	*/
	void phasespace(Ptr<MatchboxPhasespace>::ptr ps) { thePhasespace = ps; }

	/**
	* Set the XComb object to be used in the next call to
	* generateKinematics() and dSigHatDR().
	*/
	virtual void setXComb(tStdXCombPtr xc);

	/**
	* Return true, if the XComb steering this matrix element
	* should keep track of the random numbers used to generate
	* the last phase space point
	*/
	virtual bool keepRandomNumbers() const { return true; }

	/**
	* Generate incoming parton momenta. This default
	* implementation performs the standard mapping
	* from x1,x2 -> tau,y making 1/tau flat; incoming
	* parton momenta are stored in meMomenta()[0,1],
	* only massless partons are supported so far;
	* return the Jacobian of the mapping
	*/
	double generateIncomingPartons(const double* r1, const double* r2);

	/**
	* Generate internal degrees of freedom given nDim() uniform random
	* numbers in the interval ]0,1[. To help the phase space generator,
	* the 'dSigHatDR' should be a smooth function of these numbers,
	* although this is not strictly necessary. The return value should
	* be true of the generation succeeded. If so the generated momenta
	* should be stored in the meMomenta() vector. Derived classes
	* must call this method once internal degrees of freedom are setup
	* and finally return the result of this method.
	*/
	virtual bool generateKinematics(const double * r);

	/**
	* The number of internal degreed of freedom used in the matrix
	* element.
	*/
	virtual int nDim() const;

	/**
	* The number of internal degrees of freedom used in the matrix
	* element for generating a Born phase space point
	*/
	virtual int nDimBorn() const;

	/**
	* Return true, if this matrix element will generate momenta for the
	* incoming partons itself. The matrix element is required to store
	* the incoming parton momenta in meMomenta()[0,1]. No mapping in
	* tau and y is performed by the PartonExtractor object, if a
	* derived class returns true here. The phase space jacobian is to
	* include a factor 1/(x1 x2).
	*/
	virtual bool haveX1X2() const {
	return
	(phasespace() ? phasespace()->haveX1X2() : false) \|\|
	diagrams().front()->partons().size() == 3;
	}

	/**
	* Return true, if this matrix element expects
	* the incoming partons in their center-of-mass system
	*/
	virtual bool wantCMS() const {
	return
	(phasespace() ? phasespace()->wantCMS() : true) &&
	diagrams().front()->partons().size() != 3; }

	/**
	* Return the meMomenta as generated at the last
	* phase space point.
	*/
	const vector<Lorentz5Momentum>& lastMEMomenta() const { return meMomenta(); }

	/**
	* Access the meMomenta.
	*/
	vector<Lorentz5Momentum>& lastMEMomenta() { return meMomenta(); }

	//@}

	/** @name Scale choices, couplings and PDFs */
	//@{

	/**
	* Set the scale choice object
	*/
	void scaleChoice(Ptr<MatchboxScaleChoice>::ptr sc) { theScaleChoice = sc; }

	/**
	* Return the scale choice object
	*/
	Ptr<MatchboxScaleChoice>::tptr scaleChoice() const { return theScaleChoice; }

	/**
	* Set scales and alphaS
	*/
	void setScale() const;

	/**
	* Return the scale associated with the phase space point provided
	* by the last call to setKinematics().
	*/
	virtual Energy2 scale() const { return lastScale(); }

	/**
	* Return the renormalization scale for the last generated phasespace point.
	*/
	virtual Energy2 factorizationScale() const;

	/**
	* Get the factorization scale factor
	*/
	virtual double factorizationScaleFactor() const;

	/**
	* Return the (QCD) renormalization scale for the last generated phasespace point.
	*/
	virtual Energy2 renormalizationScale() const;

	/**
	* Get the renormalization scale factor
	*/
	virtual double renormalizationScaleFactor() const;

	/**
	* Return the QED renormalization scale for the last generated phasespace point.
	*/
	virtual Energy2 renormalizationScaleQED() const;

	/**
	* Set veto scales on the particles at the given
	* SubProcess which has been generated using this
	* matrix element.
	*/
	virtual void setVetoScales(tSubProPtr) const;

	/**
	* Return true, if fixed couplings are used.
	*/
	bool fixedCouplings() const;

	/**
	* Return true, if fixed couplings are used.
	*/
	bool fixedQEDCouplings() const;

	/**
	* Return the value of \f$\alpha_S\f$ associated with the phase
	* space point provided by the last call to setKinematics(). This
	* versions returns SM().alphaS(scale()).
	*/
	virtual double alphaS() const { return lastAlphaS(); }

	/**
	* Return the value of \f$\alpha_EM\f$ associated with the phase
	* space point provided by the last call to setKinematics(). This
	* versions returns SM().alphaEM(scale()).
	*/
	virtual double alphaEM() const { return lastAlphaEM(); }

	/**
	* Return true, if this matrix element provides the PDF
	* weight for the first incoming parton itself.
	*/
	virtual bool havePDFWeight1() const {
	return diagrams().front()->partons()[0]->coloured();
	}

	/**
	* Return true, if this matrix element provides the PDF
	* weight for the second incoming parton itself.
	*/
	virtual bool havePDFWeight2() const {
	return diagrams().front()->partons()[1]->coloured();
	}

	/**
	* Set the PDF weight.
	*/
	void getPDFWeight(Energy2 factorizationScale = ZERO) const;

	/**
	* Supply the PDF weight for the first incoming parton.
	*/
	double pdf1(Energy2 factorizationScale = ZERO,
	double xEx = 1.) const;

	/**
	* Supply the PDF weight for the second incoming parton.
	*/
	double pdf2(Energy2 factorizationScale = ZERO,
	double xEx = 1.) const;

	//@}

	/** @name Amplitude information and matrix element evaluation */
	//@{

	/**
	* Return the amplitude.
	*/
	Ptr<MatchboxAmplitude>::tptr matchboxAmplitude() const { return theAmplitude; }

	/**
	* Set the amplitude.
	*/
	void matchboxAmplitude(Ptr<MatchboxAmplitude>::ptr amp) { theAmplitude = amp; }

	/**
	* Return the matrix element for the kinematical configuation
	* previously provided by the last call to setKinematics(), suitably
	* scaled by sHat() to give a dimension-less number.
	*/
	virtual double me2() const;

	/**
	* Return the symmetry factor for identical final state particles.
	*/
	virtual double finalStateSymmetry() const;

	/**
	* Return the normalizing factor for the matrix element averaged
	* over quantum numbers and including running couplings.
	*/
	double me2Norm(unsigned int addAlphaS = 0) const;

	/**
	* Return the matrix element squared differential in the variables
	* given by the last call to generateKinematics().
	*/
	virtual CrossSection dSigHatDR() const;

	//@}

	/** @name One-loop corrections */
	//@{

	/**
	* Return the one-loop/tree interference.
	*/
	virtual double oneLoopInterference() const;

	/**
	* Return true, if this matrix element is capable of calculating
	* one-loop (QCD) corrections.
	*/
	virtual bool haveOneLoop() const;

	/**
	* Return true, if this matrix element only provides
	* one-loop (QCD) corrections.
	*/
	virtual bool onlyOneLoop() const;

	/**
	* Return true, if one loop corrections have been calculated in
	* dimensional reduction. Otherwise conventional dimensional
	* regularization is assumed. Note that renormalization is always
	* assumed to be MSbar.
	*/
	virtual bool isDR() const;

	/**
	* Return true, if one loop corrections are given in the conventions
	* of the integrated dipoles.
	*/
	virtual bool isCS() const;

	/**
	* Return true, if one loop corrections are given in the conventions
	* of BDK.
	*/
	virtual bool isBDK() const;

	/**
	* Return true, if one loop corrections are given in the conventions
	* of everything expanded.
	*/
	virtual bool isExpanded() const;

	/**
	* Return the value of the dimensional regularization
	* parameter. Note that renormalization scale dependence is fully
	* restored in DipoleIOperator.
	*/
	virtual Energy2 mu2() const;

	/**
	* If defined, return the coefficient of the pole in epsilon^2
	*/
	virtual double oneLoopDoublePole() const;

	/**
	* If defined, return the coefficient of the pole in epsilon
	*/
	virtual double oneLoopSinglePole() const;

	/**
	* Return true, if cancellationn of epsilon poles should be checked.
	*/
	bool checkPoles() const;

	/**
	* Simple histogram for accuracy checks
	*/
	struct AccuracyHistogram {

	/**
	* The lower bound
	*/
	double lower;

	/**
	* The upper bound
	*/
	double upper;

	/**
	* The bins, indexed by upper bound.
	*/
	map<double,double> bins;

	/**
	* The number of points of same sign
	*/
	unsigned long sameSign;

	/**
	* The number of points of opposite sign
	*/
	unsigned long oppositeSign;

	/**
	* The number of points being nan or inf
	*/
	unsigned long nans;

	/**
	* The overflow
	*/
	unsigned long overflow;

	/**
	* The underflow
	*/
	unsigned long underflow;

	/**
	* Constructor
	*/
	AccuracyHistogram(double low = -40.,
	double up = 0.,
	unsigned int nbins = 80);

	/**
	* Book two values to be checked for numerical compatibility
	*/
	void book(double a, double b);

	/**
	* Write to file.
	*/
	void dump(const std::string& prefix,
	const cPDVector& proc) const;

	/**
	* Write to persistent ostream
	*/
	void persistentOutput(PersistentOStream&) const;

	/**
	* Read from persistent istream
	*/
	void persistentInput(PersistentIStream&);

	};

	/**
	* Perform the check of epsilon pole cancellation.
	*/
	void logPoles() const;

	/**
	* Return the virtual corrections
	*/
	const vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() const {
	return theVirtuals;
	}

	/**
	* Return the virtual corrections
	*/
	vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() {
	return theVirtuals;
	}

	/**
	* Instruct this matrix element to include one-loop corrections
	*/
	void doOneLoop() { theOneLoop = true; }

	/**
	* Return true, if this matrix element includes one-loop corrections
	*/
	bool oneLoop() const { return theOneLoop; }

	/**
	* Instruct this matrix element to include one-loop corrections but
	* no Born contributions
	*/
	void doOneLoopNoBorn() { theOneLoop = true; theOneLoopNoBorn = true; }

	/**
	* Return true, if this matrix element includes one-loop corrections
	* but no Born contributions
	*/
	bool oneLoopNoBorn() const { return theOneLoopNoBorn \|\| onlyOneLoop(); }

	//@}

	/** @name Dipole subtraction */
	//@{

	/**
	* If this matrix element is considered a real
	* emission matrix element, return all subtraction
	* dipoles needed given a set of subtraction terms
	* and underlying Born matrix elements to choose
	* from.
	*/
	vector<Ptr<SubtractionDipole>::ptr>
	getDipoles(const vector<Ptr<SubtractionDipole>::ptr>&,
	const vector<Ptr<MatchboxMEBase>::ptr>&) const;

	/**
	* If this matrix element is considered a real emission matrix
	* element, but actually neglecting a subclass of the contributing
	* diagrams, return true if the given emitter-emission-spectator
	* configuration should not be considered when setting up
	* subtraction dipoles.
	*/
	virtual bool noDipole(int,int,int) const { return false; }

	/**
	* If this matrix element is considered an underlying Born matrix
	* element in the context of a subtracted real emission, but
	* actually neglecting a subclass of the contributing diagrams,
	* return true if the given emitter-spectator configuration
	* should not be considered when setting up subtraction dipoles.
	*/
	virtual bool noDipole(int,int) const { return false; }

	/**
	* Return the colour correlated matrix element squared with
	* respect to the given two partons as appearing in mePartonData(),
	* suitably scaled by sHat() to give a dimension-less number.
	*/
	virtual double colourCorrelatedME2(pair<int,int>) const;

	/**
	* Return the colour correlated matrix element squared in the
	* large-N approximation with respect to the given two partons as
	* appearing in mePartonData(), suitably scaled by sHat() to give a
	* dimension-less number.
	*/
	virtual double largeNColourCorrelatedME2(pair<int,int> ij,
	Ptr<ColourBasis>::tptr largeNBasis) const;

	/**
	* Return the colour and spin correlated matrix element squared for
	* the gluon indexed by the first argument using the given
	* correlation tensor.
	*/
	virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
	const SpinCorrelationTensor& c) const;

	//@}

	/** @name Caching and diagnostic information */
	//@{

	/**
	* Inform this matrix element that a new phase space
	* point is about to be generated, so all caches should
	* be flushed.
	*/
	virtual void flushCaches();

	/**
	* Return true, if verbose
	*/
	bool verbose() const;

	/**
	* Return true, if verbose
	*/
	bool initVerbose() const;

	/**
	* Dump the setup to an ostream
	*/
	void print(ostream&) const;

	/**
	* Print debug information on the last event
	*/
	virtual void printLastEvent(ostream&) const;

	/**
	* Write out diagnostic information for
	* generateKinematics
	*/
	void logGenerateKinematics(const double * r) const;

	/**
	* Write out diagnostic information for
	* setting scales
	*/
	void logSetScale() const;

	/**
	* Write out diagnostic information for
	* pdf evaluation
	*/
	void logPDFWeight() const;

	/**
	* Write out diagnostic information for
	* me2 evaluation
	*/
	void logME2() const;

	/**
	* Write out diagnostic information
	* for dsigdr evaluation
	*/
	void logDSigHatDR() const;

	//@}

	/** @name Reweight objects */
	//@{

	/**
	* Insert a reweight object
	*/
	void addReweight(Ptr<MatchboxReweightBase>::ptr rw) { theReweights.push_back(rw); }

	/**
	* Return the reweights
	*/
	const vector<Ptr<MatchboxReweightBase>::ptr>& reweights() const { return theReweights; }

	/**
	* Access the reweights
	*/
	vector<Ptr<MatchboxReweightBase>::ptr>& reweights() { return theReweights; }

	//@}

	/** @name Methods used to setup MatchboxMEBase objects */
	//@{

	/**
	* Return true if this object needs to be initialized before all
	* other objects (except those for which this function also returns
	* true). This default version always returns false, but subclasses
	* may override it to return true.
	*/
	virtual bool preInitialize() const { return true; }

	/**
	* Clone this matrix element.
	*/
	Ptr<MatchboxMEBase>::ptr cloneMe() const {
	return dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(clone());
	}

	/**
	* Clone the dependencies, using a given prefix.
	*/
	void cloneDependencies(const std::string& prefix = "");

	/**
	* Prepare an xcomb
	*/
	void prepareXComb(MatchboxXCombData&) const;

	/**
	* For the given event generation setup return a xcomb object
	* appropriate to this matrix element.
	*/
	virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc,
	tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
	tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
	const PBPair & newPartonBins, tCutsPtr newCuts,
	const DiagramVector & newDiagrams, bool mir,
	const PartonPairVec& allPBins,
	tStdXCombPtr newHead = tStdXCombPtr(),
	tMEPtr newME = tMEPtr());

	/**
	* For the given event generation setup return a dependent xcomb object
	* appropriate to this matrix element.
	*/
	virtual StdXCombPtr makeXComb(tStdXCombPtr newHead,
	const PBPair & newPartonBins,
	const DiagramVector & newDiagrams,
	tMEPtr newME = tMEPtr());

	//@}

	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();

	/**
	* Finalize this object. Called in the run phase just after a
	* run has ended. Used eg. to write out statistics.
	*/
	virtual void dofinish();
	//@}

	private:

	/**
	* The factory which produced this matrix element
	*/
	Ptr<MatchboxFactory>::tcptr theFactory;

	/**
	* The phase space generator to be used.
	*/
	Ptr<MatchboxPhasespace>::ptr thePhasespace;

	/**
	* The amplitude to be used
	*/
	Ptr<MatchboxAmplitude>::ptr theAmplitude;

	/**
	* The scale choice object
	*/
	Ptr<MatchboxScaleChoice>::ptr theScaleChoice;

	/**
	* The virtual corrections.
	*/
	vector<Ptr<MatchboxInsertionOperator>::ptr> theVirtuals;

	/**
	* A vector of reweight objects the sum of which
	* should be applied to reweight this matrix element
	*/
	vector<Ptr<MatchboxReweightBase>::ptr> theReweights;

	private:

	/**
	* The subprocess to be considered.
	*/
	Process theSubprocess;

	/**
	* True, if this matrix element includes one-loop corrections
	*/
	bool theOneLoop;

	/**
	* True, if this matrix element includes one-loop corrections
	* but no Born contributions
	*/
	bool theOneLoopNoBorn;

	/**
	* The process index, if this is an OLP handled matrix element
	*/
	vector<int> theOLPProcess;

	/**
	* Histograms of epsilon^2 pole cancellation
	*/
	mutable map<cPDVector,AccuracyHistogram> epsilonSquarePoleHistograms;

	/**
	* Histograms of epsilon pole cancellation
	*/
	mutable map<cPDVector,AccuracyHistogram> epsilonPoleHistograms;

	private:

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

	};

	inline PersistentOStream& operator<<(PersistentOStream& os,
	const MatchboxMEBase::AccuracyHistogram& h) {
	h.persistentOutput(os);
	return os;
	}

	inline PersistentIStream& operator>>(PersistentIStream& is,
	MatchboxMEBase::AccuracyHistogram& h) {
	h.persistentInput(is);
	return is;
	}

	}

	#endif /* HERWIG_MatchboxMEBase_H */

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 12:30 PM (1 d, 20 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4022778
Default Alt Text: (79 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions