No OneTemporary
Actions

Size

214 KB

Subscribers

None

View Options

	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,1214 +1,1218 @@
	// -- 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/Base/MatchboxAmplitude.h"
	#include "Herwig/MatrixElement/Matchbox/Base/MatchboxReweightBase.h"
	#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
	#include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h"
	#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>::tptr factory() const;

	/**
	* Set the factory which produced this matrix element
	*/
	void factory(Ptr<MatchboxFactory>::tptr 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 true, if this amplitude is capable of consistently filling
	* the rho matrices for the spin correllations
	*/
	virtual bool canFillRhoMatrix() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->canFillRhoMatrix();
	return false;
	}

	/**
	* construct the spin information for the interaction
	*/
	virtual void constructVertex(tSubProPtr) {}

	/**
	* construct the spin information for the interaction
	*/
	virtual void constructVertex(tSubProPtr sub, const ColourLines* cl);

	/**
	* 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 the vector that contains the PDG ids of
	* the light flavours, which are contained in the
	* jet particle group.
	*/
	virtual vector<int> getNLightJetVec() const;

	/**
	* Return the vector that contains the PDG ids of
	* the heavy flavours, which are contained in the
	* jet particle group.
	*/
	virtual vector<int> getNHeavyJetVec() const;

	/**
	* Return the vector that contains the PDG ids of
	* the light flavours, which are contained in the
	* proton particle group.
	*/
	virtual vector<int> getNLightProtonVec() 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 true, if colour and spin correlated matrix elements should
	* be ordered from the OLP
	*/
	virtual bool needsOLPCorrelators() const {
	return matchboxAmplitude() ? matchboxAmplitude()->needsOLPCorrelators() : true;
	}

	/**
	* 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(5,0);
	theOLPProcess[pType] = id;
	}

	/**
	* Return true, if this is a real emission matrix element which does
	* not require colour correlators.
	*/
	bool noCorrelations() const {
	return theNoCorrelations;
	}

	/**
	* Indicate that this is a real emission matrix element which does
	* not require colour correlators.
	*/
	void needsNoCorrelations() {
	theNoCorrelations = true;
	}

	/**
	* Indicate that this is a virtual matrix element which does
	* require colour correlators.
	*/
	void needsCorrelations() {
	theNoCorrelations = false;
	}

	//@}

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

	/**
	* Set the typed and momenta of the incoming and outgoing partons to
	* be used in subsequent calls to me() and colourGeometries()
	* according to the associated XComb object. If the function is
	* overridden in a sub class the new function must call the base
	* class one first.
	*/
	virtual void setKinematics();

	/**
	* Clear the information previously provided by a call to
	* setKinematics(...).
	*/
	virtual void clearKinematics();

	/**
	* 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(); }


	/**
	* leg size
	*/

	int legsize() const {return int(meMomenta().size());}

	//@}

	/** @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(Energy2 ren=ZERO,Energy2 fac=ZERO) const;

	/**
	* Indicate that this matrix element is running alphas by itself.
	*/
	virtual bool hasRunningAlphaS() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->hasRunningAlphaS();
	return false;
	}

	/**
	* Indicate that this matrix element is running alphaew by itself.
	*/
	virtual bool hasRunningAlphaEW() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->hasRunningAlphaEW();
	return false;
	}

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

	/**
	* Return the shower scale for the last generated phasespace point.
	*/
	virtual Energy2 showerScale() 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 true, if this matrix element provides the PDF
	* weight for the second incoming parton itself.
	*/
	virtual bool havePDFWeight2() const;

	/**
	* 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., double xFactor = 1.) const;

	/**
	* Supply the PDF weight for the second incoming parton.
	*/
	double pdf2(Energy2 factorizationScale = ZERO,
	double xEx = 1., double xFactor = 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 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 largeNME2(Ptr<ColourBasis>::tptr largeNBasis) 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;

	/**
	* Always same prefactor
	**/
	CrossSection prefactor()const;

	/**
	* Born part of the cross section
	**/
	CrossSection dSigHatDRB() const ;

	/**
	* Virtual corrections of the cross section
	**/
	CrossSection dSigHatDRV() const ;

	/**
	* Insertion operators of the cross section
	**/
	CrossSection dSigHatDRI() const ;

	/**
	* If diffAlpha is not 1 and the matrix element has insertion operators
	* this routine adds the difference between the insertion operator calculated
	* with an alpha-Parameter to the insertion operator without alpha-parameter.
	*/
	CrossSection dSigHatDRAlphaDiff(double alpha) 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 the amplitude is DRbar renormalized, otherwise
	* MSbar is assumed.
	*/
	virtual bool isDRbar() 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& folder, 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; }

	+ /**
	+ * Instruct this matrix element not to include one-loop corrections
	+ */
	+
	void noOneLoop() { theOneLoop = false; }

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

	void noOneLoopNoBorn() { theOneLoop = false; theOneLoopNoBorn = false; }

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

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

	/**
	* Return true, if this matrix element includes one-loop corrections
	* but no actual loop contributions
	*/
	bool oneLoopNoLoops() const { return theOneLoopNoLoops; }

	//@}

	/** @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>&,bool slim=false) 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;

	/**
	* Return the spin correlated matrix element squared for
	* the vector boson indexed by the first argument using the given
	* correlation tensor.
	*/
	virtual double spinCorrelatedME2(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; }

	/**
	* Return the theMerger.
	*/
	const MergerBasePtr merger() const;

	/**
	* Return the theMerger.
	*/
	MergerBasePtr merger() ;

	/**
	* Set the theMerger.
	*/
	void merger(MergerBasePtr v);

	/**
	* Fill the projectors
	*/
	virtual void fillProjectors();

	///
	pair<bool,bool> clustersafe(int emit,int emis,int spec);

	//@}

	/** @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 = ""){cloneDependencies(prefix,false);};

	void cloneDependencies(const std::string& prefix ,bool slim);

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

	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();

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

	/**
	* True, if this matrix element includes one-loop corrections
	* but no actual loop contributions (e.g. finite collinear terms)
	*/
	bool theOneLoopNoLoops;

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

	/**
	* True, if this is a real emission matrix element which does
	* not require colour correlators.
	*/
	bool theNoCorrelations;

	/**
	* Flag which pdfs should be included.
	*/
	mutable pair<bool,bool> theHavePDFs;

	/**
	* True, if already checked for which PDFs to include.
	*/
	mutable bool checkedPDFs;

	/**
	* The first cluster node
	*/

	MergerBasePtr theMerger;


	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 */
	diff --git a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
	--- a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
	+++ b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
	@@ -1,1373 +1,1347 @@
	// -- C++ --
	//
	// SubtractionDipole.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 SubtractionDipole class.
	//

	#include "SubtractionDipole.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Repository/Repository.h"
	#include "ThePEG/StandardModel/StandardModelBase.h"
	#include "ThePEG/Utilities/Rebinder.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDF/PartonBin.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
	#include "Herwig/MatrixElement/Matchbox/Utility/DiagramDrawer.h"

	#include <iterator>
	using std::ostream_iterator;

	using namespace Herwig;

	SubtractionDipole::SubtractionDipole()
	: MEBase(), theSplitting(false), theApply(true), theSubtractionTest(false),
	theIgnoreCuts(false),
	theRealEmitter(-1), theRealEmission(-1), theRealSpectator(-1),
	lastRealEmissionKey(realEmissionKey(cPDVector(),-1,-1,-1)),
	lastUnderlyingBornKey(underlyingBornKey(cPDVector(),-1,-1)),
	theBornEmitter(-1), theBornSpectator(-1),
	theLastSubtractionScale(ZERO), theLastSplittingScale(ZERO),
	theLastSubtractionPt(ZERO), theLastSplittingPt(ZERO),
	theLastSubtractionZ(0.0), theLastSplittingZ(0.0),
	theRealShowerSubtraction(false), theVirtualShowerSubtraction(false),
	theLoopSimSubtraction(false), theRealEmissionScales(false),
	theShowerHardScale(ZERO), theShowerScale(ZERO),
	theIsInShowerPhasespace(false), theIsAboveCutoff(false) {}

	SubtractionDipole::~SubtractionDipole() {}

	double SubtractionDipole::alpha() const{
	return factory()->alphaParameter();
	}

	void SubtractionDipole::clearBookkeeping() {
	theRealEmitter = -1;
	theRealEmission = -1;
	theRealSpectator = -1;
	theBornEmitter = -1;
	theBornSpectator = -1;
	theMergingMap.clear();
	theSplittingMap.clear();
	theIndexMap.clear();
	theUnderlyingBornDiagrams.clear();
	theRealEmissionDiagrams.clear();
	theBornToRealDiagrams.clear();
	theRealToBornDiagrams.clear();
	}

	void SubtractionDipole::setupBookkeeping(const map<Ptr<DiagramBase>::ptr,SubtractionDipole::MergeInfo>& mergeInfo,bool slim) {

	theMergingMap.clear();
	theSplittingMap.clear();

	theUnderlyingBornDiagrams.clear();
	theRealEmissionDiagrams.clear();

	theBornToRealDiagrams.clear();
	theRealToBornDiagrams.clear();

	int xemitter = -1;
	int xspectator = -1;
	map<int,int> mergeLegs;
	map<int,int> remapLegs;
	map<int,int> realBornMap;
	map<int,int> bornRealMap;

	set<Ptr<DiagramBase>::cptr> usedDiagrams;

	for ( map<Ptr<DiagramBase>::ptr,MergeInfo>::const_iterator mit = mergeInfo.begin();
	mit != mergeInfo.end(); ++mit ) {

	DiagramVector::const_iterator bd =
	theUnderlyingBornME->diagrams().end();

	// work out the most similar underlying Born diagram
	map<int,int> xRemapLegs;
	int nomapScore = 0;
	for ( DiagramVector::const_iterator b =
	theUnderlyingBornME->diagrams().begin();
	b != theUnderlyingBornME->diagrams().end(); ++b ) {
	map<int,int> theRemapLegs;
	if ( mit->second.diagram->isSame(*b,theRemapLegs) &&
	usedDiagrams.find(*b) == usedDiagrams.end() ) {
	int theNomapScore = 0;
	for ( map<int,int>::const_iterator m = theRemapLegs.begin();
	m != theRemapLegs.end(); ++m )
	if ( m->first == m->second )
	theNomapScore += 1;
	if ( theNomapScore >= nomapScore ) {
	nomapScore = theNomapScore;
	xRemapLegs = theRemapLegs;
	bd = b;
	}
	}
	}

	// no underlying Born
	if ( bd == theUnderlyingBornME->diagrams().end() )
	continue;

	// as we deal with one splitting only we now mark this diagram as used
	// since we fixed the overall remapping of the process from the first
	// occurence, see below. TODO: This confuses this code even more, and
	// clearly calls for a cleanup. This is just grown historically and got
	// messed up with experiencing different processes and setups.
	usedDiagrams.insert(*bd);

	if ( xemitter == -1 ) {

	xemitter = mit->second.emitter;
	mergeLegs = mit->second.mergeLegs;
	remapLegs = xRemapLegs;

	assert(remapLegs.find(xemitter) != remapLegs.end());
	xemitter = remapLegs[xemitter];

	// work out the leg remapping real -> born
	for ( map<int,int>::const_iterator k = mergeLegs.begin();
	k != mergeLegs.end(); ++k ) {
	assert(remapLegs.find(k->second) != remapLegs.end());
	realBornMap[k->first] = remapLegs[k->second];
	}

	// work out the leg remapping born -> real
	for ( map<int,int>::const_iterator k = realBornMap.begin();
	k != realBornMap.end(); ++k ) {
	bornRealMap[k->second] = k->first;
	}

	// work out the spectator
	assert(mergeLegs.find(realSpectator()) != mergeLegs.end());
	assert(remapLegs.find(mergeLegs[realSpectator()]) != remapLegs.end());
	xspectator = realBornMap[realSpectator()];

	}

	RealEmissionKey realKey = realEmissionKey((*mit->first).partons(),realEmitter(),realEmission(),realSpectator());
	UnderlyingBornKey bornKey = underlyingBornKey((**bd).partons(),xemitter,xspectator);
	if ( theMergingMap.find(realKey) == theMergingMap.end() )
	theMergingMap.insert(make_pair(realKey,make_pair(bornKey,realBornMap)));
	RealEmissionInfo realInfo = make_pair(realKey,bornRealMap);
	bool gotit = false;
	typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spIterator;
	pair<spIterator,spIterator> range = theSplittingMap.equal_range(bornKey);
	for ( ; range.first != range.second; ++range.first )
	if ( range.first->second == realInfo ) {
	gotit = true;
	break;
	}
	if ( !gotit )
	theSplittingMap.insert(make_pair(bornKey,realInfo));
	theUnderlyingBornDiagrams[process(realKey)].push_back(*bd);
	theRealEmissionDiagrams[process(bornKey)].push_back(mit->first);
	theBornToRealDiagrams[*bd] = mit->first;
	theRealToBornDiagrams[mit->first] = *bd;

	}


	if (slim) {
	theIndexMap.clear();
	theSplittingMap.clear();
	theBornToRealDiagrams.clear();
	theRealEmissionDiagrams.clear();
	}

	if ( theSplittingMap.empty() )
	return;

	theIndexMap.clear();

	for ( multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator s =
	theSplittingMap.begin(); s != theSplittingMap.end(); ++s ) {
	theIndexMap[process(s->first)] = make_pair(emitter(s->first),spectator(s->first));
	}

	}

	void SubtractionDipole::subtractionBookkeeping() {
	/*
	if ( theMergingMap.empty() )
	setupBookkeeping();
	*/
	assert(!theMergingMap.empty());
	lastRealEmissionKey =
	realEmissionKey(lastHeadXComb().mePartonData(),realEmitter(),realEmission(),realSpectator());
	map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
	theMergingMap.find(lastRealEmissionKey);
	if ( k == theMergingMap.end() ) {
	theApply = false;
	return;
	}
	theApply = true;
	lastUnderlyingBornKey = k->second.first;
	bornEmitter(emitter(lastUnderlyingBornKey));
	bornSpectator(spectator(lastUnderlyingBornKey));
	}

	void SubtractionDipole::splittingBookkeeping() {
	/*
	if ( theMergingMap.empty() )
	setupBookkeeping();
	*/
	assert(!theMergingMap.empty());
	map<cPDVector,pair<int,int> >::const_iterator esit =
	theIndexMap.find(lastHeadXComb().mePartonData());
	if ( esit == theIndexMap.end() ) {
	theApply = false;
	return;
	}
	theApply = true;
	pair<int,int> es = esit->second;
	bornEmitter(es.first);
	bornSpectator(es.second);
	lastUnderlyingBornKey = underlyingBornKey(lastHeadXComb().mePartonData(),bornEmitter(),bornSpectator());
	typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
	pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
	assert(kr.first != kr.second);
	lastRealEmissionInfo = kr.first;
	for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
	if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
	break;
	assert(lastRealEmissionInfo != kr.second);
	lastRealEmissionKey = lastRealEmissionInfo->second.first;
	realEmitter(emitter(lastRealEmissionKey));
	realEmission(emission(lastRealEmissionKey));
	realSpectator(spectator(lastRealEmissionKey));
	}

	StdXCombPtr SubtractionDipole::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& allBins,
	tStdXCombPtr newHead,
	tMEPtr newME) {

	if ( !newME )
	newME = this;

	if ( !splitting() ) {
	return
	underlyingBornME()->makeXComb(newMaxEnergy, inc,
	newEventHandler, newSubProcessHandler,
	newExtractor, newCKKW,
	newPartonBins, newCuts,
	newDiagrams, mir, allBins,
	newHead, newME);
	}

	return
	realEmissionME()->makeXComb(newMaxEnergy, inc,
	newEventHandler, newSubProcessHandler,
	newExtractor, newCKKW,
	newPartonBins, newCuts,
	newDiagrams, mir, allBins,
	newHead, newME);

	}

	StdXCombPtr SubtractionDipole::makeXComb(tStdXCombPtr newHead,
	const PBPair & newPartonBins,
	const DiagramVector & newDiagrams,
	tMEPtr newME) {

	if ( !newME )
	newME = this;

	if ( !splitting() ) {
	return
	underlyingBornME()->makeXComb(newHead, newPartonBins,
	newDiagrams, newME);
	}

	return
	realEmissionME()->makeXComb(newHead, newPartonBins,
	newDiagrams, newME);

	}

	StdXCombPtr SubtractionDipole::makeBornXComb(tStdXCombPtr realXC) {

	const cPDVector& proc = const_cast<const StandardXComb&>(*realXC).mePartonData();

	lastRealEmissionKey =
	realEmissionKey(proc,realEmitter(),realEmission(),realSpectator());
	map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
	theMergingMap.find(lastRealEmissionKey);

	if ( k == theMergingMap.end() )
	return StdXCombPtr();

	PartonPairVec pbs = realXC->pExtractor()->getPartons(realXC->maxEnergy(),
	realXC->particles(),
	*(realXC->cuts()));

	DiagramVector bornDiags = underlyingBornDiagrams(proc);
	assert(!bornDiags.empty());

	PartonPairVec::iterator ppit = pbs.begin();
	for ( ; ppit != pbs.end(); ++ppit ) {
	if ( ppit->first->parton() == bornDiags.front()->partons()[0] &&
	ppit->second->parton() == bornDiags.front()->partons()[1] )
	break;
	}

	assert(ppit != pbs.end());

	return
	underlyingBornME()->makeXComb(realXC,*ppit,bornDiags,this);

	}

	vector<StdXCombPtr> SubtractionDipole::makeRealXCombs(tStdXCombPtr bornXC) {

	const cPDVector& proc = const_cast<const StandardXComb&>(*bornXC).mePartonData();

	map<cPDVector,pair<int,int> >::const_iterator esit = theIndexMap.find(proc);
	if ( esit == theIndexMap.end() )
	return vector<StdXCombPtr>();
	pair<int,int> es = esit->second;
	bornEmitter(es.first);
	bornSpectator(es.second);
	lastUnderlyingBornKey = underlyingBornKey(proc,bornEmitter(),bornSpectator());

	if ( theSplittingMap.find(lastUnderlyingBornKey) == theSplittingMap.end() )
	return vector<StdXCombPtr>();

	PartonPairVec pbs = bornXC->pExtractor()->getPartons(bornXC->maxEnergy(),
	bornXC->particles(),
	*(bornXC->cuts()));

	DiagramVector realDiags = realEmissionDiagrams(proc);
	assert(!realDiags.empty());

	vector<StdXCombPtr> res;

	map<cPDVector,DiagramVector> realProcs;

	for ( MEBase::DiagramVector::const_iterator d = realDiags.begin();
	d != realDiags.end(); ++d ) {
	realProcs[(*d).partons()].push_back(d);
	}

	for ( map<cPDVector,DiagramVector>::const_iterator pr =
	realProcs.begin(); pr != realProcs.end(); ++pr ) {

	PartonPairVec::iterator ppit = pbs.begin();
	for ( ; ppit != pbs.end(); ++ppit ) {
	if ( ppit->first->parton() == pr->second.front()->partons()[0] &&
	ppit->second->parton() == pr->second.front()->partons()[1] )
	break;
	}

	assert(ppit != pbs.end());

	StdXCombPtr rxc =
	realEmissionME()->makeXComb(bornXC,*ppit,pr->second,this);

	res.push_back(rxc);

	}

	return res;

	}

	const MEBase::DiagramVector& SubtractionDipole::underlyingBornDiagrams(const cPDVector& real) const {
	static DiagramVector empty;
	map<cPDVector,DiagramVector>::const_iterator k = theUnderlyingBornDiagrams.find(real);
	if (k == theUnderlyingBornDiagrams.end() )
	return empty;
	return k->second;
	}

	tcDiagPtr SubtractionDipole::underlyingBornDiagram(tcDiagPtr realDiag) const {
	map<tcDiagPtr,tcDiagPtr>::const_iterator it = theRealToBornDiagrams.find(realDiag);
	assert(it != theRealToBornDiagrams.end());
	return it->second;
	}

	const MEBase::DiagramVector& SubtractionDipole::realEmissionDiagrams(const cPDVector& born) const {
	static DiagramVector empty;
	map<cPDVector,DiagramVector>::const_iterator k = theRealEmissionDiagrams.find(born);
	if ( k == theRealEmissionDiagrams.end() )
	return empty;
	return k->second;
	}

	tcDiagPtr SubtractionDipole::realEmissionDiagram(tcDiagPtr bornDiag) const {
	map<tcDiagPtr,tcDiagPtr>::const_iterator it = theBornToRealDiagrams.find(bornDiag);
	assert(it != theBornToRealDiagrams.end());
	return it->second;
	}

	void SubtractionDipole::getDiagrams() const {
	if ( splitting() ) {
	realEmissionME()->diagrams();
	useDiagrams(realEmissionME());
	} else {
	underlyingBornME()->diagrams();
	useDiagrams(underlyingBornME());
	}
	}

	Selector<MEBase::DiagramIndex> SubtractionDipole::diagrams(const DiagramVector & dv) const {
	Ptr<MatchboxMEBase>::tcptr me =
	splitting() ?
	realEmissionME() :
	underlyingBornME();
	if ( me->phasespace() ) {
	me->phasespace()->setXComb(lastXCombPtr());
	me->phasespace()->fillDiagramWeights();
	}
	return
	me->diagrams(dv);
	}

	Selector<const ColourLines *>
	SubtractionDipole::colourGeometries(tcDiagPtr diag) const {
	return
	splitting() ?
	realEmissionME()->colourGeometries(diag) :
	underlyingBornME()->colourGeometries(diag);
	}

	const ColourLines &
	SubtractionDipole::selectColourGeometry(tcDiagPtr diag) const {
	return
	splitting() ?
	realEmissionME()->selectColourGeometry(diag) :
	underlyingBornME()->selectColourGeometry(diag);
	}

	void SubtractionDipole::flushCaches() {
	theUnderlyingBornME->flushCaches();
	theRealEmissionME->flushCaches();
	for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
	reweights().begin(); r != reweights().end(); ++r ) {
	(**r).flushCaches();
	}
	}

	void SubtractionDipole::setXComb(tStdXCombPtr xc) {
	if ( !xc ) {
	theApply = false;
	return;
	} else {
	theApply = true;
	}
	lastMatchboxXComb(xc);
	MEBase::setXComb(xc);
	if ( splitting() ) {
	realEmissionME()->setXComb(xc);
	underlyingBornME()->setXComb(xc->head());
	splittingBookkeeping();
	} else {
	realEmissionME()->setXComb(xc->head());
	underlyingBornME()->setXComb(xc);
	subtractionBookkeeping();
	}
	if ( !apply() )
	return;
	}

	void SubtractionDipole::setKinematics() {
	MEBase::setKinematics();
	if ( splitting() )
	realEmissionME()->setKinematics();
	else
	underlyingBornME()->setKinematics();
	}

	bool SubtractionDipole::generateKinematics(const double * r) {
	if ( lastXCombPtr()->kinematicsGenerated() )
	return true;
	if ( splitting() ) {
	if ( !generateRadiationKinematics(r) )
	return false;
	if( ! realEmissionME()->lastXCombPtr()->setIncomingPartons())
	return false;
	realEmissionME()->setScale();
	double jac = jacobian();
	jac *= pow(underlyingBornME()->lastXComb().lastSHat() / realEmissionME()->lastXComb().lastSHat(),
	realEmissionME()->lastXComb().mePartonData().size()-4.);
	jacobian(jac);
	assert(lastXCombPtr() == realEmissionME()->lastXCombPtr());
	lastXCombPtr()->didGenerateKinematics();
	return true;
	}
	- if ( !generateTildeKinematics() )
	- {cout<<"\nxx 1 xx"; return false;}
	+ if ( !generateTildeKinematics() ){ return false;}
	if( ! underlyingBornME()->lastXCombPtr()->setIncomingPartons() )
	- {cout<<"\nxx 2 xx"; return false;}
	+ return false;
	underlyingBornME()->setScale();
	assert(lastXCombPtr() == underlyingBornME()->lastXCombPtr());
	if( ! underlyingBornME()->lastXCombPtr()->setIncomingPartons() )
	- {cout<<"\nxx 3 xx"; return false;}
	+ return false;
	// need to have the scale and x's available for checking shower phase space
	if ( showerApproximation() &&
	lastXCombPtr()->willPassCuts() )
	showerApproximation()->getShowerVariables();
	lastXCombPtr()->didGenerateKinematics();
	return true;
	}

	int SubtractionDipole::nDim() const {
	if ( !splitting() )
	return underlyingBornME()->nDim();
	return underlyingBornME()->nDim() + nDimRadiation();
	}

	void SubtractionDipole::clearKinematics() {
	MEBase::clearKinematics();
	if ( splitting() )
	realEmissionME()->clearKinematics();
	else
	underlyingBornME()->clearKinematics();
	}

	void SubtractionDipole::tildeKinematics(Ptr<TildeKinematics>::tptr tk) {
	theTildeKinematics = tk;
	}

	bool SubtractionDipole::generateTildeKinematics() {

	assert(!splitting());

	Ptr<TildeKinematics>::tptr kinematics = theTildeKinematics;
	if ( showerApproximation() ) {
	showerApproximation()->setBornXComb(lastXCombPtr());
	showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
	showerApproximation()->setDipole(this);
	showerApproximation()->checkCutoff();
	if ( showerApproximation()->showerTildeKinematics() &&
	isAboveCutoff() &&
	realShowerSubtraction() )
	kinematics = showerApproximation()->showerTildeKinematics();
	}

	if ( !kinematics ) {
	jacobian(0.0);
	return false;
	}

	kinematics->prepare(lastHeadXCombPtr(),lastXCombPtr());

	if ( !kinematics->doMap() ) {
	jacobian(0.0);
	return false;
	}

	theLastSubtractionScale = kinematics->lastScale();
	theLastSubtractionPt = kinematics->lastPt();
	theLastSubtractionZ = kinematics->lastZ();

	meMomenta().resize(lastHeadXComb().meMomenta().size() - 1);

	assert(mergingMap().find(lastRealEmissionKey) != mergingMap().end());
	map<int,int>& trans = theMergingMap[lastRealEmissionKey].second;

	int n = lastHeadXComb().meMomenta().size();
	for ( int k = 0; k < n; ++k ) {
	if ( k == realEmitter() \|\| k == realEmission() \|\| k == realSpectator() )
	continue;
	meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
	if ( kinematics->doesTransform() && k > 1 )
	meMomenta()[trans[k]] = kinematics->transform(meMomenta()[trans[k]]);
	}

	meMomenta()[bornEmitter()] =
	const_cast<const TildeKinematics&>(*kinematics).bornEmitterMomentum();
	meMomenta()[bornSpectator()] =
	const_cast<const TildeKinematics&>(*kinematics).bornSpectatorMomentum();

	cPDVector::const_iterator pd = mePartonData().begin();
	vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
	for ( ; pd != mePartonData().end(); ++pd, ++p ) {
	p->setMass((**pd).hardProcessMass());
	p->rescaleRho();
	}

	jacobian(realEmissionME()->lastXComb().jacobian());

	logGenerateTildeKinematics();

	return true;

	}

	void SubtractionDipole::invertedTildeKinematics(Ptr<InvertedTildeKinematics>::tptr itk) {
	theInvertedTildeKinematics = itk;
	}

	int SubtractionDipole::nDimRadiation() const {
	return invertedTildeKinematics() ?
	invertedTildeKinematics()->nDimRadiation() :
	0;
	}

	bool SubtractionDipole::generateRadiationKinematics(const double * r) {

	assert(splitting());

	Ptr<InvertedTildeKinematics>::tptr kinematics = theInvertedTildeKinematics;
	if ( showerApproximation() ) {
	showerApproximation()->setBornXComb(lastHeadXCombPtr());
	showerApproximation()->setRealXComb(lastXCombPtr());
	showerApproximation()->setDipole(this);
	if ( showerApproximation()->showerInvertedTildeKinematics() ) {
	kinematics = showerApproximation()->showerInvertedTildeKinematics();
	}
	}

	if ( !kinematics ) {
	jacobian(0.0);
	return false;
	}

	kinematics->prepare(lastXCombPtr(),lastHeadXCombPtr());

	if ( !kinematics->doMap(r) ) {
	jacobian(0.0);
	return false;
	}

	theLastSplittingScale = kinematics->lastScale();
	theLastSplittingPt = kinematics->lastPt();
	theLastSplittingZ = kinematics->lastZ();

	meMomenta().resize(lastHeadXComb().meMomenta().size() + 1);

	assert(splittingMap().find(lastUnderlyingBornKey) != splittingMap().end());
	map<int,int>& trans = const_cast<map<int,int>&>(lastRealEmissionInfo->second.second);

	int n = lastHeadXComb().meMomenta().size();
	for ( int k = 0; k < n; ++k ) {
	if ( k == bornEmitter() \|\| k == bornSpectator() )
	continue;
	meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
	if ( kinematics->doesTransform() && k > 1 )
	meMomenta()[trans[k]] = kinematics->transform(meMomenta()[trans[k]]);
	}

	meMomenta()[realEmitter()] =
	const_cast<const InvertedTildeKinematics&>(*kinematics).realEmitterMomentum();
	meMomenta()[realEmission()] =
	const_cast<const InvertedTildeKinematics&>(*kinematics).realEmissionMomentum();
	meMomenta()[realSpectator()] =
	const_cast<const InvertedTildeKinematics&>(*kinematics).realSpectatorMomentum();

	cPDVector::const_iterator pd = mePartonData().begin();
	vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
	for ( ; pd != mePartonData().end(); ++pd, ++p ) {
	p->setMass((**pd).hardProcessMass());
	p->rescaleRho();
	}

	jacobian(underlyingBornME()->lastXComb().jacobian() *
	kinematics->jacobian());

	logGenerateRadiationKinematics(r);

	return true;

	}

	void SubtractionDipole::ptCut(Energy cut) {
	theInvertedTildeKinematics->ptCut(cut);
	}

	CrossSection SubtractionDipole::dSigHatDR(Energy2 factorizationScale) const {

	double pdfweight = 1.;

	double jac = jacobian();

	if ( splitting() && jac == 0.0 ) {
	lastMECrossSection(ZERO);
	return ZERO;
	}

	if ( factorizationScale == ZERO ) {
	factorizationScale = underlyingBornME()->lastScale();
	}

	if ( havePDFWeight1() ) {
	pdfweight *= realEmissionME()->pdf1(factorizationScale);
	}

	if ( havePDFWeight2() ) {
	pdfweight *= realEmissionME()->pdf2(factorizationScale);
	}

	lastMEPDFWeight(pdfweight);

	bool needTheDipole = true;
	CrossSection shower = ZERO;

	double lastThetaMu = 1.0;

	double showerFactor = 1.;

	if ( showerApproximation() ) {
	assert(!splitting());
	showerApproximation()->setBornXComb(lastXCombPtr());
	showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
	showerApproximation()->setDipole(const_cast<SubtractionDipole*>(this));
	if ( !isAboveCutoff() ) {
	showerApproximation()->wasBelowCutoff();
	lastThetaMu = 0.0;
	} else {
	lastThetaMu = 1.0;
	}
	if ( lastThetaMu > 0.0 && isInShowerPhasespace() ) {
	if ( realShowerSubtraction() )
	shower = showerApproximation()->dSigHatDR()*lastThetaMu;
	if ( virtualShowerSubtraction() \|\| loopSimSubtraction() )
	shower = -showerApproximation()->dSigHatDR()*lastThetaMu;
	if ( virtualShowerSubtraction() &&
	isAboveCutoff() &&
	showerApproximation()->showerTildeKinematics() ) {
	// map shower to dipole kinematics; we are always above the
	// cutoff in this case
	showerFactor *=
	showerApproximation()->showerTildeKinematics()->jacobianRatio();
	}
	shower *= showerFactor;
	}
	if ( realShowerSubtraction() && lastThetaMu == 1.0 )
	needTheDipole = false;
	if ( virtualShowerSubtraction() && lastThetaMu == 0.0 )
	needTheDipole = false;
	if ( factory()->loopSimCorrections() \|\|
	factory()->meCorrectionsOnly() )
	needTheDipole = false;
	}

	double xme2 = 0.0;

	if ( needTheDipole )
	xme2 = me2();

	if ( factory()->loopSimCorrections() \|\|
	factory()->meCorrectionsOnly() ) {

	assert(showerApproximation());
	xme2 = realEmissionME()->me2() * showerApproximation()->channelWeight();

	double rws =
	pow(underlyingBornME()->lastXComb().lastAlphaS()/
	realEmissionME()->lastXComb().lastAlphaS(),
	realEmissionME()->orderInAlphaS());

	xme2 *= rws;

	double rwe =
	pow(underlyingBornME()->lastXComb().lastAlphaEM()/
	realEmissionME()->lastXComb().lastAlphaEM(),
	underlyingBornME()->orderInAlphaEW());

	xme2 *= rwe;

	}

	if ( realShowerSubtraction() )
	xme2 *= 1. - lastThetaMu;
	if ( virtualShowerSubtraction() \|\| loopSimSubtraction() )
	xme2 *= lastThetaMu;

	double coupl = lastMECouplings();
	coupl *= underlyingBornME()->lastXComb().lastAlphaS();
	lastMECouplings(coupl);

	CrossSection res =
	sqr(hbarc) * jac * pdfweight * xme2 /
	(2. * realEmissionME()->lastXComb().lastSHat());

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

	lastMECrossSection(-res-shower);

	logDSigHatDR(jac);

	return lastMECrossSection();
	}


	bool SubtractionDipole::aboveAlpha() const{return theTildeKinematics->aboveAlpha();}



	-
	-
	-CrossSection SubtractionDipole::ps(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const {
	- double jac = jacobian();
	+CrossSection SubtractionDipole::prefactor(Energy2 factorizationScale)const{
	+
	+ const double jac = jacobian();
	assert( factorizationScale != ZERO );
	assert (! splitting());
	if(!theRealEmissionME->clustersafe(realEmitter(),realEmission(),realSpectator()).second \|\| jac == 0.0 ) {
	lastMECrossSection(ZERO);
	lastME2(0.0);
	return ZERO;
	}
	double pdfweight = 1.;
	if ( havePDFWeight1() ) pdfweight *= realEmissionME()->pdf1(factorizationScale);
	if ( havePDFWeight2() ) pdfweight *= realEmissionME()->pdf2(factorizationScale);
	+

	- double ccme2 =underlyingBornME()->me2()*
	- underlyingBornME()->largeNColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),largeNBasis)
	- /underlyingBornME()->largeNME2(largeNBasis);
	+ return sqr(hbarc) * jac * pdfweight / (2. * realEmissionME()->lastXComb().lastSHat());

	- //double ccme2 =underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
	-
	- return sqr(hbarc) * jac * pdfweight * (me2Avg(ccme2)) / (2. * realEmissionME()->lastXComb().lastSHat());
	}




	+CrossSection SubtractionDipole::ps(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const {
	+
	+ double ccme2 =underlyingBornME()->me2()*
	+ underlyingBornME()->
	+ largeNColourCorrelatedME2(
	+ make_pair(bornEmitter(),bornSpectator()),largeNBasis)/
	+ underlyingBornME()->largeNME2(largeNBasis);
	+
	+ return prefactor(factorizationScale) * me2Avg(ccme2);
	+}
	+
	+
	+
	+
	pair<CrossSection,CrossSection> SubtractionDipole::dipandPs(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const {
	- double jac = jacobian();
	- assert( factorizationScale != ZERO );
	- assert (! splitting());
	- if(!theRealEmissionME->clustersafe(realEmitter(),realEmission(),realSpectator()).second \|\| jac == 0.0 ) {
	- lastMECrossSection(ZERO);
	- lastME2(0.0);
	- return make_pair(ZERO,ZERO);
	- }

	- double pdfweight = 1.;
	- if ( havePDFWeight1() ) pdfweight *= realEmissionME()->pdf1(factorizationScale);
	- if ( havePDFWeight2() ) pdfweight *= realEmissionME()->pdf2(factorizationScale);
	-
	+ CrossSection factor= prefactor(factorizationScale);
	+
	double ccme2 =underlyingBornME()->me2()*
	underlyingBornME()->
	largeNColourCorrelatedME2(
	make_pair(bornEmitter(),bornSpectator()),largeNBasis)/
	underlyingBornME()->largeNME2(largeNBasis);

	double ps = me2Avg(ccme2);
	double dip = me2();

	- CrossSection factor= sqr(hbarc) * jac * pdfweight / (2. * realEmissionME()->lastXComb().lastSHat());
	-
	-
	return make_pair(factordip,factorps);
	}

	CrossSection SubtractionDipole::dip(Energy2 factorizationScale) const {
	- double jac = jacobian();
	- assert( factorizationScale != ZERO );
	- assert (! splitting());
	- if(!theRealEmissionME->clustersafe(realEmitter(),realEmission(),realSpectator()).second \|\| jac == 0.0 ) {
	- lastMECrossSection(ZERO);
	- lastME2(0.0);
	- return ZERO;
	- }
	-
	- double pdfweight = 1.;
	- if ( havePDFWeight1() ) pdfweight *= realEmissionME()->pdf1(factorizationScale);
	- if ( havePDFWeight2() ) pdfweight *= realEmissionME()->pdf2(factorizationScale);
	-
	+ CrossSection factor= prefactor(factorizationScale);
	double dip = me2();
	-
	- CrossSection factor= sqr(hbarc) * jac * pdfweight / (2. * realEmissionME()->lastXComb().lastSHat());
	-
	-
	return factor*dip;
	}





	bool SubtractionDipole::clustersafe()const {

	return (theRealEmissionME->clustersafe(realEmitter(),realEmission(),realSpectator()).second);

	}


	-bool SubtractionDipole::clustersafe(){
	-
	- return (theRealEmissionME->clustersafe(realEmitter(),realEmission(),realSpectator()).second);
	-
	-}


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

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

	os << " subtraction '" << name() << "'\n for real emission '"
	<< theRealEmissionME->name() << "'\n using underlying Born '"
	<< theUnderlyingBornME->name() << "'\n";

	os << " tilde kinematics are '"
	<< (theTildeKinematics ? theTildeKinematics->name() : "")
	<< " '\n inverted tilde kinematics are '"
	<< (theInvertedTildeKinematics ? theInvertedTildeKinematics->name() : "") << "'\n";

	os << " the following subtraction mappings have been found:\n";

	for ( map<RealEmissionKey,UnderlyingBornInfo>::const_iterator m =
	theMergingMap.begin(); m != theMergingMap.end(); ++m ) {
	os << " " << process(m->second.first)[0]->PDGName() << " "
	<< process(m->second.first)[1]->PDGName() << " -> ";
	for ( cPDVector::const_iterator p = process(m->second.first).begin() + 2;
	p != process(m->second.first).end(); ++p ) {
	os << (**p).PDGName() << " ";
	}
	os << "[" << emitter(m->second.first) << "," << spectator(m->second.first) << "] <=> ";
	os << process(m->first)[0]->PDGName() << " "
	<< process(m->first)[1]->PDGName() << " -> ";
	for ( cPDVector::const_iterator p = process(m->first).begin() + 2;
	p != process(m->first).end(); ++p ) {
	os << (**p).PDGName() << " ";
	}
	os << "[(" << emitter(m->first) << "," << emission(m->first) << ")," << spectator(m->first) << "]\n"
	<< " non-dipole momenta ( ";
	for ( map<int,int>::const_iterator k = m->second.second.begin();
	k != m->second.second.end(); ++k ) {
	if ( k->first == spectator(m->first) )
	continue;
	os << k->second << " ";
	}
	os << ") <=> ( ";
	for ( map<int,int>::const_iterator k = m->second.second.begin();
	k != m->second.second.end(); ++k ) {
	if ( k->first == spectator(m->first) )
	continue;
	os << k->first << " ";
	}
	os << ")\n";
	}

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

	os << flush;

	}

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

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

	os << " for dipole '" << name() << "' applying ["
	<< bornEmitter() << "," << bornSpectator() << "] <=> [("
	<< realEmitter() << "," << realEmission() << ")," << realSpectator() << "]\n"
	<< " evaluated the cross section/nb " << (lastMECrossSection()/nanobarn) << "\n"
	<< " with subtraction parameters x[0] = " << subtractionParameters()[0]
	<< " x[1] = " << subtractionParameters()[1] << "\n";

	os << " the last real emission event was:\n";
	realEmissionME()->printLastEvent(os);

	os << " the last underlying Born event was:\n";
	underlyingBornME()->printLastEvent(os);


	os << "--- end SubtractionDipole last event information -------------------------------\n";

	os << flush;

	}

	void SubtractionDipole::logME2() const {

	if ( !realEmissionME()->verbose() &&
	!underlyingBornME()->verbose() )
	return;

	tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
	tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();

	generator()->log() << "'" << name() << "' evaluated me2 using\n"
	<< "Born XComb " << bornxc << " real XComb " << realxc << "\n";

	generator()->log() << "subtraction parameters: ";
	copy(subtractionParameters().begin(),subtractionParameters().end(),
	ostream_iterator<double>(generator()->log()," "));
	generator()->log() << "\n";

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

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

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

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

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

	pit = realxc->meMomenta().begin();
	dit = realxc->mePartonData().begin();

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

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

	}

	void SubtractionDipole::logDSigHatDR(double effectiveJac) const {

	if ( !realEmissionME()->verbose() &&
	!underlyingBornME()->verbose() )
	return;

	tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
	tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();

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

	}

	void SubtractionDipole::logGenerateTildeKinematics() const {

	if ( !realEmissionME()->verbose() &&
	!underlyingBornME()->verbose() )
	return;

	generator()->log() << "'" << name() << "' generating tilde kinematics.\n"
	<< "configuration: [" << bornEmitter() << ","
	<< bornSpectator() << "] => "
	<< "[(" << realEmitter() << "," << realEmission() << "),"
	<< realSpectator() << "]\n"
	<< "with real xcomb " << lastHeadXCombPtr() << " born xcomb "
	<< lastXCombPtr() << "\n"
	<< "from real emission phase space point:\n";
	Lorentz5Momentum rSum;
	vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
	cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
	size_t count = 0;
	for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr ) {
	generator()->log() << (**dr).PDGName() << " : "
	<< (*pr/GeV) << "\n";
	if ( count < 2 ) {
	rSum -= *pr;
	} else {
	rSum += *pr;
	}
	++count;

	}
	generator()->log() << "sum : " << (rSum/GeV) << "\n";

	generator()->log() << "subtraction parameters: ";
	copy(subtractionParameters().begin(),subtractionParameters().end(),
	ostream_iterator<double>(generator()->log()," "));
	generator()->log() << "\n"
	<< "with scale/GeV = " << (theLastSubtractionScale/GeV)
	<< "and pt/GeV = " << (theLastSubtractionPt/GeV) << "\n";

	generator()->log() << "generated tilde kinematics:\n";
	pr = lastXComb().meMomenta().begin();
	dr = lastXComb().mePartonData().begin();
	count = 0;
	Lorentz5Momentum bSum;
	for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr ) {
	generator()->log() << (**dr).PDGName() << " : "
	<< (*pr/GeV) << "\n";
	if ( count < 2 ) {
	bSum -= *pr;
	} else {
	bSum += *pr;
	}
	++count;
	}
	generator()->log() << "sum : " << (bSum/GeV) << "\n";

	generator()->log() << "Jacobian = " << jacobian() << "\n" << flush;

	}


	void SubtractionDipole::logGenerateRadiationKinematics(const double * r) const {

	if ( !realEmissionME()->verbose() &&
	!underlyingBornME()->verbose() )
	return;

	generator()->log() << "'" << name() << "' generating radiation kinematics.\n"
	<< "configuration: [" << bornEmitter() << ","
	<< bornSpectator() << "] => "
	<< "[(" << realEmitter() << "," << realEmission() << "),"
	<< realSpectator() << "]\n"
	<< "with born xcomb " << lastHeadXCombPtr() << " real xcomb "
	<< lastXCombPtr() << "\n"
	<< "from random numbers:\n";
	copy(r,r+nDimRadiation(),ostream_iterator<double>(generator()->log()," "));
	generator()->log() << "\n";
	generator()->log() << "and born phase space point:\n";
	vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
	cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
	for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr )
	generator()->log() << (**dr).PDGName() << " : "
	<< (*pr/GeV) << "\n";

	generator()->log() << "subtraction parameters: ";
	copy(subtractionParameters().begin(),subtractionParameters().end(),
	ostream_iterator<double>(generator()->log()," "));
	generator()->log() << "\n" << flush;

	generator()->log() << "scales: scale/GeV = " << (theLastSplittingScale/GeV)
	<< " pt/GeV = " << (theLastSplittingPt/GeV) << "\n" << flush;

	generator()->log() << "generated real emission kinematics:\n";
	pr = lastXComb().meMomenta().begin();
	dr = lastXComb().mePartonData().begin();
	for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr )
	generator()->log() << (**dr).PDGName() << " : "
	<< (*pr/GeV) << "\n";

	generator()->log() << "Jacobian = "
	<< jacobian() << " = "
	<< underlyingBornME()->lastXComb().jacobian()
	<< "\|Born * "
	<< invertedTildeKinematics()->jacobian()
	<< "\|Radiation\n" << flush;

	}

	void SubtractionDipole::doinit() {
	MEBase::doinit();
	if ( underlyingBornME() ) {
	theUnderlyingBornME->init();
	}
	if ( realEmissionME() ) {
	theRealEmissionME->init();
	}
	if ( tildeKinematics() ) {
	theTildeKinematics->init();
	}
	if ( invertedTildeKinematics() ) {
	theInvertedTildeKinematics->init();
	}
	if ( showerApproximation() ) {
	theShowerApproximation->init();
	}
	for ( vector<Ptr<SubtractionDipole>::tptr>::iterator p = thePartners.begin();
	p != thePartners.end(); ++p ) {
	(**p).init();
	}
	for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
	theReweights.begin(); rw != theReweights.end(); ++rw ) {
	(**rw).init();
	}
	}

	void SubtractionDipole::doinitrun() {
	MEBase::doinitrun();
	if ( underlyingBornME() ) {
	theUnderlyingBornME->initrun();
	}
	if ( realEmissionME() ) {
	theRealEmissionME->initrun();
	}
	if ( tildeKinematics() ) {
	theTildeKinematics->initrun();
	}
	if ( invertedTildeKinematics() ) {
	theInvertedTildeKinematics->initrun();
	}
	if ( showerApproximation() ) {
	theShowerApproximation->initrun();
	}
	for ( vector<Ptr<SubtractionDipole>::tptr>::iterator p = thePartners.begin();
	p != thePartners.end(); ++p ) {
	(**p).initrun();
	}
	for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
	theReweights.begin(); rw != theReweights.end(); ++rw ) {
	(**rw).initrun();
	}
	}

	void SubtractionDipole::cloneDependencies(const std::string& prefix,bool slim) {

	if ( underlyingBornME() ) {
	Ptr<MatchboxMEBase>::ptr myUnderlyingBornME = underlyingBornME()->cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << myUnderlyingBornME->name();
	if ( ! (generator()->preinitRegister(myUnderlyingBornME,pname.str()) ) )
	throw Exception() << "SubtractionDipole::cloneDependencies(): Matrix element " << pname.str() << " already existing." << Exception::runerror;
	myUnderlyingBornME->cloneDependencies(pname.str(),slim);
	underlyingBornME(myUnderlyingBornME);
	}

	if ( realEmissionME()&& !slim ) {
	Ptr<MatchboxMEBase>::ptr myRealEmissionME = realEmissionME()->cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << myRealEmissionME->name();
	if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
	throw Exception() << "SubtractionDipole::cloneDependencies(): Matrix element " << pname.str() << " already existing." << Exception::runerror;
	myRealEmissionME->cloneDependencies(pname.str());
	realEmissionME(myRealEmissionME);
	}

	if ( tildeKinematics() ) {
	Ptr<TildeKinematics>::ptr myTildeKinematics = tildeKinematics()->cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << myTildeKinematics->name();
	if ( ! (generator()->preinitRegister(myTildeKinematics,pname.str()) ) )
	throw Exception() << "SubtractionDipole::cloneDependencies(): Tilde kinematics " << pname.str() << " already existing." << Exception::runerror;
	myTildeKinematics->dipole(this);
	tildeKinematics(myTildeKinematics);
	}

	if ( invertedTildeKinematics()&& !slim ) {
	Ptr<InvertedTildeKinematics>::ptr myInvertedTildeKinematics = invertedTildeKinematics()->cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << myInvertedTildeKinematics->name();
	if ( ! (generator()->preinitRegister(myInvertedTildeKinematics,pname.str()) ) )
	throw Exception() << "SubtractionDipole::cloneDependencies(): Inverted tilde kinematics " << pname.str() << " already existing." << Exception::runerror;
	myInvertedTildeKinematics->dipole(this);
	invertedTildeKinematics(myInvertedTildeKinematics);
	}

	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 Exception() << "SubtractionDipole::cloneDependencies(): Reweight " << pname.str() << " already existing." << Exception::runerror;
	myReweight->cloneDependencies(pname.str());
	*rw = myReweight;
	}

	}

	void SubtractionDipole::constructVertex(tSubProPtr sub) {
	if ( splitting() )
	realEmissionME()->constructVertex(sub);
	else
	underlyingBornME()->constructVertex(sub);
	}

	void SubtractionDipole::constructVertex(tSubProPtr sub, const ColourLines* cl) {
	if ( splitting() )
	realEmissionME()->constructVertex(sub,cl);
	else
	underlyingBornME()->constructVertex(sub,cl);
	}

	void SubtractionDipole::generateSubCollision(SubProcess & sub) {
	if ( splitting() )
	realEmissionME()->generateSubCollision(sub);
	else
	underlyingBornME()->generateSubCollision(sub);
	}

	void SubtractionDipole::persistentOutput(PersistentOStream & os) const {
	os << theLastXComb << theSplitting << theApply << theSubtractionTest
	<< theIgnoreCuts << theRealEmissionME << theUnderlyingBornME
	<< thePartners << theTildeKinematics << theInvertedTildeKinematics
	<< theReweights << theRealEmitter << theRealEmission << theRealSpectator
	<< theSubtractionParameters << theMergingMap << theSplittingMap
	<< theIndexMap << theUnderlyingBornDiagrams << theRealEmissionDiagrams
	<< theBornToRealDiagrams << theRealToBornDiagrams
	<< lastRealEmissionKey << lastUnderlyingBornKey
	<< theBornEmitter << theBornSpectator << ounit(theLastSubtractionScale,GeV)
	<< ounit(theLastSplittingScale,GeV) << ounit(theLastSubtractionPt,GeV)
	<< ounit(theLastSplittingPt,GeV) << theLastSubtractionZ
	<< theLastSplittingZ << theShowerApproximation
	<< theRealShowerSubtraction << theVirtualShowerSubtraction
	<< theLoopSimSubtraction << theRealEmissionScales << theFactory
	<< ounit(theShowerHardScale,GeV) << ounit(theShowerScale,GeV)
	<< theShowerParameters << theIsInShowerPhasespace << theIsAboveCutoff;
	}

	void SubtractionDipole::persistentInput(PersistentIStream & is, int) {
	is >> theLastXComb >> theSplitting >> theApply >> theSubtractionTest
	>> theIgnoreCuts >> theRealEmissionME >> theUnderlyingBornME
	>> thePartners >> theTildeKinematics >> theInvertedTildeKinematics
	>> theReweights >> theRealEmitter >> theRealEmission >> theRealSpectator
	>> theSubtractionParameters >> theMergingMap >> theSplittingMap
	>> theIndexMap >> theUnderlyingBornDiagrams >> theRealEmissionDiagrams
	>> theBornToRealDiagrams >> theRealToBornDiagrams
	>> lastRealEmissionKey >> lastUnderlyingBornKey
	>> theBornEmitter >> theBornSpectator >> iunit(theLastSubtractionScale,GeV)
	>> iunit(theLastSplittingScale,GeV) >> iunit(theLastSubtractionPt,GeV)
	>> iunit(theLastSplittingPt,GeV) >> theLastSubtractionZ
	>> theLastSplittingZ >> theShowerApproximation
	>> theRealShowerSubtraction >> theVirtualShowerSubtraction
	>> theLoopSimSubtraction >> theRealEmissionScales >> theFactory
	>> iunit(theShowerHardScale,GeV) >> iunit(theShowerScale,GeV)
	>> theShowerParameters >> theIsInShowerPhasespace >> theIsAboveCutoff;
	lastMatchboxXComb(theLastXComb);
	typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
	pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
	lastRealEmissionInfo = kr.first;
	for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
	if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
	break;
	}

	Ptr<MatchboxFactory>::tptr SubtractionDipole::factory() const {
	return theFactory;
	}

	void SubtractionDipole::factory(Ptr<MatchboxFactory>::tptr f) {
	theFactory = f;
	}


	void SubtractionDipole::Init() {

	static ClassDocumentation<SubtractionDipole> documentation
	("SubtractionDipole represents a dipole subtraction "
	"term in the formalism of Catani and Seymour.");

	}

	// * 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).
	DescribeAbstractClass<SubtractionDipole,MEBase>
	describeSubtractionDipole("Herwig::SubtractionDipole", "Herwig.so");
	diff --git a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h
	--- a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h
	+++ b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h
	@@ -1,1269 +1,1271 @@
	// -- C++ --
	//
	// SubtractionDipole.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_SubtractionDipole_H
	#define HERWIG_SubtractionDipole_H
	//
	// This is the declaration of the SubtractionDipole class.
	//

	#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.fh"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.fh"

	#include "ThePEG/MatrixElement/MEBase.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
	#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief SubtractionDipole represents a dipole subtraction
	* term in the formalism of Catani and Seymour.
	*
	*/
	class SubtractionDipole:
	public MEBase, public LastMatchboxXCombInfo {

	public:

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

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

	public:

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

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

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

	/**
	* A helpre struct to communicate diagram merging and remapping
	* information
	*/
	struct MergeInfo {

	/**
	* The merged emitter
	*/
	int emitter;

	/**
	* The Born diagram
	*/
	Ptr<Tree2toNDiagram>::ptr diagram;

	/**
	* The merging map
	*/
	map<int,int> mergeLegs;

	};

	/**
	* Return true, if this dipole can possibly handle the indicated
	* emitter.
	*/
	virtual bool canHandleEmitter(const cPDVector& partons, int emitter) const = 0;

	/**
	* Return true, if this dipole can possibly handle the indicated
	* splitting.
	*/
	virtual bool canHandleSplitting(const cPDVector& partons, int emitter, int emission) const = 0;

	/**
	* Return true, if this dipole can possibly handle the indicated
	* spectator.
	*/
	virtual bool canHandleSpectator(const cPDVector& partons, int spectator) const = 0;

	/**
	* Return true, if this dipole applies to the selected
	* configuration.
	*/
	virtual bool canHandle(const cPDVector& partons,
	int emitter, int emission, int spectator) const = 0;

	/**
	* Return true, if this dipole is symmetric with respect to emitter
	* and emission.
	*/
	virtual bool isSymmetric() const { return false; }

	/**
	* If this is a dependent matrix element in a ME group, return true,
	* if it applies to the process set in lastXComb()
	*/
	virtual bool apply() const { return theApply; }

	/**
	* Clear the bookkeeping
	*/
	void clearBookkeeping();

	/**
	* Setup bookkeeping maps.
	*/
	void setupBookkeeping(const map<Ptr<DiagramBase>::ptr,MergeInfo>& mergeInfo,bool slim);

	/**
	* Get bookkeeping information for the given
	* real emission diagram
	*/
	void subtractionBookkeeping();

	/**
	* Determine bookkeeping information for
	* the underlying Born process supplied through
	* the lastHeadXComb() object.
	*/
	void splittingBookkeeping();

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

	/**
	* Create a dependent xcomb object for the underlying
	* Born process, given a XComb driving the real emission
	*/
	StdXCombPtr makeBornXComb(tStdXCombPtr realXC);

	/**
	* Create dependent xcomb objects for the real emission process,
	* given a XComb driving the underlying Born
	*/
	vector<StdXCombPtr> makeRealXCombs(tStdXCombPtr bornXC);

	/**
	* Return true, if bookkeeping did not find a non-trivial setup.
	*/
	bool empty() const { return theSplittingMap.empty()&&theMergingMap.empty(); }

	/**
	* Return the emitter as referred to by the real emission
	* matrix element.
	*/
	int realEmitter() const { return theRealEmitter; }

	/**
	* Set the emitter as referred to by the real emission
	* matrix element.
	*/
	void realEmitter(int id) { theRealEmitter = id; }

	/**
	* Return the emission as referred to by the real emission
	* matrix element.
	*/
	int realEmission() const { return theRealEmission; }

	/**
	* Set the emission as referred to by the real emission
	* matrix element.
	*/
	void realEmission(int id) { theRealEmission = id; }

	/**
	* Return the spectator as referred to by the real emission
	* matrix element.
	*/
	int realSpectator() const { return theRealSpectator; }

	/**
	* Set the spectator as referred to by the real emission
	* matrix element.
	*/
	void realSpectator(int id) { theRealSpectator = id; }

	/**
	* Return the emitter as referred to by the underlying
	* Born process.
	*/
	int bornEmitter() const { return theBornEmitter; }

	/**
	* Set the emitter as referred to by the underlying
	* Born process.
	*/
	void bornEmitter(int id) { theBornEmitter = id; }

	/**
	* Return the spectator as referred to by the underlying
	* Born process.
	*/
	int bornSpectator() const { return theBornSpectator; }

	/**
	* Set the spectator as referred to by the underlying
	* Born process.
	*/
	void bornSpectator(int id) { theBornSpectator = id; }

	/**
	* Define the real emission key type
	*/
	typedef pair<pair<cPDVector,int>,pair<int,int> > RealEmissionKey;

	/**
	* Create a real emission key
	*/
	static RealEmissionKey realEmissionKey(const cPDVector& proc,
	int em, int emm, int sp) {
	return make_pair(make_pair(proc,emm),make_pair(em,sp));
	}

	/**
	* Return the diagram of a real emission key
	*/
	static const cPDVector& process(const RealEmissionKey& key) {
	return key.first.first;
	}

	/**
	* Return the emission id of a real emission key
	*/
	static int emission(const RealEmissionKey& key) {
	return key.first.second;
	}

	/**
	* Return the emitter id of a real emission key
	*/
	static int emitter(const RealEmissionKey& key) {
	return key.second.first;
	}

	/**
	* Return the spectator id of a real emission key
	*/
	static int spectator(const RealEmissionKey& key) {
	return key.second.second;
	}

	/**
	* Define the underlying Born key type
	*/
	typedef pair<cPDVector,pair<int,int> > UnderlyingBornKey;

	/**
	* Create a underlying Born key
	*/
	static UnderlyingBornKey underlyingBornKey(const cPDVector& proc,
	int em, int sp) {
	return make_pair(proc,make_pair(em,sp));
	}

	/**
	* Return the diagram of a underlying Born key
	*/
	static const cPDVector& process(const UnderlyingBornKey& key) {
	return key.first;
	}

	/**
	* Return the emitter id of a underlying Born key
	*/
	static int emitter(const UnderlyingBornKey& key) {
	return key.second.first;
	}

	/**
	* Return the spectator id of a underlying Born key
	*/
	static int spectator(const UnderlyingBornKey& key) {
	return key.second.second;
	}

	/**
	* Define real emission key and index dictionary
	* for partons not involved in the given dipole.
	*/
	typedef pair<RealEmissionKey,map<int,int> > RealEmissionInfo;

	/**
	* Define underlying Born key and index dictionary
	* for partons not involved in the given dipole.
	*/
	typedef pair<UnderlyingBornKey,map<int,int> > UnderlyingBornInfo;

	/**
	* Return the merging map
	*/
	const map<RealEmissionKey,UnderlyingBornInfo>& mergingMap() const { return theMergingMap; }

	/**
	* Return the splitting map
	*/
	const multimap<UnderlyingBornKey,RealEmissionInfo>& splittingMap() const { return theSplittingMap; }

	/**
	* Return the underlying Born diagrams to be considered
	* for the given real emission process.
	*/
	const DiagramVector& underlyingBornDiagrams(const cPDVector& real) const;

	/**
	* Find the underlying Born diagram for the given real emission diagram
	*/
	tcDiagPtr underlyingBornDiagram(tcDiagPtr realDiag) const;

	/**
	* Return the real emission diagrams to be considered
	* for the given Born process.
	*/
	const DiagramVector& realEmissionDiagrams(const cPDVector& born) const;

	/**
	* Find the real emission diagram for the given underlying Born diagram
	*/
	tcDiagPtr realEmissionDiagram(tcDiagPtr bornDiag) const;

	/**
	* Add all possible diagrams with the add() function.
	*/
	virtual void getDiagrams() 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; }

	/**
	* 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 & dv) const;

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

	/**
	* Select a ColpurLines geometry. The default version returns a
	* colour geometry selected among the ones returned from
	* colourGeometries(tcDiagPtr).
	*/
	virtual const ColourLines &
	selectColourGeometry(tcDiagPtr diag) const;

	/**
	* Return the order in \f$\alpha_S\f$ in which this matrix element
	* is given.
	*/
	virtual unsigned int orderInAlphaS() const { return realEmissionME()->orderInAlphaS(); }

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

	//@}

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

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

	/**
	* Set the typed and momenta of the incoming and outgoing partons to
	* be used in subsequent calls to me() and colourGeometries()
	* according to the associated XComb object. If the function is
	* overridden in a sub class the new function must call the base
	* class one first.
	*/
	virtual void setKinematics();

	/**
	* 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.
	*/
	virtual bool generateKinematics(const double * r);

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

	/**
	* Return true, if this matrix element expects
	* the incoming partons in their center-of-mass system
	*/
	virtual bool wantCMS () const { return realEmissionME()->wantCMS(); }

	/**
	* Clear the information previously provided by a call to
	* setKinematics(...).
	*/
	virtual void clearKinematics();

	/**
	* If this is a dependent matrix element in a ME group, return true,
	* if cuts should be ignored.
	*/
	virtual bool ignoreCuts() const { return theIgnoreCuts; }

	/**
	* Indicate that cuts should be ignored
	*/
	void doIgnoreCuts(bool is = true) { theIgnoreCuts = is; }

	//@}

	/** @name Tilde kinematics */
	//@{

	/**
	* Return the TildeKinematics object used
	*/
	Ptr<TildeKinematics>::tcptr tildeKinematics() const { return theTildeKinematics; }

	/**
	* Set the TildeKinematics object used
	*/
	void tildeKinematics(Ptr<TildeKinematics>::tptr);

	/**
	* Generate the tilde kinematics from real emission
	* kinematics accessible through the XComb's
	* head object and store it in meMomenta(). This default
	* implemenation uses the tildeKinematics() object.
	*/
	virtual bool generateTildeKinematics();

	/**
	* Return the InvertedTildeKinematics object used
	*/
	Ptr<InvertedTildeKinematics>::tcptr invertedTildeKinematics() const { return theInvertedTildeKinematics; }

	/**
	* Set the InvertedTildeKinematics object used
	*/
	void invertedTildeKinematics(Ptr<InvertedTildeKinematics>::tptr);

	/**
	* Return the number of additional random numbers
	* needed to generate real emission kinematics off
	* the tilde kinematics previously supplied through
	* the XComb object. This default implementation
	* returns invertedTildeKinematics()->nDimRadiation()
	*/
	virtual int nDimRadiation() const;

	/**
	* Generate the real emission kinematics
	* off the Born kinematics accessible through the XComb's
	* head object and store it in meMomenta(); store
	* the single particle phasespace in units of lastHeadXComb()->lastSHat()
	* in jacobian(). This default
	* implemenation uses the invertedTildeKinematics() object
	*/
	virtual bool generateRadiationKinematics(const double *);

	/**
	* Set a pt cut when splitting
	*/
	void ptCut(Energy cut);

	/**
	* Return the relevant dipole scale
	*/
	Energy lastDipoleScale() const {
	return splitting() ? theLastSplittingScale : theLastSubtractionScale;
	}

	/**
	* Return the relevant pt
	*/
	Energy lastPt() const {
	return splitting() ? theLastSplittingPt : theLastSubtractionPt;
	}

	/**
	* Return the relevant momentum fractions
	*/
	double lastZ() const {
	return splitting() ? theLastSplittingZ : theLastSubtractionZ;
	}

	/**
	* Return true, if this dipole acts in splitting mode.
	*/
	bool splitting() const { return theSplitting; }

	/**
	* Switch on splitting mode for this dipole.
	*/
	void doSplitting() { theSplitting = true; }

	/**
	* Switch off splitting mode for this dipole.
	*/
	void doSubtraction() { theSplitting = false; }

	/**
	* Return the subtraction parameters.
	*/
	const vector<double>& subtractionParameters() const { return theSubtractionParameters; }

	/**
	* Access the subtraction parameters.
	*/
	vector<double>& subtractionParameters() { return theSubtractionParameters; }

	/**
	* Return the shower hard scale encountered
	*/
	Energy showerHardScale() const { return theShowerHardScale; }

	/**
	* Set the shower hard scale encountered
	*/
	void showerHardScale(Energy s) { theShowerHardScale = s; }

	/**
	* Return the shower evolution scale encountered
	*/
	Energy showerScale() const { return theShowerScale; }

	/**
	* Set the shower evolution scale encountered
	*/
	void showerScale(Energy s) { theShowerScale = s; }

	/**
	* Return the shower splitting variables encountered
	*/
	const vector<double>& showerParameters() const { return theShowerParameters; }

	/**
	* Access the shower splitting variables encountered
	*/
	vector<double>& showerParameters() { return theShowerParameters; }

	/**
	* Return true, if this configuration is in the shower phase space
	*/
	bool isInShowerPhasespace() const { return theIsInShowerPhasespace; }

	/**
	* Indicate whether this configuration is in the shower phase space
	*/
	void isInShowerPhasespace(bool yes) { theIsInShowerPhasespace = yes; }

	/**
	* Return true, if this configuration is above the shower infrared cutoff
	*/
	bool isAboveCutoff() const { return theIsAboveCutoff; }

	/**
	* Indicate whether this configuration is above the shower infrared cutoff
	*/
	void isAboveCutoff(bool yes) { theIsAboveCutoff = yes; }

	//@}

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

	/**
	* Return true, if scales should be calculated from real emission kinematics
	*/
	bool realEmissionScales() const { return theRealEmissionScales; }

	/**
	* Switch on or off that scales should be calculated from real emission kinematics
	*/
	void doRealEmissionScales(bool on = true) { theRealEmissionScales = on; }

	/**
	* Return the scale associated with the phase space point provided
	* by the last call to setKinematics().
	*/
	virtual Energy2 scale() const {
	return realEmissionScales() ?
	realEmissionME()->scale() :
	underlyingBornME()->scale();
	}

	/**
	* 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 realEmissionScales() ?
	realEmissionME()->alphaS() :
	underlyingBornME()->alphaS();
	}

	/**
	* 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 realEmissionScales() ?
	realEmissionME()->alphaEM() :
	underlyingBornME()->alphaEM();
	}

	/**
	* Return true, if this matrix element provides the PDF
	* weight for the first incoming parton itself.
	*/
	virtual bool havePDFWeight1() const { return realEmissionME()->havePDFWeight1(); }

	/**
	* Return true, if this matrix element provides the PDF
	* weight for the second incoming parton itself.
	*/
	virtual bool havePDFWeight2() const { return realEmissionME()->havePDFWeight2(); }

	/**
	* How to sample the z-distribution.
	* FlatZ = 1
	* OneOverZ = 2
	* OneOverOneMinusZ = 3
	* OneOverZOneMinusZ = 4
	*/

	virtual int samplingZ() const {return 4;}
	//@}

	/** @name Matrix elements and evaluation */
	//@{

	/**
	* Return the real emission matrix element
	*/
	Ptr<MatchboxMEBase>::tcptr realEmissionME() const {
	return theRealEmissionME;
	}

	/**
	* Return the real emission matrix element
	*/
	Ptr<MatchboxMEBase>::tptr realEmissionME() {
	return theRealEmissionME;
	}

	/**
	* Set the real emission matrix element
	*/
	void realEmissionME(Ptr<MatchboxMEBase>::tptr me) { theRealEmissionME = me; }

	/**
	* Return the underlying Born matrix element
	*/
	Ptr<MatchboxMEBase>::tcptr underlyingBornME() const {
	return theUnderlyingBornME;
	}

	/**
	* Return the underlying Born matrix element
	*/
	Ptr<MatchboxMEBase>::tptr underlyingBornME() {
	return theUnderlyingBornME;
	}

	/**
	* Set the underlying Born matrix element
	*/
	void underlyingBornME(Ptr<MatchboxMEBase>::tptr me) { theUnderlyingBornME = me; }

	/**
	* Set the dipoles which have been found along with this dipole
	*/
	void partnerDipoles(const vector<Ptr<SubtractionDipole>::tptr>& p) {
	thePartners = p;
	}

	/**
	* Return the dipoles which have been found along with this dipole
	*/
	const vector<Ptr<SubtractionDipole>::tptr>& partnerDipoles() const {
	return thePartners;
	}

	/**
	* Return the matrix element averaged over spin correlations.
	*/
	virtual double me2Avg(double ccme2) const = 0;

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

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

	/**
	* Return true if this dipole is safe for clustering
	**/
	bool clustersafe() const;

	- bool clustersafe();
	-
	+
	+ /// calculate the general prefactor for merging.
	+ CrossSection prefactor(Energy2 factorizationScale)const;
	+
	/**
	* Calculate the parton shower approximation for this dipole.
	**/

	CrossSection ps(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const;

	/**
	* Calculate the dipole with clusterfsafe flag.
	**/

	CrossSection dip(Energy2 factorizationScale) const;



	/**
	* Calculate the dipole dSigDR and the parton shower approximation for this dipole.
	**/

	pair<CrossSection,CrossSection> dipandPs(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const;

	//@}

	/** @name Methods relevant to matching */
	//@{

	/**
	* Set the shower approximation.
	*/
	void showerApproximation(Ptr<ShowerApproximation>::tptr app) {
	theShowerApproximation = app;
	}

	/**
	* Return the shower approximation.
	*/
	Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }

	/**
	* Indicate that the shower real emission contribution should be subtracted.
	*/
	void doRealShowerSubtraction() { theRealShowerSubtraction = true; }

	/**
	* Return true, if the shower real emission contribution should be subtracted.
	*/
	bool realShowerSubtraction() const { return theRealShowerSubtraction; }

	/**
	* Indicate that the shower virtual contribution should be subtracted.
	*/
	void doVirtualShowerSubtraction() { theVirtualShowerSubtraction = true; }

	/**
	* Return true, if the shower virtual contribution should be subtracted.
	*/
	bool virtualShowerSubtraction() const { return theVirtualShowerSubtraction; }

	/**
	* Indicate that the loopsim matched virtual contribution should be subtracted.
	*/
	void doLoopSimSubtraction() { theLoopSimSubtraction = true; }

	/**
	* Return true, if the loopsim matched virtual contribution should be subtracted.
	*/
	bool loopSimSubtraction() const { return theLoopSimSubtraction; }

	//@}

	/** @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();

	/**
	* Indicate that the subtraction is being tested.
	*/
	void doTestSubtraction() { theSubtractionTest = true; }

	/**
	* Return true, if the subtraction is being tested.
	*/
	bool testSubtraction() const { return theSubtractionTest; }

	/**
	* Return true, if verbose
	*/
	bool verbose() const { return realEmissionME()->verbose() \|\| underlyingBornME()->verbose(); }

	/**
	* 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
	* generateTildeKinematics
	*/
	void logGenerateTildeKinematics() const;

	/**
	* Write out diagnostic information for
	* generateRadiationKinematics
	*/
	void logGenerateRadiationKinematics(const double * r) const;

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

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

	//@}

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

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

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

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

	//@}

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

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

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

	void cloneDependencies(const std::string& prefix , bool slim);

	//@}

	/** @name Methods required to setup the event record */
	//@{

	/**
	* construct the spin information for the interaction
	*/
	virtual void constructVertex(tSubProPtr sub);

	/**
	* construct the spin information for the interaction
	*/
	virtual void constructVertex(tSubProPtr sub, const ColourLines* cl);

	/**
	* Comlete a SubProcess object using the internal degrees of freedom
	* generated in the last generateKinematics() (and possible other
	* degrees of freedom which was intergated over in dSigHatDR(). This
	* default version does nothing. Will be made purely virtual in the
	* future.
	*/
	virtual void generateSubCollision(SubProcess & sub);

	/**
	* Alpha parameter as in Nagy
	* (http://arxiv.org/pdf/hep-ph/0307268v2.pdf) to restrict dipole
	* phase space
	*/
	double alpha() const;

	/*
	* True if phase space point is above the alpha cut for this dipole.
	*/

	bool aboveAlpha() const;

	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Standard Interfaced functions. */
	//@{

	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();

	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();
	//@}

	private:

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

	/**
	* Wether or not this dipole acts in splitting mode.
	*/
	bool theSplitting;

	/**
	* True, if should apply to process in the xcomb.
	*/
	bool theApply;

	/**
	* True, if the subtraction is being tested.
	*/
	bool theSubtractionTest;

	/**
	* True if cuts should be ignored
	*/
	bool theIgnoreCuts;

	/**
	* The real emission matrix element to be considered
	*/
	Ptr<MatchboxMEBase>::ptr theRealEmissionME;

	/**
	* The underlying Born matrix element
	*/
	Ptr<MatchboxMEBase>::ptr theUnderlyingBornME;

	/**
	* The dipoles which have been found along with this dipole
	*/
	vector<Ptr<SubtractionDipole>::tptr> thePartners;

	/**
	* The TildeKinematics to be used.
	*/
	Ptr<TildeKinematics>::ptr theTildeKinematics;

	/**
	* The InvertedTildeKinematics to be used.
	*/
	Ptr<InvertedTildeKinematics>::ptr theInvertedTildeKinematics;

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

	/**
	* The emitter as referred to by the real emission
	* matrix element.
	*/
	int theRealEmitter;

	/**
	* The emission as referred to by the real emission
	* matrix element.
	*/
	int theRealEmission;

	/**
	* The spectator as referred to by the real emission
	* matrix element.
	*/
	int theRealSpectator;

	/**
	* The subtraction parameters
	*/
	vector<double> theSubtractionParameters;

	/**
	* Map real emission diagrams to underlying Born diagrams
	* and tilde emitter/spectator.
	*/
	map<RealEmissionKey,UnderlyingBornInfo> theMergingMap;

	/**
	* Map underlying Born diagrams and tilde emitter/spectator
	* to real emission diagram containing the splitting.
	*/
	multimap<UnderlyingBornKey,RealEmissionInfo> theSplittingMap;

	/**
	* Map underlying Born diagrams to emitter/spectator pairs
	*/
	map<cPDVector,pair<int,int> > theIndexMap;

	/**
	* Map real emission processes to Born diagrams
	*/
	map<cPDVector,DiagramVector> theUnderlyingBornDiagrams;

	/**
	* Map Born processes to real emission diagrams
	*/
	map<cPDVector,DiagramVector> theRealEmissionDiagrams;

	/**
	* Map underlying Born diagrams to real emission diagrams.
	*/
	map<tcDiagPtr,tcDiagPtr> theBornToRealDiagrams;

	/**
	* Map real emission diagrams to underlying Born diagrams.
	*/
	map<tcDiagPtr,tcDiagPtr> theRealToBornDiagrams;

	/**
	* The last real emission key encountered
	*/
	RealEmissionKey lastRealEmissionKey;

	/**
	* The last underlying Born key encountered
	*/
	UnderlyingBornKey lastUnderlyingBornKey;

	/**
	* The last real emission info encountered
	*/
	multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator lastRealEmissionInfo;

	/**
	* The emitter as referred to by the underlying Born
	* matrix element.
	*/
	int theBornEmitter;

	/**
	* The spectator as referred to by the underlying Born
	* matrix element.
	*/
	int theBornSpectator;

	/**
	* The last scale as generated from the tilde mapping
	*/
	Energy theLastSubtractionScale;

	/**
	* The last scale as generated from the splitting mapping
	*/
	Energy theLastSplittingScale;

	/**
	* The last pt as generated from the tilde mapping
	*/
	Energy theLastSubtractionPt;

	/**
	* The last pt as generated from the splitting mapping
	*/
	Energy theLastSplittingPt;

	/**
	* The last z as generated from the tilde mapping
	*/
	double theLastSubtractionZ;

	/**
	* The last z as generated from the splitting mapping
	*/
	double theLastSplittingZ;

	/**
	* The shower approximation.
	*/
	Ptr<ShowerApproximation>::ptr theShowerApproximation;

	/**
	* True, if the shower real emission contribution should be subtracted.
	*/
	bool theRealShowerSubtraction;

	/**
	* True, if the shower virtual contribution should be subtracted.
	*/
	bool theVirtualShowerSubtraction;

	/**
	* True, if the loopsim matched virtual contribution should be subtracted.
	*/
	bool theLoopSimSubtraction;

	/**
	* True, if scales should be calculated from real emission kinematics
	*/
	bool theRealEmissionScales;

	/**
	* Return the shower hard scale encountered
	*/
	Energy theShowerHardScale;

	/**
	* The shower evolution scale encountered
	*/
	Energy theShowerScale;

	/**
	* The shower splitting variables encountered
	*/
	vector<double> theShowerParameters;

	/**
	* True, if this configuration is in the shower phase space
	*/
	bool theIsInShowerPhasespace;

	/**
	* True, if this configuration is above the shower infrared cutoff
	*/
	bool theIsAboveCutoff;

	private:

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

	};

	}

	#endif /* HERWIG_SubtractionDipole_H */
	diff --git a/Shower/Dipole/Merging/Makefile.am b/Shower/Dipole/Merging/Makefile.am
	--- a/Shower/Dipole/Merging/Makefile.am
	+++ b/Shower/Dipole/Merging/Makefile.am
	@@ -1,16 +1,15 @@
	noinst_LTLIBRARIES = libHwDipoleShowerMerging.la

	libHwDipoleShowerMerging_la_SOURCES = \
	Node.cc \
	Node.fh \
	Node.h \
	MScale.cc \
	MScale.h \
	MergingFactory.cc \
	MergingFactory.fh \
	MergingFactory.h \
	- Merger.fh\
	Merger.h\
	Merger.cc \
	MergingReweight.h \
	MergingReweight.cc
	diff --git a/Shower/Dipole/Merging/Merger.fh b/Shower/Dipole/Merging/Merger.fh
	deleted file mode 100644
	--- a/Shower/Dipole/Merging/Merger.fh
	+++ /dev/null
	@@ -1,18 +0,0 @@
	-// -- C++ --
	-//
	-// This is the forward declaration of the Merger class.
	-//
	-#ifndef HERWIG_Merger_FH
	-#define HERWIG_Merger_FH
	-
	-#include "ThePEG/Config/Pointers.h"
	-
	-namespace Herwig {
	-using namespace ThePEG;
	-
	-class Merger;
	-ThePEG_DECLARE_POINTERS(Merger,MergerPtr);
	-
	-}
	-
	-#endif
	diff --git a/Shower/Dipole/Merging/Merger.h b/Shower/Dipole/Merging/Merger.h
	--- a/Shower/Dipole/Merging/Merger.h
	+++ b/Shower/Dipole/Merging/Merger.h
	@@ -1,407 +1,407 @@
	- // -- C++ --
	+ /// -- C++ --
	//
	- // Merger.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	- // Copyright (C) 2002-2007 The Herwig Collaboration
	+ /// Merger.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	+ /// Copyright (C) 2002-2007 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.
	+ /// 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_Merger_H
	#define HERWIG_Merger_H
	//
	- // This is the declaration of the Merger class.
	+ /// This is the declaration of the Merger class.
	//
	#include "MergingFactory.fh"
	#include "Node.fh"



	#include "ThePEG/Handlers/HandlerBase.h"
	#include "Herwig/Shower/Dipole/DipoleShowerHandler.h"
	//#include "Herwig/Shower/Dipole/Base/DipoleSplittingGenerator.h"
	#include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h"

	#include "Herwig/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.h"

	#include "ThePEG/Cuts/JetFinder.h"
	#include "ThePEG/Cuts/Cuts.h"



	namespace Herwig {

	using namespace ThePEG;

	typedef Ptr<Node>::ptr NPtr;
	typedef vector<Ptr<Node>::ptr> NPtrVec;
	//definition of a history step
	struct HistoryStep {
	HistoryStep(){}
	- HistoryStep(NPtr cn,double w ,Energy sc){
	+ HistoryStep(NPtr cn, double w , Energy sc){
	node=cn;
	weight=w;
	scale=sc;
	}
	- // the cluster node defining this step
	- // containing the full information
	+ /// the cluster node defining this step
	+ /// containing the full information
	NPtr node;
	- // current sudakov weight of the history
	+ /// current sudakov weight of the history
	double weight;
	- // current scale of the history
	+ /// current scale of the history
	Energy scale;
	};

	typedef vector< HistoryStep > History;

	- typedef multimap<DipoleIndex,Ptr<DipoleSplittingGenerator>::ptr> GeneratorMap2;
	+ typedef multimap<DipoleIndex, Ptr<DipoleSplittingGenerator>::ptr> GeneratorMap2;

	/**
	* \ingroup DipoleShower
	* \author Johannes Bellm
	*
	* \brief Merger handles the Merger ....... //TODO .
	*
	* @see \ref MergerInterfaces "The interfaces"
	* defined for Merger.
	*/
	class Merger: public MergerBase {

	friend class MergingFactory;
	friend class Node;
	friend class MScale;

	public:

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

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

	public:

	- //define the ME region for a particle vector.
	- bool matrixElementRegion(PVector incoming,PVector outgoing,
	- Energy winnerScale=ZERO,
	+ // define the ME region for a particle vector.
	+ bool matrixElementRegion(PVector incoming, PVector outgoing,
	+ Energy winnerScale=ZERO,
	Energy cutscale=ZERO);
	- // return the current merging scale,
	- // gets smeared around the central merging scale in generate kinematics.
	+ /// return the current merging scale,
	+ /// gets smeared around the central merging scale in generate kinematics.
	Energy mergingScale()const{return theMergePt;}
	- // return the current merging pt (should be unified with mergingScale)
	+ /// return the current merging pt (should be unified with mergingScale)
	Energy mergePt()const {return theMergePt;}
	- // legsize of highest process with NLO corrections
	+ /// legsize of highest process with NLO corrections
	int M()const;
	- // legsize of the highest LO merged process
	+ /// legsize of the highest LO merged process
	int N()const;
	- // legsize of the production process
	+ /// legsize of the production process
	int N0()const{return theN0;}
	- // cross section of as given by the merging
	+ /// cross section of as given by the merging
	CrossSection MergingDSigDR() ;
	- // *** virtual functions of the base class **///
	- // set the current xcomb, called from ME
	- void setXComb(Ptr<MatchboxMEBase>::ptr,tStdXCombPtr);
	- // set kinematics, called from ME
	+ /// *** virtual functions of the base class **///
	+ /// set the current xcomb, called from ME
	+ void setXComb(Ptr<MatchboxMEBase>::ptr, tStdXCombPtr);
	+ /// set kinematics, called from ME
	void setKinematics(Ptr<MatchboxMEBase>::ptr);
	- // clear kinematics, called from ME
	+ /// clear kinematics, called from ME
	void clearKinematics(Ptr<MatchboxMEBase>::ptr);
	- // generate kinematics, called from ME
	- bool generateKinematics(Ptr<MatchboxMEBase>::ptr,const double *);
	- // fill the projector the subprocess is build from
	+ /// generate kinematics, called from ME
	+ bool generateKinematics(Ptr<MatchboxMEBase>::ptr, const double *);
	+ /// fill the projector the subprocess is build from
	void fillProjectors(Ptr<MatchboxMEBase>::ptr);
	- // return true true if a given ME phase space point is TODO
	- pair<bool,bool> clusterSafe(Ptr<MatchboxMEBase>::ptr,int,int,int);
	- // return the current maximum legs, the shower should veto
	+ /// return true true if a given ME phase space point is TODO
	+ pair<bool, bool> clusterSafe(Ptr<MatchboxMEBase>::ptr, int, int, int);
	+ /// return the current maximum legs, the shower should veto
	size_t maxLegs() const {return theCurrentMaxLegs;}
	- // set the current ME
	+ /// set the current ME
	void setME(Ptr<MatchboxMEBase>::ptr me){theCurrentME=me;}

	protected:
	- // the merging factory needs to set the legsize of the production process
	+ /// the merging factory needs to set the legsize of the production process
	void N0(int n){ theN0=n;}
	- // return the large-N basis (TODO: implement check if born ME works with the choice)
	+ /// return the large-N basis (TODO: implement check if born ME works with the choice)
	Ptr<ColourBasis>::ptr largeNBasis(){return theLargeNBasis;}
	- // smeare the merging pt
	+ /// smeare the merging pt
	void smeareMergePt(){theMergePt=centralMergePt()(1.+0.(-1.+2.UseRandom::rnd())smear());}
	- // true if the phase space for initial emissions should not be restricted in z.
	+ /// true if the phase space for initial emissions should not be restricted in z.
	bool openInitialStateZ(){return theOpenInitialStateZ;}
	- // set the current renormalisation scale
	+ /// set the current renormalisation scale
	Energy renormscale() { return therenormscale;}

	private:
	- // calculate a single sudakov step for a given dipole
	- double singlesudakov(Dipole,Energy,Energy,pair<bool,bool>);
	- // calculate the sudakov supression for a clusternode between
	- // the current running scale and next scale
	- bool dosudakov(NPtr Born,Energy running, Energy next, double& sudakov0_n);
	- // cleanup
	+ /// calculate a single sudakov step for a given dipole
	+ double singlesudakov(Dipole, Energy, Energy, pair<bool, bool>);
	+ /// calculate the sudakov supression for a clusternode between
	+ /// the current running scale and next scale
	+ bool dosudakov(NPtr Born, Energy running, Energy next, double& sudakov0_n);
	+ /// cleanup
	void cleanup(NPtr);
	- // return true if the cluster node has the matching number of
	- // legs to the current projector stage
	+ /// return true if the cluster node has the matching number of
	+ /// legs to the current projector stage
	bool isProjectorStage( NPtr , int );
	- // calculate the staring scale:
	- // if Node is part of the production process, calculate according to the
	- // scale choice object in the merging scale objekt, else
	- // return max(scale as scalechoice , min(Mass(i,j)))
	+ /// calculate the staring scale:
	+ /// if Node is part of the production process, calculate according to the
	+ /// scale choice object in the merging scale objekt, else
	+ /// return max(scale as scalechoice , min(Mass(i, j)))
	Energy CKKW_StartScale(NPtr);
	- // prepare the sudakov calculation
	- void CKKW_PrepareSudakov(NPtr,Energy);
	- // number of active flavours as given by the shower
	+ /// prepare the sudakov calculation
	+ void CKKW_PrepareSudakov(NPtr, Energy);
	+ /// number of active flavours as given by the shower
	double Nf(Energy scale)const{return DSH()->Nf(scale);}
	- // pointer to the factory
	+ /// pointer to the factory
	Ptr<MergingFactory>::ptr treefactory();
	- // map from ME to first clusternode
	- map<Ptr<MatchboxMEBase>::ptr,NPtr> firstNodeMap();
	- // set the current merging pt, smeared in generate kinematics
	+ /// map from ME to first clusternode
	+ map<Ptr<MatchboxMEBase>::ptr, NPtr> firstNodeMap();
	+ /// set the current merging pt, smeared in generate kinematics
	void mergePt(Energy x) {theMergePt = x;}
	- // return the central merging pt
	+ /// return the central merging pt
	Energy centralMergePt() {return theCentralMergePt;}

	private:
	- // calculate the history weighted born cross section
	+ /// calculate the history weighted born cross section
	CrossSection MergingDSigDRBornStandard(NPtr Node);
	- // calculate the history weighted born cross section
	- // add the difference of IPK with and without alpha parameter
	- // subtract the dipoles above the alpha parameter
	+ /// calculate the history weighted born cross section
	+ /// add the difference of IPK with and without alpha parameter
	+ /// subtract the dipoles above the alpha parameter
	CrossSection MergingDSigDRBornGamma(NPtr Node);
	- // calculate the history weighted virtual contribution
	+ /// calculate the history weighted virtual contribution
	CrossSection MergingDSigDRVirtualStandard(NPtr Node);
	- // calculate the history weighted real contribution
	- // splitted into 3 differnt contibutions
	+ /// calculate the history weighted real contribution
	+ /// splitted into 3 differnt contibutions
	CrossSection MergingDSigDRRealStandard(NPtr Node);
	- // calculate the history weighted real contribution
	- // all dipoles above:
	- // N*(R rnd(i)-Dip_i) history_i U(\phi^n_i)
	+ /// calculate the history weighted real contribution
	+ /// all dipoles above:
	+ /// N*(R rnd(i)-Dip_i) history_i U(\phi^n_i)
	CrossSection MergingDSigDRRealAllAbove(NPtr Node);
	- // calculate the history weighted real contribution
	- // not all dipoles above:
	- // (R - sum PS_i) history_rnd U(\phi^n+1)
	+ /// calculate the history weighted real contribution
	+ /// not all dipoles above:
	+ /// (R - sum PS_i) history_rnd U(\phi^n+1)
	CrossSection MergingDSigDRRealBelowSubReal(NPtr Node);
	- // calculate the history weighted real contribution
	- // not all dipoles above:
	- // rnd(i)-> N*(PS_i - Dip_i) history_i U(\phi^n_i)
	+ /// calculate the history weighted real contribution
	+ /// not all dipoles above:
	+ /// rnd(i)-> N*(PS_i - Dip_i) history_i U(\phi^n_i)
	CrossSection MergingDSigDRRealBelowSubInt(NPtr Node);
	- // max legssize the shower should veto for LO
	+ /// max legssize the shower should veto for LO
	size_t maxLegsLO() const {return N0()+N();}
	- // Calculate the LO partonic cross section.
	- // if diffalpha != 1, add the difference of IPK(1)-IPK(diffalpha)
	- CrossSection TreedSigDR(Energy startscale,NPtr,double diffalpha=1.);
	- // fill the projecting xcomb
	+ /// Calculate the LO partonic cross section.
	+ /// if diffalpha != 1, add the difference of IPK(1)-IPK(diffalpha)
	+ CrossSection TreedSigDR(Energy startscale, NPtr, double diffalpha=1.);
	+ /// fill the projecting xcomb
	Energy fillProjector(int);
	- // fill the history, including calculation of sudakov supression
	+ /// fill the history, including calculation of sudakov supression
	void fillHistory(Energy, NPtr, NPtr );
	- // calculate the pdf ratio for the given clusternode
	- double pdfratio(NPtr,Energy,Energy,int);
	- // return the pdf-ratio reweight for the history
	+ /// calculate the pdf ratio for the given clusternode
	+ double pdfratio(NPtr, Energy, Energy, int);
	+ /// return the pdf-ratio reweight for the history
	double pdfReweight();
	- // return the alpha_s reweight for the history
	+ /// return the alpha_s reweight for the history
	double alphaReweight();
	- // max legssize the shower should veto for NLO
	+ /// max legssize the shower should veto for NLO
	size_t maxLegsNLO()const {return N0()+M();}
	- // calculate the virtual contribution.
	- CrossSection LoopdSigDR(Energy startscale,NPtr);
	- // calculate alpha_s expansion of the pdf-ratios
	+ /// calculate the virtual contribution.
	+ CrossSection LoopdSigDR(Energy startscale, NPtr);
	+ /// calculate alpha_s expansion of the pdf-ratios
	double sumpdfReweightUnlops();
	- // calculate alpha_s expansion of the alpha_s-ratios, including K_g
	+ /// calculate alpha_s expansion of the alpha_s-ratios, including K_g
	double sumalphaReweightUnlops();
	- // calculate alpha_s expansion of the sudakov exponents
	+ /// calculate alpha_s expansion of the sudakov exponents
	double sumfillHistoryUnlops();
	- // calculate alpha_s expansion of the single step alpha_s-ratio, including K_g
	- double alphasUnlops( Energy next,Energy fixedScale);
	- // calculate alpha_s expansion of the single step pdf-ratio
	- double pdfUnlops(NPtr,int,Energy,Energy,double,int,Energy);
	- // calculate alpha_s expansion of the single step sudakov exponent
	- bool doUNLOPS(NPtr Born,Energy running, Energy next,Energy fixedScale, double& UNLOPS);
	- // calculate alpha_s expansion of the single dipole sudakov exponent
	- double singleUNLOPS(Dipole,Energy,Energy,Energy,pair<bool,bool>);
	+ /// calculate alpha_s expansion of the single step alpha_s-ratio, including K_g
	+ double alphasUnlops( Energy next, Energy fixedScale);
	+ /// calculate alpha_s expansion of the single step pdf-ratio
	+ double pdfUnlops(NPtr, int, Energy, Energy, double, int, Energy);
	+ /// calculate alpha_s expansion of the single step sudakov exponent
	+ bool doUNLOPS(NPtr Born, Energy running, Energy next, Energy fixedScale, double& UNLOPS);
	+ /// calculate alpha_s expansion of the single dipole sudakov exponent
	+ double singleUNLOPS(Dipole, Energy, Energy, Energy, pair<bool, bool>);
	//alpha_s as given in the shower
	double as(Energy q){return DSH()->as(q);}
	//return the dipole shower handler
	Ptr<DipoleShowerHandler>::ptr DSH(){return theDipoleShowerHandler;}
	//return the const dipole shower handler
	Ptr<DipoleShowerHandler>::ptr DSH()const{return theDipoleShowerHandler;}
	- // set the current renormalisation scale
	+ /// set the current renormalisation scale
	void renormscale(Energy x) { therenormscale = x;}
	- // insert map from ME to first clusternode
	- void firstNodeMap(Ptr<MatchboxMEBase>::ptr,NPtr);
	- // the gamma parameter to subtract dipoles above a alpha parameter
	- // and subtract the corresponding IPK operator
	+ /// insert map from ME to first clusternode
	+ void firstNodeMap(Ptr<MatchboxMEBase>::ptr, NPtr);
	+ /// the gamma parameter to subtract dipoles above a alpha parameter
	+ /// and subtract the corresponding IPK operator
	double gamma()const{return theGamma;}
	- // history choice: weighted history choice
	+ /// history choice: weighted history choice
	int chooseHistory()const {return theChooseHistory;}
	- // the smearing factor for the merging scale
	+ /// the smearing factor for the merging scale
	double smear()const{return theSmearing;}
	- // flag to tell if ME region shoulcd be defined by jet algorithm
	- // currently not implemented
	+ /// flag to tell if ME region shoulcd be defined by jet algorithm
	+ /// currently not implemented
	bool MERegionByJetAlg(){return defMERegionByJetAlg;}
	- // return the large-N colour basis
	+ /// return the large-N colour basis
	void largeNBasis(Ptr<ColourBasis>::ptr x){theLargeNBasis=x;}



	private:

	- // calculate the history expansion
	+ /// calculate the history expansion
	bool Unlopsweights;
	- // use CMW scheme
	+ /// use CMW scheme
	bool theCMWScheme;
	- // true if current point should be projected
	+ /// true if current point should be projected
	bool projected;
	- // true if LO cross sections should be unitarised
	+ /// true if LO cross sections should be unitarised
	bool isUnitarized;
	- // true if NLO contributions should be unitarised
	+ /// true if NLO contributions should be unitarised
	bool isNLOUnitarized;
	- // define ME region by jet algorithm
	+ /// define ME region by jet algorithm
	bool defMERegionByJetAlg;
	- // no z-restricions on initial state emissions in clustering
	+ /// no z-restricions on initial state emissions in clustering
	bool theOpenInitialStateZ;
	- // history weight choice
	+ /// history weight choice
	int theChooseHistory;
	- // legsize of production process
	+ /// legsize of production process
	int theN0;
	- // calculate only the N particle contribution
	+ /// calculate only the N particle contribution
	int theOnlyN;
	- // the current maxlegs (either LO or NLO maxlegs)
	+ /// the current maxlegs (either LO or NLO maxlegs)
	size_t theCurrentMaxLegs;
	- // renormalisation scale factor of the ME
	+ /// renormalisation scale factor of the ME
	double xiRenME;
	- // factorisation scale factor of the ME
	+ /// factorisation scale factor of the ME
	double xiFacME;
	- // renormalisation scale factor of the PS
	+ /// renormalisation scale factor of the PS
	double xiRenSh;
	- // factorisation scale factor of the PS
	+ /// factorisation scale factor of the PS
	double xiFacSh;
	- // starting scale factor of the PS
	+ /// starting scale factor of the PS
	double xiQSh;
	- // current weight and weight of clustered born
	- double weight,weightCB;
	- // subtract the dipole contribution above a given gamma
	+ /// current weight and weight of clustered born
	+ double weight, weightCB;
	+ /// subtract the dipole contribution above a given gamma
	double theGamma;
	- // if ME region defined by jet algorithm, this is the y cut for ee
	+ /// if ME region defined by jet algorithm, this is the y cut for ee
	double ee_ycut;
	- // if ME region defined by jet algorithm, this is the d cut for pp
	+ /// if ME region defined by jet algorithm, this is the d cut for pp
	double pp_dcut;
	- // smearing factor for merging scale
	+ /// smearing factor for merging scale
	double theSmearing;
	- // current renomalistion scale
	+ /// current renomalistion scale
	Energy therenormscale;
	- // cutoff for real emission contribution
	+ /// cutoff for real emission contribution
	Energy theIRSafePT;
	- // current merging scale
	+ /// current merging scale
	Energy theMergePt;
	- // central merging scale
	+ /// central merging scale
	Energy theCentralMergePt;
	- // current cluster histoy including sudakov weights
	+ /// current cluster histoy including sudakov weights
	History history;
	- // if ME region defined by jet algorithm, this is the jetfinder
	+ /// if ME region defined by jet algorithm, this is the jetfinder
	Ptr<JetFinder>::ptr theMergingJetFinder;
	- // pointer to the large-N basis
	+ /// pointer to the large-N basis
	Ptr<ColourBasis>::ptr theLargeNBasis;
	- // current ME
	+ /// current ME
	Ptr<MatchboxMEBase>::ptr theCurrentME;
	- // Tilde kinematics pointers, only to use lastPt(emitter,emission,spectator)
	+ /// Tilde kinematics pointers, only to use lastPt(emitter, emission, spectator)
	Ptr<FFLightTildeKinematics>::ptr FFLTK;
	Ptr<FILightTildeKinematics>::ptr FILTK;
	Ptr<IFLightTildeKinematics>::ptr IFLTK;
	Ptr<IILightTildeKinematics>::ptr IILTK;
	Ptr<FFMassiveTildeKinematics>::ptr FFMTK;
	Ptr<FIMassiveTildeKinematics>::ptr FIMTK;
	Ptr<IFMassiveTildeKinematics>::ptr IFMTK;
	//pointer to the shower handler
	Ptr<DipoleShowerHandler>::ptr theDipoleShowerHandler;
	- // pointer to the MergingFactory
	+ /// pointer to the MergingFactory
	Ptr<MergingFactory>::ptr theTreeFactory;
	- // map from ME to first Node
	- map<Ptr<MatchboxMEBase>::ptr,NPtr> theFirstNodeMap;
	+ /// map from ME to first Node
	+ map<Ptr<MatchboxMEBase>::ptr, NPtr> theFirstNodeMap;

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

	//@}


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


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


	private:



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

	};

	}



	#endif /* HERWIG_Merger_H */
	diff --git a/Shower/Dipole/Merging/MergingFactory.h b/Shower/Dipole/Merging/MergingFactory.h
	--- a/Shower/Dipole/Merging/MergingFactory.h
	+++ b/Shower/Dipole/Merging/MergingFactory.h
	@@ -1,184 +1,185 @@
	- // -- C++ --
	+ /// -- C++ --
	//
	- // MergingFactory.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	- // Copyright (C) 2002-2012 The Herwig Collaboration
	+ /// MergingFactory.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.
	+ /// 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_MergingFactory_H
	#define HERWIG_MergingFactory_H
	//
	- // This is the declaration of the MergingFactory class.
	+ /// This is the declaration of the MergingFactory class.
	//

	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
	#include "Node.fh"
	#include "Merger.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup Matchbox
	* \author Johannes Bellm
	*
	* \brief MergingFactory automatically sets up a NLO
	* QCD merging carried out in dipole subtraction.
	*
	* @see \ref MergingFactoryInterfaces "The interfaces"
	* defined for MergeboxFactory.
	*/
	class MergingFactory : public MatchboxFactory {
	public:

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

	/**
	* The destructor.
	*/
	virtual ~MergingFactory();
	//@}
	- // main method to setup the ME vector
	+ /// main method to setup the ME vector
	virtual void setup();
	- //fill all amplitudes, stored in pureMEsMap
	+ /// fill all amplitudes, stored in pureMEsMap
	void fill_amplitudes();
	/// prepare the Born and virtual matrix elements.
	void prepare_BV(int i);
	- //prepare the real emission matrix elements.
	+ /// prepare the real emission matrix elements.
	void prepare_R(int i);
	- // push the born contributions to the ME vector.
	- void pushB(Ptr<MatchboxMEBase>::ptr,int);
	+ /// push the born contributions to the ME vector.
	+ void pushB(Ptr<MatchboxMEBase>::ptr, int);
	//push the virtual contributions to the ME vector.
	- void pushV(Ptr<MatchboxMEBase>::ptr,int);
	- // push the real contributions to the ME vector.
	- void pushProR(Ptr<MatchboxMEBase>::ptr,int);
	- //order matrix elements form one loop provider.
	+ void pushV(Ptr<MatchboxMEBase>::ptr, int);
	+ /// push the real contributions to the ME vector.
	+ void pushProR(Ptr<MatchboxMEBase>::ptr, int);
	+ /// order matrix elements form one loop provider.
	void orderOLPs();
	- // Debugging: push only multiplicities to the ME vector
	- // in range of specified mulltiplicity.
	+ /// Debugging: push only multiplicities to the ME vector
	+ /// in range of specified mulltiplicity.
	int onlymulti()const {
	return theonlymulti==-1?-1:(theonlymulti+processMap.find(0)->second.size());}
	- // calculate only unlops weights.
	+ /// calculate only unlops weights.
	bool onlyUnlopsweights() const {return theonlyUnlopsweights;}
	- // pointer to the merging helper.
	+ /// pointer to the merging helper.
	Ptr<Merger>::ptr MH(){return theMergingHelper;}
	- // maximal NLO mulitplicity: 0=NLO corrections to the productio process.
	+ /// maximal NLO mulitplicity: 0=NLO corrections to the productio process.
	int M()const {return theM-1;}
	- // leg size of highest multiplicity.
	+ /// leg size of highest multiplicity.
	int N()const {return theN;}
	- // Return the Map of matrix elements to be considered
	- // (the Key is the number of additional jets)
	+ /// Return the Map of matrix elements to be considered
	+ /// (the Key is the number of additional jets)
	const map<int, vector<Ptr<MatchboxMEBase>::ptr> >& pureMEsMap() const {
	return thePureMEsMap;}
	- // Access the Map of matrix elements to be considered
	- // (the Key is the number of additional jets)
	+ /// Access the Map of matrix elements to be considered
	+ /// (the Key is the number of additional jets)
	map<int, vector<Ptr<MatchboxMEBase>::ptr> >& pureMEsMap() {
	return thePureMEsMap;
	}
	//Parse a process description
	virtual vector<string> parseProcess(string);

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

	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
	+ /**
	+ * Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	private:

	- // Calculate only virtual and real contributions.
	+ /// Calculate only virtual and real contributions.
	bool theonlyNLOParts;
	- // Calculate only virtual contributions.
	+ /// Calculate only virtual contributions.
	bool theonlyvirtualNLOParts;
	- // Calculate only real contributions.
	+ /// Calculate only real contributions.
	bool theonlyrealNLOParts;
	- // Calculate only expanded histories contributions.
	+ /// Calculate only expanded histories contributions.
	bool theonlyUnlopsweights;
	- // unitarize virtual and real contributions.
	+ /// unitarize virtual and real contributions.
	bool theunitarizeNLOParts;
	- // Calculate born contributions.
	+ /// Calculate born contributions.
	bool calc_born;
	- // Calculate virtual contributions.
	+ /// Calculate virtual contributions.
	bool calc_virtual;
	- // Calculate real contributions.
	+ /// Calculate real contributions.
	bool calc_real;
	- // unitarise the LO contributions.
	+ /// unitarise the LO contributions.
	bool unitarized;
	- // unitarise the NLO contributions.
	+ /// unitarise the NLO contributions.
	bool NLOunitarized;
	- // did run setup.
	+ /// did run setup.
	bool ransetup;
	- // Debugging: push only multiplicities to the ME vector
	- // in range of specified mulltiplicity.
	+ /// Debugging: push only multiplicities to the ME vector
	+ /// in range of specified mulltiplicity.
	int theonlymulti;
	- // calculate only the specified subprocess with no.
	+ /// calculate only the specified subprocess with no.
	int theonlysub;
	- // cut the subprocesses in equal size pieces.
	+ /// cut the subprocesses in equal size pieces.
	int divideSub;
	- // interface to calculate every # subprocess.
	+ /// interface to calculate every # subprocess.
	int divideSubNumber;
	- // maximal legsize for NLO corrections.
	+ /// maximal legsize for NLO corrections.
	int theM;
	- // maximal legsize for LO contributions.
	+ /// maximal legsize for LO contributions.
	int theN;
	- // Prefix for subtraction data.
	+ /// Prefix for subtraction data.
	string theSubtractionData;
	- // map for processes.
	+ /// map for processes.
	map< int, vector<string> > processMap;
	//The matrix elements: int = number of additional jets
	map< int, vector<Ptr<MatchboxMEBase>::ptr> > thePureMEsMap;
	- // map for virtual contributions
	+ /// map for virtual contributions
	map<int, vector<Ptr<MatchboxInsertionOperator>::ptr> > theVirtualsMap;
	- // the merging helper
	+ /// the merging helper
	Ptr<Merger>::ptr theMergingHelper;

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

	};

	}

	#endif /* HERWIG_MergingFactory_H */
	diff --git a/Shower/Dipole/Merging/Node.h b/Shower/Dipole/Merging/Node.h
	--- a/Shower/Dipole/Merging/Node.h
	+++ b/Shower/Dipole/Merging/Node.h
	@@ -1,245 +1,242 @@
	// -- C++ --
	#ifndef Herwig_Node_H
	#define Herwig_Node_H
	//
	// This is the declaration of the Node class.
	//
	//#include "Node.fh"
	#include "MergingFactory.fh"
	#include "Merger.h"


	#include "ThePEG/Interface/Interfaced.h"
	#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
	#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
	#include "Herwig/Shower/Dipole/Base/DipoleEventRecord.h"
	#include "ThePEG/MatrixElement/MEBase.h"
	#include "ThePEG/Config/Pointers.h"
	#include <vector>

	namespace Herwig {
	using namespace ThePEG;
	/**
	* Here is the documentation of the Node class.
	*
	* @see \ref NodeInterfaces "The interfaces"
	* defined for Node.
	*/

	/**
	- * Define the SafeClusterMap type map<pair<pair<emitter,emmision>,spectator >
	- * ,pair<first-clustering,second-clustering> >
	+ * Define the SafeClusterMap type map<pair<pair<emitter, emmision>, spectator >
	+ * , pair<first-clustering, second-clustering> >
	*/
	- typedef map<pair<pair<int,int>,int >,pair<bool,bool> > SafeClusterMap;
	+ typedef map<pair<pair<int, int>, int >, pair<bool, bool> > SafeClusterMap;

	class Node : public Interfaced {
	public:

	/** @name Standard constructors and destructors. */
	//@{

	- /**
	- * The default constructor. Do not use!
	- */
	+ ///The default constructor. Do not use!
	Node();
	// another constructor for first nodes
	- Node(Ptr<MatchboxMEBase>::ptr nodeME,int cutstage,Ptr<Merger>::ptr mh);
	+ Node(Ptr<MatchboxMEBase>::ptr nodeME, int cutstage, Ptr<Merger>::ptr mh);
	// another constructor for underlying nodes
	- Node(Ptr<Node>::ptr deephead,
	- Ptr<Node>::ptr head,
	- Ptr<SubtractionDipole>::ptr dipol,
	- Ptr<MatchboxMEBase>::ptr nodeME,
	+ Node(Ptr<Node>::ptr deephead,
	+ Ptr<Node>::ptr head,
	+ Ptr<SubtractionDipole>::ptr dipol,
	+ Ptr<MatchboxMEBase>::ptr nodeME,
	int cutstage);
	-
	- /**
	- * The destructor.
	- */
	+ /// The destructor.
	virtual ~Node();
	//@}

	public:
	// get children from vector<Ptr<MatchboxMEBase>::ptr>
	void birth(vector<Ptr<MatchboxMEBase>::ptr> vec);
	- // recursive setXComb. proStage is the number of clusterings
	- // before the projectors get filled.
	+ /// recursive setXComb. proStage is the number of clusterings
	+ /// before the projectors get filled.
	void setXComb(tStdXCombPtr xc);
	- // calculate the dipole and ps approximation
	- pair<CrossSection,CrossSection> calcDipandPS(Energy scale);
	- // calculate the ps approximation
	+ /// calculate the dipole and ps approximation
	+ pair<CrossSection, CrossSection> calcDipandPS(Energy scale);
	+ /// calculate the ps approximation
	CrossSection calcPs(Energy scale);
	- // calculate the dipole
	+ /// calculate the dipole
	CrossSection calcDip(Energy scale);
	- // recursive flush caches and clean up XCombs.
	+ /// recursive flush caches and clean up XCombs.
	void flushCaches();
	- // recursive clearKinematics
	+ /// recursive clearKinematics
	void clearKinematics();
	- // recursive setKinematics
	+ /// recursive setKinematics
	void setKinematics();
	- // recursive generateKinematics using tilde kinematics of the dipoles
	- bool generateKinematics(const double *r, int stage,Energy2 shat);
	- // generate the kinamatics of the first node
	+ /// recursive generateKinematics using tilde kinematics of the dipoles
	+ bool generateKinematics(const double *r, int stage, Energy2 shat);
	+ /// generate the kinamatics of the first node
	void firstgenerateKinematics(const double *r, int stage);
	//return the ME
	const Ptr<MatchboxMEBase>::ptr nodeME()const;
	//return the node ME
	Ptr<MatchboxMEBase>::ptr nodeME();
	//return the parent Node
	Ptr<Node>::ptr parent()const {return theparent;}
	- // vector of children nodes created in birth
	+ /// vector of children nodes created in birth
	vector< Ptr<Node>::ptr > children()const {return thechildren;}
	//pick a random child (flat)
	Ptr<Node>::ptr randomChild();
	- // true if all children show scales above pt
	+ /// true if all children show scales above pt
	bool allAbove(Energy pt);
	- // true if the node is in the history of other.
	+ /// true if the node is in the history of other.
	bool isInHistoryOf(Ptr<Node>::ptr other);
	- // legsize of the node ME
	+ /// legsize of the node ME
	int legsize() const;
	- // set the first node (first men). only use in factory
	+ /// set the first node (first men). only use in factory
	void deepHead(Ptr<Node>::ptr deephead) {theDeepHead = deephead;}
	- // return the first node
	+ /// return the first node
	Ptr<Node>::ptr deepHead() const {return theDeepHead;}
	- // insert nodes to projector vector
	+ /// insert nodes to projector vector
	void Projector(double a, Ptr<Node>::ptr pro) {
	- pair<double,Ptr<Node>::ptr> p;
	+ pair<double, Ptr<Node>::ptr> p;
	p.first = a;
	p.second = pro;
	theProjectors.push_back(p);
	}
	- // insert nodes to projector vector
	+ /// insert nodes to projector vector
	vector< pair <double , Ptr<Node>::ptr > > Projector() {return theProjectors;}
	- // returns the dipol of the node.
	+ /// returns the dipol of the node.
	Ptr<SubtractionDipole>::ptr dipole() const;
	- // set the xcomb of the node
	+ /// set the xcomb of the node
	void xcomb(StdXCombPtr xc) { thexcomb = xc;}
	- // return the xcomb
	+ /// return the xcomb
	StdXCombPtr xcomb() const {return thexcomb;}
	- // return the current running pt
	+ /// return the current running pt
	Energy runningPt(){return theRunningPt;}
	- // set the current running pt
	+ /// set the current running pt
	void runningPt(Energy x){theRunningPt=x;}
	- // return the cut stage to cut on merging pt in generate kinematics
	+ /// return the cut stage to cut on merging pt in generate kinematics
	int cutStage() const {return theCutStage;}
	- // get the clustersafe map for this node
	+ /// get the clustersafe map for this node
	SafeClusterMap clusterSafe() const {return clustersafer;}
	- // get a vector of the next nodes, ordered in pt (and in parton shower phace space)
	- vector<Ptr<Node>::ptr> getNextOrderedNodes(bool normal=true,double hardscalefactor=1.);
	+ /// get a vector of the next nodes, ordered in pt (and in parton shower phace space)
	+ vector<Ptr<Node>::ptr> getNextOrderedNodes(bool normal=true, double hardscalefactor=1.);
	//true if the node is in shower history for a given pt
	bool inShowerPS(Energy hardpt);
	//get the history
	- Ptr<Node>::ptr getHistory(bool normal=true,double hardscalefactor=1.);
	+ Ptr<Node>::ptr getHistory(bool normal=true, double hardscalefactor=1.);
	//true if node correspond to a subtracted real.
	bool subtractedReal() {return theSubtractedReal;}
	- // set if node correspont to a subtracted real.
	+ /// set if node correspont to a subtracted real.
	void subtractedReal(bool x) { theSubtractedReal = x;}
	//true if node correspond to a virtual contribution.
	bool virtualContribution() { return theVirtualContribution ;}
	- // set if node correspont to a virtual contribution.
	+ /// set if node correspont to a virtual contribution.
	void virtualContribution(bool x) {theVirtualContribution = x;}
	//pointer to the merging helper
	Ptr<Merger>::ptr MH()const{return theMergingHelper;}
	- // set the merging helper
	+ /// set the merging helper
	void MH(Ptr<Merger>::ptr a){theMergingHelper=a;}

	Energy pT()const{return dipole()->lastPt();}

	private:
	//the xcomb of the node
	StdXCombPtr thexcomb;
	- // the Matrixelement representing this node.
	+ /// the Matrixelement representing this node.
	Ptr<MatchboxMEBase>::ptr thenodeMEPtr;
	- // the dipol used to substract
	- // and generate kinematics using tilde kinematics
	+ /// the dipol used to substract
	+ /// and generate kinematics using tilde kinematics
	Ptr<SubtractionDipole>::ptr thedipol;
	- // vector of the children node
	+ /// vector of the children node
	vector< Ptr<Node>::ptr > thechildren;
	- // the parent node
	+ /// the parent node
	Ptr<Node>::ptr theparent;
	- // The nodes of the projection stage.
	- vector< pair <double,Ptr<Node>::ptr> > theProjectors;
	- // The godfather node of whole tree.(Firstnode)
	+ /// The nodes of the projection stage.
	+ vector< pair <double, Ptr<Node>::ptr> > theProjectors;
	+ /// The godfather node of whole tree.(Firstnode)
	Ptr<Node>::ptr theDeepHead;
	- // The CutStage is number of clusterings which are possible without
	- // introducing a merging scale to cut away singularities.
	- // -> subtracted MEs have the CutStage 1.
	- // -> virtual and normal tree level ME get 0.
	+ /**
	+ * The CutStage is number of clusterings which are possible without
	+ * introducing a merging scale to cut away singularities.
	+ * -> subtracted MEs have the CutStage 1.
	+ * -> virtual and normal tree level ME get 0.
	+ */
	int theCutStage;
	- // For [[Emitter,Emission],Spectator] the mapped pair gives
	- // information if the first and the second cluster is safe.
	+ /// For [[Emitter, Emission], Spectator] the mapped pair gives
	+ /// information if the first and the second cluster is safe.
	SafeClusterMap clustersafer;
	- // the current running pt
	+ /// the current running pt
	Energy theRunningPt;
	- // tell if node belongs to an ordered history
	+ /// tell if node belongs to an ordered history
	bool isOrdered;
	- // flag to tell if node is subtracted real
	+ /// flag to tell if node is subtracted real
	bool theSubtractedReal;
	- // flag to tell if node is virtual contribution
	+ /// flag to tell if node is virtual contribution
	bool theVirtualContribution;
	- // the merging helper
	+ /// the merging helper
	Ptr<Merger>::ptr theMergingHelper;

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


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



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

	};

	}

	#endif /* Herwig_Node_H */





	diff --git a/src/Matchbox/DiagonalCKM.in b/src/Matchbox/DiagonalCKM.in
	--- a/src/Matchbox/DiagonalCKM.in
	+++ b/src/Matchbox/DiagonalCKM.in
	@@ -1,7 +1,7 @@
	# -- ThePEG-repository --

	set /Herwig/MatrixElements/Matchbox/Factory:QuarkFlavourDiagonal Yes
	-set /Herwig/Merging/MFactory:QuarkFlavourDiagonal Yes
	+set /Herwig/Merging/MergingFactory:QuarkFlavourDiagonal Yes

	set /Herwig/Vertices/FFWMatchboxVertex:Diagonal Yes

	diff --git a/src/Matchbox/EECollider.in b/src/Matchbox/EECollider.in
	--- a/src/Matchbox/EECollider.in
	+++ b/src/Matchbox/EECollider.in
	@@ -1,28 +1,28 @@
	# -- ThePEG-repository --

	read Matchbox/GenericCollider.in

	cd /Herwig/EventHandlers

	create ThePEG::FixedCMSLuminosity Luminosity FixedCMSLuminosity.so
	set Luminosity:Energy 91.2
	set EventHandler:LuminosityFunction Luminosity
	set EventHandler:BeamA /Herwig/Particles/e-
	set EventHandler:BeamB /Herwig/Particles/e+
	set /Herwig/Shower/ShowerHandler:MPIHandler NULL
	set /Herwig/Shower/PowhegShowerHandler:MPIHandler NULL
	set /Herwig/DipoleShower/DipoleShowerHandler:MPIHandler NULL
	cd /Herwig/MatrixElements/Matchbox
	set Factory:FirstPerturbativePDF No
	set Factory:SecondPerturbativePDF No
	set Factory:PartonExtractor /Herwig/Partons/EEExtractor

	-set /Herwig/Merging/MFactory:FirstPerturbativePDF No
	-set /Herwig/Merging/MFactory:SecondPerturbativePDF No
	-set /Herwig/Merging/MFactory:PartonExtractor /Herwig/Partons/EEExtractor
	+set /Herwig/Merging/MergingFactory:FirstPerturbativePDF No
	+set /Herwig/Merging/MergingFactory:SecondPerturbativePDF No
	+set /Herwig/Merging/MergingFactory:PartonExtractor /Herwig/Partons/EEExtractor

	set /Herwig/Shower/PartnerFinder:QEDPartner IIandFF

	cd /Herwig/Cuts
	set Cuts:MHatMin 22.36*GeV

	diff --git a/src/Matchbox/EPCollider.in b/src/Matchbox/EPCollider.in
	--- a/src/Matchbox/EPCollider.in
	+++ b/src/Matchbox/EPCollider.in
	@@ -1,42 +1,42 @@
	# -- ThePEG-repository --

	read Matchbox/GenericCollider.in

	cd /Herwig/EventHandlers
	create ThePEG::LuminosityFunction Luminosity
	set EventHandler:BeamA /Herwig/Particles/p+
	set Luminosity:BeamEMaxA 920.*GeV
	set EventHandler:BeamB /Herwig/Particles/e-
	set Luminosity:BeamEMaxB 30.*GeV
	set EventHandler:LuminosityFunction Luminosity
	set /Herwig/Shower/ShowerHandler:MPIHandler NULL
	set /Herwig/Shower/PowhegShowerHandler:MPIHandler NULL
	set /Herwig/DipoleShower/DipoleShowerHandler:MPIHandler NULL
	cd /Herwig/MatrixElements/Matchbox
	set Factory:FirstPerturbativePDF No
	set Factory:SecondPerturbativePDF Yes
	set Factory:PartonExtractor /Herwig/Partons/DISExtractor

	-set /Herwig/Merging/MFactory:FirstPerturbativePDF No
	-set /Herwig/Merging/MFactory:SecondPerturbativePDF Yes
	-set /Herwig/Merging/MFactory:PartonExtractor /Herwig/Partons/DISExtractor
	+set /Herwig/Merging/MergingFactory:FirstPerturbativePDF No
	+set /Herwig/Merging/MergingFactory:SecondPerturbativePDF Yes
	+set /Herwig/Merging/MergingFactory:PartonExtractor /Herwig/Partons/DISExtractor


	do /Herwig/Shower/SplittingGenerator:DeleteFinalSplitting e-->e-,gamma; /Herwig/Shower/LtoLGammaSudakov
	do /Herwig/Shower/SplittingGenerator:DeleteFinalSplitting mu-->mu-,gamma; /Herwig/Shower/LtoLGammaSudakov
	do /Herwig/Shower/SplittingGenerator:DeleteFinalSplitting tau-->tau-,gamma; /Herwig/Shower/LtoLGammaSudakov
	set /Herwig/Shower/KinematicsReconstructor:ReconstructionOption Colour

	cd /Herwig/Cuts

	set Cuts:ScaleMin 1.0*GeV
	set Cuts:X1Min 1.0e-5
	set Cuts:X2Min 1.0e-5
	set Cuts:MHatMin 0.0*GeV
	set NeutralCurrentCut:MinW2 0
	set NeutralCurrentCut:MaxW2 1e10
	set ChargedCurrentCut:MinW2 0
	set ChargedCurrentCut:MaxW2 1e10
	insert Cuts:TwoCuts[0] NeutralCurrentCut
	insert Cuts:TwoCuts[1] ChargedCurrentCut

	diff --git a/src/Matchbox/FiveFlavourNoBMassScheme.in b/src/Matchbox/FiveFlavourNoBMassScheme.in
	--- a/src/Matchbox/FiveFlavourNoBMassScheme.in
	+++ b/src/Matchbox/FiveFlavourNoBMassScheme.in
	@@ -1,102 +1,102 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox

	do Factory:StartParticleGroup p
	insert Factory:ParticleGroup 0 /Herwig/Particles/b
	insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/c
	insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/s
	insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/d
	insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/u
	insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
	insert Factory:ParticleGroup 0 /Herwig/Particles/g
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup pbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/b
	insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/c
	insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/s
	insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/d
	insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/u
	insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
	insert Factory:ParticleGroup 0 /Herwig/Particles/g
	do Factory:EndParticleGroup


	cd /Herwig/Merging/



	-do MFactory:StartParticleGroup p
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/b
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/bbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/c
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/cbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/s
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/sbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/d
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/dbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/u
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/ubar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/g
	-do MFactory:EndParticleGroup
	+do MergingFactory:StartParticleGroup p
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/b
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/bbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/c
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/cbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/s
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/sbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/d
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/dbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/u
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/ubar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/g
	+do MergingFactory:EndParticleGroup

	-do MFactory:StartParticleGroup pbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/b
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/bbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/c
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/cbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/s
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/sbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/d
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/dbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/u
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/ubar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/g
	-do MFactory:EndParticleGroup
	+do MergingFactory:StartParticleGroup pbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/b
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/bbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/c
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/cbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/s
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/sbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/d
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/dbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/u
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/ubar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/g
	+do MergingFactory:EndParticleGroup

	cd /Herwig/Particles

	do b:UnsetHardProcessMass
	do bbar:UnsetHardProcessMass
	do c:UnsetHardProcessMass
	do cbar:UnsetHardProcessMass
	set b:NominalMass 0*GeV
	set bbar:NominalMass 0*GeV
	set c:NominalMass 0*GeV
	set cbar:NominalMass 0*GeV

	cd /Herwig/DipoleShower/Kernels

	cp IFgx2qqxDipoleKernel IFgx2bbbarxDipoleKernel
	set IFgx2bbbarxDipoleKernel:Flavour /Herwig/Particles/b
	insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2bbbarxDipoleKernel

	cp IFgx2qqxDipoleKernel IFgx2bbarbxDipoleKernel
	set IFgx2bbarbxDipoleKernel:Flavour /Herwig/Particles/bbar
	insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2bbarbxDipoleKernel

	cp IFMgx2qqxDipoleKernel IFMgx2bbbarxDipoleKernel
	set IFMgx2bbbarxDipoleKernel:Flavour /Herwig/Particles/b
	insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFMgx2bbbarxDipoleKernel

	cp IFMgx2qqxDipoleKernel IFMgx2bbarbxDipoleKernel
	set IFMgx2bbarbxDipoleKernel:Flavour /Herwig/Particles/bbar
	insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFMgx2bbarbxDipoleKernel

	cp IIgx2qqxDipoleKernel IIgx2bbbarxDipoleKernel
	set IIgx2bbbarxDipoleKernel:Flavour /Herwig/Particles/b
	insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2bbbarxDipoleKernel

	cp IIgx2qqxDipoleKernel IIgx2bbarbxDipoleKernel
	set IIgx2bbarbxDipoleKernel:Flavour /Herwig/Particles/bbar
	insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2bbarbxDipoleKernel

	diff --git a/src/Matchbox/FiveFlavourScheme.in b/src/Matchbox/FiveFlavourScheme.in
	--- a/src/Matchbox/FiveFlavourScheme.in
	+++ b/src/Matchbox/FiveFlavourScheme.in
	@@ -1,69 +1,69 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox

	do Factory:StartParticleGroup p
	insert Factory:ParticleGroup 0 /Herwig/Particles/b
	insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/c
	insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/s
	insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/d
	insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/u
	insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
	insert Factory:ParticleGroup 0 /Herwig/Particles/g
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup pbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/b
	insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/c
	insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/s
	insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/d
	insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/u
	insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
	insert Factory:ParticleGroup 0 /Herwig/Particles/g
	do Factory:EndParticleGroup


	cd /Herwig/Merging/



	-do MFactory:StartParticleGroup p
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/b
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/bbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/c
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/cbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/s
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/sbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/d
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/dbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/u
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/ubar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/g
	-do MFactory:EndParticleGroup
	+do MergingFactory:StartParticleGroup p
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/b
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/bbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/c
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/cbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/s
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/sbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/d
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/dbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/u
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/ubar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/g
	+do MergingFactory:EndParticleGroup

	-do MFactory:StartParticleGroup pbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/b
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/bbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/c
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/cbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/s
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/sbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/d
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/dbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/u
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/ubar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/g
	-do MFactory:EndParticleGroup
	+do MergingFactory:StartParticleGroup pbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/b
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/bbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/c
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/cbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/s
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/sbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/d
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/dbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/u
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/ubar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/g
	+do MergingFactory:EndParticleGroup

	cd /Herwig/Particles

	set b:HardProcessMass 0*GeV
	set bbar:HardProcessMass 0*GeV
	diff --git a/src/Matchbox/FourFlavourScheme.in b/src/Matchbox/FourFlavourScheme.in
	--- a/src/Matchbox/FourFlavourScheme.in
	+++ b/src/Matchbox/FourFlavourScheme.in
	@@ -1,77 +1,77 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox

	do Factory:StartParticleGroup p
	insert Factory:ParticleGroup 0 /Herwig/Particles/c
	insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/s
	insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/d
	insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/u
	insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
	insert Factory:ParticleGroup 0 /Herwig/Particles/g
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup pbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/c
	insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/s
	insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/d
	insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/u
	insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
	insert Factory:ParticleGroup 0 /Herwig/Particles/g
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup j
	insert Factory:ParticleGroup 0 /Herwig/Particles/c
	insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/s
	insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/d
	insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/u
	insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
	insert Factory:ParticleGroup 0 /Herwig/Particles/g
	do Factory:EndParticleGroup

	cd /Herwig/Merging/
	-do MFactory:StartParticleGroup p
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/c
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/cbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/s
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/sbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/d
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/dbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/u
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/ubar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/g
	-do MFactory:EndParticleGroup
	+do MergingFactory:StartParticleGroup p
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/c
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/cbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/s
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/sbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/d
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/dbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/u
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/ubar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/g
	+do MergingFactory:EndParticleGroup

	-do MFactory:StartParticleGroup pbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/c
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/cbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/s
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/sbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/d
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/dbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/u
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/ubar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/g
	-do MFactory:EndParticleGroup
	+do MergingFactory:StartParticleGroup pbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/c
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/cbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/s
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/sbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/d
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/dbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/u
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/ubar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/g
	+do MergingFactory:EndParticleGroup


	-do MFactory:StartParticleGroup j
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/c
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/cbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/s
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/sbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/d
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/dbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/u
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/ubar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/g
	-do MFactory:EndParticleGroup
	+do MergingFactory:StartParticleGroup j
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/c
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/cbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/s
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/sbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/d
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/dbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/u
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/ubar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/g
	+do MergingFactory:EndParticleGroup
	diff --git a/src/Matchbox/GoSam-GoSam.in b/src/Matchbox/GoSam-GoSam.in
	--- a/src/Matchbox/GoSam-GoSam.in
	+++ b/src/Matchbox/GoSam-GoSam.in
	@@ -1,17 +1,17 @@
	# -- ThePEG-repository --

	# GoSam as the one and only matrix element provider

	cd /Herwig/MatrixElements/Matchbox/Amplitudes

	clear /Herwig/MatrixElements/Matchbox/Factory:Amplitudes
	insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 GoSam

	-clear /Herwig/Merging/MFactory:Amplitudes
	-insert /Herwig/Merging/MFactory:Amplitudes 0 GoSam
	+clear /Herwig/Merging/MergingFactory:Amplitudes
	+insert /Herwig/Merging/MergingFactory:Amplitudes 0 GoSam

	cd /Herwig/Particles

	do b:UnsetHardProcessMass
	do bbar:UnsetHardProcessMass

	diff --git a/src/Matchbox/HJets.in b/src/Matchbox/HJets.in
	--- a/src/Matchbox/HJets.in
	+++ b/src/Matchbox/HJets.in
	@@ -1,13 +1,13 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox

	clear Factory:Amplitudes
	-clear /Herwig/Merging/MFactory:Amplitudes
	+clear /Herwig/Merging/MergingFactory:Amplitudes
	library HJets.so
	read HJetsProcesses.in

	cd /Herwig/Particles

	set b:HardProcessMass 0*GeV
	set bbar:HardProcessMass 0*GeV
	diff --git a/src/Matchbox/IdentifiedBs.in b/src/Matchbox/IdentifiedBs.in
	--- a/src/Matchbox/IdentifiedBs.in
	+++ b/src/Matchbox/IdentifiedBs.in
	@@ -1,33 +1,33 @@
	# -- ThePEG-repository --

	read Matchbox/FourFlavourScheme.in

	cd /Herwig/MatrixElements/Matchbox

	do Factory:StartParticleGroup j
	insert Factory:ParticleGroup 0 /Herwig/Particles/c
	insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/s
	insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/d
	insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/u
	insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
	insert Factory:ParticleGroup 0 /Herwig/Particles/g
	do Factory:EndParticleGroup

	cd /Herwig/Merging

	-do MFactory:StartParticleGroup j
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/c
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/cbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/s
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/sbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/d
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/dbar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/u
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/ubar
	-insert MFactory:ParticleGroup 0 /Herwig/Particles/g
	-do MFactory:EndParticleGroup
	+do MergingFactory:StartParticleGroup j
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/c
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/cbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/s
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/sbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/d
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/dbar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/u
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/ubar
	+insert MergingFactory:ParticleGroup 0 /Herwig/Particles/g
	+do MergingFactory:EndParticleGroup


	diff --git a/src/Matchbox/IncreaseVerbosity.in b/src/Matchbox/IncreaseVerbosity.in
	--- a/src/Matchbox/IncreaseVerbosity.in
	+++ b/src/Matchbox/IncreaseVerbosity.in
	@@ -1,13 +1,13 @@
	# -- ThePEG-repository --

	#=============================================================================
	# ::::: Verbosity options for Matchbox:::::
	#
	# Settings increase the verbosity of the Matchbox interface
	#=============================================================================

	set /Herwig/Samplers/Sampler:Verbose On
	set /Herwig/MatrixElements/Matchbox/Factory:InitVerbose On
	-set /Herwig/Merging/MFactory:InitVerbose On
	+set /Herwig/Merging/MergingFactory:InitVerbose On
	cd /

	diff --git a/src/Matchbox/MadGraph-GoSam.in b/src/Matchbox/MadGraph-GoSam.in
	--- a/src/Matchbox/MadGraph-GoSam.in
	+++ b/src/Matchbox/MadGraph-GoSam.in
	@@ -1,24 +1,24 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox/Amplitudes

	clear /Herwig/MatrixElements/Matchbox/Factory:Amplitudes
	insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 GenericProcesses

	set GenericProcesses:TreeLevelAmplitude MadGraph
	set GenericProcesses:OneLoopAmplitude GoSam

	-clear /Herwig/Merging/MFactory:Amplitudes
	-insert /Herwig/Merging/MFactory:Amplitudes 0 GenericProcesses
	+clear /Herwig/Merging/MergingFactory:Amplitudes
	+insert /Herwig/Merging/MergingFactory:Amplitudes 0 GenericProcesses


	create ColorFull::TraceBasis /Herwig/Merging/LargeNColourBasis2
	set /Herwig/Merging/LargeNColourBasis2:LargeN On
	set /Herwig/Merging/Merger:LargeNBasis LargeNColourBasis2


	cd /Herwig/Particles

	do b:UnsetHardProcessMass
	do bbar:UnsetHardProcessMass

	diff --git a/src/Matchbox/MadGraph-MadGraph.in b/src/Matchbox/MadGraph-MadGraph.in
	--- a/src/Matchbox/MadGraph-MadGraph.in
	+++ b/src/Matchbox/MadGraph-MadGraph.in
	@@ -1,23 +1,23 @@
	# -- ThePEG-repository --

	# MadGraph as the one and only matrix element provider

	cd /Herwig/MatrixElements/Matchbox/Amplitudes

	clear /Herwig/MatrixElements/Matchbox/Factory:Amplitudes
	insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 MadGraph


	-clear /Herwig/Merging/MFactory:Amplitudes
	-insert /Herwig/Merging/MFactory:Amplitudes 0 MadGraph
	+clear /Herwig/Merging/MergingFactory:Amplitudes
	+insert /Herwig/Merging/MergingFactory:Amplitudes 0 MadGraph


	create ColorFull::TraceBasis /Herwig/Merging/LargeNColourBasis2
	set /Herwig/Merging/LargeNColourBasis2:LargeN On
	set /Herwig/Merging/Merger:LargeNBasis /Herwig/Merging/LargeNColourBasis2

	cd /Herwig/Particles

	do b:UnsetHardProcessMass
	do bbar:UnsetHardProcessMass

	diff --git a/src/Matchbox/MadGraph-NJet.in b/src/Matchbox/MadGraph-NJet.in
	--- a/src/Matchbox/MadGraph-NJet.in
	+++ b/src/Matchbox/MadGraph-NJet.in
	@@ -1,23 +1,23 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox/Amplitudes

	clear /Herwig/MatrixElements/Matchbox/Factory:Amplitudes
	insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 GenericProcesses

	set GenericProcesses:TreeLevelAmplitude MadGraph
	set GenericProcesses:OneLoopAmplitude NJet

	-clear /Herwig/Merging/MFactory:Amplitudes
	-insert /Herwig/Merging/MFactory:Amplitudes 0 GenericProcesses
	+clear /Herwig/Merging/MergingFactory:Amplitudes
	+insert /Herwig/Merging/MergingFactory:Amplitudes 0 GenericProcesses


	create ColorFull::TraceBasis /Herwig/Merging/LargeNColourBasis2
	set /Herwig/Merging/LargeNColourBasis2:LargeN On
	set /Herwig/Merging/Merger:LargeNBasis LargeNColourBasis2

	cd /Herwig/Particles

	do b:UnsetHardProcessMass
	do bbar:UnsetHardProcessMass

	diff --git a/src/Matchbox/MadGraph-OpenLoops.in b/src/Matchbox/MadGraph-OpenLoops.in
	--- a/src/Matchbox/MadGraph-OpenLoops.in
	+++ b/src/Matchbox/MadGraph-OpenLoops.in
	@@ -1,22 +1,22 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox/Amplitudes

	clear /Herwig/MatrixElements/Matchbox/Factory:Amplitudes
	insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 GenericProcesses

	set GenericProcesses:TreeLevelAmplitude MadGraph
	set GenericProcesses:OneLoopAmplitude OpenLoops

	-clear /Herwig/Merging/MFactory:Amplitudes
	-insert /Herwig/Merging/MFactory:Amplitudes 0 GenericProcesses
	+clear /Herwig/Merging/MergingFactory:Amplitudes
	+insert /Herwig/Merging/MergingFactory:Amplitudes 0 GenericProcesses

	create ColorFull::TraceBasis /Herwig/Merging/LargeNColourBasis2
	set /Herwig/Merging/LargeNColourBasis2:LargeN On
	set /Herwig/Merging/Merger:LargeNBasis /Herwig/Merging/LargeNColourBasis2

	cd /Herwig/Particles

	do b:UnsetHardProcessMass
	do bbar:UnsetHardProcessMass

	diff --git a/src/Matchbox/Merging.in b/src/Matchbox/Merging.in
	--- a/src/Matchbox/Merging.in
	+++ b/src/Matchbox/Merging.in
	@@ -1,236 +1,236 @@
	cd /Herwig/Merging

	clear /Herwig/EventHandlers/EventHandler:SubProcessHandlers
	-insert /Herwig/EventHandlers/EventHandler:SubProcessHandlers[0] MFactory
	-set /Herwig/Cuts/MatchboxJetMatcher:Factory MFactory
	-set /Herwig/MatrixElements/Matchbox/MECorrectionHandler:Factory MFactory
	+insert /Herwig/EventHandlers/EventHandler:SubProcessHandlers[0] MergingFactory
	+set /Herwig/Cuts/MatchboxJetMatcher:Factory MergingFactory
	+set /Herwig/MatrixElements/Matchbox/MECorrectionHandler:Factory MergingFactory


	-set MFactory:Cuts /Herwig/Cuts/Cuts
	-set MFactory:ScaleChoice /Herwig/Merging/MScale
	-set Merger:MFactory MFactory
	+set MergingFactory:Cuts /Herwig/Cuts/Cuts
	+set MergingFactory:ScaleChoice /Herwig/Merging/MScale
	+set Merger:MergingFactory MergingFactory
	set Merger:DipoleShowerHandler /Herwig/DipoleShower/DipoleShowerHandler
	set /Herwig/DipoleShower/DipoleShowerHandler:MergingHelper Merger
	set MScale:MergingHelper Merger
	set Merger:MergingJetFinder /Herwig/Cuts/JetFinder
	-set MFactory:MergingHelper Merger
	+set MergingFactory:MergingHelper Merger


	cd /Herwig/Analysis
	create ThePEG::RivetAnalysis Rivet RivetAnalysis.so



	set /Herwig/Shower/ShowerHandler:HardScaleProfile NULL
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleProfile NULL

	set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF
	set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF



	set /Herwig/Couplings/NLOAlphaS:min_active_flavours 3
	set /Herwig/Couplings/NLOAlphaS:max_active_flavours 5
	set /Herwig/Couplings/NLOAlphaS:input_alpha_s 0.118
	set /Herwig/Couplings/NLOAlphaS:input_scale 91.18*GeV


	set /Herwig/Model:QCD/RunningAlphaS /Herwig/Couplings/NLOAlphaS
	set /Herwig/DipoleShower/DipoleShowerHandler:GlobalAlphaS /Herwig/Couplings/NLOAlphaS




	cd /Herwig/Generators
	set EventGenerator:EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler

	set /Herwig/DipoleShower/DipoleShowerHandler:MaxPtIsMuF Yes

	create Herwig::MonacoSampler /Herwig/Samplers/Monaco
	set /Herwig/Samplers/Sampler:BinSampler /Herwig/Samplers/Monaco


	cd /Herwig/Particles

	do c:UnsetHardProcessMass
	do cbar:UnsetHardProcessMass
	#do b:UnsetHardProcessMass
	#do bbar:UnsetHardProcessMass
	set c:NominalMass 0*GeV
	set cbar:NominalMass 0*GeV


	##################################################
	## intrinsic pt
	##################################################
	set /Herwig/DipoleShower/IntrinsicPtGenerator:ValenceIntrinsicPtScale 2.0*GeV
	set /Herwig/DipoleShower/IntrinsicPtGenerator:SeaIntrinsicPtScale 2.0*GeV

	##################################################
	## Dipole shower tune settings
	##################################################

	set /Herwig/Particles/g:ConstituentMass 0.7899848*GeV

	##################################################
	## cutoffs
	##################################################
	set /Herwig/DipoleShower/Kinematics/FFLightKinematics:IRCutoff 0.45*GeV
	set /Herwig/DipoleShower/Kinematics/FFMassiveKinematics:IRCutoff 0.45*GeV

	set /Herwig/DipoleShower/Kinematics/FILightKinematics:IRCutoff 0.45*GeV
	set /Herwig/DipoleShower/Kinematics/FIMassiveKinematics:IRCutoff 0.45*GeV

	set /Herwig/DipoleShower/Kinematics/IFLightKinematics:IRCutoff 0.45*GeV
	set /Herwig/DipoleShower/Kinematics/IFMassiveKinematics:IRCutoff 0.45*GeV

	set /Herwig/DipoleShower/Kinematics/IILightKinematics:IRCutoff 0.45*GeV

	## needs to be synchronized with settings above
	set /Herwig/MatrixElements/Matchbox/DipoleMatching:FFPtCut 0.45*GeV
	set /Herwig/MatrixElements/Matchbox/DipoleMatching:FIPtCut 0.45*GeV
	set /Herwig/MatrixElements/Matchbox/DipoleMatching:IIPtCut 0.45*GeV

	##################################################
	## hadronization parameters -- complete?
	##################################################
	cd /Herwig/Hadronization

	set ClusterFissioner:ClMaxLight 3.30
	set ClusterFissioner:ClPowLight 2.50
	set ClusterFissioner:PSplitLight 1.29
	set ClusterDecayer:ClDirLight 1
	set ClusterDecayer:ClSmrLight 3.118342

	set ClusterFissioner:ClMaxCharm 3.11*GeV
	set ClusterFissioner:ClPowCharm 1.62
	set ClusterFissioner:PSplitCharm 2.54
	set ClusterDecayer:ClDirCharm 1
	set ClusterDecayer:ClSmrCharm 0.
	set HadronSelector:SingleHadronLimitCharm 0.0

	set ClusterFissioner:ClMaxBottom 1.38*GeV
	set ClusterFissioner:ClPowBottom 0.62
	set ClusterFissioner:PSplitBottom 0.20
	set ClusterDecayer:ClDirBottom 1
	set ClusterDecayer:ClSmrBottom 0.
	set HadronSelector:SingleHadronLimitBottom 0.12

	set HadronSelector:PwtUquark 1.0
	set HadronSelector:PwtDquark 1.0
	set HadronSelector:PwtSquark 1.09
	set HadronSelector:PwtCquark 1.0
	set HadronSelector:PwtBquark 1.0
	set HadronSelector:PwtDIquark 0.66
	set HadronSelector:SngWt 1.0
	set HadronSelector:DecWt 1.0


	cd /










	cd /Herwig/Samplers

	set Monaco:EnhancementFactor 1.2
	set Monaco:InitialPoints 5000
	set Monaco:LuminosityMapperBins 8
	set Monaco:NIterations 4
	set Monaco:RemapChannelDimension Yes
	set Monaco:RemapperMinSelection 0.0001
	set Monaco:RemapperPoints 10000
	set Monaco:UseAllIterations No

	set Sampler:UpdateAfter 1000
	set Sampler:AddUpSamplers Off
	set Sampler:GlobalMaximumWeight Off
	set Sampler:FlatSubprocesses Off
	set Sampler:MinSelection 0.000001
	set Sampler:AlmostUnweighted On
	set Sampler:BinSampler:Kappa 1.
	set Sampler:RunCombinationData Off
	set Sampler:Verbose On

	#######
	# CMW scheme
	#######

	cd /Herwig/DipoleShower/Kernels
	set FFgx2ggxDipoleKernel:CMWScheme On
	set FFqx2qgxDipoleKernel:CMWScheme On
	set FFgx2ddxDipoleKernel:CMWScheme On
	set FFgx2uuxDipoleKernel:CMWScheme On
	set FFgx2ccxDipoleKernel:CMWScheme On
	set FFgx2ssxDipoleKernel:CMWScheme On
	set FFgx2bbxDipoleKernel:CMWScheme On
	set FFMgx2ggxDipoleKernel:CMWScheme On
	set FFMdx2dgxDipoleKernel:CMWScheme On
	set FFMux2ugxDipoleKernel:CMWScheme On
	set FFMcx2cgxDipoleKernel:CMWScheme On
	set FFMsx2sgxDipoleKernel:CMWScheme On
	set FFMbx2bgxDipoleKernel:CMWScheme On
	set FFMtx2tgxDipoleKernel:CMWScheme On
	set FFMgx2ddxDipoleKernel:CMWScheme On
	set FFMgx2uuxDipoleKernel:CMWScheme On
	set FFMgx2ccxDipoleKernel:CMWScheme On
	set FFMgx2ssxDipoleKernel:CMWScheme On
	set FFMgx2bbxDipoleKernel:CMWScheme On
	set FIgx2ggxDipoleKernel:CMWScheme On
	set FIqx2qgxDipoleKernel:CMWScheme On
	set FIgx2ddxDipoleKernel:CMWScheme On
	set FIgx2uuxDipoleKernel:CMWScheme On
	set FIgx2ccxDipoleKernel:CMWScheme On
	set FIgx2ssxDipoleKernel:CMWScheme On
	set FIgx2bbxDipoleKernel:CMWScheme On
	set FIMdx2dgxDipoleKernel:CMWScheme On
	set FIMux2ugxDipoleKernel:CMWScheme On
	set FIMcx2cgxDipoleKernel:CMWScheme On
	set FIMsx2sgxDipoleKernel:CMWScheme On
	set FIMbx2bgxDipoleKernel:CMWScheme On
	set FIMtx2tgxDipoleKernel:CMWScheme On
	set FIMgx2ddxDipoleKernel:CMWScheme On
	set FIMgx2uuxDipoleKernel:CMWScheme On
	set FIMgx2ccxDipoleKernel:CMWScheme On
	set FIMgx2ssxDipoleKernel:CMWScheme On
	set FIMgx2bbxDipoleKernel:CMWScheme On
	#set FIMgx2ttxDipoleKernel:CMWScheme On
	set IFgx2ggxDipoleKernel:CMWScheme On
	set IFqx2qgxDipoleKernel:CMWScheme On
	set IFqx2gqxDipoleKernel:CMWScheme On
	set IFgx2ddbarxDipoleKernel:CMWScheme On
	set IFgx2dbardxDipoleKernel:CMWScheme On
	set IFgx2uubarxDipoleKernel:CMWScheme On
	set IFgx2ubaruxDipoleKernel:CMWScheme On
	set IFgx2ccbarxDipoleKernel:CMWScheme On
	set IFgx2cbarcxDipoleKernel:CMWScheme On
	set IFgx2ssbarxDipoleKernel:CMWScheme On
	set IFgx2sbarsxDipoleKernel:CMWScheme On
	set IFMgx2ggxDipoleKernel:CMWScheme On
	set IFMqx2qgxDipoleKernel:CMWScheme On
	set IFMqx2gqxDipoleKernel:CMWScheme On
	set IFMgx2ddbarxDipoleKernel:CMWScheme On
	set IFMgx2dbardxDipoleKernel:CMWScheme On
	set IFMgx2uubarxDipoleKernel:CMWScheme On
	set IFMgx2ubaruxDipoleKernel:CMWScheme On
	set IFMgx2ccbarxDipoleKernel:CMWScheme On
	set IFMgx2cbarcxDipoleKernel:CMWScheme On
	set IFMgx2ssbarxDipoleKernel:CMWScheme On
	set IFMgx2sbarsxDipoleKernel:CMWScheme On
	set IIgx2ggxDipoleKernel:CMWScheme On
	set IIqx2qgxDipoleKernel:CMWScheme On
	set IIqx2gqxDipoleKernel:CMWScheme On
	set IIgx2ddbarxDipoleKernel:CMWScheme On
	set IIgx2dbardxDipoleKernel:CMWScheme On
	set IIgx2uubarxDipoleKernel:CMWScheme On
	set IIgx2ubaruxDipoleKernel:CMWScheme On
	set IIgx2ccbarxDipoleKernel:CMWScheme On
	set IIgx2cbarcxDipoleKernel:CMWScheme On
	set IIgx2ssbarxDipoleKernel:CMWScheme On
	set IIgx2sbarsxDipoleKernel:CMWScheme On
	diff --git a/src/Matchbox/MuDown.in b/src/Matchbox/MuDown.in
	--- a/src/Matchbox/MuDown.in
	+++ b/src/Matchbox/MuDown.in
	@@ -1,22 +1,22 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox

	set Factory:RenormalizationScaleFactor 0.5
	set Factory:FactorizationScaleFactor 0.5

	-set /Herwig/Merging/MFactory:RenormalizationScaleFactor 0.5
	-set /Herwig/Merging/MFactory:FactorizationScaleFactor 0.5
	+set /Herwig/Merging/MergingFactory:RenormalizationScaleFactor 0.5
	+set /Herwig/Merging/MergingFactory:FactorizationScaleFactor 0.5

	set MEMatching:RenormalizationScaleFactor 0.5
	set MEMatching:FactorizationScaleFactor 0.5

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 0.5
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 0.5

	set /Herwig/Shower/ShowerHandler:RenormalizationScaleFactor 0.5
	set /Herwig/Shower/ShowerHandler:FactorizationScaleFactor 0.5

	set /Herwig/Shower/PowhegShowerHandler:RenormalizationScaleFactor 0.5
	set /Herwig/Shower/PowhegShowerHandler:FactorizationScaleFactor 0.5

	diff --git a/src/Matchbox/MuUp.in b/src/Matchbox/MuUp.in
	--- a/src/Matchbox/MuUp.in
	+++ b/src/Matchbox/MuUp.in
	@@ -1,24 +1,24 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox

	set Factory:RenormalizationScaleFactor 2.0
	set Factory:FactorizationScaleFactor 2.0

	-set /Herwig/Merging/MFactory:RenormalizationScaleFactor 2.0
	-set /Herwig/Merging/MFactory:FactorizationScaleFactor 2.0
	+set /Herwig/Merging/MergingFactory:RenormalizationScaleFactor 2.0
	+set /Herwig/Merging/MergingFactory:FactorizationScaleFactor 2.0

	set MEMatching:RenormalizationScaleFactor 2.0
	set MEMatching:FactorizationScaleFactor 2.0

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 2.0
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 2.0

	set /Herwig/Shower/ShowerHandler:RenormalizationScaleFactor 2.0
	set /Herwig/Shower/ShowerHandler:FactorizationScaleFactor 2.0

	set /Herwig/Shower/PowhegShowerHandler:RenormalizationScaleFactor 2.0
	set /Herwig/Shower/PowhegShowerHandler:FactorizationScaleFactor 2.0



	diff --git a/src/Matchbox/NJet-NJet.in b/src/Matchbox/NJet-NJet.in
	--- a/src/Matchbox/NJet-NJet.in
	+++ b/src/Matchbox/NJet-NJet.in
	@@ -1,11 +1,11 @@
	# -- ThePEG-repository --

	# NJet as the one and only matrix element provider

	cd /Herwig/MatrixElements/Matchbox/Amplitudes

	clear /Herwig/MatrixElements/Matchbox/Factory:Amplitudes
	insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 NJet

	-clear /Herwig/Merging/MFactory:Amplitudes
	-insert /Herwig/Merging/MFactory:Amplitudes 0 NJets
	+clear /Herwig/Merging/MergingFactory:Amplitudes
	+insert /Herwig/Merging/MergingFactory:Amplitudes 0 NJets
	diff --git a/src/Matchbox/PPCollider.in b/src/Matchbox/PPCollider.in
	--- a/src/Matchbox/PPCollider.in
	+++ b/src/Matchbox/PPCollider.in
	@@ -1,30 +1,30 @@
	# -- ThePEG-repository --

	read Matchbox/GenericCollider.in

	cd /Herwig/EventHandlers

	create ThePEG::FixedCMSLuminosity Luminosity FixedCMSLuminosity.so
	set Luminosity:Energy 13000.0
	set EventHandler:LuminosityFunction Luminosity
	set EventHandler:BeamA /Herwig/Particles/p+
	set EventHandler:BeamB /Herwig/Particles/p+
	cd /Herwig/MatrixElements/Matchbox
	set Factory:FirstPerturbativePDF Yes
	set Factory:SecondPerturbativePDF Yes
	set Factory:PartonExtractor /Herwig/Partons/QCDExtractor

	cd /Herwig/Merging
	-set MFactory:FirstPerturbativePDF Yes
	-set MFactory:SecondPerturbativePDF Yes
	-set MFactory:PartonExtractor /Herwig/Partons/QCDExtractor
	+set MergingFactory:FirstPerturbativePDF Yes
	+set MergingFactory:SecondPerturbativePDF Yes
	+set MergingFactory:PartonExtractor /Herwig/Partons/QCDExtractor


	cd /Herwig/Cuts

	set Cuts:ScaleMin 2.0*GeV
	set Cuts:X1Min 1.0e-5
	set Cuts:X2Min 1.0e-5
	set Cuts:MHatMin 0.0*GeV


	diff --git a/src/Matchbox/StandardModelLike.in b/src/Matchbox/StandardModelLike.in
	--- a/src/Matchbox/StandardModelLike.in
	+++ b/src/Matchbox/StandardModelLike.in
	@@ -1,22 +1,22 @@
	# -- ThePEG-repository --

	set /Herwig/MatrixElements/Matchbox/Factory:EnforceChargeConservation Yes
	set /Herwig/MatrixElements/Matchbox/Factory:EnforceColourConservation Yes

	set /Herwig/MatrixElements/Matchbox/Factory:EnforceLeptonNumberConservation Yes
	set /Herwig/MatrixElements/Matchbox/Factory:EnforceQuarkNumberConservation Yes

	set /Herwig/MatrixElements/Matchbox/Factory:LeptonFlavourDiagonal Yes




	-set /Herwig/Merging/MFactory:EnforceChargeConservation Yes
	-set /Herwig/Merging/MFactory:EnforceColourConservation Yes
	+set /Herwig/Merging/MergingFactory:EnforceChargeConservation Yes
	+set /Herwig/Merging/MergingFactory:EnforceColourConservation Yes

	-set /Herwig/Merging/MFactory:EnforceLeptonNumberConservation Yes
	-set /Herwig/Merging/MFactory:EnforceQuarkNumberConservation Yes
	+set /Herwig/Merging/MergingFactory:EnforceLeptonNumberConservation Yes
	+set /Herwig/Merging/MergingFactory:EnforceQuarkNumberConservation Yes

	-set /Herwig/Merging/MFactory:LeptonFlavourDiagonal Yes
	+set /Herwig/Merging/MergingFactory:LeptonFlavourDiagonal Yes


	diff --git a/src/Matchbox/VBFNLO.in b/src/Matchbox/VBFNLO.in
	--- a/src/Matchbox/VBFNLO.in
	+++ b/src/Matchbox/VBFNLO.in
	@@ -1,17 +1,17 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox
	clear Factory:Amplitudes
	insert Factory:Amplitudes 0 /Herwig/MatrixElements/Matchbox/Amplitudes/VBFNLO


	-clear /Herwig/Merging/MFactory:Amplitudes
	-insert /Herwig/Merging/MFactory:Amplitudes 0 /Herwig/MatrixElements/Matchbox/Amplitudes/VBFNLO
	+clear /Herwig/Merging/MergingFactory:Amplitudes
	+insert /Herwig/Merging/MergingFactory:Amplitudes 0 /Herwig/MatrixElements/Matchbox/Amplitudes/VBFNLO

	cd /Herwig/Particles

	set b:HardProcessMass 0*GeV
	set bbar:HardProcessMass 0*GeV
	set tau+:HardProcessMass 0*GeV
	set tau-:HardProcessMass 0*GeV

	diff --git a/src/Matchbox/VBFNLOPhasespace.in b/src/Matchbox/VBFNLOPhasespace.in
	--- a/src/Matchbox/VBFNLOPhasespace.in
	+++ b/src/Matchbox/VBFNLOPhasespace.in
	@@ -1,9 +1,9 @@
	# -- ThePEG-repository --

	cd /Herwig/MatrixElements/Matchbox/Phasespace
	set /Herwig/MatrixElements/Matchbox/Factory:Phasespace VBFNLOPhasespace

	-insert /Herwig/Merging/MFactory:Phasespace VBFNLOPhasespace
	+insert /Herwig/Merging/MergingFactory:Phasespace VBFNLOPhasespace

	set /Herwig/MatrixElements/Matchbox/MECorrectionHandler:Phasespace InvertibleLabframePhasespace

	diff --git a/src/Merging/setupfiles/setupLEP1 b/src/Merging/setupfiles/setupLEP1
	--- a/src/Merging/setupfiles/setupLEP1
	+++ b/src/Merging/setupfiles/setupLEP1
	@@ -1,13 +1,13 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set EEMFactory:RenormalizationScaleFactor 1.
	-set EEMFactory:FactorizationScaleFactor 1.
	+set EEMergingFactory:RenormalizationScaleFactor 1.
	+set EEMergingFactory:FactorizationScaleFactor 1.


	set Merger:chooseHistory 8
	set Merger:MergingScale 3.5.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLEP10 b/src/Merging/setupfiles/setupLEP10
	--- a/src/Merging/setupfiles/setupLEP10
	+++ b/src/Merging/setupfiles/setupLEP10
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 2.
	-set EEMFactory:RenormalizationScaleFactor 1.
	-set EEMFactory:FactorizationScaleFactor 1.
	+set EEMergingFactory:RenormalizationScaleFactor 1.
	+set EEMergingFactory:FactorizationScaleFactor 1.

	-set EEMFactory:onlyabove 0
	+set EEMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 5.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLEP2 b/src/Merging/setupfiles/setupLEP2
	--- a/src/Merging/setupfiles/setupLEP2
	+++ b/src/Merging/setupfiles/setupLEP2
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set EEMFactory:RenormalizationScaleFactor 1.
	-set EEMFactory:FactorizationScaleFactor 1.
	+set EEMergingFactory:RenormalizationScaleFactor 1.
	+set EEMergingFactory:FactorizationScaleFactor 1.

	-set EEMFactory:onlyabove 0
	+set EEMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 5.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLEP3 b/src/Merging/setupfiles/setupLEP3
	--- a/src/Merging/setupfiles/setupLEP3
	+++ b/src/Merging/setupfiles/setupLEP3
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set EEMFactory:RenormalizationScaleFactor 2.
	-set EEMFactory:FactorizationScaleFactor 2.
	+set EEMergingFactory:RenormalizationScaleFactor 2.
	+set EEMergingFactory:FactorizationScaleFactor 2.

	-set EEMFactory:onlyabove 0
	+set EEMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 6.5.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLEP4 b/src/Merging/setupfiles/setupLEP4
	--- a/src/Merging/setupfiles/setupLEP4
	+++ b/src/Merging/setupfiles/setupLEP4
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set EEMFactory:RenormalizationScaleFactor 2.
	-set EEMFactory:FactorizationScaleFactor 2.
	+set EEMergingFactory:RenormalizationScaleFactor 2.
	+set EEMergingFactory:FactorizationScaleFactor 2.

	-set EEMFactory:onlyabove 0
	+set EEMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 5.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLEP5 b/src/Merging/setupfiles/setupLEP5
	--- a/src/Merging/setupfiles/setupLEP5
	+++ b/src/Merging/setupfiles/setupLEP5
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set EEMFactory:RenormalizationScaleFactor 0.5
	-set EEMFactory:FactorizationScaleFactor 0.5
	+set EEMergingFactory:RenormalizationScaleFactor 0.5
	+set EEMergingFactory:FactorizationScaleFactor 0.5

	-set EEMFactory:onlyabove 0
	+set EEMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 5.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLEP6 b/src/Merging/setupfiles/setupLEP6
	--- a/src/Merging/setupfiles/setupLEP6
	+++ b/src/Merging/setupfiles/setupLEP6
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 2.
	-set EEMFactory:RenormalizationScaleFactor 1.
	-set EEMFactory:FactorizationScaleFactor 1.
	+set EEMergingFactory:RenormalizationScaleFactor 1.
	+set EEMergingFactory:FactorizationScaleFactor 1.

	-set EEMFactory:onlyabove 0
	+set EEMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 5.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLEP7 b/src/Merging/setupfiles/setupLEP7
	--- a/src/Merging/setupfiles/setupLEP7
	+++ b/src/Merging/setupfiles/setupLEP7
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 0.5
	-set EEMFactory:RenormalizationScaleFactor 1.
	-set EEMFactory:FactorizationScaleFactor 1.
	+set EEMergingFactory:RenormalizationScaleFactor 1.
	+set EEMergingFactory:FactorizationScaleFactor 1.

	-set EEMFactory:onlyabove 0
	+set EEMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 5.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLEP8 b/src/Merging/setupfiles/setupLEP8
	--- a/src/Merging/setupfiles/setupLEP8
	+++ b/src/Merging/setupfiles/setupLEP8
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 2.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 2.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set EEMFactory:RenormalizationScaleFactor 1.
	-set EEMFactory:FactorizationScaleFactor 1.
	+set EEMergingFactory:RenormalizationScaleFactor 1.
	+set EEMergingFactory:FactorizationScaleFactor 1.

	-set EEMFactory:onlyabove 0
	+set EEMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 5.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLEP9 b/src/Merging/setupfiles/setupLEP9
	--- a/src/Merging/setupfiles/setupLEP9
	+++ b/src/Merging/setupfiles/setupLEP9
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 0.5
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 0.5
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set EEMFactory:RenormalizationScaleFactor 1.
	-set EEMFactory:FactorizationScaleFactor 1.
	+set EEMergingFactory:RenormalizationScaleFactor 1.
	+set EEMergingFactory:FactorizationScaleFactor 1.

	-set EEMFactory:onlyabove 0
	+set EEMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 5.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLHC1 b/src/Merging/setupfiles/setupLHC1
	--- a/src/Merging/setupfiles/setupLHC1
	+++ b/src/Merging/setupfiles/setupLHC1
	@@ -1,13 +1,13 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set PPMFactory:RenormalizationScaleFactor 1.
	-set PPMFactory:FactorizationScaleFactor 1.
	+set PPMergingFactory:RenormalizationScaleFactor 1.
	+set PPMergingFactory:FactorizationScaleFactor 1.


	set Merger:chooseHistory 8
	set Merger:MergingScale 10.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLHC10 b/src/Merging/setupfiles/setupLHC10
	--- a/src/Merging/setupfiles/setupLHC10
	+++ b/src/Merging/setupfiles/setupLHC10
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 2.
	-set PPMFactory:RenormalizationScaleFactor 1.
	-set PPMFactory:FactorizationScaleFactor 1.
	+set PPMergingFactory:RenormalizationScaleFactor 1.
	+set PPMergingFactory:FactorizationScaleFactor 1.

	-set PPMFactory:onlyabove 0
	+set PPMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 20.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLHC2 b/src/Merging/setupfiles/setupLHC2
	--- a/src/Merging/setupfiles/setupLHC2
	+++ b/src/Merging/setupfiles/setupLHC2
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set PPMFactory:RenormalizationScaleFactor 1.
	-set PPMFactory:FactorizationScaleFactor 1.
	+set PPMergingFactory:RenormalizationScaleFactor 1.
	+set PPMergingFactory:FactorizationScaleFactor 1.

	-set PPMFactory:onlyabove 0
	+set PPMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 20.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLHC3 b/src/Merging/setupfiles/setupLHC3
	--- a/src/Merging/setupfiles/setupLHC3
	+++ b/src/Merging/setupfiles/setupLHC3
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set PPMFactory:RenormalizationScaleFactor 2.
	-set PPMFactory:FactorizationScaleFactor 2.
	+set PPMergingFactory:RenormalizationScaleFactor 2.
	+set PPMergingFactory:FactorizationScaleFactor 2.

	-set PPMFactory:onlyabove 0
	+set PPMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 30.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLHC4 b/src/Merging/setupfiles/setupLHC4
	--- a/src/Merging/setupfiles/setupLHC4
	+++ b/src/Merging/setupfiles/setupLHC4
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set PPMFactory:RenormalizationScaleFactor 2.
	-set PPMFactory:FactorizationScaleFactor 2.
	+set PPMergingFactory:RenormalizationScaleFactor 2.
	+set PPMergingFactory:FactorizationScaleFactor 2.

	-set PPMFactory:onlyabove 0
	+set PPMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 20.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLHC5 b/src/Merging/setupfiles/setupLHC5
	--- a/src/Merging/setupfiles/setupLHC5
	+++ b/src/Merging/setupfiles/setupLHC5
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set PPMFactory:RenormalizationScaleFactor 0.5
	-set PPMFactory:FactorizationScaleFactor 0.5
	+set PPMergingFactory:RenormalizationScaleFactor 0.5
	+set PPMergingFactory:FactorizationScaleFactor 0.5

	-set PPMFactory:onlyabove 0
	+set PPMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 20.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLHC6 b/src/Merging/setupfiles/setupLHC6
	--- a/src/Merging/setupfiles/setupLHC6
	+++ b/src/Merging/setupfiles/setupLHC6
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 2.
	-set PPMFactory:RenormalizationScaleFactor 1.
	-set PPMFactory:FactorizationScaleFactor 1.
	+set PPMergingFactory:RenormalizationScaleFactor 1.
	+set PPMergingFactory:FactorizationScaleFactor 1.

	-set PPMFactory:onlyabove 0
	+set PPMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 20.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLHC7 b/src/Merging/setupfiles/setupLHC7
	--- a/src/Merging/setupfiles/setupLHC7
	+++ b/src/Merging/setupfiles/setupLHC7
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 1.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 0.5
	-set PPMFactory:RenormalizationScaleFactor 1.
	-set PPMFactory:FactorizationScaleFactor 1.
	+set PPMergingFactory:RenormalizationScaleFactor 1.
	+set PPMergingFactory:FactorizationScaleFactor 1.

	-set PPMFactory:onlyabove 0
	+set PPMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 20.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLHC8 b/src/Merging/setupfiles/setupLHC8
	--- a/src/Merging/setupfiles/setupLHC8
	+++ b/src/Merging/setupfiles/setupLHC8
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 2.
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 2.
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set PPMFactory:RenormalizationScaleFactor 1.
	-set PPMFactory:FactorizationScaleFactor 1.
	+set PPMergingFactory:RenormalizationScaleFactor 1.
	+set PPMergingFactory:FactorizationScaleFactor 1.

	-set PPMFactory:onlyabove 0
	+set PPMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 20.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.
	diff --git a/src/Merging/setupfiles/setupLHC9 b/src/Merging/setupfiles/setupLHC9
	--- a/src/Merging/setupfiles/setupLHC9
	+++ b/src/Merging/setupfiles/setupLHC9
	@@ -1,14 +1,14 @@
	cd /Herwig/Merging

	set /Herwig/DipoleShower/DipoleShowerHandler:RenormalizationScaleFactor 0.5
	set /Herwig/DipoleShower/DipoleShowerHandler:FactorizationScaleFactor 0.5
	set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleFactor 1.
	-set PPMFactory:RenormalizationScaleFactor 1.
	-set PPMFactory:FactorizationScaleFactor 1.
	+set PPMergingFactory:RenormalizationScaleFactor 1.
	+set PPMergingFactory:FactorizationScaleFactor 1.

	-set PPMFactory:onlyabove 0
	+set PPMergingFactory:onlyabove 0

	set Merger:chooseHistory 8
	set Merger:MergingScale 20.*GeV
	set Merger:MergingScaleSmearing 0.1
	set Merger:gamma 1.

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 4:30 PM (21 h, 22 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023451
Default Alt Text: (214 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions