No OneTemporary
Actions

Size

99 KB

Subscribers

None

View Options

	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,525 +1,525 @@
	// -- C++ --
	//
	// QTildeMatching.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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 < max(theQTildeSudakov->pT2min(),sqr(safeCut()) ))
	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) >=
	max(theQTildeSudakov->pT2min(),sqr(safeCut()));
	else
	return sqr(z(1.-z)qtilde) - sqr((1.-z)mu) - zsqr(Qg) >=
	max(theQTildeSudakov->pT2min(),sqr(safeCut()));
	}
	if ( dipole()->bornEmitter() < 2 ) {
	return
	sqr((1.-z)qtilde) - zsqr(Qg) >=
	max(theQTildeSudakov->pT2min(),sqr(safeCut()));
	}
	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() /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()]);
	}
	// allow small violations (numerical inaccuracies)
	if ( z <= 0 && z >= -1e-6 ) {
	z = std::numeric_limits<double>::epsilon();
	} else if ( z >= 1 && z <= 1+1e-6 ) {
	z = 1-std::numeric_limits<double>::epsilon();
	}

	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);
	}
	if ( qtilde2 < ZERO ) {
	qtilde2 = ZERO;
	}

	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();
	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
	+ static Reference<QTildeMatching,PartnerFinder> interfaceQTildeFinder
	("QTildeFinder",
	"Set the partner finder to calculate hard scales.",
	&QTildeMatching::theQTildeFinder, false, false, true, false, false);
	interfaceQTildeFinder.rank(-1);

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

	static Reference<QTildeMatching,ShowerHandler> interfaceShowerHandler
	("ShowerHandler",
	"The QTilde shower handler to use.",
	&QTildeMatching::theShowerHandler, false, false, true, true, false);
	interfaceShowerHandler.rank(-1);

	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);
	interfaceCorrectForXZMismatch.rank(-1);

	}

	diff --git a/MatrixElement/Matchbox/Matching/QTildeMatching.h b/MatrixElement/Matchbox/Matching/QTildeMatching.h
	--- a/MatrixElement/Matchbox/Matching/QTildeMatching.h
	+++ b/MatrixElement/Matchbox/Matching/QTildeMatching.h
	@@ -1,201 +1,201 @@
	// -- C++ --
	//
	// QTildeMatching.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef Herwig_QTildeMatching_H
	#define Herwig_QTildeMatching_H
	//
	// This is the declaration of the QTildeMatching class.
	//

	#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
	#include "Herwig/Shower/ShowerHandler.h"
	-#include "Herwig/Shower/QTilde/Default/QTildeFinder.h"
	+#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
	#include "Herwig/Shower/QTilde/Default/QTildeSudakov.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief QTildeMatching implements NLO matching with the default shower.
	*
	*/
	class QTildeMatching: public Herwig::ShowerApproximation {

	public:

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

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

	public:

	/**
	* Return the shower approximation to the real emission cross
	* section for the given pair of Born and real emission
	* configurations.
	*/
	virtual CrossSection dSigHatDR() const;

	/**
	* Return the shower approximation splitting kernel for the given
	* pair of Born and real emission configurations in units of the
	* Born center of mass energy squared, and including a weight to
	* project onto the splitting given by the dipole used.
	*/
	virtual double me2() const;

	/**
	* Determine if the configuration is below or above the cutoff.
	*/
	virtual void checkCutoff();

	/**
	* Determine all kinematic variables which are not provided by the
	* dipole kinematics; store all shower variables in the respective
	* dipole object for later use.
	*/
	virtual void getShowerVariables();

	protected:

	/**
	* Return true, if the shower was able to generate an emission
	* leading from the given Born to the given real emission process.
	*/
	virtual bool isInShowerPhasespace() const;

	/**
	* Return true, if the shower emission leading from the given Born
	* to the given real emission process would have been generated
	* above the shower's infrared cutoff.
	*/
	virtual bool isAboveCutoff() const;

	/**
	* Calculate qtilde^2 and z for the splitting considered
	*/
	void calculateShowerVariables() const;

	/**
	* Return the splitting function as a function of the kinematic
	* variables
	*/
	double splitFn(const pair<Energy2,double>&) const;

	public:

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

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

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

	protected:

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

	protected:

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

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

	private:

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

	/**
	* The shower handler to be used
	*/
	Ptr<ShowerHandler>::ptr theShowerHandler;

	/**
	* The qtilde partner finder for calculating the hard scales
	*/
	- Ptr<QTildeFinder>::ptr theQTildeFinder;
	+ Ptr<PartnerFinder>::ptr theQTildeFinder;

	/**
	* The qtilde Sudakov to access the cutoff
	*/
	Ptr<QTildeSudakov>::ptr theQTildeSudakov;

	/**
	* True, if PDF weight should be corrected for z/x mismatch at the
	* hard phase space boundary
	*/
	bool theCorrectForXZMismatch;

	};

	}

	#endif /* Herwig_QTildeMatching_H */
	diff --git a/Shower/QTilde/Base/PartnerFinder.cc b/Shower/QTilde/Base/PartnerFinder.cc
	--- a/Shower/QTilde/Base/PartnerFinder.cc
	+++ b/Shower/QTilde/Base/PartnerFinder.cc
	@@ -1,536 +1,757 @@
	// -- C++ --
	//
	// PartnerFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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 PartnerFinder class.
	//

	#include "PartnerFinder.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/Utilities/DescribeClass.h"

	using namespace Herwig;

	-DescribeAbstractClass<PartnerFinder,Interfaced>
	+DescribeClass<PartnerFinder,Interfaced>
	describePartnerFinder ("Herwig::PartnerFinder","HwShower.so");

	// some useful functions to avoid using #define
	namespace {

	// return bool if final-state particle
	inline bool FS(const tShowerParticlePtr a) {
	return a->isFinalState();
	}

	// return colour line pointer
	inline Ptr<ThePEG::ColourLine>::transient_pointer
	CL(const tShowerParticlePtr a, unsigned int index=0) {
	return a->colourInfo()->colourLines().empty() ? ThePEG::tColinePtr() :
	const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->colourLines()[index]);
	}

	// return colour line size
	inline unsigned int
	CLSIZE(const tShowerParticlePtr a) {
	return a->colourInfo()->colourLines().size();
	}

	inline Ptr<ThePEG::ColourLine>::transient_pointer
	ACL(const tShowerParticlePtr a, unsigned int index=0) {
	return a->colourInfo()->antiColourLines().empty() ? ThePEG::tColinePtr() :
	const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->antiColourLines()[index]);
	}

	inline unsigned int
	ACLSIZE(const tShowerParticlePtr a) {
	return a->colourInfo()->antiColourLines().size();
	}
	}

	void PartnerFinder::persistentOutput(PersistentOStream & os) const {
	- os << partnerMethod_ << QEDPartner_ << scaleChoice_;
	+ os << partnerMethod_ << QEDPartner_ << scaleChoice_
	+ << _finalFinalConditions << _initialFinalDecayConditions
	+ << _initialInitialConditions;
	}

	void PartnerFinder::persistentInput(PersistentIStream & is, int) {
	- is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_;
	+ is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_
	+ >> _finalFinalConditions >> _initialFinalDecayConditions
	+ >>_initialInitialConditions;
	}

	void PartnerFinder::Init() {

	static ClassDocumentation<PartnerFinder> documentation
	("This class is responsible for finding the partners for each interaction types ",
	"and within the evolution scale range specified by the ShowerVariables ",
	"then to determine the initial evolution scales for each pair of partners.");

	static Switch<PartnerFinder,int> interfacePartnerMethod
	("PartnerMethod",
	"Choice of partner finding method for gluon evolution.",
	&PartnerFinder::partnerMethod_, 0, false, false);
	static SwitchOption interfacePartnerMethodRandom
	(interfacePartnerMethod,
	"Random",
	"Choose partners of a gluon randomly.",
	0);
	static SwitchOption interfacePartnerMethodMaximum
	(interfacePartnerMethod,
	"Maximum",
	"Choose partner of gluon with largest angle.",
	1);

	static Switch<PartnerFinder,int> interfaceQEDPartner
	("QEDPartner",
	"Control of which particles to use as the partner for QED radiation",
	&PartnerFinder::QEDPartner_, 0, false, false);
	static SwitchOption interfaceQEDPartnerAll
	(interfaceQEDPartner,
	"All",
	"Consider all possible choices which give a positive contribution"
	" in the soft limit.",
	0);
	static SwitchOption interfaceQEDPartnerIIandFF
	(interfaceQEDPartner,
	"IIandFF",
	"Only allow initial-initial or final-final combinations",
	1);
	static SwitchOption interfaceQEDPartnerIF
	(interfaceQEDPartner,
	"IF",
	"Only allow initial-final combinations",
	2);

	static Switch<PartnerFinder,int> interfaceScaleChoice
	("ScaleChoice",
	"The choice of the evolution scales",
	&PartnerFinder::scaleChoice_, 0, false, false);
	static SwitchOption interfaceScaleChoicePartner
	(interfaceScaleChoice,
	"Partner",
	"Scale of all interactions is that of the evolution partner",
	0);
	static SwitchOption interfaceScaleChoiceDifferent
	(interfaceScaleChoice,
	"Different",
	"Allow each interaction to have different scales",
	1);

	+
	+ static Switch<PartnerFinder,unsigned int> interfaceFinalFinalConditions
	+ ("FinalFinalConditions",
	+ "The initial conditions for the shower of a final-final colour connection",
	+ &PartnerFinder::_finalFinalConditions, 0, false, false);
	+ static SwitchOption interfaceFinalFinalConditionsSymmetric
	+ (interfaceFinalFinalConditions,
	+ "Symmetric",
	+ "The symmetric choice",
	+ 0);
	+ static SwitchOption interfaceFinalFinalConditionsColoured
	+ (interfaceFinalFinalConditions,
	+ "Coloured",
	+ "Maximal radiation from the coloured particle",
	+ 1);
	+ static SwitchOption interfaceFinalFinalConditionsAntiColoured
	+ (interfaceFinalFinalConditions,
	+ "AntiColoured",
	+ "Maximal emission from the anticoloured particle",
	+ 2);
	+ static SwitchOption interfaceFinalFinalConditionsRandom
	+ (interfaceFinalFinalConditions,
	+ "Random",
	+ "Randomly selected maximal emission from one of the particles",
	+ 3);
	+
	+ static Switch<PartnerFinder,unsigned int> interfaceInitialFinalDecayConditions
	+ ("InitialFinalDecayConditions",
	+ "The initial conditions for the shower of an initial-final"
	+ " decay colour connection.",
	+ &PartnerFinder::_initialFinalDecayConditions, 0, false, false);
	+ static SwitchOption interfaceInitialFinalDecayConditionsSymmetric
	+ (interfaceInitialFinalDecayConditions,
	+ "Symmetric",
	+ "The symmetric choice",
	+ 0);
	+ static SwitchOption interfaceInitialFinalDecayConditionsMaximal
	+ (interfaceInitialFinalDecayConditions,
	+ "Maximal",
	+ "Maximal radiation from the decay product",
	+ 1);
	+ static SwitchOption interfaceInitialFinalDecayConditionsSmooth
	+ (interfaceInitialFinalDecayConditions,
	+ "Smooth",
	+ "Smooth matching in the soft limit",
	+ 2);
	+
	+ static Switch<PartnerFinder,unsigned int> interfaceInitialInitialConditions
	+ ("InitialInitialConditions",
	+ "The initial conditions for the shower of an initial-initial"
	+ " colour connection.",
	+ &PartnerFinder::_initialInitialConditions, 0, false, false);
	+ static SwitchOption interfaceInitialInitialConditionsSymmetric
	+ (interfaceInitialInitialConditions,
	+ "Symmetric",
	+ "The symmetric choice",
	+ 0);
	+ static SwitchOption interfaceInitialInitialConditionsMaximiseB
	+ (interfaceInitialInitialConditions,
	+ "MaximiseB",
	+ "Maximal radiation from parton b",
	+ 1);
	+ static SwitchOption interfaceInitialInitialConditionsMaximiseC
	+ (interfaceInitialInitialConditions,
	+ "MaximiseC",
	+ "Maximal radiation from parton c",
	+ 2);
	}

	void PartnerFinder::setInitialEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	ShowerInteraction type,
	const bool setPartners) {
	// clear the existing partners
	for(ShowerParticleVector::const_iterator cit = particles.begin();
	cit != particles.end(); ++cit) (*cit)->clearPartners();
	// set them
	if(type==ShowerInteraction::QCD) {
	setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
	}
	else if(type==ShowerInteraction::QED) {
	setInitialQEDEvolutionScales(particles,isDecayCase,setPartners);
	}
	else {
	setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
	setInitialQEDEvolutionScales(particles,isDecayCase,false);
	}
	// print out for debugging
	if(Debug::level>=10) {
	for(ShowerParticleVector::const_iterator cit = particles.begin();
	cit != particles.end(); ++cit) {
	generator()->log() << "Particle: " << **cit << "\n";
	if(!(**cit).partner()) continue;
	generator()->log() << "Primary partner: " << (*cit).partner() << "\n";
	for(vector<ShowerParticle::EvolutionPartner>::const_iterator it= (**cit).partners().begin();
	it!=(**cit).partners().end();++it) {
	generator()->log() << static_cast<long>(it->type) << " "
	<< it->weight << " "
	<< it->scale/GeV << " "
	<< *(it->partner)
	<< "\n";
	}
	}
	generator()->log() << flush;
	}
	}

	void PartnerFinder::setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	const bool setPartners) {
	// set the partners and the scales
	if(setPartners) {
	// Loop over particles and consider only coloured particles which don't
	// have already their colour partner fixed and that don't have children
	// (the latter requirement is relaxed in the case isDecayCase is true).
	// Build a map which has as key one of these particles (i.e. a pointer
	// to a ShowerParticle object) and as a corresponding value the vector
	// of all its possible normal candidate colour partners, defined as follows:
	// --- have colour, and no children (this is not required in the case
	// isDecayCase is true);
	// --- if both are initial (incoming) state particles, then the (non-null) colourLine()
	// of one of them must match the (non-null) antiColourLine() of the other.
	// --- if one is an initial (incoming) state particle and the other is
	// a final (outgoing) state particle, then both must have the
	// same (non-null) colourLine() or the same (non-null) antiColourLine();
	// Notice that this definition exclude the special case of baryon-violating
	// processes (as in R-parity Susy), which will show up as particles
	// without candidate colour partners, and that we will be treated a part later
	// (this means that no modifications of the following loop is needed!)
	ShowerParticleVector::const_iterator cit, cjt;
	for(cit = particles.begin(); cit != particles.end(); ++cit) {
	// Skip colourless particles
	if(!(*cit)->data().coloured()) continue;
	// find the partners
	vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners =
	findQCDPartners(*cit,particles);
	// must have a partner
	if(partners.empty()) {
	throw Exception() << "`Failed to make colour connections in "
	<< "PartnerFinder::setQCDInitialEvolutionScales"
	<< (**cit)
	<< Exception::eventerror;
	}
	// Calculate the evolution scales for all possible pairs of of particles
	vector<pair<Energy,Energy> > scales;
	for(unsigned int ix=0;ix< partners.size();++ix) {
	scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
	isDecayCase));
	}
	// In the case of more than one candidate colour partners,
	// there are now two approaches to choosing the partner. The
	// first method is based on two assumptions:
	// 1) the choice of which is the colour partner is done
	// randomly between the available candidates.
	// 2) the choice of which is the colour partner is done
	// independently from each particle: in other words,
	// if for a particle "i" its selected colour partner is
	// the particle "j", then the colour partner of "j"
	// does not have to be necessarily "i".
	// The second method always chooses the furthest partner
	// for hard gluons and gluinos.
	// store the choice
	int position(-1);
	// random choice
	if( partnerMethod_ == 0 ) {
	// random choice of partner
	position = UseRandom::irnd(partners.size());
	}
	// take the one with largest angle
	else if (partnerMethod_ == 1 ) {
	if ((*cit)->perturbative() == 1 &&
	(*cit)->dataPtr()->iColour()==PDT::Colour8 ) {
	assert(partners.size()==2);
	// Determine largest angle
	double maxAngle(0.);
	for(unsigned int ix=0;ix<partners.size();++ix) {
	double angle = (*cit)->momentum().vect().
	angle(partners[ix].second->momentum().vect());
	if(angle>maxAngle) {
	maxAngle = angle;
	position = ix;
	}
	}
	}
	else
	position = UseRandom::irnd(partners.size());
	}
	else
	assert(false);
	// set the evolution partner
	(*cit)->partner(partners[position].second);
	for(unsigned int ix=0;ix<partners.size();++ix) {
	(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
	1.,partners[ix].first,
	scales[ix].first));
	}
	// set scales for all interactions to that of the partner, default
	Energy scale = scales[position].first;
	for(unsigned int ix=0;ix<partners.size();++ix) {
	if(partners[ix].first==ShowerPartnerType::QCDColourLine) {
	(**cit).scales().QCD_c =
	(**cit).scales().QCD_c_noAO =
	(scaleChoice_==0 ? scale : scales[ix].first);
	}
	else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
	(**cit).scales().QCD_ac =
	(**cit).scales().QCD_ac_noAO =
	(scaleChoice_==0 ? scale : scales[ix].first);
	}
	else
	assert(false);
	}
	}
	}
	// primary partner set, set the others and do the scale
	else {
	for(ShowerParticleVector::const_iterator cit = particles.begin();
	cit != particles.end(); ++cit) {
	// Skip colourless particles
	if(!(*cit)->data().coloured()) continue;
	// find the partners
	vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners =
	findQCDPartners(*cit,particles);
	// must have a partner
	if(partners.empty()) {
	throw Exception() << "`Failed to make colour connections in "
	<< "PartnerFinder::setQCDInitialEvolutionScales"
	<< (**cit)
	<< Exception::eventerror;
	}
	// Calculate the evolution scales for all possible pairs of of particles
	vector<pair<Energy,Energy> > scales;
	int position(-1);
	for(unsigned int ix=0;ix< partners.size();++ix) {
	if(partners[ix].second) position = ix;
	scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
	isDecayCase));
	}
	assert(position>=0);
	for(unsigned int ix=0;ix<partners.size();++ix) {
	(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
	1.,partners[ix].first,
	scales[ix].first));
	}
	// set scales for all interactions to that of the partner, default
	Energy scale = scales[position].first;
	for(unsigned int ix=0;ix<partners.size();++ix) {
	if(partners[ix].first==ShowerPartnerType::QCDColourLine) {
	(**cit).scales().QCD_c =
	(**cit).scales().QCD_c_noAO =
	(scaleChoice_==0 ? scale : scales[ix].first);
	}
	else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
	(**cit).scales().QCD_ac =
	(**cit).scales().QCD_ac_noAO =
	(scaleChoice_==0 ? scale : scales[ix].first);
	}
	else {
	assert(false);
	}
	}
	}
	}
	}

	void PartnerFinder::setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	const bool setPartners) {
	// loop over all the particles
	for(ShowerParticleVector::const_iterator cit = particles.begin();
	cit != particles.end(); ++cit) {
	// not charged or photon continue
	if(!(**cit).dataPtr()->charged()) continue;
	// find the potential partners
	vector<pair<double,tShowerParticlePtr> > partners = findQEDPartners(*cit,particles,isDecayCase);
	if(partners.empty()) {
	throw Exception() << "Failed to find partner in "
	<< "PartnerFinder::setQEDInitialEvolutionScales"
	<< (**cit) << Exception::eventerror;
	}
	// calculate the probabilities
	double prob(0.);
	for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
	// normalise
	for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
	// set the partner if required
	int position(-1);
	// use QCD partner if set
	if(!setPartners&&(*cit)->partner()) {
	for(unsigned int ix=0;ix<partners.size();++ix) {
	if((*cit)->partner()==partners[ix].second) {
	position = ix;
	break;
	}
	}
	}
	// set the partner
	if(setPartners\|\|!(*cit)->partner()\|\|position<0) {
	prob = UseRandom::rnd();
	for(unsigned int ix=0;ix<partners.size();++ix) {
	if(partners[ix].first>prob) {
	position = ix;
	break;
	}
	prob -= partners[ix].first;
	}
	if(position>=0&&(setPartners\|\|!(*cit)->partner())) {
	(*cit)->partner(partners[position].second);
	}
	}
	// must have a partner
	if(position<0) throw Exception() << "Failed to find partner in "
	<< "PartnerFinder::setQEDInitialEvolutionScales"
	<< (**cit) << Exception::eventerror;
	// Calculate the evolution scales for all possible pairs of of particles
	vector<pair<Energy,Energy> > scales;
	for(unsigned int ix=0;ix< partners.size();++ix) {
	scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
	isDecayCase));
	}
	// store all the possible partners
	for(unsigned int ix=0;ix<partners.size();++ix) {
	(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
	partners[ix].first,
	ShowerPartnerType::QED,
	scales[ix].first));
	}
	// set scales
	(**cit).scales().QED = scales[position].first;
	(**cit).scales().QED_noAO = scales[position].first;
	}
	}

	pair<Energy,Energy> PartnerFinder::
	calculateInitialEvolutionScales(const ShowerPPair &particlePair,
	const bool isDecayCase) {
	bool FS1=FS(particlePair.first),FS2= FS(particlePair.second);

	- if(FS1 && FS2)
	- return calculateFinalFinalScales(particlePair);
	+ if(FS1 && FS2) {
	+ bool colouredFirst =
	+ particlePair.first->colourLine()&&
	+ particlePair.first->colourLine()==particlePair.second->antiColourLine();
	+ return calculateFinalFinalScales(particlePair.first->momentum(),
	+ particlePair.second->momentum(),
	+ colouredFirst);
	+ }
	else if(FS1 && !FS2) {
	ShowerPPair a(particlePair.second, particlePair.first);
	- pair<Energy,Energy> rval = calculateInitialFinalScales(a,isDecayCase);
	+ pair<Energy,Energy> rval = calculateInitialFinalScales(a.first->momentum(),
	+ a.second->momentum(),
	+ isDecayCase);
	return pair<Energy,Energy>(rval.second,rval.first);
	}
	else if(!FS1 &&FS2)
	- return calculateInitialFinalScales(particlePair,isDecayCase);
	+ return calculateInitialFinalScales(particlePair.first->momentum(),particlePair.second->momentum(),isDecayCase);
	else
	- return calculateInitialInitialScales(particlePair);
	+ return calculateInitialInitialScales(particlePair.first->momentum(),particlePair.second->momentum());
	}

	vector< pair<ShowerPartnerType, tShowerParticlePtr> >
	PartnerFinder::findQCDPartners(tShowerParticlePtr particle,
	const ShowerParticleVector &particles) {
	vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners;
	ShowerParticleVector::const_iterator cjt;
	for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
	if(!(cjt)->data().coloured() \|\| particle==cjt) continue;
	// one initial-state and one final-state particle
	if(FS(particle) != FS(*cjt)) {
	// loop over all the colours of both particles
	for(unsigned int ix=0; ix<CLSIZE(particle); ++ix) {
	for(unsigned int jx=0; jx<CLSIZE(*cjt); ++jx) {
	if((CL(particle,ix) && CL(particle,ix)==CL(*cjt,jx))) {
	partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
	}
	}
	}
	//loop over all the anti-colours of both particles
	for(unsigned int ix=0; ix<ACLSIZE(particle); ++ix) {
	for(unsigned int jx=0; jx<ACLSIZE(*cjt); ++jx) {
	if((ACL(particle,ix) && ACL(particle,ix)==ACL(*cjt,jx))) {
	partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
	}
	}
	}
	}
	// two initial-state or two final-state particles
	else {
	//loop over the colours of the first particle and the anti-colours of the other
	for(unsigned int ix=0; ix<CLSIZE(particle); ++ix){
	for(unsigned int jx=0; jx<ACLSIZE(*cjt); ++jx){
	if(CL(particle,ix) && CL(particle,ix)==ACL(*cjt,jx)) {
	partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
	}
	}
	}
	//loop over the anti-colours of the first particle and the colours of the other
	for(unsigned int ix=0; ix<ACLSIZE(particle); ++ix){
	for(unsigned int jx=0; jx<CLSIZE(*cjt); jx++){
	if(ACL(particle,ix) && ACL(particle,ix)==CL(*cjt,jx)) {
	partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
	}
	}
	}
	}
	}
	// if we haven't found any partners look for RPV
	if (partners.empty()) {
	// special for RPV
	tColinePtr col = CL(particle);
	if(FS(particle)&&col&&col->sourceNeighbours().first) {
	tColinePair cpair = col->sourceNeighbours();
	for(cjt=particles.begin();cjt!=particles.end();++cjt) {
	if(( FS(cjt) && ( CL(cjt) == cpair.first \|\| CL(*cjt) == cpair.second))\|\|
	(!FS(cjt) && (ACL(cjt) == cpair.first \|\| ACL(*cjt) == cpair.second ))) {
	partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
	}
	}
	}
	else if(col&&col->sinkNeighbours().first) {
	tColinePair cpair = col->sinkNeighbours();
	for(cjt=particles.begin();cjt!=particles.end();++cjt) {
	if(( FS(cjt) && (ACL(cjt) == cpair.first \|\| ACL(*cjt) == cpair.second))\|\|
	(!FS(cjt) && ( CL(cjt) == cpair.first \|\| CL(*cjt) == cpair.second))) {
	partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
	}
	}
	}
	col = ACL(particle);
	if(FS(particle)&&col&&col->sinkNeighbours().first) {
	tColinePair cpair = col->sinkNeighbours();
	for(cjt=particles.begin();cjt!=particles.end();++cjt) {
	if(( FS(cjt) && (ACL(cjt) == cpair.first \|\| ACL(*cjt) == cpair.second))\|\|
	(!FS(cjt) && ( CL(cjt) == cpair.first \|\| CL(*cjt) == cpair.second ))) {
	partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
	}
	}
	}
	else if(col&&col->sourceNeighbours().first) {
	tColinePair cpair = col->sourceNeighbours();
	for(cjt=particles.begin();cjt!=particles.end();++cjt) {
	if(( FS(cjt) && ( CL(cjt) == cpair.first \|\| CL(*cjt) == cpair.second))\|\|
	(!FS(cjt) && (ACL(cjt) == cpair.first \|\|ACL(*cjt) == cpair.second))) {
	partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
	}
	}
	}
	}
	// return the partners
	return partners;
	}

	vector< pair<double, tShowerParticlePtr> >
	PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
	const ShowerParticleVector &particles,
	const bool isDecayCase) {
	vector< pair<double, tShowerParticlePtr> > partners;
	ShowerParticleVector::const_iterator cjt;
	double pcharge = particle->id()==ParticleID::gamma ? 1 : double(particle->data().iCharge());
	vector< pair<double, tShowerParticlePtr> > photons;
	for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
	if(particle == *cjt) continue;
	if((*cjt).id()==ParticleID::gamma) photons.push_back(make_pair(1.,cjt));
	if(!(*cjt)->data().charged() ) continue;
	double charge = pchargedouble((cjt)->data().iCharge());
	if( FS(particle) != FS(cjt) ) charge =-1.;
	if( QEDPartner_ != 0 && !isDecayCase ) {
	// only include II and FF as requested
	if( QEDPartner_ == 1 && FS(particle) != FS(*cjt) )
	continue;
	// only include IF is requested
	else if(QEDPartner_ == 2 && FS(particle) == FS(*cjt) )
	continue;
	}
	if(particle->id()==ParticleID::gamma) charge = -abs(charge);
	// only keep positive dipoles
	if(charge<0.) partners.push_back(make_pair(-charge,*cjt));
	}
	if(particle->id()==ParticleID::gamma&& partners.empty()) {
	return photons;
	}
	return partners;
	}
	+
	+pair<Energy,Energy>
	+PartnerFinder::calculateFinalFinalScales(
	+ const Lorentz5Momentum & p1,
	+ const Lorentz5Momentum & p2,
	+ bool colouredFirst)
	+{
	+ static const double eps=1e-7;
	+ // Using JHEP 12(2003)045 we find that we need ktilde = 1/2(1+b-c+lambda)
	+ // ktilde = qtilde^2/Q^2 therefore qtilde = sqrt(ktilde*Q^2)
	+ // find momenta in rest frame of system
	+ // calculate quantities for the scales
	+ Energy2 Q2 = (p1+p2).m2();
	+ double b = p1.mass2()/Q2;
	+ double c = p2.mass2()/Q2;
	+ if(b<0.) {
	+ if(b<-eps) {
	+ throw Exception() << "Negative Mass squared b = " << b
	+ << "in PartnerFinder::calculateFinalFinalScales()"
	+ << Exception::eventerror;
	+ }
	+ b = 0.;
	+ }
	+ if(c<0.) {
	+ if(c<-eps) {
	+ throw Exception() << "Negative Mass squared c = " << c
	+ << "in PartnerFinder::calculateFinalFinalScales()"
	+ << Exception::eventerror;
	+ }
	+ c = 0.;
	+ }
	+ // KMH & PR - 16 May 2008 - swapped lambda calculation from
	+ // double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)),
	+ // which should be identical for p1 & p2 onshell in their COM
	+ // but in the inverse construction for the Nason method, this
	+ // was not the case, leading to misuse.
	+ double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c))
	+ (sqrt(b)-1.-sqrt(c))(sqrt(c)-1.-sqrt(b)));
	+ // symmetric case
	+ unsigned int iopt=finalFinalConditions();
	+ Energy firstQ,secondQ;
	+ if(iopt==0) {
	+ firstQ = sqrt(0.5Q2(1.+b-c+lam));
	+ secondQ = sqrt(0.5Q2(1.-b+c+lam));
	+ }
	+ // assymetric choice
	+ else {
	+ double kappab,kappac;
	+ // calculate kappa with coloured line getting maximum
	+ if((iopt==1&&colouredFirst)\|\| // first particle coloured+maximal for coloured
	+ (iopt==2&&!colouredFirst)\|\| // first particle anticoloured+maximal for acoloured
	+ (iopt==3&&UseRandom::rndbool(0.5))) { // random choice
	+ kappab=4.(1.-2.sqrt(c)-b+c);
	+ kappac=c+0.25*sqr(1.-b-c+lam)/(kappab-b);
	+ }
	+ else {
	+ kappac=4.(1.-2.sqrt(b)-c+b);
	+ kappab=b+0.25*sqr(1.-b-c+lam)/(kappac-c);
	+ }
	+ // calculate the scales
	+ firstQ = sqrt(Q2*kappab);
	+ secondQ = sqrt(Q2*kappac);
	+ }
	+ return pair<Energy,Energy>(firstQ, secondQ);
	+}
	+
	+
	+pair<Energy,Energy>
	+PartnerFinder::calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
	+ const bool isDecayCase) {
	+ if(!isDecayCase) {
	+ // In this case from JHEP 12(2003)045 we find the conditions
	+ // ktilde_b = (1+c) and ktilde_c = (1+2c)
	+ // We also find that c = m_c^2/Q^2. The process is a+b->c where
	+ // particle a is not colour connected (considered as a colour singlet).
	+ // Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and
	+ // q_c = sqrt(Q^2+2 m_c^2)
	+ // We also assume that the first particle in the pair is the initial
	+ // state particle and the second is the final state one c
	+ Energy2 mc2 = sqr(pc.mass());
	+ Energy2 Q2 = -(pb-pc).m2();
	+ return pair<Energy,Energy>(sqrt(Q2+mc2), sqrt(Q2+2*mc2));
	+ }
	+ else {
	+ // In this case from JHEP 12(2003)045 we find, for the decay
	+ // process b->c+a(neutral), the condition
	+ // (ktilde_b-1)(ktilde_c-c)=(1/4)sqr(1-a+c+lambda).
	+ // We also assume that the first particle in the pair is the initial
	+ // state particle (b) and the second is the final state one (c).
	+ // - We find maximal phase space coverage through emissions from
	+ // c if we set ktilde_c = 4.*(sqr(1.-sqrt(a))-c)
	+ // - We find the most 'symmetric' way to populate the phase space
	+ // occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda)
	+ // - We find the most 'smooth' way to populate the phase space
	+ // occurs for...
	+ Energy2 mb2(sqr(pb.mass()));
	+ double a=(pb-pc).m2()/mb2;
	+ double c=sqr(pc.mass())/mb2;
	+ double lambda = 1. + aa + cc - 2.a - 2.c - 2.ac;
	+ lambda = sqrt(max(lambda,0.));
	+ double PROD = 0.25*sqr(1. - a + c + lambda);
	+ double ktilde_b, ktilde_c,cosi(0.);
	+ switch(initialFinalDecayConditions()) {
	+ case 0: // the 'symmetric' choice
	+ ktilde_c = 0.5*(1-a+c+lambda) + c ;
	+ ktilde_b = 1.+PROD/(ktilde_c-c) ;
	+ break;
	+ case 1: // the 'maximal' choice
	+ ktilde_c = 4.0*(sqr(1.-sqrt(a))-c);
	+ ktilde_b = 1.+PROD/(ktilde_c-c) ;
	+ break;
	+ case 2: // the 'smooth' choice
	+ // c is a problem if very small here use 1GeV as minimum
	+ c = max(c,1.*GeV2/mb2);
	+ cosi = (sqr(1-sqrt(c))-a)/lambda;
	+ ktilde_b = 2.0/(1.0-cosi);
	+ ktilde_c = (1.0-a+c+lambda)(1.0+c-a-lambdacosi)/(2.0*(1.0+cosi));
	+ break;
	+ default:
	+ throw Exception() << "Invalid option for decay shower's phase space"
	+ << " PartnerFinder::calculateInitialFinalScales"
	+ << Exception::abortnow;
	+ }
	+ return pair<Energy,Energy>(sqrt(mb2ktilde_b),sqrt(mb2ktilde_c));
	+ }
	+}
	+
	+pair<Energy,Energy>
	+PartnerFinder::calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2) {
	+ // This case is quite simple. From JHEP 12(2003)045 we find the condition
	+ // that ktilde_b = ktilde_c = 1. In this case we have the process
	+ // b+c->a so we need merely boost to the CM frame of the two incoming
	+ // particles and then qtilde is equal to the energy in that frame
	+ Energy Q = sqrt((p1+p2).m2());
	+ if(_initialInitialConditions==1) {
	+ return pair<Energy,Energy>(sqrt(2.0)Q,sqrt(0.5)Q);
	+ } else if(_initialInitialConditions==2) {
	+ return pair<Energy,Energy>(sqrt(0.5)Q,sqrt(2.0)Q);
	+ } else {
	+ return pair<Energy,Energy>(Q,Q);
	+ }
	+}
	diff --git a/Shower/QTilde/Base/PartnerFinder.h b/Shower/QTilde/Base/PartnerFinder.h
	--- a/Shower/QTilde/Base/PartnerFinder.h
	+++ b/Shower/QTilde/Base/PartnerFinder.h
	@@ -1,208 +1,297 @@
	// -- C++ --
	//
	// PartnerFinder.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_PartnerFinder_H
	#define HERWIG_PartnerFinder_H
	//
	// This is the declaration of the PartnerFinder class.
	//
	#include "Herwig/Shower/QTilde/ShowerConfig.h"
	#include "ThePEG/Interface/Interfaced.h"
	#include "PartnerFinder.fh"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* typedef of a pair of particle for calculating the evolution scales
	*/
	typedef pair<tShowerParticlePtr,tShowerParticlePtr> ShowerPPair;

	/** \ingroup Shower
	*
	* This class is responsible of two related tasks:
	* - it finds the partners
	* - for each pair of partners (and interaction therefore)
	* it sets the initial evolution scales of both of them.
	*
	* In general the finding of the partners is performed by this class but
	* the calculation of the initial evolution scales should be implemented
	* for different shower evolution models in classes inheriting from this one.
	* Notice that a given particle has, in general, a different partner
	* for each different interaction; however, given a partner, its
	* initial evolution scale, Q, is purely a kinematical relationship
	* between the pair, without dependence on the dynamics (i.e. type of interaction).
	*
	* @see \ref PartnerFinderInterfaces "The interfaces"
	* defined for PartnerFinder.
	*/
	class PartnerFinder: public Interfaced {

	public:

	/**
	* The default constructor.
	*/
	- PartnerFinder() : partnerMethod_(0), QEDPartner_(0), scaleChoice_(0) {}
	+ PartnerFinder() : partnerMethod_(0), QEDPartner_(0), scaleChoice_(0),
	+ _finalFinalConditions(0),
	+ _initialFinalDecayConditions(0),
	+ _initialInitialConditions(0) {}

	/**
	* Given in input a collection of particles (ShowerParticle objects),
	* each of these methods set the initial evolution scales of those particles,
	* between the ones given in input, that do not have yet their
	* evolution scale set.
	* The input collection of particles can be either the full collection of
	* showering particles (kept in the main class ShowerHandler,
	* in the case isDecayCase is false, or simply, in the case isDecayCase
	* is true, the decaying particle and its decay products.
	* The methods returns true, unless something wrong (inconsistencies,
	* or undefined values) happens.
	*
	* These methods are virtual but in most cases inheriting classes should not
	* need to overide them as they simply find the relevant partner and call
	* one of the calculateScale members to calculate the scale.
	*/
	//@{
	/**
	* Set the initial scales
	* @param particles The particles to be considered
	* @param isDecayCase Whether or not this is a decay
	* @param setPartners Whether to set the colour partners or just the scales
	*/
	virtual void setInitialEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	ShowerInteraction,
	const bool setPartners=true);
	//@}
	+protected:
	+
	+ /** @name Clone Methods. */
	+ //@{
	+ /**
	+ * Make a simple clone of this object.
	+ * @return a pointer to the new object.
	+ */
	+ virtual IBPtr clone() const {return new_ptr(*this);}
	+
	+ /** Make a clone of this object, possibly modifying the cloned object
	+ * to make it sane.
	+ * @return a pointer to the new object.
	+ */
	+ virtual IBPtr fullclone() const {return new_ptr(*this);}
	+ //@}
	+

	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:

	/**
	+ * Access function for the initial conditions for the shower
	+ */
	+ //@{
	+ /**
	+ * Initial conditions for the shower of a final-final colour connection
	+ * - 0 is the symmetric choice
	+ * - 1 is maximal emmision from the coloured particle
	+ * - 2 is maximal emmision from the anticoloured particle
	+ * - 3 is randomly selected maximal emmision
	+ */
	+ unsigned int finalFinalConditions() const
	+ {return _finalFinalConditions;}
	+
	+ /**
	+ * Initial conditions for the shower of an initial-final decay colour connection
	+ * - 0 is the symmetric choice
	+ * - 1 is maximal emission from the decay product
	+ * - 2 is the smooth choice
	+ */
	+ unsigned int initialFinalDecayConditions() const
	+ {return _initialFinalDecayConditions;}
	+ //@}
	+
	+protected:
	+
	+ /**
	* Members to set the scales for different interactions
	*/
	//@{
	/**
	* Set initial scales for a QCD interaction
	*/
	virtual void setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	const bool setPartners=true);

	/**
	* Set initial scales for a QED interaction
	*/
	virtual void setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	const bool setPartners=true);
	//@}

	/**
	* Find the QCD partners
	* @param particle The particle to find the partners for
	* @param particles The full set of particles to search
	*/
	vector< pair<ShowerPartnerType, tShowerParticlePtr> >
	findQCDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles);

	/**
	* Find the QED partners
	* @param particle The particle to find the partners for
	* @param particles The full set of particles to search
	*/
	vector< pair<double, tShowerParticlePtr> >
	findQEDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles,
	const bool isDecayCase);

	+public:
	/**
	* Given a pair of particles, supposedly partners w.r.t. an interaction,
	* this method returns their initial evolution scales as a pair.
	* If something wrong happens, it returns the null (ZERO,ZERO) pair.
	* This method is used by the above setXXXInitialEvolutionScales
	* methods.
	* These methods must be overiden in inheriting classes
	*/
	//@{
	/**
	* General method to calculate the initial evolution scales
	*/
	- virtual pair<Energy,Energy> calculateInitialEvolutionScales(const ShowerPPair &,
	- const bool isDecayCase);
	+ pair<Energy,Energy> calculateInitialEvolutionScales(const ShowerPPair &,
	+ const bool isDecayCase);

	/**
	- * Calculate the initial evolution scales for two final-state particles
	+ * Calculate the initial evolution scales given momenta
	*/
	- virtual pair<Energy,Energy> calculateFinalFinalScales(const ShowerPPair &)=0;
	+ pair<Energy,Energy> calculateFinalFinalScales(const Lorentz5Momentum & p1,
	+ const Lorentz5Momentum & p2,
	+ bool colouredfirst);

	/**
	- * Calculate the initial evolution scales for two initial-state particles
	+ * Calculate the initial evolution scales given momenta
	*/
	- virtual pair<Energy,Energy> calculateInitialInitialScales(const ShowerPPair &)=0;
	+ pair<Energy,Energy> calculateInitialInitialScales(const Lorentz5Momentum& p1,
	+ const Lorentz5Momentum& p2);

	/**
	- * Calculate the initial evolution scales for one initial
	- * and one final-state particles
	+ * Calculate the initial evolution scales given momenta
	*/
	- virtual pair<Energy,Energy> calculateInitialFinalScales(const ShowerPPair &,
	- const bool isDecayCase)=0;
	+ pair<Energy,Energy> calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
	+ const bool isDecayCase);
	+
	+
	//@}

	private:

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

	private:

	/**
	* Method for choosing colour partner
	*/
	int partnerMethod_;

	/**
	* Choice for the QED radiation partner
	*/
	int QEDPartner_;

	/**
	* Choice of the scale
	*/
	int scaleChoice_;
	+
	+ /**
	+ * Flags controlling the initial conditions for the shower
	+ */
	+ //@{
	+ /**
	+ * Initial conditions for the shower with a final-final colour
	+ * connection
	+ */
	+ unsigned int _finalFinalConditions;
	+
	+ /**
	+ * Initial conditions for the shower with an initial-final decay colour
	+ * connection. This is done according to the top decay colour
	+ * connection calculation in JHEP12(2003)_045. The options act as follows:
	+ * 0: This is the default 'symmetric' choice which more or less divides
	+ * the phase space evenly between the parent and its charged child.
	+ * 1: This 'maximal' choice maximises the phase space available for
	+ * gluons emitted from the charged child.
	+ * 2: This (experimental) 'smooth' choice does not suffer from
	+ * a discontinuity at the boundary between the region populated by
	+ * emissions from the charged child and the region populated by emissions
	+ * from the parent. This does, however, mean that the phase space
	+ * available for emissions from the charged child is fairly minimal.
	+ */
	+ unsigned int _initialFinalDecayConditions;
	+
	+ /**
	+ * Initial conditions for the shower with an initial-initial colour
	+ * connection. This is done according to the colour connection
	+ * calculation in JHEP12(2003)_045. The options act as follows:
	+ * 0: This is the default 'symmetric' choice which more or less divides
	+ * the phase space evenly between the two incoming partons.
	+ * 1: This increases the phase space for emission from "parton b".
	+ * 2: This increases the phase space for emission from "parton c".
	+ */
	+ unsigned int _initialInitialConditions;
	+ //@}
	};

	}

	#endif /* HERWIG_PartnerFinder_H */
	diff --git a/Shower/QTilde/Default/QTildeFinder.cc b/Shower/QTilde/Default/QTildeFinder.cc
	deleted file mode 100644
	--- a/Shower/QTilde/Default/QTildeFinder.cc
	+++ /dev/null
	@@ -1,273 +0,0 @@
	-// -- C++ --
	-//
	-// QTildeFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	-// Copyright (C) 2002-2017 The Herwig Collaboration
	-//
	-// Herwig is licenced under version 3 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 QTildeFinder class.
	-//
	-
	-#include "QTildeFinder.h"
	-#include "ThePEG/Interface/ClassDocumentation.h"
	-#include "ThePEG/Persistency/PersistentOStream.h"
	-#include "ThePEG/Persistency/PersistentIStream.h"
	-#include "ThePEG/Interface/Switch.h"
	-#include "ThePEG/Repository/EventGenerator.h"
	-#include "ThePEG/EventRecord/Event.h"
	-#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
	-#include "ThePEG/Repository/UseRandom.h"
	-#include "ThePEG/Utilities/DescribeClass.h"
	-
	-using namespace Herwig;
	-
	-DescribeClass<QTildeFinder,Herwig::PartnerFinder>
	-describeQTildeFinder ("Herwig::QTildeFinder","HwShower.so");
	-
	-void QTildeFinder::persistentOutput(PersistentOStream & os) const {
	- os << _finalFinalConditions << _initialFinalDecayConditions
	- << _initialInitialConditions;
	-}
	-
	-void QTildeFinder::persistentInput(PersistentIStream & is, int) {
	- is >> _finalFinalConditions >> _initialFinalDecayConditions
	- >>_initialInitialConditions;
	-}
	-
	-void QTildeFinder::Init() {
	-
	- static ClassDocumentation<QTildeFinder> documentation
	- ("This class is responsible for finding the partners for each interaction types ",
	- "and within the evolution scale range specified by the ShowerVariables ",
	- "then to determine the initial evolution scales for each pair of partners.");
	-
	- static Switch<QTildeFinder,unsigned int> interfaceFinalFinalConditions
	- ("FinalFinalConditions",
	- "The initial conditions for the shower of a final-final colour connection",
	- &QTildeFinder::_finalFinalConditions, 0, false, false);
	- static SwitchOption interfaceFinalFinalConditionsSymmetric
	- (interfaceFinalFinalConditions,
	- "Symmetric",
	- "The symmetric choice",
	- 0);
	- static SwitchOption interfaceFinalFinalConditionsColoured
	- (interfaceFinalFinalConditions,
	- "Coloured",
	- "Maximal radiation from the coloured particle",
	- 1);
	- static SwitchOption interfaceFinalFinalConditionsAntiColoured
	- (interfaceFinalFinalConditions,
	- "AntiColoured",
	- "Maximal emission from the anticoloured particle",
	- 2);
	- static SwitchOption interfaceFinalFinalConditionsRandom
	- (interfaceFinalFinalConditions,
	- "Random",
	- "Randomly selected maximal emission from one of the particles",
	- 3);
	-
	- static Switch<QTildeFinder,unsigned int> interfaceInitialFinalDecayConditions
	- ("InitialFinalDecayConditions",
	- "The initial conditions for the shower of an initial-final"
	- " decay colour connection.",
	- &QTildeFinder::_initialFinalDecayConditions, 0, false, false);
	- static SwitchOption interfaceInitialFinalDecayConditionsSymmetric
	- (interfaceInitialFinalDecayConditions,
	- "Symmetric",
	- "The symmetric choice",
	- 0);
	- static SwitchOption interfaceInitialFinalDecayConditionsMaximal
	- (interfaceInitialFinalDecayConditions,
	- "Maximal",
	- "Maximal radiation from the decay product",
	- 1);
	- static SwitchOption interfaceInitialFinalDecayConditionsSmooth
	- (interfaceInitialFinalDecayConditions,
	- "Smooth",
	- "Smooth matching in the soft limit",
	- 2);
	-
	- static Switch<QTildeFinder,unsigned int> interfaceInitialInitialConditions
	- ("InitialInitialConditions",
	- "The initial conditions for the shower of an initial-initial"
	- " colour connection.",
	- &QTildeFinder::_initialInitialConditions, 0, false, false);
	- static SwitchOption interfaceInitialInitialConditionsSymmetric
	- (interfaceInitialInitialConditions,
	- "Symmetric",
	- "The symmetric choice",
	- 0);
	- static SwitchOption interfaceInitialInitialConditionsMaximiseB
	- (interfaceInitialInitialConditions,
	- "MaximiseB",
	- "Maximal radiation from parton b",
	- 1);
	- static SwitchOption interfaceInitialInitialConditionsMaximiseC
	- (interfaceInitialInitialConditions,
	- "MaximiseC",
	- "Maximal radiation from parton c",
	- 2);
	-}
	-
	-pair<Energy,Energy> QTildeFinder::
	-calculateInitialFinalScales(const ShowerPPair &ppair, const bool isDecayCase) {
	- return
	- calculateInitialFinalScales(ppair.first->momentum(),ppair.second->momentum(),isDecayCase);
	-}
	-
	-pair<Energy,Energy> QTildeFinder::
	-calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
	- const bool isDecayCase) {
	- if(!isDecayCase) {
	- // In this case from JHEP 12(2003)045 we find the conditions
	- // ktilde_b = (1+c) and ktilde_c = (1+2c)
	- // We also find that c = m_c^2/Q^2. The process is a+b->c where
	- // particle a is not colour connected (considered as a colour singlet).
	- // Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and
	- // q_c = sqrt(Q^2+2 m_c^2)
	- // We also assume that the first particle in the pair is the initial
	- // state particle and the second is the final state one c
	- Energy2 mc2 = sqr(pc.mass());
	- Energy2 Q2 = -(pb-pc).m2();
	- return pair<Energy,Energy>(sqrt(Q2+mc2), sqrt(Q2+2*mc2));
	- }
	- else {
	- // In this case from JHEP 12(2003)045 we find, for the decay
	- // process b->c+a(neutral), the condition
	- // (ktilde_b-1)(ktilde_c-c)=(1/4)sqr(1-a+c+lambda).
	- // We also assume that the first particle in the pair is the initial
	- // state particle (b) and the second is the final state one (c).
	- // - We find maximal phase space coverage through emissions from
	- // c if we set ktilde_c = 4.*(sqr(1.-sqrt(a))-c)
	- // - We find the most 'symmetric' way to populate the phase space
	- // occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda)
	- // - We find the most 'smooth' way to populate the phase space
	- // occurs for...
	- Energy2 mb2(sqr(pb.mass()));
	- double a=(pb-pc).m2()/mb2;
	- double c=sqr(pc.mass())/mb2;
	- double lambda = 1. + aa + cc - 2.a - 2.c - 2.ac;
	- lambda = sqrt(max(lambda,0.));
	- double PROD = 0.25*sqr(1. - a + c + lambda);
	- double ktilde_b, ktilde_c,cosi(0.);
	- switch(initialFinalDecayConditions()) {
	- case 0: // the 'symmetric' choice
	- ktilde_c = 0.5*(1-a+c+lambda) + c ;
	- ktilde_b = 1.+PROD/(ktilde_c-c) ;
	- break;
	- case 1: // the 'maximal' choice
	- ktilde_c = 4.0*(sqr(1.-sqrt(a))-c);
	- ktilde_b = 1.+PROD/(ktilde_c-c) ;
	- break;
	- case 2: // the 'smooth' choice
	- // c is a problem if very small here use 1GeV as minimum
	- c = max(c,1.*GeV2/mb2);
	- cosi = (sqr(1-sqrt(c))-a)/lambda;
	- ktilde_b = 2.0/(1.0-cosi);
	- ktilde_c = (1.0-a+c+lambda)(1.0+c-a-lambdacosi)/(2.0*(1.0+cosi));
	- break;
	- default:
	- throw Exception() << "Invalid option for decay shower's phase space"
	- << " QTildeFinder::calculateInitialFinalScales"
	- << Exception::abortnow;
	- }
	- return pair<Energy,Energy>(sqrt(mb2ktilde_b),sqrt(mb2ktilde_c));
	- }
	-}
	-
	-pair<Energy,Energy> QTildeFinder::
	-calculateInitialInitialScales(const ShowerPPair &ppair) {
	- return
	- calculateInitialInitialScales(ppair.first->momentum(),
	- ppair.second->momentum());
	-}
	-
	-pair<Energy,Energy> QTildeFinder::
	-calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2) {
	- // This case is quite simple. From JHEP 12(2003)045 we find the condition
	- // that ktilde_b = ktilde_c = 1. In this case we have the process
	- // b+c->a so we need merely boost to the CM frame of the two incoming
	- // particles and then qtilde is equal to the energy in that frame
	- Energy Q = sqrt((p1+p2).m2());
	- if(_initialInitialConditions==1) {
	- return pair<Energy,Energy>(sqrt(2.0)Q,sqrt(0.5)Q);
	- } else if(_initialInitialConditions==2) {
	- return pair<Energy,Energy>(sqrt(0.5)Q,sqrt(2.0)Q);
	- } else {
	- return pair<Energy,Energy>(Q,Q);
	- }
	-}
	-
	-pair<Energy,Energy> QTildeFinder::
	-calculateFinalFinalScales(const ShowerPPair & pp) {
	- bool colouredFirst =
	- pp.first->colourLine()&&
	- pp.first->colourLine()==pp.second->antiColourLine();
	- return calculateFinalFinalScales(pp.first->momentum(),pp.second->momentum(),
	- colouredFirst);
	-}
	-
	-pair<Energy,Energy> QTildeFinder::
	-calculateFinalFinalScales(Lorentz5Momentum p1, Lorentz5Momentum p2,
	- bool colouredFirst) {
	- static const double eps=1e-7;
	- // Using JHEP 12(2003)045 we find that we need ktilde = 1/2(1+b-c+lambda)
	- // ktilde = qtilde^2/Q^2 therefore qtilde = sqrt(ktilde*Q^2)
	- // find momenta in rest frame of system
	- // calculate quantities for the scales
	- Energy2 Q2 = (p1+p2).m2();
	- double b = p1.mass2()/Q2;
	- double c = p2.mass2()/Q2;
	- if(b<0.) {
	- if(b<-eps) {
	- throw Exception() << "Negative Mass squared b = " << b
	- << "in QTildeFinder::calculateFinalFinalScales()"
	- << Exception::eventerror;
	- }
	- b = 0.;
	- }
	- if(c<0.) {
	- if(c<-eps) {
	- throw Exception() << "Negative Mass squared c = " << c
	- << "in QTildeFinder::calculateFinalFinalScales()"
	- << Exception::eventerror;
	- }
	- c = 0.;
	- }
	- // KMH & PR - 16 May 2008 - swapped lambda calculation from
	- // double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)),
	- // which should be identical for p1 & p2 onshell in their COM
	- // but in the inverse construction for the Nason method, this
	- // was not the case, leading to misuse.
	- double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c))
	- (sqrt(b)-1.-sqrt(c))(sqrt(c)-1.-sqrt(b)));
	- // symmetric case
	- unsigned int iopt=finalFinalConditions();
	- Energy firstQ,secondQ;
	- if(iopt==0) {
	- firstQ = sqrt(0.5Q2(1.+b-c+lam));
	- secondQ = sqrt(0.5Q2(1.-b+c+lam));
	- }
	- // assymetric choice
	- else {
	- double kappab,kappac;
	- // calculate kappa with coloured line getting maximum
	- if((iopt==1&&colouredFirst)\|\| // first particle coloured+maximal for coloured
	- (iopt==2&&!colouredFirst)\|\| // first particle anticoloured+maximal for acoloured
	- (iopt==3&&UseRandom::rndbool(0.5))) { // random choice
	- kappab=4.(1.-2.sqrt(c)-b+c);
	- kappac=c+0.25*sqr(1.-b-c+lam)/(kappab-b);
	- }
	- else {
	- kappac=4.(1.-2.sqrt(b)-c+b);
	- kappab=b+0.25*sqr(1.-b-c+lam)/(kappac-c);
	- }
	- // calculate the scales
	- firstQ = sqrt(Q2*kappab);
	- secondQ = sqrt(Q2*kappac);
	- }
	- return pair<Energy,Energy>(firstQ, secondQ);
	-}
	diff --git a/Shower/QTilde/Default/QTildeFinder.h b/Shower/QTilde/Default/QTildeFinder.h
	deleted file mode 100644
	--- a/Shower/QTilde/Default/QTildeFinder.h
	+++ /dev/null
	@@ -1,206 +0,0 @@
	-// -- C++ --
	-//
	-// QTildeFinder.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	-// Copyright (C) 2002-2017 The Herwig Collaboration
	-//
	-// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	-// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	-//
	-#ifndef HERWIG_QTildeFinder_H
	-#define HERWIG_QTildeFinder_H
	-//
	-// This is the declaration of the QTildeFinder class.
	-//
	-
	-#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
	-#include "Herwig/Shower/QTilde/ShowerConfig.h"
	-#include "ThePEG/Interface/Interfaced.h"
	-
	-namespace Herwig {
	-using namespace ThePEG;
	-
	-/** \ingroup Shower
	- *
	- * The QTildeFinder class is responsible for finding the partners and
	- * setting the initial evolution scales for the shower evolution described
	- * in JHEP 0312:045,2003.
	- *
	- * @see \ref QTildeFinderInterfaces "The interfaces"
	- * defined for QTildeFinder.
	- */
	-class QTildeFinder: public PartnerFinder {
	-
	-public:
	-
	- /**
	- * The default constructor.
	- */
	- QTildeFinder() : _finalFinalConditions(0),
	- _initialFinalDecayConditions(0),
	- _initialInitialConditions(0) {}
	-
	- /** @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();
	-
	-public:
	-
	- /**
	- * Calculate the initial evolution scales given momenta
	- */
	- pair<Energy,Energy> calculateFinalFinalScales(Lorentz5Momentum p1, Lorentz5Momentum p2,
	- bool colouredfirst);
	-
	- /**
	- * Calculate the initial evolution scales given momenta
	- */
	- pair<Energy,Energy> calculateInitialInitialScales(const Lorentz5Momentum& p1,
	- const Lorentz5Momentum& p2);
	-
	- /**
	- * Calculate the initial evolution scales given momenta
	- */
	- pair<Energy,Energy> calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
	- const bool isDecayCase);
	-
	-protected:
	-
	- /**
	- * Given a pair of particles, supposedly partners w.r.t. an interaction,
	- * this method returns their initial evolution scales as a pair.
	- * If something wrong happens, it returns the null (ZERO,ZERO) pair.
	- * This method is used by the above setXXXInitialEvolutionScales
	- * methods.
	- */
	- //@{
	- /**
	- * Calculate the initial evolution scales for two final-state particles
	- */
	- virtual pair<Energy,Energy> calculateFinalFinalScales(const ShowerPPair &);
	-
	- /**
	- * Calculate the initial evolution scales for two initial-state particles
	- */
	- virtual pair<Energy,Energy> calculateInitialInitialScales(const ShowerPPair &);
	-
	- /**
	- * Calculate the initial evolution scales for one initial
	- * and one final-state particles
	- */
	- virtual pair<Energy,Energy> calculateInitialFinalScales(const ShowerPPair &,
	- const bool isDecayCase);
	- //@}
	-
	- /**
	- * Access function for the initial conditions for the shower
	- */
	- //@{
	- /**
	- * Initial conditions for the shower of a final-final colour connection
	- * - 0 is the symmetric choice
	- * - 1 is maximal emmision from the coloured particle
	- * - 2 is maximal emmision from the anticoloured particle
	- * - 3 is randomly selected maximal emmision
	- */
	- unsigned int finalFinalConditions() const
	- {return _finalFinalConditions;}
	-
	- /**
	- * Initial conditions for the shower of an initial-final decay colour connection
	- * - 0 is the symmetric choice
	- * - 1 is maximal emission from the decay product
	- * - 2 is the smooth choice
	- */
	- unsigned int initialFinalDecayConditions() const
	- {return _initialFinalDecayConditions;}
	- //@}
	-
	-protected:
	-
	- /** @name Clone Methods. */
	- //@{
	- /**
	- * Make a simple clone of this object.
	- * @return a pointer to the new object.
	- */
	- virtual IBPtr clone() const {return new_ptr(*this);}
	-
	- /** Make a clone of this object, possibly modifying the cloned object
	- * to make it sane.
	- * @return a pointer to the new object.
	- */
	- virtual IBPtr fullclone() const {return new_ptr(*this);}
	- //@}
	-
	-private:
	-
	- /**
	- * The assignment operator is private and must never be called.
	- * In fact, it should not even be implemented.
	- */
	- QTildeFinder & operator=(const QTildeFinder &);
	-
	-private:
	-
	- /**
	- * Flags controlling the initial conditions for the shower
	- */
	- //@{
	- /**
	- * Initial conditions for the shower with a final-final colour
	- * connection
	- */
	- unsigned int _finalFinalConditions;
	-
	- /**
	- * Initial conditions for the shower with an initial-final decay colour
	- * connection. This is done according to the top decay colour
	- * connection calculation in JHEP12(2003)_045. The options act as follows:
	- * 0: This is the default 'symmetric' choice which more or less divides
	- * the phase space evenly between the parent and its charged child.
	- * 1: This 'maximal' choice maximises the phase space available for
	- * gluons emitted from the charged child.
	- * 2: This (experimental) 'smooth' choice does not suffer from
	- * a discontinuity at the boundary between the region populated by
	- * emissions from the charged child and the region populated by emissions
	- * from the parent. This does, however, mean that the phase space
	- * available for emissions from the charged child is fairly minimal.
	- */
	- unsigned int _initialFinalDecayConditions;
	-
	- /**
	- * Initial conditions for the shower with an initial-initial colour
	- * connection. This is done according to the colour connection
	- * calculation in JHEP12(2003)_045. The options act as follows:
	- * 0: This is the default 'symmetric' choice which more or less divides
	- * the phase space evenly between the two incoming partons.
	- * 1: This increases the phase space for emission from "parton b".
	- * 2: This increases the phase space for emission from "parton c".
	- */
	- unsigned int _initialInitialConditions;
	- //@}
	-};
	-
	-}
	-
	-#endif /* HERWIG_QTildeFinder_H */
	diff --git a/Shower/QTilde/Makefile.am b/Shower/QTilde/Makefile.am
	--- a/Shower/QTilde/Makefile.am
	+++ b/Shower/QTilde/Makefile.am
	@@ -1,43 +1,42 @@
	SUBDIRS = Matching
	pkglib_LTLIBRARIES = HwShower.la
	HwShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 25:0:0
	HwShower_la_SOURCES = \
	Couplings/ShowerAlphaQCD.h Couplings/ShowerAlphaQCD.cc \
	Couplings/ShowerAlphaQED.h Couplings/ShowerAlphaQED.cc\
	QTildeShowerHandler.h QTildeShowerHandler.fh QTildeShowerHandler.cc \
	SplittingFunctions/HalfHalfOneSplitFn.h SplittingFunctions/HalfHalfOneSplitFn.cc\
	SplittingFunctions/OneOneOneSplitFn.h SplittingFunctions/OneOneOneSplitFn.cc\
	SplittingFunctions/OneOneOneMassiveSplitFn.h SplittingFunctions/OneOneOneMassiveSplitFn.cc\
	SplittingFunctions/ZeroZeroOneSplitFn.h SplittingFunctions/ZeroZeroOneSplitFn.cc\
	SplittingFunctions/OneHalfHalfSplitFn.h SplittingFunctions/OneHalfHalfSplitFn.cc\
	SplittingFunctions/HalfOneHalfSplitFn.h SplittingFunctions/HalfOneHalfSplitFn.cc\
	SplittingFunctions/CMWOneOneOneSplitFn.h SplittingFunctions/CMWOneOneOneSplitFn.cc\
	SplittingFunctions/CMWHalfHalfOneSplitFn.h SplittingFunctions/CMWHalfHalfOneSplitFn.cc\
	Default/QTildeSudakov.cc Default/QTildeSudakov.h\
	Default/Decay_QTildeShowerKinematics1to2.cc \
	Default/Decay_QTildeShowerKinematics1to2.h \
	Default/IS_QTildeShowerKinematics1to2.cc Default/IS_QTildeShowerKinematics1to2.h \
	Default/FS_QTildeShowerKinematics1to2.cc Default/FS_QTildeShowerKinematics1to2.h \
	-Default/QTildeFinder.cc Default/QTildeFinder.h\
	Default/QTildeReconstructor.cc Default/QTildeReconstructor.h Default/QTildeReconstructor.tcc \
	Base/KinematicsReconstructor.cc \
	Base/KinematicsReconstructor.h \
	Base/KinematicsReconstructor.fh \
	Base/HardTree.cc Base/HardTree.h Base/HardTree.fh \
	Base/HardBranching.h Base/HardBranching.fh Base/HardBranching.cc\
	Base/PartnerFinder.h Base/PartnerFinder.fh Base/PartnerFinder.cc \
	Base/ShowerVeto.h Base/ShowerVeto.fh Base/ShowerVeto.cc \
	Base/FullShowerVeto.h Base/FullShowerVeto.fh Base/FullShowerVeto.cc \
	SplittingFunctions/SplittingGenerator.cc SplittingFunctions/SplittingGenerator.h\
	SplittingFunctions/SplittingGenerator.fh \
	Base/ShowerTree.h Base/ShowerTree.fh Base/ShowerTree.cc \
	ShowerConfig.h ShowerConfig.cc \
	Base/Branching.h \
	Base/ShowerParticle.cc Base/ShowerParticle.fh Base/ShowerParticle.h \
	Base/ShowerKinematics.fh Base/ShowerKinematics.h Base/ShowerKinematics.cc \
	Base/ShowerBasis.fh Base/ShowerBasis.h Base/ShowerBasis.cc \
	Base/ShowerProgenitor.fh Base/ShowerProgenitor.h \
	Base/SudakovFormFactor.cc Base/SudakovFormFactor.h Base/SudakovFormFactor.fh \
	SplittingFunctions/SplittingFunction.h SplittingFunctions/SplittingFunction.fh \
	SplittingFunctions/SplittingFunction.cc \
	Base/ShowerVertex.cc Base/ShowerVertex.fh Base/ShowerVertex.h
	diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in
	--- a/src/defaults/Shower.in
	+++ b/src/defaults/Shower.in
	@@ -1,310 +1,310 @@
	# -- ThePEG-repository --

	############################################################
	# Setup of default parton shower
	#
	# Useful switches for users are marked near the top of
	# this file.
	#
	# Don't edit this file directly, but reset the switches
	# in your own input files!
	############################################################
	library HwMPI.so
	library HwShower.so
	library HwMatching.so

	mkdir /Herwig/Shower
	cd /Herwig/Shower

	create Herwig::QTildeShowerHandler ShowerHandler
	newdef ShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
	newdef ShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
	# use LO PDFs for Shower, can be changed later
	newdef ShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
	newdef ShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
	newdef ShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
	newdef ShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF

	#####################################
	# initial setup, don't change these!
	#####################################
	create Herwig::SplittingGenerator SplittingGenerator
	create Herwig::ShowerAlphaQCD AlphaQCD
	create Herwig::ShowerAlphaQED AlphaQED
	-create Herwig::QTildeFinder PartnerFinder
	+create Herwig::PartnerFinder PartnerFinder
	newdef PartnerFinder:PartnerMethod 1
	newdef PartnerFinder:ScaleChoice 1
	create Herwig::QTildeReconstructor KinematicsReconstructor
	newdef KinematicsReconstructor:ReconstructionOption Colour3
	newdef KinematicsReconstructor:InitialStateReconOption SofterFraction
	newdef KinematicsReconstructor:InitialInitialBoostOption LongTransBoost

	newdef /Herwig/Partons/RemnantDecayer:AlphaS AlphaQCD
	newdef /Herwig/Partons/RemnantDecayer:AlphaEM AlphaQED

	newdef ShowerHandler:PartnerFinder PartnerFinder
	newdef ShowerHandler:KinematicsReconstructor KinematicsReconstructor
	newdef ShowerHandler:SplittingGenerator SplittingGenerator
	newdef ShowerHandler:SpinCorrelations Yes
	newdef ShowerHandler:SoftCorrelations Singular

	##################################################################
	# Intrinsic pT
	#
	# Recommended:
	# 1.9 GeV for Tevatron W/Z production.
	# 2.1 GeV for LHC W/Z production at 10 TeV
	# 2.2 GeV for LHC W/Z production at 14 TeV
	#
	# Set all parameters to 0 to disable
	##################################################################
	newdef ShowerHandler:IntrinsicPtGaussian 1.3*GeV
	newdef ShowerHandler:IntrinsicPtBeta 0
	newdef ShowerHandler:IntrinsicPtGamma 0*GeV
	newdef ShowerHandler:IntrinsicPtIptmax 0*GeV
	#############################################################
	# Set up truncated shower handler.
	#############################################################

	create Herwig::PowhegShowerHandler PowhegShowerHandler
	set PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
	set PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
	newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
	newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
	newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
	newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
	newdef PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
	newdef PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
	newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
	newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
	newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
	newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
	newdef PowhegShowerHandler:PartnerFinder PartnerFinder
	newdef PowhegShowerHandler:KinematicsReconstructor KinematicsReconstructor
	newdef PowhegShowerHandler:SplittingGenerator SplittingGenerator
	newdef PowhegShowerHandler:Interactions QCDandQED
	newdef PowhegShowerHandler:SpinCorrelations Yes
	newdef PowhegShowerHandler:SoftCorrelations Singular
	newdef PowhegShowerHandler:IntrinsicPtGaussian 1.3*GeV
	newdef PowhegShowerHandler:IntrinsicPtBeta 0
	newdef PowhegShowerHandler:IntrinsicPtGamma 0*GeV
	newdef PowhegShowerHandler:IntrinsicPtIptmax 0*GeV
	newdef PowhegShowerHandler:ReconstructionOption OffShell5

	#############################################################
	# End of interesting user servicable section.
	#
	# Anything that follows below should only be touched if you
	# know what you're doing.
	#
	# Really.
	#############################################################
	#
	# a few default values
	newdef ShowerHandler:MECorrMode 1
	newdef ShowerHandler:ReconstructionOption OffShell5
	newdef AlphaQCD:ScaleFactor 1.0
	newdef AlphaQCD:NPAlphaS 2
	newdef AlphaQCD:Qmin 0.935
	newdef AlphaQCD:NumberOfLoops 2
	newdef AlphaQCD:InputOption 1
	newdef AlphaQCD:AlphaMZ 0.126234
	#
	#
	# Lets set up all the splittings
	create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn
	set QtoQGammaSplitFn:InteractionType QED
	set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral
	set QtoQGammaSplitFn:AngularOrdered Yes

	create Herwig::HalfHalfOneSplitFn QtoQGSplitFn
	newdef QtoQGSplitFn:InteractionType QCD
	newdef QtoQGSplitFn:ColourStructure TripletTripletOctet
	set QtoQGSplitFn:AngularOrdered Yes

	create Herwig::OneOneOneSplitFn GtoGGSplitFn
	newdef GtoGGSplitFn:InteractionType QCD
	newdef GtoGGSplitFn:ColourStructure OctetOctetOctet
	set GtoGGSplitFn:AngularOrdered Yes

	create Herwig::OneOneOneMassiveSplitFn WtoWGammaSplitFn
	newdef WtoWGammaSplitFn:InteractionType QED
	newdef WtoWGammaSplitFn:ColourStructure ChargedChargedNeutral
	set WtoWGammaSplitFn:AngularOrdered Yes

	create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn
	newdef GtoQQbarSplitFn:InteractionType QCD
	newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet
	set GtoQQbarSplitFn:AngularOrdered Yes

	create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn
	newdef GammatoQQbarSplitFn:InteractionType QED
	newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged
	set GammatoQQbarSplitFn:AngularOrdered Yes

	create Herwig::HalfOneHalfSplitFn QtoGQSplitFn
	newdef QtoGQSplitFn:InteractionType QCD
	newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet
	set QtoGQSplitFn:AngularOrdered Yes

	create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn
	newdef QtoGammaQSplitFn:InteractionType QED
	newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged
	set QtoGammaQSplitFn:AngularOrdered Yes
	#
	# Now the Sudakovs
	create Herwig::QTildeSudakov SudakovCommon
	newdef SudakovCommon:Alpha AlphaQCD
	newdef SudakovCommon:cutoffKinScale 0.0*GeV
	newdef SudakovCommon:PDFmax 1.0
	newdef SudakovCommon:CutOffOption pT
	newdef SudakovCommon:pTmin 1.222798*GeV

	cp SudakovCommon QtoQGSudakov
	newdef QtoQGSudakov:SplittingFunction QtoQGSplitFn
	newdef QtoQGSudakov:PDFmax 1.9

	cp SudakovCommon QtoQGammaSudakov
	set QtoQGammaSudakov:SplittingFunction QtoQGammaSplitFn
	set QtoQGammaSudakov:Alpha AlphaQED
	set QtoQGammaSudakov:PDFmax 1.9

	cp QtoQGammaSudakov LtoLGammaSudakov
	# Technical parameter to stop evolution.
	set LtoLGammaSudakov:pTmin 0.000001

	cp SudakovCommon GtoGGSudakov
	newdef GtoGGSudakov:SplittingFunction GtoGGSplitFn
	newdef GtoGGSudakov:PDFmax 2.0

	cp SudakovCommon WtoWGammaSudakov
	newdef WtoWGammaSudakov:SplittingFunction WtoWGammaSplitFn
	set WtoWGammaSudakov:Alpha AlphaQED

	cp SudakovCommon GtoQQbarSudakov
	newdef GtoQQbarSudakov:SplittingFunction GtoQQbarSplitFn
	newdef GtoQQbarSudakov:PDFmax 120.0

	cp SudakovCommon GammatoQQbarSudakov
	newdef GammatoQQbarSudakov:SplittingFunction GammatoQQbarSplitFn
	set GammatoQQbarSudakov:Alpha AlphaQED
	newdef GammatoQQbarSudakov:PDFmax 120.0

	cp SudakovCommon GtobbbarSudakov
	newdef GtobbbarSudakov:SplittingFunction GtoQQbarSplitFn
	newdef GtobbbarSudakov:PDFmax 40000.0

	cp SudakovCommon GtoccbarSudakov
	newdef GtoccbarSudakov:SplittingFunction GtoQQbarSplitFn
	newdef GtoccbarSudakov:PDFmax 2000.0

	cp SudakovCommon QtoGQSudakov
	newdef QtoGQSudakov:SplittingFunction QtoGQSplitFn

	cp SudakovCommon QtoGammaQSudakov
	newdef QtoGammaQSudakov:SplittingFunction QtoGammaQSplitFn
	set QtoGammaQSudakov:Alpha AlphaQED

	cp SudakovCommon utoGuSudakov
	newdef utoGuSudakov:SplittingFunction QtoGQSplitFn
	newdef utoGuSudakov:PDFFactor OverOneMinusZ
	newdef utoGuSudakov:PDFmax 5.0

	cp SudakovCommon dtoGdSudakov
	newdef dtoGdSudakov:SplittingFunction QtoGQSplitFn
	newdef dtoGdSudakov:PDFFactor OverOneMinusZ


	#
	# Now add the final splittings
	#
	do SplittingGenerator:AddFinalSplitting u->u,g; QtoQGSudakov
	do SplittingGenerator:AddFinalSplitting d->d,g; QtoQGSudakov
	do SplittingGenerator:AddFinalSplitting s->s,g; QtoQGSudakov
	do SplittingGenerator:AddFinalSplitting c->c,g; QtoQGSudakov
	do SplittingGenerator:AddFinalSplitting b->b,g; QtoQGSudakov
	do SplittingGenerator:AddFinalSplitting t->t,g; QtoQGSudakov
	#
	do SplittingGenerator:AddFinalSplitting g->g,g; GtoGGSudakov
	#
	do SplittingGenerator:AddFinalSplitting g->u,ubar; GtoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting g->d,dbar; GtoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting g->s,sbar; GtoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting g->c,cbar; GtoccbarSudakov
	do SplittingGenerator:AddFinalSplitting g->b,bbar; GtobbbarSudakov
	do SplittingGenerator:AddFinalSplitting g->t,tbar; GtoQQbarSudakov
	#
	do SplittingGenerator:AddFinalSplitting gamma->u,ubar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->d,dbar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->s,sbar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->c,cbar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->b,bbar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->t,tbar; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->e-,e+; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->mu-,mu+; GammatoQQbarSudakov
	do SplittingGenerator:AddFinalSplitting gamma->tau-,tau+; GammatoQQbarSudakov
	#
	do SplittingGenerator:AddFinalSplitting u->u,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddFinalSplitting d->d,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddFinalSplitting s->s,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddFinalSplitting c->c,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddFinalSplitting b->b,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddFinalSplitting t->t,gamma; QtoQGammaSudakov

	do SplittingGenerator:AddFinalSplitting e-->e-,gamma; LtoLGammaSudakov
	do SplittingGenerator:AddFinalSplitting mu-->mu-,gamma; LtoLGammaSudakov
	do SplittingGenerator:AddFinalSplitting tau-->tau-,gamma; LtoLGammaSudakov

	do SplittingGenerator:AddFinalSplitting W+->W+,gamma; WtoWGammaSudakov
	#
	# Now lets add the initial splittings. Remember the form a->b,c; means
	# that the current particle b is given and we backward branch to new
	# particle a which is initial state and new particle c which is final state
	#
	do SplittingGenerator:AddInitialSplitting u->u,g; QtoQGSudakov
	do SplittingGenerator:AddInitialSplitting d->d,g; QtoQGSudakov
	do SplittingGenerator:AddInitialSplitting s->s,g; QtoQGSudakov
	do SplittingGenerator:AddInitialSplitting c->c,g; QtoQGSudakov
	do SplittingGenerator:AddInitialSplitting b->b,g; QtoQGSudakov
	do SplittingGenerator:AddInitialSplitting u->u,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddInitialSplitting d->d,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddInitialSplitting s->s,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddInitialSplitting c->c,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddInitialSplitting b->b,gamma; QtoQGammaSudakov
	do SplittingGenerator:AddInitialSplitting t->t,gamma; QtoQGammaSudakov

	do SplittingGenerator:AddInitialSplitting g->g,g; GtoGGSudakov
	#
	do SplittingGenerator:AddInitialSplitting g->d,dbar; GtoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting g->u,ubar; GtoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting g->s,sbar; GtoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting g->c,cbar; GtoccbarSudakov
	do SplittingGenerator:AddInitialSplitting g->b,bbar; GtobbbarSudakov
	#
	do SplittingGenerator:AddInitialSplitting gamma->d,dbar; GammatoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting gamma->u,ubar; GammatoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting gamma->s,sbar; GammatoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting gamma->c,cbar; GammatoQQbarSudakov
	do SplittingGenerator:AddInitialSplitting gamma->b,bbar; GammatoQQbarSudakov
	#
	do SplittingGenerator:AddInitialSplitting d->g,d; dtoGdSudakov
	do SplittingGenerator:AddInitialSplitting u->g,u; utoGuSudakov
	do SplittingGenerator:AddInitialSplitting s->g,s; QtoGQSudakov
	do SplittingGenerator:AddInitialSplitting c->g,c; QtoGQSudakov
	do SplittingGenerator:AddInitialSplitting b->g,b; QtoGQSudakov
	do SplittingGenerator:AddInitialSplitting dbar->g,dbar; dtoGdSudakov
	do SplittingGenerator:AddInitialSplitting ubar->g,ubar; utoGuSudakov
	do SplittingGenerator:AddInitialSplitting sbar->g,sbar; QtoGQSudakov
	do SplittingGenerator:AddInitialSplitting cbar->g,cbar; QtoGQSudakov
	do SplittingGenerator:AddInitialSplitting bbar->g,bbar; QtoGQSudakov
	#
	do SplittingGenerator:AddInitialSplitting d->gamma,d; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting u->gamma,u; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting s->gamma,s; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting c->gamma,c; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting b->gamma,b; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting dbar->gamma,dbar; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting ubar->gamma,ubar; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting sbar->gamma,sbar; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting cbar->gamma,cbar; QtoGammaQSudakov
	do SplittingGenerator:AddInitialSplitting bbar->gamma,bbar; QtoGammaQSudakov

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 5:27 PM (12 h, 6 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023638
Default Alt Text: (99 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions