No OneTemporary
Actions

Size

114 KB

Subscribers

None

View Options

	diff --git a/DipoleShower/Kinematics/FFMassiveKinematics.cc b/DipoleShower/Kinematics/FFMassiveKinematics.cc
	--- a/DipoleShower/Kinematics/FFMassiveKinematics.cc
	+++ b/DipoleShower/Kinematics/FFMassiveKinematics.cc
	@@ -1,350 +1,343 @@
	// -- C++ --
	//
	// FFMassiveKinematics.cc 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.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the FFMassiveKinematics class.
	//

	#include "FFMassiveKinematics.h"
	#include "ThePEG/Interface/ClassDocumentation.h"

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
	#include "Herwig++/DipoleShower/Kernels/DipoleSplittingKernel.h"

	// TODO: remove after verification
	// only for checking for NaN or inf
	#include <gsl/gsl_math.h>

	using namespace Herwig;

	FFMassiveKinematics::FFMassiveKinematics()
	: DipoleSplittingKinematics() {}

	FFMassiveKinematics::~FFMassiveKinematics() {}

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

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

	pair<double,double> FFMassiveKinematics::kappaSupport(const DipoleSplittingInfo&) const {
	return make_pair(0.0,1.0);
	}

	pair<double,double> FFMassiveKinematics::xiSupport(const DipoleSplittingInfo& split) const {
	double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
	if ( split.index().emitterData()->id() == ParticleID::g ) {
	if ( split.emissionData()->id() != ParticleID::g )
	return make_pair(0.5(1.-c),0.5(1.+c));
	double b = log((1.+c)/(1.-c));
	return make_pair(-b,b);
	}
	return make_pair(-log(0.5(1.+c)),-log(0.5(1.-c)));
	}

	Energy FFMassiveKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const {
	return (pEmitter+pSpectator).m();
	}

	Energy FFMassiveKinematics::ptMax(Energy dScale,
	double, double,
	const DipoleIndex& ind,
	const DipoleSplittingKernel& split) const {
	double mui = split.emitter(ind)->mass() / dScale;
	double mu = split.emission(ind)->mass() / dScale;
	double muj = split.spectator(ind)->mass() / dScale;
	double mui2 = sqr( mui ), mu2 = sqr( mu ), muj2 = sqr( muj );
	-
	- if(mui>(1.+mu-muj)&&
	- ( mui2mui2 + mu2mu2 + sqr(1.-muj)sqr(1.-muj) - 2.( mui2mu2+mui2sqr(1.-muj)+mu2*sqr(1.-muj) ))>0
	- ) {
	- return rootOfKallen( mui2, mu2, sqr(1.-muj) ) / ( 2.-2.sqrt(muj2) ) dScale;
	- }
	- else
	- return ZERO;
	+ return rootOfKallen( mui2, mu2, sqr(1.-muj) ) / ( 2.-2.sqrt(muj2) ) dScale;
	}

	Energy FFMassiveKinematics::QMax(Energy dScale,
	double, double,
	const DipoleIndex& ind,
	const DipoleSplittingKernel&) const {
	assert(false && "implementation missing");
	double Muj = ind.spectatorData()->mass() / dScale;
	return dScale * ( 1.-2.*Muj+sqr(Muj) );
	}

	Energy FFMassiveKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
	// from Martin's thesis
	double z = split.lastZ();
	Energy mi = split.emitterData()->mass();
	Energy m = split.emissionData()->mass();
	Energy2 pt2 = z(1.-z)sqr(scale) - (1-z)sqr(mi) - zsqr(m);
	assert(pt2 >= ZERO);
	return sqrt(pt2);
	}

	Energy FFMassiveKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
	// from Martin's thesis
	double z = split.lastZ();
	Energy mi = split.emitterData()->mass();
	Energy m = split.emissionData()->mass();
	Energy2 Q2 = (sqr(scale) + (1-z)sqr(mi) + zsqr(m))/(z*(1.-z));
	return sqrt(Q2);
	}

	double FFMassiveKinematics::ptToRandom(Energy pt, Energy,
	double,double,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const {
	return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
	}

	// own kinematics close to Dinsdale,Ternick,Weinzierl
	bool FFMassiveKinematics::generateSplitting(double kappa, double xi, double rphi,
	DipoleSplittingInfo& info,
	const DipoleSplittingKernel&) {

	Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);

	if ( pt > info.hardPt() \|\| pt < IRCutoff() ) {
	jacobian(0.0);
	return false;
	}

	double z;
	double mapZJacobian;

	if ( info.index().emitterData()->id() == ParticleID::g ) {
	if ( info.emissionData()->id() != ParticleID::g ) {
	z = xi;
	mapZJacobian = 1.;
	} else {
	z = exp(xi)/(1.+exp(xi));
	mapZJacobian = z*(1.-z);
	}
	} else {
	z = 1.-exp(-xi);
	mapZJacobian = 1.-z;
	}

	// masses
	double mui2 = sqr( info.emitterData()->mass() / info.scale() );
	double mu2 = sqr( info.emissionData()->mass() / info.scale() );
	double muj2 = sqr( info.spectatorData()->mass() / info.scale() );
	double Mui2 = 0.;
	if ( info.emitterData()->id() + info.emissionData()->id() == 0 ) Mui2 = 0.; // gluon
	else Mui2 = mui2; // (anti)quark
	double Muj2 = muj2;

	// new: 2011-08-31
	// 2011-11-08: this does happen
	if( sqrt(mui2)+sqrt(mu2)+sqrt(muj2) > 1. ){
	jacobian(0.0);
	return false;
	}

	double bar = 1.-mui2-mu2-muj2;
	double y = ( sqr( pt / info.scale() ) + sqr(1.-z)mui2 + zz*mu2 ) /
	(z(1.-z)bar);

	// do this here for simplicity
	Energy ptmax1 = rootOfKallen( mui2, mu2, sqr(1.-sqrt(muj2)) ) /
	( 2.-2.sqrt(muj2) ) info.scale();
	Energy auxHardPt = ptmax1 > info.hardPt() ? info.hardPt() : ptmax1;
	// 2011-11-09
	assert(ptmax1>=info.hardPt());

	// phasespace constraint to incorporate ptMax
	double zp1 = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) +
	rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
	sqrt( 1.-sqr(pt/auxHardPt) ) ) /
	( 2.*sqr(1.-sqrt(muj2)) );
	double zm1 = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) -
	rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
	sqrt( 1.-sqr(pt/auxHardPt) ) ) /
	( 2.*sqr(1.-sqrt(muj2)) );
	if ( z > zp1 \|\| z < zm1 \|\|
	pt > auxHardPt) {
	jacobian(0.0);
	return false;
	}

	// kinematic phasespace boundaries for (y,z)
	// same as in Dittmaier hep-ph/9904440v2 (equivalent to CS)
	double ym = 2.sqrt(mui2)sqrt(mu2)/bar;
	double yp = 1. - 2.sqrt(muj2)(1.-sqrt(muj2))/bar;
	if ( y < ym \|\| y > yp ) {
	jacobian(0.0);
	return false;
	}

	double zm = ( (2.mui2+bary)(1.-y) - sqrt(yy-ymym)sqrt(sqr(2.muj2+bar-bary)-4.*muj2) ) /
	( 2.(1.-y)(mui2+mu2+bar*y) );
	double zp = ( (2.mui2+bary)(1.-y) + sqrt(yy-ymym)sqrt(sqr(2.muj2+bar-bary)-4.*muj2) ) /
	( 2.(1.-y)(mui2+mu2+bar*y) );

	if ( z < zm \|\| z > zp ) {
	jacobian(0.0);
	return false;
	}

	double phi = 2.Constants::pirphi;

	jacobian( 2. * mapZJacobian * (1.-y) *
	log(0.5 * generator()->maximumCMEnergy()/IRCutoff()) *
	bar / rootOfKallen(1.,Mui2,Muj2) );

	lastPt(pt);
	lastZ(z);
	lastPhi(phi);

	if ( theMCCheck )
	theMCCheck->book(1.,1.,info.scale(),info.hardPt(),pt,z,jacobian());

	return true;

	}

	// kinematics close to Dinsdale,Ternick,Weinzierl
	// revised 2011-08-22
	// revised 2011-11-06
	void FFMassiveKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) {

	double z = dInfo.lastZ();
	Energy pt = dInfo.lastPt();

	// masses
	double mui2 = sqr( dInfo.emitterData()->mass() / dInfo.scale() );
	double mu2 = sqr( dInfo.emissionData()->mass() / dInfo.scale() );
	double muj2 = sqr( dInfo.spectatorData()->mass() / dInfo.scale() );

	double y = ( sqr( pt / dInfo.scale() ) + sqr(1.-z)mui2 + zz*mu2 ) /
	(z(1.-z)(1.-mui2-mu2-muj2));

	Energy2 sbar = sqr(dInfo.scale()) *(1.-mui2-mu2-muj2);

	// CMF: particle energies
	Energy Ei = ( sbar(1.-(1.-z)(1.-y)) + 2.sqr(dInfo.scale())mui2 ) / (2.*dInfo.scale());
	Energy E = ( sbar(1.- z (1.-y)) + 2.sqr(dInfo.scale())mu2 ) / (2.*dInfo.scale());
	Energy Ej = ( sbar(1.- y ) + 2.sqr(dInfo.scale())muj2 ) / (2.dInfo.scale());
	// CMF: momenta in z-direction (axis of pEmitter & pSpectator)
	Energy qi3 = (2.EiEj-z(1.-y)sbar ) / 2./sqrt(EjEj-sqr(dInfo.scale())muj2);
	Energy q3 = (2.E Ej-(1.-z)(1.-y)sbar) / 2./sqrt(EjEj-sqr(dInfo.scale())muj2);
	Energy qj3 = sqrt( sqr(Ej) - sqr(dInfo.scale())*muj2 );

	// get z axis in the dipole's CMF which is parallel to pSpectator
	Boost toCMF = (pEmitter+pSpectator).findBoostToCM();
	Lorentz5Momentum pjAux = pSpectator; pjAux.boost(toCMF);
	ThreeVector<double> pjAxis = pjAux.vect().unit();

	// set the momenta in this special reference frame
	// note that pt might in some cases differ from the physical pt!
	// phi is defined exactly as in getKt
	Energy ptResc = sqrt( sqr(Ei)-sqr(dInfo.scale())*mui2-sqr(qi3) );
	Lorentz5Momentum em ( ptResccos(dInfo.lastPhi()), -ptRescsin(dInfo.lastPhi()), qi3, Ei );
	Lorentz5Momentum emm ( -ptResccos(dInfo.lastPhi()), ptRescsin(dInfo.lastPhi()), q3, E );
	Lorentz5Momentum spe ( 0.GeV, 0.GeV, qj3, Ej );

	// output the mismatch between pt and physical pt
	// ofstream output1("ptDiffOnPtAxis-uub-m.dat",ofstream::app);
	// ofstream output2("ptDiffOnCosAxis-uub-m.dat",ofstream::app);
	// if( abs(dInfo.spectatorData()->id())==5 && dInfo.emitterData()->id()+dInfo.emissionData()->id()==0 &&
	// abs(dInfo.emitterData()->id())==1 ) {
	// output1 << pt/dInfo.scale() << " " << abs(ptResc-pt)/(ptResc+pt) << " " << endl;
	// output2 << em.vect().unit()*emm.vect().unit() << " " << abs(ptResc-pt)/(ptResc+pt) << " " << endl;
	// }
	// output1.close(); output2.close();

	// rotate back
	em.rotateUz (pjAxis);
	emm.rotateUz(pjAxis);
	spe.rotateUz(pjAxis);

	// boost back
	em.boost (-toCMF);
	emm.boost(-toCMF);
	spe.boost(-toCMF);

	// mass shells, rescale energy
	em.setMass(dInfo.scale()*sqrt(mui2));
	em.rescaleEnergy();
	emm.setMass(dInfo.scale()*sqrt(mu2));
	emm.rescaleEnergy();
	spe.setMass(dInfo.scale()*sqrt(muj2));
	spe.rescaleEnergy();

	// book
	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	// TODO: remove
	// 2011-11-09: never occurred
	if(em.t()/GeV>=0. && emm.t()/GeV>=0. && spe.t()/GeV>=0.);
	else cout << "FFMassiveKinematics::generateKinematics momenta corrupt" << endl;

	// 2011-11-03 LEP run with full masses:
	// x,y,t no problem
	// z order > 5.e-7 happend 41 times in 10000 runs
	// maximum was 2.5e-6

	// double order=2.4e-6;
	// Lorentz5Momentum pDif=em+emm+spe-pEmitter-pSpectator, pSum=em+emm+spe+pEmitter+pSpectator;
	// if(abs(pDif.x()/(pSum.x()==ZERO ? 1.*GeV : pSum.x())) <order \|\| (pDif-pSum).x()==ZERO);
	// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: x " <<
	// abs(pDif.x()/(pSum.x()==ZERO ? 1.*GeV : pSum.x())) << endl <<
	// " " << em/GeV << " " << emm/GeV << " " << spe/GeV << endl <<
	// " " << pEmitter/GeV << " " << pSpectator/GeV << endl;
	// if(abs(pDif.y()/(pSum.y()==ZERO ? 1.*GeV : pSum.y())) <order \|\| (pDif-pSum).y()==ZERO);
	// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: y " <<
	// abs(pDif.y()/(pSum.y()==ZERO ? 1.*GeV : pSum.y())) << endl;
	// if(abs(pDif.z()/(pSum.z()==ZERO ? 1.*GeV : pSum.z())) <order \| (pDif-pSum).z()==ZERO);
	// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: z " <<
	// abs(pDif.z()/(pSum.z()==ZERO ? 1.*GeV : pSum.z())) << endl <<
	// " " << em/GeV << " " << emm/GeV << " " << spe/GeV << endl <<
	// " " << pEmitter/GeV << " " << pSpectator/GeV << endl;
	// if(abs(pDif.t()/(pSum.t()==ZERO ? 1.*GeV : pSum.t())) <order \|\| (pDif-pSum).t()==ZERO);
	// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: t " <<
	// // abs(pDif.t()/(pSum.t()==ZERO ? 1.*GeV : pSum.t())) << endl;
	// cout << endl;

	}

	// If needed, insert default implementations of function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void FFMassiveKinematics::persistentOutput(PersistentOStream & ) const {
	}

	void FFMassiveKinematics::persistentInput(PersistentIStream & , int) {
	}

	ClassDescription<FFMassiveKinematics> FFMassiveKinematics::initFFMassiveKinematics;
	// Definition of the static class description member.

	void FFMassiveKinematics::Init() {

	static ClassDocumentation<FFMassiveKinematics> documentation
	("FFMassiveKinematics implements massive splittings "
	"off a final-final dipole.");

	}

	diff --git a/DipoleShower/Kinematics/FFMassiveKinematics.h b/DipoleShower/Kinematics/FFMassiveKinematics.h
	--- a/DipoleShower/Kinematics/FFMassiveKinematics.h
	+++ b/DipoleShower/Kinematics/FFMassiveKinematics.h
	@@ -1,261 +1,262 @@
	// -- C++ --
	//
	// FFMassiveKinematics.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.
	//
	#ifndef HERWIG_FFMassiveKinematics_H
	#define HERWIG_FFMassiveKinematics_H
	//
	// This is the declaration of the FFMassiveKinematics class.
	//

	#include "DipoleSplittingKinematics.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer
	*
	* \brief FFMassiveKinematics implements massive splittings
	* off a final-final dipole.
	*
	*/
	class FFMassiveKinematics: public DipoleSplittingKinematics {

	public:

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

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

	public:

	/**
	* Return the boundaries in between the evolution
	* variable random number is to be sampled; the lower
	* cuoff is assumed to correspond to the infrared cutoff.
	*/
	virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;

	/**
	* Return the boundaries in between the momentum
	* fraction random number is to be sampled.
	*/
	virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;

	/**
	* Return the boundaries on the momentum fraction
	*/
	virtual pair<double,double> zBoundaries(Energy,
	const DipoleSplittingInfo&,
	const DipoleSplittingKernel&) const {
	return make_pair(0.0,1.0);
	}

	/**
	* Return the dipole scale associated to the
	* given pair of emitter and spectator. This
	* should be the invariant mass or absolute value
	* final/final or initial/initial and the absolute
	* value of the momentum transfer for intial/final or
	* final/initial dipoles.
	*/
	virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const;

	/**
	* Return the maximum pt for the given dipole scale.
	*/
	virtual Energy ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex& ind,
	const DipoleSplittingKernel& split) const;

	/**
	* Return the maximum virtuality for the given dipole scale.
	*/
	virtual Energy QMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex& dIndex,
	const DipoleSplittingKernel& split) const;

	/**
	* Return the pt given a virtuality.
	*/
	virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the virtuality given a pt.
	*/
	virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the random number associated to
	* the given pt.
	*/
	virtual double ptToRandom(Energy pt, Energy dScale,
	double emX, double specX,
	const DipoleIndex& dIndex,
	const DipoleSplittingKernel& split) const;

	/**
	* Generate splitting variables given three random numbers
	* and the momentum fractions of the emitter and spectator.
	* Return true on success.
	*/
	virtual bool generateSplitting(double kappa, double xi, double phi,
	DipoleSplittingInfo& dIndex,
	const DipoleSplittingKernel& split);

	/**
	* Generate the full kinematics given emitter and
	* spectator momentum and a previously completeted
	* DipoleSplittingInfo object.
	*/
	virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo);

	public:

	/**
	* Triangular / Kallen function
	*/
	template <class T>
	inline double rootOfKallen (T a, T b, T c) const {
	- return sqrt( aa + bb + cc - 2.( ab+ac+b*c ) ); }
	+ double sres=aa + bb + cc - 2.( ab+ac+b*c );
	+ return sres>0.?sqrt( sres ):0.; }

	/**
	* Perform a rotation on both momenta such that the first one will
	* point along the (positive) z axis. Rotate back to the original
	* reference frame by applying rotateUz(returnedVector) to each momentum.
	*/
	ThreeVector<double> rotateToZ (Lorentz5Momentum& pTarget, Lorentz5Momentum& p1){
	ThreeVector<double> oldAxis = pTarget.vect().unit();
	double ct = oldAxis.z(); double st = sqrt( 1.-sqr(ct) ); // cos,sin(theta)
	double cp = oldAxis.x()/st; double sp = oldAxis.y()/st; // cos,sin(phi)
	pTarget.setZ( pTarget.vect().mag() ); pTarget.setX( 0.GeV ); pTarget.setY( 0.GeV );
	Lorentz5Momentum p1old = p1;
	p1.setX( spp1old.x() - cpp1old.y() );
	p1.setY( ctcpp1old.x() + ctspp1old.y() - st*p1old.z() );
	p1.setZ( stcpp1old.x() + stspp1old.y() + ct*p1old.z() );
	return oldAxis;
	}

	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 static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<FFMassiveKinematics> initFFMassiveKinematics;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the FFMassiveKinematics class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::FFMassiveKinematics>
	: public ClassTraitsBase<Herwig::FFMassiveKinematics> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::FFMassiveKinematics"; }
	/**
	* The name of a file containing the dynamic library where the class
	* FFMassiveKinematics is implemented. It may also include several, space-separated,
	* libraries if the class FFMassiveKinematics depends on other classes (base classes
	* excepted). In this case the listed libraries will be dynamically
	* linked in the order they are specified.
	*/
	static string library() { return "HwDipoleShower.so"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_FFMassiveKinematics_H */
	diff --git a/DipoleShower/Kinematics/FIMassiveKinematics.cc b/DipoleShower/Kinematics/FIMassiveKinematics.cc
	--- a/DipoleShower/Kinematics/FIMassiveKinematics.cc
	+++ b/DipoleShower/Kinematics/FIMassiveKinematics.cc
	@@ -1,282 +1,279 @@
	// -- C++ --
	//
	// FIMassiveKinematics.cc 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.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the FIMassiveKinematics class.
	//

	#include "FIMassiveKinematics.h"
	#include "ThePEG/Interface/ClassDocumentation.h"


	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
	#include "Herwig++/DipoleShower/Kernels/DipoleSplittingKernel.h"

	using namespace Herwig;

	FIMassiveKinematics::FIMassiveKinematics()
	: DipoleSplittingKinematics() {}

	FIMassiveKinematics::~FIMassiveKinematics() {}

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

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

	pair<double,double> FIMassiveKinematics::kappaSupport(const DipoleSplittingInfo&) const {
	return make_pair(0.0,1.0);
	}

	pair<double,double> FIMassiveKinematics::xiSupport(const DipoleSplittingInfo& split) const {
	double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
	if ( split.index().emitterData()->id() == ParticleID::g ) {
	if ( split.emissionData()->id() != ParticleID::g )
	return make_pair(0.5(1.-c),0.5(1.+c));
	double b = log((1.+c)/(1.-c));
	return make_pair(-b,b);
	}
	return make_pair(-log(0.5(1.+c)),-log(0.5(1.-c)));
	}

	// sbar
	Energy FIMassiveKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const {
	return sqrt(2.(pEmitterpSpectator));
	}

	Energy FIMassiveKinematics::ptMax(Energy dScale,
	double, double specX,
	const DipoleIndex& ind,
	const DipoleSplittingKernel& split) const {
	Energy mi = split.emitter(ind)->mass(), m = split.emission(ind)->mass();
	Energy2 mi2 = sqr(mi), m2 = sqr(m);
	// Energy2 Mi2 = split.emitter(int)->id() + split.emission(int)->id() == 0 ?
	// 0.*GeV2 : mi2;
	Energy2 Mi2 = mi2 == m2 ? 0.*GeV2 : mi2;

	// s^star/x
	Energy2 s = sqr(dScale) * (1.-specX)/specX + Mi2;
	Energy roots = sqrt(s);
	- if(roots>mi+m)
	- return .5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
	- else
	- return ZERO;
	+ return .5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
	}

	// what is this? in FF it is given by y+dScale = sqrt( 2qiq / bar )->max
	Energy FIMassiveKinematics::QMax(Energy dScale,
	double, double specX,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const {
	assert(false && "implementation missing");
	cout << "FIMassiveKinematics::QMax called.\n" << flush;
	// this is sqrt( 2qi*q ) -> max;
	return dScale * sqrt((1.-specX)/specX);
	}

	Energy FIMassiveKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
	// from Martin's thesis
	double z = split.lastZ();
	Energy mi = split.emitterData()->mass();
	Energy m = split.emissionData()->mass();
	Energy2 pt2 = z(1.-z)sqr(scale) - (1-z)sqr(mi) - zsqr(m);
	assert(pt2 >= ZERO);
	return sqrt(pt2);
	}

	Energy FIMassiveKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
	// from Martin's thesis
	double z = split.lastZ();
	Energy mi = split.emitterData()->mass();
	Energy m = split.emissionData()->mass();
	Energy2 Q2 = (sqr(scale) + (1-z)sqr(mi) + zsqr(m))/(z*(1.-z));
	return sqrt(Q2);
	}

	double FIMassiveKinematics::ptToRandom(Energy pt, Energy,
	double,double,
	const DipoleIndex&,
	const DipoleSplittingKernel&) const {
	return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
	}

	bool FIMassiveKinematics::generateSplitting(double kappa, double xi, double rphi,
	DipoleSplittingInfo& info,
	const DipoleSplittingKernel&) {

	if ( info.spectatorX() < xMin() ) {
	jacobian(0.0);
	return false;
	}

	Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);

	if ( pt > info.hardPt() \|\| pt < IRCutoff() ) {
	jacobian(0.0);
	return false;
	}

	double z;
	double mapZJacobian;

	if ( info.index().emitterData()->id() == ParticleID::g ) {
	if ( info.emissionData()->id() != ParticleID::g ) {
	z = xi;
	mapZJacobian = 1.;
	} else {
	z = exp(xi)/(1.+exp(xi));
	mapZJacobian = z*(1.-z);
	}
	} else {
	z = 1.-exp(-xi);
	mapZJacobian = 1.-z;
	}

	// double s = z*(1.-z);

	// double xs = info.spectatorX();

	// double x = 1. / ( 1. + sqr(pt/info.scale()) / s );

	// double zp = 0.5*(1.+sqrt(1.-sqr(pt/info.hardPt())));
	// double zm = 0.5*(1.-sqrt(1.-sqr(pt/info.hardPt())));

	Energy2 mi2 = sqr(info.emitterData()->mass());
	Energy2 m2 = sqr(info.emissionData()->mass());
	Energy2 Mi2 = info.emitterData()->id()+info.emissionData()->id() == 0 ?
	0.*GeV2 : mi2;

	// s^star/x
	Energy2 s = sqr(info.scale()) * (1.-info.spectatorX())/info.spectatorX() + Mi2;

	double xs = info.spectatorX();
	double x = 1. / ( 1. +
	( sqr(pt) + (1.-z)mi2 + zm2 - z(1.-z)Mi2 ) /
	( z(1.-z)s ) );

	double zm1 = .5( 1.+(mi2-m2)/s - rootOfKallen(s/s,mi2/s,m2/s)
	sqrt( 1.-sqr(pt/info.hardPt()) ) );
	double zp1 = .5( 1.+(mi2-m2)/s + rootOfKallen(s/s,mi2/s,m2/s)
	sqrt( 1.-sqr(pt/info.hardPt()) ) );

	if ( // pt < IRCutoff() \|\|
	// pt > info.hardPt() \|\|
	z > zp1 \|\| z < zm1 \|\|
	x < xs ) {
	jacobian(0.0);
	return false;
	}

	// additional purely kinematic constraints
	double mui2 = x*mi2/sqr(info.scale());
	double mu2 = x*m2/sqr(info.scale());
	double Mui2 = x*Mi2/sqr(info.scale());
	double xp = 1. + Mui2 - sqr(sqrt(mui2)+sqrt(mu2));
	double root = sqr(1.-x+Mui2-mui2-mu2)-4.mui2mu2;
	if( root < 0. && root>-1e-10 )
	root = 0.;
	else if (root <0. ) {
	jacobian(0.0);
	return false;
	}

	root = sqrt(root);
	double zm = .5*( 1.-x+Mui2+mui2-mui2 - root ) / (1.-x+Mui2);
	double zp = .5*( 1.-x+Mui2+mui2-mui2 + root ) / (1.-x+Mui2);
	if (x > xp \|\|
	z > zp \|\| z < zm ) {
	jacobian(0.0);
	return false;
	}

	double phi = 2.Constants::pirphi;

	jacobian( 2. * mapZJacobian * log(0.5 * generator()->maximumCMEnergy()/IRCutoff()));

	lastPt(pt);
	lastZ(z);
	lastPhi(phi);
	lastSpectatorZ(x);

	if ( theMCCheck )
	theMCCheck->book(1.,info.spectatorX(),info.scale(),info.hardPt(),pt,z,jacobian());

	return true;

	}

	void FIMassiveKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo) {

	Energy pt = dInfo.lastPt();
	double z = dInfo.lastZ();

	Lorentz5Momentum kt =
	getKt (pSpectator, pEmitter, pt, dInfo.lastPhi(),true);

	Energy2 mi2 = sqr(dInfo.emitterData()->mass());
	Energy2 m2 = sqr(dInfo.emissionData()->mass());
	Energy2 Mi2 = dInfo.emitterData()->id() + dInfo.emissionData()->id() == 0 ?
	0.*GeV2 : mi2;

	double xInv = ( 1. +
	(ptpt+(1.-z)mi2+zm2-z(1.-z)*Mi2) /
	(z(1.-z)sqr(dInfo.scale())) );

	Lorentz5Momentum em = z*pEmitter +
	(sqr(pt)+mi2-zzMi2)/(zsqr(dInfo.scale()))pSpectator + kt;
	em.setMass(sqrt(mi2));
	em.rescaleEnergy();

	Lorentz5Momentum emm = (1.-z)*pEmitter +
	(ptpt+m2-sqr(1.-z)Mi2)/((1.-z)sqr(dInfo.scale()))pSpectator - kt;
	emm.setMass(sqrt(m2));
	emm.rescaleEnergy();

	Lorentz5Momentum spe = xInv*pSpectator;
	spe.setMass(ZERO);
	spe.rescaleEnergy();

	emitterMomentum(em);
	emissionMomentum(emm);
	spectatorMomentum(spe);

	}

	// If needed, insert default implementations of function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void FIMassiveKinematics::persistentOutput(PersistentOStream & ) const {
	}

	void FIMassiveKinematics::persistentInput(PersistentIStream & , int) {
	}

	ClassDescription<FIMassiveKinematics> FIMassiveKinematics::initFIMassiveKinematics;
	// Definition of the static class description member.

	void FIMassiveKinematics::Init() {

	static ClassDocumentation<FIMassiveKinematics> documentation
	("FIMassiveKinematics implements massless splittings "
	"off a final-initial dipole.");

	}

	diff --git a/DipoleShower/Kinematics/FIMassiveKinematics.h b/DipoleShower/Kinematics/FIMassiveKinematics.h
	--- a/DipoleShower/Kinematics/FIMassiveKinematics.h
	+++ b/DipoleShower/Kinematics/FIMassiveKinematics.h
	@@ -1,244 +1,246 @@
	// -- C++ --
	//
	// FILightKinematics.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.
	//
	#ifndef HERWIG_FILightKinematics_H
	#define HERWIG_FILightKinematics_H
	//
	// This is the declaration of the FILightKinematics class.
	//

	#include "DipoleSplittingKinematics.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer, Martin Stoll
	*
	* \brief FIMassiveKinematics implements massless splittings
	* off a final-initial dipole.
	*
	*/
	class FIMassiveKinematics: public DipoleSplittingKinematics {

	public:

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

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

	public:

	/**
	* Return the boundaries in between the evolution
	* variable random number is to be sampled; the lower
	* cuoff is assumed to correspond to the infrared cutoff.
	*/
	virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;

	/**
	* Return the boundaries in between the momentum
	* fraction random number is to be sampled.
	*/
	virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;

	/**
	* Return the boundaries on the momentum fraction
	*/
	virtual pair<double,double> zBoundaries(Energy,
	const DipoleSplittingInfo&,
	const DipoleSplittingKernel&) const {
	return make_pair(0.0,1.0);
	}

	/**
	* Return the dipole scale associated to the
	* given pair of emitter and spectator. This
	* should be the invariant mass or absolute value
	* final/final or initial/initial and the absolute
	* value of the momentum transfer for intial/final or
	* final/initial dipoles.
	*/
	virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator) const;

	/**
	* Return the maximum pt for the given dipole scale.
	*/
	virtual Energy ptMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex& dIndex,
	const DipoleSplittingKernel&) const;

	/**
	* Return the maximum virtuality for the given dipole scale.
	*/
	virtual Energy QMax(Energy dScale,
	double emX, double specX,
	const DipoleIndex& dIndex,
	const DipoleSplittingKernel&) const;

	/**
	* Return the pt given a virtuality.
	*/
	virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the virtuality given a pt.
	*/
	virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;

	/**
	* Return the random number associated to
	* the given pt.
	*/
	virtual double ptToRandom(Energy pt, Energy dScale,
	double emX, double specX,
	const DipoleIndex& dIndex,
	const DipoleSplittingKernel&) const;

	/**
	* Generate splitting variables given three random numbers
	* and the momentum fractions of the emitter and spectator.
	* Return true on success.
	*/
	virtual bool generateSplitting(double kappa, double xi, double phi,
	DipoleSplittingInfo& dIndex,
	const DipoleSplittingKernel&);

	/**
	* Generate the full kinematics given emitter and
	* spectator momentum and a previously completeted
	* DipoleSplittingInfo object.
	*/
	virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
	const Lorentz5Momentum& pSpectator,
	const DipoleSplittingInfo& dInfo);

	public:

	/**
	* Triangular / Kallen function
	*/
	template <class T>
	inline double rootOfKallen (T a, T b, T c) const {
	- return sqrt( aa + bb + cc - 2.( ab+ac+b*c ) ); }
	+ double sres=aa + bb + cc - 2.( ab+ac+b*c );
	+ return sres>0.?sqrt( sres ):0.; }
	+

	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 static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<FIMassiveKinematics> initFIMassiveKinematics;

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the FIMassiveKinematics class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::FIMassiveKinematics>
	: public ClassTraitsBase<Herwig::FIMassiveKinematics> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::FIMassiveKinematics"; }
	/**
	* The name of a file containing the dynamic library where the class
	* FIMassiveKinematics is implemented. It may also include several, space-separated,
	* libraries if the class FIMassiveKinematics depends on other classes (base classes
	* excepted). In this case the listed libraries will be dynamically
	* linked in the order they are specified.
	*/
	static string library() { return "HwDipoleShower.so"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_FIMassiveKinematics_H */
	diff --git a/MatrixElement/Matchbox/Matching/DipoleMatching.cc b/MatrixElement/Matchbox/Matching/DipoleMatching.cc
	--- a/MatrixElement/Matchbox/Matching/DipoleMatching.cc
	+++ b/MatrixElement/Matchbox/Matching/DipoleMatching.cc
	@@ -1,129 +1,139 @@
	// -- C++ --
	//
	// DipoleMatching.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 DipoleMatching class.
	//

	#include "DipoleMatching.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"


	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"

	using namespace Herwig;

	DipoleMatching::DipoleMatching() {}

	DipoleMatching::~DipoleMatching() {}

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

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

	CrossSection DipoleMatching::dSigHatDR() const {

	double xme2 = 0.;

	pair<int,int> ij(dipole()->bornEmitter(),
	dipole()->bornSpectator());
	double ccme2 =
	dipole()->underlyingBornME()->largeNColourCorrelatedME2(ij,theLargeNBasis);

	+ if(ccme2==0.)return 0.*nanobarn;
	+
	+ double lnme2=dipole()->underlyingBornME()->largeNME2(theLargeNBasis);
	+ if(lnme2==0){
	+ generator()->log() <<"\nQTildeMatching: ";
	+ generator()->log() <<"\n LargeNME2 is ZERO, while largeNColourCorrelatedME2 is not ZERO." ;
	+ generator()->log() <<"\n This is too seriuos.\n" ;
	+ generator()->log() << Exception::runerror;
	+ }
	+
	+
	ccme2 *=
	- dipole()->underlyingBornME()->me2() /
	- dipole()->underlyingBornME()->largeNME2(theLargeNBasis);
	+ dipole()->underlyingBornME()->me2() /lnme2;

	xme2 = dipole()->me2Avg(ccme2);

	xme2 /= dipole()->underlyingBornME()->lastXComb().lastAlphaS();
	double bornPDF = bornPDFWeight(dipole()->underlyingBornME()->lastScale());
	if ( bornPDF == 0.0 )
	return ZERO;
	xme2 *= bornPDF;

	if ( profileScales() )
	xme2 *= profileScales()->hardScaleProfile(dipole()->showerHardScale(),dipole()->lastPt());

	CrossSection res =
	sqr(hbarc) *
	realXComb()->jacobian() *
	subtractionScaleWeight() *
	xme2 /
	(2. * realXComb()->lastSHat());

	return res;

	}

	double DipoleMatching::me2() const {
	throw Exception() << "DipoleMatching::me2(): Not intented to use. Disable the ShowerApproximationGenerator."
	<< Exception::runerror;
	return 0.;
	}

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


	void DipoleMatching::persistentOutput(PersistentOStream & os) const {
	os << theShowerHandler;
	}

	void DipoleMatching::persistentInput(PersistentIStream & is, int) {
	is >> theShowerHandler;
	}

	void DipoleMatching::doinit() {
	ShowerApproximation::doinit();
	if ( theShowerHandler ) {
	if ( theShowerHandler->scaleFactorOption() < 2 ) {
	hardScaleFactor(theShowerHandler->hardScaleFactor());
	factorizationScaleFactor(theShowerHandler->factorizationScaleFactor());
	renormalizationScaleFactor(theShowerHandler->renormalizationScaleFactor());
	}
	profileScales(theShowerHandler->profileScales());
	restrictPhasespace(theShowerHandler->restrictPhasespace());
	hardScaleIsMuF(theShowerHandler->hardScaleIsMuF());
	}
	}

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

	void DipoleMatching::Init() {

	static ClassDocumentation<DipoleMatching> documentation
	("DipoleMatching implements NLO matching with the dipole shower.");

	static Reference<DipoleMatching,ShowerHandler> interfaceShowerHandler
	("ShowerHandler",
	"",
	&DipoleMatching::theShowerHandler, false, false, true, true, false);

	}

	diff --git a/MatrixElement/Matchbox/Matching/QTildeMatching.cc b/MatrixElement/Matchbox/Matching/QTildeMatching.cc
	--- a/MatrixElement/Matchbox/Matching/QTildeMatching.cc
	+++ b/MatrixElement/Matchbox/Matching/QTildeMatching.cc
	@@ -1,501 +1,511 @@
	// -- C++ --
	//
	// QTildeMatching.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 QTildeMatching class.
	//

	#include "QTildeMatching.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"


	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
	#include "Herwig++/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"

	using namespace Herwig;

	QTildeMatching::QTildeMatching()
	: theCorrectForXZMismatch(true) {}

	QTildeMatching::~QTildeMatching() {}

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

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

	void QTildeMatching::checkCutoff() {
	if ( showerTildeKinematics() ) {
	showerTildeKinematics()->
	prepare(realCXComb(),bornCXComb());
	showerTildeKinematics()->dipole(dipole());
	showerTildeKinematics()->getShowerVariables();
	}
	}

	void QTildeMatching::getShowerVariables() {

	// already filled from checkCutoff in this case
	if ( showerTildeKinematics() )
	return;

	// get the shower variables
	calculateShowerVariables();

	// check for the cutoff
	dipole()->isAboveCutoff(isAboveCutoff());

	// get the hard scale
	dipole()->showerHardScale(hardScale());

	// check for phase space
	dipole()->isInShowerPhasespace(isInShowerPhasespace());

	}

	bool QTildeMatching::isInShowerPhasespace() const {

	assert((theQTildeSudakov->cutOffOption() == 0 \|\| theQTildeSudakov->cutOffOption() == 2) &&
	"implementation only provided for default and pt cutoff");

	Energy qtildeHard = ZERO;

	Energy qtilde = dipole()->showerScale();
	assert(!dipole()->showerParameters().empty());
	double z = dipole()->showerParameters()[0];

	// FF
	if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) {
	qtildeHard =
	theQTildeFinder->
	calculateFinalFinalScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
	bornCXComb()->meMomenta()[dipole()->bornSpectator()],
	bornCXComb()->mePartonData()[dipole()->bornEmitter()]->iColour() == PDT::Colour3).first;
	}

	// FI
	if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) {
	qtildeHard =
	theQTildeFinder->
	calculateInitialFinalScales(bornCXComb()->meMomenta()[dipole()->bornSpectator()],
	bornCXComb()->meMomenta()[dipole()->bornEmitter()],false).second;
	}

	// IF
	if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) {
	qtildeHard =
	theQTildeFinder->
	calculateInitialFinalScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
	bornCXComb()->meMomenta()[dipole()->bornSpectator()],false).first;
	if ( z < (dipole()->bornEmitter() == 0 ? bornCXComb()->lastX1() : bornCXComb()->lastX2()) )
	return false;
	}

	// II
	if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2 ) {
	qtildeHard =
	theQTildeFinder->
	calculateInitialInitialScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
	bornCXComb()->meMomenta()[dipole()->bornSpectator()]).first;
	if ( z < (dipole()->bornEmitter() == 0 ? bornCXComb()->lastX1() : bornCXComb()->lastX2()) )
	return false;
	}


	Energy Qg = theQTildeSudakov->kinScale();
	Energy2 pt2 = ZERO;
	if ( dipole()->bornEmitter() > 1 ) {
	Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
	if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
	pt2 = sqr(z(1.-z)qtilde) - sqr(mu);
	else
	pt2 = sqr(z(1.-z)qtilde) - sqr((1.-z)mu) - zsqr(Qg);
	}
	if ( dipole()->bornEmitter() < 2 ) {
	pt2 = sqr((1.-z)qtilde) - zsqr(Qg);
	}

	if ( pt2 < theQTildeSudakov->pT2min() )
	return false;

	bool hardVeto = restrictPhasespace() && sqrt(pt2) >= dipole()->showerHardScale();
	return qtilde <= qtildeHard && !hardVeto;

	}

	bool QTildeMatching::isAboveCutoff() const {

	assert((theQTildeSudakov->cutOffOption() == 0 \|\| theQTildeSudakov->cutOffOption() == 2) &&
	"implementation only provided for default and pt cutoff");

	Energy qtilde = dipole()->showerScale();
	assert(!dipole()->showerParameters().empty());
	double z = dipole()->showerParameters()[0];
	Energy Qg = theQTildeSudakov->kinScale();
	if ( dipole()->bornEmitter() > 1 ) {
	Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
	if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
	return sqr(z(1.-z)qtilde) - sqr(mu) >= theQTildeSudakov->pT2min();
	else
	return sqr(z(1.-z)qtilde) - sqr((1.-z)mu) - zsqr(Qg) >= theQTildeSudakov->pT2min();
	}
	if ( dipole()->bornEmitter() < 2 ) {
	return
	sqr((1.-z)qtilde) - zsqr(Qg) >= theQTildeSudakov->pT2min();
	}
	return false;
	}

	CrossSection QTildeMatching::dSigHatDR() const {

	assert(!dipole()->showerParameters().empty());
	pair<Energy2,double> vars =
	make_pair(sqr(dipole()->showerScale()),
	dipole()->showerParameters()[0]);

	pair<int,int> ij(dipole()->bornEmitter(),
	dipole()->bornSpectator());
	double ccme2 =
	dipole()->underlyingBornME()->largeNColourCorrelatedME2(ij,theLargeNBasis);
	-
	+
	+ if(ccme2==0.)return 0.*nanobarn;
	+
	+ double lnme2=dipole()->underlyingBornME()->largeNME2(theLargeNBasis);
	+ if(lnme2==0){
	+ generator()->log() <<"\nQTildeMatching: ";
	+ generator()->log() <<"\n largeNME2 is ZERO, while largeNColourCorrelatedME2 is not ZERO." ;
	+ generator()->log() <<"\n This is too seriuos.\n" ;
	+ generator()->log() << Exception::runerror;
	+ }
	+
	ccme2 *=
	- dipole()->underlyingBornME()->me2() /
	- dipole()->underlyingBornME()->largeNME2(theLargeNBasis);
	+ dipole()->underlyingBornME()->me2() /lnme2;
	+

	Energy2 prop = ZERO;
	if ( dipole()->bornEmitter() > 1 ) {
	prop =
	(realCXComb()->meMomenta()[dipole()->realEmitter()] +
	realCXComb()->meMomenta()[dipole()->realEmission()]).m2()
	- bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2();
	} else {
	prop =
	2.vars.second(realCXComb()->meMomenta()[dipole()->realEmitter()]*
	realCXComb()->meMomenta()[dipole()->realEmission()]);
	}

	// note alphas included downstream from subtractionScaleWeight()
	double xme2 = -8.Constants::piccme2splitFn(vars)realXComb()->lastSHat()/prop;
	xme2 *=
	pow(realCXComb()->lastSHat() / bornCXComb()->lastSHat(),
	bornCXComb()->mePartonData().size()-4.);

	double bornPDF = bornPDFWeight(dipole()->underlyingBornME()->lastScale());
	if ( bornPDF == 0.0 )
	return ZERO;

	xme2 *= bornPDF;

	xme2 *= dipole()->realEmissionME()->finalStateSymmetry() /
	dipole()->underlyingBornME()->finalStateSymmetry();

	// take care of mismatch between z and x as we are approaching the
	// hard phase space boundary
	// TODO get rid of this useless scale option business and simplify PDF handling in here
	if ( dipole()->bornEmitter() < 2 && theCorrectForXZMismatch ) {
	Energy2 emissionScale = ZERO;
	if ( emissionScaleInSubtraction() == showerScale ) {
	emissionScale = showerFactorizationScale();
	} else if ( emissionScaleInSubtraction() == realScale ) {
	emissionScale = dipole()->realEmissionME()->lastScale();
	} else if ( emissionScaleInSubtraction() == bornScale ) {
	emissionScale = dipole()->underlyingBornME()->lastScale();
	}
	double xzMismatch =
	dipole()->subtractionParameters()[0] / dipole()->showerParameters()[0];
	double realCorrectedPDF =
	dipole()->bornEmitter() == 0 ?
	dipole()->realEmissionME()->pdf1(emissionScale,theExtrapolationX,
	xzMismatch) :
	dipole()->realEmissionME()->pdf2(emissionScale,theExtrapolationX,
	xzMismatch);
	double realPDF =
	dipole()->bornEmitter() == 0 ?
	dipole()->realEmissionME()->pdf1(emissionScale,theExtrapolationX,1.0) :
	dipole()->realEmissionME()->pdf2(emissionScale,theExtrapolationX,1.0);
	if ( realPDF == 0.0 \|\| realCorrectedPDF == 0.0 )
	return ZERO;
	xme2 *= realCorrectedPDF / realPDF;
	}

	Energy qtilde = sqrt(vars.first);
	double z = vars.second;
	Energy2 pt2 = ZERO;
	Energy Qg = theQTildeSudakov->kinScale();
	if ( dipole()->bornEmitter() > 1 ) {
	Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
	if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
	pt2 = sqr(z(1.-z)qtilde) - sqr(mu);
	else
	pt2 = sqr(z(1.-z)qtilde) - sqr((1.-z)mu) - zsqr(Qg);
	}
	if ( dipole()->bornEmitter() < 2 ) {
	pt2 = sqr((1.-z)qtilde) - zsqr(Qg);
	}
	assert(pt2 >= ZERO);

	if ( profileScales() )
	xme2 *= profileScales()->hardScaleProfile(dipole()->showerHardScale(),sqrt(pt2));

	CrossSection res =
	sqr(hbarc) *
	realXComb()->jacobian() *
	subtractionScaleWeight() *
	xme2 /
	(2. * realXComb()->lastSHat());

	return res;

	}

	double QTildeMatching::me2() const {
	throw Exception() << "QTildeMatching::me2(): Not intented to use. Disable the ShowerApproximationGenerator."
	<< Exception::runerror;
	return 0.;
	}

	void QTildeMatching::calculateShowerVariables() const {

	Lorentz5Momentum n;
	Energy2 Q2 = ZERO;

	const Lorentz5Momentum& pb = bornCXComb()->meMomenta()[dipole()->bornEmitter()];
	const Lorentz5Momentum& pc = bornCXComb()->meMomenta()[dipole()->bornSpectator()];

	if ( dipole()->bornEmitter() > 1 ) {
	Q2 = (pb+pc).m2();
	} else {
	Q2 = -(pb-pc).m2();
	}

	if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) {
	double b = sqr(bornCXComb()->meMomenta()[dipole()->bornEmitter()].m())/Q2;
	double c = sqr(bornCXComb()->meMomenta()[dipole()->bornSpectator()].m())/Q2;
	double lambda = sqrt(1.+sqr(b)+sqr(c)-2.b-2.c-2.bc);
	n = (1.-0.5(1.-b+c-lambda))pc - 0.5(1.-b+c-lambda)pb;
	}

	if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) {
	n = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
	}

	if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) {
	double c = sqr(bornCXComb()->meMomenta()[dipole()->bornSpectator()].m())/Q2;
	n = (1.+c)pc - cpb;
	}

	if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2 ) {
	n = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
	}

	// the light-cone condition is numerically not very stable, so we
	// explicitly push it on the light-cone here
	n.setMass(ZERO);
	n.rescaleEnergy();

	double z = 0.0;

	if ( dipole()->bornEmitter() > 1 ) {
	z = 1. -
	(n*realCXComb()->meMomenta()[dipole()->realEmission()])/
	(n*bornCXComb()->meMomenta()[dipole()->bornEmitter()]);
	} else {
	z = 1. -
	(n*realCXComb()->meMomenta()[dipole()->realEmission()])/
	(n*realCXComb()->meMomenta()[dipole()->realEmitter()]);
	}

	Energy2 qtilde2 = ZERO;
	Energy2 q2 = ZERO;

	if ( dipole()->bornEmitter() > 1 ) {
	q2 =
	(realCXComb()->meMomenta()[dipole()->realEmitter()] + realCXComb()->meMomenta()[dipole()->realEmission()]).m2();
	qtilde2 = (q2 - bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2())/(z*(1.-z));
	} else {
	q2 =
	-(realCXComb()->meMomenta()[dipole()->realEmitter()] - realCXComb()->meMomenta()[dipole()->realEmission()]).m2();
	qtilde2 = (q2 + bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2())/(1.-z);
	}

	assert(qtilde2 >= ZERO && z >= 0.0 && z <= 1.0);

	dipole()->showerScale(sqrt(qtilde2));
	dipole()->showerParameters().resize(1);
	dipole()->showerParameters()[0] = z;

	}

	double QTildeMatching::splitFn(const pair<Energy2,double>& vars) const {
	const Energy2& qtilde2 = vars.first;
	const double& z = vars.second;
	double Nc = SM().Nc();

	// final state branching
	if ( dipole()->bornEmitter() > 1 ) {
	// final state quark quark branching
	if ( abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 7 ) {
	Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
	return
	((sqr(Nc)-1.)/(2.Nc))(1+sqr(z)-2.sqr(m)/(zqtilde2))/(1.-z);
	}
	// final state gluon branching
	if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g ) {
	if ( realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
	// ATTENTION the factor 2 here is intentional as it cancels to the 1/2
	// stemming from the large-N colour correlator
	return 2.Nc(z/(1.-z)+(1.-z)/z+z*(1.-z));
	}
	if ( abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
	Energy m = realCXComb()->mePartonData()[dipole()->realEmission()]->hardProcessMass();
	return (1./2.)(1.-2.z(1.-z)+2.sqr(m)/(z(1.-z)qtilde2));
	}
	}
	// final state squark branching
	if ((abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) > 1000000 &&
	abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 1000007) \|\|
	(abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) > 2000000 &&
	abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 2000007)){
	Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
	return ((sqr(Nc)-1.)/Nc)(z-sqr(m)/(zqtilde2))/(1.-z);
	}
	// final state gluino branching
	if (bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == 1000021){
	Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
	return Nc(1.+sqr(z)-2.sqr(m)/(z*qtilde2))/(1.-z);
	}
	}

	// initial state branching
	if ( dipole()->bornEmitter() < 2 ) {
	// g/g
	if ( realCXComb()->mePartonData()[dipole()->realEmitter()]->id() == ParticleID::g &&
	realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
	// see above for factor of 2
	return 2.Nc(z/(1.-z)+(1.-z)/z+z*(1.-z));
	}
	// q/q
	if ( abs(realCXComb()->mePartonData()[dipole()->realEmitter()]->id()) < 7 &&
	realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
	return
	((sqr(Nc)-1.)/(2.Nc))(1+sqr(z))/(1.-z);
	}
	// g/q
	if ( realCXComb()->mePartonData()[dipole()->realEmitter()]->id() == ParticleID::g &&
	abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
	return (1./2.)(1.-2.z*(1.-z));
	}
	// q/g
	if ( abs(realCXComb()->mePartonData()[dipole()->realEmitter()]->id()) < 7 &&
	abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
	return
	((sqr(Nc)-1.)/(2.Nc))(1+sqr(1.-z))/z;
	}
	}

	return 0.0;

	}

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

	void QTildeMatching::doinit() {
	assert(theShowerHandler && theQTildeFinder && theQTildeSudakov);
	theShowerHandler->init();
	theQTildeFinder->init();
	theQTildeSudakov->init();
	if ( theShowerHandler->scaleFactorOption() < 2 ) {
	hardScaleFactor(theShowerHandler->hardScaleFactor());
	factorizationScaleFactor(theShowerHandler->factorizationScaleFactor());
	renormalizationScaleFactor(theShowerHandler->renormalizationScaleFactor());
	}
	profileScales(theShowerHandler->profileScales());
	restrictPhasespace(theShowerHandler->restrictPhasespace());
	hardScaleIsMuF(theShowerHandler->hardScaleIsMuF());
	ShowerApproximation::doinit();
	}

	void QTildeMatching::doinitrun() {
	assert(theShowerHandler && theQTildeFinder && theQTildeSudakov);
	theShowerHandler->initrun();
	theQTildeFinder->initrun();
	theQTildeSudakov->initrun();
	ShowerApproximation::doinitrun();
	}


	void QTildeMatching::persistentOutput(PersistentOStream & os) const {
	os << theQTildeFinder << theQTildeSudakov
	<< theShowerHandler << theCorrectForXZMismatch;
	}

	void QTildeMatching::persistentInput(PersistentIStream & is, int) {
	is >> theQTildeFinder >> theQTildeSudakov
	>> theShowerHandler >> theCorrectForXZMismatch;
	}


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

	void QTildeMatching::Init() {

	static ClassDocumentation<QTildeMatching> documentation
	("QTildeMatching implements NLO matching with the default shower.");

	static Reference<QTildeMatching,QTildeFinder> interfaceQTildeFinder
	("QTildeFinder",
	"Set the partner finder to calculate hard scales.",
	&QTildeMatching::theQTildeFinder, false, false, true, false, false);

	static Reference<QTildeMatching,QTildeSudakov> interfaceQTildeSudakov
	("QTildeSudakov",
	"Set the partner finder to calculate hard scales.",
	&QTildeMatching::theQTildeSudakov, false, false, true, false, false);

	static Reference<QTildeMatching,ShowerHandler> interfaceShowerHandler
	("ShowerHandler",
	"",
	&QTildeMatching::theShowerHandler, false, false, true, true, false);

	static Switch<QTildeMatching,bool> interfaceCorrectForXZMismatch
	("CorrectForXZMismatch",
	"Correct for x/z mismatch near hard phase space boundary.",
	&QTildeMatching::theCorrectForXZMismatch, true, false, false);
	static SwitchOption interfaceCorrectForXZMismatchYes
	(interfaceCorrectForXZMismatch,
	"Yes",
	"Include the correction factor.",
	true);
	static SwitchOption interfaceCorrectForXZMismatchNo
	(interfaceCorrectForXZMismatch,
	"No",
	"Do not include the correction factor.",
	false);

	}

	diff --git a/Shower/Couplings/ShowerAlphaQCD.cc b/Shower/Couplings/ShowerAlphaQCD.cc
	--- a/Shower/Couplings/ShowerAlphaQCD.cc
	+++ b/Shower/Couplings/ShowerAlphaQCD.cc
	@@ -1,377 +1,372 @@
	// -- C++ --
	//
	// ShowerAlphaQCD.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ShowerAlphaQCD class.
	//
	#include "ShowerAlphaQCD.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Command.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Config/Constants.h"

	using namespace Herwig;

	DescribeClass<ShowerAlphaQCD,ShowerAlpha>
	describeShowerAlphaQCD("Herwig::ShowerAlphaQCD","HwShower.so");

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

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

	void ShowerAlphaQCD::persistentOutput(PersistentOStream & os) const {
	os << _asType << _asMaxNP << ounit(_qmin,GeV) << _nloop << _lambdaopt << _thresopt
	<< ounit(_lambdain,GeV) << _alphain << _inopt
	<< _tolerance << _maxtry << _alphamin
	<< ounit(_thresholds,GeV) << ounit(_lambda,GeV);
	}

	void ShowerAlphaQCD::persistentInput(PersistentIStream & is, int) {
	is >> _asType >> _asMaxNP >> iunit(_qmin,GeV) >> _nloop >> _lambdaopt >> _thresopt
	>> iunit(_lambdain,GeV) >> _alphain >> _inopt
	>> _tolerance >> _maxtry >> _alphamin
	>> iunit(_thresholds,GeV) >> iunit(_lambda,GeV);
	}

	void ShowerAlphaQCD::Init() {

	static ClassDocumentation<ShowerAlphaQCD> documentation
	("This (concrete) class describes the QCD alpha running.");

	static Switch<ShowerAlphaQCD, int> intAsType
	("NPAlphaS",
	"Behaviour of AlphaS in the NP region",
	&ShowerAlphaQCD::_asType, 1, false, false);
	static SwitchOption intAsTypeZero
	(intAsType, "Zero","zero below Q_min", 1);
	static SwitchOption intAsTypeConst
	(intAsType, "Const","const as(qmin) below Q_min", 2);
	static SwitchOption intAsTypeLin
	(intAsType, "Linear","growing linearly below Q_min", 3);
	static SwitchOption intAsTypeQuad
	(intAsType, "Quadratic","growing quadratically below Q_min", 4);
	static SwitchOption intAsTypeExx1
	(intAsType, "Exx1", "quadratic from AlphaMaxNP down to as(Q_min)", 5);
	static SwitchOption intAsTypeExx2
	(intAsType, "Exx2", "const = AlphaMaxNP below Q_min", 6);

	// default such that as(qmin) = 1 in the current parametrization.
	// min = Lambda3
	static Parameter<ShowerAlphaQCD,Energy> intQmin
	("Qmin", "Q < Qmin is treated with NP parametrization as of (unit [GeV])",
	&ShowerAlphaQCD::_qmin, GeV, 0.630882GeV, 0.330445GeV,
	100.0*GeV,false,false,false);

	static Parameter<ShowerAlphaQCD,double> interfaceAlphaMaxNP
	("AlphaMaxNP",
	"Max value of alpha in NP region, only relevant if NPAlphaS = 5,6",
	&ShowerAlphaQCD::_asMaxNP, 1.0, 0., 100.0,
	false, false, Interface::limited);

	static Parameter<ShowerAlphaQCD,unsigned int> interfaceNumberOfLoops
	("NumberOfLoops",
	"The number of loops to use in the alpha_S calculation",
	&ShowerAlphaQCD::_nloop, 3, 1, 3,
	false, false, Interface::limited);

	static Switch<ShowerAlphaQCD,bool> interfaceLambdaOption
	("LambdaOption",
	"Option for the calculation of the Lambda used in the simulation from the input"
	" Lambda_MSbar",
	&ShowerAlphaQCD::_lambdaopt, false, false, false);
	static SwitchOption interfaceLambdaOptionfalse
	(interfaceLambdaOption,
	"Same",
	"Use the same value",
	false);
	static SwitchOption interfaceLambdaOptionConvert
	(interfaceLambdaOption,
	"Convert",
	"Use the conversion to the Herwig scheme from NPB349, 635",
	true);

	static Parameter<ShowerAlphaQCD,Energy> interfaceLambdaQCD
	("LambdaQCD",
	"Input value of Lambda_MSBar",
	&ShowerAlphaQCD::_lambdain, MeV, 0.208364GeV, 100.0MeV, 500.0*MeV,
	false, false, Interface::limited);

	static Parameter<ShowerAlphaQCD,double> interfaceAlphaMZ
	("AlphaMZ",
	"The input value of the strong coupling at the Z mass ",
	&ShowerAlphaQCD::_alphain, 0.118, 0.1, 0.2,
	false, false, Interface::limited);

	static Switch<ShowerAlphaQCD,bool> interfaceInputOption
	("InputOption",
	"Option for inputing the initial value of the coupling",
	&ShowerAlphaQCD::_inopt, true, false, false);
	static SwitchOption interfaceInputOptionAlphaMZ
	(interfaceInputOption,
	"AlphaMZ",
	"Use the value of alpha at MZ to calculate the coupling",
	true);
	static SwitchOption interfaceInputOptionLambdaQCD
	(interfaceInputOption,
	"LambdaQCD",
	"Use the input value of Lambda to calculate the coupling",
	false);

	static Parameter<ShowerAlphaQCD,double> interfaceTolerance
	("Tolerance",
	"The tolerance for discontinuities in alphaS at thresholds.",
	&ShowerAlphaQCD::_tolerance, 1e-10, 1e-20, 1e-4,
	false, false, Interface::limited);

	static Parameter<ShowerAlphaQCD,unsigned int> interfaceMaximumIterations
	("MaximumIterations",
	"The maximum number of iterations for the Newton-Raphson method to converge.",
	&ShowerAlphaQCD::_maxtry, 100, 10, 1000,
	false, false, Interface::limited);

	static Switch<ShowerAlphaQCD,bool> interfaceThresholdOption
	("ThresholdOption",
	"Whether to use the consistuent or normal masses for the thresholds",
	&ShowerAlphaQCD::_thresopt, true, false, false);
	static SwitchOption interfaceThresholdOptionCurrent
	(interfaceThresholdOption,
	"Current",
	"Use the current masses",
	true);
	static SwitchOption interfaceThresholdOptionConstituent
	(interfaceThresholdOption,
	"Constituent",
	"Use the constitent masses.",
	false);

	static Command<ShowerAlphaQCD> interfaceValue
	("Value",
	"",
	&ShowerAlphaQCD::value, false);

	}

	void ShowerAlphaQCD::doinit() {
	ShowerAlpha::doinit();
	// calculate the value of 5-flavour lambda
	// evaluate the initial
	// value of Lambda from alphas if needed using Newton-Raphson
	- if(_inopt)
	- {_lambda[2]=computeLambda(getParticleData(ParticleID::Z0)->mass(),_alphain,5);}
	+ if(_inopt) {
	+ _lambda[2]=computeLambda(getParticleData(ParticleID::Z0)->mass(),_alphain,5);
	+ }
	// otherwise it was an input parameter
	else{_lambda[2]=_lambdain;}
	// convert lambda to the Monte Carlo scheme if needed
	using Constants::pi;
	- if(_lambdaopt){_lambda[2] =exp(0.5(67.-3.*sqr(pi)-50./3.)/23.)/sqrt(2.);}
	+ if(_lambdaopt) _lambda[2] =exp(0.5(67.-3.*sqr(pi)-50./3.)/23.)/sqrt(2.);
	// compute the threshold matching
	// top threshold
	for(int ix=1;ix<4;++ix) {
	if(_thresopt)
	_thresholds[ix]=getParticleData(ix+3)->mass();
	else
	_thresholds[ix]=getParticleData(ix+3)->constituentMass();
	}
	// compute 6 flavour lambda by matching at top mass using Newton Raphson
	_lambda[3]=computeLambda(_thresholds[3],alphaS(_thresholds[3],_lambda[2],5),6);
	// bottom threshold
	// compute 4 flavour lambda by matching at bottom mass using Newton Raphson
	_lambda[1]=computeLambda(_thresholds[2],alphaS(_thresholds[2],_lambda[2],5),4);
	// charm threshold
	// compute 3 flavour lambda by matching at charm mass using Newton Raphson
	_lambda[0]=computeLambda(_thresholds[1],alphaS(_thresholds[1],_lambda[1],4),3);
	// final threshold is qmin
	_thresholds[0]=_qmin;
	// compute the maximum value of as
	if ( _asType < 5 ) _alphamin = value(sqr(_qmin)+1.0e-8*sqr(MeV)); // approx as = 1
	else _alphamin = max(_asMaxNP, value(sqr(_qmin)+1.0e-8*sqr(MeV)));
	// check consistency lambda_3 < qmin
	if(_lambda[0]>_qmin)
	Throw<InitException>() << "The value of Qmin is less than Lambda_3 in"
	<< " ShowerAlphaQCD::doinit " << Exception::abortnow;
	-
	}

	double ShowerAlphaQCD::value(const Energy2 scale) const {
	pair<short,Energy> nflam;
	Energy q = sqrt(scale);
	double val(0.);
	// special handling if the scale is less than Qmin
	if (q < _qmin) {
	nflam = getLamNfTwoLoop(_qmin);
	double val0 = alphaS(_qmin, nflam.second, nflam.first);
	switch (_asType) {
	case 1:
	// flat, zero; the default type with no NP effects.
	val = 0.;
	break;
	case 2:
	// flat, non-zero alpha_s = alpha_s(q2min).
	val = val0;
	break;
	case 3:
	// linear in q
	val = val0*q/_qmin;
	break;
	case 4:
	// quadratic in q
	val = val0*sqr(q/_qmin);
	break;
	case 5:
	// quadratic in q, starting off at asMaxNP, ending on as(qmin)
	val = (val0 - _asMaxNP)*sqr(q/_qmin) + _asMaxNP;
	break;
	case 6:
	// just asMaxNP and constant
	val = _asMaxNP;
	break;
	}
	} else {
	// the 'ordinary' case
	nflam = getLamNfTwoLoop(q);
	val = alphaS(q, nflam.second, nflam.first);
	}
	return scaleFactor() * val;
	}

	double ShowerAlphaQCD::overestimateValue() const {
	return scaleFactor() * _alphamin;
	}

	double ShowerAlphaQCD::ratio(const Energy2 scale) const {
	pair<short,Energy> nflam;
	Energy q = sqrt(scale);
	double val(0.);
	// special handling if the scale is less than Qmin
	if (q < _qmin) {
	nflam = getLamNfTwoLoop(_qmin);
	double val0 = alphaS(_qmin, nflam.second, nflam.first);
	switch (_asType) {
	case 1:
	// flat, zero; the default type with no NP effects.
	val = 0.;
	break;
	case 2:
	// flat, non-zero alpha_s = alpha_s(q2min).
	val = val0;
	break;
	case 3:
	// linear in q
	val = val0*q/_qmin;
	break;
	case 4:
	// quadratic in q
	val = val0*sqr(q/_qmin);
	break;
	case 5:
	// quadratic in q, starting off at asMaxNP, ending on as(qmin)
	val = (val0 - _asMaxNP)*sqr(q/_qmin) + _asMaxNP;
	break;
	case 6:
	// just asMaxNP and constant
	val = _asMaxNP;
	break;
	}
	} else {
	// the 'ordinary' case
	nflam = getLamNfTwoLoop(q);
	val = alphaS(q, nflam.second, nflam.first);
	}
	// denominator
	return val/_alphamin;
	}

	string ShowerAlphaQCD::value (string scale) {
	istringstream readscale(scale);
	double inScale; readscale >> inScale;
	Energy theScale = inScale * GeV;
	initialize();
	ostringstream showvalue ("");
	showvalue << "alpha_s (" << theScale/GeV << " GeV) = "
	<< value (sqr(theScale));
	return showvalue.str();
	}

	pair<short, Energy> ShowerAlphaQCD::getLamNfTwoLoop(Energy q) const {
	short nf = 6;
	// get lambda and nf according to the thresholds
	if (q < _thresholds[1]) nf = 3;
	else if (q < _thresholds[2]) nf = 4;
	else if (q < _thresholds[3]) nf = 5;
	return pair<short,Energy>(nf, _lambda[nf-3]);
	}

	Energy ShowerAlphaQCD::computeLambda(Energy match,
	double alpha,
	unsigned int nflav) const {
	Energy lamtest=200.0*MeV;
	double xtest;
	unsigned int ntry=0;
	- do
	- {
	- ++ntry;
	- xtest=log(sqr(match/lamtest));
	- xtest+= (alpha-alphaS(match,lamtest,nflav))/derivativealphaS(match,lamtest,nflav);
	- lamtest=match/exp(0.5*xtest);
	- }
	+ do {
	+ ++ntry;
	+ xtest = log(sqr(match/lamtest));
	+ xtest += (alpha-alphaS(match,lamtest,nflav))/derivativealphaS(match,lamtest,nflav);
	+ Energy newLambda = match/exp(0.5*xtest);
	+ lamtest = newLambda<match ? newLambda : 0.5*(lamtest+match);
	+ }
	while(abs(alpha-alphaS(match,lamtest,nflav)) > _tolerance && ntry < _maxtry);
	return lamtest;
	}

	-double ShowerAlphaQCD::derivativealphaS(Energy q, Energy lam, int nf) const
	-{
	+double ShowerAlphaQCD::derivativealphaS(Energy q, Energy lam, int nf) const {
	using Constants::pi;
	double lx = log(sqr(q/lam));
	double b0 = 11. - 2./3.*nf;
	double b1 = 51. - 19./3.*nf;
	double b2 = 2857. - 5033./9.nf + 325./27.sqr(nf);
	if(_nloop==1)
	return -4.pi/(b0sqr(lx));
	else if(_nloop==2)
	- return -4.pi/(b0sqr(lx))(1.-2.b1/sqr(b0)/lx(1.-2.log(lx)));
	+ return -4.pi/(b0sqr(lx))(1.+2.b1/sqr(b0)/lx(1.-2.log(lx)));
	else
	- return -4.pi/(b0sqr(lx))*(1.
	- - 2.b1/sqr(b0)/lx(1.-2.*log(lx))
	- + 4.sqr(b1)/(sqr(sqr(b0))sqr(lx))*
	- (1.
	- - 2.*log(lx)
	- + 3.(sqr(log(lx) - 0.5)+b2b0/(8.*sqr(b1))-1.25)
	- )
	- );
	+ return -4.pi/(b0sqr(lx))*
	+ (1. + 2.b1/sqr(b0)/lx(1.-2.*log(lx))
	+ + 4.sqr(b1)/(sqr(sqr(b0))sqr(lx))(1. - 2.log(lx)
	+ + 3.(sqr(log(lx) - 0.5)+b2b0/(8.*sqr(b1))-1.25)));
	}

	double ShowerAlphaQCD::alphaS(Energy q, Energy lam, int nf) const {
	using Constants::pi;
	double lx(log(sqr(q/lam)));
	double b0 = 11. - 2./3.*nf;
	double b1 = 51. - 19./3.*nf;
	double b2 = 2857. - 5033./9.nf + 325./27.sqr(nf);
	// one loop
	if(_nloop==1)
	{return 4.pi/(b0lx);}
	// two loop
	else if(_nloop==2) {
	return 4.pi/(b0lx)(1.-2.b1/sqr(b0)*log(lx)/lx);
	}
	// three loop
	else
	{return 4.pi/(b0lx)(1.-2.b1/sqr(b0)*log(lx)/lx +
	4.sqr(b1)/(sqr(sqr(b0))sqr(lx))*
	(sqr(log(lx) - 0.5) + b2b0/(8.sqr(b1)) - 5./4.));}
	}

	diff --git a/Tests/Rivet/Templates/Hadron-Dipole.in b/Tests/Rivet/Templates/Hadron-Dipole.in
	--- a/Tests/Rivet/Templates/Hadron-Dipole.in
	+++ b/Tests/Rivet/Templates/Hadron-Dipole.in
	@@ -1,42 +1,55 @@
	##################################################
	## Collider type
	##################################################
	read Matchbox/PPCollider.in

	##################################################
	## Matrix element library selection
	##################################################
	read Matchbox/MadGraph-OpenLoops.in
	##################################################
	## Matching and shower selection
	## Please also see flavour scheme settings
	## towards the end of the input file.
	##################################################
	read Matchbox/MCatNLO-DipoleShower.in
	##################################################
	## PDF choice
	##################################################
	read Matchbox/FiveFlavourScheme.in
	## required for dipole shower in five flavour scheme
	# read Matchbox/FiveFlavourNoBMassScheme.in
	read Matchbox/MMHT2014.in
	##################################################
	# Create the Herwig analysis
	##################################################
	create ThePEG::RivetAnalysis /Herwig/Analysis/RivetAnalysis RivetAnalysis.so
	insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
	set /Herwig/Analysis/RivetAnalysis:Debug No

	##################################################
	## Save the generator
	##################################################
	do /Herwig/MatrixElements/Matchbox/Factory:ProductionMode

	## Select the process
	cd /Herwig/MatrixElements/Matchbox
	+
	+do Factory:StartParticleGroup pnob
	+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
	+
	${process}

	read ${parameterFile}

	cd /Herwig/Generators
	saverun ${runname} EventGenerator
	diff --git a/Tests/python/make_input_files.py b/Tests/python/make_input_files.py
	--- a/Tests/python/make_input_files.py
	+++ b/Tests/python/make_input_files.py
	@@ -1,494 +1,495 @@
	#! /usr/bin/env python
	import logging,sys,os
	from string import strip, Template

	import sys
	if sys.version_info[:3] < (2,4,0):
	print "rivet scripts require Python version >= 2.4.0... exiting"
	sys.exit(1)

	if __name__ == "__main__":
	import logging
	from optparse import OptionParser, OptionGroup
	parser = OptionParser(usage="%prog name [...]")

	(opts, args) = parser.parse_args()
	## Check args
	if len(args) != 1:
	logging.error("Must specify at least input file")
	sys.exit(1)

	name=args[0]

	collider=""
	# select the template to load
	# collider
	if(name.find("BFactory")==0) :
	collider="BFactory"
	elif(name.find("LEP")==0) :
	collider="LEP"
	elif(name.find("DIS")==0) :
	collider="DIS"
	elif(name.find("TVT")==0) :
	collider="TVT"
	elif(name.find("LHC")==0) :
	collider="LHC"
	simulation=""
	if(name.find("Matchbox")>0) :
	simulation="Matchbox"
	elif(name.find("Dipole")>0) :
	simulation="Dipole"
	elif(name.find("Powheg")>0) :
	simulation="Powheg"

	if(collider=="") :
	logging.error("Can\'t find collider")
	sys.exit(1)

	# find the template
	if(simulation=="") :
	if(collider.find("TVT")>=0 or collider.find("LHC")>=0) :
	templateName="Hadron.in"
	elif(collider.find("BFactory")<0) :
	templateName= "%s.in" % (collider)
	else :
	templateName= "LEP.in"
	else :
	if(collider.find("TVT")>=0 or collider.find("LHC")>=0) :
	templateName= "Hadron-%s.in" % (simulation)
	elif(collider.find("BFactory")<0) :
	templateName= "%s-%s.in" % (collider,simulation)
	else :
	templateName= "LEP-%s.in" % (simulation)
	with open(os.path.join("Rivet/Templates",templateName), 'r') as f:
	templateText = f.read()
	template = Template( templateText )

	# work out the name of the parameter file
	if(simulation=="") :
	istart = 1
	else :
	istart = 2
	nameSplit=name.split("-")
	parameterName=nameSplit[istart]
	for i in range(istart+1,len(nameSplit)) :
	parameterName += "-%s" % nameSplit[i]

	# work out the process and parameters
	process=""
	bscheme="read Matchbox/FiveFlavourScheme.in"
	# Bfactory
	if(collider=="BFactory") :
	if(simulation=="") :
	process = "set /Herwig/MatrixElements/MEee2gZ2qq:MaximumFlavour 4"
	if(parameterName=="10.58") :
	process += "\ncreate Herwig::MEee2VectorMeson /Herwig/MatrixElements/MEUpsilon HwMELepton.so\nset /Herwig/MatrixElements/MEUpsilon:VectorMeson /Herwig/Particles/Upsilon(4S)\nset /Herwig/MatrixElements/MEUpsilon:Coupling 0.0004151809\ninsert /Herwig/MatrixElements/SimpleEE:MatrixElements 0 /Herwig/MatrixElements/MEUpsilon"
	elif(simulation=="Powheg") :
	process = "set /Herwig/MatrixElements/PowhegMEee2gZ2qq:MaximumFlavour 4"
	if(parameterName=="10.58") :
	process += "\ncreate Herwig::MEee2VectorMeson /Herwig/MatrixElements/MEUpsilon HwMELepton.so\nset /Herwig/MatrixElements/MEUpsilon:VectorMeson /Herwig/Particles/Upsilon(4S)\nset /Herwig/MatrixElements/MEUpsilon:Coupling 0.0004151809\ninsert /Herwig/MatrixElements/SimpleEE:MatrixElements 0 /Herwig/MatrixElements/MEUpsilon"
	elif(simulation=="Matchbox" or simulation == "Dipole" ) :
	process = "do Factory:Process e- e+ -> u ubar\ndo Factory:Process e- e+ -> d dbar\ndo Factory:Process e- e+ -> c cbar\ndo Factory:Process e- e+ -> s sbar"
	if(parameterName=="10.58") :
	process += "\ninsert /Herwig/Generators/EventGenerator:EventHandler:SubProcessHandlers 0 /Herwig/MatrixElements/SimpleEE\ncreate Herwig::MEee2VectorMeson /Herwig/MatrixElements/MEUpsilon HwMELepton.so\nset /Herwig/MatrixElements/MEUpsilon:VectorMeson /Herwig/Particles/Upsilon(4S)\nset /Herwig/MatrixElements/MEUpsilon:Coupling 0.0004151809\ninsert /Herwig/MatrixElements/SimpleEE:MatrixElements 0 /Herwig/MatrixElements/MEUpsilon\n"
	# DIS
	elif(collider=="DIS") :
	if(simulation=="") :
	if(parameterName.find("NoME")>=0) :
	process = "set /Herwig/Shower/Evolver:MECorrMode 0"
	parameterName=parameterName.replace("NoME-","")
	else :
	process = ""
	elif(simulation=="Powheg") :
	process = ""
	elif(simulation=="Matchbox" or simulation == "Dipole" ) :
	if(parameterName.find("e-")>=0) :
	process="do Factory:Process e- p -> e- j"
	else :
	process="do Factory:Process e+ p -> e+ j"
	# LEP
	elif(collider=="LEP") :
	+ print parameterName
	if(simulation=="") :
	process=""
	- if(parameterName=="-10") :
	+ if(parameterName=="10") :
	process="set /Herwig/MatrixElements/MEee2gZ2qq:MaximumFlavour 4"
	elif(simulation=="Powheg") :
	process=""
	- if(parameterName=="-10") :
	+ if(parameterName=="10") :
	process="set /Herwig/MatrixElements/PowhegMEee2gZ2qq:MaximumFlavour 4"
	elif(simulation=="Matchbox" or simulation == "Dipole" ) :
	- if(parameterName=="-10") :
	+ if(parameterName=="10") :
	process="do Factory:Process e- e+ -> u ubar\ndo Factory:Process e- e+ -> d dbar\ndo Factory:Process e- e+ -> c cbar\ndo Factory:Process e- e+ -> s sbar"
	else :
	process="do Factory:Process e- e+ -> j j"
	# TVT
	elif(collider=="TVT") :
	process="set /Herwig/Generators/EventGenerator:EventHandler:BeamB /Herwig/Particles/pbar-\n"
	if(parameterName.find("Run-II")>=0) :
	process+="set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 1960.0\n"
	elif(parameterName.find("Run-I")>=0) :
	process+="set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 1800.0\n"

	if(simulation=="") :
	if(parameterName.find("Run-I-WZ")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Electron\ninsert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
	elif(parameterName.find("Run-I-W")>=0 or parameterName.find("Run-II-W")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Electron\n"
	elif(parameterName.find("Run-I-Z")>=0 or parameterName.find("Run-II-Z-e")>=0) :
	process +="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
	elif(parameterName.find("Run-II-Z-mu")>=0) :
	process +="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Muon\n"
	elif(simulation=="Powheg") :
	if(parameterName.find("Run-I-WZ")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff\nset PowhegMEqq2W2ff:Process Electron\ninsert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
	elif(parameterName.find("Run-I-W")>=0 or parameterName.find("Run-II-W")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff\nset PowhegMEqq2W2ff:Process Electron\n"
	elif(parameterName.find("Run-I-Z")>=0 or parameterName.find("Run-II-Z-e")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
	elif(parameterName.find("Run-II-Z-mu")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Muon\n"
	elif(simulation=="Matchbox" or simulation == "Dipole" ) :
	if(parameterName.find("Run-I-WZ")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p pbar e+ e-\ndo Factory:Process p pbar e+ nu\ndo Factory:Process p pbar e- nu\n"
	process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120GeV\n"
	elif(parameterName.find("Run-I-W")>=0 or parameterName.find("Run-II-W")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p pbar e+ nu\ndo Factory:Process p pbar e- nu\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	elif(parameterName.find("Run-I-Z")>=0 or parameterName.find("Run-II-Z-e")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p pbar e+ e-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120GeV\n"
	elif(parameterName.find("Run-II-Z-mu")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p pbar mu+ mu-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120GeV\n"


	# LHC
	elif(collider=="LHC") :
	if(parameterName.find("-7-")>=0) :
	process="set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 7000.0\n"
	elif(parameterName.find("-8-")>=0) :
	process="set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 8000.0\n"
	else :
	process="set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 7000.0\n"

	if(simulation=="") :
	if(parameterName.find("WZ")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEPP2VV\nset MEPP2VV:Process WZ\n"
	process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+;\n"
	process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_ebar,e-; /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
	process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("WW-emu")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEPP2VV\nset MEPP2VV:Process WW\n"
	process+="set /Herwig/Particles/W+:Synchronized 0\n"
	process+="set /Herwig/Particles/W-:Synchronized 0\n"
	process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+;\n"
	process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("WW-ll")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEPP2VV\nset MEPP2VV:Process WW\n"
	process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+; /Herwig/Particles/W+/W+->nu_tau,tau+;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("ZZ-ll")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEPP2VV\nset MEPP2VV:Process ZZ\n"
	process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("ZZ-lv")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEPP2VV\nset MEPP2VV:Process ZZ\n"
	process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+; /Herwig/Particles/Z0/Z0->nu_e,nu_ebar; /Herwig/Particles/Z0/Z0->nu_mu,nu_mubar; /Herwig/Particles/Z0/Z0->nu_tau,nu_taubar;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("LHC-W-Z-e")>0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
	process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Electron\n"
	elif(parameterName.find("LHC-W-Z-mu")>0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Muon\n"
	process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Muon\n"
	elif(parameterName.find("W-e")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Electron\n"
	elif(parameterName.find("W-mu")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Muon\n"
	elif(parameterName.find("Z-e")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
	elif(parameterName.find("Z-mu")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Muon\n"
	elif(parameterName.find("Z-LowMass-e")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
	process+="set /Herwig/Cuts/QCDCuts:MHatMin 20.GeV\nset /Herwig/Cuts/MassCut:MinM 20.GeV\nset /Herwig/Cuts/MassCut:MaxM 70.*GeV\n"
	elif(parameterName.find("Z-MedMass-e")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Electron\n"
	process+="set /Herwig/Cuts/QCDCuts:MHatMin 40.GeV\nset /Herwig/Cuts/MassCut:MinM 40.GeV\nset /Herwig/Cuts/MassCut:MaxM 130.*GeV\n"
	elif(parameterName.find("Z-LowMass-mu")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Muon\n"
	process+="set /Herwig/Cuts/QCDCuts:MHatMin 20.GeV\nset /Herwig/Cuts/MassCut:MinM 10.GeV\nset /Herwig/Cuts/MassCut:MaxM 70.*GeV\n"
	elif(parameterName.find("W-Jet")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEWJet\nset MEWJet:WDecay Electron\n"
	if(parameterName.find("W-Jet-1-e")>=0) :
	process+="set /Herwig/Cuts/WBosonKtCut:MinKT 100.0*GeV\n"
	parameterName=parameterName.replace("W-Jet-1-e","W-Jet-e")
	elif(parameterName.find("W-Jet-2-e")>=0) :
	process+="set /Herwig/Cuts/WBosonKtCut:MinKT 190.0*GeV\n"
	parameterName=parameterName.replace("W-Jet-2-e","W-Jet-e")
	elif(parameterName.find("W-Jet-3-e")>=0) :
	process+="set /Herwig/Cuts/WBosonKtCut:MinKT 270.0*GeV\n"
	parameterName=parameterName.replace("W-Jet-3-e","W-Jet-e")
	elif(parameterName.find("Z-Jet")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] MEZJet\nset MEZJet:ZDecay Electron\n"
	if(parameterName.find("Z-Jet-1-e")>=0) :
	process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 100.0*GeV\n"
	parameterName=parameterName.replace("Z-Jet-1-e","Z-Jet-e")
	elif(parameterName.find("Z-Jet-2-e")>=0) :
	process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 190.0*GeV\n"
	parameterName=parameterName.replace("Z-Jet-2-e","Z-Jet-e")
	elif(parameterName.find("Z-Jet-3-e")>=0) :
	process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 270.0*GeV\n"
	parameterName=parameterName.replace("Z-Jet-3-e","Z-Jet-e")
	else :
	logging.error(" Process %s not supported for internal matrix elements" % name)
	sys.exit(1)
	elif(simulation=="Powheg") :
	if(parameterName.find("WZ")>=0) :
	process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
	process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV\nset PowhegMEPP2VV:Process WZ\n"
	process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+;\n"
	process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_ebar,e-; /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
	process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("WW-emu")>=0) :
	process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
	process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV\nset PowhegMEPP2VV:Process WW\n"
	process+="set /Herwig/Particles/W+:Synchronized 0\n"
	process+="set /Herwig/Particles/W-:Synchronized 0\n"
	process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+;\n"
	process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("WW-ll")>=0) :
	process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
	process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV\nset PowhegMEPP2VV:Process WW\n"
	process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+; /Herwig/Particles/W+/W+->nu_tau,tau+;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("ZZ-ll")>=0) :
	process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
	process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV\nset PowhegMEPP2VV:Process ZZ\n"
	process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("ZZ-lv")>=0) :
	process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
	process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEPP2VV\nset PowhegMEPP2VV:Process ZZ\n"
	process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+; /Herwig/Particles/Z0/Z0->nu_e,nu_ebar; /Herwig/Particles/Z0/Z0->nu_mu,nu_mubar; /Herwig/Particles/Z0/Z0->nu_tau,nu_taubar;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("LHC-W-Z-e")>0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff\nset PowhegMEqq2W2ff:Process Electron\n"
	elif(parameterName.find("LHC-W-Z-mu")>0) :
	process+="insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff\nset MEqq2gZ2ff:Process Muon\n"
	process+="insert SimpleQCD:MatrixElements[0] MEqq2W2ff\nset MEqq2W2ff:Process Muon\n"
	elif(parameterName.find("W-e")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff\nset PowhegMEqq2W2ff:Process Electron\n"
	elif(parameterName.find("W-mu")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2W2ff\nset PowhegMEqq2W2ff:Process Muon\n"
	elif(parameterName.find("Z-e")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
	elif(parameterName.find("Z-mu")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Muon\n"
	elif(parameterName.find("Z-LowMass-e")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
	process+="set /Herwig/Cuts/QCDCuts:MHatMin 20.GeV\nset /Herwig/Cuts/MassCut:MinM 20.GeV\nset /Herwig/Cuts/MassCut:MaxM 70.*GeV\n"
	elif(parameterName.find("Z-MedMass-e")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Electron\n"
	process+="set /Herwig/Cuts/QCDCuts:MHatMin 40.GeV\nset /Herwig/Cuts/MassCut:MinM 40.GeV\nset /Herwig/Cuts/MassCut:MaxM 130.*GeV\n"
	elif(parameterName.find("Z-LowMass-mu")>=0) :
	process+="insert SimpleQCD:MatrixElements[0] PowhegMEqq2gZ2ff\nset PowhegMEqq2gZ2ff:Process Muon\n"
	process+="set /Herwig/Cuts/QCDCuts:MHatMin 20.GeV\nset /Herwig/Cuts/MassCut:MinM 10.GeV\nset /Herwig/Cuts/MassCut:MaxM 70.*GeV\n"
	else :
	logging.error(" Process %s not supported for internal POWHEG matrix elements" % name)
	sys.exit(1)

	elif(simulation=="Matchbox" or simulation == "Dipole" ) :
	if(parameterName.find("WZ")>=0) :
	process+="set /Herwig/Particles/W+:HardProcessWidth 0.*GeV\n"
	process+="set /Herwig/Particles/W-:HardProcessWidth 0.*GeV\n"
	process+="set /Herwig/Particles/Z0:HardProcessWidth 0.*GeV\n"
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p W+ Z0\ndo Factory:Process p p W- Z0\n"
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 171.6*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n\n"
	process+="set /Herwig/Cuts/MassCut:MinM 66GeV\nset /Herwig/Cuts/MassCut:MaxM 116GeV\n"
	process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+;\n"
	process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_ebar,e-; /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
	process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("WW-emu")>=0) :
	process+="set /Herwig/Particles/W+:HardProcessWidth 0.*GeV\n"
	process+="set /Herwig/Particles/W-:HardProcessWidth 0.*GeV\n"
	process+="set /Herwig/Particles/Z0:HardProcessWidth 0.*GeV\n"
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process pnob pnob W+ W-\n"
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 160.8*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
	process+="set /Herwig/Particles/W+:Synchronized 0\n"
	process+="set /Herwig/Particles/W-:Synchronized 0\n"
	process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+;\n"
	process+="do /Herwig/Particles/W-:SelectDecayModes /Herwig/Particles/W-/W-->nu_mubar,mu-;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("WW-ll")>=0) :
	process+="set /Herwig/Particles/W+:HardProcessWidth 0.*GeV\n"
	process+="set /Herwig/Particles/W-:HardProcessWidth 0.*GeV\n"
	process+="set /Herwig/Particles/Z0:HardProcessWidth 0.*GeV\n"
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process pnob pnob W+ W-\n"
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 160.8*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
	process+="do /Herwig/Particles/W+:SelectDecayModes /Herwig/Particles/W+/W+->nu_e,e+; /Herwig/Particles/W+/W+->nu_mu,mu+; /Herwig/Particles/W+/W+->nu_tau,tau+;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("ZZ-ll")>=0) :
	process+="set /Herwig/Particles/W+:HardProcessWidth 0.*GeV\n"
	process+="set /Herwig/Particles/W-:HardProcessWidth 0.*GeV\n"
	process+="set /Herwig/Particles/Z0:HardProcessWidth 0.*GeV\n"
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p Z0 Z0\n"
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 182.2*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
	process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("ZZ-lv")>=0) :
	process+="set /Herwig/Particles/W+:HardProcessWidth 0.*GeV\n"
	process+="set /Herwig/Particles/W-:HardProcessWidth 0.*GeV\n"
	process+="set /Herwig/Particles/Z0:HardProcessWidth 0.*GeV\n"
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p Z0 Z0\n"
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 182.2*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
	process+="do /Herwig/Particles/Z0:SelectDecayModes /Herwig/Particles/Z0/Z0->e-,e+; /Herwig/Particles/Z0/Z0->mu-,mu+; /Herwig/Particles/Z0/Z0->tau-,tau+; /Herwig/Particles/Z0/Z0->nu_e,nu_ebar; /Herwig/Particles/Z0/Z0->nu_mu,nu_mubar; /Herwig/Particles/Z0/Z0->nu_tau,nu_taubar;\n"
	process+="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	process+="insert /Herwig/EventHandlers/LHCHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	elif(parameterName.find("LHC-W-Z-e")>0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\n"
	process+="do Factory:Process p p e+ e-\ndo Factory:Process p p e+ nu\ndo Factory:Process p p e- nu\n"
	process+="set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120GeV\n"
	elif(parameterName.find("LHC-W-Z-mu")>0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\n"
	process+="do Factory:Process p p mu+ mu-\ndo Factory:Process p p mu+ nu\ndo Factory:Process p p mu- nu\n"
	process+="set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120GeV\n"
	elif(parameterName.find("W-e")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ nu\ndo Factory:Process p p e- nu\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	elif(parameterName.find("W-mu")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p mu+ nu\ndo Factory:Process p p mu- nu\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	elif(parameterName.find("Z-e")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	elif(parameterName.find("Z-mu")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p mu+ mu-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120GeV\n"
	elif(parameterName.find("Z-LowMass-e")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 20GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 70GeV\n"
	elif(parameterName.find("Z-MedMass-e")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 40GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 130GeV\n"
	elif(parameterName.find("Z-LowMass-mu")>=0) :
	process+="set Factory:OrderInAlphaS 0\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p mu+ mu-\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale\n"
	process+="set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 20GeV\nset /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 70GeV\n"
	elif(parameterName.find("W-Jet")>=0) :
	process+="set Factory:OrderInAlphaS 1\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ nu j\ndo Factory:Process p p e- nu j\n\n"
	process+="set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/HTScale\n"
	process+="set /Herwig/Cuts/QCDCuts:JetFinder /Herwig/Cuts/JetFinder\n"
	process+="insert /Herwig/Cuts/QCDCuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
	process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
	if(parameterName.find("W-Jet-1-e")>=0) :
	process+="set /Herwig/Cuts/FirstJet:PtMin 100.*GeV\n"
	parameterName=parameterName.replace("W-Jet-1-e","W-Jet-e")
	elif(parameterName.find("W-Jet-2-e")>=0) :
	process+="set /Herwig/Cuts/FirstJet:PtMin 190.0*GeV\n"
	parameterName=parameterName.replace("W-Jet-2-e","W-Jet-e")
	elif(parameterName.find("W-Jet-3-e")>=0) :
	process+="set /Herwig/Cuts/FirstJet:PtMin 270.0*GeV\n"
	parameterName=parameterName.replace("W-Jet-3-e","W-Jet-e")
	elif(parameterName.find("Z-Jet")>=0) :
	process+="set Factory:OrderInAlphaS 1\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e- j\n"
	process+="set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/HTScale\n"
	process+="set /Herwig/Cuts/QCDCuts:JetFinder /Herwig/Cuts/JetFinder\n"
	process+="insert /Herwig/Cuts/QCDCuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
	process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
	if(parameterName.find("Z-Jet-1-e")>=0) :
	process+="set /Herwig/Cuts/FirstJet:PtMin 100.*GeV\n"
	parameterName=parameterName.replace("Z-Jet-1-e","Z-Jet-e")
	elif(parameterName.find("Z-Jet-2-e")>=0) :
	process+="set /Herwig/Cuts/FirstJet:PtMin 190.0*GeV\n"
	parameterName=parameterName.replace("Z-Jet-2-e","Z-Jet-e")
	elif(parameterName.find("Z-Jet-3-e")>=0) :
	process+="set /Herwig/Cuts/FirstJet:PtMin 270.0*GeV\n"
	parameterName=parameterName.replace("Z-Jet-3-e","Z-Jet-e")
	elif(parameterName.find("Z-bb")>=0) :
	bscheme="read Matchbox/FourFlavourScheme.in"
	process+="set /Herwig/Particles/b:HardProcessMass 4.2GeV\nset /Herwig/Particles/bbar:HardProcessMass 4.2GeV\n"
	process+="set Factory:OrderInAlphaS 2\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e- b bbar\n"
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 91.2*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
	process+="set /Herwig/Cuts/MassCut:MinM 66GeV\nset /Herwig/Cuts/MassCut:MaxM 116GeV\n"
	process+="set /Herwig/Cuts/QCDCuts:JetFinder /Herwig/Cuts/JetFinder\n"
	process+="insert /Herwig/Cuts/QCDCuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
	process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
	process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/SecondJet\n"
	process+="set /Herwig/Cuts/FirstJet:PtMin 18.*GeV\n"
	process+="set /Herwig/Cuts/SecondJet:PtMin 15.*GeV\n"
	elif(parameterName.find("Z-b")>=0) :
	process+="set Factory:OrderInAlphaS 1\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ e- b\ndo Factory:Process p p e+ e- bbar\n"
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 91.2*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
	process+="set /Herwig/Cuts/MassCut:MinM 66GeV\nset /Herwig/Cuts/MassCut:MaxM 116GeV\n"
	process+="set /Herwig/Cuts/QCDCuts:JetFinder /Herwig/Cuts/JetFinder\n"
	process+="insert /Herwig/Cuts/QCDCuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
	process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
	process+="set /Herwig/Cuts/FirstJet:PtMin 15.*GeV\n"
	elif(parameterName.find("W-b")>=0) :
	bscheme="read Matchbox/FourFlavourScheme.in"
	process += "set /Herwig/Particles/b:HardProcessMass 4.2GeV\nset /Herwig/Particles/bbar:HardProcessMass 4.2GeV\n"
	process += "set Factory:OrderInAlphaS 2\nset Factory:OrderInAlphaEW 2\ndo Factory:Process p p e+ nu b bbar\ndo Factory:Process p p e- nu b bbar\n"
	process += "set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 80.4*GeV\nset Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
	process+="set /Herwig/Cuts/QCDCuts:JetFinder /Herwig/Cuts/JetFinder\n"
	process+="insert /Herwig/Cuts/QCDCuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
	process+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
	process+="set /Herwig/Cuts/FirstJet:PtMin 30.*GeV\n"
	else :
	logging.error(" Process %s not supported for Matchbox matrix elements" % name)
	sys.exit(1)
	# write the file
	if(simulation=="Matchbox" or simulation =="Dipole") :
	with open(os.path.join("Rivet",name+".in") ,'w') as f:
	f.write( template.substitute({ 'process' : process,
	'runname' : name,
	'bscheme' : bscheme,
	'parameterFile' : os.path.join(collider,collider+"-"+parameterName+".in") }))
	else :
	with open(os.path.join("Rivet",name+".in") ,'w') as f:
	f.write( template.substitute({ 'process' : process,
	'runname' : name,
	'parameterFile' : os.path.join(collider,collider+"-"+parameterName+".in") }))

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:21 PM (1 d, 13 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805128
Default Alt Text: (114 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions