No OneTemporary
Actions

Size

282 KB

Subscribers

None

View Options

This file is larger than 256 KB, so syntax highlighting was skipped.

	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,538 +1,538 @@
	// -- 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"
	#include "Herwig/Shower/QTilde/Kinematics/KinematicHelpers.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 {

	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()]).first;
	+ bornCXComb()->meMomenta()[dipole()->bornSpectator()],0).first; //0 is the key to symmetric choice for final final scale
	}

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


	Energy2 pt2 = ZERO;

	const vector<Energy> & masses = theQTildeSudakov->virtualMasses({{
	bornCXComb()->mePartonData()[dipole()->bornEmitter() ],
	realCXComb()->mePartonData()[dipole()->realEmitter() ],
	realCXComb()->mePartonData()[dipole()->realEmission()]
	}});

	const Energy2 m22 = sqr(masses[2]);

	if ( dipole()->bornEmitter() > 1 ) {
	const Energy2 m02 = sqr(masses[0]);
	const Energy2 m12 = sqr(masses[1]);

	pt2 = QTildeKinematics::pT2_FSR(sqr(qtilde),z,m02,m12,m22);
	}
	else {
	pt2 = QTildeKinematics::pT2_ISR(sqr(qtilde),z,m22);
	}

	if ( pt2 < max(theQTildeSudakov->pT2min(),sqr(safeCut()) ))
	return false;

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

	}

	bool QTildeMatching::isAboveCutoff() const {
	Energy qtilde = dipole()->showerScale();
	assert(!dipole()->showerParameters().empty());
	double z = dipole()->showerParameters()[0];


	const vector<Energy> & masses = theQTildeSudakov->virtualMasses({{
	bornCXComb()->mePartonData()[dipole()->bornEmitter() ],
	realCXComb()->mePartonData()[dipole()->realEmitter() ],
	realCXComb()->mePartonData()[dipole()->realEmission()]
	}});

	const Energy2 m22 = sqr(masses[2]);

	if ( dipole()->bornEmitter() > 1 ) {
	const Energy2 m02 = sqr(masses[0]);
	const Energy2 m12 = sqr(masses[1]);

	const Energy2 pt2 = QTildeKinematics::pT2_FSR(sqr(qtilde),z,m02,m12,m22);

	return pt2 >= max(theQTildeSudakov->pT2min(),sqr(safeCut()));
	}
	else {
	const Energy2 pt2 = QTildeKinematics::pT2_ISR(sqr(qtilde),z,m22);
	return pt2 >= 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;

	const vector<Energy> & masses = theQTildeSudakov->virtualMasses({{
	bornCXComb()->mePartonData()[dipole()->bornEmitter() ],
	realCXComb()->mePartonData()[dipole()->realEmitter() ],
	realCXComb()->mePartonData()[dipole()->realEmission()]
	}});

	const Energy2 m22 = sqr(masses[2]);
	if ( dipole()->bornEmitter() > 1 ) {
	const Energy2 m02 = sqr(masses[0]);
	const Energy2 m12 = sqr(masses[1]);

	pt2 = QTildeKinematics::pT2_FSR(sqr(qtilde),z,m02,m12,m22);
	}
	else {
	pt2 = QTildeKinematics::pT2_ISR(sqr(qtilde),z,m22);
	}
	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,PartnerFinder> interfaceQTildeFinder
	("QTildeFinder",
	"Set the partner finder to calculate hard scales.",
	&QTildeMatching::theQTildeFinder, false, false, true, false, false);
	interfaceQTildeFinder.rank(-1);

	static Reference<QTildeMatching,SudakovFormFactor> 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/Shower/QTilde/Base/PartnerFinder.cc b/Shower/QTilde/Base/PartnerFinder.cc
	--- a/Shower/QTilde/Base/PartnerFinder.cc
	+++ b/Shower/QTilde/Base/PartnerFinder.cc
	@@ -1,694 +1,733 @@
	// -- 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;

	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 size_t
	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 size_t
	ACLSIZE(const tShowerParticlePtr a) {
	return a->colourInfo()->antiColourLines().size();
	}
	}

	void PartnerFinder::persistentOutput(PersistentOStream & os) const {
	os << partnerMethod_ << QEDPartner_ << scaleChoice_;
	}

	void PartnerFinder::persistentInput(PersistentIStream & is, int) {
	is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_;
	}

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

	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 if(type==ShowerInteraction::EW) {
	setInitialEWEvolutionScales(particles,isDecayCase,false);
	}
	else if(type==ShowerInteraction::QEDQCD) {
	setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
	setInitialQEDEvolutionScales(particles,isDecayCase,false);
	}
	else if(type==ShowerInteraction::ALL) {
	setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
	setInitialQEDEvolutionScales(particles,isDecayCase,false);
	setInitialEWEvolutionScales(particles,isDecayCase,false);
	}
	else
	assert(false);
	// \todo EW scales here
	// 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) {
	// 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!

	for ( const auto & sp : particles ) {
	// Skip colourless particles
	if(!sp->data().coloured()) continue;
	// find the partners
	auto partners = findQCDPartners(sp,particles);
	// must have a partner
	if(partners.empty()) {
	throw Exception() << "`Failed to make colour connections in "
	<< "PartnerFinder::setQCDInitialEvolutionScales"
	<< *sp
	<< Exception::eventerror;
	}
	// Calculate the evolution scales for all possible pairs of of particles
	vector<pair<Energy,Energy> > scales;
	int position = -1;
	for(size_t ix=0; ix< partners.size(); ++ix) {
	scales.push_back(
	- calculateInitialEvolutionScales(
	+ calculateInitialEvolutionScales( //QCD shower
	ShowerPPair(sp, partners[ix].second),
	- isDecayCase
	- )
	- );
	+ isDecayCase, 0)); //choices: 0->symmetric 1->colored 'maximal' 2->anti-coloured maximal
	if (!setPartners && partners[ix].second) position = ix;
	}
	assert(setPartners \|\| position >= 0);

	// set partners if required
	if (setPartners) {
	// 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.
	// 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 (sp->perturbative() == 1 &&
	sp->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 = sp->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
	sp->partner(partners[position].second);
	}


	// primary partner set, set the others and do the scale
	for(size_t ix=0; ix<partners.size(); ++ix) {
	sp->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) {
	sp->scales().QCD_c =
	sp->scales().QCD_c_noAO =
	(scaleChoice_==0 ? scale : scales[ix].first);
	}
	else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
	sp->scales().QCD_ac =
	sp->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(const auto & sp : particles) {
	// not charged or photon continue
	if(!sp->dataPtr()->charged()) continue;
	// find the potential partners
	vector<pair<double,tShowerParticlePtr> > partners = findQEDPartners(sp,particles,isDecayCase);
	if(partners.empty()) {
	throw Exception() << "Failed to find partner in "
	<< "PartnerFinder::setQEDInitialEvolutionScales"
	<< *sp << 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&&sp->partner()) {
	for(unsigned int ix=0;ix<partners.size();++ix) {
	if(sp->partner()==partners[ix].second) {
	position = ix;
	break;
	}
	}
	}
	// set the partner
	if(setPartners\|\|!sp->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\|\|!sp->partner())) {
	sp->partner(partners[position].second);
	}
	}
	// must have a partner
	if(position<0) throw Exception() << "Failed to find partner in "
	<< "PartnerFinder::setQEDInitialEvolutionScales"
	<< *sp << 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) {
	+ for(unsigned int ix=0;ix< partners.size();++ix) { //QED shower
	scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp,partners[ix].second),
	- isDecayCase));
	+ isDecayCase, 0)); //choices: 0->symmetric 1->colored 'maximal' 2->anti-coloured maximal
	}
	// store all the possible partners
	for(unsigned int ix=0;ix<partners.size();++ix) {
	sp->addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
	partners[ix].first,
	ShowerPartnerType::QED,
	scales[ix].first));
	}
	// set scales
	sp->scales().QED = scales[position].first;
	sp->scales().QED_noAO = scales[position].first;
	}
	}

	pair<Energy,Energy> PartnerFinder::
	calculateInitialEvolutionScales(const ShowerPPair &particlePair,
	- const bool isDecayCase) {
	+ const bool isDecayCase, int key) {
	bool FS1=FS(particlePair.first),FS2= FS(particlePair.second);
	-
	- if(FS1 && FS2)
	- return calculateFinalFinalScales(particlePair.first->momentum(),
	- particlePair.second->momentum());
	+ if(FS1 && FS2){
	+ return calculateFinalFinalScales(particlePair.first->momentum(),particlePair.second->momentum(), key);
	+ }
	else if(FS1 && !FS2) {
	pair<Energy,Energy> rval = calculateInitialFinalScales(particlePair.second->momentum(),
	particlePair.first->momentum(),
	isDecayCase);
	return { rval.second, rval.first };
	}
	else if(!FS1 &&FS2)
	return calculateInitialFinalScales(particlePair.first->momentum(),particlePair.second->momentum(),isDecayCase);
	else
	return calculateInitialInitialScales(particlePair.first->momentum(),particlePair.second->momentum());
	}

	vector< pair<ShowerPartnerType, tShowerParticlePtr> >
	PartnerFinder::findQCDPartners(tShowerParticlePtr particle,
	const ShowerParticleVector &particles) {
	vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners;
	for(const auto & sp : particles) {
	if(!sp->data().coloured() \|\| particle==sp) continue;
	// one initial-state and one final-state particle
	if(FS(particle) != FS(sp)) {
	// loop over all the colours of both particles
	for(size_t ix=0; ix<CLSIZE(particle); ++ix) {
	for(size_t jx=0; jx<CLSIZE(sp); ++jx) {
	if((CL(particle,ix) && CL(particle,ix)==CL(sp,jx))) {
	partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
	}
	}
	}
	//loop over all the anti-colours of both particles
	for(size_t ix=0; ix<ACLSIZE(particle); ++ix) {
	for(size_t jx=0; jx<ACLSIZE(sp); ++jx) {
	if((ACL(particle,ix) && ACL(particle,ix)==ACL(sp,jx))) {
	partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
	}
	}
	}
	}
	// 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(size_t ix=0; ix<CLSIZE(particle); ++ix){
	for(size_t jx=0; jx<ACLSIZE(sp); ++jx){
	if(CL(particle,ix) && CL(particle,ix)==ACL(sp,jx)) {
	partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
	}
	}
	}
	//loop over the anti-colours of the first particle and the colours of the other
	for(size_t ix=0; ix<ACLSIZE(particle); ++ix){
	for(size_t jx=0; jx<CLSIZE(sp); jx++){
	if(ACL(particle,ix) && ACL(particle,ix)==CL(sp,jx)) {
	partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
	}
	}
	}
	}
	}
	// 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(const auto & sp : particles) {
	if(( FS(sp) && ( CL(sp) == cpair.first \|\| CL(sp) == cpair.second))\|\|
	(!FS(sp) && (ACL(sp) == cpair.first \|\| ACL(sp) == cpair.second ))) {
	partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
	}
	}
	}
	else if(col&&col->sinkNeighbours().first) {
	tColinePair cpair = col->sinkNeighbours();
	for(const auto & sp : particles) {
	if(( FS(sp) && (ACL(sp) == cpair.first \|\| ACL(sp) == cpair.second))\|\|
	(!FS(sp) && ( CL(sp) == cpair.first \|\| CL(sp) == cpair.second))) {
	partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
	}
	}
	}
	col = ACL(particle);
	if(FS(particle)&&col&&col->sinkNeighbours().first) {
	tColinePair cpair = col->sinkNeighbours();
	for(const auto & sp : particles) {
	if(( FS(sp) && (ACL(sp) == cpair.first \|\| ACL(sp) == cpair.second))\|\|
	(!FS(sp) && ( CL(sp) == cpair.first \|\| CL(sp) == cpair.second ))) {
	partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
	}
	}
	}
	else if(col&&col->sourceNeighbours().first) {
	tColinePair cpair = col->sourceNeighbours();
	for(const auto & sp : particles) {
	if(( FS(sp) && ( CL(sp) == cpair.first \|\| CL(sp) == cpair.second))\|\|
	(!FS(sp) && (ACL(sp) == cpair.first \|\|ACL(sp) == cpair.second))) {
	partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
	}
	}
	}
	}
	// return the partners
	return partners;
	}

	vector< pair<double, tShowerParticlePtr> >
	PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
	const ShowerParticleVector & particles,
	const bool isDecayCase) {
	vector< pair<double, tShowerParticlePtr> > partners;
	const double pcharge =
	particle->id()==ParticleID::gamma ? 1 : double(particle->data().iCharge());
	vector< pair<double, tShowerParticlePtr> > photons;
	for (const auto & sp : particles) {
	if (particle == sp) continue;
	if (sp->id()==ParticleID::gamma) photons.push_back(make_pair(1.,sp));
	if (!sp->data().charged() ) continue;
	double charge = pcharge*double((sp)->data().iCharge());
	if ( FS(particle) != FS(sp) ) charge *=-1.;
	if ( QEDPartner_ != 0 && !isDecayCase ) {
	// only include II and FF as requested
	if ( QEDPartner_ == 1 && FS(particle) != FS(sp) )
	continue;
	// only include IF is requested
	else if (QEDPartner_ == 2 && FS(particle) == FS(sp) )
	continue;
	}
	if (particle->id()==ParticleID::gamma) charge = -abs(charge);
	// only keep positive dipoles
	if (charge<0.) partners.push_back({ -charge, sp });
	}
	if (particle->id()==ParticleID::gamma && partners.empty())
	return photons;
	return partners;
	}

	void PartnerFinder::setInitialEWEvolutionScales(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) {
	// if not weakly interacting continue
	if( !weaklyInteracting( (**cit).dataPtr()))
	continue;
	// find the potential partners
	vector<pair<double,tShowerParticlePtr> > partners = findEWPartners(*cit,particles,isDecayCase);
	if(partners.empty()) {
	throw Exception() << "Failed to find partner in "
	<< "PartnerFinder::setEWInitialEvolutionScales"
	<< (**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::setEWInitialEvolutionScales"
	<< (**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) {
	+ for(unsigned int ix=0;ix< partners.size();++ix) { //EW shower
	scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
	- isDecayCase));
	+ isDecayCase, 2)); //choices: 0->symmetric 1->colored 'maximal' 2->anti-coloured maximal
	}
	// 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::EW,
	scales[ix].first));
	}
	// set scales
	- (**cit).scales().EW = scales[position].first;
	+ (**cit).scales().EW_Z = scales[position].first;
	+ (**cit).scales().EW_W = scales[position].first;
	(**cit).scales().EW_noAO = scales[position].first;
	}
	}

	vector< pair<double, tShowerParticlePtr> >
	PartnerFinder::findEWPartners(tShowerParticlePtr particle,
	const ShowerParticleVector &particles,
	const bool isDecayCase ) {
	vector< pair<double, tShowerParticlePtr> > partners;
	ShowerParticleVector::const_iterator cjt;
	for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
	if(!weaklyInteracting((*cjt)->dataPtr()) \|\|
	particle == *cjt) continue;
	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;
	}
	double charge = 1.;
	// only keep positive dipoles
	if(charge>0.) partners.push_back(make_pair(charge,*cjt));
	}
	return partners;
	}

	pair<Energy,Energy>
	PartnerFinder::calculateFinalFinalScales(
	const Lorentz5Momentum & p1,
	- const Lorentz5Momentum & p2)
	+ const Lorentz5Momentum & p2, int key)
	{
	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.
	const 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
	- return { sqrt(0.5Q2(1.+b-c+lam)), sqrt(0.5Q2(1.-b+c+lam)) };
	+
	+ Energy firstQ, secondQ;
	+ double kappab(0.), kappac(0.);
	+ //key = 0; // symmetric case pre-selection
	+ switch(key)
	+ {
	+ case 0: // symmetric case
	+ firstQ = sqrt(0.5Q2(1.+b-c+lam));
	+ secondQ = sqrt(0.5Q2(1.-b+c+lam));
	+ break;
	+ case 1: // assymetric choices - colored
	+ kappac=4.(1.-2.sqrt(b)-c+b);
	+ kappab=b+0.25*sqr(1.-b-c+lam)/(kappac-c);
	+ firstQ = sqrt(Q2*kappab);
	+ secondQ = sqrt(Q2*kappac);
	+ break;
	+ case 2: // assymetric choices - anticolored
	+ kappab=4.(1.-2.sqrt(c)-b+c);
	+ kappac=c+0.25*sqr(1.-b-c+lam)/(kappab-b);
	+ firstQ = sqrt(Q2*kappab);
	+ secondQ = sqrt(Q2*kappac);
	+ break;
	+ }
	+ // calculate the scales
	+ 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
	const Energy2 mc2 = sqr(pc.mass());
	const Energy2 Q2 = -(pb-pc).m2();
	return { 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.));
	const double PROD = 0.25*sqr(1. - a + c + lambda);
	- const double ktilde_c = 0.5*(1-a+c+lambda) + c ;
	- const double ktilde_b = 1.+PROD/(ktilde_c-c) ;
	+ int key = 0;
	+ double ktilde_b, ktilde_c, cosi = 0.;
	+ switch (key)
	+ {
	+ 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 = 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;
	+ }
	return { 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
	const Energy Q = sqrt((p1+p2).m2());
	return {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,256 +1,256 @@
	// -- 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) {}

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

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

	/**
	* Set initial scales for a EW interaction
	*/
	virtual void setInitialEWEvolutionScales(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:
	/**
	* Find the EW partners
	* @param particle The particle to find the partners for
	* @param particles The full set of particles to search
	*/
	vector< pair<double, tShowerParticlePtr> >
	findEWPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles,
	const bool isDecayCase);

	/**
	* 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
	*/
	pair<Energy,Energy> calculateInitialEvolutionScales(const ShowerPPair &,
	- const bool isDecayCase);
	+ const bool isDecayCase, int key);

	/**
	* Calculate the initial evolution scales given momenta
	*/
	pair<Energy,Energy> calculateFinalFinalScales(const Lorentz5Momentum & p1,
	- const Lorentz5Momentum & p2);
	+ const Lorentz5Momentum & p2, int key);

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

	/**
	* Find weakling interacting particles
	*/
	bool weaklyInteracting(tcPDPtr pd) {
	long id = abs(pd->id());
	return ( id==ParticleID::Wplus \|\| id ==ParticleID::Z0 \|\| (id>=1 && id<=6 ) \|\| (id>=11 && id<=16));
	}

	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	- PartnerFinder & operator=(const PartnerFinder &) = delete;
	+ 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_;

	};

	}

	#endif /* HERWIG_PartnerFinder_H */
	diff --git a/Shower/QTilde/Base/ShowerParticle.cc b/Shower/QTilde/Base/ShowerParticle.cc
	--- a/Shower/QTilde/Base/ShowerParticle.cc
	+++ b/Shower/QTilde/Base/ShowerParticle.cc
	@@ -1,579 +1,585 @@
	// -- C++ --
	//
	// ShowerParticle.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 ShowerParticle class.
	//

	#include "ShowerParticle.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include <ostream>

	using namespace Herwig;

	PPtr ShowerParticle::clone() const {
	return new_ptr(*this);
	}

	PPtr ShowerParticle::fullclone() const {
	return new_ptr(*this);
	}

	ClassDescription<ShowerParticle> ShowerParticle::initShowerParticle;
	// Definition of the static class description member.

	void ShowerParticle::vetoEmission(ShowerPartnerType, Energy scale) {
	scales_.QED = min(scale,scales_.QED );
	scales_.QED_noAO = min(scale,scales_.QED_noAO );
	scales_.QCD_c = min(scale,scales_.QCD_c );
	scales_.QCD_c_noAO = min(scale,scales_.QCD_c_noAO );
	scales_.QCD_ac = min(scale,scales_.QCD_ac );
	scales_.QCD_ac_noAO = min(scale,scales_.QCD_ac_noAO);
	- scales_.EW = min(scale,scales_.EW );
	+ scales_.EW_Z = min(scale,scales_.EW_Z );
	+ scales_.EW_W = min(scale,scales_.EW_W );
	if(spinInfo()) spinInfo()->undecay();
	}

	void ShowerParticle::addPartner(EvolutionPartner in ) {
	partners_.push_back(in);
	}

	namespace {

	LorentzRotation boostToShower(Lorentz5Momentum & porig, tShowerBasisPtr basis) {
	LorentzRotation output;
	assert(basis);
	if(basis->frame()==ShowerBasis::BackToBack) {
	// we are doing the evolution in the back-to-back frame
	// work out the boostvector
	Boost boostv(-(basis->pVector()+basis->nVector()).boostVector());
	// momentum of the parton
	Lorentz5Momentum ptest(basis->pVector());
	// construct the Lorentz boost
	output = LorentzRotation(boostv);
	ptest *= output;
	Axis axis(ptest.vect().unit());
	// now rotate so along the z axis as needed for the splitting functions
	if(axis.perp2()>1e-10) {
	double sinth(sqrt(1.-sqr(axis.z())));
	output.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	}
	else if(axis.z()<0.) {
	output.rotate(Constants::pi,Axis(1.,0.,0.));
	}
	- porig *= output;
	+ porig = output*basis->pVector();
	+ porig.setX(ZERO);
	+ porig.setY(ZERO);
	}
	else {
	output = LorentzRotation(-basis->pVector().boostVector());
	- porig *= output;
	+ porig = output*basis->pVector();
	+ porig.setX(ZERO);
	+ porig.setY(ZERO);
	+ porig.setZ(ZERO);
	}
	return output;
	}

	RhoDMatrix bosonMapping(ShowerParticle & particle,
	const Lorentz5Momentum & porig,
	VectorSpinPtr vspin,
	const LorentzRotation & rot,
	Helicity::Direction dir) {
	// rotate the original basis
	vector<LorentzPolarizationVector> sbasis;
	for(unsigned int ix=0;ix<3;++ix) {
	sbasis.push_back(vspin->getProductionBasisState(ix));
	sbasis.back().transform(rot);
	}
	// splitting basis
	vector<LorentzPolarizationVector> fbasis;
	bool massless(particle.id()==ParticleID::g\|\|particle.id()==ParticleID::gamma);
	VectorWaveFunction wave(porig,particle.dataPtr(),dir==outgoing ? incoming : outgoing);
	for(unsigned int ix=0;ix<3;++ix) {
	if(massless&&ix==1) {
	fbasis.push_back(LorentzPolarizationVector());
	}
	else {
	wave.reset(ix,vector_phase);
	fbasis.push_back(wave.wave());
	}
	}
	// work out the mapping
	RhoDMatrix mapping=RhoDMatrix(PDT::Spin1,false);
	for(unsigned int ix=0;ix<3;++ix) {
	for(unsigned int iy=0;iy<3;++iy) {
	mapping(ix,iy)= -conj(sbasis[ix].dot(fbasis[iy]));
	if(particle.id()<0) mapping(ix,iy)=conj(mapping(ix,iy));
	}
	}
	return mapping;
	}

	RhoDMatrix fermionMapping(ShowerParticle & particle,
	const Lorentz5Momentum & porig,
	FermionSpinPtr fspin,
	const LorentzRotation & rot) {
	// extract the original basis states
	vector<LorentzSpinor<SqrtEnergy> > sbasis;
	for(unsigned int ix=0;ix<2;++ix) {
	sbasis.push_back(fspin->getProductionBasisState(ix));
	sbasis.back().transform(rot);
	}
	// calculate the states in the splitting basis
	vector<LorentzSpinor<SqrtEnergy> > fbasis;
	SpinorWaveFunction wave(porig,particle.dataPtr(),
	particle.id()>0 ? incoming : outgoing);
	for(unsigned int ix=0;ix<2;++ix) {
	wave.reset(ix);
	fbasis.push_back(wave.dimensionedWave());
	}
	RhoDMatrix mapping=RhoDMatrix(PDT::Spin1Half,false);
	for(unsigned int ix=0;ix<2;++ix) {
	if(fbasis[0].s2()==complex<SqrtEnergy>()) {
	mapping(ix,0) = sbasis[ix].s3()/fbasis[0].s3();
	mapping(ix,1) = sbasis[ix].s2()/fbasis[1].s2();
	}
	else {
	mapping(ix,0) = sbasis[ix].s2()/fbasis[0].s2();
	mapping(ix,1) = sbasis[ix].s3()/fbasis[1].s3();
	}
	}
	return mapping;
	}

	FermionSpinPtr createFermionSpinInfo(ShowerParticle & particle,
	const Lorentz5Momentum & porig,
	const LorentzRotation & rot,
	Helicity::Direction dir) {
	// calculate the splitting basis for the branching
	// and rotate back to construct the basis states
	LorentzRotation rinv = rot.inverse();
	SpinorWaveFunction wave;
	if(particle.id()>0)
	wave=SpinorWaveFunction(porig,particle.dataPtr(),incoming);
	else
	wave=SpinorWaveFunction(porig,particle.dataPtr(),outgoing);
	FermionSpinPtr fspin = new_ptr(FermionSpinInfo(particle.momentum(),dir==outgoing));
	for(unsigned int ix=0;ix<2;++ix) {
	wave.reset(ix);
	LorentzSpinor<SqrtEnergy> basis = wave.dimensionedWave();
	basis.transform(rinv);
	fspin->setBasisState(ix,basis);
	}
	particle.spinInfo(fspin);
	return fspin;
	}

	VectorSpinPtr createVectorSpinInfo(ShowerParticle & particle,
	const Lorentz5Momentum & porig,
	const LorentzRotation & rot,
	Helicity::Direction dir) {
	// calculate the splitting basis for the branching
	// and rotate back to construct the basis states
	LorentzRotation rinv = rot.inverse();
	bool massless(particle.id()==ParticleID::g\|\|particle.id()==ParticleID::gamma);
	VectorWaveFunction wave(porig,particle.dataPtr(),dir);
	VectorSpinPtr vspin = new_ptr(VectorSpinInfo(particle.momentum(),dir==outgoing));
	for(unsigned int ix=0;ix<3;++ix) {
	LorentzPolarizationVector basis;
	if(massless&&ix==1) {
	basis = LorentzPolarizationVector();
	}
	else {
	wave.reset(ix,vector_phase);
	basis = wave.wave();
	}
	basis *= rinv;
	vspin->setBasisState(ix,basis);
	}
	particle.spinInfo(vspin);
	vspin-> DMatrix() = RhoDMatrix(PDT::Spin1);
	vspin->rhoMatrix() = RhoDMatrix(PDT::Spin1);
	if(massless) {
	vspin-> DMatrix()(0,0) = 0.5;
	vspin->rhoMatrix()(0,0) = 0.5;
	vspin-> DMatrix()(2,2) = 0.5;
	vspin->rhoMatrix()(2,2) = 0.5;
	}
	return vspin;
	}

	}

	RhoDMatrix ShowerParticle::extractRhoMatrix(bool forward) {
	// get the spin density matrix and the mapping
	RhoDMatrix mapping;
	SpinPtr inspin;
	bool needMapping = getMapping(inspin,mapping);
	// set the decayed flag
	inspin->decay();
	// get the spin density matrix
	RhoDMatrix rho = forward ? inspin->rhoMatrix() : inspin->DMatrix();
	// map to the shower basis if needed
	if(needMapping) {
	RhoDMatrix rhop(rho.iSpin(),false);
	for(int ixb=0;ixb<rho.iSpin();++ixb) {
	for(int iyb=0;iyb<rho.iSpin();++iyb) {
	if(mapping(iyb,ixb)==0.)continue;
	for(int iya=0;iya<rho.iSpin();++iya) {
	if(rho(iya,iyb)==0.)continue;
	for(int ixa=0;ixa<rho.iSpin();++ixa) {
	rhop(ixa,ixb) += rho(iya,iyb)mapping(iya,ixa)conj(mapping(iyb,ixb));
	}
	}
	}
	}
	rhop.normalize();
	rho = rhop;
	}
	return rho;
	}

	bool ShowerParticle::getMapping(SpinPtr & output, RhoDMatrix & mapping) {
	// if the particle is not from the hard process
	if(!this->perturbative()) {
	// mapping is the identity
	output=this->spinInfo();
	mapping=RhoDMatrix(this->dataPtr()->iSpin());
	if(output) {
	return false;
	}
	else {
	- Lorentz5Momentum porig=momentum();
	+ Lorentz5Momentum porig;
	LorentzRotation rot = boostToShower(porig,showerBasis());
	Helicity::Direction dir = this->isFinalState() ? outgoing : incoming;
	if(this->dataPtr()->iSpin()==PDT::Spin0) {
	assert(false);
	}
	else if(this->dataPtr()->iSpin()==PDT::Spin1Half) {
	output = createFermionSpinInfo(*this,porig,rot,dir);
	}
	else if(this->dataPtr()->iSpin()==PDT::Spin1) {
	output = createVectorSpinInfo(*this,porig,rot,dir);
	}
	else {
	assert(false);
	}
	return false;
	}
	}
	// if particle is final-state and is from the hard process
	else if(this->isFinalState()) {
	assert(this->perturbative()==1 \|\| this->perturbative()==2);
	// get transform to shower frame
	- Lorentz5Momentum porig=momentum();
	+ Lorentz5Momentum porig;
	LorentzRotation rot = boostToShower(porig,showerBasis());
	// the rest depends on the spin of the particle
	PDT::Spin spin(this->dataPtr()->iSpin());
	mapping=RhoDMatrix(spin,false);
	// do the spin dependent bit
	if(spin==PDT::Spin0) {
	ScalarSpinPtr sspin=dynamic_ptr_cast<ScalarSpinPtr>(this->spinInfo());
	if(!sspin) {
	ScalarWaveFunction::constructSpinInfo(this,outgoing,true);
	}
	output=this->spinInfo();
	return false;
	}
	else if(spin==PDT::Spin1Half) {
	FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
	// spin info exists get information from it
	if(fspin) {
	output=fspin;
	mapping = fermionMapping(*this,porig,fspin,rot);
	return true;
	}
	// spin info does not exist create it
	else {
	output = createFermionSpinInfo(*this,porig,rot,outgoing);
	return false;
	}
	}
	else if(spin==PDT::Spin1) {
	VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
	// spin info exists get information from it
	if(vspin) {
	output=vspin;
	mapping = bosonMapping(*this,porig,vspin,rot,outgoing);
	return true;
	}
	else {
	output = createVectorSpinInfo(*this,porig,rot,outgoing);
	return false;
	}
	}
	// not scalar/fermion/vector
	else
	assert(false);
	}
	// incoming to hard process
	else if(this->perturbative()==1 && !this->isFinalState()) {
	// get the basis vectors
	// get transform to shower frame
	- Lorentz5Momentum porig=momentum();
	+ Lorentz5Momentum porig;
	LorentzRotation rot = boostToShower(porig,showerBasis());
	porig *= this->x();
	// the rest depends on the spin of the particle
	PDT::Spin spin(this->dataPtr()->iSpin());
	mapping=RhoDMatrix(spin);
	// do the spin dependent bit
	if(spin==PDT::Spin0) {
	cerr << "testing spin 0 not yet implemented " << endl;
	assert(false);
	}
	// spin-1/2
	else if(spin==PDT::Spin1Half) {
	FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
	// spin info exists get information from it
	if(fspin) {
	output=fspin;
	mapping = fermionMapping(*this,porig,fspin,rot);
	return true;
	}
	// spin info does not exist create it
	else {
	output = createFermionSpinInfo(*this,porig,rot,incoming);
	return false;
	}
	}
	// spin-1
	else if(spin==PDT::Spin1) {
	VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
	// spinInfo exists map it
	if(vspin) {
	output=vspin;
	mapping = bosonMapping(*this,porig,vspin,rot,incoming);
	return true;
	}
	// create the spininfo
	else {
	output = createVectorSpinInfo(*this,porig,rot,incoming);
	return false;
	}
	}
	assert(false);
	}
	// incoming to decay
	else if(this->perturbative() == 2 && !this->isFinalState()) {
	// get the basis vectors
	- //Lorentz5Momentum porig;
	- //LorentzRotation rot=boostToShower(porig,showerBasis());
	+ Lorentz5Momentum porig;
	+ LorentzRotation rot=boostToShower(porig,showerBasis());
	// the rest depends on the spin of the particle
	PDT::Spin spin(this->dataPtr()->iSpin());
	mapping=RhoDMatrix(spin);
	// do the spin dependent bit
	if(spin==PDT::Spin0) {
	cerr << "testing spin 0 not yet implemented " << endl;
	assert(false);
	}
	// spin-1/2
	else if(spin==PDT::Spin1Half) {
	// FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
	// // spin info exists get information from it
	// if(fspin) {
	// output=fspin;
	// mapping = fermionMapping(*this,porig,fspin,rot);
	// return true;
	// // spin info does not exist create it
	// else {
	// output = createFermionSpinInfo(*this,porig,rot,incoming);
	// return false;
	// }
	// }
	assert(false);
	}
	// // spin-1
	// else if(spin==PDT::Spin1) {
	// VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
	// // spinInfo exists map it
	// if(vspin) {
	// output=vspin;
	// mapping = bosonMapping(*this,porig,vspin,rot);
	// return true;
	// }
	// // create the spininfo
	// else {
	// output = createVectorSpinInfo(*this,porig,rot,incoming);
	// return false;
	// }
	// }
	// assert(false);
	assert(false);
	}
	else
	assert(false);
	return true;
	}

	void ShowerParticle::constructSpinInfo(bool timeLike) {
	// now construct the required spininfo and calculate the basis states
	PDT::Spin spin(dataPtr()->iSpin());
	if(spin==PDT::Spin0) {
	ScalarWaveFunction::constructSpinInfo(this,outgoing,timeLike);
	}
	// calculate the basis states and construct the SpinInfo for a spin-1/2 particle
	else if(spin==PDT::Spin1Half) {
	// outgoing particle
	if(id()>0) {
	vector<LorentzSpinorBar<SqrtEnergy> > stemp;
	SpinorBarWaveFunction::calculateWaveFunctions(stemp,this,outgoing);
	SpinorBarWaveFunction::constructSpinInfo(stemp,this,outgoing,timeLike);
	}
	// outgoing antiparticle
	else {
	vector<LorentzSpinor<SqrtEnergy> > stemp;
	SpinorWaveFunction::calculateWaveFunctions(stemp,this,outgoing);
	SpinorWaveFunction::constructSpinInfo(stemp,this,outgoing,timeLike);
	}
	}
	// calculate the basis states and construct the SpinInfo for a spin-1 particle
	else if(spin==PDT::Spin1) {
	bool massless(id()==ParticleID::g\|\|id()==ParticleID::gamma);
	vector<Helicity::LorentzPolarizationVector> vtemp;
	VectorWaveFunction::calculateWaveFunctions(vtemp,this,outgoing,massless,vector_phase);
	VectorWaveFunction::constructSpinInfo(vtemp,this,outgoing,timeLike,massless);
	}
	else {
	throw Exception() << "Spins higher than 1 are not yet implemented in "
	<< "FS_QtildaShowerKinematics1to2::constructVertex() "
	<< Exception::runerror;
	}
	}

	void ShowerParticle::initializeDecay() {
	Lorentz5Momentum p, n, ppartner, pcm;
	assert(perturbative()!=1);
	// this is for the initial decaying particle
	if(perturbative()==2) {
	ShowerBasisPtr newBasis;
	p = momentum();
	Lorentz5Momentum ppartner(partner()->momentum());
	// removed to make inverse recon work properly
	//if(partner->thePEGBase()) ppartner=partner->thePEGBase()->momentum();
	pcm=ppartner;
	Boost boost(p.findBoostToCM());
	pcm.boost(boost);
	n = Lorentz5Momentum( ZERO,0.5p.mass()pcm.vect().unit());
	n.boost( -boost);
	newBasis = new_ptr(ShowerBasis());
	newBasis->setBasis(p,n,ShowerBasis::Rest);
	showerBasis(newBasis,false);
	}
	else {
	showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0])->showerBasis(),true);
	}
	}

	void ShowerParticle::initializeInitialState(PPtr parent) {
	// For the time being we are considering only 1->2 branching
	Lorentz5Momentum p, n, pthis, pcm;
	assert(perturbative()!=2);
	if(perturbative()==1) {
	ShowerBasisPtr newBasis;
	// find the partner and its momentum
	if(!partner()) return;
	if(partner()->isFinalState()) {
	Lorentz5Momentum pa = -momentum()+partner()->momentum();
	Lorentz5Momentum pb = momentum();
	Energy scale=parent->momentum().t();
	Lorentz5Momentum pbasis(ZERO,parent->momentum().vect().unit()*scale);
	Axis axis(pa.vect().unit());
	LorentzRotation rot;
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	if(axis.perp2()>1e-20) {
	rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	rot.rotateX(Constants::pi);
	}
	if(abs(1.-pa.e()/pa.vect().mag())>1e-6) rot.boostZ( pa.e()/pa.vect().mag());
	pb *= rot;
	if(pb.perp2()/GeV2>1e-20) {
	Boost trans = -1./pb.e()*pb.vect();
	trans.setZ(0.);
	rot.boost(trans);
	}
	pbasis *=rot;
	rot.invert();
	n = rot*Lorentz5Momentum(ZERO,-pbasis.vect());
	p = rot*Lorentz5Momentum(ZERO, pbasis.vect());
	}
	else {
	pcm = parent->momentum();
	p = Lorentz5Momentum(ZERO, pcm.vect());
	n = Lorentz5Momentum(ZERO, -pcm.vect());
	}
	newBasis = new_ptr(ShowerBasis());
	newBasis->setBasis(p,n,ShowerBasis::BackToBack);
	showerBasis(newBasis,false);
	}
	else {
	showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(children()[0])->showerBasis(),true);
	}
	}

	void ShowerParticle::initializeFinalState() {
	// set the basis vectors
	Lorentz5Momentum p,n;
	if(perturbative()!=0) {
	ShowerBasisPtr newBasis;
	// find the partner() and its momentum
	if(!partner()) return;
	Lorentz5Momentum ppartner(partner()->momentum());
	// momentum of the emitting particle
	p = momentum();
	Lorentz5Momentum pcm;
	// if the partner is a final-state particle then the reference
	// vector is along the partner in the rest frame of the pair
	if(partner()->isFinalState()) {
	Boost boost=(p + ppartner).findBoostToCM();
	pcm = ppartner;
	pcm.boost(boost);
	n = Lorentz5Momentum(ZERO,pcm.vect());
	n.boost( -boost);
	}
	else if(!partner()->isFinalState()) {
	// if the partner is an initial-state particle then the reference
	// vector is along the partner which should be massless
	if(perturbative()==1) {
	n = Lorentz5Momentum(ZERO,ppartner.vect());
	}
	// if the partner is an initial-state decaying particle then the reference
	// vector is along the backwards direction in rest frame of decaying particle
	else {
	Boost boost=ppartner.findBoostToCM();
	pcm = p;
	pcm.boost(boost);
	n = Lorentz5Momentum( ZERO, -pcm.vect());
	n.boost( -boost);
	}
	}
	newBasis = new_ptr(ShowerBasis());
	newBasis->setBasis(p,n,ShowerBasis::BackToBack);
	showerBasis(newBasis,false);
	}
	else if(initiatesTLS()) {
	showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0]->children()[0])->showerBasis(),true);
	}
	else {
	showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0])->showerBasis(),true);
	}
	}

	void ShowerParticle::setShowerMomentum(bool timeLike) {
	Energy m = this->mass() > ZERO ? this->mass() : this->data().mass();
	// calculate the momentum of the assuming on-shell
	Energy2 pt2 = sqr(this->showerParameters().pt);
	double alpha = timeLike ? this->showerParameters().alpha : this->x();
	double beta = 0.5(sqr(m) + pt2 - sqr(alpha)showerBasis()->pVector().m2())/(alpha*showerBasis()->p_dot_n());
	Lorentz5Momentum porig=showerBasis()->sudakov2Momentum(alpha,beta,
	this->showerParameters().ptx,
	this->showerParameters().pty);
	porig.setMass(m);
	this->set5Momentum(porig);
	}
	diff --git a/Shower/QTilde/Base/ShowerParticle.h b/Shower/QTilde/Base/ShowerParticle.h
	--- a/Shower/QTilde/Base/ShowerParticle.h
	+++ b/Shower/QTilde/Base/ShowerParticle.h
	@@ -1,535 +1,537 @@
	// -- C++ --
	//
	// ShowerParticle.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_ShowerParticle_H
	#define HERWIG_ShowerParticle_H
	//
	// This is the declaration of the ShowerParticle class.
	//

	#include "ThePEG/EventRecord/Particle.h"
	#include "Herwig/Shower/QTilde/SplittingFunctions/SplittingFunction.fh"
	#include "Herwig/Shower/QTilde/ShowerConfig.h"
	#include "Herwig/Shower/QTilde/Kinematics/ShowerBasis.h"
	#include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.h"
	#include "ShowerParticle.fh"
	#include <iosfwd>

	namespace Herwig {

	using namespace ThePEG;

	/** \ingroup Shower
	* This class represents a particle in the showering process.
	* It inherits from the Particle class of ThePEG and has some
	* specifics information useful only during the showering process.
	*
	* Notice that:
	* - for forward evolution, it is clear what is meant by parent/child;
	* for backward evolution, however, it depends whether we want
	* to keep a physical picture or a Monte-Carlo effective one.
	* In the former case, an incoming particle (emitting particle)
	* splits into an emitted particle and the emitting particle after
	* the emission: the latter two are then children of the
	* emitting particle, the parent. In the Monte-Carlo effective
	* picture, we have that the particle close to the hard subprocess,
	* with higher (space-like) virtuality, splits into an emitted particle
	* and the emitting particle at lower virtuality: the latter two are,
	* in this case, the children of the first one, the parent. However we
	* choose a more physical picture where the new emitting particle is the
	* parented of the emitted final-state particle and the original emitting
	* particle.
	* - the pointer to a SplitFun object is set only in the case
	* that the particle has undergone a shower emission. This is similar to
	* the case of the decay of a normal Particle where
	* the pointer to a Decayer object is set only in the case
	* that the particle has undergone to a decay.
	* In the case of particle connected directly to the hard subprocess,
	* there is no pointer to the hard subprocess, but there is a method
	* isFromHardSubprocess() which returns true only in this case.
	*
	* @see Particle
	* @see ShowerConfig
	* @see ShowerKinematics
	*/
	class ShowerParticle: public Particle {

	public:

	/**
	* Struct for all the info on an evolution partner
	*/
	struct EvolutionPartner {

	/**
	* Constructor
	*/
	EvolutionPartner(tShowerParticlePtr p,double w, ShowerPartnerType t,
	Energy s) : partner(p), weight(w), type(t), scale(s)
	{}

	/**
	* The partner
	*/
	tShowerParticlePtr partner;

	/**
	* Weight
	*/
	double weight;

	/**
	* Type
	*/
	ShowerPartnerType type;

	/**
	* The assoicated evolution scale
	*/
	Energy scale;
	};

	/**
	* Struct to store the evolution scales
	*/
	struct EvolutionScales {

	/**
	* Constructor
	*/
	EvolutionScales() : QED(),QCD_c(),QCD_ac(),
	QED_noAO(),QCD_c_noAO(),QCD_ac_noAO()
	{}

	/**
	* QED scale
	*/
	Energy QED;

	/**
	* QCD colour scale
	*/
	Energy QCD_c;

	/**
	* QCD anticolour scale
	*/
	Energy QCD_ac;

	/**
	* QED scale
	*/
	Energy QED_noAO;

	/**
	* QCD colour scale
	*/
	Energy QCD_c_noAO;

	/**
	* QCD anticolour scale
	*/
	Energy QCD_ac_noAO;

	/**
	- * EW scale
	+ * EW scales
	*/
	- Energy EW;
	+ Energy EW_Z;
	+ Energy EW_W;

	/**
	* EW scale
	*/
	Energy EW_noAO;

	};


	/** @name Construction and descruction functions. */
	//@{

	/**
	* Standard Constructor. Note that the default constructor is
	* private - there is no particle without a pointer to a
	* ParticleData object.
	* @param x the ParticleData object
	* @param fs Whether or not the particle is an inital or final-state particle
	* @param tls Whether or not the particle initiates a time-like shower
	*/
	ShowerParticle(tcEventPDPtr x, bool fs, bool tls=false)
	: Particle(x), _isFinalState(fs),
	_perturbative(0), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
	_vMass(ZERO), _thePEGBase() {}

	/**
	* Copy constructor from a ThePEG Particle
	* @param x ThePEG particle
	* @param pert Where the particle came from
	* @param fs Whether or not the particle is an inital or final-state particle
	* @param tls Whether or not the particle initiates a time-like shower
	*/
	ShowerParticle(const Particle & x, unsigned int pert, bool fs, bool tls=false)
	: Particle(x), _isFinalState(fs),
	_perturbative(pert), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
	_vMass(ZERO), _thePEGBase(&x) {}
	//@}

	public:

	/**
	* Set a preliminary momentum for the particle
	*/
	void setShowerMomentum(bool timelike);

	/**
	* Construct the spin info object for a shower particle
	*/
	void constructSpinInfo(bool timelike);

	/**
	* Perform any initial calculations needed after the branching has been selected
	*/
	void initializeDecay();

	/**
	* Perform any initial calculations needed after the branching has been selected
	* @param parent The beam particle
	*/
	void initializeInitialState(PPtr parent);

	/**
	* Perform any initial calculations needed after the branching has been selected
	*/
	void initializeFinalState();

	/**
	* Access/Set various flags about the state of the particle
	*/
	//@{
	/**
	* Access the flag that tells if the particle is final state
	* or initial state.
	*/
	bool isFinalState() const { return _isFinalState; }

	/**
	* Access the flag that tells if the particle is initiating a
	* time like shower when it has been emitted in an initial state shower.
	*/
	bool initiatesTLS() const { return _initiatesTLS; }

	/**
	* Access the flag which tells us where the particle came from
	* This is 0 for a particle produced in the shower, 1 if the particle came
	* from the hard sub-process and 2 is it came from a decay.
	*/
	unsigned int perturbative() const { return _perturbative; }
	//@}

	/**
	* Set/Get the momentum fraction for initial-state particles
	*/
	//@{
	/**
	* For an initial state particle get the fraction of the beam momentum
	*/
	void x(double x) { _x = x; }

	/**
	* For an initial state particle set the fraction of the beam momentum
	*/
	double x() const { return _x; }
	//@}

	/**
	* Set/Get methods for the ShowerKinematics objects
	*/
	//@{
	/**
	* Access/ the ShowerKinematics object.
	*/
	const ShoKinPtr & showerKinematics() const { return _showerKinematics; }


	/**
	* Set the ShowerKinematics object.
	*/
	void showerKinematics(const ShoKinPtr in) { _showerKinematics = in; }
	//@}

	/**
	* Set/Get methods for the ShowerBasis objects
	*/
	//@{
	/**
	* Access/ the ShowerBasis object.
	*/
	const ShowerBasisPtr & showerBasis() const { return _showerBasis; }


	/**
	* Set the ShowerBasis object.
	*/
	void showerBasis(const ShowerBasisPtr in, bool copy) {
	if(!copy)
	_showerBasis = in;
	else {
	_showerBasis = new_ptr(ShowerBasis());
	_showerBasis->setBasis(in->pVector(),in->nVector(),in->frame());
	}
	}
	//@}

	/**
	* Members relating to the initial evolution scale and partner for the particle
	*/
	//@{
	/**
	* Veto emission at a given scale
	*/
	void vetoEmission(ShowerPartnerType type, Energy scale);

	/**
	* Access to the evolution scales
	*/
	const EvolutionScales & scales() const {return scales_;}

	/**
	* Access to the evolution scales
	*/
	EvolutionScales & scales() {return scales_;}

	/**
	* Return the virtual mass\f$
	*/
	Energy virtualMass() const { return _vMass; }

	/**
	* Set the virtual mass
	*/
	void virtualMass(Energy mass) { _vMass = mass; }

	/**
	* Return the partner
	*/
	tShowerParticlePtr partner() const { return _partner; }

	/**
	* Set the partner
	*/
	void partner(const tShowerParticlePtr partner) { _partner = partner; }

	/**
	* Get the possible partners
	*/
	vector<EvolutionPartner> & partners() { return partners_; }

	/**
	* Add a possible partners
	*/
	void addPartner(EvolutionPartner in );

	/**
	* Clear the evolution partners
	*/
	void clearPartners() { partners_.clear(); }

	/**
	* Return the progenitor of the shower
	*/
	tShowerParticlePtr progenitor() const { return _progenitor; }

	/**
	* Set the progenitor of the shower
	*/
	void progenitor(const tShowerParticlePtr progenitor) { _progenitor = progenitor; }
	//@}


	/**
	* Members to store and provide access to variables for a specific
	* shower evolution scheme
	*/
	//@{
	struct Parameters {
	Parameters() : alpha(1.), beta(), ptx(), pty(), pt() {}
	double alpha;
	double beta;
	Energy ptx;
	Energy pty;
	Energy pt;
	};


	/**
	* Set the vector containing dimensionless variables
	*/
	Parameters & showerParameters() { return _parameters; }
	//@}

	/**
	* If this particle came from the hard process get a pointer to ThePEG particle
	* it came from
	*/
	const tcPPtr thePEGBase() const { return _thePEGBase; }

	public:

	/**
	* Extract the rho matrix including mapping needed in the shower
	*/
	RhoDMatrix extractRhoMatrix(bool forward);

	/**
	* For a particle which came from the hard process get the spin density and
	* the mapping required to the basis used in the Shower
	* @param rho The \f$\rho\f$ matrix
	* @param mapping The mapping
	* @param showerkin The ShowerKinematics object
	*/
	bool getMapping(SpinPtr &, RhoDMatrix & map);

	protected:

	/**
	* Standard clone function.
	*/
	virtual PPtr clone() const;

	/**
	* Standard clone function.
	*/
	virtual PPtr fullclone() const;

	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<ShowerParticle> initShowerParticle;

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

	private:

	/**
	* Whether the particle is in the final or initial state
	*/
	bool _isFinalState;

	/**
	* Whether the particle came from
	*/
	unsigned int _perturbative;

	/**
	* Does a particle produced in the backward shower initiate a time-like shower
	*/
	bool _initiatesTLS;

	/**
	* Dimensionless parameters
	*/
	Parameters _parameters;

	/**
	* The beam energy fraction for particle's in the initial state
	*/
	double _x;

	/**
	* The shower kinematics for the particle
	*/
	ShoKinPtr _showerKinematics;

	/**
	* The shower basis for the particle
	*/
	ShowerBasisPtr _showerBasis;

	/**
	* Storage of the evolution scales
	*/
	EvolutionScales scales_;

	/**
	* Virtual mass
	*/
	Energy _vMass;

	/**
	* Partners
	*/
	tShowerParticlePtr _partner;

	/**
	* Pointer to ThePEG Particle this ShowerParticle was created from
	*/
	const tcPPtr _thePEGBase;

	/**
	* Progenitor
	*/
	tShowerParticlePtr _progenitor;

	/**
	* Partners
	*/
	vector<EvolutionPartner> partners_;

	};

	inline ostream & operator<<(ostream & os, const ShowerParticle::EvolutionScales & es) {
	os << "Scales: QED=" << es.QED / GeV
	<< " QCD_c=" << es.QCD_c / GeV
	<< " QCD_ac=" << es.QCD_ac / GeV
	- << " EW=" << es.EW / GeV
	+ << " EW_Z=" << es.EW_Z / GeV
	+ << " EW_W=" << es.EW_W / GeV
	<< " QED_noAO=" << es.QED_noAO / GeV
	<< " QCD_c_noAO=" << es.QCD_c_noAO / GeV
	<< " QCD_ac_noAO=" << es.QCD_ac_noAO / GeV
	<< " EW_noAO=" << es.EW_noAO / GeV
	<< '\n';
	return os;
	}

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	/** @endcond */

	}

	#endif /* HERWIG_ShowerParticle_H */
	diff --git a/Shower/QTilde/QTildeShowerHandler.cc b/Shower/QTilde/QTildeShowerHandler.cc
	--- a/Shower/QTilde/QTildeShowerHandler.cc
	+++ b/Shower/QTilde/QTildeShowerHandler.cc
	@@ -1,3221 +1,3225 @@
	// -- C++ --
	//
	// QTildeShowerHandler.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 QTildeShowerHandler class.
	//

	#include "QTildeShowerHandler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Parameter.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 "ThePEG/Utilities/EnumIO.h"
	#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
	#include "Herwig/PDF/MPIPDF.h"
	#include "Herwig/PDF/MinBiasPDF.h"
	#include "Herwig/Shower/QTilde/Base/ShowerTree.h"
	#include "Herwig/Shower/QTilde/Base/HardTree.h"
	#include "Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.h"
	#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
	#include "Herwig/PDF/HwRemDecayer.h"
	#include "Herwig/Shower/QTilde/Base/ShowerVertex.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "Herwig/Shower/RealEmissionProcess.h"
	#include "Herwig/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h"
	#include "Herwig/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h"
	#include "Herwig/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h"


	using namespace Herwig;

	bool QTildeShowerHandler::_hardEmissionWarn = true;
	bool QTildeShowerHandler::_missingTruncWarn = true;

	QTildeShowerHandler::QTildeShowerHandler() :
	_maxtry(100), _meCorrMode(1), _reconOpt(0),
	_hardVetoReadOption(false),
	_iptrms(ZERO), _beta(0.), _gamma(ZERO), _iptmax(),
	_limitEmissions(0), _initialenhance(1.), _finalenhance(1.),
	_nReWeight(100), _reWeight(false),
	interaction_(ShowerInteraction::QEDQCD),
	_trunc_Mode(true), _hardEmission(1),
	_softOpt(2), _hardPOWHEG(false), muPt(ZERO)
	{}

	QTildeShowerHandler::~QTildeShowerHandler() {}

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

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

	void QTildeShowerHandler::persistentOutput(PersistentOStream & os) const {
	os << _splittingGenerator << _maxtry
	<< _meCorrMode << _hardVetoReadOption
	<< _limitEmissions << _softOpt << _hardPOWHEG
	<< ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
	<< _vetoes << _fullShowerVetoes << _nReWeight << _reWeight
	<< _trunc_Mode << _hardEmission << _reconOpt
	<< ounit(muPt,GeV) << oenum(interaction_)
	<< _reconstructor << _partnerfinder;
	}

	void QTildeShowerHandler::persistentInput(PersistentIStream & is, int) {
	is >> _splittingGenerator >> _maxtry
	>> _meCorrMode >> _hardVetoReadOption
	>> _limitEmissions >> _softOpt >> _hardPOWHEG
	>> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
	>> _vetoes >> _fullShowerVetoes >> _nReWeight >> _reWeight
	>> _trunc_Mode >> _hardEmission >> _reconOpt
	>> iunit(muPt,GeV) >> ienum(interaction_)
	>> _reconstructor >> _partnerfinder;
	}


	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<QTildeShowerHandler,ShowerHandler>
	describeHerwigQTildeShowerHandler("Herwig::QTildeShowerHandler", "HwShower.so");

	void QTildeShowerHandler::Init() {

	static ClassDocumentation<QTildeShowerHandler> documentation
	("TheQTildeShowerHandler class is the main class"
	" for the angular-ordered parton shower",
	"The Shower evolution was performed using an algorithm described in "
	"\\cite{Marchesini:1983bm,Marchesini:1987cf,Gieseke:2003rz,Bahr:2008pv}.",
	"%\\cite{Marchesini:1983bm}\n"
	"\\bibitem{Marchesini:1983bm}\n"
	" G.~Marchesini and B.~R.~Webber,\n"
	" ``Simulation Of QCD Jets Including Soft Gluon Interference,''\n"
	" Nucl.\\ Phys.\\ B {\\bf 238}, 1 (1984).\n"
	" %%CITATION = NUPHA,B238,1;%%\n"
	"%\\cite{Marchesini:1987cf}\n"
	"\\bibitem{Marchesini:1987cf}\n"
	" G.~Marchesini and B.~R.~Webber,\n"
	" ``Monte Carlo Simulation of General Hard Processes with Coherent QCD\n"
	" Radiation,''\n"
	" Nucl.\\ Phys.\\ B {\\bf 310}, 461 (1988).\n"
	" %%CITATION = NUPHA,B310,461;%%\n"
	"%\\cite{Gieseke:2003rz}\n"
	"\\bibitem{Gieseke:2003rz}\n"
	" S.~Gieseke, P.~Stephens and B.~Webber,\n"
	" ``New formalism for QCD parton showers,''\n"
	" JHEP {\\bf 0312}, 045 (2003)\n"
	" [arXiv:hep-ph/0310083].\n"
	" %%CITATION = JHEPA,0312,045;%%\n"
	);

	static Reference<QTildeShowerHandler,SplittingGenerator>
	interfaceSplitGen("SplittingGenerator",
	"A reference to the SplittingGenerator object",
	&Herwig::QTildeShowerHandler::_splittingGenerator,
	false, false, true, false);

	static Parameter<QTildeShowerHandler,unsigned int> interfaceMaxTry
	("MaxTry",
	"The maximum number of attempts to generate the shower from a"
	" particular ShowerTree",
	&QTildeShowerHandler::_maxtry, 100, 1, 100000,
	false, false, Interface::limited);

	static Parameter<QTildeShowerHandler,unsigned int> interfaceNReWeight
	("NReWeight",
	"The number of attempts for the shower when reweighting",
	&QTildeShowerHandler::_nReWeight, 100, 10, 10000,
	false, false, Interface::limited);

	static Switch<QTildeShowerHandler, unsigned int> ifaceMECorrMode
	("MECorrMode",
	"Choice of the ME Correction Mode",
	&QTildeShowerHandler::_meCorrMode, 1, false, false);
	static SwitchOption on
	(ifaceMECorrMode,"HardPlusSoft","hard+soft on", 1);
	static SwitchOption hard
	(ifaceMECorrMode,"Hard","only hard on", 2);
	static SwitchOption soft
	(ifaceMECorrMode,"Soft","only soft on", 3);

	static Switch<QTildeShowerHandler, bool> ifaceHardVetoReadOption
	("HardVetoReadOption",
	"Apply read-in scale veto to all collisions or just the primary one?",
	&QTildeShowerHandler::_hardVetoReadOption, false, false, false);
	static SwitchOption AllCollisions
	(ifaceHardVetoReadOption,
	"AllCollisions",
	"Read-in pT veto applied to primary and secondary collisions.",
	false);
	static SwitchOption PrimaryCollision
	(ifaceHardVetoReadOption,
	"PrimaryCollision",
	"Read-in pT veto applied to primary but not secondary collisions.",
	true);

	static Parameter<QTildeShowerHandler, Energy> ifaceiptrms
	("IntrinsicPtGaussian",
	"RMS of intrinsic pT of Gaussian distribution:\n"
	"2(1-Beta)exp(-sqr(intrinsicpT/RMS))/sqr(RMS)",
	&QTildeShowerHandler::_iptrms, GeV, ZERO, ZERO, 1000000.0*GeV,
	false, false, Interface::limited);

	static Parameter<QTildeShowerHandler, double> ifacebeta
	("IntrinsicPtBeta",
	"Proportion of inverse quadratic distribution in generating intrinsic pT.\n"
	"(1-Beta) is the proportion of Gaussian distribution",
	&QTildeShowerHandler::_beta, 0, 0, 1,
	false, false, Interface::limited);

	static Parameter<QTildeShowerHandler, Energy> ifacegamma
	("IntrinsicPtGamma",
	"Parameter for inverse quadratic:\n"
	"2BetaGamma/(sqr(Gamma)+sqr(intrinsicpT))",
	&QTildeShowerHandler::_gamma,GeV, ZERO, ZERO, 100000.0*GeV,
	false, false, Interface::limited);

	static Parameter<QTildeShowerHandler, Energy> ifaceiptmax
	("IntrinsicPtIptmax",
	"Upper bound on intrinsic pT for inverse quadratic",
	&QTildeShowerHandler::_iptmax,GeV, ZERO, ZERO, 100000.0*GeV,
	false, false, Interface::limited);

	static RefVector<QTildeShowerHandler,ShowerVeto> ifaceVetoes
	("Vetoes",
	"The vetoes to be checked during showering",
	&QTildeShowerHandler::_vetoes, -1,
	false,false,true,true,false);

	static RefVector<QTildeShowerHandler,FullShowerVeto> interfaceFullShowerVetoes
	("FullShowerVetoes",
	"The vetos to be appliede on the full final state of the shower",
	&QTildeShowerHandler::_fullShowerVetoes, -1, false, false, true, false, false);

	static Switch<QTildeShowerHandler,unsigned int> interfaceLimitEmissions
	("LimitEmissions",
	"Limit the number and type of emissions for testing",
	&QTildeShowerHandler::_limitEmissions, 0, false, false);
	static SwitchOption interfaceLimitEmissionsNoLimit
	(interfaceLimitEmissions,
	"NoLimit",
	"Allow an arbitrary number of emissions",
	0);
	static SwitchOption interfaceLimitEmissionsOneInitialStateEmission
	(interfaceLimitEmissions,
	"OneInitialStateEmission",
	"Allow one emission in the initial state and none in the final state",
	1);
	static SwitchOption interfaceLimitEmissionsOneFinalStateEmission
	(interfaceLimitEmissions,
	"OneFinalStateEmission",
	"Allow one emission in the final state and none in the initial state",
	2);
	static SwitchOption interfaceLimitEmissionsHardOnly
	(interfaceLimitEmissions,
	"HardOnly",
	"Only allow radiation from the hard ME correction",
	3);
	static SwitchOption interfaceLimitEmissionsOneEmission
	(interfaceLimitEmissions,
	"OneEmission",
	"Allow one emission in either the final state or initial state, but not both",
	4);

	static Switch<QTildeShowerHandler,bool> interfaceTruncMode
	("TruncatedShower", "Include the truncated shower?",
	&QTildeShowerHandler::_trunc_Mode, 1, false, false);
	static SwitchOption interfaceTruncMode0
	(interfaceTruncMode,"No","Truncated Shower is OFF", 0);
	static SwitchOption interfaceTruncMode1
	(interfaceTruncMode,"Yes","Truncated Shower is ON", 1);

	static Switch<QTildeShowerHandler,int> interfaceHardEmission
	("HardEmission",
	"Whether to use ME corrections or POWHEG for the hardest emission",
	&QTildeShowerHandler::_hardEmission, 0, false, false);
	static SwitchOption interfaceHardEmissionNone
	(interfaceHardEmission,
	"None",
	"No Corrections",
	0);
	static SwitchOption interfaceHardEmissionMECorrection
	(interfaceHardEmission,
	"MECorrection",
	"Old fashioned ME correction",
	1);
	static SwitchOption interfaceHardEmissionPOWHEG
	(interfaceHardEmission,
	"POWHEG",
	"Powheg style hard emission",
	2);

	static Switch<QTildeShowerHandler,ShowerInteraction> interfaceInteractions
	("Interactions",
	"The interactions to be used in the shower",
	&QTildeShowerHandler::interaction_, ShowerInteraction::QEDQCD, false, false);
	static SwitchOption interfaceInteractionsQCD
	(interfaceInteractions,
	"QCD",
	"Only QCD radiation",
	ShowerInteraction::QCD);
	static SwitchOption interfaceInteractionsQED
	(interfaceInteractions,
	"QED",
	"Only QEd radiation",
	ShowerInteraction::QED);
	static SwitchOption interfaceInteractionEWOnly
	(interfaceInteractions,
	"EWOnly",
	"Only EW",
	ShowerInteraction::EW);
	static SwitchOption interfaceInteractionQEDQCD
	(interfaceInteractions,
	"QEDQCD",
	"QED and QCD",
	ShowerInteraction::QEDQCD);
	static SwitchOption interfaceInteractionALL
	(interfaceInteractions,
	"ALL",
	"QED, QCD and EW",
	ShowerInteraction::ALL);

	static Switch<QTildeShowerHandler,unsigned int> interfaceReconstructionOption
	("ReconstructionOption",
	"Treatment of the reconstruction of the transverse momentum of "
	"a branching from the evolution scale.",
	&QTildeShowerHandler::_reconOpt, 0, false, false);
	static SwitchOption interfaceReconstructionOptionCutOff
	(interfaceReconstructionOption,
	"CutOff",
	"Use the cut-off masses in the calculation",
	0);
	static SwitchOption interfaceReconstructionOptionOffShell5
	(interfaceReconstructionOption,
	"OffShell5",
	"Try and preserve q2 but if pt negative just zero it",
	5);

	static Switch<QTildeShowerHandler,unsigned int> interfaceSoftCorrelations
	("SoftCorrelations",
	"Option for the treatment of soft correlations in the parton shower",
	&QTildeShowerHandler::_softOpt, 2, false, false);
	static SwitchOption interfaceSoftCorrelationsNone
	(interfaceSoftCorrelations,
	"No",
	"No soft correlations",
	0);
	static SwitchOption interfaceSoftCorrelationsFull
	(interfaceSoftCorrelations,
	"Full",
	"Use the full eikonal",
	1);
	static SwitchOption interfaceSoftCorrelationsSingular
	(interfaceSoftCorrelations,
	"Singular",
	"Use original Webber-Marchisini form",
	2);

	static Switch<QTildeShowerHandler,bool> interfaceHardPOWHEG
	("HardPOWHEG",
	"Treatment of powheg emissions which are too hard to have a shower interpretation",
	&QTildeShowerHandler::_hardPOWHEG, false, false, false);
	static SwitchOption interfaceHardPOWHEGAsShower
	(interfaceHardPOWHEG,
	"AsShower",
	"Still interpret as shower emissions",
	false);
	static SwitchOption interfaceHardPOWHEGRealEmission
	(interfaceHardPOWHEG,
	"RealEmission",
	"Generate shower from the real emmission configuration",
	true);

	static Reference<QTildeShowerHandler,KinematicsReconstructor> interfaceKinematicsReconstructor
	("KinematicsReconstructor",
	"Reference to the KinematicsReconstructor object",
	&QTildeShowerHandler::_reconstructor, false, false, true, false, false);

	static Reference<QTildeShowerHandler,PartnerFinder> interfacePartnerFinder
	("PartnerFinder",
	"Reference to the PartnerFinder object",
	&QTildeShowerHandler::_partnerfinder, false, false, true, false, false);

	}

	tPPair QTildeShowerHandler::cascade(tSubProPtr sub,
	XCPtr xcomb) {
	// use me for reference in tex file etc
	useMe();
	prepareCascade(sub);
	// set things up in the base class
	resetWeights();
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	// start of the try block for the whole showering process
	unsigned int countFailures=0;
	while (countFailures<maxtry()) {
	try {
	decay_.clear();
	done_.clear();
	PerturbativeProcessPtr hard;
	DecayProcessMap decay;
	splitHardProcess(firstInteraction() ? tagged() :
	tPVector(currentSubProcess()->outgoing().begin(),
	currentSubProcess()->outgoing().end()),
	hard,decay);
	ShowerTree::constructTrees(hard_,decay_,hard,decay);
	// if no hard process
	if(!hard_) throw Exception() << "Shower starting with a decay"
	<< "is not implemented"
	<< Exception::runerror;
	// perform the shower for the hard process
	showerHardProcess(hard_,xcomb);
	done_.push_back(hard_);
	hard_->updateAfterShower(decay_);
	// if no decaying particles to shower break out of the loop
	if(decay_.empty()) break;
	// shower the decay products
	while(!decay_.empty()) {
	// find particle whose production process has been showered
	ShowerDecayMap::iterator dit = decay_.begin();
	while(!dit->second->parent()->hasShowered() && dit!=decay_.end()) ++dit;
	assert(dit!=decay_.end());
	// get the particle
	ShowerTreePtr decayingTree = dit->second;
	// remove it from the multimap
	decay_.erase(dit);
	// make sure the particle has been decayed
	QTildeShowerHandler::decay(decayingTree,decay_);
	// now shower the decay
	showerDecay(decayingTree);
	done_.push_back(decayingTree);
	decayingTree->updateAfterShower(decay_);
	}
	// suceeded break out of the loop
	break;
	}
	catch (KinematicsReconstructionVeto) {
	resetWeights();
	++countFailures;
	}
	catch ( ... ) {
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	throw;
	}
	}
	// if loop exited because of too many tries, throw event away
	if (countFailures >= maxtry()) {
	resetWeights();
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	throw Exception() << "Too many tries for main while loop "
	<< "in QTildeShowerHandler::cascade()."
	<< Exception::eventerror;
	}
	//enter the particles in the event record
	fillEventRecord();
	// clear storage
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	// non hadronic case return
	if (!isResolvedHadron(incomingBeams().first ) &&
	!isResolvedHadron(incomingBeams().second) )
	return incomingBeams();
	// remake the remnants (needs to be after the colours are sorted
	// out in the insertion into the event record)
	if ( firstInteraction() ) return remakeRemnant(sub->incoming());
	//Return the new pair of incoming partons. remakeRemnant is not
	//necessary here, because the secondary interactions are not yet
	//connected to the remnants.
	return make_pair(findFirstParton(sub->incoming().first ),
	findFirstParton(sub->incoming().second));
	}

	void QTildeShowerHandler::fillEventRecord() {
	// create a new step
	StepPtr pstep = newStep();
	assert(!done_.empty());
	assert(done_[0]->isHard());
	// insert the steps
	for(unsigned int ix=0;ix<done_.size();++ix) {
	done_[ix]->fillEventRecord(pstep,doISR(),doFSR());
	}
	}

	HardTreePtr QTildeShowerHandler::generateCKKW(ShowerTreePtr ) const {
	return HardTreePtr();
	}

	void QTildeShowerHandler::doinit() {
	ShowerHandler::doinit();
	// check on the reweighting
	for(unsigned int ix=0;ix<_fullShowerVetoes.size();++ix) {
	if(_fullShowerVetoes[ix]->behaviour()==1) {
	_reWeight = true;
	break;
	}
	}
	if(_reWeight && maximumTries()<_nReWeight) {
	throw Exception() << "Reweight being performed in the shower but the number of attempts for the"
	<< "shower is less than that for the reweighting.\n"
	<< "Maximum number of attempt for the shower "
	<< fullName() << ":MaxTry is " << maximumTries() << "\nand for reweighting is "
	<< fullName() << ":NReWeight is " << _nReWeight << "\n"
	<< "we recommend the number of attempts is 10 times the number for reweighting\n"
	<< Exception::runerror;
	}
	ShowerTree::_vmin2 = vMin();
	ShowerTree::_spaceTime = includeSpaceTime();
	}

	void QTildeShowerHandler::doinitrun() {
	ShowerHandler::doinitrun();
	ShowerTree::_vmin2 = vMin();
	ShowerTree::_spaceTime = includeSpaceTime();
	}

	void QTildeShowerHandler::generateIntrinsicpT(vector<ShowerProgenitorPtr> particlesToShower) {
	_intrinsic.clear();
	if ( !ipTon() \|\| !doISR() ) return;
	// don't do anything for the moment for secondary scatters
	if( !firstInteraction() ) return;
	// generate intrinsic pT
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	// only consider initial-state particles
	if(particlesToShower[ix]->progenitor()->isFinalState()) continue;
	if(!particlesToShower[ix]->progenitor()->dataPtr()->coloured()) continue;
	Energy ipt;
	if(UseRandom::rnd() > _beta) {
	ipt=_iptrms*sqrt(-log(UseRandom::rnd()));
	}
	else {
	ipt=_gamma*sqrt(pow(1.+sqr(_iptmax/_gamma), UseRandom::rnd())-1.);
	}
	pair<Energy,double> pt = make_pair(ipt,UseRandom::rnd(Constants::twopi));
	_intrinsic[particlesToShower[ix]] = pt;
	}
	}

	void QTildeShowerHandler::setupMaximumScales(const vector<ShowerProgenitorPtr> & p,
	XCPtr xcomb) {
	// let POWHEG events radiate freely
	if(_hardEmission==2&&hardTree()) {
	vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
	for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(Constants::MaxEnergy);
	return;
	}
	// return if no vetos
	if (!restrictPhasespace()) return;
	// find out if hard partonic subprocess.
	bool isPartonic(false);
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit = _currenttree->incomingLines().begin();
	Lorentz5Momentum pcm;
	for(; cit!=currentTree()->incomingLines().end(); ++cit) {
	pcm += cit->first->progenitor()->momentum();
	isPartonic \|= cit->first->progenitor()->coloured();
	}
	// find minimum pt from hard process, the maximum pt from all outgoing
	// coloured lines (this is simpler and more general than
	// 2stu/(s^2+t^2+u^2)). Maximum scale for scattering processes will
	// be transverse mass.
	Energy ptmax = generator()->maximumCMEnergy();
	// general case calculate the scale
	if ( !hardScaleIsMuF() \|\| (hardVetoReadOption()&&!firstInteraction()) ) {
	// scattering process
	if(currentTree()->isHard()) {
	assert(xcomb);
	// coloured incoming particles
	if (isPartonic) {
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt = currentTree()->outgoingLines().begin();
	for(; cjt!=currentTree()->outgoingLines().end(); ++cjt) {
	if (cjt->first->progenitor()->coloured())
	ptmax = min(ptmax,cjt->first->progenitor()->momentum().mt());
	}
	}
	if (ptmax == generator()->maximumCMEnergy() ) ptmax = pcm.m();
	if(hardScaleIsMuF()&&hardVetoReadOption()&&
	!firstInteraction()) {
	ptmax=min(ptmax,sqrt(xcomb->lastShowerScale()));
	}
	}
	// decay, incoming() is the decaying particle.
	else {
	ptmax = currentTree()->incomingLines().begin()->first
	->progenitor()->momentum().mass();
	}
	}
	// hepeup.SCALUP is written into the lastXComb by the
	// LesHouchesReader itself - use this by user's choice.
	// Can be more general than this.
	else {
	if(currentTree()->isHard()) {
	assert(xcomb);
	ptmax = sqrt( xcomb->lastShowerScale() );
	}
	else {
	ptmax = currentTree()->incomingLines().begin()->first
	->progenitor()->momentum().mass();
	}
	}
	ptmax *= hardScaleFactor();
	// set maxHardPt for all progenitors. For partonic processes this
	// is now the max pt in the FS, for non-partonic processes or
	// processes with no coloured FS the invariant mass of the IS
	vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
	for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(ptmax);
	}

	void QTildeShowerHandler::setupHardScales(const vector<ShowerProgenitorPtr> & p,
	XCPtr xcomb) {
	if ( hardScaleIsMuF() &&
	(!hardVetoReadOption() \|\| firstInteraction()) ) {
	Energy hardScale = ZERO;
	if(currentTree()->isHard()) {
	assert(xcomb);
	hardScale = sqrt( xcomb->lastShowerScale() );
	}
	else {
	hardScale = currentTree()->incomingLines().begin()->first
	->progenitor()->momentum().mass();
	}
	hardScale *= hardScaleFactor();
	vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
	for (; ckt != p.end(); ckt++) (*ckt)->hardScale(hardScale);
	muPt = hardScale;
	}
	}

	void QTildeShowerHandler::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) {
	_hardme = HwMEBasePtr();
	// extract the matrix element
	tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(xcomb);
	if(lastXC) {
	_hardme = dynamic_ptr_cast<HwMEBasePtr>(lastXC->matrixElement());
	}
	_decayme = HwDecayerBasePtr();
	// set the current tree
	currentTree(hard);
	hardTree(HardTreePtr());
	// work out the type of event
	currentTree()->xcombPtr(dynamic_ptr_cast<StdXCombPtr>(xcomb));
	currentTree()->identifyEventType();
	checkFlags();
	// generate the showering
	doShowering(true,xcomb);
	}

	RealEmissionProcessPtr QTildeShowerHandler::hardMatrixElementCorrection(bool hard) {
	// set the initial enhancement factors for the soft correction
	_initialenhance = 1.;
	_finalenhance = 1.;
	// see if we can get the correction from the matrix element
	// or decayer
	RealEmissionProcessPtr real;
	if(hard) {
	if(_hardme&&_hardme->hasMECorrection()) {
	_hardme->initializeMECorrection(_currenttree->perturbativeProcess(),
	_initialenhance,_finalenhance);
	if(hardMEC())
	real =
	_hardme->applyHardMatrixElementCorrection(_currenttree->perturbativeProcess());
	}
	}
	else {
	if(_decayme&&_decayme->hasMECorrection()) {
	_decayme->initializeMECorrection(_currenttree->perturbativeProcess(),
	_initialenhance,_finalenhance);
	if(hardMEC())
	real = _decayme->applyHardMatrixElementCorrection(_currenttree->perturbativeProcess());
	}
	}
	return real;
	}

	ShowerParticleVector QTildeShowerHandler::createTimeLikeChildren(tShowerParticlePtr, IdList ids) {
	// Create the ShowerParticle objects for the two children of
	// the emitting particle; set the parent/child relationship
	// if same as definition create particles, otherwise create cc
	ShowerParticleVector children;
	for(unsigned int ix=0;ix<2;++ix) {
	children.push_back(new_ptr(ShowerParticle(ids[ix+1],true)));
	if(children[ix]->id()==_progenitor->id()&&!ids[ix+1]->stable()&&abs(ids[ix+1]->id())!=ParticleID::tauminus)
	children[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
	else
	children[ix]->set5Momentum(Lorentz5Momentum(ids[ix+1]->mass()));
	}
	return children;
	}

	bool QTildeShowerHandler::timeLikeShower(tShowerParticlePtr particle,
	ShowerInteraction type,
	Branching fb, bool first) {
	// don't do anything if not needed
	if(_limitEmissions == 1 \|\| hardOnly() \|\|
	( _limitEmissions == 2 && _nfs != 0) \|\|
	( _limitEmissions == 4 && _nfs + _nis != 0) ) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	// generate the emission
	ShowerParticleVector children;
	// generate the emission
	if(!fb.kinematics)
	fb = selectTimeLikeBranching(particle,type,HardBranchingPtr());
	// no emission, return
	if(!fb.kinematics) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	Branching fc[2] = {Branching(),Branching()};

	assert(fb.kinematics);
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	// check highest pT
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the children
	children = createTimeLikeChildren(particle,fb.ids);
	// update the children
	particle->showerKinematics()->
	updateChildren(particle, children,fb.type);
	// update number of emissions
	++_nfs;
	if(_limitEmissions!=0) {
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}
	// select branchings for children
	fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
	fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
	// shower the first particle
	if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	if(_reconOpt>=1)
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	// branching has happened
	if(first&&!children.empty())
	particle->showerKinematics()->resetChildren(particle,children);
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}

	bool
	QTildeShowerHandler::spaceLikeShower(tShowerParticlePtr particle, PPtr beam,
	ShowerInteraction type) {
	//using the pdf's associated with the ShowerHandler assures, that
	//modified pdf's are used for the secondary interactions via
	//CascadeHandler::resetPDFs(...)
	tcPDFPtr pdf;
	if(firstPDF().particle() == _beam)
	pdf = firstPDF().pdf();
	if(secondPDF().particle() == _beam)
	pdf = secondPDF().pdf();
	Energy freeze = pdfFreezingScale();
	// don't do anything if not needed
	if(_limitEmissions == 2 \|\| hardOnly() \|\|
	( _limitEmissions == 1 && _nis != 0 ) \|\|
	( _limitEmissions == 4 && _nis + _nfs != 0 ) ) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	Branching bb;
	// generate branching
	while (true) {
	bb=_splittingGenerator->chooseBackwardBranching(*particle,beam,
	_initialenhance,
	_beam,type,
	pdf,freeze);
	// return if no emission
	if(!bb.kinematics) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	// if not vetoed break
	if(!spaceLikeVetoed(bb,particle)) break;
	// otherwise reset scale and continue
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	}
	// assign the splitting function and shower kinematics
	particle->showerKinematics(bb.kinematics);
	if(bb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(bb.kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// particles as in Sudakov form factor
	tcPDPtr part[2]={bb.ids[0],bb.ids[2]};
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent = new_ptr(ShowerParticle(part[0],false));
	ShowerParticlePtr otherChild = new_ptr(ShowerParticle(part[1],true,true));
	ShowerParticleVector theChildren;
	theChildren.push_back(particle);
	theChildren.push_back(otherChild);
	//this updates the evolution scale
	particle->showerKinematics()->
	updateParent(newParent, theChildren,bb.type);
	// update the history if needed
	_currenttree->updateInitialStateShowerProduct(_progenitor,newParent);
	_currenttree->addInitialStateBranching(particle,newParent,otherChild);
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	++_nis;
	bool emitted = _limitEmissions==0 ?
	spaceLikeShower(newParent,beam,type) : false;
	if(newParent->spinInfo()) newParent->spinInfo()->develop();
	// now reconstruct the momentum
	if(!emitted) {
	if(_intrinsic.find(_progenitor)==_intrinsic.end()) {
	bb.kinematics->updateLast(newParent,ZERO,ZERO);
	}
	else {
	pair<Energy,double> kt=_intrinsic[_progenitor];
	bb.kinematics->updateLast(newParent,
	kt.first*cos(kt.second),
	kt.first*sin(kt.second));
	}
	}
	particle->showerKinematics()->
	updateChildren(newParent, theChildren,bb.type);
	if(_limitEmissions!=0) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}
	// perform the shower of the final-state particle
	timeLikeShower(otherChild,type,Branching(),true);
	updateHistory(otherChild);
	if(theChildren[1]->spinInfo()) theChildren[1]->spinInfo()->develop();
	// return the emitted
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}

	void QTildeShowerHandler::showerDecay(ShowerTreePtr decay) {
	// work out the type of event
	currentTree()->xcombPtr(StdXCombPtr());
	currentTree()->identifyEventType();
	_decayme = HwDecayerBasePtr();
	_hardme = HwMEBasePtr();
	// find the decayer
	// try the normal way if possible
	tDMPtr dm = decay->incomingLines().begin()->first->original() ->decayMode();
	if(!dm) dm = decay->incomingLines().begin()->first->copy() ->decayMode();
	if(!dm) dm = decay->incomingLines().begin()->first->progenitor()->decayMode();
	// otherwise make a string and look it up
	if(!dm) {
	string tag = decay->incomingLines().begin()->first->original()->dataPtr()->name()
	+ "->";
	OrderedParticles outgoing;
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) {
	if(abs(decay->incomingLines().begin()->first->original()->id()) == ParticleID::t &&
	abs(it->first->original()->id())==ParticleID::Wplus &&
	decay->treelinks().size() == 1) {
	ShowerTreePtr Wtree = decay->treelinks().begin()->first;
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	it2=Wtree->outgoingLines().begin();it2!=Wtree->outgoingLines().end();++it2) {
	outgoing.insert(it2->first->original()->dataPtr());
	}
	}
	else {
	outgoing.insert(it->first->original()->dataPtr());
	}
	}
	for(OrderedParticles::const_iterator it=outgoing.begin(); it!=outgoing.end();++it) {
	if(it!=outgoing.begin()) tag += ",";
	tag +=(**it).name();
	}
	tag += ";";
	dm = findDecayMode(tag);
	}
	if(dm) _decayme = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
	// set the ShowerTree to be showered
	currentTree(decay);
	decay->applyTransforms();
	hardTree(HardTreePtr());
	// generate the showering
	doShowering(false,XCPtr());
	// if no vetos
	// force calculation of spin correlations
	SpinPtr spInfo = decay->incomingLines().begin()->first->progenitor()->spinInfo();
	if(spInfo) {
	if(!spInfo->developed()) spInfo->needsUpdate();
	spInfo->develop();
	}
	}

	bool QTildeShowerHandler::spaceLikeDecayShower(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass,ShowerInteraction type,
	Branching fb) {
	// don't do anything if not needed
	if(_limitEmissions == 1 \|\| hardOnly() \|\|
	( _limitEmissions == 3 && _nis != 0) \|\|
	( _limitEmissions == 4 && _nfs + _nis != 0) ) {
	return false;
	}
	// generate the emission
	ShowerParticleVector children;
	// generate the emission
	if(!fb.kinematics)
	fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,
	HardBranchingPtr());
	// no emission, return
	if(!fb.kinematics) return false;
	Branching fc[2];
	if(particle->virtualMass()==ZERO)
	particle->virtualMass(_progenitor->progenitor()->mass());
	fc[0] = Branching();
	fc[1] = Branching();
	assert(fb.kinematics);
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the ShowerParticle objects for the two children
	children = createTimeLikeChildren(particle,fb.ids);
	// updateChildren the children
	particle->showerKinematics()->
	updateChildren(particle, children, fb.type);
	// select branchings for children
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,
	type,HardBranchingPtr());
	fc[1] = selectTimeLikeBranching (children[1],type,HardBranchingPtr());
	// old default
	++_nis;
	// shower the first particle
	_currenttree->updateInitialStateShowerProduct(_progenitor,children[0]);
	_currenttree->addInitialStateBranching(particle,children[0],children[1]);
	if(fc[0].kinematics) spaceLikeDecayShower(children[0],maxScales,minmass,type,Branching());
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],true);
	updateHistory(children[1]);
	// TODO NEED AN UPDATE HERE FOR RECONOPT!=0
	// branching has happened
	return true;
	}

	vector<ShowerProgenitorPtr> QTildeShowerHandler::setupShower(bool hard) {
	RealEmissionProcessPtr real;
	// generate hard me if needed
	if(_hardEmission==1) {
	real = hardMatrixElementCorrection(hard);
	if(real&&!real->outgoing().empty()) setupMECorrection(real);
	}
	// generate POWHEG hard emission if needed
	else if(_hardEmission==2)
	hardestEmission(hard);
	// set the initial colour partners
	setEvolutionPartners(hard,interaction_,false);
	// get the particles to be showered
	vector<ShowerProgenitorPtr> particlesToShower =
	currentTree()->extractProgenitors();
	// return the answer
	return particlesToShower;
	}

	void QTildeShowerHandler::setEvolutionPartners(bool hard,ShowerInteraction type,
	bool clear) {
	// match the particles in the ShowerTree and hardTree
	if(hardTree() && !hardTree()->connect(currentTree()))
	throw Exception() << "Can't match trees in "
	<< "QTildeShowerHandler::setEvolutionPartners()"
	<< Exception::eventerror;
	// extract the progenitors
	vector<ShowerParticlePtr> particles =
	currentTree()->extractProgenitorParticles();
	// clear the partners if needed
	if(clear) {
	for(unsigned int ix=0;ix<particles.size();++ix) {
	particles[ix]->partner(ShowerParticlePtr());
	particles[ix]->clearPartners();
	}
	}
	// sort out the colour partners
	if(hardTree()) {
	// find the partner
	for(unsigned int ix=0;ix<particles.size();++ix) {
	tShowerParticlePtr partner = hardTree()->particles()[particles[ix]]->branchingParticle()->partner();
	if(!partner) continue;
	for(map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	it=hardTree()->particles().begin();
	it!=hardTree()->particles().end();++it) {
	if(it->second->branchingParticle()==partner) {
	particles[ix]->partner(it->first);
	break;
	}
	}
	if(!particles[ix]->partner())
	throw Exception() << "Can't match partners in "
	<< "QTildeShowerHandler::setEvolutionPartners()"
	<< Exception::eventerror;
	}
	}

	// Set the initial evolution scales
	partnerFinder()->
	setInitialEvolutionScales(particles,!hard,interaction_,!_hardtree);
	if(hardTree() && _hardPOWHEG) {
	bool tooHard=false;
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit=hardTree()->particles().end();
	for(unsigned int ix=0;ix<particles.size();++ix) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	mit = hardTree()->particles().find(particles[ix]);
	Energy hardScale(ZERO);
	ShowerPartnerType type(ShowerPartnerType::Undefined);
	// final-state
	if(particles[ix]->isFinalState()) {
	if(mit!= eit && !mit->second->children().empty()) {
	hardScale = mit->second->scale();
	type = mit->second->type();
	}
	}
	// initial-state
	else {
	if(mit!= eit && mit->second->parent()) {
	hardScale = mit->second->parent()->scale();
	type = mit->second->parent()->type();
	}
	}
	if(type!=ShowerPartnerType::Undefined) {
	if(type==ShowerPartnerType::QED) {
	tooHard \|= particles[ix]->scales().QED_noAO<hardScale;
	}
	else if(type==ShowerPartnerType::QCDColourLine) {
	tooHard \|= particles[ix]->scales().QCD_c_noAO<hardScale;
	}
	else if(type==ShowerPartnerType::QCDAntiColourLine) {
	tooHard \|= particles[ix]->scales().QCD_ac_noAO<hardScale;
	}
	- else if(type==ShowerPartnerType::EW) {
	- tooHard \|= particles[ix]->scales().EW<hardScale;
	+ else if(type==ShowerPartnerType::EW &&
	+ (particles[1]->id()==ParticleID::Z0
	+ \|\| particles[2]->id()==ParticleID::Z0)) {
	+ tooHard \|= particles[ix]->scales().EW_Z<hardScale;
	+ }
	+ else if(type==ShowerPartnerType::EW &&
	+ (abs(particles[1]->id())==ParticleID::Wplus
	+ \|\| abs(particles[2]->id())==ParticleID::Wplus)) {
	+ tooHard \|= particles[ix]->scales().EW_W<hardScale;
	}
	}
	}
	if(tooHard) convertHardTree(hard,type);
	}
	}

	void QTildeShowerHandler::updateHistory(tShowerParticlePtr particle) {
	if(!particle->children().empty()) {
	ShowerParticleVector theChildren;
	for(unsigned int ix=0;ix<particle->children().size();++ix) {
	ShowerParticlePtr part = dynamic_ptr_cast<ShowerParticlePtr>
	(particle->children()[ix]);
	theChildren.push_back(part);
	}
	// update the history if needed
	if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
	_currenttree->updateFinalStateShowerProduct(_progenitor,
	particle,theChildren);
	_currenttree->addFinalStateBranching(particle,theChildren);
	for(unsigned int ix=0;ix<theChildren.size();++ix)
	updateHistory(theChildren[ix]);
	}
	}

	bool QTildeShowerHandler::startTimeLikeShower(ShowerInteraction type) {
	// initialize basis vectors etc
	if(!progenitor()->progenitor()->partner()) return false;
	progenitor()->progenitor()->initializeFinalState();
	if(hardTree()) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit=hardTree()->particles().end(),
	mit = hardTree()->particles().find(progenitor()->progenitor());
	if( mit != eit && !mit->second->children().empty() ) {
	bool output=truncatedTimeLikeShower(progenitor()->progenitor(),
	mit->second ,type,Branching(),true);
	if(output) updateHistory(progenitor()->progenitor());
	return output;
	}
	}
	// do the shower
	bool output = hardOnly() ? false :
	timeLikeShower(progenitor()->progenitor() ,type,Branching(),true) ;
	if(output) updateHistory(progenitor()->progenitor());
	return output;
	}

	bool QTildeShowerHandler::startSpaceLikeShower(PPtr parent, ShowerInteraction type) {
	// initialise the basis vectors
	if(!progenitor()->progenitor()->partner()) return false;
	progenitor()->progenitor()->initializeInitialState(parent);
	if(hardTree()) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit =hardTree()->particles().end(),
	mit = hardTree()->particles().find(progenitor()->progenitor());
	if( mit != eit && mit->second->parent() ) {
	return truncatedSpaceLikeShower( progenitor()->progenitor(),
	parent, mit->second->parent(), type );
	}
	}
	// perform the shower
	return hardOnly() ? false :
	spaceLikeShower(progenitor()->progenitor(),parent,type);
	}

	bool QTildeShowerHandler::
	startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
	Energy minimumMass,ShowerInteraction type) {
	// set up the particle basis vectors
	if(!progenitor()->progenitor()->partner()) return false;
	progenitor()->progenitor()->initializeDecay();
	if(hardTree()) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit =hardTree()->particles().end(),
	mit = hardTree()->particles().find(progenitor()->progenitor());
	if( mit != eit && mit->second->parent() ) {
	HardBranchingPtr branch=mit->second;
	while(branch->parent()) branch=branch->parent();
	return truncatedSpaceLikeDecayShower(progenitor()->progenitor(),maxScales,
	minimumMass, branch ,type, Branching());
	}
	}
	// perform the shower
	return hardOnly() ? false :
	spaceLikeDecayShower(progenitor()->progenitor(),maxScales,minimumMass,type,Branching());
	}

	bool QTildeShowerHandler::timeLikeVetoed(const Branching & fb,
	ShowerParticlePtr particle) {
	// work out type of interaction
	ShowerInteraction type = convertInteraction(fb.type);
	// check whether emission was harder than largest pt of hard subprocess
	if ( restrictPhasespace() && fb.kinematics->pT() > _progenitor->maxHardPt() )
	return true;
	// soft matrix element correction veto
	if( softMEC()) {
	if(_hardme && _hardme->hasMECorrection()) {
	if(_hardme->softMatrixElementVeto(particle,
	_progenitor->progenitor(),
	particle->isFinalState(),
	_progenitor->highestpT(),
	fb.ids, fb.kinematics->z(),
	fb.kinematics->scale(),
	fb.kinematics->pT()))
	return true;
	}
	else if(_decayme && _decayme->hasMECorrection()) {
	if(_decayme->softMatrixElementVeto(particle,
	_progenitor->progenitor(),
	particle->isFinalState(),
	_progenitor->highestpT(),
	fb.ids, fb.kinematics->z(),
	fb.kinematics->scale(),
	fb.kinematics->pT()))
	return true;
	}
	}
	// veto on maximum pt
	if(fb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
	// general vetos
	if (fb.kinematics && !_vetoes.empty()) {
	bool vetoed=false;
	for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
	v != _vetoes.end(); ++v) {
	bool test = (**v).vetoTimeLike(_progenitor,particle,fb,currentTree());
	switch((**v).vetoType()) {
	case ShowerVeto::Emission:
	vetoed \|= test;
	break;
	case ShowerVeto::Shower:
	if(test) throw VetoShower();
	break;
	case ShowerVeto::Event:
	if(test) throw Veto();
	break;
	}
	}
	if(vetoed) return true;
	}
	if ( firstInteraction() &&
	profileScales() ) {
	double weight =
	profileScales()->
	hardScaleProfile(_progenitor->hardScale(),fb.kinematics->pT());
	if ( UseRandom::rnd() > weight )
	return true;
	}
	return false;
	}

	bool QTildeShowerHandler::spaceLikeVetoed(const Branching & bb,
	ShowerParticlePtr particle) {
	// work out type of interaction
	ShowerInteraction type = convertInteraction(bb.type);
	// check whether emission was harder than largest pt of hard subprocess
	if (restrictPhasespace() && bb.kinematics->pT() > _progenitor->maxHardPt())
	return true;
	// apply the soft correction
	if( softMEC() && _hardme && _hardme->hasMECorrection() ) {
	if(_hardme->softMatrixElementVeto(particle,
	_progenitor->progenitor(),
	particle->isFinalState(),
	_progenitor->highestpT(),
	bb.ids, bb.kinematics->z(),
	bb.kinematics->scale(),
	bb.kinematics->pT()))
	return true;
	}
	// the more general vetos

	// check vs max pt for the shower
	if(bb.kinematics->pT()>_progenitor->maximumpT(type)) return true;

	if (!_vetoes.empty()) {
	bool vetoed=false;
	for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
	v != _vetoes.end(); ++v) {
	bool test = (**v).vetoSpaceLike(_progenitor,particle,bb,currentTree());
	switch ((**v).vetoType()) {
	case ShowerVeto::Emission:
	vetoed \|= test;
	break;
	case ShowerVeto::Shower:
	if(test) throw VetoShower();
	break;
	case ShowerVeto::Event:
	if(test) throw Veto();
	break;
	}
	}
	if (vetoed) return true;
	}
	if ( firstInteraction() &&
	profileScales() ) {
	double weight =
	profileScales()->
	hardScaleProfile(_progenitor->hardScale(),bb.kinematics->pT());
	if ( UseRandom::rnd() > weight )
	return true;
	}
	return false;
	}

	bool QTildeShowerHandler::spaceLikeDecayVetoed( const Branching & fb,
	ShowerParticlePtr particle) {
	// work out type of interaction
	ShowerInteraction type = convertInteraction(fb.type);
	// apply the soft correction
	if( softMEC() && _decayme && _decayme->hasMECorrection() ) {
	if(_decayme->softMatrixElementVeto(particle,
	_progenitor->progenitor(),
	particle->isFinalState(),
	_progenitor->highestpT(),
	fb.ids, fb.kinematics->z(),
	fb.kinematics->scale(),
	fb.kinematics->pT()))
	return true;
	}
	// veto on hardest pt in the shower
	if(fb.kinematics->pT()> _progenitor->maximumpT(type)) return true;
	// general vetos
	if (!_vetoes.empty()) {
	bool vetoed=false;
	for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
	v != _vetoes.end(); ++v) {
	bool test = (**v).vetoSpaceLike(_progenitor,particle,fb,currentTree());
	switch((**v).vetoType()) {
	case ShowerVeto::Emission:
	vetoed \|= test;
	break;
	case ShowerVeto::Shower:
	if(test) throw VetoShower();
	break;
	case ShowerVeto::Event:
	if(test) throw Veto();
	break;
	}
	if (vetoed) return true;
	}
	}
	return false;
	}

	void QTildeShowerHandler::hardestEmission(bool hard) {
	HardTreePtr ISRTree;
	// internal POWHEG in production or decay
	if( (( _hardme && _hardme->hasPOWHEGCorrection()!=0 ) \|\|
	( _decayme && _decayme->hasPOWHEGCorrection()!=0 ) ) ) {
	RealEmissionProcessPtr real;
	unsigned int type(0);
	// production
	if(_hardme) {
	assert(hard);
	real = _hardme->generateHardest( currentTree()->perturbativeProcess(),
	interaction_);
	type = _hardme->hasPOWHEGCorrection();
	}
	// decay
	else {
	assert(!hard);
	real = _decayme->generateHardest( currentTree()->perturbativeProcess() );
	type = _decayme->hasPOWHEGCorrection();
	}
	if(real) {
	// set up ther hard tree
	if(!real->outgoing().empty()) _hardtree = new_ptr(HardTree(real));
	// set up the vetos
	currentTree()->setVetoes(real->pT(),type);
	}
	// store initial state POWHEG radiation
	if(_hardtree && _hardme && _hardme->hasPOWHEGCorrection()==1)
	ISRTree = _hardtree;
	}
	else if (hard) {
	// Get minimum pT cutoff used in shower approximation
	Energy maxpt = 1.*GeV;

	if ( currentTree()->showerApproximation() ) {
	int colouredIn = 0;
	int colouredOut = 0;
	for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
	= currentTree()->outgoingLines().begin();
	it != currentTree()->outgoingLines().end(); ++it ) {
	if( it->second->coloured() ) ++colouredOut;
	}
	for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
	= currentTree()->incomingLines().begin();
	it != currentTree()->incomingLines().end(); ++it ) {
	if( it->second->coloured() ) ++colouredIn;
	}
	if ( currentTree()->showerApproximation()->ffPtCut() == currentTree()->showerApproximation()->fiPtCut() &&
	currentTree()->showerApproximation()->ffPtCut() == currentTree()->showerApproximation()->iiPtCut() )
	maxpt = currentTree()->showerApproximation()->ffPtCut();
	else if ( colouredIn == 2 && colouredOut == 0 )
	maxpt = currentTree()->showerApproximation()->iiPtCut();
	else if ( colouredIn == 0 && colouredOut > 1 )
	maxpt = currentTree()->showerApproximation()->ffPtCut();
	else if ( colouredIn == 2 && colouredOut == 1 )
	maxpt = min(currentTree()->showerApproximation()->iiPtCut(), currentTree()->showerApproximation()->fiPtCut());
	else if ( colouredIn == 1 && colouredOut > 1 )
	maxpt = min(currentTree()->showerApproximation()->ffPtCut(), currentTree()->showerApproximation()->fiPtCut());
	else
	maxpt = min(min(currentTree()->showerApproximation()->iiPtCut(), currentTree()->showerApproximation()->fiPtCut()),
	currentTree()->showerApproximation()->ffPtCut());
	}

	// Generate hardtree from born and real emission subprocesses
	_hardtree = generateCKKW(currentTree());

	// Find transverse momentum of hardest emission
	if (_hardtree){
	for(set<HardBranchingPtr>::iterator it=_hardtree->branchings().begin();
	it!=_hardtree->branchings().end();++it) {
	if ((it)->parent() && (it)->status()==HardBranching::Incoming)
	maxpt=(*it)->branchingParticle()->momentum().perp();
	if ((it)->children().size()==2 && (it)->status()==HardBranching::Outgoing){
	if ((*it)->branchingParticle()->id()!=21 &&
	abs((*it)->branchingParticle()->id())>5 ){
	if ((*it)->children()[0]->branchingParticle()->id()==21 \|\|
	abs((*it)->children()[0]->branchingParticle()->id())<6)
	maxpt=(*it)->children()[0]->branchingParticle()->momentum().perp();
	else if ((*it)->children()[1]->branchingParticle()->id()==21 \|\|
	abs((*it)->children()[1]->branchingParticle()->id())<6)
	maxpt=(*it)->children()[1]->branchingParticle()->momentum().perp();
	}
	else {
	if ( abs((*it)->branchingParticle()->id())<6){
	if (abs((*it)->children()[0]->branchingParticle()->id())<6)
	maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp();
	else
	maxpt = (*it)->children()[0]->branchingParticle()->momentum().perp();
	}
	else maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp();
	}
	}
	}
	}


	// Hardest (pt) emission should be the first powheg emission.
	maxpt=min(sqrt(lastXCombPtr()->lastShowerScale()),maxpt);

	// set maximum pT for subsequent emissions from S events
	if ( currentTree()->isPowhegSEvent() ) {
	for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
	= currentTree()->outgoingLines().begin();
	it != currentTree()->outgoingLines().end(); ++it ) {
	if( ! it->second->coloured() ) continue;
	it->first->maximumpT(maxpt, ShowerInteraction::QCD );
	}
	for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
	= currentTree()->incomingLines().begin();
	it != currentTree()->incomingLines().end(); ++it ) {
	if( ! it->second->coloured() ) continue;
	it->first->maximumpT(maxpt, ShowerInteraction::QCD );
	}
	}
	}
	else
	_hardtree = generateCKKW(currentTree());

	// if hard me doesn't have a FSR powheg
	// correction use decay powheg correction
	if (_hardme && _hardme->hasPOWHEGCorrection()<2) {
	addFSRUsingDecayPOWHEG(ISRTree);
	}
	// connect the trees
	if(_hardtree) {
	connectTrees(currentTree(),_hardtree,hard);
	}
	}

	void QTildeShowerHandler::addFSRUsingDecayPOWHEG(HardTreePtr ISRTree) {
	// check for intermediate colour singlet resonance
	const ParticleVector inter = _hardme->subProcess()->intermediates();
	if (inter.size()!=1 \|\| inter[0]->momentum().m2()/GeV2 < 0 \|\|
	inter[0]->dataPtr()->iColour()!=PDT::Colour0) {
	return;
	}

	// ignore cases where outgoing particles are not coloured
	map<ShowerProgenitorPtr, tShowerParticlePtr > out = currentTree()->outgoingLines();
	if (out.size() != 2 \|\|
	out. begin()->second->dataPtr()->iColour()==PDT::Colour0 \|\|
	out.rbegin()->second->dataPtr()->iColour()==PDT::Colour0) {
	return;
	}

	// look up decay mode
	tDMPtr dm;
	string tag;
	string inParticle = inter[0]->dataPtr()->name() + "->";
	vector<string> outParticles;
	outParticles.push_back(out.begin ()->first->progenitor()->dataPtr()->name());
	outParticles.push_back(out.rbegin()->first->progenitor()->dataPtr()->name());
	for (int it=0; it<2; ++it){
	tag = inParticle + outParticles[it] + "," + outParticles[(it+1)%2] + ";";
	dm = generator()->findDecayMode(tag);
	if(dm) break;
	}

	// get the decayer
	HwDecayerBasePtr decayer;
	if(dm) decayer = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
	// check if decayer has a FSR POWHEG correction
	if (!decayer \|\| decayer->hasPOWHEGCorrection()<2) {
	return;
	}
	// generate the hardest emission
	// create RealEmissionProcess
	PPtr in = new_ptr(*inter[0]);
	RealEmissionProcessPtr newProcess(new_ptr(RealEmissionProcess()));
	newProcess->bornIncoming().push_back(in);
	newProcess->bornOutgoing().push_back(out.begin ()->first->progenitor());
	newProcess->bornOutgoing().push_back(out.rbegin()->first->progenitor());
	// generate the FSR
	newProcess = decayer->generateHardest(newProcess);
	HardTreePtr FSRTree;
	if(newProcess) {
	// set up ther hard tree
	if(!newProcess->outgoing().empty()) FSRTree = new_ptr(HardTree(newProcess));
	// set up the vetos
	currentTree()->setVetoes(newProcess->pT(),2);
	}
	if(!FSRTree) return;

	// if there is no ISRTree make _hardtree from FSRTree
	if (!ISRTree){
	vector<HardBranchingPtr> inBranch,hardBranch;
	for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit =currentTree()->incomingLines().begin();
	cit!=currentTree()->incomingLines().end();++cit ) {
	inBranch.push_back(new_ptr(HardBranching(cit->second,SudakovPtr(),
	HardBranchingPtr(),
	HardBranching::Incoming)));
	inBranch.back()->beam(cit->first->original()->parents()[0]);
	hardBranch.push_back(inBranch.back());
	}
	if(inBranch[0]->branchingParticle()->dataPtr()->coloured()) {
	inBranch[0]->colourPartner(inBranch[1]);
	inBranch[1]->colourPartner(inBranch[0]);
	}
	for(set<HardBranchingPtr>::iterator it=FSRTree->branchings().begin();
	it!=FSRTree->branchings().end();++it) {
	if((**it).branchingParticle()->id()!=in->id())
	hardBranch.push_back(*it);
	}
	hardBranch[2]->colourPartner(hardBranch[3]);
	hardBranch[3]->colourPartner(hardBranch[2]);
	HardTreePtr newTree = new_ptr(HardTree(hardBranch,inBranch,
	ShowerInteraction::QCD));
	_hardtree = newTree;
	}

	// Otherwise modify the ISRTree to include the emission in FSRTree
	else {
	vector<tShowerParticlePtr> FSROut, ISROut;
	set<HardBranchingPtr>::iterator itFSR, itISR;
	// get outgoing particles
	for(itFSR =FSRTree->branchings().begin();
	itFSR!=FSRTree->branchings().end();++itFSR){
	if ((**itFSR).status()==HardBranching::Outgoing)
	FSROut.push_back((*itFSR)->branchingParticle());
	}
	for(itISR =ISRTree->branchings().begin();
	itISR!=ISRTree->branchings().end();++itISR){
	if ((**itISR).status()==HardBranching::Outgoing)
	ISROut.push_back((*itISR)->branchingParticle());
	}

	// find COM frame formed by outgoing particles
	LorentzRotation eventFrameFSR, eventFrameISR;
	eventFrameFSR = ((FSROut[0]->momentum()+FSROut[1]->momentum()).findBoostToCM());
	eventFrameISR = ((ISROut[0]->momentum()+ISROut[1]->momentum()).findBoostToCM());

	// find rotation between ISR and FSR frames
	int j=0;
	if (ISROut[0]->id()!=FSROut[0]->id()) j=1;
	eventFrameISR.rotateZ( (eventFrameFSR*FSROut[0]->momentum()).phi()-
	(eventFrameISR*ISROut[j]->momentum()).phi() );
	eventFrameISR.rotateY( (eventFrameFSR*FSROut[0]->momentum()).theta()-
	(eventFrameISR*ISROut[j]->momentum()).theta() );
	eventFrameISR.invert();

	for (itFSR=FSRTree->branchings().begin();
	itFSR!=FSRTree->branchings().end();++itFSR){
	if ((**itFSR).branchingParticle()->id()==in->id()) continue;
	for (itISR =ISRTree->branchings().begin();
	itISR!=ISRTree->branchings().end();++itISR){
	if ((**itISR).status()==HardBranching::Incoming) continue;
	if ((**itFSR).branchingParticle()->id()==
	(**itISR).branchingParticle()->id()){
	// rotate FSRTree particle to ISRTree event frame
	(*itISR).branchingParticle()->setMomentum(eventFrameISR
	eventFrameFSR*
	(**itFSR).branchingParticle()->momentum());
	(**itISR).branchingParticle()->rescaleMass();
	// add the children of the FSRTree particles to the ISRTree
	if(!(**itFSR).children().empty()){
	(itISR).addChild((itFSR).children()[0]);
	(itISR).addChild((itFSR).children()[1]);
	// rotate momenta to ISRTree event frame
	(*itISR).children()[0]->branchingParticle()->setMomentum(eventFrameISR
	eventFrameFSR*
	(**itFSR).children()[0]->branchingParticle()->momentum());
	(*itISR).children()[1]->branchingParticle()->setMomentum(eventFrameISR
	eventFrameFSR*
	(**itFSR).children()[1]->branchingParticle()->momentum());
	}
	}
	}
	}
	_hardtree = ISRTree;
	}
	}

	bool QTildeShowerHandler::truncatedTimeLikeShower(tShowerParticlePtr particle,
	HardBranchingPtr branch,
	ShowerInteraction type,
	Branching fb, bool first) {
	// select a branching if we don't have one
	if(!fb.kinematics)
	fb = selectTimeLikeBranching(particle,type,branch);
	// must be an emission, the forced one it not a truncated one
	assert(fb.kinematics);
	ShowerParticleVector children;
	Branching fc[2] = {Branching(),Branching()};
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the children
	children = createTimeLikeChildren(particle,fb.ids);
	// update the children
	particle->showerKinematics()->
	updateChildren(particle, children,fb.type);
	// select branchings for children
	if(!fc[0].kinematics) {
	// select branching for first particle
	if(!fb.hard && fb.iout ==1 )
	fc[0] = selectTimeLikeBranching(children[0],type,branch);
	else if(fb.hard && !branch->children()[0]->children().empty() )
	fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]);
	else
	fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
	}
	// select branching for the second particle
	if(!fc[1].kinematics) {
	// select branching for first particle
	if(!fb.hard && fb.iout ==2 )
	fc[1] = selectTimeLikeBranching(children[1],type,branch);
	else if(fb.hard && !branch->children()[1]->children().empty() )
	fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]);
	else
	fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
	}
	// shower the first particle
	if(fc[0].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 1)
	truncatedTimeLikeShower(children[0],branch,type,fc[0],false);
	// hard emission and subsquent hard emissions
	else if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 2)
	truncatedTimeLikeShower(children[1],branch,type,fc[1],false);
	// hard emission and subsquent hard emissions
	else if(fb.hard && !branch->children()[1]->children().empty() )
	truncatedTimeLikeShower(children[1],branch->children()[1],type,fc[1],false);
	else
	timeLikeShower(children[1],type,fc[1],false);
	}
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	if(first&&!children.empty())
	particle->showerKinematics()->resetChildren(particle,children);
	if(particle->spinInfo()) particle->spinInfo()->develop();
	// TODO NEED AN UPDATE HERE FOR RECONOPT!=0 ?
	return true;
	}

	bool QTildeShowerHandler::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
	HardBranchingPtr branch,
	ShowerInteraction type) {
	tcPDFPtr pdf;
	if(firstPDF().particle() == beamParticle())
	pdf = firstPDF().pdf();
	if(secondPDF().particle() == beamParticle())
	pdf = secondPDF().pdf();
	Energy freeze = pdfFreezingScale();
	Branching bb;
	// parameters of the force branching
	double z(0.);
	HardBranchingPtr timelike;
	for( unsigned int ix = 0; ix < branch->children().size(); ++ix ) {
	if( branch->children()[ix]->status() ==HardBranching::Outgoing) {
	timelike = branch->children()[ix];
	}
	if( branch->children()[ix]->status() ==HardBranching::Incoming )
	z = branch->children()[ix]->z();
	}
	// generate truncated branching
	tcPDPtr part[2];
	if(z>=0.&&z<=1.) {
	while (true) {
	if( !isTruncatedShowerON() \|\| hardOnly() ) break;
	bb = splittingGenerator()->chooseBackwardBranching( *particle,
	beam, 1., beamParticle(),
	type , pdf,freeze);
	if( !bb.kinematics \|\| bb.kinematics->scale() < branch->scale() ) {
	bb = Branching();
	break;
	}
	// particles as in Sudakov form factor
	part[0] = bb.ids[0];
	part[1] = bb.ids[2];
	double zsplit = bb.kinematics->z();
	// apply the vetos for the truncated shower
	// if doesn't carry most of momentum
	ShowerInteraction type2 = convertInteraction(bb.type);
	if(type2==branch->sudakov()->interactionType() &&
	zsplit < 0.5) {
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	continue;
	}
	// others
	if( part[0]->id() != particle->id() \|\| // if particle changes type
	bb.kinematics->pT() > progenitor()->maximumpT(type2) \|\| // pt veto
	bb.kinematics->scale() < branch->scale()) { // angular ordering veto
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	continue;
	}
	// and those from the base class
	if(spaceLikeVetoed(bb,particle)) {
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	continue;
	}
	break;
	}
	}
	if( !bb.kinematics ) {
	//do the hard emission
	ShoKinPtr kinematics = new_ptr(IS_QTildeShowerKinematics1to2(
	branch->scale(), z, branch->phi(),
	branch->children()[0]->pT(), branch->sudakov()
	));
	// assign the splitting function and shower kinematics
	particle->showerKinematics( kinematics );
	if(kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent =
	new_ptr( ShowerParticle( branch->branchingParticle()->dataPtr(), false ) );
	ShowerParticlePtr otherChild =
	new_ptr( ShowerParticle( timelike->branchingParticle()->dataPtr(),
	true, true ) );
	ShowerParticleVector theChildren;
	theChildren.push_back( particle );
	theChildren.push_back( otherChild );
	particle->showerKinematics()->
	updateParent( newParent, theChildren, branch->type());
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
	currentTree()->addInitialStateBranching( particle, newParent, otherChild );
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	bool emitted=false;
	if(!hardOnly()) {
	if( branch->parent() ) {
	emitted = truncatedSpaceLikeShower( newParent, beam, branch->parent() , type);
	}
	else {
	emitted = spaceLikeShower( newParent, beam , type);
	}
	}
	if( !emitted ) {
	if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
	kinematics->updateLast( newParent, ZERO, ZERO );
	}
	else {
	pair<Energy,double> kt = intrinsicpT()[progenitor()];
	kinematics->updateLast( newParent,
	kt.first*cos( kt.second ),
	kt.first*sin( kt.second ) );
	}
	}
	particle->showerKinematics()->
	updateChildren( newParent, theChildren,bb.type);
	if(hardOnly()) return true;
	// perform the shower of the final-state particle
	if( timelike->children().empty() ) {
	timeLikeShower( otherChild , type,Branching(),true);
	}
	else {
	truncatedTimeLikeShower( otherChild, timelike , type,Branching(), true);
	}
	updateHistory(otherChild);
	// return the emitted
	return true;
	}
	// assign the splitting function and shower kinematics
	particle->showerKinematics( bb.kinematics );
	if(bb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(bb.kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent = new_ptr( ShowerParticle( part[0], false ) );
	ShowerParticlePtr otherChild = new_ptr( ShowerParticle( part[1], true, true ) );
	ShowerParticleVector theChildren;
	theChildren.push_back( particle );
	theChildren.push_back( otherChild );
	particle->showerKinematics()->
	updateParent( newParent, theChildren, bb.type);
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
	currentTree()->addInitialStateBranching( particle, newParent, otherChild );
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	bool emitted = truncatedSpaceLikeShower( newParent, beam, branch,type);
	// now reconstruct the momentum
	if( !emitted ) {
	if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
	bb.kinematics->updateLast( newParent, ZERO, ZERO );
	}
	else {
	pair<Energy,double> kt = intrinsicpT()[ progenitor() ];
	bb.kinematics->updateLast( newParent,
	kt.first*cos( kt.second ),
	kt.first*sin( kt.second ) );
	}
	}
	particle->showerKinematics()->
	updateChildren( newParent, theChildren, bb.type);
	// perform the shower of the final-state particle
	timeLikeShower( otherChild , type,Branching(),true);
	updateHistory(otherChild);
	// return the emitted
	return true;
	}

	bool QTildeShowerHandler::
	truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass, HardBranchingPtr branch,
	ShowerInteraction type, Branching fb) {
	// select a branching if we don't have one
	if(!fb.kinematics)
	fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,branch);
	// must be an emission, the forced one it not a truncated one
	assert(fb.kinematics);
	ShowerParticleVector children;
	Branching fc[2]={Branching(),Branching()};
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the ShowerParticle objects for the two children
	children = createTimeLikeChildren(particle,fb.ids);
	// updateChildren the children
	particle->showerKinematics()->
	updateChildren(particle, children, fb.type);
	// select branchings for children
	if(!fc[0].kinematics) {
	if(children[0]->id()==particle->id()) {
	// select branching for first particle
	if(!fb.hard)
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,branch);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,
	branch->children()[0]);
	else
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,
	HardBranchingPtr());
	}
	else {
	// select branching for first particle
	if(fb.hard && !branch->children()[0]->children().empty() )
	fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]);
	else
	fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
	}
	}
	// select branching for the second particle
	if(!fc[1].kinematics) {
	if(children[1]->id()==particle->id()) {
	// select branching for first particle
	if(!fb.hard)
	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,branch);
	else if(fb.hard && ! branch->children()[1]->children().empty() )
	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,
	branch->children()[1]);
	else
	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,
	HardBranchingPtr());
	}
	else {
	if(fb.hard && !branch->children()[1]->children().empty() )
	fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]);
	else
	fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
	}
	}
	// old default
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),children[0]);
	currentTree()->addInitialStateBranching(particle,children[0],children[1]);
	// shower the first particle
	if(fc[0].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[0]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[0]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[0]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	}
	// shower the second particle
	if(fc[1].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[1]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[1]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[1]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[1],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[1],false);
	}
	}
	updateHistory(children[1]);
	// TODO DO WE NEED A CHECK IF RECONPT !=0
	return true;
	}

	void QTildeShowerHandler::connectTrees(ShowerTreePtr showerTree,
	HardTreePtr hardTree, bool hard ) {
	ShowerParticleVector particles;
	// find the Sudakovs
	for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
	cit!=hardTree->branchings().end();++cit) {
	// Sudakovs for ISR
	if((cit).parent()&&(cit).status()==HardBranching::Incoming) {
	++_nis;
	array<long,3> br;
	br[0] = (**cit).parent()->branchingParticle()->id();
	br[1] = (**cit). branchingParticle()->id();
	br[2] = (*cit).parent()->children()[0]==cit ?
	(**cit).parent()->children()[1]->branchingParticle()->id() :
	(**cit).parent()->children()[0]->branchingParticle()->id();
	BranchingList branchings = splittingGenerator()->initialStateBranchings();
	if(br[1]<0&&br[0]==br[1]) {
	br[0] = abs(br[0]);
	br[1] = abs(br[1]);
	}
	else if(br[1]<0) {
	br[1] = -br[1];
	br[2] = -br[2];
	}
	long index = abs(br[1]);
	SudakovPtr sudakov;
	for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
	cjt != branchings.upper_bound(index); ++cjt ) {
	IdList ids = cjt->second.particles;
	if(ids[0]->id()==br[0]&&ids[1]->id()==br[1]&&ids[2]->id()==br[2]) {
	sudakov=cjt->second.sudakov;
	break;
	}
	}
	if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
	<< "QTildeShowerHandler::connectTrees() for ISR"
	<< Exception::runerror;
	(**cit).parent()->sudakov(sudakov);
	}
	// Sudakovs for FSR
	else if(!(**cit).children().empty()) {
	++_nfs;
	array<long,3> br;
	br[0] = (**cit) .branchingParticle()->id();
	br[1] = (**cit).children()[0]->branchingParticle()->id();
	br[2] = (**cit).children()[1]->branchingParticle()->id();
	BranchingList branchings = splittingGenerator()->finalStateBranchings();
	if(br[0]<0) {
	br[0] = abs(br[0]);
	br[1] = abs(br[1]);
	br[2] = abs(br[2]);
	}
	long index = br[0];
	SudakovPtr sudakov;
	for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
	cjt != branchings.upper_bound(index); ++cjt ) {
	IdList ids = cjt->second.particles;
	if(ids[0]->id()==br[0]&&ids[1]->id()==br[1]&&ids[2]->id()==br[2]) {
	sudakov=cjt->second.sudakov;
	break;
	}
	}
	if(!sudakov) {
	throw Exception() << "Can't find Sudakov for the hard emission in "
	<< "QTildeShowerHandler::connectTrees()"
	<< Exception::runerror;
	}
	(**cit).sudakov(sudakov);
	}
	}
	// calculate the evolution scale
	for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
	cit!=hardTree->branchings().end();++cit) {
	particles.push_back((*cit)->branchingParticle());
	}
	partnerFinder()->
	setInitialEvolutionScales(particles,!hard,interaction_,true);
	hardTree->partnersSet(true);
	// inverse reconstruction
	if(hard) {
	kinematicsReconstructor()->
	deconstructHardJets(hardTree,interaction_);
	}
	else
	kinematicsReconstructor()->
	deconstructDecayJets(hardTree,interaction_);
	// now reset the momenta of the showering particles
	vector<ShowerProgenitorPtr> particlesToShower=showerTree->extractProgenitors();
	// match them
	map<ShowerProgenitorPtr,HardBranchingPtr> partners;
	for(set<HardBranchingPtr>::const_iterator bit=hardTree->branchings().begin();
	bit!=hardTree->branchings().end();++bit) {
	Energy2 dmin( 1e30*GeV2 );
	ShowerProgenitorPtr partner;
	for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
	pit!=particlesToShower.end();++pit) {
	if(partners.find(*pit)!=partners.end()) continue;
	if( (bit).branchingParticle()->id() != (pit).progenitor()->id() ) continue;
	if( (**bit).branchingParticle()->isFinalState() !=
	(**pit).progenitor()->isFinalState() ) continue;
	if( (**pit).progenitor()->isFinalState() ) {
	Energy2 dtest =
	sqr( (pit).progenitor()->momentum().x() - (bit).showerMomentum().x() ) +
	sqr( (pit).progenitor()->momentum().y() - (bit).showerMomentum().y() ) +
	sqr( (pit).progenitor()->momentum().z() - (bit).showerMomentum().z() ) +
	sqr( (pit).progenitor()->momentum().t() - (bit).showerMomentum().t() );
	// add mass difference for identical particles (e.g. Z0 Z0 production)
	dtest += 1e10sqr((pit).progenitor()->momentum().m()-(*bit).showerMomentum().m());
	if( dtest < dmin ) {
	partner = *pit;
	dmin = dtest;
	}
	}
	else {
	// ensure directions are right
	if((pit).progenitor()->momentum().z()/(bit).showerMomentum().z()>ZERO) {
	partner = *pit;
	break;
	}
	}
	}
	if(!partner) throw Exception() << "Failed to match shower and hard trees in QTildeShowerHandler::hardestEmission"
	<< Exception::eventerror;
	partners[partner] = *bit;
	}
	for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
	pit!=particlesToShower.end();++pit) {
	HardBranchingPtr partner = partners[*pit];
	if((**pit).progenitor()->dataPtr()->stable()) {
	(**pit).progenitor()->set5Momentum(partner->showerMomentum());
	(**pit).copy()->set5Momentum(partner->showerMomentum());
	}
	else {
	Lorentz5Momentum oldMomentum = (**pit).progenitor()->momentum();
	Lorentz5Momentum newMomentum = partner->showerMomentum();
	LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
	(**pit).progenitor()->transform(boost);
	(**pit).copy() ->transform(boost);
	boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
	(**pit).progenitor()->transform(boost);
	(**pit).copy() ->transform(boost);
	}
	}
	// correction boosts for daughter trees
	for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	tit = showerTree->treelinks().begin();
	tit != showerTree->treelinks().end();++tit) {
	ShowerTreePtr decayTree = tit->first;
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit = decayTree->incomingLines().begin();
	// reset the momentum of the decay particle
	Lorentz5Momentum oldMomentum = cit->first->progenitor()->momentum();
	Lorentz5Momentum newMomentum = tit->second.second->momentum();
	LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
	decayTree->transform(boost,true);
	boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
	decayTree->transform(boost,true);
	}
	}

	void QTildeShowerHandler::doShowering(bool hard,XCPtr xcomb) {
	// zero number of emissions
	_nis = _nfs = 0;
	// if MC@NLO H event and limited emissions
	// indicate both final and initial state emission
	if ( currentTree()->isMCatNLOHEvent() && _limitEmissions != 0 ) {
	_nis = _nfs = 1;
	}
	// extract particles to shower
	vector<ShowerProgenitorPtr> particlesToShower(setupShower(hard));
	// check if we should shower
	bool colCharge = false;
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(particlesToShower[ix]->progenitor()->dataPtr()->coloured() \|\|
	particlesToShower[ix]->progenitor()->dataPtr()->charged()) {
	colCharge = true;
	break;
	}
	}
	if(!colCharge) {
	_currenttree->hasShowered(true);
	return;
	}
	// setup the maximum scales for the shower
	if (restrictPhasespace()) setupMaximumScales(particlesToShower,xcomb);
	// set the hard scales for the profiles
	setupHardScales(particlesToShower,xcomb);
	// specific stuff for hard processes and decays
	Energy minmass(ZERO), mIn(ZERO);
	// hard process generate the intrinsic p_T once and for all
	if(hard) {
	generateIntrinsicpT(particlesToShower);
	}
	// decay compute the minimum mass of the final-state
	else {
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(particlesToShower[ix]->progenitor()->isFinalState()) {
	- if(particlesToShower[ix]->progenitor()->dataPtr()->stable()){
	- auto dm= ShowerHandler::currentHandler()->retConstituentMasses()?
	- particlesToShower[ix]->progenitor()->dataPtr()->constituentMass():
	- particlesToShower[ix]->progenitor()->dataPtr()->mass();
	- minmass += dm;
	- }else
	+ if(particlesToShower[ix]->progenitor()->dataPtr()->stable())
	+ minmass += particlesToShower[ix]->progenitor()->dataPtr()->constituentMass();
	+ else
	minmass += particlesToShower[ix]->progenitor()->mass();
	}
	else {
	mIn = particlesToShower[ix]->progenitor()->mass();
	}
	}
	// throw exception if decay can't happen
	if ( minmass > mIn ) {
	throw Exception() << "QTildeShowerHandler.cc: Mass of decaying particle is "
	<< "below constituent masses of decay products."
	<< Exception::eventerror;
	}
	}
	// setup for reweighted
	bool reWeighting = _reWeight && hard && ShowerHandler::currentHandler()->firstInteraction();
	double eventWeight=0.;
	unsigned int nTryReWeight(0);
	// create random particle vector (only need to do once)
	vector<ShowerProgenitorPtr> tmp;
	unsigned int nColouredIncoming = 0;
	while(particlesToShower.size()>0){
	unsigned int xx=UseRandom::irnd(particlesToShower.size());
	tmp.push_back(particlesToShower[xx]);
	particlesToShower.erase(particlesToShower.begin()+xx);
	}
	particlesToShower=tmp;
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(!particlesToShower[ix]->progenitor()->isFinalState() &&
	particlesToShower[ix]->progenitor()->coloured()) ++nColouredIncoming;
	}
	bool switchRecon = hard && nColouredIncoming !=1;
	// main shower loop
	unsigned int ntry(0);
	bool reconstructed = false;
	do {
	// clear results of last attempt if needed
	if(ntry!=0) {
	currentTree()->clear();
	setEvolutionPartners(hard,interaction_,true);
	_nis = _nfs = 0;
	// if MC@NLO H event and limited emissions
	// indicate both final and initial state emission
	if ( currentTree()->isMCatNLOHEvent() && _limitEmissions != 0 ) {
	_nis = _nfs = 1;
	}
	for(unsigned int ix=0; ix<particlesToShower.size();++ix) {
	SpinPtr spin = particlesToShower[ix]->progenitor()->spinInfo();
	if(spin && spin->decayVertex() &&
	dynamic_ptr_cast<tcSVertexPtr>(spin->decayVertex())) {
	spin->decayVertex(VertexPtr());
	}
	}
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(particlesToShower[ix]->progenitor()->isFinalState() \|\|
	(hard && !particlesToShower[ix]->progenitor()->isFinalState())) {
	if(particlesToShower[ix]->progenitor()->spinInfo())
	particlesToShower[ix]->progenitor()->spinInfo()->reset();
	}
	}
	}
	// loop over particles
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	// extract the progenitor
	progenitor(particlesToShower[ix]);
	// final-state radiation
	if(progenitor()->progenitor()->isFinalState()) {
	if(!doFSR()) continue;
	// perform shower
	progenitor()->hasEmitted(startTimeLikeShower(interaction_));
	}
	// initial-state radiation
	else {
	if(!doISR()) continue;
	// hard process
	if(hard) {
	// get the PDF
	setBeamParticle(_progenitor->beam());
	if(!beamParticle()) {
	throw Exception() << "Incorrect type of beam particle in "
	<< "QTildeShowerHandler::doShowering(). "
	<< "This should not happen for conventional choices but may happen if you have used a"
	<< " non-default choice and have not changed the create ParticleData line in the input files"
	<< " for this particle to create BeamParticleData."
	<< Exception::runerror;
	}
	// perform the shower
	// set the beam particle
	tPPtr beamparticle=progenitor()->original();
	if(!beamparticle->parents().empty())
	beamparticle=beamparticle->parents()[0];
	// generate the shower
	progenitor()->hasEmitted(startSpaceLikeShower(beamparticle,
	interaction_));
	}
	// decay
	else {
	// skip colour and electrically neutral particles
	if(!progenitor()->progenitor()->dataPtr()->coloured() &&
	!progenitor()->progenitor()->dataPtr()->charged()) {
	progenitor()->hasEmitted(false);
	continue;
	}
	// perform shower
	// set the scales correctly. The current scale is the maximum scale for
	// emission not the starting scale
	ShowerParticle::EvolutionScales maxScales(progenitor()->progenitor()->scales());
	progenitor()->progenitor()->scales() = ShowerParticle::EvolutionScales();
	if(progenitor()->progenitor()->dataPtr()->charged()) {
	progenitor()->progenitor()->scales().QED = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QED_noAO = progenitor()->progenitor()->mass();
	}
	if(progenitor()->progenitor()->hasColour()) {
	progenitor()->progenitor()->scales().QCD_c = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QCD_c_noAO = progenitor()->progenitor()->mass();
	}
	if(progenitor()->progenitor()->hasAntiColour()) {
	progenitor()->progenitor()->scales().QCD_ac = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QCD_ac_noAO = progenitor()->progenitor()->mass();
	}
	// perform the shower
	progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass,
	interaction_));
	}
	}
	}
	// do the kinematic reconstruction, checking if it worked
	reconstructed = hard ?
	kinematicsReconstructor()->
	reconstructHardJets (currentTree(),intrinsicpT(),interaction_,
	switchRecon && ntry>maximumTries()/2) :
	kinematicsReconstructor()->
	reconstructDecayJets(currentTree(),interaction_);
	if(!reconstructed) continue;
	// apply vetos on the full shower
	for(vector<FullShowerVetoPtr>::const_iterator it=_fullShowerVetoes.begin();
	it!=_fullShowerVetoes.end();++it) {
	int veto = (**it).applyVeto(currentTree());
	if(veto<0) continue;
	// veto the shower
	if(veto==0) {
	reconstructed = false;
	break;
	}
	// veto the shower and reweight
	else if(veto==1) {
	reconstructed = false;
	break;
	}
	// veto the event
	else if(veto==2) {
	throw Veto();
	}
	}
	if(reWeighting) {
	if(reconstructed) eventWeight += 1.;
	reconstructed=false;
	++nTryReWeight;
	if(nTryReWeight==_nReWeight) {
	reWeighting = false;
	if(eventWeight==0.) throw Veto();
	}
	}
	}
	while(!reconstructed&&maximumTries()>++ntry);
	// check if failed to generate the shower
	if(ntry==maximumTries()) {
	if(hard)
	throw ShowerHandler::ShowerTriesVeto(ntry);
	else
	throw Exception() << "Failed to generate the shower after "
	<< ntry << " attempts in QTildeShowerHandler::showerDecay()"
	<< Exception::eventerror;
	}
	// handle the weights and apply any reweighting required
	if(nTryReWeight>0) {
	tStdEHPtr seh = dynamic_ptr_cast<tStdEHPtr>(generator()->currentEventHandler());
	static bool first = true;
	if(seh) {
	seh->reweight(eventWeight/double(nTryReWeight));
	}
	else if(first) {
	generator()->log() << "Reweighting the shower only works with internal Herwig7 processes"
	<< "Presumably you are showering Les Houches Events. These will not be"
	<< "reweighted\n";
	first = false;
	}
	}
	// tree has now showered
	_currenttree->hasShowered(true);
	hardTree(HardTreePtr());
	}

	void QTildeShowerHandler:: convertHardTree(bool hard,ShowerInteraction type) {
	map<ColinePtr,ColinePtr> cmap;
	// incoming particles
	for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit=currentTree()->incomingLines().begin();cit!=currentTree()->incomingLines().end();++cit) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	mit = hardTree()->particles().find(cit->first->progenitor());
	// put the colour lines in the map
	ShowerParticlePtr oldParticle = cit->first->progenitor();
	ShowerParticlePtr newParticle = mit->second->branchingParticle();
	ColinePtr cLine = oldParticle-> colourLine();
	ColinePtr aLine = oldParticle->antiColourLine();
	if(newParticle->colourLine() &&
	cmap.find(newParticle-> colourLine())==cmap.end())
	cmap[newParticle-> colourLine()] = cLine;
	if(newParticle->antiColourLine() &&
	cmap.find(newParticle->antiColourLine())==cmap.end())
	cmap[newParticle->antiColourLine()] = aLine;
	// check whether or not particle emits
	bool emission = mit->second->parent();
	if(emission) {
	if(newParticle->colourLine()) {
	ColinePtr ctemp = newParticle-> colourLine();
	ctemp->removeColoured(newParticle);
	}
	if(newParticle->antiColourLine()) {
	ColinePtr ctemp = newParticle->antiColourLine();
	ctemp->removeAntiColoured(newParticle);
	}
	newParticle = mit->second->parent()->branchingParticle();
	}
	// get the new colour lines
	ColinePtr newCLine,newALine;
	// sort out colour lines
	if(newParticle->colourLine()) {
	ColinePtr ctemp = newParticle-> colourLine();
	ctemp->removeColoured(newParticle);
	if(cmap.find(ctemp)!=cmap.end()) {
	newCLine = cmap[ctemp];
	}
	else {
	newCLine = new_ptr(ColourLine());
	cmap[ctemp] = newCLine;
	}
	}
	// and anticolour lines
	if(newParticle->antiColourLine()) {
	ColinePtr ctemp = newParticle->antiColourLine();
	ctemp->removeAntiColoured(newParticle);
	if(cmap.find(ctemp)!=cmap.end()) {
	newALine = cmap[ctemp];
	}
	else {
	newALine = new_ptr(ColourLine());
	cmap[ctemp] = newALine;
	}
	}
	// remove colour lines from old particle
	if(aLine) {
	aLine->removeAntiColoured(cit->first->copy());
	aLine->removeAntiColoured(cit->first->progenitor());
	}
	if(cLine) {
	cLine->removeColoured(cit->first->copy());
	cLine->removeColoured(cit->first->progenitor());
	}
	// add particle to colour lines
	if(newCLine) newCLine->addColoured (newParticle);
	if(newALine) newALine->addAntiColoured(newParticle);
	// insert new particles
	cit->first->copy(newParticle);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,false)));
	cit->first->progenitor(sp);
	currentTree()->incomingLines()[cit->first]=sp;
	cit->first->perturbative(!emission);
	// and the emitted particle if needed
	if(emission) {
	ShowerParticlePtr newOut = mit->second->parent()->children()[1]->branchingParticle();
	if(newOut->colourLine()) {
	ColinePtr ctemp = newOut-> colourLine();
	ctemp->removeColoured(newOut);
	assert(cmap.find(ctemp)!=cmap.end());
	cmap[ctemp]->addColoured (newOut);
	}
	if(newOut->antiColourLine()) {
	ColinePtr ctemp = newOut->antiColourLine();
	ctemp->removeAntiColoured(newOut);
	assert(cmap.find(ctemp)!=cmap.end());
	cmap[ctemp]->addAntiColoured(newOut);
	}
	ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true));
	ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout));
	out->perturbative(false);
	currentTree()->outgoingLines().insert(make_pair(out,sout));
	}
	if(hard) {
	// sort out the value of x
	if(mit->second->beam()->momentum().z()>ZERO) {
	sp->x(newParticle->momentum(). plus()/mit->second->beam()->momentum(). plus());
	}
	else {
	sp->x(newParticle->momentum().minus()/mit->second->beam()->momentum().minus());
	}
	}
	}
	// outgoing particles
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cit=currentTree()->outgoingLines().begin();cit!=currentTree()->outgoingLines().end();++cit) {
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,
	tShowerParticlePtr> >::const_iterator tit;
	for(tit = currentTree()->treelinks().begin();
	tit != currentTree()->treelinks().end();++tit) {
	if(tit->second.first && tit->second.second==cit->first->progenitor())
	break;
	}
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	mit = hardTree()->particles().find(cit->first->progenitor());
	if(mit==hardTree()->particles().end()) continue;
	// put the colour lines in the map
	ShowerParticlePtr oldParticle = cit->first->progenitor();
	ShowerParticlePtr newParticle = mit->second->branchingParticle();
	ShowerParticlePtr newOut;
	ColinePtr cLine = oldParticle-> colourLine();
	ColinePtr aLine = oldParticle->antiColourLine();
	if(newParticle->colourLine() &&
	cmap.find(newParticle-> colourLine())==cmap.end())
	cmap[newParticle-> colourLine()] = cLine;
	if(newParticle->antiColourLine() &&
	cmap.find(newParticle->antiColourLine())==cmap.end())
	cmap[newParticle->antiColourLine()] = aLine;
	// check whether or not particle emits
	bool emission = !mit->second->children().empty();
	if(emission) {
	if(newParticle->colourLine()) {
	ColinePtr ctemp = newParticle-> colourLine();
	ctemp->removeColoured(newParticle);
	}
	if(newParticle->antiColourLine()) {
	ColinePtr ctemp = newParticle->antiColourLine();
	ctemp->removeAntiColoured(newParticle);
	}
	newParticle = mit->second->children()[0]->branchingParticle();
	newOut = mit->second->children()[1]->branchingParticle();
	if(newParticle->id()!=oldParticle->id()&&newParticle->id()==newOut->id())
	swap(newParticle,newOut);
	}
	// get the new colour lines
	ColinePtr newCLine,newALine;
	// sort out colour lines
	if(newParticle->colourLine()) {
	ColinePtr ctemp = newParticle-> colourLine();
	ctemp->removeColoured(newParticle);
	if(cmap.find(ctemp)!=cmap.end()) {
	newCLine = cmap[ctemp];
	}
	else {
	newCLine = new_ptr(ColourLine());
	cmap[ctemp] = newCLine;
	}
	}
	// and anticolour lines
	if(newParticle->antiColourLine()) {
	ColinePtr ctemp = newParticle->antiColourLine();
	ctemp->removeAntiColoured(newParticle);
	if(cmap.find(ctemp)!=cmap.end()) {
	newALine = cmap[ctemp];
	}
	else {
	newALine = new_ptr(ColourLine());
	cmap[ctemp] = newALine;
	}
	}
	// remove colour lines from old particle
	if(aLine) {
	aLine->removeAntiColoured(cit->first->copy());
	aLine->removeAntiColoured(cit->first->progenitor());
	}
	if(cLine) {
	cLine->removeColoured(cit->first->copy());
	cLine->removeColoured(cit->first->progenitor());
	}
	// special for unstable particles
	if(newParticle->id()==oldParticle->id() &&
	(tit!=currentTree()->treelinks().end()\|\|!oldParticle->dataPtr()->stable())) {
	Lorentz5Momentum oldMomentum = oldParticle->momentum();
	Lorentz5Momentum newMomentum = newParticle->momentum();
	LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
	if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false);
	oldParticle->transform(boost);
	boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
	oldParticle->transform(boost);
	if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false);
	newParticle=oldParticle;
	}
	// add particle to colour lines
	if(newCLine) newCLine->addColoured (newParticle);
	if(newALine) newALine->addAntiColoured(newParticle);
	// insert new particles
	cit->first->copy(newParticle);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,true)));
	cit->first->progenitor(sp);
	currentTree()->outgoingLines()[cit->first]=sp;
	cit->first->perturbative(!emission);
	// and the emitted particle if needed
	if(emission) {
	if(newOut->colourLine()) {
	ColinePtr ctemp = newOut-> colourLine();
	ctemp->removeColoured(newOut);
	assert(cmap.find(ctemp)!=cmap.end());
	cmap[ctemp]->addColoured (newOut);
	}
	if(newOut->antiColourLine()) {
	ColinePtr ctemp = newOut->antiColourLine();
	ctemp->removeAntiColoured(newOut);
	assert(cmap.find(ctemp)!=cmap.end());
	cmap[ctemp]->addAntiColoured(newOut);
	}
	ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true));
	ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout));
	out->perturbative(false);
	currentTree()->outgoingLines().insert(make_pair(out,sout));
	}
	// update any decay products
	if(tit!=currentTree()->treelinks().end())
	currentTree()->updateLink(tit->first,make_pair(cit->first,sp));
	}
	// reset the tree
	currentTree()->resetShowerProducts();
	// reextract the particles and set the colour partners
	vector<ShowerParticlePtr> particles =
	currentTree()->extractProgenitorParticles();
	// clear the partners
	for(unsigned int ix=0;ix<particles.size();++ix) {
	particles[ix]->partner(ShowerParticlePtr());
	particles[ix]->clearPartners();
	}
	// clear the tree
	hardTree(HardTreePtr());
	// Set the initial evolution scales
	partnerFinder()->
	setInitialEvolutionScales(particles,!hard,type,!_hardtree);
	}

	Branching QTildeShowerHandler::selectTimeLikeBranching(tShowerParticlePtr particle,
	ShowerInteraction type,
	HardBranchingPtr branch) {
	Branching fb;
	unsigned int iout=0;
	while (true) {
	// break if doing truncated shower and no truncated shower needed
	if(branch && (!isTruncatedShowerON()\|\|hardOnly())) break;
	fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
	// no emission break
	if(!fb.kinematics) break;
	// special for truncated shower
	if(branch) {
	// check haven't evolved too far
	if(fb.kinematics->scale() < branch->scale()) {
	fb=Branching();
	break;
	}
	// find the truncated line
	iout=0;
	if(fb.ids[1]->id()!=fb.ids[2]->id()) {
	if(fb.ids[1]->id()==particle->id()) iout=1;
	else if (fb.ids[2]->id()==particle->id()) iout=2;
	}
	else if(fb.ids[1]->id()==particle->id()) {
	if(fb.kinematics->z()>0.5) iout=1;
	else iout=2;
	}
	// apply the vetos for the truncated shower
	// no flavour changing branchings
	if(iout==0) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
	// only if same interaction for forced branching
	ShowerInteraction type2 = convertInteraction(fb.type);
	// and evolution
	if(type2==branch->sudakov()->interactionType()) {
	if(zsplit < 0.5 \|\| // hardest line veto
	fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// pt veto
	if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// standard vetos for all emissions
	if(timeLikeVetoed(fb,particle)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	continue;
	}
	// special for already decayed particles
	// don't allow flavour changing branchings
	bool vetoDecay = false;
	for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,
	tShowerParticlePtr> >::const_iterator tit = currentTree()->treelinks().begin();
	tit != currentTree()->treelinks().end();++tit) {
	if(tit->second.first == progenitor()) {
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	it = currentTree()->outgoingLines().find(progenitor());
	if(it!=currentTree()->outgoingLines().end() && particle == it->second &&
	fb.ids[0]!=fb.ids[1] && fb.ids[1]!=fb.ids[2]) {
	vetoDecay = true;
	break;
	}
	}
	}
	if(vetoDecay) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	continue;
	}
	break;
	}
	// normal case
	if(!branch) {
	if(fb.kinematics) fb.hard = false;
	return fb;
	}
	// truncated emission
	if(fb.kinematics) {
	fb.hard = false;
	fb.iout = iout;
	return fb;
	}
	// otherwise need to return the hard emission
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin = new_ptr(FS_QTildeShowerKinematics1to2(
	branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT(),
	branch->sudakov()
	));

	IdList idlist(3);
	idlist[0] = particle->dataPtr();
	idlist[1] = branch->children()[0]->branchingParticle()->dataPtr();
	idlist[2] = branch->children()[1]->branchingParticle()->dataPtr();
	fb = Branching( showerKin, idlist, branch->sudakov(),branch->type() );
	fb.hard = true;
	fb.iout=0;
	// return it
	return fb;
	}

	Branching QTildeShowerHandler::selectSpaceLikeDecayBranching(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass,ShowerInteraction type,
	HardBranchingPtr branch) {
	Branching fb;
	unsigned int iout=0;
	while (true) {
	// break if doing truncated shower and no truncated shower needed
	if(branch && (!isTruncatedShowerON()\|\|hardOnly())) break;
	// select branching
	fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass,
	_initialenhance,type);
	// return if no radiation
	if(!fb.kinematics) break;
	// special for truncated shower
	if(branch) {
	// check haven't evolved too far
	if(fb.kinematics->scale() < branch->scale()) {
	fb=Branching();
	break;
	}
	// find the truncated line
	iout=0;
	if(fb.ids[1]->id()!=fb.ids[2]->id()) {
	if(fb.ids[1]->id()==particle->id()) iout=1;
	else if (fb.ids[2]->id()==particle->id()) iout=2;
	}
	else if(fb.ids[1]->id()==particle->id()) {
	if(fb.kinematics->z()>0.5) iout=1;
	else iout=2;
	}
	// apply the vetos for the truncated shower
	// no flavour changing branchings
	if(iout==0) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	ShowerInteraction type2 = convertInteraction(fb.type);
	double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
	if(type2==branch->sudakov()->interactionType()) {
	if(zsplit < 0.5 \|\| // hardest line veto
	fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// pt veto
	if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// if not vetoed break
	if(spaceLikeDecayVetoed(fb,particle)) {
	// otherwise reset scale and continue
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	break;
	}
	// normal case
	if(!branch) {
	if(fb.kinematics) fb.hard = false;
	return fb;
	}
	// truncated emission
	if(fb.kinematics) {
	fb.hard = false;
	fb.iout = iout;
	return fb;
	}
	// otherwise need to return the hard emission
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin = new_ptr(Decay_QTildeShowerKinematics1to2(
	branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT(),
	branch->sudakov()));
	IdList idlist(3);
	idlist[0] = particle->dataPtr();
	idlist[1] = branch->children()[0]->branchingParticle()->dataPtr();
	idlist[2] = branch->children()[1]->branchingParticle()->dataPtr();
	// create the branching
	fb = Branching( showerKin, idlist, branch->sudakov(),ShowerPartnerType::QCDColourLine );
	fb.hard=true;
	fb.iout=0;
	// return it
	return fb;
	}

	void QTildeShowerHandler::checkFlags() {
	string error = "Inconsistent hard emission set-up in QTildeShowerHandler::showerHardProcess(). ";
	if ( ( currentTree()->isMCatNLOSEvent() \|\| currentTree()->isMCatNLOHEvent() ) ) {
	if (_hardEmission ==2 )
	throw Exception() << error
	<< "Cannot generate POWHEG matching with MC@NLO shower "
	<< "approximation. Add 'set QTildeShowerHandler:HardEmission 0' to input file."
	<< Exception::runerror;
	if ( canHandleMatchboxTrunc() )
	throw Exception() << error
	<< "Cannot use truncated qtilde shower with MC@NLO shower "
	<< "approximation. Set LHCGenerator:EventHandler"
	<< ":CascadeHandler to '/Herwig/Shower/ShowerHandler' or "
	<< "'/Herwig/Shower/Dipole/DipoleShowerHandler'."
	<< Exception::runerror;
	}
	else if ( ((currentTree()->isPowhegSEvent() \|\| currentTree()->isPowhegHEvent()) ) &&
	_hardEmission != 2){
	if ( canHandleMatchboxTrunc())
	throw Exception() << error
	<< "Unmatched events requested for POWHEG shower "
	<< "approximation. Set QTildeShowerHandler:HardEmission to "
	<< "'POWHEG'."
	<< Exception::runerror;
	else if (_hardEmissionWarn) {
	_hardEmissionWarn = false;
	_hardEmission=2;
	throw Exception() << error
	<< "Unmatched events requested for POWHEG shower "
	<< "approximation. Changing QTildeShowerHandler:HardEmission from "
	<< _hardEmission << " to 2"
	<< Exception::warning;
	}
	}

	if ( currentTree()->isPowhegSEvent() \|\| currentTree()->isPowhegHEvent()) {
	if (currentTree()->showerApproximation()->needsTruncatedShower() &&
	!canHandleMatchboxTrunc() )
	throw Exception() << error
	<< "Current shower handler cannot generate truncated shower. "
	<< "Set Generator:EventHandler:CascadeHandler to "
	<< "'/Herwig/Shower/PowhegShowerHandler'."
	<< Exception::runerror;
	}
	else if ( currentTree()->truncatedShower() && _missingTruncWarn) {
	_missingTruncWarn=false;
	throw Exception() << "Warning: POWHEG shower approximation used without "
	<< "truncated shower. Set Generator:EventHandler:"
	<< "CascadeHandler to '/Herwig/Shower/PowhegShowerHandler' and "
	<< "'MEMatching:TruncatedShower Yes'."
	<< Exception::warning;
	}
	// else if ( !dipme && _hardEmissionMode > 1 &&
	// firstInteraction())
	// throw Exception() << error
	// << "POWHEG matching requested for LO events. Include "
	// << "'set Factory:ShowerApproximation MEMatching' in input file."
	// << Exception::runerror;
	}


	tPPair QTildeShowerHandler::remakeRemnant(tPPair oldp){
	// get the parton extractor
	PartonExtractor & pex = *lastExtractor();
	// get the new partons
	tPPair newp = make_pair(findFirstParton(oldp.first ),
	findFirstParton(oldp.second));
	// if the same do nothing
	if(newp == oldp) return oldp;
	// Creates the new remnants and returns the new PartonBinInstances
	// ATTENTION Broken here for very strange configuration
	PBIPair newbins = pex.newRemnants(oldp, newp, newStep());
	newStep()->addIntermediate(newp.first);
	newStep()->addIntermediate(newp.second);
	// return the new partons
	return newp;
	}

	PPtr QTildeShowerHandler::findFirstParton(tPPtr seed) const{
	if(seed->parents().empty()) return seed;
	tPPtr parent = seed->parents()[0];
	//if no parent there this is a loose end which will
	//be connected to the remnant soon.
	if(!parent \|\| parent == incomingBeams().first \|\|
	parent == incomingBeams().second ) return seed;
	else return findFirstParton(parent);
	}

	void QTildeShowerHandler::decay(ShowerTreePtr tree, ShowerDecayMap & decay) {
	// must be one incoming particle
	assert(tree->incomingLines().size()==1);
	// apply any transforms
	tree->applyTransforms();
	// if already decayed return
	if(!tree->outgoingLines().empty()) return;
	// now we need to replace the particle with a new copy after the shower
	// find particle after the shower
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	tit = tree->parent()->treelinks().find(tree);
	assert(tit!=tree->parent()->treelinks().end());
	ShowerParticlePtr newparent=tit->second.second;
	PerturbativeProcessPtr newProcess = new_ptr(PerturbativeProcess());
	newProcess->incoming().push_back(make_pair(newparent,PerturbativeProcessPtr()));
	DecayProcessMap decayMap;
	ShowerHandler::decay(newProcess,decayMap);
	ShowerTree::constructTrees(tree,decay,newProcess,decayMap);
	}

	namespace {

	ShowerProgenitorPtr
	findFinalStateLine(ShowerTreePtr tree, long id, Lorentz5Momentum momentum) {
	map<ShowerProgenitorPtr,tShowerParticlePtr>::iterator partner;
	Energy2 dmin(1e30*GeV2);
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::iterator
	cit =tree->outgoingLines().begin(); cit!=tree->outgoingLines().end(); ++cit) {
	if(cit->second->id()!=id) continue;
	Energy2 test =
	sqr(cit->second->momentum().x()-momentum.x())+
	sqr(cit->second->momentum().y()-momentum.y())+
	sqr(cit->second->momentum().z()-momentum.z())+
	sqr(cit->second->momentum().t()-momentum.t());
	if(test<dmin) {
	dmin = test;
	partner = cit;
	}
	}
	return partner->first;
	}

	ShowerProgenitorPtr
	findInitialStateLine(ShowerTreePtr tree, long id, Lorentz5Momentum momentum) {
	map<ShowerProgenitorPtr,ShowerParticlePtr>::iterator partner;
	Energy2 dmin(1e30*GeV2);
	for(map<ShowerProgenitorPtr,ShowerParticlePtr>::iterator
	cit =tree->incomingLines().begin(); cit!=tree->incomingLines().end(); ++cit) {
	if(cit->second->id()!=id) continue;
	Energy2 test =
	sqr(cit->second->momentum().x()-momentum.x())+
	sqr(cit->second->momentum().y()-momentum.y())+
	sqr(cit->second->momentum().z()-momentum.z())+
	sqr(cit->second->momentum().t()-momentum.t());
	if(test<dmin) {
	dmin = test;
	partner = cit;
	}
	}
	return partner->first;
	}

	void fixSpectatorColours(PPtr newSpect,ShowerProgenitorPtr oldSpect,
	ColinePair & cline,ColinePair & aline, bool reconnect) {
	cline.first = oldSpect->progenitor()->colourLine();
	cline.second = newSpect->colourLine();
	aline.first = oldSpect->progenitor()->antiColourLine();
	aline.second = newSpect->antiColourLine();
	if(!reconnect) return;
	if(cline.first) {
	cline.first ->removeColoured(oldSpect->copy());
	cline.first ->removeColoured(oldSpect->progenitor());
	cline.second->removeColoured(newSpect);
	cline.first ->addColoured(newSpect);
	}
	if(aline.first) {
	aline.first ->removeAntiColoured(oldSpect->copy());
	aline.first ->removeAntiColoured(oldSpect->progenitor());
	aline.second->removeAntiColoured(newSpect);
	aline.first ->addAntiColoured(newSpect);
	}
	}

	void fixInitialStateEmitter(ShowerTreePtr tree, PPtr newEmit,PPtr emitted, ShowerProgenitorPtr emitter,
	ColinePair cline,ColinePair aline,double x) {
	// sort out the colours
	if(emitted->dataPtr()->iColour()==PDT::Colour8) {
	// emitter
	if(cline.first && cline.first == emitter->progenitor()->antiColourLine() &&
	cline.second !=newEmit->antiColourLine()) {
	// sort out not radiating line
	ColinePtr col = emitter->progenitor()->colourLine();
	if(col) {
	col->removeColoured(emitter->copy());
	col->removeColoured(emitter->progenitor());
	newEmit->colourLine()->removeColoured(newEmit);
	col->addColoured(newEmit);
	}
	}
	else if(aline.first && aline.first == emitter->progenitor()->colourLine() &&
	aline.second !=newEmit->colourLine()) {
	// sort out not radiating line
	ColinePtr anti = emitter->progenitor()->antiColourLine();
	if(anti) {
	anti->removeAntiColoured(emitter->copy());
	anti->removeAntiColoured(emitter->progenitor());
	newEmit->colourLine()->removeAntiColoured(newEmit);
	anti->addAntiColoured(newEmit);
	}
	}
	else
	assert(false);
	// emitted
	if(cline.first && cline.second==emitted->colourLine()) {
	cline.second->removeColoured(emitted);
	cline.first->addColoured(emitted);
	}
	else if(aline.first && aline.second==emitted->antiColourLine()) {
	aline.second->removeAntiColoured(emitted);
	aline.first->addAntiColoured(emitted);
	}
	else
	assert(false);
	}
	else {
	if(emitter->progenitor()->antiColourLine() ) {
	ColinePtr col = emitter->progenitor()->antiColourLine();
	col->removeAntiColoured(emitter->copy());
	col->removeAntiColoured(emitter->progenitor());
	if(newEmit->antiColourLine()) {
	newEmit->antiColourLine()->removeAntiColoured(newEmit);
	col->addAntiColoured(newEmit);
	}
	else if (emitted->colourLine()) {
	emitted->colourLine()->removeColoured(emitted);
	col->addColoured(emitted);
	}
	else
	assert(false);
	}
	if(emitter->progenitor()->colourLine() ) {
	ColinePtr col = emitter->progenitor()->colourLine();
	col->removeColoured(emitter->copy());
	col->removeColoured(emitter->progenitor());
	if(newEmit->colourLine()) {
	newEmit->colourLine()->removeColoured(newEmit);
	col->addColoured(newEmit);
	}
	else if (emitted->antiColourLine()) {
	emitted->antiColourLine()->removeAntiColoured(emitted);
	col->addAntiColoured(emitted);
	}
	else
	assert(false);
	}
	}
	// update the emitter
	emitter->copy(newEmit);
	ShowerParticlePtr sp = new_ptr(ShowerParticle(*newEmit,1,false));
	sp->x(x);
	emitter->progenitor(sp);
	tree->incomingLines()[emitter]=sp;
	emitter->perturbative(false);
	// add emitted
	sp=new_ptr(ShowerParticle(*emitted,1,true));
	ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(emitter->original(),emitted,sp));
	gluon->perturbative(false);
	tree->outgoingLines().insert(make_pair(gluon,sp));
	}

	void fixFinalStateEmitter(ShowerTreePtr tree, PPtr newEmit,PPtr emitted, ShowerProgenitorPtr emitter,
	ColinePair cline,ColinePair aline) {
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
	// special case if decayed
	for(tit = tree->treelinks().begin(); tit != tree->treelinks().end();++tit) {
	if(tit->second.first && tit->second.second==emitter->progenitor())
	break;
	}
	// sort out the colour lines
	if(cline.first && cline.first == emitter->progenitor()->antiColourLine() &&
	cline.second !=newEmit->antiColourLine()) {
	// sort out not radiating line
	ColinePtr col = emitter->progenitor()->colourLine();
	if(col) {
	col->removeColoured(emitter->copy());
	col->removeColoured(emitter->progenitor());
	newEmit->colourLine()->removeColoured(newEmit);
	col->addColoured(newEmit);
	}
	}
	else if(aline.first && aline.first == emitter->progenitor()->colourLine() &&
	aline.second !=newEmit->colourLine()) {
	// sort out not radiating line
	ColinePtr anti = emitter->progenitor()->antiColourLine();
	if(anti) {
	anti->removeAntiColoured(emitter->copy());
	anti->removeAntiColoured(emitter->progenitor());
	newEmit->colourLine()->removeAntiColoured(newEmit);
	anti->addAntiColoured(newEmit);
	}
	}
	else
	assert(false);
	// update the emitter
	emitter->copy(newEmit);
	ShowerParticlePtr sp = new_ptr(ShowerParticle(*newEmit,1,true));
	emitter->progenitor(sp);
	tree->outgoingLines()[emitter]=sp;
	emitter->perturbative(false);
	// update for decaying particles
	if(tit!=tree->treelinks().end())
	tree->updateLink(tit->first,make_pair(emitter,sp));
	// add the emitted particle
	// sort out the colour
	if(cline.first && cline.second==emitted->antiColourLine()) {
	cline.second->removeAntiColoured(emitted);
	cline.first->addAntiColoured(emitted);
	}
	else if(aline.first && aline.second==emitted->colourLine()) {
	aline.second->removeColoured(emitted);
	aline.first->addColoured(emitted);
	}
	else
	assert(false);
	sp=new_ptr(ShowerParticle(*emitted,1,true));
	ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(emitter->original(),
	emitted,sp));
	gluon->perturbative(false);
	tree->outgoingLines().insert(make_pair(gluon,sp));
	}

	}

	void QTildeShowerHandler::setupMECorrection(RealEmissionProcessPtr real) {
	assert(real);
	// II emission
	if(real->emitter() < real->incoming().size() &&
	real->spectator() < real->incoming().size()) {
	// recoiling system
	for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt= currentTree()->outgoingLines().begin();
	cjt != currentTree()->outgoingLines().end();++cjt ) {
	cjt->first->progenitor()->transform(real->transformation());
	cjt->first->copy()->transform(real->transformation());
	}
	// the the radiating system
	ShowerProgenitorPtr emitter,spectator;
	unsigned int iemit = real->emitter();
	unsigned int ispect = real->spectator();
	int ig = int(real->emitted())-int(real->incoming().size());
	emitter = findInitialStateLine(currentTree(),
	real->bornIncoming()[iemit]->id(),
	real->bornIncoming()[iemit]->momentum());
	spectator = findInitialStateLine(currentTree(),
	real->bornIncoming()[ispect]->id(),
	real->bornIncoming()[ispect]->momentum());
	// sort out the colours
	ColinePair cline,aline;
	fixSpectatorColours(real->incoming()[ispect],spectator,cline,aline,true);
	// update the spectator
	spectator->copy(real->incoming()[ispect]);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->incoming()[ispect],1,false)));
	sp->x(ispect ==0 ? real->x().first :real->x().second);
	spectator->progenitor(sp);
	currentTree()->incomingLines()[spectator]=sp;
	spectator->perturbative(true);
	// now for the emitter
	fixInitialStateEmitter(currentTree(),real->incoming()[iemit],real->outgoing()[ig],
	emitter,cline,aline,iemit ==0 ? real->x().first :real->x().second);
	}
	// FF emission
	else if(real->emitter() >= real->incoming().size() &&
	real->spectator() >= real->incoming().size()) {
	assert(real->outgoing()[real->emitted()-real->incoming().size()]->id()==ParticleID::g);
	// find the emitter and spectator in the shower tree
	ShowerProgenitorPtr emitter,spectator;
	int iemit = int(real->emitter())-int(real->incoming().size());
	emitter = findFinalStateLine(currentTree(),
	real->bornOutgoing()[iemit]->id(),
	real->bornOutgoing()[iemit]->momentum());
	int ispect = int(real->spectator())-int(real->incoming().size());
	spectator = findFinalStateLine(currentTree(),
	real->bornOutgoing()[ispect]->id(),
	real->bornOutgoing()[ispect]->momentum());
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
	// first the spectator
	// special case if decayed
	for(tit = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) {
	if(tit->second.first && tit->second.second==spectator->progenitor())
	break;
	}
	// sort out the colours
	ColinePair cline,aline;
	fixSpectatorColours(real->outgoing()[ispect],spectator,cline,aline,true);
	// update the spectator
	spectator->copy(real->outgoing()[ispect]);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->outgoing()[ispect],1,true)));
	spectator->progenitor(sp);
	currentTree()->outgoingLines()[spectator]=sp;
	spectator->perturbative(true);
	// update for decaying particles
	if(tit!=currentTree()->treelinks().end())
	currentTree()->updateLink(tit->first,make_pair(spectator,sp));
	// now the emitting particle
	int ig = int(real->emitted())-int(real->incoming().size());
	fixFinalStateEmitter(currentTree(),real->outgoing()[iemit],
	real->outgoing()[ig],
	emitter,cline,aline);
	}
	// IF emission
	else {
	// scattering process
	if(real->incoming().size()==2) {
	ShowerProgenitorPtr emitter,spectator;
	unsigned int iemit = real->emitter();
	unsigned int ispect = real->spectator();
	int ig = int(real->emitted())-int(real->incoming().size());
	ColinePair cline,aline;
	// incoming spectator
	if(ispect<2) {
	spectator = findInitialStateLine(currentTree(),
	real->bornIncoming()[ispect]->id(),
	real->bornIncoming()[ispect]->momentum());
	fixSpectatorColours(real->incoming()[ispect],spectator,cline,aline,true);
	// update the spectator
	spectator->copy(real->incoming()[ispect]);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->incoming()[ispect],1,false)));
	sp->x(ispect ==0 ? real->x().first :real->x().second);
	spectator->progenitor(sp);
	currentTree()->incomingLines()[spectator]=sp;
	spectator->perturbative(true);
	}
	// outgoing spectator
	else {
	spectator = findFinalStateLine(currentTree(),
	real->bornOutgoing()[ispect-real->incoming().size()]->id(),
	real->bornOutgoing()[ispect-real->incoming().size()]->momentum());
	// special case if decayed
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
	for(tit = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) {
	if(tit->second.first && tit->second.second==spectator->progenitor())
	break;
	}
	fixSpectatorColours(real->outgoing()[ispect-real->incoming().size()],spectator,cline,aline,true);
	// update the spectator
	spectator->copy(real->outgoing()[ispect-real->incoming().size()]);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->outgoing()[ispect-real->incoming().size()],1,true)));
	spectator->progenitor(sp);
	currentTree()->outgoingLines()[spectator]=sp;
	spectator->perturbative(true);
	// update for decaying particles
	if(tit!=currentTree()->treelinks().end())
	currentTree()->updateLink(tit->first,make_pair(spectator,sp));
	}
	// incoming emitter
	if(iemit<2) {
	emitter = findInitialStateLine(currentTree(),
	real->bornIncoming()[iemit]->id(),
	real->bornIncoming()[iemit]->momentum());
	fixInitialStateEmitter(currentTree(),real->incoming()[iemit],real->outgoing()[ig],
	emitter,aline,cline,iemit ==0 ? real->x().first :real->x().second);
	}
	// outgoing emitter
	else {
	emitter = findFinalStateLine(currentTree(),
	real->bornOutgoing()[iemit-real->incoming().size()]->id(),
	real->bornOutgoing()[iemit-real->incoming().size()]->momentum());
	fixFinalStateEmitter(currentTree(),real->outgoing()[iemit-real->incoming().size()],
	real->outgoing()[ig],emitter,aline,cline);
	}
	}
	// decay process
	else {
	assert(real->spectator()==0);
	unsigned int iemit = real->emitter()-real->incoming().size();
	int ig = int(real->emitted())-int(real->incoming().size());
	ColinePair cline,aline;
	// incoming spectator
	ShowerProgenitorPtr spectator = findInitialStateLine(currentTree(),
	real->bornIncoming()[0]->id(),
	real->bornIncoming()[0]->momentum());
	fixSpectatorColours(real->incoming()[0],spectator,cline,aline,false);
	// find the emitter
	ShowerProgenitorPtr emitter =
	findFinalStateLine(currentTree(),
	real->bornOutgoing()[iemit]->id(),
	real->bornOutgoing()[iemit]->momentum());
	// recoiling system
	for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt= currentTree()->outgoingLines().begin();
	cjt != currentTree()->outgoingLines().end();++cjt ) {
	if(cjt->first==emitter) continue;
	cjt->first->progenitor()->transform(real->transformation());
	cjt->first->copy()->transform(real->transformation());
	}
	// sort out the emitter
	fixFinalStateEmitter(currentTree(),real->outgoing()[iemit],
	real->outgoing()[ig],emitter,aline,cline);
	}
	}
	// clean up the shower tree
	_currenttree->resetShowerProducts();
	}
	diff --git a/Shower/QTilde/SplittingFunctions/SplittingFunction.cc b/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
	--- a/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
	+++ b/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
	@@ -1,1043 +1,1066 @@
	// -- C++ --
	//
	// SplittingFunction.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 SplittingFunction class.
	//

	#include "SplittingFunction.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Utilities/EnumIO.h"
	#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
	#include "ThePEG/Utilities/DescribeClass.h"

	using namespace Herwig;

	DescribeAbstractClass<SplittingFunction,Interfaced>
	describeSplittingFunction ("Herwig::SplittingFunction","");

	void SplittingFunction::Init() {

	static ClassDocumentation<SplittingFunction> documentation
	("The SplittingFunction class is the based class for 1->2 splitting functions"
	" in Herwig");

	static Switch<SplittingFunction,ColourStructure> interfaceColourStructure
	("ColourStructure",
	"The colour structure for the splitting function",
	&SplittingFunction::_colourStructure, Undefined, false, false);
	static SwitchOption interfaceColourStructureTripletTripletOctet
	(interfaceColourStructure,
	"TripletTripletOctet",
	"3 -> 3 8",
	TripletTripletOctet);
	static SwitchOption interfaceColourStructureOctetOctetOctet
	(interfaceColourStructure,
	"OctetOctetOctet",
	"8 -> 8 8",
	OctetOctetOctet);
	static SwitchOption interfaceColourStructureOctetTripletTriplet
	(interfaceColourStructure,
	"OctetTripletTriplet",
	"8 -> 3 3bar",
	OctetTripletTriplet);
	static SwitchOption interfaceColourStructureTripletOctetTriplet
	(interfaceColourStructure,
	"TripletOctetTriplet",
	"3 -> 8 3",
	TripletOctetTriplet);
	static SwitchOption interfaceColourStructureSextetSextetOctet
	(interfaceColourStructure,
	"SextetSextetOctet",
	"6 -> 6 8",
	SextetSextetOctet);

	static SwitchOption interfaceColourStructureChargedChargedNeutral
	(interfaceColourStructure,
	"ChargedChargedNeutral",
	"q -> q 0",
	ChargedChargedNeutral);
	static SwitchOption interfaceColourStructureNeutralChargedCharged
	(interfaceColourStructure,
	"NeutralChargedCharged",
	"0 -> q qbar",
	NeutralChargedCharged);
	static SwitchOption interfaceColourStructureChargedNeutralCharged
	(interfaceColourStructure,
	"ChargedNeutralCharged",
	"q -> 0 q",
	ChargedNeutralCharged);
	static SwitchOption interfaceColourStructureEW
	(interfaceColourStructure,
	"EW",
	"q -> q W/Z",
	EW);

	static Switch<SplittingFunction,ShowerInteraction>
	interfaceInteractionType
	("InteractionType",
	"Type of the interaction",
	&SplittingFunction::_interactionType,
	ShowerInteraction::UNDEFINED, false, false);
	static SwitchOption interfaceInteractionTypeQCD
	(interfaceInteractionType,
	"QCD","QCD",ShowerInteraction::QCD);
	static SwitchOption interfaceInteractionTypeQED
	(interfaceInteractionType,
	"QED","QED",ShowerInteraction::QED);
	static SwitchOption interfaceInteractionTypeEW
	(interfaceInteractionType,
	"EW","EW",ShowerInteraction::EW);

	static Switch<SplittingFunction,bool> interfaceAngularOrdered
	("AngularOrdered",
	"Whether or not this interaction is angular ordered, "
	"normally only g->q qbar and gamma-> f fbar are the only ones which aren't.",
	&SplittingFunction::angularOrdered_, true, false, false);
	static SwitchOption interfaceAngularOrderedYes
	(interfaceAngularOrdered,
	"Yes",
	"Interaction is angular ordered",
	true);
	static SwitchOption interfaceAngularOrderedNo
	(interfaceAngularOrdered,
	"No",
	"Interaction isn't angular ordered",
	false);

	static Switch<SplittingFunction,unsigned int> interfaceScaleChoice
	("ScaleChoice",
	"The scale choice to be used",
	&SplittingFunction::scaleChoice_, 2, false, false);
	static SwitchOption interfaceScaleChoicepT
	(interfaceScaleChoice,
	"pT",
	"pT of the branching",
	0);
	static SwitchOption interfaceScaleChoiceQ2
	(interfaceScaleChoice,
	"Q2",
	"Q2 of the branching",
	1);
	static SwitchOption interfaceScaleChoiceFromAngularOrdering
	(interfaceScaleChoice,
	"FromAngularOrdering",
	"If angular order use pT, otherwise Q2",
	2);

	}

	void SplittingFunction::persistentOutput(PersistentOStream & os) const {
	os << oenum(_interactionType)
	<< oenum(_colourStructure) << _colourFactor
	<< angularOrdered_ << scaleChoice_;
	}

	void SplittingFunction::persistentInput(PersistentIStream & is, int) {
	is >> ienum(_interactionType)
	>> ienum(_colourStructure) >> _colourFactor
	>> angularOrdered_ >> scaleChoice_;
	}

	void SplittingFunction::colourConnection(tShowerParticlePtr parent,
	tShowerParticlePtr first,
	tShowerParticlePtr second,
	ShowerPartnerType partnerType,
	const bool back) const {
	if(_colourStructure==TripletTripletOctet) {
	if(!back) {
	ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// ensure input consistency
	assert(( cparent.first && !cparent.second &&
	partnerType==ShowerPartnerType::QCDColourLine) \|\|
	( !cparent.first && cparent.second &&
	partnerType==ShowerPartnerType::QCDAntiColourLine));
	// q -> q g
	if(cparent.first) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(second);
	newline->addColoured ( first);
	newline->addAntiColoured (second);
	}
	// qbar -> qbar g
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.second->addAntiColoured(second);
	newline->addColoured(second);
	newline->addAntiColoured(first);
	}
	// Set progenitor
	first->progenitor(parent->progenitor());
	second->progenitor(parent->progenitor());
	}
	else {
	ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// ensure input consistency
	assert(( cfirst.first && !cfirst.second &&
	partnerType==ShowerPartnerType::QCDColourLine) \|\|
	( !cfirst.first && cfirst.second &&
	partnerType==ShowerPartnerType::QCDAntiColourLine));
	// q -> q g
	if(cfirst.first) {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.first->addAntiColoured(second);
	newline->addColoured(second);
	newline->addColoured(parent);
	}
	// qbar -> qbar g
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.second->addColoured(second);
	newline->addAntiColoured(second);
	newline->addAntiColoured(parent);
	}
	// Set progenitor
	parent->progenitor(first->progenitor());
	second->progenitor(first->progenitor());
	}
	}
	else if(_colourStructure==OctetOctetOctet) {
	if(!back) {
	ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// ensure input consistency
	assert(cparent.first&&cparent.second);
	// ensure first gluon is hardest
	if( first->id()==second->id() && parent->showerKinematics()->z()<0.5 )
	swap(first,second);
	// colour line radiates
	if(partnerType==ShowerPartnerType::QCDColourLine) {
	// The colour line is radiating
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(second);
	cparent.second->addAntiColoured(first);
	newline->addColoured(first);
	newline->addAntiColoured(second);
	}
	// anti colour line radiates
	else if(partnerType==ShowerPartnerType::QCDAntiColourLine) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(first);
	cparent.second->addAntiColoured(second);
	newline->addColoured(second);
	newline->addAntiColoured(first);
	}
	else
	assert(false);
	}
	else {
	ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// ensure input consistency
	assert(cfirst.first&&cfirst.second);
	// The colour line is radiating
	if(partnerType==ShowerPartnerType::QCDColourLine) {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.first->addAntiColoured(second);
	cfirst.second->addAntiColoured(parent);
	newline->addColoured(parent);
	newline->addColoured(second);
	}
	// anti colour line radiates
	else if(partnerType==ShowerPartnerType::QCDAntiColourLine) {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.first->addColoured(parent);
	cfirst.second->addColoured(second);
	newline->addAntiColoured(second);
	newline->addAntiColoured(parent);
	}
	else
	assert(false);
	}
	}
	else if(_colourStructure == OctetTripletTriplet) {
	if(!back) {
	ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// ensure input consistency
	assert(cparent.first&&cparent.second);
	cparent.first ->addColoured ( first);
	cparent.second->addAntiColoured(second);
	// Set progenitor
	first->progenitor(parent->progenitor());
	second->progenitor(parent->progenitor());
	}
	else {
	ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// ensure input consistency
	assert(( cfirst.first && !cfirst.second) \|\|
	(!cfirst.first && cfirst.second));
	// g -> q qbar
	if(cfirst.first) {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.first->addColoured(parent);
	newline->addAntiColoured(second);
	newline->addAntiColoured(parent);
	}
	// g -> qbar q
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cfirst.second->addAntiColoured(parent);
	newline->addColoured(second);
	newline->addColoured(parent);
	}
	// Set progenitor
	parent->progenitor(first->progenitor());
	second->progenitor(first->progenitor());
	}
	}
	else if(_colourStructure == TripletOctetTriplet) {
	if(!back) {
	ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// ensure input consistency
	assert(( cparent.first && !cparent.second) \|\|
	(!cparent.first && cparent.second));
	// q -> g q
	if(cparent.first) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(first);
	newline->addColoured (second);
	newline->addAntiColoured( first);
	}
	// qbar -> g qbar
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.second->addAntiColoured(first);
	newline->addColoured ( first);
	newline->addAntiColoured(second);
	}
	// Set progenitor
	first->progenitor(parent->progenitor());
	second->progenitor(parent->progenitor());
	}
	else {
	ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// ensure input consistency
	assert(cfirst.first&&cfirst.second);
	// q -> g q
	if(parent->id()>0) {
	cfirst.first ->addColoured(parent);
	cfirst.second->addColoured(second);
	}
	else {
	cfirst.first ->addAntiColoured(second);
	cfirst.second->addAntiColoured(parent);
	}
	// Set progenitor
	parent->progenitor(first->progenitor());
	second->progenitor(first->progenitor());
	}
	}
	else if(_colourStructure==SextetSextetOctet) {
	//make sure we're not doing backward evolution
	assert(!back);

	//make sure something sensible
	assert(parent->colourLine() \|\| parent->antiColourLine());

	//get the colour lines or anti-colour lines
	bool isAntiColour=true;
	ColinePair cparent;
	if(parent->colourLine()) {
	cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[0]),
	const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[1]));
	isAntiColour=false;
	}
	else {
	cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[0]),
	const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[1]));
	}

	//check for sensible input
	// assert(cparent.first && cparent.second);

	// sextet has 2 colour lines
	if(!isAntiColour) {
	//pick at random which of the colour topolgies to take
	double topology = UseRandom::rnd();
	if(topology < 0.25) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(second);
	cparent.second->addColoured(first);
	newline->addColoured(first);
	newline->addAntiColoured(second);
	}
	else if(topology >=0.25 && topology < 0.5) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(first);
	cparent.second->addColoured(second);
	newline->addColoured(first);
	newline->addAntiColoured(second);
	}
	else if(topology >= 0.5 && topology < 0.75) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(second);
	cparent.second->addColoured(first);
	newline->addColoured(first);
	newline->addAntiColoured(second);
	}
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addColoured(first);
	cparent.second->addColoured(second);
	newline->addColoured(first);
	newline->addAntiColoured(second);
	}
	}
	// sextet has 2 anti-colour lines
	else {
	double topology = UseRandom::rnd();
	if(topology < 0.25){
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addAntiColoured(second);
	cparent.second->addAntiColoured(first);
	newline->addAntiColoured(first);
	newline->addColoured(second);
	}
	else if(topology >=0.25 && topology < 0.5) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addAntiColoured(first);
	cparent.second->addAntiColoured(second);
	newline->addAntiColoured(first);
	newline->addColoured(second);
	}
	else if(topology >= 0.5 && topology < 0.75) {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addAntiColoured(second);
	cparent.second->addAntiColoured(first);
	newline->addAntiColoured(first);
	newline->addColoured(second);
	}
	else {
	ColinePtr newline=new_ptr(ColourLine());
	cparent.first->addAntiColoured(first);
	cparent.second->addAntiColoured(second);
	newline->addAntiColoured(first);
	newline->addColoured(second);
	}
	}
	}
	else if(_colourStructure == ChargedChargedNeutral) {
	if(!parent->data().coloured()) return;
	if(!back) {
	ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// q -> q g
	if(cparent.first) {
	cparent.first->addColoured(first);
	}
	// qbar -> qbar g
	if(cparent.second) {
	cparent.second->addAntiColoured(first);
	}
	}
	else {
	ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// q -> q g
	if(cfirst.first) {
	cfirst.first->addColoured(parent);
	}
	// qbar -> qbar g
	if(cfirst.second) {
	cfirst.second->addAntiColoured(parent);
	}
	}
	}
	else if(_colourStructure == ChargedNeutralCharged) {
	if(!parent->data().coloured()) return;
	if(!back) {
	ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// q -> q g
	if(cparent.first) {
	cparent.first->addColoured(second);
	}
	// qbar -> qbar g
	if(cparent.second) {
	cparent.second->addAntiColoured(second);
	}
	}
	else {
	if (second->dataPtr()->iColour()==PDT::Colour3 ) {
	ColinePtr newline=new_ptr(ColourLine());
	newline->addColoured(second);
	newline->addColoured(parent);
	}
	else if (second->dataPtr()->iColour()==PDT::Colour3bar ) {
	ColinePtr newline=new_ptr(ColourLine());
	newline->addAntiColoured(second);
	newline->addAntiColoured(parent);
	}
	}
	}
	else if(_colourStructure == NeutralChargedCharged ) {
	if(!back) {
	if(first->dataPtr()->coloured()) {
	ColinePtr newline=new_ptr(ColourLine());
	if(first->dataPtr()->iColour()==PDT::Colour3) {
	newline->addColoured (first );
	newline->addAntiColoured(second);
	}
	else if (first->dataPtr()->iColour()==PDT::Colour3bar) {
	newline->addColoured (second);
	newline->addAntiColoured(first );
	}
	else
	assert(false);
	}
	}
	else {
	ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// gamma -> q qbar
	if(cfirst.first) {
	cfirst.first->addAntiColoured(second);
	}
	// gamma -> qbar q
	else if(cfirst.second) {
	cfirst.second->addColoured(second);
	}
	else
	assert(false);
	}
	}
	else if(_colourStructure == EW) {
	if(!parent->data().coloured()) return;
	if(!back) {
	ColinePair cparent = ColinePair(parent->colourLine(),
	parent->antiColourLine());
	// q -> q g
	if(cparent.first) {
	cparent.first->addColoured(first);
	}
	// qbar -> qbar g
	if(cparent.second) {
	cparent.second->addAntiColoured(first);
	}
	}
	else {
	ColinePair cfirst = ColinePair(first->colourLine(),
	first->antiColourLine());
	// q -> q g
	if(cfirst.first) {
	cfirst.first->addColoured(parent);
	}
	// qbar -> qbar g
	if(cfirst.second) {
	cfirst.second->addAntiColoured(parent);
	}
	}
	}
	else {
	assert(false);
	}
	}

	void SplittingFunction::doinit() {
	Interfaced::doinit();
	assert(_interactionType!=ShowerInteraction::UNDEFINED);
	assert((_colourStructure>0&&_interactionType==ShowerInteraction::QCD) \|\|
	(_colourStructure<0&&(_interactionType==ShowerInteraction::QED \|\|
	_interactionType==ShowerInteraction::EW)) );
	if(_colourFactor>0.) return;
	// compute the colour factors if need
	if(_colourStructure==TripletTripletOctet) {
	_colourFactor = 4./3.;
	}
	else if(_colourStructure==OctetOctetOctet) {
	_colourFactor = 3.;
	}
	else if(_colourStructure==OctetTripletTriplet) {
	_colourFactor = 0.5;
	}
	else if(_colourStructure==TripletOctetTriplet) {
	_colourFactor = 4./3.;
	}
	else if(_colourStructure==SextetSextetOctet) {
	_colourFactor = 10./3.;
	}
	else if(_colourStructure<0) {
	_colourFactor = 1.;
	}
	else {
	assert(false);
	}
	}

	bool SplittingFunction::checkColours(const IdList & ids) const {
	if(_colourStructure==TripletTripletOctet) {
	if(ids[0]!=ids[1]) return false;
	if((ids[0]->iColour()==PDT::Colour3\|\|ids[0]->iColour()==PDT::Colour3bar) &&
	ids[2]->iColour()==PDT::Colour8) return true;
	return false;
	}
	else if(_colourStructure==OctetOctetOctet) {
	for(unsigned int ix=0;ix<3;++ix) {
	if(ids[ix]->iColour()!=PDT::Colour8) return false;
	}
	return true;
	}
	else if(_colourStructure==OctetTripletTriplet) {
	if(ids[0]->iColour()!=PDT::Colour8) return false;
	if(ids[1]->iColour()==PDT::Colour3&&ids[2]->iColour()==PDT::Colour3bar)
	return true;
	if(ids[1]->iColour()==PDT::Colour3bar&&ids[2]->iColour()==PDT::Colour3)
	return true;
	return false;
	}
	else if(_colourStructure==TripletOctetTriplet) {
	if(ids[0]!=ids[2]) return false;
	if((ids[0]->iColour()==PDT::Colour3\|\|ids[0]->iColour()==PDT::Colour3bar) &&
	ids[1]->iColour()==PDT::Colour8) return true;
	return false;
	}
	else if(_colourStructure==SextetSextetOctet) {
	if(ids[0]!=ids[1]) return false;
	if((ids[0]->iColour()==PDT::Colour6 \|\| ids[0]->iColour()==PDT::Colour6bar) &&
	ids[2]->iColour()==PDT::Colour8) return true;
	return false;
	}
	else if(_colourStructure==ChargedChargedNeutral) {
	if(ids[0]!=ids[1]) return false;
	if(ids[2]->iCharge()!=0) return false;
	if(ids[0]->iCharge()==ids[1]->iCharge()) return true;
	return false;
	}
	else if(_colourStructure==ChargedNeutralCharged) {
	if(ids[0]!=ids[2]) return false;
	if(ids[1]->iCharge()!=0) return false;
	if(ids[0]->iCharge()==ids[2]->iCharge()) return true;
	return false;
	}
	else if(_colourStructure==NeutralChargedCharged) {
	if(ids[1]->id()!=-ids[2]->id()) return false;
	if(ids[0]->iCharge()!=0) return false;
	if(ids[1]->iCharge()==-ids[2]->iCharge()) return true;
	return false;
	}
	else {
	assert(false);
	}
	return false;
	}

	namespace {

	bool hasColour(tPPtr p) {
	PDT::Colour colour = p->dataPtr()->iColour();
	return colour==PDT::Colour3 \|\| colour==PDT::Colour8 \|\| colour == PDT::Colour6;
	}

	bool hasAntiColour(tPPtr p) {
	PDT::Colour colour = p->dataPtr()->iColour();
	return colour==PDT::Colour3bar \|\| colour==PDT::Colour8 \|\| colour == PDT::Colour6bar;
	}

	}

	void SplittingFunction::evaluateFinalStateScales(ShowerPartnerType partnerType,
	Energy scale, double z,
	tShowerParticlePtr parent,
	tShowerParticlePtr emitter,
	tShowerParticlePtr emitted) {
	// identify emitter and emitted
	double zEmitter = z, zEmitted = 1.-z;
	bool bosonSplitting(false);
	// special for g -> gg, particle highest z is emitter
	if(emitter->id() == emitted->id() && emitter->id() == parent->id() &&
	zEmitted > zEmitter) {
	swap(zEmitted,zEmitter);
	swap( emitted, emitter);
	}
	// otherwise if particle ID same
	else if(emitted->id()==parent->id()) {
	swap(zEmitted,zEmitter);
	swap( emitted, emitter);
	}
	// no real emitter/emitted
	else if(emitter->id()!=parent->id()) {
	bosonSplitting = true;
	}
	// may need to add angularOrder flag here
	// now the various scales
	// QED
	if(partnerType==ShowerPartnerType::QED) {
	assert(colourStructure()==ChargedChargedNeutral \|\|
	colourStructure()==ChargedNeutralCharged \|\|
	colourStructure()==NeutralChargedCharged );
	// normal case
	if(!bosonSplitting) {
	assert(colourStructure()==ChargedChargedNeutral \|\|
	colourStructure()==ChargedNeutralCharged );
	// set the scales
	// emitter
	emitter->scales().QED = zEmitter*scale;
	emitter->scales().QED_noAO = scale;
	emitter->scales().QCD_c = min(scale,parent->scales().QCD_c );
	emitter->scales().QCD_c_noAO = min(scale,parent->scales().QCD_c_noAO );
	emitter->scales().QCD_ac = min(scale,parent->scales().QCD_ac );
	emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO);
	- emitter->scales().EW = min(scale,parent->scales().EW );
	+ emitter->scales().EW_Z = min(scale,parent->scales().EW_Z );
	+ emitter->scales().EW_W = min(scale,parent->scales().EW_W );
	emitter->scales().EW_noAO = min(scale,parent->scales().EW_noAO );
	// emitted
	emitted->scales().QED = zEmitted*scale;
	emitted->scales().QED_noAO = scale;
	emitted->scales().QCD_c = ZERO;
	emitted->scales().QCD_c_noAO = ZERO;
	emitted->scales().QCD_ac = ZERO;
	emitted->scales().QCD_ac_noAO = ZERO;
	- emitted->scales().EW = min(scale,parent->scales().EW );
	+ emitted->scales().EW_Z = min(scale,parent->scales().EW_Z );
	+ emitted->scales().EW_W = min(scale,parent->scales().EW_W );
	emitted->scales().EW_noAO = min(scale,parent->scales().EW_noAO );
	}
	// gamma -> f fbar
	else {
	assert(colourStructure()==NeutralChargedCharged );
	// emitter
	emitter->scales().QED = zEmitter*scale;
	emitter->scales().QED_noAO = scale;
	if(hasColour(emitter)) {
	emitter->scales().QCD_c = zEmitter*scale;
	emitter->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(emitter)) {
	emitter->scales().QCD_ac = zEmitter*scale;
	emitter->scales().QCD_ac_noAO = scale;
	}
	- emitter->scales().EW = zEmitter*scale;
	+ emitter->scales().EW_Z = zEmitter*scale;
	+ emitter->scales().EW_W = zEmitter*scale;
	emitter->scales().EW_noAO = scale;
	// emitted
	emitted->scales().QED = zEmitted*scale;
	emitted->scales().QED_noAO = scale;
	if(hasColour(emitted)) {
	emitted->scales().QCD_c = zEmitted*scale;
	emitted->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(emitted)) {
	emitted->scales().QCD_ac = zEmitted*scale;
	emitted->scales().QCD_ac_noAO = scale;
	}
	- emitted->scales().EW = zEmitted*scale;
	+ emitted->scales().EW_Z = zEmitted*scale;
	+ emitted->scales().EW_W = zEmitted*scale;
	emitted->scales().EW_noAO = scale;
	}
	}
	// QCD
	else if (partnerType==ShowerPartnerType::QCDColourLine \|\|
	partnerType==ShowerPartnerType::QCDAntiColourLine) {
	// normal case eg q -> q g and g -> g g
	if(!bosonSplitting) {
	emitter->scales().QED = min(scale,parent->scales().QED );
	emitter->scales().QED_noAO = min(scale,parent->scales().QED_noAO);
	- emitter->scales().EW = min(scale,parent->scales().EW );
	+ emitter->scales().EW_Z = min(scale,parent->scales().EW_Z );
	+ emitter->scales().EW_W = min(scale,parent->scales().EW_W );
	emitter->scales().EW_noAO = min(scale,parent->scales().EW_noAO);
	if(partnerType==ShowerPartnerType::QCDColourLine) {
	emitter->scales().QCD_c = zEmitter*scale;
	emitter->scales().QCD_c_noAO = scale;
	emitter->scales().QCD_ac = min(zEmitter*scale,parent->scales().QCD_ac );
	emitter->scales().QCD_ac_noAO = min( scale,parent->scales().QCD_ac_noAO);
	}
	else {
	emitter->scales().QCD_c = min(zEmitter*scale,parent->scales().QCD_c );
	emitter->scales().QCD_c_noAO = min( scale,parent->scales().QCD_c_noAO );
	emitter->scales().QCD_ac = zEmitter*scale;
	emitter->scales().QCD_ac_noAO = scale;
	}
	// emitted
	emitted->scales().QED = ZERO;
	emitted->scales().QED_noAO = ZERO;
	emitted->scales().QCD_c = zEmitted*scale;
	emitted->scales().QCD_c_noAO = scale;
	emitted->scales().QCD_ac = zEmitted*scale;
	emitted->scales().QCD_ac_noAO = scale;
	- emitted->scales().EW = min(scale,parent->scales().EW );
	+ emitted->scales().EW_Z = min(scale,parent->scales().EW_Z );
	+ emitted->scales().EW_W = min(scale,parent->scales().EW_W );
	emitted->scales().EW_noAO = min(scale,parent->scales().EW_noAO);
	}
	// g -> q qbar
	else {
	// emitter
	if(emitter->dataPtr()->charged()) {
	emitter->scales().QED = zEmitter*scale;
	emitter->scales().QED_noAO = scale;
	}
	- emitter->scales().EW = zEmitter*scale;
	+ emitter->scales().EW_Z = zEmitter*scale;
	+ emitter->scales().EW_W = zEmitter*scale;
	emitter->scales().EW_noAO = scale;
	emitter->scales().QCD_c = zEmitter*scale;
	emitter->scales().QCD_c_noAO = scale;
	emitter->scales().QCD_ac = zEmitter*scale;
	emitter->scales().QCD_ac_noAO = scale;
	// emitted
	if(emitted->dataPtr()->charged()) {
	emitted->scales().QED = zEmitted*scale;
	emitted->scales().QED_noAO = scale;
	}
	- emitted->scales().EW = zEmitted*scale;
	+ emitted->scales().EW_Z = zEmitted*scale;
	+ emitted->scales().EW_W = zEmitted*scale;
	emitted->scales().EW_noAO = scale;
	emitted->scales().QCD_c = zEmitted*scale;
	emitted->scales().QCD_c_noAO = scale;
	emitted->scales().QCD_ac = zEmitted*scale;
	emitted->scales().QCD_ac_noAO = scale;
	}
	}
	else if(partnerType==ShowerPartnerType::EW) {
	// EW
	- emitter->scales().EW = zEmitter*scale;
	+ emitter->scales().EW_Z = zEmitter*scale;
	+ emitter->scales().EW_W = zEmitter*scale;
	emitter->scales().EW_noAO = scale;
	- emitted->scales().EW = zEmitted*scale;
	+ emitted->scales().EW_Z = zEmitted*scale;
	+ emitted->scales().EW_W = zEmitted*scale;
	emitted->scales().EW_noAO = scale;
	// QED
	// W radiation AO
	if(emitted->dataPtr()->charged()) {
	emitter->scales().QED = zEmitter*scale;
	emitter->scales().QED_noAO = scale;
	emitted->scales().QED = zEmitted*scale;
	emitted->scales().QED_noAO = scale;
	}
	// Z don't
	else {
	emitter->scales().QED = min(scale,parent->scales().QED );
	emitter->scales().QED_noAO = min(scale,parent->scales().QED_noAO);
	emitted->scales().QED = ZERO;
	emitted->scales().QED_noAO = ZERO;
	}
	// QCD
	emitter->scales().QCD_c = min(scale,parent->scales().QCD_c );
	emitter->scales().QCD_c_noAO = min(scale,parent->scales().QCD_c_noAO );
	emitter->scales().QCD_ac = min(scale,parent->scales().QCD_ac );
	emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO);
	emitted->scales().QCD_c = ZERO;
	emitted->scales().QCD_c_noAO = ZERO;
	emitted->scales().QCD_ac = ZERO;
	emitted->scales().QCD_ac_noAO = ZERO;
	}
	else
	assert(false);
	}

	void SplittingFunction::evaluateInitialStateScales(ShowerPartnerType partnerType,
	Energy scale, double z,
	tShowerParticlePtr parent,
	tShowerParticlePtr spacelike,
	tShowerParticlePtr timelike) {
	// scale for time-like child
	Energy AOScale = (1.-z)*scale;
	// QED
	if(partnerType==ShowerPartnerType::QED) {
	if(parent->id()==spacelike->id()) {
	// parent
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c );
	parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO );
	parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac );
	parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
	// timelike
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	timelike->scales().QCD_c = ZERO;
	timelike->scales().QCD_c_noAO = ZERO;
	timelike->scales().QCD_ac = ZERO;
	timelike->scales().QCD_ac_noAO = ZERO;
	}
	else if(parent->id()==timelike->id()) {
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	if(hasColour(parent)) {
	parent ->scales().QCD_c = scale;
	parent ->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(parent)) {
	parent ->scales().QCD_ac = scale;
	parent ->scales().QCD_ac_noAO = scale;
	}
	// timelike
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	if(hasColour(timelike)) {
	timelike->scales().QCD_c = AOScale;
	timelike->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(timelike)) {
	timelike->scales().QCD_ac = AOScale;
	timelike->scales().QCD_ac_noAO = scale;
	}
	}
	else {
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	parent ->scales().QCD_c = ZERO ;
	parent ->scales().QCD_c_noAO = ZERO ;
	parent ->scales().QCD_ac = ZERO ;
	parent ->scales().QCD_ac_noAO = ZERO ;
	// timelike
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	if(hasColour(timelike)) {
	timelike->scales().QCD_c = min(AOScale,spacelike->scales().QCD_ac );
	timelike->scales().QCD_c_noAO = min( scale,spacelike->scales().QCD_ac_noAO);
	}
	if(hasAntiColour(timelike)) {
	timelike->scales().QCD_ac = min(AOScale,spacelike->scales().QCD_c );
	timelike->scales().QCD_ac_noAO = min( scale,spacelike->scales().QCD_c_noAO );
	}
	}
	}
	// QCD
	else if (partnerType==ShowerPartnerType::QCDColourLine \|\|
	partnerType==ShowerPartnerType::QCDAntiColourLine) {
	// timelike
	if(timelike->dataPtr()->charged()) {
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	}
	if(hasColour(timelike)) {
	timelike->scales().QCD_c = AOScale;
	timelike->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(timelike)) {
	timelike->scales().QCD_ac = AOScale;
	timelike->scales().QCD_ac_noAO = scale;
	}
	if(parent->id()==spacelike->id()) {
	parent ->scales().QED = min(scale,spacelike->scales().QED );
	parent ->scales().QED_noAO = min(scale,spacelike->scales().QED_noAO );
	parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c );
	parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO );
	parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac );
	parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
	}
	else {
	if(parent->dataPtr()->charged()) {
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	}
	if(hasColour(parent)) {
	parent ->scales().QCD_c = scale;
	parent ->scales().QCD_c_noAO = scale;
	}
	if(hasAntiColour(parent)) {
	parent ->scales().QCD_ac = scale;
	parent ->scales().QCD_ac_noAO = scale;
	}
	}
	}
	else if(partnerType==ShowerPartnerType::EW) {
	if(abs(spacelike->id())!=ParticleID::Wplus &&
	spacelike->id() !=ParticleID::Z0 ) {
	// QCD scales
	parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c );
	parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO );
	parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac );
	parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
	timelike->scales().QCD_c = ZERO;
	timelike->scales().QCD_c_noAO = ZERO;
	timelike->scales().QCD_ac = ZERO;
	timelike->scales().QCD_ac_noAO = ZERO;
	// QED scales
	if(timelike->id()==ParticleID::Z0) {
	parent ->scales().QED = min(scale,spacelike->scales().QED );
	parent ->scales().QED_noAO = min(scale,spacelike->scales().QED_noAO );
	timelike->scales().QED = ZERO;
	timelike->scales().QED_noAO = ZERO;
	}
	else {
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	}
	// EW scales
	- parent ->scales().EW = scale;
	+ if(timelike->id()==ParticleID::Z0) {
	+ parent ->scales().EW_Z = scale;
	+ timelike->scales().EW_Z = AOScale;
	+ }
	+ if(timelike->id()==ParticleID::Wplus){
	+ parent ->scales().EW_W = scale;
	+ timelike->scales().EW_W = AOScale;
	+ }
	parent ->scales().EW_noAO = scale;
	- timelike->scales().EW = AOScale;
	timelike->scales().EW_noAO = scale;
	}
	else assert(false);
	}
	else
	assert(false);
	}

	void SplittingFunction::evaluateDecayScales(ShowerPartnerType partnerType,
	Energy scale, double z,
	tShowerParticlePtr parent,
	tShowerParticlePtr spacelike,
	tShowerParticlePtr timelike) {
	assert(parent->id()==spacelike->id());
	// angular-ordered scale for 2nd child
	Energy AOScale = (1.-z)*scale;
	// QED
	if(partnerType==ShowerPartnerType::QED) {
	// timelike
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	timelike->scales().QCD_c = ZERO;
	timelike->scales().QCD_c_noAO = ZERO;
	timelike->scales().QCD_ac = ZERO;
	timelike->scales().QCD_ac_noAO = ZERO;
	- timelike->scales().EW = ZERO;
	+ timelike->scales().EW_Z = ZERO;
	+ timelike->scales().EW_W = ZERO;
	timelike->scales().EW_noAO = ZERO;
	// spacelike
	spacelike->scales().QED = scale;
	spacelike->scales().QED_noAO = scale;
	- spacelike->scales().EW = max(scale,parent->scales().EW );
	+ spacelike->scales().EW_Z = max(scale,parent->scales().EW_Z );
	+ spacelike->scales().EW_W = max(scale,parent->scales().EW_W );
	spacelike->scales().EW_noAO = max(scale,parent->scales().EW_noAO );
	}
	// QCD
	else if(partnerType==ShowerPartnerType::QCDColourLine \|\|
	partnerType==ShowerPartnerType::QCDAntiColourLine) {
	// timelike
	timelike->scales().QED = ZERO;
	timelike->scales().QED_noAO = ZERO;
	timelike->scales().QCD_c = AOScale;
	timelike->scales().QCD_c_noAO = scale;
	timelike->scales().QCD_ac = AOScale;
	timelike->scales().QCD_ac_noAO = scale;
	- timelike->scales().EW = ZERO;
	+ timelike->scales().EW_Z = ZERO;
	+ timelike->scales().EW_W = ZERO;
	timelike->scales().EW_noAO = ZERO;
	// spacelike
	spacelike->scales().QED = max(scale,parent->scales().QED );
	spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO );
	- spacelike->scales().EW = max(scale,parent->scales().EW );
	+ spacelike->scales().EW_Z = max(scale,parent->scales().EW_Z );
	+ spacelike->scales().EW_W = max(scale,parent->scales().EW_W );
	spacelike->scales().EW_noAO = max(scale,parent->scales().EW_noAO );
	}
	else if(partnerType==ShowerPartnerType::EW) {
	// EW
	- timelike->scales().EW = AOScale;
	+ timelike->scales().EW_Z = AOScale;
	+ timelike->scales().EW_W = AOScale;
	timelike->scales().EW_noAO = scale;
	- spacelike->scales().EW = max(scale,parent->scales().EW );
	+ spacelike->scales().EW_Z = max(scale,parent->scales().EW_Z );
	+ spacelike->scales().EW_W = max(scale,parent->scales().EW_W );
	spacelike->scales().EW_noAO = max(scale,parent->scales().EW_noAO );
	// QCD
	timelike->scales().QCD_c = ZERO;
	timelike->scales().QCD_c_noAO = ZERO;
	timelike->scales().QCD_ac = ZERO;
	timelike->scales().QCD_ac_noAO = ZERO;
	- timelike->scales().EW = ZERO;
	+ timelike->scales().EW_Z = ZERO;
	+ timelike->scales().EW_W = ZERO;
	timelike->scales().EW_noAO = ZERO;
	// QED
	timelike->scales().QED = ZERO;
	timelike->scales().QED_noAO = ZERO;
	spacelike->scales().QED = max(scale,parent->scales().QED );
	spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO );
	}
	else
	assert(false);
	spacelike->scales().QCD_c = max(scale,parent->scales().QCD_c );
	spacelike->scales().QCD_c_noAO = max(scale,parent->scales().QCD_c_noAO );
	spacelike->scales().QCD_ac = max(scale,parent->scales().QCD_ac );
	spacelike->scales().QCD_ac_noAO = max(scale,parent->scales().QCD_ac_noAO);
	}
	diff --git a/Shower/QTilde/SplittingFunctions/SplittingGenerator.cc b/Shower/QTilde/SplittingFunctions/SplittingGenerator.cc
	--- a/Shower/QTilde/SplittingFunctions/SplittingGenerator.cc
	+++ b/Shower/QTilde/SplittingFunctions/SplittingGenerator.cc
	@@ -1,647 +1,696 @@
	// -- C++ --
	//
	// SplittingGenerator.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 SplittingGenerator class.
	//

	#include "SplittingGenerator.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Command.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Utilities/StringUtils.h"
	#include "ThePEG/Repository/Repository.h"
	#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
	#include "Herwig/Shower/ShowerHandler.h"
	#include "ThePEG/Utilities/Rebinder.h"
	#include <cassert>
	#include "ThePEG/Utilities/DescribeClass.h"

	using namespace Herwig;

	namespace {

	bool checkInteraction(ShowerInteraction allowed,
	ShowerInteraction splitting) {
	if(allowed==ShowerInteraction::ALL)
	return true;
	else if(allowed==ShowerInteraction::QEDQCD &&
	(splitting==ShowerInteraction::QED \|\|
	splitting==ShowerInteraction::QCD ))
	return true;
	else if(allowed == splitting)
	return true;
	else
	return false;
	}

	}

	DescribeClass<SplittingGenerator,Interfaced>
	describeSplittingGenerator ("Herwig::SplittingGenerator","");

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

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

	void SplittingGenerator::persistentOutput(PersistentOStream & os) const {
	os << _bbranchings << _fbranchings << _deTuning;
	}

	void SplittingGenerator::persistentInput(PersistentIStream & is, int) {
	is >> _bbranchings >> _fbranchings >> _deTuning;
	}

	void SplittingGenerator::Init() {

	static ClassDocumentation<SplittingGenerator> documentation
	("There class is responsible for initializing the Sudakov form factors ",
	"and generating splittings.");

	static Command<SplittingGenerator> interfaceAddSplitting
	("AddFinalSplitting",
	"Adds another splitting to the list of splittings considered "
	"in the shower. Command is a->b,c; Sudakov",
	&SplittingGenerator::addFinalSplitting);

	static Command<SplittingGenerator> interfaceAddInitialSplitting
	("AddInitialSplitting",
	"Adds another splitting to the list of initial splittings to consider "
	"in the shower. Command is a->b,c; Sudakov. Here the particle a is the "
	"particle that is PRODUCED by the splitting. b is the initial state "
	"particle that is splitting in the shower.",
	&SplittingGenerator::addInitialSplitting);

	static Command<SplittingGenerator> interfaceDeleteSplitting
	("DeleteFinalSplitting",
	"Deletes a splitting from the list of splittings considered "
	"in the shower. Command is a->b,c; Sudakov",
	&SplittingGenerator::deleteFinalSplitting);

	static Command<SplittingGenerator> interfaceDeleteInitialSplitting
	("DeleteInitialSplitting",
	"Deletes a splitting from the list of initial splittings to consider "
	"in the shower. Command is a->b,c; Sudakov. Here the particle a is the "
	"particle that is PRODUCED by the splitting. b is the initial state "
	"particle that is splitting in the shower.",
	&SplittingGenerator::deleteInitialSplitting);

	static Parameter<SplittingGenerator,double> interfaceDetuning
	("Detuning",
	"The Detuning parameter to make the veto algorithm less efficient to improve the weight variations",
	&SplittingGenerator::_deTuning, 1.0, 1.0, 10.0,
	false, false, Interface::limited);

	}

	string SplittingGenerator::addSplitting(string arg, bool final) {
	string partons = StringUtils::car(arg);
	string sudakov = StringUtils::cdr(arg);
	vector<tPDPtr> products;
	string::size_type next = partons.find("->");
	if(next == string::npos)
	return "Error: Invalid string for splitting " + arg;
	if(partons.find(';') == string::npos)
	return "Error: Invalid string for splitting " + arg;
	tPDPtr parent = Repository::findParticle(partons.substr(0,next));
	partons = partons.substr(next+2);
	do {
	next = min(partons.find(','), partons.find(';'));
	tPDPtr pdp = Repository::findParticle(partons.substr(0,next));
	partons = partons.substr(next+1);
	if(pdp) products.push_back(pdp);
	else return "Error: Could not create splitting from " + arg;
	} while(partons[0] != ';' && partons.size());
	SudakovPtr s;
	s = dynamic_ptr_cast<SudakovPtr>(Repository::TraceObject(sudakov));
	if(!s) return "Error: Could not load Sudakov " + sudakov + '\n';
	IdList ids;
	ids.push_back(parent);
	for(vector<tPDPtr>::iterator it = products.begin(); it!=products.end(); ++it)
	ids.push_back(*it);
	// check splitting can handle this
	if(!s->splittingFn()->accept(ids))
	return "Error: Sudakov " + sudakov + " SplittingFunction can't handle particles\n";
	// add to map
	addToMap(ids,s,final);
	return "";
	}

	string SplittingGenerator::deleteSplitting(string arg, bool final) {
	string partons = StringUtils::car(arg);
	string sudakov = StringUtils::cdr(arg);
	vector<tPDPtr> products;
	string::size_type next = partons.find("->");
	if(next == string::npos)
	return "Error: Invalid string for splitting " + arg;
	if(partons.find(';') == string::npos)
	return "Error: Invalid string for splitting " + arg;
	tPDPtr parent = Repository::findParticle(partons.substr(0,next));
	partons = partons.substr(next+2);
	do {
	next = min(partons.find(','), partons.find(';'));
	tPDPtr pdp = Repository::findParticle(partons.substr(0,next));
	partons = partons.substr(next+1);
	if(pdp) products.push_back(pdp);
	else return "Error: Could not create splitting from " + arg;
	} while(partons[0] != ';' && partons.size());
	SudakovPtr s;
	s = dynamic_ptr_cast<SudakovPtr>(Repository::TraceObject(sudakov));
	if(!s) return "Error: Could not load Sudakov " + sudakov + '\n';
	IdList ids;
	ids.push_back(parent);
	for(vector<tPDPtr>::iterator it = products.begin(); it!=products.end(); ++it)
	ids.push_back(*it);
	// check splitting can handle this
	if(!s->splittingFn()->accept(ids))
	return "Error: Sudakov " + sudakov + " SplittingFunction can't handle particles\n";
	// delete from map
	deleteFromMap(ids,s,final);
	return "";
	}

	void SplittingGenerator::addToMap(const IdList &ids, const SudakovPtr &s, bool final) {
	if(!final) {
	// search if the branching was already included.
	auto binsert =BranchingInsert(abs(ids[1]->id()),BranchingElement(s,ids));
	// get the range of already inserted splittings.
	auto eqrange=_bbranchings.equal_range(binsert.first);
	for(auto it = eqrange.first; it != eqrange.second; ++it){
	if((*it).second == binsert.second)
	throw Exception()<<"SplittingGenerator: Trying to insert existing splitting.\n"
	<< Exception::setuperror;
	}
	_bbranchings.insert(binsert);
	s->addSplitting(ids);
	}
	else {
	// search if the branching was already included.
	auto binsert =BranchingInsert(abs(ids[0]->id()),BranchingElement(s,ids));
	// get the range of already inserted splittings.
	auto eqrange=_fbranchings.equal_range(binsert.first);
	for(auto it = eqrange.first; it != eqrange.second; ++it){
	if((*it).second ==binsert.second)
	throw Exception()<<"SplittingGenerator: Trying to insert existing splitting.\n"
	<< Exception::setuperror;
	}

	_fbranchings.insert(binsert);
	s->addSplitting(ids);
	}
	}

	void SplittingGenerator::deleteFromMap(const IdList &ids,
	const SudakovPtr &s, bool final) {
	bool didRemove=false;
	if(!final) {
	pair<BranchingList::iterator,BranchingList::iterator>
	range = _bbranchings.equal_range(abs(ids[1]->id()));
	for(BranchingList::iterator it=range.first;
	it!=range.second&&it!=_bbranchings.end();++it) {
	if(it->second.sudakov==s&&it->second.particles==ids) {
	BranchingList::iterator it2=it;
	--it;
	_bbranchings.erase(it2);
	didRemove=true;
	}
	}
	s->removeSplitting(ids);
	}
	else {
	pair<BranchingList::iterator,BranchingList::iterator>
	range = _fbranchings.equal_range(abs(ids[0]->id()));
	for(BranchingList::iterator it=range.first;
	it!=range.second&&it!=_fbranchings.end();++it) {
	if(it->second.sudakov==s&&it->second.particles==ids) {
	BranchingList::iterator it2 = it;
	--it;
	_fbranchings.erase(it2);
	didRemove=true;
	}
	}
	s->removeSplitting(ids);
	}
	if (!didRemove)
	throw Exception()<<"SplittingGenerator: Try to remove non existing splitting.\n"
	<< Exception::setuperror;
	}

	Branching SplittingGenerator::chooseForwardBranching(ShowerParticle &particle,
	double enhance,
	ShowerInteraction type) const {
	RhoDMatrix rho;
	bool rhoCalc(false);
	Energy newQ = ZERO;
	ShoKinPtr kinematics = ShoKinPtr();
	ShowerPartnerType partnerType(ShowerPartnerType::Undefined);
	SudakovPtr sudakov = SudakovPtr();
	IdList ids;
	// First, find the eventual branching, corresponding to the highest scale.
	long index = abs(particle.data().id());
	// if no branchings return empty branching struct
	if( _fbranchings.find(index) == _fbranchings.end() )
	return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined);
	// otherwise select branching
	for(BranchingList::const_iterator cit = _fbranchings.lower_bound(index);
	cit != _fbranchings.upper_bound(index); ++cit) {
	// check either right interaction or doing both
	if(!checkInteraction(type,cit->second.sudakov->interactionType())) continue;
	if(!rhoCalc) {
	rho = particle.extractRhoMatrix(true);
	rhoCalc = true;
	}
	// whether or not this interaction should be angular ordered
	bool angularOrdered = cit->second.sudakov->splittingFn()->angularOrdered();
	ShoKinPtr newKin;
	ShowerPartnerType type;
	IdList particles = particle.id()!=cit->first ? cit->second.conjugateParticles : cit->second.particles;
	// work out which starting scale we need
	if(cit->second.sudakov->interactionType()==ShowerInteraction::QED) {
	type = ShowerPartnerType::QED;
	Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
	newKin = cit->second.sudakov->
	generateNextTimeBranching(startingScale,particles,rho,enhance,
	_deTuning);
	}
	else if(cit->second.sudakov->interactionType()==ShowerInteraction::QCD) {
	// special for octets
	if(particle.dataPtr()->iColour()==PDT::Colour8) {
	// octet -> octet octet
	if(cit->second.sudakov->splittingFn()->colourStructure()==OctetOctetOctet) {
	type = ShowerPartnerType::QCDColourLine;
	Energy startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
	newKin= cit->second.sudakov->
	generateNextTimeBranching(startingScale,particles,rho,0.5*enhance,
	_deTuning);
	startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
	ShoKinPtr newKin2 = cit->second.sudakov->
	generateNextTimeBranching(startingScale,particles,rho,0.5*enhance,
	_deTuning);
	// pick the one with the highest scale
	if( ( newKin && newKin2 && newKin2->scale() > newKin->scale()) \|\|
	(!newKin && newKin2) ) {
	newKin = newKin2;
	type = ShowerPartnerType::QCDAntiColourLine;
	}
	}
	// other g -> q qbar
	else {
	Energy startingScale = angularOrdered ?
	max(particle.scales().QCD_c , particle.scales().QCD_ac ) :
	max(particle.scales().QCD_c_noAO, particle.scales().QCD_ac_noAO);
	newKin= cit->second.sudakov->
	generateNextTimeBranching(startingScale,particles,rho,enhance,
	_deTuning);
	type = UseRandom::rndbool() ?
	ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
	}
	}
	// everything else q-> qg etc
	else {
	Energy startingScale;
	if(particle.hasColour()) {
	type = ShowerPartnerType::QCDColourLine;
	startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
	}
	else {
	type = ShowerPartnerType::QCDAntiColourLine;
	startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
	}
	newKin= cit->second.sudakov->
	generateNextTimeBranching(startingScale,particles,rho,enhance,
	_deTuning);
	}
	}
	else if(cit->second.sudakov->interactionType()==ShowerInteraction::EW) {
	type = ShowerPartnerType::EW;
	- Energy startingScale = particle.scales().EW;
	+ // -----------------------------------------------------------------//
	+ Energy startingScale;
	+ if(particles[1]->id()==ParticleID::Z0 \|\| particles[2]->id()==ParticleID::Z0){
	+ startingScale = particle.scales().EW_Z; //ref: symmetry choice
	+ //cerr<<"-->HINT: Z0\n";
	+ }
	+ else if(abs(particles[1]->id())==ParticleID::Wplus
	+ \|\| abs(particles[2]->id())==ParticleID::Wplus){
	+ startingScale = particle.scales().EW_W; //ref: anticolor choice
	+ //cerr<<"-->HINT: W+/W-\n";
	+ }
	+ else{ //ensurance!
	+ startingScale = particle.scales().EW_Z; //ref: symmetry choice
	+ cerr<<"suspicious activity in EW shower: P[0]: "
	+ <<particles[0]->id()<<" P[1]: "<<particles[1]->id()<<" P[2]:"<<particles[2]->id()<<"\n";
	+ }
	+ // -----------------------------------------------------------------//
	newKin = cit->second.sudakov->
	generateNextTimeBranching(startingScale,particles,rho,enhance,_deTuning);
	}
	// shouldn't be anything else
	else
	assert(false);
	// if no kinematics contine
	if(!newKin) continue;
	// select highest scale
	if( newKin->scale() > newQ ) {
	kinematics = newKin;
	newQ = newKin->scale();
	ids = particles;
	sudakov = cit->second.sudakov;
	partnerType = type;
	}
	}
	// return empty branching if nothing happened
	if(!kinematics) {
	if ( particle.spinInfo() ) particle.spinInfo()->undecay();
	return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
	ShowerPartnerType::Undefined);
	}
	// if not hard generate phi
	kinematics->phi(sudakov->generatePhiForward(particle,ids,kinematics,rho));
	// and return it
	return Branching(kinematics, ids,sudakov,partnerType);
	}

	Branching SplittingGenerator::
	chooseDecayBranching(ShowerParticle &particle,
	const ShowerParticle::EvolutionScales & stoppingScales,
	Energy minmass, double enhance,
	ShowerInteraction interaction) const {
	RhoDMatrix rho(particle.dataPtr()->iSpin());
	Energy newQ = Constants::MaxEnergy;
	ShoKinPtr kinematics;
	SudakovPtr sudakov;
	ShowerPartnerType partnerType(ShowerPartnerType::Undefined);
	IdList ids;
	// First, find the eventual branching, corresponding to the lowest scale.
	long index = abs(particle.data().id());
	// if no branchings return empty branching struct
	if(_fbranchings.find(index) == _fbranchings.end())
	return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined);
	// otherwise select branching
	for(BranchingList::const_iterator cit = _fbranchings.lower_bound(index);
	cit != _fbranchings.upper_bound(index); ++cit) {
	// check interaction doesn't change flavour
	if(cit->second.particles[1]->id()!=index&&cit->second.particles[2]->id()!=index) continue;
	// check either right interaction or doing both
	if(!checkInteraction(interaction,cit->second.sudakov->interactionType())) continue;
	// whether or not this interaction should be angular ordered
	bool angularOrdered = cit->second.sudakov->splittingFn()->angularOrdered();
	ShoKinPtr newKin;
	IdList particles = particle.id()!=cit->first ? cit->second.conjugateParticles : cit->second.particles;
	ShowerPartnerType type;
	// work out which starting scale we need
	if(cit->second.sudakov->interactionType()==ShowerInteraction::QED) {
	type = ShowerPartnerType::QED;
	Energy stoppingScale = angularOrdered ? stoppingScales.QED : stoppingScales.QED_noAO;
	Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
	if(startingScale < stoppingScale ) {
	newKin = cit->second.sudakov->
	generateNextDecayBranching(startingScale,stoppingScale,minmass,particles,rho,
	enhance,_deTuning);
	}
	}
	else if(cit->second.sudakov->interactionType()==ShowerInteraction::QCD) {
	// special for octets
	if(particle.dataPtr()->iColour()==PDT::Colour8) {
	// octet -> octet octet
	if(cit->second.sudakov->splittingFn()->colourStructure()==OctetOctetOctet) {
	Energy stoppingColour = angularOrdered ? stoppingScales.QCD_c : stoppingScales.QCD_c_noAO;
	Energy stoppingAnti = angularOrdered ? stoppingScales.QCD_ac : stoppingScales.QCD_ac_noAO;
	Energy startingColour = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
	Energy startingAnti = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
	type = ShowerPartnerType::QCDColourLine;
	if(startingColour<stoppingColour) {
	newKin= cit->second.sudakov->
	generateNextDecayBranching(startingColour,stoppingColour,minmass,
	particles,rho,0.5*enhance,_deTuning);
	}
	ShoKinPtr newKin2;
	if(startingAnti<stoppingAnti) {
	newKin2 = cit->second.sudakov->
	generateNextDecayBranching(startingAnti,stoppingAnti,minmass,
	particles,rho,0.5*enhance,_deTuning);
	}
	// pick the one with the lowest scale
	if( (newKin&&newKin2&&newKin2->scale()<newKin->scale()) \|\|
	(!newKin&&newKin2) ) {
	newKin = newKin2;
	type = ShowerPartnerType::QCDAntiColourLine;
	}
	}
	// other
	else {
	assert(false);
	}
	}
	// everything else
	else {
	Energy startingScale,stoppingScale;
	if(particle.hasColour()) {
	type = ShowerPartnerType::QCDColourLine;
	stoppingScale = angularOrdered ? stoppingScales.QCD_c : stoppingScales.QCD_c_noAO;
	startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
	}
	else {
	type = ShowerPartnerType::QCDAntiColourLine;
	stoppingScale = angularOrdered ? stoppingScales.QCD_ac : stoppingScales.QCD_ac_noAO;
	startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
	}
	if(startingScale < stoppingScale ) {
	newKin = cit->second.sudakov->
	generateNextDecayBranching(startingScale,stoppingScale,minmass,particles,rho,
	enhance,_deTuning);
	}
	}
	}
	else if(cit->second.sudakov->interactionType()==ShowerInteraction::EW) {
	type = ShowerPartnerType::EW;
	- Energy stoppingScale = stoppingScales.EW;
	- Energy startingScale = particle.scales().EW;
	+ // -----------------------------------------------------------------//
	+ Energy stoppingScale, startingScale;
	+ if(particles[1]->id()==ParticleID::Z0 \|\| particles[2]->id()==ParticleID::Z0){
	+ stoppingScale = stoppingScales.EW_Z;
	+ startingScale = particle.scales().EW_Z; //ref: symmetry choice
	+ //cerr<<"-->HINT: Z0\n";
	+ }
	+ else if(abs(particles[1]->id())==ParticleID::Wplus \|\| abs(particles[2]->id())==ParticleID::Wplus){
	+ stoppingScale = stoppingScales.EW_W;
	+ startingScale = particle.scales().EW_W; //ref: anticolor choice
	+ //cerr<<"-->HINT: W+/W-\n";
	+ }
	+ else{ //ensurance!
	+ stoppingScale = stoppingScales.EW_Z;
	+ startingScale = particle.scales().EW_Z; //ref: symmetry choice
	+ cerr<<"suspicious activity in EW shower: P[0]: "
	+ <<particles[0]->id()<<" P[1]: "<<particles[1]->id()<<" P[2]:"<<particles[2]->id()<<"\n";
	+ }
	+ // -----------------------------------------------------------------//
	if(startingScale < stoppingScale ) {
	newKin = cit->second.sudakov->
	generateNextDecayBranching(startingScale,stoppingScale,minmass,particles,rho,enhance,_deTuning);
	}
	}
	// shouldn't be anything else
	else
	assert(false);
	if(!newKin) continue;
	// select highest scale
	if(newKin->scale() < newQ ) {
	newQ = newKin->scale();
	ids = particles;
	kinematics=newKin;
	sudakov=cit->second.sudakov;
	partnerType = type;
	}
	}
	// return empty branching if nothing happened
	if(!kinematics) return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
	ShowerPartnerType::Undefined);
	// and generate phi
	kinematics->phi(sudakov->generatePhiDecay(particle,ids,kinematics,rho));
	// and return it
	return Branching(kinematics, ids,sudakov,partnerType);
	}

	Branching SplittingGenerator::
	chooseBackwardBranching(ShowerParticle &particle,PPtr,
	double enhance,
	Ptr<BeamParticleData>::transient_const_pointer beam,
	ShowerInteraction type,
	tcPDFPtr pdf, Energy freeze) const {
	RhoDMatrix rho;
	bool rhoCalc(false);
	Energy newQ=ZERO;
	ShoKinPtr kinematics=ShoKinPtr();
	ShowerPartnerType partnerType(ShowerPartnerType::Undefined);
	SudakovPtr sudakov;
	IdList ids;
	// First, find the eventual branching, corresponding to the highest scale.
	long index = abs(particle.id());
	// if no possible branching return
	if(_bbranchings.find(index) == _bbranchings.end())
	return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined);
	// otherwise select branching
	for(BranchingList::const_iterator cit = _bbranchings.lower_bound(index);
	cit != _bbranchings.upper_bound(index); ++cit ) {
	// check either right interaction or doing both
	if(!checkInteraction(type,cit->second.sudakov->interactionType())) continue;
	// setup the PDF
	cit->second.sudakov->setPDF(pdf,freeze);
	//calc rho as needed
	if(!rhoCalc) {
	rho = particle.extractRhoMatrix(false);
	rhoCalc = true;
	}
	// whether or not this interaction should be angular ordered
	bool angularOrdered = cit->second.sudakov->splittingFn()->angularOrdered();
	ShoKinPtr newKin;
	IdList particles = particle.id()!=cit->first ? cit->second.conjugateParticles : cit->second.particles;
	ShowerPartnerType type;
	if(cit->second.sudakov->interactionType()==ShowerInteraction::QED) {
	type = ShowerPartnerType::QED;
	Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
	newKin=cit->second.sudakov->
	generateNextSpaceBranching(startingScale,particles, particle.x(),rho,enhance,
	beam,_deTuning);
	}
	else if(cit->second.sudakov->interactionType()==ShowerInteraction::QCD) {
	// special for octets
	if(particle.dataPtr()->iColour()==PDT::Colour8) {
	// octet -> octet octet
	if(cit->second.sudakov->splittingFn()->colourStructure()==OctetOctetOctet) {
	type = ShowerPartnerType::QCDColourLine;
	Energy startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
	newKin = cit->second.sudakov->
	generateNextSpaceBranching(startingScale,particles, particle.x(),rho,0.5*enhance,
	beam,_deTuning);
	startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
	ShoKinPtr newKin2 = cit->second.sudakov->
	generateNextSpaceBranching(startingScale,particles, particle.x(),rho,
	0.5*enhance,beam,_deTuning);
	// pick the one with the highest scale
	if( (newKin&&newKin2&&newKin2->scale()>newKin->scale()) \|\|
	(!newKin&&newKin2) ) {
	newKin = newKin2;
	type = ShowerPartnerType::QCDAntiColourLine;
	}
	}
	else {
	Energy startingScale = angularOrdered ?
	max(particle.scales().QCD_c , particle.scales().QCD_ac ) :
	max(particle.scales().QCD_c_noAO, particle.scales().QCD_ac_noAO);
	type = UseRandom::rndbool() ?
	ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
	newKin=cit->second.sudakov->
	generateNextSpaceBranching(startingScale,particles, particle.x(),rho,enhance,beam,_deTuning);
	}
	}
	// everything else
	else {
	Energy startingScale;
	if(particle.hasColour()) {
	type = ShowerPartnerType::QCDColourLine;
	startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
	}
	else {
	type = ShowerPartnerType::QCDAntiColourLine;
	startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
	}
	newKin=cit->second.sudakov->
	generateNextSpaceBranching(startingScale,particles,particle.x(),rho,enhance,beam,_deTuning);
	}
	}
	else if(cit->second.sudakov->interactionType()==ShowerInteraction::EW) {
	type = ShowerPartnerType::EW;
	- Energy startingScale = particle.scales().EW;
	+ // -----------------------------------------------------------------//
	+ Energy startingScale;
	+ if(particles[1]->id()==ParticleID::Z0 \|\| particles[2]->id()==ParticleID::Z0){
	+ startingScale = particle.scales().EW_Z; //ref: symmetry choice
	+ //cerr<<"-->HINT: Z0\n";
	+ }
	+ else if(abs(particles[1]->id())==ParticleID::Wplus
	+ \|\| abs(particles[2]->id())==ParticleID::Wplus){
	+ startingScale = particle.scales().EW_W; //ref: anticolor choice
	+ //cerr<<"-->HINT: W+/W-\n";
	+ }
	+ else{ //ensurance!
	+ startingScale = particle.scales().EW_Z; //ref: symmetry choice
	+ cerr<<"suspicious activity in EW shower: P[0]: "
	+ <<particles[0]->id()<<" P[1]: "<<particles[1]->id()<<" P[2]:"<<particles[2]->id()<<"\n";
	+ }
	+ // -----------------------------------------------------------------//
	newKin=cit->second.sudakov->
	generateNextSpaceBranching(startingScale,particles,particle.x(),rho,enhance,beam,_deTuning);
	}
	// shouldn't be anything else
	else
	assert(false);
	// if no kinematics contine
	if(!newKin) continue;
	// select highest scale
	if(newKin->scale() > newQ) {
	newQ = newKin->scale();
	kinematics=newKin;
	ids = particles;
	sudakov=cit->second.sudakov;
	partnerType = type;
	}
	}
	// return empty branching if nothing happened
	if(!kinematics) {
	if ( particle.spinInfo() ) particle.spinInfo()->undecay();
	return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
	ShowerPartnerType::Undefined);
	}
	// initialize the ShowerKinematics
	// and generate phi
	kinematics->phi(sudakov->generatePhiBackward(particle,ids,kinematics,rho));
	// return the answer
	return Branching(kinematics, ids,sudakov,partnerType);
	}

	void SplittingGenerator::rebind(const TranslationMap & trans) {
	BranchingList::iterator cit;
	for(cit=_fbranchings.begin();cit!=_fbranchings.end();++cit) {
	(cit->second).sudakov=trans.translate((cit->second).sudakov);
	for(unsigned int ix=0;ix<(cit->second).particles.size();++ix) {
	(cit->second).particles[ix]=trans.translate((cit->second).particles[ix]);
	}
	for(unsigned int ix=0;ix<(cit->second).conjugateParticles.size();++ix) {
	(cit->second).conjugateParticles[ix]=trans.translate((cit->second).conjugateParticles[ix]);
	}
	}
	for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit) {
	(cit->second).sudakov=trans.translate((cit->second).sudakov);
	for(unsigned int ix=0;ix<(cit->second).particles.size();++ix) {
	(cit->second).particles[ix]=trans.translate((cit->second).particles[ix]);
	}
	for(unsigned int ix=0;ix<(cit->second).conjugateParticles.size();++ix) {
	(cit->second).conjugateParticles[ix]=trans.translate((cit->second).conjugateParticles[ix]);
	}
	}
	Interfaced::rebind(trans);
	}

	IVector SplittingGenerator::getReferences() {
	IVector ret = Interfaced::getReferences();
	BranchingList::iterator cit;
	for(cit=_fbranchings.begin();cit!=_fbranchings.end();++cit) {
	ret.push_back((cit->second).sudakov);
	for(unsigned int ix=0;ix<(cit->second).particles.size();++ix)
	ret.push_back(const_ptr_cast<tPDPtr>((cit->second).particles[ix]));
	for(unsigned int ix=0;ix<(cit->second).conjugateParticles.size();++ix)
	ret.push_back(const_ptr_cast<tPDPtr>((cit->second).conjugateParticles[ix]));
	}
	for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit) {
	ret.push_back((cit->second).sudakov);
	for(unsigned int ix=0;ix<(cit->second).particles.size();++ix)
	ret.push_back(const_ptr_cast<tPDPtr>((cit->second).particles[ix]));
	for(unsigned int ix=0;ix<(cit->second).conjugateParticles.size();++ix)
	ret.push_back(const_ptr_cast<tPDPtr>((cit->second).conjugateParticles[ix]));
	}
	return ret;
	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 3:44 PM (1 d, 14 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4022654
Default Alt Text: (282 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions