No OneTemporary
Actions

Size

245 KB

Subscribers

None

View Options

	diff --git a/Decay/Perturbative/SMTopDecayer.cc b/Decay/Perturbative/SMTopDecayer.cc
	--- a/Decay/Perturbative/SMTopDecayer.cc
	+++ b/Decay/Perturbative/SMTopDecayer.cc
	@@ -1,1202 +1,1202 @@
	// -- C++ --
	//
	// SMTopDecayer.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the SMTopDecayer class.
	//

	#include "SMTopDecayer.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "Herwig++/Decay/DecayVertex.h"
	#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	#include "Herwig++/PDT/ThreeBodyAllOn1IntegralCalculator.h"

	using namespace Herwig;
	using namespace ThePEG::Helicity;

	SMTopDecayer::SMTopDecayer()
	: _wquarkwgt(6,0.),_wleptonwgt(3,0.), _xg_sampling(1.5),
	_initialenhance(1.), _finalenhance(2.3), _useMEforT2(true) {
	_wleptonwgt[0] = 0.302583;
	_wleptonwgt[1] = 0.301024;
	_wleptonwgt[2] = 0.299548;
	_wquarkwgt[0] = 0.851719;
	_wquarkwgt[1] = 0.0450162;
	_wquarkwgt[2] = 0.0456962;
	_wquarkwgt[3] = 0.859839;
	_wquarkwgt[4] = 3.9704e-06;
	_wquarkwgt[5] = 0.000489657;
	generateIntermediates(true);
	}

	bool SMTopDecayer::accept(tcPDPtr parent, const tPDVector & children) const {
	if(abs(parent->id()) != ParticleID::t) return false;
	int id0(0),id1(0),id2(0);
	for(tPDVector::const_iterator it = children.begin();
	it != children.end();++it) {
	int id=(**it).id(),absid(abs(id));
	if(absid==ParticleID::b&&double(id)/double(parent->id())>0) {
	id0=id;
	}
	else {
	switch (absid) {
	case ParticleID::nu_e:
	case ParticleID::nu_mu:
	case ParticleID::nu_tau:
	id1 = id;
	break;
	case ParticleID::eminus:
	case ParticleID::muminus:
	case ParticleID::tauminus:
	id2 = id;
	break;
	case ParticleID::b:
	case ParticleID::d:
	case ParticleID::s:
	id1 = id;
	break;
	case ParticleID::u:
	case ParticleID::c:
	id2=id;
	break;
	default :
	break;
	}
	}
	}
	if(id0==0\|\|id1==0\|\|id2==0) return false;
	if(double(id1)/double(id2)>0) return false;
	return true;
	}

	ParticleVector SMTopDecayer::decay(const Particle & parent,
	const tPDVector & children) const {
	int id1(0),id2(0);
	for(tPDVector::const_iterator it = children.begin();
	it != children.end();++it) {
	int id=(**it).id(),absid=abs(id);
	if(absid == ParticleID::b && double(id)/double(parent.id())>0) continue;
	//leptons
	if(absid > 10 && absid%2==0) id1=absid;
	if(absid > 10 && absid%2==1) id2=absid;
	//quarks
	if(absid < 10 && absid%2==0) id2=absid;
	if(absid < 10 && absid%2==1) id1=absid;
	}
	unsigned int imode(0);
	if(id2 >=11 && id2<=16) imode = (id1-12)/2;
	else imode = id1+1+id2/2;
	bool cc = parent.id() == ParticleID::tbar;
	ParticleVector out(generate(true,cc,imode,parent));
	//arrange colour flow
	PPtr pparent=const_ptr_cast<PPtr>(&parent);
	out[1]->incomingColour(pparent,out[1]->id()<0);
	ParticleVector products = out[0]->children();
	if(products[0]->hasColour())
	products[0]->colourNeighbour(products[1],true);
	else if(products[0]->hasAntiColour())
	products[0]->colourNeighbour(products[1],false);
	return out;
	}

	void SMTopDecayer::persistentOutput(PersistentOStream & os) const {
	os << _wvertex << _wquarkwgt << _wleptonwgt << _wplus << _alpha
	<< _initialenhance << _finalenhance << _xg_sampling << _useMEforT2;
	}

	void SMTopDecayer::persistentInput(PersistentIStream & is, int) {
	is >> _wvertex >> _wquarkwgt >> _wleptonwgt >> _wplus >> _alpha
	>> _initialenhance >> _finalenhance >> _xg_sampling >> _useMEforT2;
	}

	ClassDescription<SMTopDecayer> SMTopDecayer::initSMTopDecayer;
	// Definition of the static class description member.

	void SMTopDecayer::Init() {

	static ClassDocumentation<SMTopDecayer> documentation
	("This is the implementation of the SMTopDecayer which "
	"decays top quarks into bottom quarks and either leptons "
	"or quark-antiquark pairs including the matrix element for top decay",
	"The matrix element correction for top decay \\cite{Hamilton:2006ms}.",
	"%\\cite{Hamilton:2006ms}\n"
	"\\bibitem{Hamilton:2006ms}\n"
	" K.~Hamilton and P.~Richardson,\n"
	" ``A simulation of QCD radiation in top quark decays,''\n"
	" JHEP {\\bf 0702}, 069 (2007)\n"
	" [arXiv:hep-ph/0612236].\n"
	" %%CITATION = JHEPA,0702,069;%%\n");

	static ParVector<SMTopDecayer,double> interfaceQuarkWeights
	("QuarkWeights",
	"Maximum weights for the hadronic decays",
	&SMTopDecayer::_wquarkwgt, 6, 1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static ParVector<SMTopDecayer,double> interfaceLeptonWeights
	("LeptonWeights",
	"Maximum weights for the semi-leptonic decays",
	&SMTopDecayer::_wleptonwgt, 3, 1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<SMTopDecayer,double> interfaceEnhancementFactor
	("InitialEnhancementFactor",
	"The enhancement factor for initial-state radiation in the shower to ensure"
	" the weight for the matrix element correction is less than one.",
	&SMTopDecayer::_initialenhance, 1.0, 1.0, 10000.0,
	false, false, Interface::limited);

	static Parameter<SMTopDecayer,double> interfaceFinalEnhancementFactor
	("FinalEnhancementFactor",
	"The enhancement factor for final-state radiation in the shower to ensure"
	" the weight for the matrix element correction is less than one",
	&SMTopDecayer::_finalenhance, 1.6, 1.0, 1000.0,
	false, false, Interface::limited);

	static Parameter<SMTopDecayer,double> interfaceSamplingTopHardMEC
	("SamplingTopHardMEC",
	"The importance sampling power for choosing an initial xg, "
	"to sample xg according to xg^-_xg_sampling",
	&SMTopDecayer::_xg_sampling, 1.5, 1.2, 2.0,
	false, false, Interface::limited);

	static Switch<SMTopDecayer,bool> interfaceUseMEForT2
	("UseMEForT2",
	"Use the matrix element correction, if available to fill the T2"
	" region for the decay shower and don't fill using the shower",
	&SMTopDecayer::_useMEforT2, true, false, false);
	static SwitchOption interfaceUseMEForT2Shower
	(interfaceUseMEForT2,
	"Shower",
	"Use the shower to fill the T2 region",
	false);
	static SwitchOption interfaceUseMEForT2ME
	(interfaceUseMEForT2,
	"ME",
	"Use the Matrix element to fill the T2 region",
	true);

	static Reference<SMTopDecayer,ShowerAlpha> interfaceCoupling
	("Coupling",
	"Pointer to the object to calculate the coupling for the correction",
	&SMTopDecayer::_alpha, false, false, true, false, false);
	}

	double SMTopDecayer::me2(const int, const Particle & inpart,
	const ParticleVector & decay,
	MEOption meopt) const {
	// spinors etc for the decaying particle
	if(meopt==Initialize) {
	// spinors and rho
	if(inpart.id()>0)
	SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho,
	const_ptr_cast<tPPtr>(&inpart),
	incoming);
	else
	SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho,
	const_ptr_cast<tPPtr>(&inpart),
	incoming);
	ME(DecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
	PDT::Spin1Half,PDT::Spin1Half));
	}
	// setup spin info when needed
	if(meopt==Terminate) {
	// for the decaying particle
	if(inpart.id()>0) {
	SpinorWaveFunction::
	constructSpinInfo(_inHalf,const_ptr_cast<tPPtr>(&inpart),incoming,true);
	SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true);
	SpinorWaveFunction ::constructSpinInfo(_outHalf ,decay[1],outgoing,true);
	SpinorBarWaveFunction::constructSpinInfo(_outHalfBar,decay[2],outgoing,true);
	}
	else {
	SpinorBarWaveFunction::
	constructSpinInfo(_inHalfBar,const_ptr_cast<tPPtr>(&inpart),incoming,true);
	SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true);
	SpinorBarWaveFunction::constructSpinInfo(_outHalfBar,decay[1],outgoing,true);
	SpinorWaveFunction ::constructSpinInfo(_outHalf ,decay[2],outgoing,true);
	}
	}

	if ( ( decay[1]->momentum() + decay[2]->momentum() ).m()
	< decay[1]->data().constituentMass() + decay[2]->data().constituentMass() )
	return 0.0;

	// spinors for the decay product
	if(inpart.id()>0) {
	SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar ,decay[0],outgoing);
	SpinorWaveFunction ::calculateWaveFunctions(_outHalf ,decay[1],outgoing);
	SpinorBarWaveFunction::calculateWaveFunctions(_outHalfBar,decay[2],outgoing);
	}
	else {
	SpinorWaveFunction ::calculateWaveFunctions(_inHalf ,decay[0],outgoing);
	SpinorBarWaveFunction::calculateWaveFunctions(_outHalfBar,decay[1],outgoing);
	SpinorWaveFunction ::calculateWaveFunctions(_outHalf ,decay[2],outgoing);
	}
	Energy2 scale(sqr(inpart.mass()));
	if(inpart.id() == ParticleID::t) {
	//Define intermediate vector wave-function for Wplus
	tcPDPtr Wplus(getParticleData(ParticleID::Wplus));
	VectorWaveFunction inter;
	unsigned int thel,bhel,fhel,afhel;
	for(thel = 0;thel<2;++thel){
	for(bhel = 0;bhel<2;++bhel){
	inter = _wvertex->evaluate(scale,1,Wplus,_inHalf[thel],
	_inHalfBar[bhel]);
	for(afhel=0;afhel<2;++afhel){
	for(fhel=0;fhel<2;++fhel){
	ME()(thel,bhel,afhel,fhel) =
	_wvertex->evaluate(scale,_outHalf[afhel],
	_outHalfBar[fhel],inter);
	}
	}
	}
	}
	}
	else if(inpart.id() == ParticleID::tbar) {
	VectorWaveFunction inter;
	tcPDPtr Wminus(getParticleData(ParticleID::Wminus));
	unsigned int tbhel,bbhel,afhel,fhel;
	for(tbhel = 0;tbhel<2;++tbhel){
	for(bbhel = 0;bbhel<2;++bbhel){
	inter = _wvertex->
	evaluate(scale,1,Wminus,_inHalf[bbhel],_inHalfBar[tbhel]);
	for(afhel=0;afhel<2;++afhel){
	for(fhel=0;fhel<2;++fhel){
	ME()(tbhel,bbhel,fhel,afhel) =
	_wvertex->evaluate(scale,_outHalf[afhel],
	_outHalfBar[fhel],inter);
	}
	}
	}
	}
	}
	double output = (ME().contract(_rho)).real();
	if(abs(decay[1]->id())<=6) output *=3.;
	return output;
	}

	void SMTopDecayer::doinit() {
	DecayIntegrator::doinit();
	//get vertices from SM object
	tcHwSMPtr hwsm = dynamic_ptr_cast<tcHwSMPtr>(standardModel());
	if(!hwsm) throw InitException() << "Must have Herwig::StandardModel in "
	<< "SMTopDecayer::doinit()";
	_wvertex = hwsm->vertexFFW();
	//initialise
	_wvertex->init();
	//set up decay modes
	_wplus = getParticleData(ParticleID::Wplus);
	DecayPhaseSpaceModePtr mode;
	DecayPhaseSpaceChannelPtr Wchannel;
	tPDVector extpart(4);
	vector<double> wgt(1,1.0);
	extpart[0] = getParticleData(ParticleID::t);
	extpart[1] = getParticleData(ParticleID::b);
	//lepton modes
	for(int i=11; i<17;i+=2) {
	extpart[2] = getParticleData(-i);
	extpart[3] = getParticleData(i+1);
	mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
	Wchannel = new_ptr(DecayPhaseSpaceChannel(mode));
	Wchannel->addIntermediate(extpart[0],0,0.0,-1,1);
	Wchannel->addIntermediate(_wplus,0,0.0,2,3);
	Wchannel->init();
	mode->addChannel(Wchannel);
	addMode(mode,_wleptonwgt[(i-11)/2],wgt);
	}
	//quark modes
	unsigned int iz=0;
	for(int ix=1;ix<6;ix+=2) {
	for(int iy=2;iy<6;iy+=2) {
	// check that the combination of particles is allowed
	if(_wvertex->allowed(-ix,iy,ParticleID::Wminus)) {
	extpart[2] = getParticleData(-ix);
	extpart[3] = getParticleData( iy);
	mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
	Wchannel = new_ptr(DecayPhaseSpaceChannel(mode));
	Wchannel->addIntermediate(extpart[0],0,0.0,-1,1);
	Wchannel->addIntermediate(_wplus,0,0.0,2,3);
	Wchannel->init();
	mode->addChannel(Wchannel);
	addMode(mode,_wquarkwgt[iz],wgt);
	++iz;
	}
	else {
	throw InitException() << "SMTopDecayer::doinit() the W vertex"
	<< "cannot handle all the quark modes"
	<< Exception::abortnow;
	}
	}
	}
	}

	void SMTopDecayer::dataBaseOutput(ofstream & os,bool header) const {
	if(header) os << "update decayers set parameters=\"";
	// parameters for the DecayIntegrator base class
	for(unsigned int ix=0;ix<_wquarkwgt.size();++ix) {
	os << "newdef " << name() << ":QuarkWeights " << ix << " "
	<< _wquarkwgt[ix] << "\n";
	}
	for(unsigned int ix=0;ix<_wleptonwgt.size();++ix) {
	os << "newdef " << name() << ":LeptonWeights " << ix << " "
	<< _wleptonwgt[ix] << "\n";
	}
	DecayIntegrator::dataBaseOutput(os,false);
	if(header) os << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
	}

	void SMTopDecayer::doinitrun() {
	DecayIntegrator::doinitrun();
	if(initialize()) {
	for(unsigned int ix=0;ix<numberModes();++ix) {
	if(ix<3) _wleptonwgt[ix ] = mode(ix)->maxWeight();
	else _wquarkwgt [ix-3] = mode(ix)->maxWeight();
	}
	}
	}

	WidthCalculatorBasePtr SMTopDecayer::threeBodyMEIntegrator(const DecayMode & dm) const {
	// identify W decay products
	int sign = dm.parent()->id() > 0 ? 1 : -1;
	int iferm(0),ianti(0);
	for(ParticleMSet::const_iterator pit=dm.products().begin();
	pit!=dm.products().end();++pit) {
	int id = (**pit).id();
	if(id*sign != ParticleID::b) {
	if (idsign > 0 ) iferm = idsign;
	else ianti = id*sign;
	}
	}
	assert(iferm!=0&&ianti!=0);
	// work out which mode we are doing
	int imode(-1);
	for(unsigned int ix=0;ix<numberModes();++ix) {
	if(mode(ix)->externalParticles(2)->id() == ianti &&
	mode(ix)->externalParticles(3)->id() == iferm ) {
	imode = ix;
	break;
	}
	}
	assert(imode>=0);
	// get the masses we need
	Energy m[3] = {mode(imode)->externalParticles(1)->mass(),
	mode(imode)->externalParticles(3)->mass(),
	mode(imode)->externalParticles(2)->mass()};
	return
	new_ptr(ThreeBodyAllOn1IntegralCalculator<SMTopDecayer>
	(3,_wplus->mass(),_wplus->width(),0.0,*this,imode,m[0],m[1],m[2]));
	}

	InvEnergy SMTopDecayer::threeBodydGammads(const int imode, const Energy2 mt2,
	const Energy2 mffb2, const Energy mb,
	const Energy mf, const Energy mfb) const {
	Energy mffb(sqrt(mffb2));
	Energy mw(_wplus->mass());
	Energy2 mw2(sqr(mw)),gw2(sqr(_wplus->width()));
	Energy mt(sqrt(mt2));
	Energy Eb = 0.5*(mt2-mffb2-sqr(mb))/mffb;
	Energy Ef = 0.5*(mffb2-sqr(mfb)+sqr(mf))/mffb;
	Energy Ebm = sqrt(sqr(Eb)-sqr(mb));
	Energy Efm = sqrt(sqr(Ef)-sqr(mf));
	Energy2 upp = sqr(Eb+Ef)-sqr(Ebm-Efm);
	Energy2 low = sqr(Eb+Ef)-sqr(Ebm+Efm);
	InvEnergy width=(dGammaIntegrand(mffb2,upp,mt,mb,mf,mfb,mw)-
	dGammaIntegrand(mffb2,low,mt,mb,mf,mfb,mw))
	/32./mt2/mt/8/pow(Constants::pi,3)/(sqr(mffb2-mw2)+mw2*gw2);
	// couplings
	width = 0.25sqr(4.Constants::pigenerator()->standardModel()->alphaEM(mt2)/
	generator()->standardModel()->sin2ThetaW());
	width = generator()->standardModel()->CKM(mode(imode)->externalParticles(0),
	*mode(imode)->externalParticles(1));
	if(abs(mode(imode)->externalParticles(2)->id())<=6) {
	width *=3.;
	if(abs(mode(imode)->externalParticles(2)->id())%2==0)
	width =generator()->standardModel()->CKM(mode(imode)->externalParticles(2),
	*mode(imode)->externalParticles(3));
	else
	width =generator()->standardModel()->CKM(mode(imode)->externalParticles(3),
	*mode(imode)->externalParticles(2));
	}
	// final spin average
	assert(!isnan(width*GeV));
	return 0.5*width;
	}

	Energy6 SMTopDecayer::dGammaIntegrand(Energy2 mffb2, Energy2 mbf2, Energy mt,
	Energy mb, Energy mf, Energy mfb, Energy mw) const {
	Energy2 mt2(sqr(mt)) ,mb2(sqr(mb)) ,mf2(sqr(mf )),mfb2(sqr(mfb )),mw2(sqr(mw ));
	Energy4 mt4(sqr(mt2)),mb4(sqr(mb2)),mf4(sqr(mf2)),mfb4(sqr(mfb2)),mw4(sqr(mw2));
	return -mbf2 * ( + 6 * mb2 * mf2 * mfb2 * mffb2 + 6 * mb2 * mt2 * mfb2 * mffb2
	+ 6 * mb2 * mt2 * mf2 * mffb2 + 12 * mb2 * mt2 * mf2 * mfb2
	- 3 * mb2 * mfb4 * mffb2 + 3 * mb2 * mf2 * mffb2 * mffb2
	- 3 * mb2 * mf4 * mffb2 - 6 * mb2 * mt2 * mfb4
	- 6 * mb2 * mt2 * mf4 - 3 * mb4 * mfb2 * mffb2
	- 3 * mb4 * mf2 * mffb2 - 6 * mb4 * mf2 * mfb2
	+ 3 * mt4 * mf4 + 3 * mb4 * mfb4
	+ 3 * mb4 * mf4 + 3 * mt4 * mfb4
	+ 3 * mb2 * mfb2 * mffb2 * mffb2 + 3 * mt2 * mfb2 * mffb2 * mffb2
	- 3 * mt2 * mfb4 * mffb2 + 3 * mt2 * mf2 * mffb2 * mffb2
	- 3 * mt2 * mf4 * mffb2 - 3 * mt4 * mfb2 * mffb2
	- 3 * mt4 * mf2 * mffb2 - 6 * mt4 * mf2 * mfb2
	+ 6 * mt2 * mf2 * mfb2 * mffb2 + 12 * mt2 * mf2 * mw4
	+ 12 * mb2 * mfb2 * mw4 + 12 * mb2 * mt2 * mw4
	+ 6 * mw2 * mt2 * mfb2 * mbf2 - 12 * mw2 * mt2 * mf2 * mffb2
	- 6 * mw2 * mt2 * mf2 * mbf2 - 12 * mw2 * mt2 * mf2 * mfb2
	- 12 * mw2 * mb2 * mfb2 * mffb2 - 6 * mw2 * mb2 * mfb2 * mbf2
	+ 6 * mw2 * mb2 * mf2 * mbf2 - 12 * mw2 * mb2 * mf2 * mfb2
	- 12 * mw2 * mb2 * mt2 * mfb2 - 12 * mw2 * mb2 * mt2 * mf2
	+ 12 * mf2 * mfb2 * mw4 + 4 * mbf2 * mbf2 * mw4
	- 6 * mfb2 * mbf2 * mw4 - 6 * mf2 * mbf2 * mw4
	- 6 * mt2 * mbf2 * mw4 - 6 * mb2 * mbf2 * mw4
	+ 12 * mw2 * mt2 * mf4 + 12 * mw2 * mt4 * mf2
	+ 12 * mw2 * mb2 * mfb4 + 12 * mw2 * mb4 * mfb2) /mw4 / 3.;
	}

	void SMTopDecayer::initializeMECorrection(ShowerTreePtr tree, double & initial,
	double & final) {
	// check the outgoing particles
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
	ShowerParticlePtr part[2];
	unsigned int ix(0);
	for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
	part[ix]=cit->first->progenitor();
	++ix;
	}
	// check the final-state particles and get the masses
	if(abs(part[0]->id())==ParticleID::Wplus&&abs(part[1]->id())==ParticleID::b) {
	_ma=part[0]->mass();
	_mc=part[1]->mass();
	}
	else if(abs(part[1]->id())==ParticleID::Wplus&&abs(part[0]->id())==ParticleID::b) {
	_ma=part[1]->mass();
	_mc=part[0]->mass();
	}
	else {
	return;
	}
	// set the top mass
	_mt=tree->incomingLines().begin()->first->progenitor()->mass();
	// set the gluon mass
	_mg=getParticleData(ParticleID::g)->constituentMass();
	// set the radiation enhancement factors
	initial = _initialenhance;
	final = _finalenhance;
	// reduced mass parameters
	_a=sqr(_ma/_mt);
	_g=sqr(_mg/_mt);
	_c=sqr(_mc/_mt);
	useMe();
	}

	void SMTopDecayer::applyHardMatrixElementCorrection(ShowerTreePtr tree) {
	// Get b and a and put them in particle vector ba in that order...
	ParticleVector ba;
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
	for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit)
	ba.push_back(cit->first->copy());
	PPtr temp;
	if(abs(ba[0]->id())!=5) swap(ba[0],ba[1]);
	// Get the starting scales for the showers $\tilde{\kappa}_{b}$
	// $\tilde{\kappa}_{c}$. These are needed in order to know the boundary
	// of the dead zone.
	+ _ktb = _ktc = -1.;
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cjt;
	- // for(cjt = tree->incomingLines().begin();
	- // cjt!= tree->incomingLines().end();++cjt) {
	- // if(abs(cjt->first->progenitor()->id())==6)
	- // _ktb=sqr(cjt->first->progenitor()->
	- // evolutionScale(true,cjt->first->progenitor()->id()>0 ?
	- // ShowerPartnerType::QCDColourLine :
	- // ShowerPartnerType::QCDAntiColourLine)/_mt);
	- // }
	- // for(cit = tree->outgoingLines().begin();
	- // cit!= tree->outgoingLines().end();++cit) {
	- // if(abs(cit->first->progenitor()->id())==5)
	- // _ktc=sqr(cit->first->progenitor()->
	- // evolutionScale(true,cit->first->progenitor()->id()>0 ?
	- // ShowerPartnerType::QCDColourLine :
	- // ShowerPartnerType::QCDAntiColourLine)/_mt);
	- // }
	- assert(false);
	+ for(cjt = tree->incomingLines().begin();
	+ cjt!= tree->incomingLines().end();++cjt) {
	+ if(abs(cjt->first->progenitor()->id())!=6) continue;
	+ if(cjt->first->progenitor()->id()>0)
	+ _ktb=sqr(cjt->first->progenitor()->scales().QCD_c /_mt);
	+ else
	+ _ktb=sqr(cjt->first->progenitor()->scales().QCD_ac/_mt);
	+ }
	+ for(cit = tree->outgoingLines().begin();
	+ cit!= tree->outgoingLines().end();++cit) {
	+ if(abs(cit->first->progenitor()->id())!=5) continue;
	+ if(cit->first->progenitor()->id()>0)
	+ _ktc=sqr(cit->first->progenitor()->scales().QCD_c /_mt);
	+ else
	+ _ktc=sqr(cit->first->progenitor()->scales().QCD_ac/_mt);
	+ }
	if (_ktb<=0.\|\|_ktc<=0.) {
	throw Exception()
	<< "SMTopDecayer::applyHardMatrixElementCorrection()"
	<< " did not set ktb,ktc"
	<< Exception::abortnow;
	}
	// Now decide if we get an emission into the dead region.
	// If there is an emission newfs stores momenta for a,c,g
	// according to NLO decay matrix element.
	vector<Lorentz5Momentum> newfs = applyHard(ba,_ktb,_ktc);
	// If there was no gluon emitted return.
	if(newfs.size()!=3) return;
	// Sanity checks to ensure energy greater than mass etc :)
	bool check = true;
	tcPDPtr gluondata=getParticleData(ParticleID::g);
	if (newfs[0].e()<ba[0]->data().constituentMass()) check = false;
	if (newfs[1].e()<ba[1]->mass()) check = false;
	if (newfs[2].e()<gluondata->constituentMass()) check = false;
	// Return if insane:
	if (!check) return;
	// Set masses in 5-vectors:
	newfs[0].setMass(ba[0]->mass());
	newfs[1].setMass(ba[1]->mass());
	newfs[2].setMass(ZERO);
	// The next part of this routine sets the colour structure.
	// To do this for decays we assume that the gluon comes from c!
	// First create new particle objects for c, a and gluon:
	PPtr newg = gluondata->produceParticle(newfs[2]);
	PPtr newc,newa;
	newc = ba[0]->data().produceParticle(newfs[0]);
	newa = ba[1];
	newa->set5Momentum(newfs[1]);
	// set the colour lines
	ColinePtr col;
	if(ba[0]->id()>0) {
	col=ba[0]->colourLine();
	col->addColoured(newg);
	newg->colourNeighbour(newc);
	}
	else {
	col=ba[0]->antiColourLine();
	col->addAntiColoured(newg);
	newg->antiColourNeighbour(newc);
	}
	// change the existing quark and antiquark
	PPtr orig;
	for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
	if(cit->first->progenitor()->id()==newc->id()) {
	// remove old particles from colour line
	if(newc->id()>0) {
	col->removeColoured(cit->first->copy());
	col->removeColoured(cit->first->progenitor());
	}
	else {
	col->removeAntiColoured(cit->first->copy());
	col->removeAntiColoured(cit->first->progenitor());
	}
	// insert new particles
	cit->first->copy(newc);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newc,2,true)));
	cit->first->progenitor(sp);
	tree->outgoingLines()[cit->first]=sp;
	cit->first->perturbative(false);
	orig=cit->first->original();
	}
	else {
	cit->first->copy(newa);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newa,2,true)));
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,
	tShowerParticlePtr> >::const_iterator tit;
	for(tit = tree->treelinks().begin();
	tit != tree->treelinks().end();++tit) {
	if(tit->second.first && tit->second.second==cit->first->progenitor())
	break;
	}
	cit->first->progenitor(sp);
	if(tit!=tree->treelinks().end())
	tree->updateLink(tit->first,make_pair(cit->first,sp));
	tree->outgoingLines()[cit->first]=sp;
	cit->first->perturbative(true);
	}
	}
	// Add the gluon to the shower:
	ShowerParticlePtr sg =new_ptr(ShowerParticle(*newg,2,true));
	ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(orig,newg,sg));
	gluon->perturbative(false);
	tree->outgoingLines().insert(make_pair(gluon,sg));
	if(!inTheDeadRegion(_xg,_xa,_ktb,_ktc)) {
	generator()->log()
	<< "SMTopDecayer::applyHardMatrixElementCorrection()\n"
	<< "Just found a point that escaped from the dead region!\n"
	<< " _xg: " << _xg << " _xa: " << _xa
	<< " newfs.size(): " << newfs.size() << endl;
	}
	tree->hardMatrixElementCorrection(true);
	}

	vector<Lorentz5Momentum> SMTopDecayer::
	applyHard(const ParticleVector &p,double ktb, double ktc) {
	// ********************************* //
	// First we see if we get a dead //
	// region event: _xa,_xg //
	// ********************************* //
	vector<Lorentz5Momentum> fs;
	// Return if there is no (NLO) gluon emission:

	double weight = getHard(ktb,ktc);
	if(weight>1.) {
	generator()->log() << "Weight greater than 1 for hard emission in "
	<< "SMTopDecayer::applyHard xg = " << _xg
	<< " xa = " << _xa << "\n";
	weight=1.;
	}
	// Accept/Reject
	if (weight<UseRandom::rnd()\|\|p.size()!= 2) return fs;
	// Drop events if getHard returned a negative weight
	// as in events that, somehow have escaped from the dead region
	// or, worse, the allowed region.
	if(weight<0.) return fs;

	// Calculate xc by momentum conservation:
	_xc = 2.-_xa-_xg;

	// ************************************ //
	// Now we get the boosts & rotations to //
	// go from lab to top rest frame with //
	// a in the +z direction. //
	// ************************************ //
	Lorentz5Momentum pa_lab,pb_lab,pc_lab,pg_lab;
	// Calculate momentum of b:
	pb_lab = p[0]->momentum() + p[1]->momentum();
	// Define/assign momenta of c,a and the gluon:
	if(abs(p[0]->id())==5) {
	pc_lab = p[0]->momentum();
	pa_lab = p[1]->momentum();
	} else {
	pc_lab = p[1]->momentum();
	pa_lab = p[0]->momentum();
	}
	// Calculate the boost to the b rest frame:
	SpinOneLorentzRotation rot0(pb_lab.findBoostToCM());
	// Calculate the rotation matrix to position a along the +z direction
	// in the rest frame of b and does a random rotation about z:
	SpinOneLorentzRotation rot1 = rotateToZ(rot0*pa_lab);
	// Calculate the boost from the b rest frame back to the lab:
	// and the inverse of the random rotation about the z-axis and the
	// rotation required to align a with +z:
	SpinOneLorentzRotation invrot = rot0.inverse()*rot1.inverse();

	// ************************************ //
	// Now we construct the momenta in the //
	// b rest frame using _xa,_xg. //
	// First we construct b, then c and g, //
	// finally we generate a by momentum //
	// conservation. //
	// ************************************ //
	Lorentz5Momentum pa_brf, pb_brf(_mt), pc_brf, pg_brf;
	// First we set the top quark to being on-shell and at rest.
	// Second we set the energies of c and g,
	pc_brf.setE(0.5_mt(2.-_xa-_xg));
	pg_brf.setE(0.5_mt_xg);
	// then their masses,
	pc_brf.setMass(_mc);
	pg_brf.setMass(ZERO);
	// Now set the z-component of c and g. For pg we simply start from
	// _xa and _xg, while for pc we assume it is equal to minus the sum
	// of the z-components of a (assumed to point in the +z direction) and g.
	double root=sqrt(_xa_xa-4._a);
	pg_brf.setZ(_mt(1.-_xa-_xg+0.5_xa*_xg-_c+_a)/root);
	pc_brf.setZ(-1.( pg_brf.z()+_mt0.5*root));
	// Now set the y-component of c and g's momenta
	pc_brf.setY(ZERO);
	pg_brf.setY(ZERO);
	// Now set the x-component of c and g's momenta
	pg_brf.setX(sqrt(sqr(pg_brf.t())-sqr(pg_brf.z())));
	pc_brf.setX(-pg_brf.x());
	// Momenta b,c,g are now set. Now we obtain a from momentum conservation,
	pa_brf = pb_brf-pc_brf-pg_brf;
	pa_brf.setMass(pa_brf.m());
	pa_brf.rescaleEnergy();

	// ************************************ //
	// Now we orient the momenta and boost //
	// them back to the original lab frame. //
	// ************************************ //
	// As in herwig6507 we assume that, in the rest frame
	// of b, we have aligned the W boson momentum in the
	// +Z direction by rot1rot0pa_lab, therefore
	// we obtain the new pa_lab by applying:
	// invrot*pa_brf.
	pa_lab = invrot*pa_brf;
	pb_lab = invrot*pb_brf;
	pc_lab = invrot*pc_brf;
	pg_lab = invrot*pg_brf;
	fs.push_back(pc_lab);
	fs.push_back(pa_lab);
	fs.push_back(pg_lab);
	return fs;
	}

	double SMTopDecayer::getHard(double ktb, double ktc) {
	// zero the variables
	_xg = 0.;
	_xa = 0.;
	_xc = 0.;
	// Get a phase space point in the dead region:
	double volume_factor = deadRegionxgxa(ktb,ktc);
	// if outside region return -1
	if(volume_factor<0) return volume_factor;
	// Compute the weight for this phase space point:
	double weight = volume_factorme(_xa,_xg)(1.+_a-_c-_xa);
	// Alpha_S and colour factors - this hard wired Alpha_S needs removing.
	weight *= (4./3.)/Constants::pi
	(_alpha->value(_mt_mt_xg(1.-_xa+_a-_c)
	/(2.-_xg-_xa-_c)));
	return weight;
	}

	bool SMTopDecayer::softMatrixElementVeto(ShowerProgenitorPtr initial,
	ShowerParticlePtr parent,Branching br) {
	// check if we need to apply the full correction
	long id[2]={abs(initial->progenitor()->id()),abs(parent->id())};
	// the initial-state correction
	if(id[0]==ParticleID::t&&id[1]==ParticleID::t)
	{
	Energy pt=br.kinematics->pT();
	// check if hardest so far
	// if not just need to remove effect of enhancement
	bool veto(false);
	// if not hardest so far
	if(pt<initial->highestpT())
	veto=!UseRandom::rndbool(1./_initialenhance);
	// if hardest so far do calculation
	else
	{
	// values of kappa and z
	double z(br.kinematics->z()),kappa(sqr(br.kinematics->scale()/_mt));
	// parameters for the translation
	double w(1.-(1.-z)(kappa-1.)),u(1.+_a-_c-(1.-z)kappa),v(sqr(u)-4._aw*z);
	// veto if outside phase space
	if(v<0.)
	veto=true;
	// otherwise calculate the weight
	else
	{
	v = sqrt(v);
	double xa((0.5(u+v)/w+0.5(u-v)/z)),xg((1.-z)*kappa);
	double f(me(xa,xg)),
	J(0.5(u+v)/sqr(w)-0.5(u-v)/sqr(z)+_asqr(w-z)/(vw*z));
	double wgt(fJ2./kappa/(1.+sqr(z)-2.*z/kappa)/_initialenhance);
	// This next `if' prevents the hardest emission from the
	// top shower ever entering the so-called T2 region of the
	// phase space if that region is to be populated by the hard MEC.
	if(_useMEforT2&&xg>xgbcut(_ktb)) wgt = 0.;
	if(wgt>1.) {
	generator()->log() << "Violation of maximum for initial-state "
	<< " soft veto in "
	<< "SMTopDecayer::softMatrixElementVeto"
	<< "xg = " << xg << " xa = " << xa
	<< "weight = " << wgt << "\n";
	wgt=1.;
	}
	// compute veto from weight
	veto = !UseRandom::rndbool(wgt);
	}
	// if not vetoed reset max
	if(!veto) initial->highestpT(pt);
	}
	// if vetoing reset the scale
	if(veto) parent->vetoEmission(br.type,br.kinematics->scale());
	// return the veto
	return veto;
	}
	// final-state correction
	else if(id[0]==ParticleID::b&&id[1]==ParticleID::b)
	{
	Energy pt=br.kinematics->pT();
	// check if hardest so far
	// if not just need to remove effect of enhancement
	bool veto(false);
	// if not hardest so far
	if(pt<initial->highestpT()) return !UseRandom::rndbool(1./_finalenhance);
	// if hardest so far do calculation
	// values of kappa and z
	double z(br.kinematics->z()),kappa(sqr(br.kinematics->scale()/_mt));
	// momentum fractions
	double xa(1.+_a-_c-z(1.-z)kappa),r(0.5(1.+_c/(1.+_a-xa))),root(sqr(xa)-4._a);
	if(root<0.) {
	generator()->log() << "Imaginary root for final-state veto in "
	<< "SMTopDecayer::softMatrixElementVeto"
	<< "\nz = " << z << "\nkappa = " << kappa
	<< "\nxa = " << xa
	<< "\nroot^2= " << root;
	parent->vetoEmission(br.type,br.kinematics->scale());
	return true;
	}
	root=sqrt(root);
	double xg((2.-xa)(1.-r)-(z-r)root);
	// xfact (below) is supposed to equal xg/(1-z).
	double xfact(zkappa/2./(z(1.-z)kappa+_c)(2.-xa-root)+root);
	// calculate the full result
	double f(me(xa,xg));
	// jacobian
	double J(z*root);
	double wgt(fJ2.kappa/(1.+sqr(z)-2._c/kappa/z)/sqr(xfact)/_finalenhance);
	if(wgt>1.) {
	generator()->log() << "Violation of maximum for final-state soft veto in "
	<< "SMTopDecayer::softMatrixElementVeto"
	<< "xg = " << xg << " xa = " << xa
	<< "weight = " << wgt << "\n";
	wgt=1.;
	}
	// compute veto from weight
	veto = !UseRandom::rndbool(wgt);
	// if not vetoed reset max
	if(!veto) initial->highestpT(pt);
	// if vetoing reset the scale
	if(veto) parent->vetoEmission(br.type,br.kinematics->scale());
	// return the veto
	return veto;
	}
	// otherwise don't veto
	else return !UseRandom::rndbool(1./_finalenhance);
	}

	double SMTopDecayer::me(double xw,double xg) {
	double prop(1.+_a-_c-xw),xg2(sqr(xg));
	double lambda=sqrt(1.+_a_a+_c_c-2._a-2._c-2._a_c);
	double denom=(1.-2_a_a+_a+_c_a+_c_c-2.*_c);
	double wgt=-_cxg2/prop+(1.-_a+_c)xg-(prop*(1 - xg)+xg2)
	+(0.5(1.+2._a+_c)sqr(prop-xg)xg+2._aprop*xg2)/denom;
	return wgt/(lambda*prop);
	}


	// This function is auxiliary to the xab function.
	double SMTopDecayer::xgbr(int toggle) {
	return 1.+togglesqrt(_a)-_c(1.-toggle*sqrt(_a))/(1.-_a);
	}

	// This function is auxiliary to the xab function.
	double SMTopDecayer::ktr(double xgb, int toggle) {
	return 2.*xgb/
	(xgb+toggle*sqrt((1.-1./_a)
	*(xgb-xgbr( 1))
	*(xgb-xgbr(-1))));
	}

	// Function xab determines xa (2*W energy fraction) for a given value
	// of xg (2*gluon energy fraction) and kappa tilde (q tilde squared over
	// m_top squared). Hence this function allows you to draw 1: the total
	// phase space volume in the xa vs xg plane 2: for a given value of
	// kappa tilde (i.e. starting evolution scale) the associated contour
	// in the xa vs xg plane (and hence the regions that either shower can
	// populate). This calculation is done assuming the emission came from
	// the top quark i.e. kappa tilde here is the q tilde squared of the TOP
	// quark divided by m_top squared.
	double SMTopDecayer::xab(double xgb, double kt, int toggle) {
	double xab;
	if(toggle==2) {
	// This applies for g==0.&&kt==ktr(a,c,0.,xgb,1).
	xab = -2._a(xgb-2.)/(1.+_a-_c-xgb);
	} else if(toggle==1) {
	// This applies for kt==1&&g==0.
	double lambda = sqrt(sqr(xgb-1.+_a+_c)-4._a_c);
	xab = (0.5/(kt-xgb))(kt(1.+_a-_c-xgb)-lambda)
	+ (0.5/(kt+xgb(1.-kt)))(kt*(1.+_a-_c-xgb)+lambda);
	} else {
	// This is the form of xab FOR _g=0.
	double ktmktrpktmktrm = ktkt - 4._a(kt-1.)xgb*xgb
	/ (sqr(1.-_a-_c-xgb)-4._a_c);
	if(fabs(kt-(2.xgb-2._g)/(xgb-sqrt(xgbxgb-4._g)))/kt>1.e-6) {
	double lambda = sqrt((sqr(1.-_a-_c-xgb)-4._a_c)*ktmktrpktmktrm);
	xab = (0.5/(kt-xgb))(kt(1.+_a-_c-xgb)-lambda)
	+ (0.5/(kt+xgb(1.-kt)))(kt*(1.+_a-_c-xgb)+lambda);
	}
	else {
	// This is the value of xa as a function of xb when kt->infinity.
	// Where we take any kt > (2.xgb-2._g)/(xgb-sqrt(xgbxgb-4._g))
	// as being effectively infinite. This kt value is actually the
	// maximum allowed value kt can have if the phase space is calculated
	// without the approximation of _g=0 (massless gluon). This formula
	// for xab below is then valid for _g=0 AND kt=infinity only.
	xab = ( 2._c+_a(xgb-2.)
	+ 3.*xgb
	- xgb(_c+xgb+sqrt(_a_a-2.(_c-xgb+1.)_a+sqr(_c+xgb-1.)))
	- 2.
	)/2./(xgb-1.);
	}
	}
	if(isnan(xab)) {
	double ktmktrpktmktrm = ( sqr(xgbkt-2.(xgb-_g))
	-ktkt(1.-1./_a)*(xgb-xgbr( 1)-_g/(1.+sqrt(_a)))
	*(xgb-xgbr(-1)-_g/(1.-sqrt(_a)))
	)/
	(xgbxgb-(1.-1./_a)(xgb-xgbr( 1)-_g/(1.+sqrt(_a)))
	*(xgb-xgbr(-1)-_g/(1.-sqrt(_a)))
	);
	double lambda = sqrt((xgb-1.+sqr(sqrt(_a)+sqrt(_c-_g)))
	(xgb-1.+sqr(sqrt(_a)-sqrt(_c-_g)))
	ktmktrpktmktrm);
	xab = (0.5/(kt-xgb+_g))(kt(1.+_a-_c+_g-xgb)-lambda)
	+ (0.5/(kt+xgb(1.-kt)-_g))(kt*(1.+_a-_c+_g-xgb)+lambda);
	if(isnan(xab))
	throw Exception() << "TopMECorrection::xab complex x_a value.\n"
	<< " xgb = " << xgb << "\n"
	<< " xab = " << xab << "\n"
	<< " toggle = " << toggle << "\n"
	<< " ktmktrpktmktrm = " << ktmktrpktmktrm
	<< Exception::eventerror;
	}
	return xab;
	}

	// xgbcut is the point along the xg axis where the upper bound on the
	// top quark (i.e. b) emission phase space goes back on itself in the
	// xa vs xg plane i.e. roughly mid-way along the xg axis in
	// the xa vs xg Dalitz plot.
	double SMTopDecayer::xgbcut(double kt) {
	double lambda2 = 1.+_a_a+_c_c-2._a-2._c-2._a_c;
	double num1 = ktkt(1.-_a-_c);
	double num2 = 2.ktsqrt(_a(ktkt_c+lambda2(kt-1.)));
	return (num1-num2)/(ktkt-4._a*(kt-1.));
	}

	double SMTopDecayer::xaccut(double kt) {
	return 1.+_a-_c-0.25*kt;
	}

	double SMTopDecayer::z(double xac, double kt,
	int toggle1, int toggle2) {
	double z = -1.0;
	if(toggle2==0) {
	z = (kt+toggle1sqrt(kt(kt-4.(1.+_a-_c-xac))))/(2.kt);
	} else if(toggle2==1) {
	z = ((1.+_a+_c-xac)+toggle1*(1.+_a-_c-xac))
	/(2.*(1.+_a-xac));
	} else if(toggle2==2) {
	z = 0.5;
	} else {
	throw Exception() << "Cannot determine z in SMTopDecayer::z()"
	<< Exception::eventerror;
	}
	return z;
	}

	double SMTopDecayer::xgc(double xac, double kt,
	int toggle1, int toggle2) {
	double tiny(1.e-6);
	double xaToMinBoundary(xacxac-4._a);
	if(xaToMinBoundary<0) {
	if(fabs(xaToMinBoundary/(1.-_a)/(1.-_a))<tiny)
	xaToMinBoundary *= -1.;
	else
	throw Exception() << "SMTopDecayer::xgc xa not in phase space!"
	<< Exception::eventerror;
	}
	return (2.-xac)(1.-0.5(1.+_c/(1.+_a-xac)))
	-(z(xac,kt,toggle1,toggle2)-0.5*(1.+_c/(1.+_a-xac)))
	*sqrt(xaToMinBoundary);
	}

	double SMTopDecayer::xginvc0(double xg , double kt) {
	// The function xg(kappa_tilde_c,xa) surely, enough, draws a
	// line of constant kappa_tilde_c in the xg, xa Dalitz plot.
	// Such a function can therefore draw the upper and lower
	// edges of the phase space for emission from c (the b-quark).
	// However, to sample the soft part of the dead zone effectively
	// we want to generate a value of xg first and THEN distribute
	// xa in the associated allowed part of the dead zone. Hence, the
	// function we want, to define the dead zone in xa for a given
	// xg, is the inverse of xg(kappa_tilde_c,xa). The full expression
	// for xg(kappa_tilde_c,xa) is complicated and, sure enough,
	// does not invert. Therefore we try to overestimate the size
	// of the dead zone initially, rejecting events which do not
	// fall exactly inside it afterwards, with the immediate aim
	// of getting an approximate version of xg(kappa_tilde_c,xa)
	// that can be inverted. We do this by simply setting c=0 i.e.
	// the b-quark mass to zero (and the gluon mass of course), in
	// the full expression xg(...). The result of inverting this
	// function is the output of this routine (a value of xa) hence
	// the name xginvc0. xginvc0 is calculated to be,
	// xginvc0 = (1./3.)(1.+a+pow((U+sqrt(4.VVV+U*U))/2.,1./3.)
	// -Vpow(2./(U+sqrt(4.VVV+U*U)),1./3.)
	// )
	// U = 2.aaa - 66.aa + 9.aktxg + 18.akt
	// - 66.a + 27.ktxgxg - 45.ktxg +18.*kt +2. ;
	// V = -1.-aa-14.a-3.ktxg+3.kt;
	// This function, as with many functions in this ME correction,
	// is plagued by cuts that have to handled carefully in numerical
	// implementation. We have analysed the cuts and hence we implement
	// it in the following way, with a series of 'if' statements.
	//
	// A useful -definition- to know in deriving the v<0 terms is
	// that tanh^-1(z) = 0.5*(log(1.+z)-log(1.-z)).
	double u,v,output;
	u = 2._a_a_a-66._a*_a
	+9.xgkt_a+18.kt*_a
	-66._a+27.xgxgkt
	-45.xgkt+18.*kt+2.;
	v = -_a_a-14._a-3.xgkt+3.*kt-1.;
	double u2=uu,v3=vv*v;
	if(v<0.) {
	if(u>0.&&(4.v3+u2)<0.) output = cos( atan(sqrt(-4.v3-u2)/u)/3.);
	else if(u>0.&&(4.v3+u2)>0.) output = cosh(atanh(sqrt( 4.v3+u2)/u)/3.);
	else output = cos(( atan(sqrt(-4.*v3-u2)/u)
	+Constants::pi)/3.);
	output = 2.sqrt(-v);
	} else {
	output = sinh(log((u+sqrt(4.v3+u2))/(2.sqrt(v3)))/3.);
	output = 2.sqrt(v);
	}
	if(isnan(output)\|\|isinf(output)) {
	throw Exception() << "TopMECorrection::xginvc0:\n"
	<< "possible numerical instability detected.\n"
	<< "\n v = " << v << " u = " << u << "\n4.v3+u2 = " << 4.v3+u2
	<< "\n_a = " << _a << " ma = " << sqrt(_a_mt_mt/GeV2)
	<< "\n_c = " << _c << " mc = " << sqrt(_c_mt_mt/GeV2)
	<< "\n_g = " << _g << " mg = " << sqrt(_g_mt_mt/GeV2)
	<< Exception::eventerror;
	}
	return ( 1.+_a +output)/3.;
	}

	double SMTopDecayer::approxDeadMaxxa(double xg,double ktb,double ktc) {
	double maxxa(0.);
	double x = min(xginvc0(xg,ktc),
	xab(xg,(2.xg-2._g)/(xg-sqrt(xgxg-4._g)),0));
	double y(-9999999999.);
	if(xg>2.*sqrt(_g)&&xg<=xgbcut(ktb)) {
	y = max(xab(xg,ktb,0),xab(xg,1.,1));
	} else if(xg>=xgbcut(ktb)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
	y = max(xab(xg,ktr(xg,1),2),xab(xg,1.,1));
	}
	if(xg>2.*sqrt(_g)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
	if(x>=y) { maxxa = x ; }
	else { maxxa = -9999999.; }
	} else {
	maxxa = -9999999.;
	}
	return maxxa;
	}

	double SMTopDecayer::approxDeadMinxa(double xg,double ktb,double ktc) {
	double minxa(0.);
	double x = min(xginvc0(xg,ktc),
	xab(xg,(2.xg-2._g)/(xg-sqrt(xgxg-4._g)),0));
	double y(-9999999999.);
	if(xg>2.*sqrt(_g)&&xg<=xgbcut(ktb)) {
	y = max(xab(xg,ktb,0),xab(xg,1.,1));
	} else if(xg>=xgbcut(ktb)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
	if(_useMEforT2) y = xab(xg,1.,1);
	else y = max(xab(xg,ktr(xg,1),2),xab(xg,1.,1));
	}
	if(xg>2.*sqrt(_g)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
	if(x>=y) { minxa = y ; }
	else { minxa = 9999999.; }
	} else {
	minxa = 9999999.;
	}
	return minxa;
	}

	// This function returns true if the phase space point (xg,xa) is in the
	// kinematically allowed phase space.
	bool SMTopDecayer::inTheAllowedRegion(double xg , double xa) {
	bool output(true);
	if(xg<2.*sqrt(_g)\|\|xg>1.-sqr(sqrt(_a)+sqrt(_c))) output = false;
	if(xa<xab(xg,1.,1)) output = false;
	if(xa>xab(xg,(2.xg-2._g)/(xg-sqrt(xgxg-4._g)),0)) output = false;
	return output;
	}

	// This function returns true if the phase space point (xg,xa) is in the
	// approximate (overestimated) dead region.
	bool SMTopDecayer::inTheApproxDeadRegion(double xg , double xa,
	double ktb, double ktc) {
	bool output(true);
	if(!inTheAllowedRegion(xg,xa)) output = false;
	if(xa<approxDeadMinxa(xg,ktb,ktc)) output = false;
	if(xa>approxDeadMaxxa(xg,ktb,ktc)) output = false;
	return output;
	}

	// This function returns true if the phase space point (xg,xa) is in the
	// dead region.
	bool SMTopDecayer::inTheDeadRegion(double xg , double xa,
	double ktb, double ktc) {
	bool output(true);
	if(!inTheApproxDeadRegion(xg,xa,ktb,ktc)) output = false;
	if(xa>xaccut(ktc)) {
	if(xg<xgc(max(xaccut(ktc),2.*sqrt(_a)),ktc, 1,2)&&
	xg>xgc(xa,ktc, 1,0)) { output = false; }
	if(xg>xgc(max(xaccut(ktc),2.*sqrt(_a)),ktc,-1,2)&&
	xg<xgc(xa,ktc,-1,0)) { output = false; }
	}
	return output;
	}

	// This function attempts to generate a phase space point in the dead
	// region and returns the associated phase space volume factor needed for
	// the associated event weight.
	double SMTopDecayer::deadRegionxgxa(double ktb,double ktc) {
	_xg=0.;
	_xa=0.;
	// Here we set limits on xg and generate a value inside the bounds.
	double xgmin(2.*sqrt(_g)),xgmax(1.-sqr(sqrt(_a)+sqrt(_c)));
	// Generate _xg.
	if(_xg_sampling==2.) {
	_xg=xgminxgmax/(xgmin+UseRandom::rnd()(xgmax-xgmin));
	} else {
	_xg=xgmin*xgmax/pow(( pow(xgmin,_xg_sampling-1.)
	+ UseRandom::rnd()*(pow(xgmax,_xg_sampling-1.)
	-pow(xgmin,_xg_sampling-1.))
	),1./(_xg_sampling-1.));
	}
	// Here we set the bounds on _xa for given _xg.
	if(_xg<xgmin\|\|xgmin>xgmax)
	throw Exception() << "TopMECorrection::deadRegionxgxa:\n"
	<< "upper xg bound is less than the lower xg bound.\n"
	<< "\n_xg = " << _xg
	<< "\n2.sqrt(_g) = " << 2.sqrt(_g)
	<< "\n_a = " << _a << " ma = " << sqrt(_a_mt_mt/GeV2)
	<< "\n_c = " << _c << " mc = " << sqrt(_c_mt_mt/GeV2)
	<< "\n_g = " << _g << " mg = " << sqrt(_g_mt_mt/GeV2)
	<< Exception::eventerror;
	double xamin(approxDeadMinxa(_xg,ktb,ktc));
	double xamax(approxDeadMaxxa(_xg,ktb,ktc));
	// Are the bounds sensible? If not return.
	if(xamax<=xamin) return -1.;
	_xa=1.+_a-(1.+_a-xamax)*pow((1.+_a-xamin)/(1.+_a-xamax),UseRandom::rnd());
	// If outside the allowed region return -1.
	if(!inTheDeadRegion(_xg,_xa,ktb,ktc)) return -1.;
	// The integration volume for the weight
	double xg_vol,xa_vol;
	if(_xg_sampling==2.) {
	xg_vol = (xgmax-xgmin)
	/ (xgmax*xgmin);
	} else {
	xg_vol = (pow(xgmax,_xg_sampling-1.)-pow(xgmin,_xg_sampling-1.))
	/ ((_xg_sampling-1.)pow(xgmaxxgmin,_xg_sampling-1.));
	}
	xa_vol = log((1.+_a-xamin)/(1.+_a-xamax));
	// Here we return the integral volume factor multiplied by the part of the
	// weight left over which is not included in the BRACES function, i.e.
	// the part of _xg^-2 which is not absorbed in the integration measure.
	return xg_volxa_volpow(_xg,_xg_sampling-2.);
	}

	LorentzRotation SMTopDecayer::rotateToZ(Lorentz5Momentum v) {
	// compute the rotation matrix
	LorentzRotation trans;
	// rotate so in z-y plane
	trans.rotateZ(-atan2(v.y(),v.x()));
	// rotate so along Z
	trans.rotateY(-acos(v.z()/v.vect().mag()));
	// generate random rotation
	double c,s,cs;
	do
	{
	c = 2.*UseRandom::rnd()-1.;
	s = 2.*UseRandom::rnd()-1.;
	cs = cc+ss;
	}
	while(cs>1.\|\|cs==0.);
	double cost=(cc-ss)/cs,sint=2.cs/cs;
	// apply random azimuthal rotation
	trans.rotateZ(atan2(sint,cost));
	return trans;
	}
	diff --git a/Shower/Base/Evolver.cc b/Shower/Base/Evolver.cc
	--- a/Shower/Base/Evolver.cc
	+++ b/Shower/Base/Evolver.cc
	@@ -1,2051 +1,2054 @@
	// -- C++ --
	//
	// Evolver.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the Evolver class.
	//
	#include "Evolver.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "ThePEG/Utilities/EnumIO.h"
	#include "ShowerKinematics.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Handlers/EventHandler.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ShowerTree.h"
	#include "ShowerProgenitor.h"
	#include "KinematicsReconstructor.h"
	#include "PartnerFinder.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "Herwig++/Shower/ShowerHandler.h"

	using namespace Herwig;

	namespace {

	void findChildren(tShowerParticlePtr parent,set<ShowerParticlePtr> & fs) {
	for(unsigned int ix=0;ix<parent->children().size();++ix) {
	tShowerParticlePtr child=
	dynamic_ptr_cast<tShowerParticlePtr>(parent->children()[ix]);
	if(child) findChildren(child,fs);
	}
	if(parent->children().empty()) {
	if(parent->isFinalState()) fs.insert(parent);
	}
	}
	}

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

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

	void Evolver::persistentOutput(PersistentOStream & os) const {
	os << _model << _splittingGenerator << _maxtry
	<< _meCorrMode << _hardVetoMode << _hardVetoRead << _hardVetoReadOption
	<< _limitEmissions
	<< ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
	<< _vetoes << _trunc_Mode << _hardEmissionMode
	<< _colourEvolutionMethod << _reconOpt
	<< interaction_<< interactions_.size();
	for(unsigned int ix=0;ix<interactions_.size();++ix)
	os << oenum(interactions_[ix]);
	}

	void Evolver::persistentInput(PersistentIStream & is, int) {
	unsigned int isize;
	is >> _model >> _splittingGenerator >> _maxtry
	>> _meCorrMode >> _hardVetoMode >> _hardVetoRead >> _hardVetoReadOption
	>> _limitEmissions
	>> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
	>> _vetoes >> _trunc_Mode >> _hardEmissionMode
	>> _colourEvolutionMethod >> _reconOpt
	>> interaction_ >> isize;
	interactions_.resize(isize);
	for(unsigned int ix=0;ix<interactions_.size();++ix)
	is >> ienum(interactions_[ix]);
	}

	void Evolver::doinit() {
	Interfaced::doinit();
	if(interaction_==0) {
	interactions_.push_back(ShowerInteraction::QCD);
	interactions_.push_back(ShowerInteraction::QED);
	}
	else if(interaction_==1) {
	interactions_.push_back(ShowerInteraction::QCD);
	}
	else if(interaction_==2) {
	interactions_.push_back(ShowerInteraction::QED);
	interactions_.push_back(ShowerInteraction::QCD);
	}
	else if(interaction_==3) {
	interactions_.push_back(ShowerInteraction::QED);
	}
	else if(interaction_==4) {
	interactions_.push_back(ShowerInteraction::Both);
	}
	}

	ClassDescription<Evolver> Evolver::initEvolver;
	// Definition of the static class description member.

	void Evolver::Init() {

	static ClassDocumentation<Evolver> documentation
	("This class is responsible for carrying out the showering,",
	"including the kinematics reconstruction, in a given scale range,"
	"including the option of the POWHEG approach to simulated next-to-leading order"
	" radiation\\cite{Nason:2004rx}.",
	"%\\cite{Nason:2004rx}\n"
	"\\bibitem{Nason:2004rx}\n"
	" P.~Nason,\n"
	" ``A new method for combining NLO QCD with shower Monte Carlo algorithms,''\n"
	" JHEP {\\bf 0411} (2004) 040\n"
	" [arXiv:hep-ph/0409146].\n"
	" %%CITATION = JHEPA,0411,040;%%\n");

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

	static Reference<Evolver,ShowerModel> interfaceShowerModel
	("ShowerModel",
	"The pointer to the object which defines the shower evolution model.",
	&Evolver::_model, false, false, true, false, false);

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

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

	static Switch<Evolver, unsigned int> ifaceHardVetoMode
	("HardVetoMode",
	"Choice of the Hard Veto Mode",
	&Evolver::_hardVetoMode, 1, false, false);
	static SwitchOption HVoff
	(ifaceHardVetoMode,"No","hard vetos off", 0);
	static SwitchOption HVon
	(ifaceHardVetoMode,"Yes","hard vetos on", 1);
	static SwitchOption HVIS
	(ifaceHardVetoMode,"Initial", "only IS emissions vetoed", 2);
	static SwitchOption HVFS
	(ifaceHardVetoMode,"Final","only FS emissions vetoed", 3);

	static Switch<Evolver, unsigned int> ifaceHardVetoRead
	("HardVetoScaleSource",
	"If hard veto scale is to be read",
	&Evolver::_hardVetoRead, 0, false, false);
	static SwitchOption HVRcalc
	(ifaceHardVetoRead,"Calculate","Calculate from hard process", 0);
	static SwitchOption HVRread
	(ifaceHardVetoRead,"Read","Read from XComb->lastScale", 1);

	static Switch<Evolver, bool> ifaceHardVetoReadOption
	("HardVetoReadOption",
	"Apply read-in scale veto to all collisions or just the primary one?",
	&Evolver::_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<Evolver, Energy> ifaceiptrms
	("IntrinsicPtGaussian",
	"RMS of intrinsic pT of Gaussian distribution:\n"
	"2(1-Beta)exp(-sqr(intrinsicpT/RMS))/sqr(RMS)",
	&Evolver::_iptrms, GeV, ZERO, ZERO, 1000000.0*GeV,
	false, false, Interface::limited);

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

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

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

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

	static Switch<Evolver,unsigned int> interfaceLimitEmissions
	("LimitEmissions",
	"Limit the number and type of emissions for testing",
	&Evolver::_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<Evolver,bool> interfaceTruncMode
	("TruncatedShower", "Include the truncated shower?",
	&Evolver::_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<Evolver,unsigned int> interfaceHardEmissionMode
	("HardEmissionMode",
	"Whether to use ME corrections or POWHEG for the hardest emission",
	&Evolver::_hardEmissionMode, 0, false, false);
	static SwitchOption interfaceHardEmissionModeMECorrection
	(interfaceHardEmissionMode,
	"MECorrection",
	"Old fashioned ME correction",
	0);
	static SwitchOption interfaceHardEmissionModePOWHEG
	(interfaceHardEmissionMode,
	"POWHEG",
	"Powheg style hard emission",
	1);

	static Switch<Evolver,int> interfaceColourEvolutionMethod
	("ColourEvolutionMethod",
	"Choice of method for choosing the colour factor in gluon evolution",
	&Evolver::_colourEvolutionMethod, 0, false, false);
	static SwitchOption interfaceColourEvolutionMethodDefault
	(interfaceColourEvolutionMethod,
	"Default",
	"Colour factor is CA for all scales",
	0);
	static SwitchOption interfaceColourEvolutionMethodHalfCA
	(interfaceColourEvolutionMethod,
	"HalfCA",
	"Only use half the normal radiation until second scale is reached",
	1);

	static Switch<Evolver,unsigned int > interfaceInteractions
	("Interactions",
	"The interactions to be used in the shower",
	&Evolver::interaction_, 1, false, false);
	static SwitchOption interfaceInteractionsQCDFirst
	(interfaceInteractions,
	"QCDFirst",
	"QCD first then QED",
	0);
	static SwitchOption interfaceInteractionsQCDOnly
	(interfaceInteractions,
	"QCDOnly",
	"Only QCD",
	1);
	static SwitchOption interfaceInteractionsQEDFirst
	(interfaceInteractions,
	"QEDFirst",
	"QED first then QCD",
	2);
	static SwitchOption interfaceInteractionsQEDOnly
	(interfaceInteractions,
	"QEDOnly",
	"Only QED",
	3);
	static SwitchOption interfaceInteractionsBothAtOnce
	(interfaceInteractions,
	"BothAtOnce",
	"Generate both at the same time",
	4);

	static Switch<Evolver,unsigned int> interfaceReconstructionOption
	("ReconstructionOption",
	"Treatment of the reconstruction of the transverse momentum of "
	"a branching from the evolution scale.",
	&Evolver::_reconOpt, 0, false, false);
	static SwitchOption interfaceReconstructionOptionCutOff
	(interfaceReconstructionOption,
	"CutOff",
	"Use the cut-off masses in the calculation",
	0);
	static SwitchOption interfaceReconstructionOptionOffShell
	(interfaceReconstructionOption,
	"OffShell",
	"Use the off-shell masses in the calculation",
	1);

	}

	void Evolver::generateIntrinsicpT(vector<ShowerProgenitorPtr> particlesToShower) {
	_intrinsic.clear();
	if ( !ipTon() \|\| !isISRadiationON() ) return;
	// don't do anything for the moment for secondary scatters
	if( !ShowerHandler::currentHandler()->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 Evolver::setupMaximumScales(vector<ShowerProgenitorPtr> & p) {
	// return if no vetos
	if (!hardVetoOn()) 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 maximum 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 = -1.0*GeV;
	// general case calculate the scale
	if (!hardVetoXComb()\|\|
	(hardVetoReadOption()&&
	!ShowerHandler::currentHandler()->firstInteraction())) {
	// scattering process
	if(currentTree()->isHard()) {
	// 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 = max(ptmax,cjt->first->progenitor()->momentum().mt());
	}
	}
	if (ptmax < ZERO) ptmax = pcm.m();
	if(hardVetoXComb()&&hardVetoReadOption()&&
	!ShowerHandler::currentHandler()->firstInteraction()) {
	ptmax=min(ptmax,sqrt(ShowerHandler::currentHandler()
	->lastXCombPtr()->lastScale()));
	}
	}
	// 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 {
	ptmax = sqrt( ShowerHandler::currentHandler()
	->lastXCombPtr()->lastScale() );
	}
	// 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 Evolver::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());
	// number of attempts if more than one interaction switched on
	unsigned int interactionTry=0;
	do {
	try {
	// generate the showering
	doShowering(true);
	// if no vetos return
	return;
	}
	catch (InteractionVeto) {
	currentTree()->clear();
	++interactionTry;
	}
	}
	while(interactionTry<=5);
	throw Exception() << "Too many tries for shower in "
	<< "Evolver::showerHardProcess()"
	<< Exception::eventerror;
	}

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

	bool Evolver::timeLikeShower(tShowerParticlePtr particle,
	ShowerInteraction::Type type,
	bool first) {
	// don't do anything if not needed
	if(_limitEmissions == 1 \|\| hardOnly() \|\|
	( _limitEmissions == 2 && _nfs != 0) \|\|
	( _limitEmissions == 4 && _nfs + _nis != 0) ) return false;
	ShowerParticleVector theChildren;
	int ntry=0;
	do {
	++ntry;
	// generate the emission
	Branching fb;
	while (true) {
	fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
	// no emission return
	if(!fb.kinematics) return false;
	// if emission OK break
	if(!timeLikeVetoed(fb,particle,type)) break;
	// otherwise reset scale and continue - SO IS involved in veto algorithm
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	}
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	// Assign the splitting function to the emitting particle.
	// For the time being we are considering only 1->2 branching
	// 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
	tcPDPtr pdata[2];
	for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
	if(particle->id()!=fb.ids[0]) {
	for(unsigned int ix=0;ix<2;++ix) {
	tPDPtr cc(pdata[ix]->CC());
	if(cc) pdata[ix]=cc;
	}
	}
	theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
	theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
	// update the children
	particle->showerKinematics()->updateChildren(particle, theChildren,
	fb.type,true);
	// update number of emissions
	++_nfs;
	if(_limitEmissions!=0) return true;
	// shower the first particle
	timeLikeShower(theChildren[0],type,false);
	// shower the second particle
	timeLikeShower(theChildren[1],type,false);
	// that's if for old approach
	if(_reconOpt==0) break;
	// branching has happened
	particle->showerKinematics()->updateParent(particle, theChildren,
	fb.type,true);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<theChildren.size();++ix)
	particle->abandonChild(theChildren[ix]);
	theChildren.clear();
	}
	}
	while(particle->virtualMass()==ZERO&&ntry<50);
	if(first)
	particle->showerKinematics()->resetChildren(particle,theChildren);
	return true;
	}

	bool
	Evolver::spaceLikeShower(tShowerParticlePtr particle, PPtr beam,
	ShowerInteraction::Type 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(ShowerHandler::currentHandler()->firstPDF().particle() == _beam)
	pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
	if(ShowerHandler::currentHandler()->secondPDF().particle() == _beam)
	pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
	Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
	// don't do anything if not needed
	if(_limitEmissions == 2 \|\| hardOnly() \|\|
	( _limitEmissions == 1 && _nis != 0 ) \|\|
	( _limitEmissions == 4 && _nis + _nfs != 0 ) ) 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) return false;
	// if not vetoed break
	if(!spaceLikeVetoed(bb,particle,type)) break;
	// otherwise reset scale and continue
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	}
	// assign the splitting function and shower kinematics
	particle->showerKinematics(bb.kinematics);
	// For the time being we are considering only 1->2 branching
	// particles as in Sudakov form factor
	tcPDPtr part[2]={getParticleData(bb.ids[0]),
	getParticleData(bb.ids[2])};
	if(particle->id()!=bb.ids[1]) {
	if(part[0]->CC()) part[0]=part[0]->CC();
	if(part[1]->CC()) part[1]=part[1]->CC();
	}
	// 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,true);
	// 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;
	// 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,true);
	if(_limitEmissions!=0) return true;
	// perform the shower of the final-state particle
	timeLikeShower(otherChild,type,true);
	// return the emitted
	return true;
	}

	void Evolver::showerDecay(ShowerTreePtr decay) {
	_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()
	+ "->";
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) {
	if(it!=decay->outgoingLines().begin()) tag += ",";
	tag += it->first->original()->dataPtr()->name();
	}
	tag += ";";
	dm = generator()->findDecayMode(tag);
	}
	if(dm) _decayme = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
	// set the ShowerTree to be showered
	currentTree(decay);
	decay->applyTransforms();
	hardTree(HardTreePtr());
	unsigned int interactionTry=0;
	do {
	try {
	// generate the showering
	doShowering(false);
	// if no vetos return
	return;
	}
	catch (InteractionVeto) {
	currentTree()->clear();
	++interactionTry;
	}
	}
	while(interactionTry<=5);
	throw Exception() << "Too many tries for QED shower in Evolver::showerDecay()"
	<< Exception::eventerror;
	}

	bool Evolver::spaceLikeDecayShower(tShowerParticlePtr particle,
	- const map<ShowerPartnerType::Type,pair<Energy,Energy> > & maxScales,
	+ const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass,ShowerInteraction::Type type) {
	Branching fb;
	while (true) {
	fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass,
	_initialenhance,type);
	// return if no radiation
	if(!fb.kinematics) return false;
	// if not vetoed break
	if(!spaceLikeDecayVetoed(fb,particle,type)) break;
	// otherwise reset scale and continue
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	}
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	// For the time being we are considering only 1->2 branching
	// 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
	tcPDPtr pdata[2];
	for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
	if(particle->id()!=fb.ids[0]) {
	for(unsigned int ix=0;ix<2;++ix) {
	tPDPtr cc(pdata[ix]->CC());
	if(cc) pdata[ix]=cc;
	}
	}
	ShowerParticleVector theChildren;
	theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
	theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
	// some code moved to updateChildren
	particle->showerKinematics()->updateChildren(particle, theChildren,
	fb.type,true);
	// In the case of splittings which involves coloured particles,
	// set properly the colour flow of the branching.
	// update the history if needed
	_currenttree->updateInitialStateShowerProduct(_progenitor,theChildren[0]);
	_currenttree->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
	// shower the first particle
	spaceLikeDecayShower(theChildren[0],maxScales,minmass,type);
	// shower the second particle
	timeLikeShower(theChildren[1],type,true);
	// branching has happened
	return true;
	}

	vector<ShowerProgenitorPtr> Evolver::setupShower(bool hard) {
	// generate POWHEG hard emission if needed
	if(_hardEmissionMode==1) hardestEmission(hard);
	// set the initial colour partners
	setEvolutionPartners(hard,_hardtree ?
	_hardtree->interaction() : interactions_[0]);
	// generate hard me if needed
	if(_hardEmissionMode==0) hardMatrixElementCorrection(hard);
	// get the particles to be showered
	vector<ShowerProgenitorPtr> particlesToShower =
	currentTree()->extractProgenitors();
	// remake the colour partners if needed
	if(_hardEmissionMode==0 && _currenttree->hardMatrixElementCorrection()) {
	setEvolutionPartners(hard,interactions_[0]);
	_currenttree->resetShowerProducts();
	}
	// return the answer
	return particlesToShower;
	}

	void Evolver::setEvolutionPartners(bool hard,ShowerInteraction::Type type) {
	// match the particles in the ShowerTree and hardTree
	if(hardTree() && !hardTree()->connect(currentTree()))
	throw Exception() << "Can't match trees in "
	<< "Evolver::setEvolutionPartners()"
	<< Exception::eventerror;
	// extract the progenitors
	vector<ShowerParticlePtr> particles =
	currentTree()->extractProgenitorParticles();
	// sort out the colour partners
	if(hardTree()) {
	// find the partner
	for(unsigned int ix=0;ix<particles.size();++ix) {
	tHardBranchingPtr partner =
	hardTree()->particles()[particles[ix]]->colourPartner();
	if(!partner) continue;
	for(map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	it=hardTree()->particles().begin();
	it!=hardTree()->particles().end();++it) {
	if(it->second==partner) particles[ix]->partner(it->first);
	}
	if(!particles[ix]->partner())
	throw Exception() << "Can't match partners in "
	<< "Evolver::setEvolutionPartners()"
	<< Exception::eventerror;
	}
	}
	// Set the initial evolution scales
	showerModel()->partnerFinder()->
	setInitialEvolutionScales(particles,!hard,type,!_hardtree);
	}

	void Evolver::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 Evolver::startTimeLikeShower(ShowerInteraction::Type type) {
	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);
	if(output) updateHistory(progenitor()->progenitor());
	return output;
	}
	}
	bool output = hardOnly() ? false :
	timeLikeShower(progenitor()->progenitor() ,type,true) ;
	if(output) updateHistory(progenitor()->progenitor());
	return output;
	}

	bool Evolver::startSpaceLikeShower(PPtr parent, ShowerInteraction::Type type) {
	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 );
	}
	}
	return hardOnly() ? false :
	spaceLikeShower(progenitor()->progenitor(),parent,type);
	}

	bool Evolver::
	-startSpaceLikeDecayShower(const map<ShowerPartnerType::Type,pair<Energy,Energy> > & maxScales,
	+startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
	Energy minimumMass,ShowerInteraction::Type type) {
	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);
	}
	}
	return hardOnly() ? false :
	spaceLikeDecayShower(progenitor()->progenitor(),maxScales,minimumMass,type);
	}

	bool Evolver::timeLikeVetoed(const Branching & fb,
	ShowerParticlePtr particle,
	ShowerInteraction::Type type) {
	// check whether emission was harder than largest pt of hard subprocess
	if ( hardVetoFS() && fb.kinematics->pT() > _progenitor->maxHardPt() )
	return true;
	// soft matrix element correction veto
	if( softMEC()) {
	if(_hardme && _hardme->hasMECorrection()) {
	if(_hardme->softMatrixElementVeto(_progenitor,particle,fb))
	return true;
	}
	else if(_decayme && _decayme->hasMECorrection()) {
	if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
	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);
	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;
	}

	bool Evolver::spaceLikeVetoed(const Branching & bb,ShowerParticlePtr particle,
	ShowerInteraction::Type type) {
	// check whether emission was harder than largest pt of hard subprocess
	if (hardVetoIS() && bb.kinematics->pT() > _progenitor->maxHardPt())
	return true;
	// apply the soft correction
	if( softMEC() && _hardme && _hardme->hasMECorrection() ) {
	if(_hardme->softMatrixElementVeto(_progenitor,particle,bb))
	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);
	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;
	}

	bool Evolver::spaceLikeDecayVetoed( const Branching & fb,
	ShowerParticlePtr particle,
	ShowerInteraction::Type type ) {
	// apply the soft correction
	if( softMEC() && _decayme && _decayme->hasMECorrection() ) {
	if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
	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);
	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 Evolver::hardestEmission(bool hard) {
	if( ( _hardme && _hardme->hasPOWHEGCorrection()) \|\|
	(_decayme && _decayme->hasPOWHEGCorrection())) {
	if(_hardme) {
	assert(hard);
	_hardtree = _hardme->generateHardest( currentTree(),interactions_ );
	}
	else {
	assert(!hard);
	_hardtree = _decayme->generateHardest( currentTree() );
	}
	if(!_hardtree) return;
	// join up the two trees
	connectTrees(currentTree(),_hardtree,hard);
	}
	else {
	_hardtree = ShowerHandler::currentHandler()->generateCKKW(currentTree());
	}
	}

	bool Evolver::truncatedTimeLikeShower(tShowerParticlePtr particle,
	HardBranchingPtr branch,
	ShowerInteraction::Type type) {
	int ntry=0;
	do {
	++ntry;
	Branching fb;
	unsigned int iout=0;
	tcPDPtr pdata[2];
	while (true) {
	// no truncated shower break
	if(!isTruncatedShowerON()\|\|hardOnly()) break;
	// generate emission
	fb=splittingGenerator()->chooseForwardBranching(*particle,1.,type);
	// no emission break
	if(!fb.kinematics) break;
	// check haven't evolved too far
	if(fb.kinematics->scale() < branch->scale()) {
	fb=Branching();
	break;
	}
	// get the particle data objects
	for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
	if(particle->id()!=fb.ids[0]) {
	for(unsigned int ix=0;ix<2;++ix) {
	tPDPtr cc(pdata[ix]->CC());
	if(cc) pdata[ix]=cc;
	}
	}
	// find the truncated line
	iout=0;
	if(pdata[0]->id()!=pdata[1]->id()) {
	if(pdata[0]->id()==particle->id()) iout=1;
	else if (pdata[1]->id()==particle->id()) iout=2;
	}
	else if(pdata[0]->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
	// and evolution
	if(type==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(type)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	// should do base class vetos as well
	if(timeLikeVetoed(fb,particle,type)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	break;
	}
	// if no branching force trunctaed emission
	if(!fb.kinematics) {
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin=
	branch->sudakov()->createFinalStateBranching(branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT());
	- // particle->evolutionScale(branch->type(),make_pair(branch->scale(),branch->scale()));
	- assert(false);
	showerKin->initialize( *particle,PPtr() );
	IdList idlist(3);
	idlist[0] = particle->id();
	idlist[1] = branch->children()[0]->branchingParticle()->id();
	idlist[2] = branch->children()[1]->branchingParticle()->id();
	fb = Branching( showerKin, idlist, branch->sudakov(),branch->type() );
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics( fb.kinematics );
	// Assign the splitting function to the emitting particle.
	// For the time being we are considering only 1->2 branching
	// 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 theChildren;
	theChildren.push_back(new_ptr(ShowerParticle(branch->children()[0]->
	branchingParticle()->dataPtr(),true)));
	theChildren.push_back(new_ptr(ShowerParticle(branch->children()[1]->
	branchingParticle()->dataPtr(),true)));
	particle->showerKinematics()->
	updateChildren(particle, theChildren,fb.type,
	type==branch->sudakov()->interactionType());
	// shower the first particle
	if( branch->children()[0]->children().empty() ) {
	if( ! hardOnly() )
	timeLikeShower(theChildren[0],type,false);
	}
	else {
	truncatedTimeLikeShower( theChildren[0],branch->children()[0],type);
	}
	// shower the second particle
	if( branch->children()[1]->children().empty() ) {
	if( ! hardOnly() )
	timeLikeShower( theChildren[1] , type,false);
	}
	else {
	truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type);
	}
	// that's if for old approach
	if(_reconOpt==0) return true;
	// branching has happened
	particle->showerKinematics()->updateParent(particle, theChildren,
	fb.type,true);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<theChildren.size();++ix)
	particle->abandonChild(theChildren[ix]);
	theChildren.clear();
	}
	else return true;
	}
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	// Assign the splitting function to the emitting particle.
	// For the time being we are considering only 1->2 branching
	// 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 theChildren;
	theChildren.push_back( new_ptr( ShowerParticle( pdata[0], true ) ) );
	theChildren.push_back( new_ptr( ShowerParticle( pdata[1], true ) ) );
	particle->showerKinematics()->
	updateChildren( particle, theChildren , fb.type, true);
	// shower the first particle
	if( iout == 1 ) truncatedTimeLikeShower( theChildren[0], branch , type );
	else timeLikeShower( theChildren[0] , type,false);
	// shower the second particle
	if( iout == 2 ) truncatedTimeLikeShower( theChildren[1], branch , type );
	else timeLikeShower( theChildren[1] , type,false);
	// that's if for old approach
	if(_reconOpt==0) return true;
	// branching has happened
	particle->showerKinematics()->updateParent(particle, theChildren,fb.type,true);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<theChildren.size();++ix)
	particle->abandonChild(theChildren[ix]);
	theChildren.clear();
	}
	else return true;
	}
	while(ntry<50);
	return false;
	}

	bool Evolver::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
	HardBranchingPtr branch,
	ShowerInteraction::Type type) {
	tcPDFPtr pdf;
	if(ShowerHandler::currentHandler()->firstPDF().particle() == beamParticle())
	pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
	if(ShowerHandler::currentHandler()->secondPDF().particle() == beamParticle())
	pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
	Energy freeze = ShowerHandler::currentHandler()->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] = getParticleData( bb.ids[0] );
	part[1] = getParticleData( bb.ids[2] );

	//is emitter anti-particle
	if( particle->id() != bb.ids[1]) {
	if( part[0]->CC() ) part[0] = part[0]->CC();
	if( part[1]->CC() ) part[1] = part[1]->CC();
	}
	double zsplit = bb.kinematics->z();
	// apply the vetos for the truncated shower
	// if doesn't carry most of momentum
	if(type==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(type) \|\| // 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,type)) {
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	continue;
	}
	break;
	}
	}
	if( !bb.kinematics ) {
	//do the hard emission
	ShoKinPtr kinematics =
	branch->sudakov()->createInitialStateBranching( branch->scale(), z, branch->phi(),
	branch->children()[0]->pT() );
	kinematics->initialize( *particle, beam );
	// assign the splitting function and shower kinematics
	particle->showerKinematics( kinematics );
	// 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(),
	type==branch->sudakov()->interactionType() );
	// 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,
	type==branch->sudakov()->interactionType() );
	if(hardOnly()) return true;
	// perform the shower of the final-state particle
	if( timelike->children().empty() ) {
	timeLikeShower( otherChild , type,true);
	}
	else {
	truncatedTimeLikeShower( otherChild, timelike , type);
	}
	// return the emitted
	return true;
	}
	// assign the splitting function and shower kinematics
	particle->showerKinematics( bb.kinematics );
	// 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, true);
	// 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, true);
	// perform the shower of the final-state particle
	timeLikeShower( otherChild , type,true);
	// return the emitted
	return true;
	}

	bool Evolver::
	truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
	- const map<ShowerPartnerType::Type,pair<Energy,Energy> > & maxScales,
	+ const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass, HardBranchingPtr branch,
	ShowerInteraction::Type type) {
	Branching fb;
	unsigned int iout=0;
	tcPDPtr pdata[2];
	while (true) {
	// no truncated shower break
	if(!isTruncatedShowerON()\|\|hardOnly()) break;
	fb=splittingGenerator()->chooseDecayBranching(*particle,maxScales,minmass,1.,type);
	// return if no radiation
	if(!fb.kinematics) break;
	// check haven't evolved too far
	if(fb.kinematics->scale() < branch->scale()) {
	fb=Branching();
	break;
	}
	// get the particle data objects
	for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
	if(particle->id()!=fb.ids[0]) {
	for(unsigned int ix=0;ix<2;++ix) {
	tPDPtr cc(pdata[ix]->CC());
	if(cc) pdata[ix]=cc;
	}
	}
	// find the truncated line
	iout=0;
	if(pdata[0]->id()!=pdata[1]->id()) {
	if(pdata[0]->id()==particle->id()) iout=1;
	else if (pdata[1]->id()==particle->id()) iout=2;
	}
	else if(pdata[0]->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();
	if(type==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(type)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	// should do base class vetos as well
	// if not vetoed break
	if(!spaceLikeDecayVetoed(fb,particle,type)) break;
	// otherwise reset scale and continue
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	}
	// if no branching insert hard emission and continue
	if(!fb.kinematics) {
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin=
	branch->sudakov()->createDecayBranching(branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT());
	// particle->evolutionScale(branch->type(),make_pair(branch->scale(),branch->scale()));
	assert(false);
	showerKin->initialize( *particle,PPtr() );
	IdList idlist(3);
	idlist[0] = particle->id();
	idlist[1] = branch->children()[0]->branchingParticle()->id();
	idlist[2] = branch->children()[1]->branchingParticle()->id();
	assert(false);
	fb = Branching( showerKin, idlist, branch->sudakov(),ShowerPartnerType::QCDColourLine );
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics( fb.kinematics );
	// Assign the splitting function to the emitting particle.
	// For the time being we are considering only 1->2 branching
	// 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 theChildren;
	theChildren.push_back(new_ptr(ShowerParticle(branch->children()[0]->
	branchingParticle()->dataPtr(),true)));
	theChildren.push_back(new_ptr(ShowerParticle(branch->children()[1]->
	branchingParticle()->dataPtr(),true)));
	particle->showerKinematics()->
	updateChildren(particle, theChildren,fb.type,
	type==branch->sudakov()->interactionType());
	if(theChildren[0]->id()==particle->id()) {
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[0]);
	currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
	// shower the space-like particle
	if( branch->children()[0]->children().empty() ) {
	if( ! hardOnly() ) spaceLikeDecayShower(theChildren[0],maxScales,minmass,type);
	}
	else {
	truncatedSpaceLikeDecayShower( theChildren[0],maxScales,minmass,
	branch->children()[0],type);
	}
	// shower the second particle
	if( branch->children()[1]->children().empty() ) {
	if( ! hardOnly() ) timeLikeShower( theChildren[1] , type,true);
	}
	else {
	truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type);
	}
	}
	else {
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[1]);
	currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
	// shower the space-like particle
	if( branch->children()[1]->children().empty() ) {
	if( ! hardOnly() ) spaceLikeDecayShower(theChildren[1],maxScales,minmass,type);
	}
	else {
	truncatedSpaceLikeDecayShower( theChildren[1],maxScales,minmass,
	branch->children()[1],type);
	}
	// shower the second particle
	if( branch->children()[0]->children().empty() ) {
	if( ! hardOnly() ) timeLikeShower( theChildren[0] , type,true);
	}
	else {
	truncatedTimeLikeShower( theChildren[0],branch->children()[0] ,type);
	}
	}
	return true;
	}
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	// For the time being we are considering only 1->2 branching
	// 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 theChildren;
	theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
	theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
	particle->showerKinematics()->updateChildren(particle, theChildren,fb.type,true);
	// In the case of splittings which involves coloured particles,
	// set properly the colour flow of the branching.
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[0]);
	currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
	// shower the first particle
	truncatedSpaceLikeDecayShower(theChildren[0],maxScales,minmass,branch,type);
	// shower the second particle
	timeLikeShower(theChildren[1],type,true);
	// branching has happened
	return true;
	}

	bool Evolver::constructDecayTree(vector<ShowerProgenitorPtr> & particlesToShower,
	ShowerInteraction::Type inter) {
	Energy ptmax(-GeV);
	// get the maximum pt is all ready a hard tree
	if(hardTree()) {
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(particlesToShower[ix]->maximumpT(inter)>ptmax&&
	particlesToShower[ix]->progenitor()->isFinalState())
	ptmax = particlesToShower[ix]->maximumpT(inter);
	}
	}
	vector<HardBranchingPtr> spaceBranchings,allBranchings;
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(particlesToShower[ix]->progenitor()->isFinalState()) {
	HardBranchingPtr newBranch;
	if(particlesToShower[ix]->hasEmitted()) {
	newBranch =
	new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
	particlesToShower[ix]->progenitor()->
	showerKinematics()->SudakovFormFactor(),
	HardBranchingPtr(),HardBranching::Outgoing));
	constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
	}
	else {
	newBranch =
	new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
	SudakovPtr(),HardBranchingPtr(),
	HardBranching::Outgoing));
	}
	allBranchings.push_back(newBranch);
	}
	else {
	HardBranchingPtr newBranch;
	if(particlesToShower[ix]->hasEmitted()) {
	newBranch =
	new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
	particlesToShower[ix]->progenitor()->
	showerKinematics()->SudakovFormFactor(),
	HardBranchingPtr(),HardBranching::Decay));
	constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
	HardBranchingPtr last=newBranch;
	do {
	for(unsigned int ix=0;ix<last->children().size();++ix) {
	if(last->children()[ix]->branchingParticle()->id()==
	particlesToShower[ix]->id()) {
	last = last->children()[ix];
	continue;
	}
	}
	}
	while(!last->children().empty());
	last->status(HardBranching::Incoming);
	spaceBranchings.push_back(newBranch);
	allBranchings .push_back(last);
	}
	else {
	newBranch =
	new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
	SudakovPtr(),HardBranchingPtr(),
	HardBranching::Incoming));
	spaceBranchings.push_back(newBranch);
	allBranchings .push_back(newBranch);
	}
	}
	}
	HardTreePtr QCDTree = new_ptr(HardTree(allBranchings,spaceBranchings,inter));
	// set the charge partners
	ShowerParticleVector particles;
	particles.push_back(spaceBranchings.back()->branchingParticle());
	for(set<HardBranchingPtr>::iterator cit=QCDTree->branchings().begin();
	cit!=QCDTree->branchings().end();++cit) {
	if((*cit)->status()==HardBranching::Outgoing)
	particles.push_back((*cit)->branchingParticle());
	}
	// get the partners
	showerModel()->partnerFinder()->setInitialEvolutionScales(particles,true,inter,true);
	// do the inverse recon
	if(!showerModel()->kinematicsReconstructor()->
	deconstructDecayJets(QCDTree,this,inter)) {
	return false;
	}
	// clear the old shower
	currentTree()->clear();
	// set the hard tree
	hardTree(QCDTree);
	// set the charge partners
	setEvolutionPartners(false,inter);
	// get the particles to be showered
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
	particlesToShower.clear();
	// incoming particles
	for(cit=currentTree()->incomingLines().begin();
	cit!=currentTree()->incomingLines().end();++cit)
	particlesToShower.push_back(((*cit).first));
	assert(particlesToShower.size()==1);
	// outgoing particles
	for(cjt=currentTree()->outgoingLines().begin();
	cjt!=currentTree()->outgoingLines().end();++cjt) {
	particlesToShower.push_back(((*cjt).first));
	if(ptmax>ZERO) particlesToShower.back()->maximumpT(ptmax,inter);
	}
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit=hardTree()->particles().end(),
	mit = hardTree()->particles().find(particlesToShower[ix]->progenitor());
	if( mit != eit) {
	if(mit->second->status()==HardBranching::Outgoing)
	particlesToShower[ix]->progenitor()->set5Momentum(mit->second->pVector());
	}
	}
	return true;
	}

	bool Evolver::constructHardTree(vector<ShowerProgenitorPtr> & particlesToShower,
	ShowerInteraction::Type inter) {
	bool noEmission = true;
	vector<HardBranchingPtr> spaceBranchings,allBranchings;
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(particlesToShower[ix]->progenitor()->isFinalState()) {
	HardBranchingPtr newBranch;
	if(particlesToShower[ix]->hasEmitted()) {
	noEmission = false;
	newBranch =
	new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
	particlesToShower[ix]->progenitor()->
	showerKinematics()->SudakovFormFactor(),
	HardBranchingPtr(),HardBranching::Outgoing));
	constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
	}
	else {
	newBranch =
	new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
	SudakovPtr(),HardBranchingPtr(),
	HardBranching::Outgoing));
	}
	allBranchings.push_back(newBranch);
	}
	else {
	HardBranchingPtr first,last;
	if(!particlesToShower[ix]->progenitor()->parents().empty()) {
	noEmission = false;
	constructSpaceLikeLine(particlesToShower[ix]->progenitor(),
	first,last,SudakovPtr(),
	particlesToShower[ix]->original()->parents()[0]);
	}
	else {
	first = new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
	SudakovPtr(),HardBranchingPtr(),
	HardBranching::Incoming));
	if(particlesToShower[ix]->original()->parents().empty())
	first->beam(particlesToShower[ix]->original());
	else
	first->beam(particlesToShower[ix]->original()->parents()[0]);
	last = first;
	}
	spaceBranchings.push_back(first);
	allBranchings.push_back(last);
	}
	}
	if(!noEmission) {
	HardTreePtr QCDTree = new_ptr(HardTree(allBranchings,spaceBranchings,
	inter));
	// set the charge partners
	ShowerParticleVector particles;
	for(set<HardBranchingPtr>::iterator cit=QCDTree->branchings().begin();
	cit!=QCDTree->branchings().end();++cit) {
	particles.push_back((*cit)->branchingParticle());
	}
	// get the partners
	showerModel()->partnerFinder()->setInitialEvolutionScales(particles,false,
	inter,true);
	// do the inverse recon
	if(!showerModel()->kinematicsReconstructor()->
	deconstructHardJets(QCDTree,this,inter))
	throw Exception() << "Can't to shower deconstruction for QED shower in"
	<< "QEDEvolver::showerHard" << Exception::eventerror;
	// set the hard tree
	hardTree(QCDTree);
	}
	// clear the old shower
	currentTree()->clear();
	// set the charge partners
	setEvolutionPartners(true,inter);
	// get the particles to be showered
	particlesToShower = currentTree()->extractProgenitors();
	// reset momenta
	if(hardTree()) {
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit=hardTree()->particles().end(),
	mit = hardTree()->particles().find(particlesToShower[ix]->progenitor());
	if( mit != eit) {
	particlesToShower[ix]->progenitor()->set5Momentum(mit->second->showerMomentum());
	}
	}
	}
	return true;
	}

	void Evolver::constructTimeLikeLine(tHardBranchingPtr branch,
	tShowerParticlePtr particle) {
	for(unsigned int ix=0;ix<particle->children().size();++ix) {
	HardBranching::Status status = branch->status();
	tShowerParticlePtr child =
	dynamic_ptr_cast<ShowerParticlePtr>(particle->children()[ix]);
	if(child->children().empty()) {
	HardBranchingPtr newBranch =
	new_ptr(HardBranching(child,SudakovPtr(),branch,status));
	branch->addChild(newBranch);
	}
	else {
	HardBranchingPtr newBranch =
	new_ptr(HardBranching(child,child->showerKinematics()->SudakovFormFactor(),
	branch,status));
	constructTimeLikeLine(newBranch,child);
	branch->addChild(newBranch);
	}
	}
	// sort out the type of interaction
	if(!branch->children().empty()) {
	if(branch->branchingParticle()->id()==ParticleID::gamma \|\|
	branch->children()[0]->branchingParticle()->id()==ParticleID::gamma \|\|
	branch->children()[1]->branchingParticle()->id()==ParticleID::gamma)
	branch->type(ShowerPartnerType::QED);
	else {
	if(branch->branchingParticle()->id()==
	branch->children()[0]->branchingParticle()->id()) {
	if(branch->branchingParticle()->dataPtr()->iColour()==PDT::Colour8) {
	tShowerParticlePtr emittor =
	branch->branchingParticle()->showerKinematics()->z()>0.5 ?
	branch->children()[0]->branchingParticle() :
	branch->children()[1]->branchingParticle();
	if(branch->branchingParticle()->colourLine()==emittor->colourLine())
	branch->type(ShowerPartnerType::QCDAntiColourLine);
	else if(branch->branchingParticle()->antiColourLine()==emittor->antiColourLine())
	branch->type(ShowerPartnerType::QCDColourLine);
	else
	assert(false);
	}
	else if(branch->branchingParticle()->colourLine()) {
	branch->type(ShowerPartnerType::QCDColourLine);
	}
	else if(branch->branchingParticle()->antiColourLine()) {
	branch->type(ShowerPartnerType::QCDAntiColourLine);
	}
	else
	assert(false);
	}
	else if(branch->branchingParticle()->id()==ParticleID::g &&
	branch->children()[0]->branchingParticle()->id()==
	-branch->children()[1]->branchingParticle()->id()) {
	if(branch->branchingParticle()->showerKinematics()->z()>0.5)
	branch->type(ShowerPartnerType::QCDAntiColourLine);
	else
	branch->type(ShowerPartnerType::QCDColourLine);

	}
	else
	assert(false);
	}
	}
	}

	void Evolver::constructSpaceLikeLine(tShowerParticlePtr particle,
	HardBranchingPtr & first,
	HardBranchingPtr & last,
	SudakovPtr sud,PPtr beam) {
	if(!particle) return;
	if(!particle->parents().empty()) {
	tShowerParticlePtr parent =
	dynamic_ptr_cast<ShowerParticlePtr>(particle->parents()[0]);
	SudakovPtr newSud=particle->showerKinematics()->SudakovFormFactor();
	constructSpaceLikeLine(parent,first,last,newSud,beam);
	}
	HardBranchingPtr newBranch =
	new_ptr(HardBranching(particle,sud,last,HardBranching::Incoming));
	newBranch->beam(beam);
	if(!first) {
	first=newBranch;
	last =newBranch;
	return;
	}
	last->addChild(newBranch);
	tShowerParticlePtr timeChild =
	dynamic_ptr_cast<ShowerParticlePtr>(particle->parents()[0]->children()[1]);
	HardBranchingPtr timeBranch;
	if(!timeChild->children().empty()) {
	timeBranch =
	new_ptr(HardBranching(timeChild,
	timeChild->showerKinematics()->SudakovFormFactor(),
	last,HardBranching::Outgoing));
	constructTimeLikeLine(timeBranch,timeChild);
	}
	else {
	timeBranch =
	new_ptr(HardBranching(timeChild,SudakovPtr(),last,HardBranching::Outgoing));
	}
	last->addChild(timeBranch);
	last=newBranch;
	}

	void Evolver::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;
	IdList br(3);
	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.second;
	if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
	sudakov=cjt->second.first;
	break;
	}
	}
	if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
	<< "Evolver::connectTrees() for ISR"
	<< Exception::runerror;
	(**cit).parent()->sudakov(sudakov);
	}
	// Sudakovs for FSR
	else if(!(**cit).children().empty()) {
	++_nfs;
	IdList br(3);
	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.second;
	if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
	sudakov=cjt->second.first;
	break;
	}
	}
	if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
	<< "Evolver::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());
	}
	showerModel()->partnerFinder()->
	setInitialEvolutionScales(particles,!hard,hardTree->interaction(),
	!hardTree->partnersSet());
	hardTree->partnersSet(true);
	// inverse reconstruction
	if(hard)
	showerModel()->kinematicsReconstructor()->
	deconstructHardJets(hardTree,ShowerHandler::currentHandler()->evolver(),
	hardTree->interaction());
	else
	showerModel()->kinematicsReconstructor()->
	deconstructDecayJets(hardTree,ShowerHandler::currentHandler()->evolver(),
	hardTree->interaction());
	// now reset the momenta of the showering particles
	vector<ShowerProgenitorPtr> particlesToShower;
	for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit=showerTree->incomingLines().begin();
	cit!=showerTree->incomingLines().end();++cit )
	particlesToShower.push_back(cit->first);
	// extract the showering particles
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cit=showerTree->outgoingLines().begin();
	cit!=showerTree->outgoingLines().end();++cit )
	particlesToShower.push_back(cit->first);
	// match them
	vector<bool> matched(particlesToShower.size(),false);
	for(set<HardBranchingPtr>::const_iterator cit=hardTree->branchings().begin();
	cit!=hardTree->branchings().end();++cit) {
	Energy2 dmin( 1e30*GeV2 );
	int iloc(-1);
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(matched[ix]) continue;
	if( (**cit).branchingParticle()->id() != particlesToShower[ix]->progenitor()->id() ) continue;
	if( (**cit).branchingParticle()->isFinalState() !=
	particlesToShower[ix]->progenitor()->isFinalState() ) continue;
	Energy2 dtest =
	sqr( particlesToShower[ix]->progenitor()->momentum().x() - (**cit).showerMomentum().x() ) +
	sqr( particlesToShower[ix]->progenitor()->momentum().y() - (**cit).showerMomentum().y() ) +
	sqr( particlesToShower[ix]->progenitor()->momentum().z() - (**cit).showerMomentum().z() ) +
	sqr( particlesToShower[ix]->progenitor()->momentum().t() - (**cit).showerMomentum().t() );
	// add mass difference for identical particles (e.g. Z0 Z0 production)
	dtest += 1e10*sqr(particlesToShower[ix]->progenitor()->momentum().m()-
	(**cit).showerMomentum().m());
	if( dtest < dmin ) {
	iloc = ix;
	dmin = dtest;
	}
	}
	if(iloc<0) throw Exception() << "Failed to match shower and hard trees in Evolver::hardestEmission"
	<< Exception::eventerror;
	particlesToShower[iloc]->progenitor()->set5Momentum((**cit).showerMomentum());
	matched[iloc] = true;
	}
	// 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());
	boost.boost (-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
	decayTree->transform(boost,true);
	}
	}

	void Evolver::doShowering(bool hard) {
	// order of the interactions
	bool showerOrder(true);
	// zero number of emissions
	_nis = _nfs = 0;
	// extract particles to shower
	vector<ShowerProgenitorPtr> particlesToShower(setupShower(hard));
	// setup the maximum scales for the shower
	setupMaximumScales(particlesToShower);
	// specifc 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()) {
	minmass += max( particlesToShower[ix]->progenitor()->mass(),
	particlesToShower[ix]->progenitor()->dataPtr()->constituentMass() );
	}
	else {
	mIn = particlesToShower[ix]->progenitor()->mass();
	}
	}
	// throw exception if decay can'ty happen
	if ( minmass > mIn ) {
	throw Exception() << "Evolver.cc: Mass of decaying particle is "
	<< "below constituent masses of decay products."
	<< Exception::eventerror;
	}
	}
	// check if interactions in right order
	if(hardTree() && hardTree()->interaction()!=interactions_[0]) {
	showerOrder = false;
	swap(interactions_[0],interactions_[1]);
	}
	// loop over possible interactions
	for(unsigned int inter=0;inter<interactions_.size();++inter) {
	// set up for second pass if required
	if(inter!=0) {
	// zero intrinsic pt so only added first time round
	intrinsicpT().clear();
	// construct the tree and throw veto if not possible
	if(!(hard ?
	constructHardTree (particlesToShower,interactions_[inter]) :
	constructDecayTree(particlesToShower,interactions_[inter])))
	throw InteractionVeto();
	}
	// 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,interactions_[inter]);
	_nis = _nfs = 0;
	}
	// generate the shower
	// pick random starting point
	unsigned int istart=UseRandom::irnd(particlesToShower.size());
	unsigned int istop = particlesToShower.size();
	// loop over particles with random starting point
	for(unsigned int ix=istart;ix<=istop;++ix) {
	if(ix==particlesToShower.size()) {
	if(istart!=0) {
	istop = istart-1;
	ix=0;
	}
	else break;
	}
	// extract the progenitor
	progenitor(particlesToShower[ix]);
	// final-state radiation
	if(progenitor()->progenitor()->isFinalState()) {
	if(!isFSRadiationON()) continue;
	// perform shower
	progenitor()->hasEmitted(startTimeLikeShower(interactions_[inter]));
	}
	// initial-state radiation
	else {
	if(!isISRadiationON()) continue;
	// hard process
	if(hard) {
	// get the PDF
	setBeamParticle(_progenitor->beam());
	assert(beamParticle());
	// 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,
	interactions_[inter]));
	}
	// 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
	- map<ShowerPartnerType::Type,pair<Energy,Energy> > maxScales;
	- // map<ShowerPartnerType::Type,pair<Energy,Energy> >
	- // maxScales = progenitor()->progenitor()->evolutionScales();
	- // // starting scale is mass
	- // Energy startScale=progenitor()->progenitor()->mass();
	- // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::iterator it=maxScales.begin();
	- // it!=maxScales.end();++it)
	- // progenitor()->progenitor()->evolutionScale(it->first,make_pair(startScale,startScale));
	- assert(false);
	+ 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,
	interactions_[inter]));
	}
	}
	}
	// do the kinematic reconstruction, checking if it worked
	reconstructed = hard ?
	showerModel()->kinematicsReconstructor()->
	reconstructHardJets (currentTree(),intrinsicpT(),interactions_[inter]) :
	showerModel()->kinematicsReconstructor()->
	reconstructDecayJets(currentTree(),interactions_[inter]);
	}
	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 Evolver::showerDecay()"
	<< Exception::eventerror;
	}
	}
	// tree has now showered
	_currenttree->hasShowered(true);
	if(!showerOrder) swap(interactions_[0],interactions_[1]);
	hardTree(HardTreePtr());
	}
	diff --git a/Shower/Base/Evolver.h b/Shower/Base/Evolver.h
	--- a/Shower/Base/Evolver.h
	+++ b/Shower/Base/Evolver.h
	@@ -1,772 +1,772 @@
	// -- C++ --
	//
	// Evolver.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_Evolver_H
	#define HERWIG_Evolver_H
	//
	// This is the declaration of the Evolver class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "Herwig++/Shower/SplittingFunctions/SplittingGenerator.h"
	#include "ShowerModel.h"
	#include "ThePEG/PDF/BeamParticleData.h"
	#include "ShowerTree.h"
	#include "ShowerProgenitor.fh"
	#include "Herwig++/Shower/ShowerHandler.fh"
	#include "Branching.h"
	#include "ShowerVeto.h"
	#include "HardTree.h"
	#include "ThePEG/Handlers/XComb.h"
	#include "Evolver.fh"
	#include "Herwig++/MatrixElement/HwMEBase.h"
	#include "Herwig++/Decay/HwDecayerBase.h"

	namespace Herwig {

	using namespace ThePEG;

	/**\ingroup Shower
	* Exception class
	* used to communicate failure of QED shower
	*/
	struct InteractionVeto {};

	/** \ingroup Shower
	* The Evolver class class performs the sohwer evolution of hard scattering
	* and decay processes in Herwig++.
	*
	* @see \ref EvolverInterfaces "The interfaces"
	* defined for Evolver.
	*/
	class Evolver: public Interfaced {

	/**
	* The ShowerHandler is a friend to set some parameters at initialisation
	*/
	friend class ShowerHandler;

	public:

	/**
	* Pointer to an XComb object
	*/
	typedef Ptr<XComb>::pointer XCPtr;

	public:

	/**
	* Default Constructor
	*/
	Evolver() : _maxtry(100), _meCorrMode(1), _hardVetoMode(1),
	_hardVetoRead(0), _reconOpt(0), _hardVetoReadOption(false),
	_iptrms(ZERO), _beta(0.), _gamma(ZERO), _iptmax(),
	_limitEmissions(0), _initialenhance(1.), _finalenhance(1.),
	interaction_(1), _trunc_Mode(true), _hardEmissionMode(0),
	_colourEvolutionMethod(0)
	{}

	/**
	* Members to perform the shower
	*/
	//@{
	/**
	* Perform the shower of the hard process
	*/
	virtual void showerHardProcess(ShowerTreePtr,XCPtr);

	/**
	* Perform the shower of a decay
	*/
	virtual void showerDecay(ShowerTreePtr);
	//@}

	/**
	* Access to the flags and shower variables
	*/
	//@{
	/**
	* Is there any showering switched on
	*/
	bool showeringON() const { return isISRadiationON() \|\| isFSRadiationON(); }

	/**
	* It returns true/false if the initial-state radiation is on/off.
	*/
	bool isISRadiationON() const { return _splittingGenerator->isISRadiationON(); }

	/**
	* It returns true/false if the final-state radiation is on/off.
	*/
	bool isFSRadiationON() const { return _splittingGenerator->isFSRadiationON(); }

	/**
	* Get the ShowerModel
	*/
	ShowerModelPtr showerModel() const {return _model;}

	/**
	* Get the SplittingGenerator
	*/
	tSplittingGeneratorPtr splittingGenerator() const { return _splittingGenerator; }
	//@}

	/**
	* Connect the Hard and Shower trees
	*/
	virtual void connectTrees(ShowerTreePtr showerTree, HardTreePtr hardTree, bool hard );

	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:

	/**
	* Perform the shower
	*/
	void doShowering(bool hard);

	/**
	* Generate the hard matrix element correction
	*/
	virtual void hardMatrixElementCorrection(bool);

	/**
	* Generate the hardest emission
	*/
	virtual void hardestEmission(bool hard);

	/**
	* Extract the particles to be showered, set the evolution scales
	* and apply the hard matrix element correction
	* @param hard Whether this is a hard process or decay
	* @return The particles to be showered
	*/
	virtual vector<ShowerProgenitorPtr> setupShower(bool hard);

	/**
	* set the colour partners
	*/
	virtual void setEvolutionPartners(bool hard,ShowerInteraction::Type);

	/**
	* Methods to perform the evolution of an individual particle, including
	* recursive calling on the products
	*/
	//@{
	/**
	* It does the forward evolution of the time-like input particle
	* (and recursively for all its radiation products).
	* accepting only emissions which conforms to the showerVariables
	* and soft matrix element correction.
	* If at least one emission has occurred then the method returns true.
	* @param particle The particle to be showered
	*/
	virtual bool timeLikeShower(tShowerParticlePtr particle, ShowerInteraction::Type,
	bool first);

	/**
	* It does the backward evolution of the space-like input particle
	* (and recursively for all its time-like radiation products).
	* accepting only emissions which conforms to the showerVariables.
	* If at least one emission has occurred then the method returns true
	* @param particle The particle to be showered
	* @param beam The beam particle
	*/
	virtual bool spaceLikeShower(tShowerParticlePtr particle,PPtr beam,
	ShowerInteraction::Type);

	/**
	* If does the forward evolution of the input on-shell particle
	* involved in a decay
	* (and recursively for all its time-like radiation products).
	* accepting only emissions which conforms to the showerVariables.
	* @param particle The particle to be showered
	* @param maxscale The maximum scale for the shower.
	* @param minimumMass The minimum mass of the final-state system
	*/
	virtual bool
	spaceLikeDecayShower(tShowerParticlePtr particle,
	- const map<ShowerPartnerType::Type,pair<Energy,Energy> > & maxScales,
	+ const ShowerParticle::EvolutionScales & maxScales,
	Energy minimumMass,ShowerInteraction::Type);

	/**
	* Truncated shower from a time-like particle
	*/
	virtual bool truncatedTimeLikeShower(tShowerParticlePtr particle,
	HardBranchingPtr branch,
	ShowerInteraction::Type type);

	/**
	* Truncated shower from a space-like particle
	*/
	virtual bool truncatedSpaceLikeShower(tShowerParticlePtr particle,PPtr beam,
	HardBranchingPtr branch,
	ShowerInteraction::Type type);

	/**
	* Truncated shower from a time-like particle
	*/
	virtual bool truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
	- const map<ShowerPartnerType::Type,pair<Energy,Energy> > & maxScales,
	+ const ShowerParticle::EvolutionScales & maxScales,
	Energy minimumMass, HardBranchingPtr branch,
	ShowerInteraction::Type type);
	//@}

	/**
	* Switches for matrix element corrections
	*/
	//@{
	/**
	* Any ME correction?
	*/
	bool MECOn() const {
	return _meCorrMode > 0 && _hardEmissionMode==0;
	}

	/**
	* Any hard ME correction?
	*/
	bool hardMEC() const {
	return (_meCorrMode == 1 \|\| _meCorrMode == 2) && _hardEmissionMode==0;
	}

	/**
	* Any soft ME correction?
	*/
	bool softMEC() const {
	return (_meCorrMode == 1 \|\| _meCorrMode > 2) && _hardEmissionMode==0;
	}
	//@}

	/**
	* Is the truncated shower on?
	*/
	bool isTruncatedShowerON() const {return _trunc_Mode;}

	/**
	* Switch for intrinsic pT
	*/
	//@{
	/**
	* Any intrinsic pT?
	*/
	bool ipTon() const {
	return _iptrms != ZERO \|\| ( _beta == 1.0 && _gamma != ZERO && _iptmax !=ZERO );
	}
	//@}

	/*@name Additional shower vetoes /
	//@{
	/**
	* Insert a veto.
	*/
	void addVeto (ShowerVetoPtr v) { _vetoes.push_back(v); }

	/**
	* Remove a veto.
	*/
	void removeVeto (ShowerVetoPtr v) {
	vector<ShowerVetoPtr>::iterator vit = find(_vetoes.begin(),_vetoes.end(),v);
	if (vit != _vetoes.end())
	_vetoes.erase(vit);
	}

	//@}

	/**
	* Switches for vetoing hard emissions
	*/
	//@{
	/**
	* Vetos on?
	*/
	bool hardVetoOn() const { return _hardVetoMode > 0; }

	/**
	* veto hard emissions in IS shower?
	*/
	bool hardVetoIS() const { return _hardVetoMode == 1 \|\| _hardVetoMode == 2; }

	/**
	* veto hard emissions in FS shower?
	*/
	bool hardVetoFS() const { return _hardVetoMode == 1 \|\| _hardVetoMode > 2; }

	/**
	* veto hard emissions according to lastScale from XComb?
	*/
	bool hardVetoXComb() const {return (_hardVetoRead == 1);}

	/**
	* Returns true if the hard veto read-in is to be applied to only
	* the primary collision and false otherwise.
	*/
	bool hardVetoReadOption() const {return _hardVetoReadOption;}
	//@}

	/**
	* Enhancement factors for radiation needed to generate the soft matrix
	* element correction.
	*/
	//@{
	/**
	* Access the enhancement factor for initial-state radiation
	*/
	double initialStateRadiationEnhancementFactor() const { return _initialenhance; }

	/**
	* Access the enhancement factor for final-state radiation
	*/
	double finalStateRadiationEnhancementFactor() const { return _finalenhance; }

	/**
	* Set the enhancement factor for initial-state radiation
	*/
	void initialStateRadiationEnhancementFactor(double in) { _initialenhance=in; }

	/**
	* Set the enhancement factor for final-state radiation
	*/
	void finalStateRadiationEnhancementFactor(double in) { _finalenhance=in; }
	//@}

	/**
	* Access to set/get the HardTree currently beinging showered
	*/
	//@{
	/**
	* The HardTree currently being showered
	*/
	tHardTreePtr hardTree() {return _hardtree;}

	/**
	* The HardTree currently being showered
	*/
	void hardTree(tHardTreePtr in) {_hardtree = in;}
	//@}

	/**
	* Access/set the beam particle for the current initial-state shower
	*/
	//@{
	/**
	* Get the beam particle data
	*/
	Ptr<BeamParticleData>::const_pointer beamParticle() const { return _beam; }

	/**
	* Set the beam particle data
	*/
	void setBeamParticle(Ptr<BeamParticleData>::const_pointer in) { _beam=in; }
	//@}

	/**
	* Set/Get the current tree being evolverd for inheriting classes
	*/
	//@{
	/**
	* Get the tree
	*/
	tShowerTreePtr currentTree() { return _currenttree; }

	/**
	* Set the tree
	*/
	void currentTree(tShowerTreePtr tree) { _currenttree=tree; }

	//@}

	/**
	* Access the maximum number of attempts to generate the shower
	*/
	unsigned int maximumTries() const { return _maxtry; }

	/**
	* Set/Get the ShowerProgenitor for the current shower
	*/
	//@{
	/**
	* Access the progenitor
	*/
	ShowerProgenitorPtr progenitor() { return _progenitor; }

	/**
	* Set the progenitor
	*/
	void progenitor(ShowerProgenitorPtr in) { _progenitor=in; }
	//@}

	/**
	* Calculate the intrinsic \f$p_T\f$.
	*/
	virtual void generateIntrinsicpT(vector<ShowerProgenitorPtr>);

	/**
	* Access to the intrinsic \f$p_T\f$ for inheriting classes
	*/
	map<tShowerProgenitorPtr,pair<Energy,double> > & intrinsicpT() { return _intrinsic; }

	/**
	* find the maximally allowed pt acc to the hard process.
	*/
	void setupMaximumScales(vector<ShowerProgenitorPtr> &);

	protected:

	/**
	* Start the shower of a timelike particle
	*/
	virtual bool startTimeLikeShower(ShowerInteraction::Type);

	/**
	* Update of the time-like stuff
	*/
	void updateHistory(tShowerParticlePtr particle);

	/**
	* Start the shower of a spacelike particle
	*/
	virtual bool startSpaceLikeShower(PPtr,ShowerInteraction::Type);

	/**
	* Start the shower of a spacelike particle
	*/
	virtual bool
	- startSpaceLikeDecayShower(const map<ShowerPartnerType::Type,pair<Energy,Energy> > & maxScales,
	+ startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
	Energy minimumMass,ShowerInteraction::Type);

	/**
	* Vetos for the timelike shower
	*/
	virtual bool timeLikeVetoed(const Branching &,ShowerParticlePtr,
	ShowerInteraction::Type);

	/**
	* Vetos for the spacelike shower
	*/
	virtual bool spaceLikeVetoed(const Branching &,ShowerParticlePtr,
	ShowerInteraction::Type);

	/**
	* Vetos for the spacelike shower
	*/
	virtual bool spaceLikeDecayVetoed(const Branching &,ShowerParticlePtr,
	ShowerInteraction::Type);

	/**
	* Only generate the hard emission, for testing only.
	*/
	bool hardOnly() const {return _limitEmissions==3;}

	/**
	* Members to construct the HardTree from the shower if needed
	*/
	//@{
	/**
	* Construct the tree for a scattering process
	*/
	bool constructHardTree(vector<ShowerProgenitorPtr> & particlesToShower,
	ShowerInteraction::Type inter);

	/**
	* Construct the tree for a decay process
	*/
	bool constructDecayTree(vector<ShowerProgenitorPtr> & particlesToShower,
	ShowerInteraction::Type inter);

	/**
	* Construct a time-like line
	*/
	void constructTimeLikeLine(tHardBranchingPtr branch,tShowerParticlePtr particle);

	/**
	* Construct a space-like line
	*/
	void constructSpaceLikeLine(tShowerParticlePtr particle,
	HardBranchingPtr & first, HardBranchingPtr & last,
	SudakovPtr sud,PPtr beam);
	//@}

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	protected:

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

	private:

	/**
	* Get the octet -> octet octet reduction factor.
	*/
	double getReductionFactor(tShowerParticlePtr particle);

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

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

	private:

	/**
	* Pointer to the model for the shower evolution model
	*/
	ShowerModelPtr _model;

	/**
	* Pointer to the splitting generator
	*/
	SplittingGeneratorPtr _splittingGenerator;

	/**
	* Maximum number of tries to generate the shower of a particular tree
	*/
	unsigned int _maxtry;

	/**
	* Matrix element correction switch
	*/
	unsigned int _meCorrMode;

	/**
	* Hard emission veto switch
	*/
	unsigned int _hardVetoMode;

	/**
	* Hard veto to be read switch
	*/
	unsigned int _hardVetoRead;

	/**
	* Control of the reconstruction option
	*/
	unsigned int _reconOpt;

	/**
	* If hard veto pT scale is being read-in this determines
	* whether the read-in value is applied to primary and
	* secondary (MPI) scatters or just the primary one, with
	* the usual computation of the veto being performed for
	* the secondary (MPI) scatters.
	*/
	bool _hardVetoReadOption;

	/**
	* rms intrinsic pT of Gaussian distribution
	*/
	Energy _iptrms;

	/**
	* Proportion of inverse quadratic intrinsic pT distribution
	*/
	double _beta;

	/**
	* Parameter for inverse quadratic: 2BetaGamma/(sqr(Gamma)+sqr(intrinsicpT))
	*/
	Energy _gamma;

	/**
	* Upper bound on intrinsic pT for inverse quadratic
	*/
	Energy _iptmax;

	/**
	* Limit the number of emissions for testing
	*/
	unsigned int _limitEmissions;

	/**
	* The progenitor of the current shower
	*/
	ShowerProgenitorPtr _progenitor;

	/**
	* Matrix element
	*/
	HwMEBasePtr _hardme;

	/**
	* Decayer
	*/
	HwDecayerBasePtr _decayme;

	/**
	* The ShowerTree currently being showered
	*/
	ShowerTreePtr _currenttree;

	/**
	* The HardTree currently being showered
	*/
	HardTreePtr _hardtree;

	/**
	* Radiation enhancement factors for use with the veto algorithm
	* if needed by the soft matrix element correction
	*/
	//@{
	/**
	* Enhancement factor for initial-state radiation
	*/
	double _initialenhance;

	/**
	* Enhancement factor for final-state radiation
	*/
	double _finalenhance;
	//@}

	/**
	* The beam particle data for the current initial-state shower
	*/
	Ptr<BeamParticleData>::const_pointer _beam;

	/**
	* Storage of the intrinsic \f$p_t\f$ of the particles
	*/
	map<tShowerProgenitorPtr,pair<Energy,double> > _intrinsic;

	/**
	* Vetoes
	*/
	vector<ShowerVetoPtr> _vetoes;

	/**
	* number of IS emissions
	*/
	unsigned int _nis;

	/**
	* Number of FS emissions
	*/
	unsigned int _nfs;

	/**
	* The option for wqhich interactions to use
	*/
	unsigned int interaction_;

	/**
	* Interactions allowed in the shower
	*/
	vector<ShowerInteraction::Type> interactions_;

	/**
	* Truncated shower switch
	*/
	bool _trunc_Mode;

	/**
	* Count of the number of truncated emissions
	*/
	unsigned int _truncEmissions;

	/**
	* Mode for the hard emissions
	*/
	unsigned int _hardEmissionMode;

	/**
	* Colour evolution method
	*/
	int _colourEvolutionMethod;
	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	/** @endcond */

	}

	#endif /* HERWIG_Evolver_H */
	diff --git a/Shower/Base/PartnerFinder.cc b/Shower/Base/PartnerFinder.cc
	--- a/Shower/Base/PartnerFinder.cc
	+++ b/Shower/Base/PartnerFinder.cc
	@@ -1,520 +1,522 @@
	// -- C++ --
	//
	// PartnerFinder.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the 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 "ShowerParticle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Utilities/Debug.h"

	using namespace Herwig;

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

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

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

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

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

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

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

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

	AbstractClassDescription<PartnerFinder> PartnerFinder::initPartnerFinder;
	// Definition of the static class description member.

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

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

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

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

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

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

	vector< pair<double, tShowerParticlePtr> >
	PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
	const ShowerParticleVector &particles) {
	vector< pair<double, tShowerParticlePtr> > partners;
	ShowerParticleVector::const_iterator cjt;
	for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
	if(!(cjt)->data().charged() \|\| particle == cjt) continue;
	double charge = double(particle->data().iCharge()(cjt)->data().iCharge());
	if( FS(particle) != FS(cjt) ) charge =-1.;
	if( QEDPartner_ != 0 ) {
	// only include II and FF as requested
	if( QEDPartner_ == 1 && FS(particle) != FS(*cjt) )
	continue;
	// ony include IF is requested
	else if(QEDPartner_ == 2 && FS(particle) == FS(*cjt) )
	continue;
	}
	// only keep positive dipoles
	if(charge<0.) partners.push_back(make_pair(-charge,*cjt));
	}
	return partners;
	}
	diff --git a/Shower/SplittingFunctions/SplittingFunction.cc b/Shower/SplittingFunctions/SplittingFunction.cc
	--- a/Shower/SplittingFunctions/SplittingFunction.cc
	+++ b/Shower/SplittingFunctions/SplittingFunction.cc
	@@ -1,922 +1,897 @@
	// -- C++ --
	//
	// SplittingFunction.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the 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/Base/ShowerParticle.h"

	using namespace Herwig;

	// Static variable needed for the type description system in ThePEG.
	DescribeAbstractClass<SplittingFunction,Interfaced>
	describeThePEGSplittingFunction("Herwig::SplittingFunction", "Herwig.so");

	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 Switch<SplittingFunction,ShowerInteraction::Type>
	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 Switch<SplittingFunction,int> interfaceSplittingColourMethod
	("SplittingColourMethod",
	"Choice of assigning colour in 8->88 splittings.",
	&SplittingFunction::_splittingColourMethod, 0, false, false);
	static SwitchOption interfaceSplittingColourMethodRandom
	(interfaceSplittingColourMethod,
	"Random",
	"Choose colour assignments randomly.",
	0);
	static SwitchOption interfaceSplittingColourMethodCorrectLines
	(interfaceSplittingColourMethod,
	"CorrectLines",
	"Choose correct lines for colour.",
	1);
	static SwitchOption interfaceSplittingColourMethodRandomRecord
	(interfaceSplittingColourMethod,
	"RandomRecord",
	"Choose colour assignments randomly and record the result.",
	2);


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

	}

	void SplittingFunction::persistentOutput(PersistentOStream & os) const {
	using namespace ShowerInteraction;
	os << oenum(_interactionType) << _interactionOrder
	<< oenum(_colourStructure) << _colourFactor
	<< angularOrdered_ << _splittingColourMethod;
	}

	void SplittingFunction::persistentInput(PersistentIStream & is, int) {
	using namespace ShowerInteraction;
	is >> ienum(_interactionType) >> _interactionOrder
	>> ienum(_colourStructure) >> _colourFactor
	>> angularOrdered_ >> _splittingColourMethod;
	}

	void SplittingFunction::colourConnection(tShowerParticlePtr parent,
	tShowerParticlePtr first,
	tShowerParticlePtr second,
	ShowerPartnerType::Type 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 {
	assert(false);
	}
	}

	void SplittingFunction::doinit() {
	Interfaced::doinit();
	assert(_interactionType!=ShowerInteraction::UNDEFINED);
	assert((_colourStructure>0&&_interactionType==ShowerInteraction::QCD) \|\|
	(_colourStructure<0&&_interactionType==ShowerInteraction::QED) );
	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 {
	tcPDPtr pd[3]={getParticleData(ids[0]),
	getParticleData(ids[1]),
	getParticleData(ids[2])};
	if(_colourStructure==TripletTripletOctet) {
	if(ids[0]!=ids[1]) return false;
	if((pd[0]->iColour()==PDT::Colour3\|\|pd[0]->iColour()==PDT::Colour3bar) &&
	pd[2]->iColour()==PDT::Colour8) return true;
	return false;
	}
	else if(_colourStructure==OctetOctetOctet) {
	for(unsigned int ix=0;ix<3;++ix) {
	if(pd[ix]->iColour()!=PDT::Colour8) return false;
	}
	return true;
	}
	else if(_colourStructure==OctetTripletTriplet) {
	if(pd[0]->iColour()!=PDT::Colour8) return false;
	if(pd[1]->iColour()==PDT::Colour3&&pd[2]->iColour()==PDT::Colour3bar)
	return true;
	if(pd[1]->iColour()==PDT::Colour3bar&&pd[2]->iColour()==PDT::Colour3)
	return true;
	return false;
	}
	else if(_colourStructure==TripletOctetTriplet) {
	if(ids[0]!=ids[2]) return false;
	if((pd[0]->iColour()==PDT::Colour3\|\|pd[0]->iColour()==PDT::Colour3bar) &&
	pd[1]->iColour()==PDT::Colour8) return true;
	return false;
	}
	else if(_colourStructure==SextetSextetOctet) {
	if(ids[0]!=ids[1]) return false;
	if((pd[0]->iColour()==PDT::Colour6 \|\| pd[0]->iColour()==PDT::Colour6bar) &&
	pd[2]->iColour()==PDT::Colour8) return true;
	return false;
	}
	else if(_colourStructure==ChargedChargedNeutral) {
	if(ids[0]!=ids[1]) return false;
	if(pd[2]->iCharge()!=0) return false;
	if(pd[0]->iCharge()==pd[1]->iCharge()) return true;
	return false;
	}
	else if(_colourStructure==ChargedNeutralCharged) {
	if(ids[0]!=ids[2]) return false;
	if(pd[1]->iCharge()!=0) return false;
	if(pd[0]->iCharge()==pd[2]->iCharge()) return true;
	return false;
	}
	else if(_colourStructure==NeutralChargedCharged) {
	if(ids[1]!=-ids[2]) return false;
	if(pd[0]->iCharge()!=0) return false;
	if(pd[1]->iCharge()==-pd[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::Type 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);
	// 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;
	}
	// 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;
	}
	// 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;
	}
	}
	}
	// QCD
	else {
	// 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);
	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;
	}
	// g -> q qbar
	else {
	// emitter
	if(emitter->dataPtr()->charged()) {
	emitter->scales().QED = zEmitter*scale;
	emitter->scales().QED_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().QCD_c = zEmitted*scale;
	emitted->scales().QCD_c_noAO = scale;
	emitted->scales().QCD_ac = zEmitted*scale;
	emitted->scales().QCD_ac_noAO = scale;
	}
	}
	}

	-void SplittingFunction::evaluateInitialStateScales(ShowerPartnerType::Type type,
	+void SplittingFunction::evaluateInitialStateScales(ShowerPartnerType::Type partnerType,
	Energy scale, double z,
	tShowerParticlePtr parent,
	- tShowerParticlePtr first,
	- tShowerParticlePtr second) {
	-
	-
	-
	-
	-
	- // emitter->scales().QED = ;
	- // emitter->scales().QED_noAO = ;
	- // emitter->scales().QCD_c = ;
	- // emitter->scales().QCD_c_noAO = ;
	- // emitter->scales().QCD_ac = ;
	- // emitter->scales().QCD_ac_noAO = ;
	- // // emitted
	- // emitted->scales().QED = ;
	- // emitted->scales().QED_noAO = ;
	- // emitted->scales().QCD_c = ;
	- // emitted->scales().QCD_c_noAO = ;
	- // emitted->scales().QCD_ac = ;
	- // emitted->scales().QCD_ac_noAO = ;
	-
	-
	-
	- // // scale for time-like child
	- // Energy AOScale = (angularOrder ? (1.-z()) : 1. )*scale;
	- // // scales for parent if same as spacelike child
	- // if(parent->id()==children[0]->id()) {
	- // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	- // it = children[0]->evolutionScales().begin();
	- // it!=children[0]->evolutionScales().end();++it) {
	- // parent->evolutionScale(it->first,make_pair(min(scale,it->second.first ),
	- // min(scale,it->second.second)));
	- // }
	- // if(partnerType==ShowerPartnerType::QED) {
	- // children[1]->evolutionScale(partnerType,make_pair(AOScale,scale));
	- // }
	- // PDT::Colour childColour = children[1]->dataPtr()->iColour();
	- // if(childColour==PDT::Colour3\|\|childColour==PDT::Colour8\|\|childColour==PDT::Colour6) {
	- // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(AOScale,scale));
	- // }
	- // if(childColour==PDT::Colour3bar\|\|childColour==PDT::Colour8\|\|childColour==PDT::Colour6bar) {
	- // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(AOScale,scale));
	- // }
	- // }
	- // // scales if parent same as timelike child
	- // else if(parent->id()==children[1]->id()) {
	- // if(parent->dataPtr()->charged()) {
	- // parent ->evolutionScale(ShowerPartnerType::QED,make_pair(scale,scale));
	- // children[1]->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale,scale));
	- // }
	- // PDT::Colour childColour = children[1]->dataPtr()->iColour();
	- // if(partnerType==ShowerPartnerType::QED) {
	- // if(childColour==PDT::Colour3\|\|childColour==PDT::Colour8\|\|childColour==PDT::Colour6) {
	- // parent ->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(scale,scale));
	- // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(AOScale,scale));
	- // }
	- // if(childColour==PDT::Colour3bar\|\|childColour==PDT::Colour8\|\|childColour==PDT::Colour6bar) {
	- // parent ->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(scale,scale));
	- // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(AOScale,scale));
	- // }
	- // }
	- // else {
	- // if(childColour==PDT::Colour3\|\|childColour==PDT::Colour8\|\|childColour==PDT::Colour6) {
	- // pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDColourLine);
	- // parent ->evolutionScale(ShowerPartnerType::QCDColourLine,
	- // make_pair(min(scale,pScale.first ),
	- // min(scale,pScale.second)));
	- // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(AOScale,scale));
	- // }
	- // if(childColour==PDT::Colour3bar\|\|childColour==PDT::Colour8\|\|childColour==PDT::Colour6bar) {
	- // pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDAntiColourLine);
	- // parent ->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
	- // make_pair(min(scale,pScale.first ),
	- // min(scale,pScale.second)));
	- // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(AOScale,scale));
	- // }
	- // }
	- // }
	- // // g -> q qbar, or gamma -> e+e-
	- // else if(children[0]->id()==-children[1]->id()) {
	- // PDT::Colour childColour = children[1]->dataPtr()->iColour();
	- // if(partnerType==ShowerPartnerType::QED) {
	- // parent ->evolutionScale(ShowerPartnerType::QED,make_pair(scale,scale));
	- // children[1]->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale,scale));
	- // if(childColour==PDT::Colour3\|\|childColour==PDT::Colour8\|\|childColour==PDT::Colour6) {
	- // pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDAntiColourLine);
	- // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,
	- // make_pair(min(AOScale,pScale.first ),
	- // min(AOScale,pScale.second)));
	- // }
	- // if(childColour==PDT::Colour3bar\|\|childColour==PDT::Colour8\|\|childColour==PDT::Colour6bar) {
	- // pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDColourLine );
	- // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
	- // make_pair(min(AOScale,pScale.first ),
	- // min(AOScale,pScale.second)));
	- // }
	- // }
	- // else {
	- // pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QED);
	- // children[1]->evolutionScale(ShowerPartnerType::QED,
	- // make_pair(min(AOScale,pScale.first ),
	- // min(AOScale,pScale.second)));
	- // if(childColour==PDT::Colour3\|\|childColour==PDT::Colour8\|\|childColour==PDT::Colour6) {
	- // pScale = children[0]->evolutionScale(ShowerPartnerType::QCDAntiColourLine);
	- // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,
	- // make_pair(min(AOScale,pScale.first ),
	- // min(AOScale,pScale.second)));
	- // }
	- // if(childColour==PDT::Colour3bar\|\|childColour==PDT::Colour8\|\|childColour==PDT::Colour6bar) {
	- // pScale = children[0]->evolutionScale(ShowerPartnerType::QCDColourLine );
	- // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
	- // make_pair(min(AOScale,pScale.first ),
	- // min(AOScale,pScale.second)));
	- // }
	- // parent->evolutionScale(ShowerPartnerType::QCDColourLine,
	- // make_pair(min(scale,pScale.first ),
	- // min(scale,pScale.second)));
	- // parent->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
	- // make_pair(min(scale,pScale.first ),
	- // min(scale,pScale.second)));
	- // }
	- // }
	- assert(false);
	+ 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 {
	+ // 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;
	+ }
	+ }
	+ }
	}

	-void SplittingFunction::evaluateDecayScales(ShowerPartnerType::Type type,
	+void SplittingFunction::evaluateDecayScales(ShowerPartnerType::Type partnerType,
	Energy scale, double z,
	tShowerParticlePtr parent,
	- tShowerParticlePtr first,
	- tShowerParticlePtr second) {
	- // // angular-ordered scale for 2nd child
	- // Energy AOScale = angularOrder ? (1.-z())*scale : scale;
	- // // QED
	- // if(partnerType==ShowerPartnerType::QED) {
	- // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	- // it = parent->evolutionScales().begin();
	- // it!=parent->evolutionScales().end();++it) {
	- // children[0]->evolutionScale(it->first,make_pair(min(scale,it->second.first ),
	- // min(scale,it->second.second)));
	- // if(it->first==ShowerPartnerType::QED) {
	- // children[1]->evolutionScale(it->first,make_pair(min(AOScale,it->second.first ),
	- // min(scale,it->second.second)));
	- // }
	- // }
	- // }
	- // // QCD
	- // else {
	- // // scales for the emitter
	- // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	- // it = parent->evolutionScales().begin();
	- // it!=parent->evolutionScales().end();++it) {
	- // children[0]->evolutionScale(it->first,make_pair(min(scale,it->second.first ),
	- // min(scale,it->second.second)));
	- // }
	- // PDT::Colour emittedColour = children[1]->dataPtr()->iColour();
	- // if(emittedColour==PDT::Colour3\|\|emittedColour==PDT::Colour8\|\|emittedColour==PDT::Colour6) {
	- // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,
	- // make_pair(AOScale,scale));
	- // }
	- // if(emittedColour==PDT::Colour3bar\|\|emittedColour==PDT::Colour8\|\|emittedColour==PDT::Colour6bar) {
	- // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
	- // make_pair(AOScale,scale));
	- // }
	- // if(children[1]->dataPtr()->charged()) {
	- // children[1]->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale,scale));
	- // }
	- // }
	- assert(false);
	+ 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;
	+ // spacelike
	+ spacelike->scales().QED = scale;
	+ spacelike->scales().QED_noAO = scale;
	+ }
	+ // QCD
	+ else {
	+ // 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;
	+ // spacelike
	+ spacelike->scales().QED = max(scale,parent->scales().QED );
	+ spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO );
	+ }
	+ 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/SplittingFunctions/SplittingGenerator.cc b/Shower/SplittingFunctions/SplittingGenerator.cc
	--- a/Shower/SplittingFunctions/SplittingGenerator.cc
	+++ b/Shower/SplittingFunctions/SplittingGenerator.cc
	@@ -1,555 +1,546 @@
	// -- C++ --
	//
	// SplittingGenerator.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the 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/Reference.h"
	#include "ThePEG/Interface/Command.h"
	#include "ThePEG/Utilities/StringUtils.h"
	#include "ThePEG/Repository/Repository.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "Herwig++/Shower/ShowerHandler.h"
	#include "ThePEG/Utilities/Rebinder.h"
	#include <cassert>

	using namespace Herwig;

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

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


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

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

	ClassDescription<SplittingGenerator> SplittingGenerator::initSplittingGenerator;
	// Definition of the static class description member.

	void SplittingGenerator::Init() {

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

	static Switch<SplittingGenerator, bool> interfaceISRMode
	("ISR",
	"Include initial-state radiation?",
	&SplittingGenerator::_isr_Mode, 1, false, false);
	static SwitchOption interfaceISRMode0
	(interfaceISRMode,"No","ISR (Initial State Radiation) is OFF", 0);
	static SwitchOption interfaceISRMode1
	(interfaceISRMode,"Yes","ISR (Initial State Radiation) is ON", 1);

	static Switch<SplittingGenerator, bool> interfaceFSRMode
	("FSR",
	"Include final-state radiation?",
	&SplittingGenerator::_fsr_Mode, 1, false, false);
	static SwitchOption interfaceFSRMode0
	(interfaceFSRMode,"No","FSR (Final State Radiation) is OFF", 0);
	static SwitchOption interfaceFSRMode1
	(interfaceFSRMode,"Yes","FSR (Final State Radiation) is ON", 1);

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

	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->id());
	for(vector<tPDPtr>::iterator it = products.begin(); it!=products.end(); ++it)
	ids.push_back((*it)->id());
	// check splitting can handle this
	if(!s->splittingFn()->accept(ids))
	return "Error: Sudakov " + sudakov + "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->id());
	for(vector<tPDPtr>::iterator it = products.begin(); it!=products.end(); ++it)
	ids.push_back((*it)->id());
	// check splitting can handle this
	if(!s->splittingFn()->accept(ids))
	return "Error: Sudakov " + sudakov + "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(isISRadiationON() && !final) {
	_bbranchings.insert(BranchingInsert(ids[1],BranchingElement(s,ids)));
	s->addSplitting(ids);
	}
	if(isFSRadiationON() && final) {
	_fbranchings.insert(BranchingInsert(ids[0],BranchingElement(s,ids)));
	s->addSplitting(ids);
	}
	}

	void SplittingGenerator::deleteFromMap(const IdList &ids,
	const SudakovPtr &s, bool final) {
	if(isISRadiationON() && !final) {
	pair<BranchingList::iterator,BranchingList::iterator>
	range = _bbranchings.equal_range(ids[1]);
	for(BranchingList::iterator it=range.first;it!=range.second&&it->first==ids[1];++it) {
	if(it->second.first==s&&it->second.second==ids)
	_bbranchings.erase(it);
	}
	s->removeSplitting(ids);
	}
	if(isFSRadiationON() && final) {
	pair<BranchingList::iterator,BranchingList::iterator>
	range = _fbranchings.equal_range(ids[0]);
	for(BranchingList::iterator it=range.first;it!=range.second&&it->first==ids[0];++it) {
	if(it->second.first==s&&it->second.second==ids)
	_fbranchings.erase(it);
	}
	s->removeSplitting(ids);
	}
	}

	Branching SplittingGenerator::chooseForwardBranching(ShowerParticle &particle,
	double enhance,
	ShowerInteraction::Type type) const {
	Energy newQ = ZERO;
	ShoKinPtr kinematics = ShoKinPtr();
	ShowerPartnerType::Type 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(type != cit->second.first->interactionType() &&
	type != ShowerInteraction::Both ) continue;
	// whether or not this interaction should be angular ordered
	bool angularOrdered = cit->second.first->splittingFn()->angularOrdered();
	ShoKinPtr newKin;
	ShowerPartnerType::Type type;
	// work out which starting scale we need
	if(cit->second.first->interactionType()==ShowerInteraction::QED) {
	type = ShowerPartnerType::QED;
	Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
	newKin = cit->second.first->
	generateNextTimeBranching(startingScale,cit->second.second,
	particle.id()!=cit->first,enhance);
	}
	else if(cit->second.first->interactionType()==ShowerInteraction::QCD) {
	// special for octets
	if(particle.dataPtr()->iColour()==PDT::Colour8) {
	// octet -> octet octet
	if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) {
	type = ShowerPartnerType::QCDColourLine;
	Energy startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
	newKin= cit->second.first->
	generateNextTimeBranching(startingScale,cit->second.second,
	particle.id()!=cit->first,0.5*enhance);
	startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
	ShoKinPtr newKin2 = cit->second.first->
	generateNextTimeBranching(startingScale,cit->second.second,
	particle.id()!=cit->first,0.5*enhance);
	// 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_c_noAO);
	newKin= cit->second.first->
	generateNextTimeBranching(startingScale, cit->second.second,
	particle.id()!=cit->first,enhance);
	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.first->
	generateNextTimeBranching(startingScale,cit->second.second,
	particle.id()!=cit->first,enhance);
	}
	}
	// 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 = cit->second.second;
	sudakov = cit->second.first;
	partnerType = type;
	}
	}
	// return empty branching if nothing happened
	if(!kinematics) return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
	ShowerPartnerType::Undefined);
	// If a branching has been selected initialize it
	kinematics->initialize(particle,PPtr());
	// and return it
	return Branching(kinematics, ids,sudakov,partnerType);
	}

	Branching SplittingGenerator::
	chooseDecayBranching(ShowerParticle &particle,
	- const map<ShowerPartnerType::Type,pair<Energy,Energy> > & stoppingScales,
	+ const ShowerParticle::EvolutionScales & stoppingScales,
	Energy minmass, double enhance,
	ShowerInteraction::Type interaction) const {
	Energy newQ = Constants::MaxEnergy;
	ShoKinPtr kinematics;
	SudakovPtr sudakov;
	ShowerPartnerType::Type 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.second[1]!=index&&cit->second.second[2]!=index) continue;
	// check either right interaction or doing both
	if(interaction != cit->second.first->interactionType() &&
	interaction != ShowerInteraction::Both ) continue;
	// whether or not this interaction should be angular ordered
	bool angularOrdered = cit->second.first->splittingFn()->angularOrdered();
	ShoKinPtr newKin;
	ShowerPartnerType::Type type;
	// work out which starting scale we need
	- // if(cit->second.first->interactionType()==ShowerInteraction::QED) {
	- // type = ShowerPartnerType::QED;
	- // map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	- // it=stoppingScales.find(type);
	- // if(it==stoppingScales.end()) continue;
	- // Energy stoppingScale = angularOrdered ? it->second.first : it->second.second;
	- // if(particle.evolutionScale(angularOrdered,type) < stoppingScale ) {
	- // newKin = cit->second.first->
	- // generateNextDecayBranching(particle.evolutionScale(angularOrdered,type),
	- // stoppingScale,minmass,cit->second.second,
	- // particle.id()!=cit->first,enhance);
	- // cerr << "testing in QED " << particle << "\n";
	- // assert(false);
	- // }
	- // }
	- // else if(cit->second.first->interactionType()==ShowerInteraction::QCD) {
	- // // special for octets
	- // if(particle.dataPtr()->iColour()==PDT::Colour8) {
	- // map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	- // it=stoppingScales.find(ShowerPartnerType::QCDColourLine);
	- // if(it==stoppingScales.end()) continue;
	- // Energy stoppingColour = angularOrdered ? it->second.first : it->second.second;
	- // it=stoppingScales.find(ShowerPartnerType::QCDAntiColourLine);
	- // Energy stoppingAnti = angularOrdered ? it->second.first : it->second.second;
	- // // octet -> octet octet
	- // if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) {
	- // type = ShowerPartnerType::QCDColourLine;
	- // newKin= cit->second.first->
	- // generateNextDecayBranching(particle.evolutionScale(angularOrdered,type),
	- // stoppingColour,minmass,
	- // cit->second.second,
	- // particle.id()!=cit->first,0.5*enhance);
	- // ShoKinPtr newKin2 = cit->second.first->
	- // generateNextDecayBranching(particle.evolutionScale(angularOrdered,
	- // ShowerPartnerType::QCDAntiColourLine),
	- // stoppingAnti,minmass,
	- // cit->second.second,
	- // particle.id()!=cit->first,0.5*enhance);
	- // // pick the one with the lowest scale
	- // if( (newKin&&newKin2&&newKin2->scale()<newKin->scale()) \|\|
	- // (!newKin&&newKin2) ) {
	- // newKin = newKin2;
	- // type = ShowerPartnerType::QCDAntiColourLine;
	- // }
	- // }
	- // // other
	- // else {
	- // Energy startingScale =
	- // min(particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDColourLine),
	- // particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDAntiColourLine));
	- // newKin = cit->second.first->
	- // generateNextDecayBranching(startingScale,
	- // max(stoppingColour,stoppingAnti),minmass,
	- // cit->second.second,
	- // particle.id()!=cit->first,enhance);
	- // type = cit->second.second[0]<0 ?
	- // ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
	- // }
	- // }
	- // // everything else
	- // else {
	- // type = particle.hasColour() ?
	- // ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
	- // map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	- // it=stoppingScales.find(type);
	- // if(it==stoppingScales.end()) continue;
	- // Energy stoppingScale = angularOrdered ? it->second.first : it->second.second;
	- // if(particle.evolutionScale(angularOrdered,type) < stoppingScale ) {
	- // newKin = cit->second.first->
	- // generateNextDecayBranching(particle.evolutionScale(angularOrdered,type),
	- // stoppingScale,minmass,cit->second.second,
	- // particle.id()!=cit->first,enhance);
	- // }
	- // }
	- // }
	- // // shouldn't be anything else
	- // else
	- // assert(false);
	- assert(false);
	+ if(cit->second.first->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.first->
	+ generateNextDecayBranching(startingScale,stoppingScale,minmass,cit->second.second,
	+ particle.id()!=cit->first,enhance);
	+ }
	+ }
	+ else if(cit->second.first->interactionType()==ShowerInteraction::QCD) {
	+ // special for octets
	+ if(particle.dataPtr()->iColour()==PDT::Colour8) {
	+ // octet -> octet octet
	+ if(cit->second.first->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.first->
	+ generateNextDecayBranching(startingColour,stoppingColour,minmass,
	+ cit->second.second,
	+ particle.id()!=cit->first,0.5*enhance);
	+ }
	+ ShoKinPtr newKin2;
	+ if(startingAnti<stoppingAnti) {
	+ newKin2 = cit->second.first->
	+ generateNextDecayBranching(startingAnti,stoppingAnti,minmass,
	+ cit->second.second,
	+ particle.id()!=cit->first,0.5*enhance);
	+ }
	+ // 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.first->
	+ generateNextDecayBranching(startingScale,stoppingScale,minmass,cit->second.second,
	+ particle.id()!=cit->first,enhance);
	+ }
	+ }
	+ }
	+ // shouldn't be anything else
	+ else
	+ assert(false);
	if(!newKin) continue;
	// select highest scale
	if(newKin->scale() < newQ ) {
	newQ = newKin->scale();
	ids = cit->second.second;
	kinematics=newKin;
	sudakov=cit->second.first;
	partnerType = type;
	}
	}
	// return empty branching if nothing happened
	if(!kinematics) return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
	ShowerPartnerType::Undefined);
	// initialize the branching
	kinematics->initialize(particle,PPtr());
	// and return it
	return Branching(kinematics, ids,sudakov,partnerType);
	}

	Branching SplittingGenerator::
	chooseBackwardBranching(ShowerParticle &particle,PPtr beamparticle,
	double enhance,
	Ptr<BeamParticleData>::transient_const_pointer beam,
	ShowerInteraction::Type type,
	tcPDFPtr pdf, Energy freeze) const {
	Energy newQ=ZERO;
	ShoKinPtr kinematics=ShoKinPtr();
	ShowerPartnerType::Type 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(type != cit->second.first->interactionType() &&
	type != ShowerInteraction::Both ) continue;
	// setup the PDF
	cit->second.first->setPDF(pdf,freeze);
	// whether or not this interaction should be angular ordered
	bool angularOrdered = cit->second.first->splittingFn()->angularOrdered();
	ShoKinPtr newKin;
	ShowerPartnerType::Type type;
	- // if(cit->second.first->interactionType()==ShowerInteraction::QED) {
	- // type = ShowerPartnerType::QED;
	- // newKin=cit->second.first->
	- // generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
	- // cit->second.second, particle.x(),
	- // particle.id()!=cit->first,enhance,beam);
	- // cerr << "testing in QED " << particle << "\n";
	- // assert(false);
	- // }
	- // else if(cit->second.first->interactionType()==ShowerInteraction::QCD) {
	- // // special for octets
	- // if(particle.dataPtr()->iColour()==PDT::Colour8) {
	- // // octet -> octet octet
	- // if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) {
	- // type = ShowerPartnerType::QCDColourLine;
	- // newKin=cit->second.first->
	- // generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
	- // cit->second.second, particle.x(),
	- // particle.id()!=cit->first,0.5*enhance,beam);
	- // ShoKinPtr newKin2 = cit->second.first->
	- // generateNextSpaceBranching(particle.evolutionScale(angularOrdered,
	- // ShowerPartnerType::QCDAntiColourLine),
	- // cit->second.second, particle.x(),
	- // particle.id()!=cit->first,0.5*enhance,beam);
	- // // pick the one with the highest scale
	- // if( (newKin&&newKin2&&newKin2->scale()>newKin->scale()) \|\|
	- // (!newKin&&newKin2) ) {
	- // newKin = newKin2;
	- // type = ShowerPartnerType::QCDAntiColourLine;
	- // }
	- // }
	- // // other
	- // else {
	- // Energy startingScale =
	- // max(particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDColourLine),
	- // particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDAntiColourLine));
	- // type = UseRandom::rndbool() ?
	- // ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
	- // newKin=cit->second.first->
	- // generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
	- // cit->second.second, particle.x(),
	- // particle.id()!=cit->first,enhance,beam);
	- // }
	- // }
	- // // everything else
	- // else {
	- // type = particle.hasColour() ?
	- // ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
	- // newKin=cit->second.first->
	- // generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
	- // cit->second.second, particle.x(),
	- // particle.id()!=cit->first,enhance,beam);
	- // }
	- // }
	- // // shouldn't be anything else
	- // else
	- // assert(false);
	- assert(false);
	+ if(cit->second.first->interactionType()==ShowerInteraction::QED) {
	+ type = ShowerPartnerType::QED;
	+ Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
	+ newKin=cit->second.first->
	+ generateNextSpaceBranching(startingScale,cit->second.second, particle.x(),
	+ particle.id()!=cit->first,enhance,beam);
	+ }
	+ else if(cit->second.first->interactionType()==ShowerInteraction::QCD) {
	+ // special for octets
	+ if(particle.dataPtr()->iColour()==PDT::Colour8) {
	+ // octet -> octet octet
	+ if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) {
	+ type = ShowerPartnerType::QCDColourLine;
	+ Energy startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
	+ newKin = cit->second.first->
	+ generateNextSpaceBranching(startingScale,cit->second.second, particle.x(),
	+ particle.id()!=cit->first,0.5*enhance,beam);
	+ startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
	+ ShoKinPtr newKin2 = cit->second.first->
	+ generateNextSpaceBranching(startingScale,cit->second.second, particle.x(),
	+ particle.id()!=cit->first,0.5*enhance,beam);
	+ // 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_c_noAO);
	+ type = UseRandom::rndbool() ?
	+ ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
	+ newKin=cit->second.first->
	+ generateNextSpaceBranching(startingScale,cit->second.second, particle.x(),
	+ particle.id()!=cit->first,enhance,beam);
	+ }
	+ }
	+ // 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.first->
	+ generateNextSpaceBranching(startingScale,cit->second.second, particle.x(),
	+ particle.id()!=cit->first,enhance,beam);
	+ }
	+ }
	+ // 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 = cit->second.second;
	sudakov=cit->second.first;
	partnerType = type;
	}
	}
	// return empty branching if nothing happened
	if(!kinematics) return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
	ShowerPartnerType::Undefined);
	// initialize the ShowerKinematics
	// and return it
	kinematics->initialize(particle,beamparticle);
	// 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).first=trans.translate((cit->second).first);}
	for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit)
	{(cit->second).first=trans.translate((cit->second).first);}
	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).first);}
	for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit)
	{ret.push_back((cit->second).first);}
	return ret;
	}

	diff --git a/Shower/SplittingFunctions/SplittingGenerator.h b/Shower/SplittingFunctions/SplittingGenerator.h
	--- a/Shower/SplittingFunctions/SplittingGenerator.h
	+++ b/Shower/SplittingFunctions/SplittingGenerator.h
	@@ -1,377 +1,378 @@
	// -- C++ --
	//
	// SplittingGenerator.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_SplittingGenerator_H
	#define HERWIG_SplittingGenerator_H
	//
	// This is the declaration of the SplittingGenerator class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "Herwig++/Shower/Base/Branching.h"
	#include "Herwig++/Shower/Base/SudakovFormFactor.h"
	#include "SplittingGenerator.fh"
	+#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "Herwig++/Shower/Base/ShowerKinematics.h"

	namespace Herwig {

	using namespace ThePEG;

	/** \ingroup Shower
	*
	* This class is responsible for creating, at the beginning of the Run,
	* all the SplittingFunction objects and the corresponding
	* SudakovFormFactor objects, and then of the generation of splittings
	* (radiation emissions) during the event.
	* Many switches are defined in this class which allowed the user to turn on/off:
	* - each type of interaction (QCD, QED, EWK,...);
	* - initial- and final-state radiation for all type of interactions;
	* - initial- and final-state radiation for each type of interaction;
	* - each type of splitting (\f$u\to ug\f$, \f$d\to dg\f$, \f$\ldots\f$,
	* \f$g\to gg\f$, \f$g\to u\bar{u}\f$, \f$\ldots\f$).
	*
	* These switches are useful mainly for debugging, but eventually can
	* also be used for a "quick and dirty" estimation of systematic errors.
	*
	* In the future it should be possible to implement in this class
	*
	* - the \f$1\to2\f$ azimuthal correlations for soft emission due to QCD coherence
	* using the ShowerParticle object provided in the input.
	* - Similarly having the \f$\rho-D\f$ matrix and the SplittingFunction pointer
	* it should be possible to implement the spin correlations.
	*
	* @see SudakovFormFactor
	* @see SplitFun
	*
	* @see \ref SplittingGeneratorInterfaces "The interfaces"
	* defined for SplittingGenerator.
	*/
	class SplittingGenerator: public Interfaced {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	SplittingGenerator() : _isr_Mode(1), _fsr_Mode(1) {}
	//@}

	public:

	/**
	* Methods to select the next branching and reconstruct the kinematics
	*/
	//@{
	/**
	* Choose a new forward branching for a time-like particle
	* The method returns:
	* - a pointer to a ShowerKinematics object, which
	* contains the information about the new scale and all other
	* kinematics variables that need to be generated simultaneously;
	* - a pointer to the SudakovFormFactor object associated
	* with the chosen emission.
	* - The PDG codes of the particles in the branching,
	* as a Branching struct.
	*
	* In the case no branching has been generated, both the returned
	* pointers are null ( ShoKinPtr() , tSudakovFFPtr() ).
	*
	* @param particle The particle to be evolved
	* @param enhance The factor by which to ehnace the emission of radiation
	* @param type The type of interaction to generate
	* @return The Branching struct for the branching
	*/
	Branching chooseForwardBranching(ShowerParticle & particle,
	double enhance,
	ShowerInteraction::Type type) const;

	/**
	* Select the next branching of a particles for the initial-state shower
	* in the particle's decay.
	* @param particle The particle being showerwed
	* @param maxscale The maximum scale
	* @param minmass Minimum mass of the particle after the branching
	* @param enhance The factor by which to ehnace the emission of radiation
	* @param type The type of interaction to generate
	* @return The Branching struct for the branching
	*/
	Branching chooseDecayBranching(ShowerParticle & particle,
	- const map<ShowerPartnerType::Type,pair<Energy,Energy> > & maxScales,
	+ const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass,double enhance,
	ShowerInteraction::Type type) const;

	/**
	* Choose a new backward branching for a space-like particle.
	* The method returns:
	* - a pointer to a ShowerKinematics object, which
	* contains the information about the new scale and all other
	* kinematics variables that need to be generated simultaneously;
	* - a pointer to the SudakovFormFactor object associated
	* with the chosen emission.
	* - The PDG codes of the particles in the branching,
	* as a Branching struct.
	*
	* In the case no branching has been generated, both the returned
	* pointers are null ( ShoKinPtr() , tSudakovFFPtr() ).
	*
	* @param particle The particle to be evolved
	* @param enhance The factor by which to ehnace the emission of radiation
	* @param beamparticle The beam particle
	* @param beam The BeamParticleData object
	* @param type The type of interaction to generate
	* @return The Branching struct for the branching
	*/
	Branching
	chooseBackwardBranching(ShowerParticle & particle,
	PPtr beamparticle,
	double enhance,
	Ptr<BeamParticleData>::transient_const_pointer beam,
	ShowerInteraction::Type type,
	tcPDFPtr , Energy ) const;
	//@}

	public:

	/**
	* Access to the switches
	*/
	//@{
	/**
	* It returns true/false if the initial-state radiation is on/off.
	*/
	bool isISRadiationON() const { return _isr_Mode; }

	/**
	* It returns true/false if the final-state radiation is on/off.
	*/
	bool isFSRadiationON() const { return _fsr_Mode; }
	//@}

	/**
	* Methods to parse the information from the input files to create the
	* branchings
	*/
	//@{
	/**
	* Add a final-state splitting
	*/
	string addFinalSplitting(string arg) { return addSplitting(arg,true); }

	/**
	* Add an initial-state splitting
	*/
	string addInitialSplitting(string arg) { return addSplitting(arg,false); }

	/**
	* Add a final-state splitting
	*/
	string deleteFinalSplitting(string arg) { return deleteSplitting(arg,true); }

	/**
	* Add an initial-state splitting
	*/
	string deleteInitialSplitting(string arg) { return deleteSplitting(arg,false); }
	//@}

	/**
	* Access to the splittings
	*/
	//@{
	/**
	* Access the final-state branchings
	*/
	const BranchingList & finalStateBranchings() const { return _fbranchings; }

	/**
	* Access the initial-state branchings
	*/
	const BranchingList & initialStateBranchings() const { return _bbranchings; }
	//@}

	public:

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

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

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

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	protected:

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Rebind pointer to other Interfaced objects. Called in the setup phase
	* after all objects used in an EventGenerator has been cloned so that
	* the pointers will refer to the cloned objects afterwards.
	* @param trans a TranslationMap relating the original objects to
	* their respective clones.
	* @throws RebindException if no cloned object was found for a given
	* pointer.
	*/
	virtual void rebind(const TranslationMap & trans)
	;

	/**
	* Return a vector of all pointers to Interfaced objects used in this
	* object.
	* @return a vector of pointers.
	*/
	virtual IVector getReferences();
	//@}

	private:

	/**
	* Add a branching to the map
	* @param ids PDG coeds of the particles in the branching
	* @param sudakov The SudakovFormFactor for the branching
	* @param final Whether this is an initial- or final-state branching
	*/
	void addToMap(const IdList & ids, const SudakovPtr & sudakov, bool final);

	/**
	* Remove a branching to the map
	* @param ids PDG coeds of the particles in the branching
	* @param sudakov The SudakovFormFactor for the branching
	* @param final Whether this is an initial- or final-state branching
	*/
	void deleteFromMap(const IdList & ids, const SudakovPtr & sudakov, bool final);

	/**
	* Obtain the reference vectors for a final-state particle
	* @param particle The particle
	* @param p The p reference vector
	* @param n The n reference vector
	*/
	void finalStateBasisVectors(ShowerParticle particle, Lorentz5Momentum & p,
	Lorentz5Momentum & n) const;

	/**
	* Add a splitting
	* @param in string to be parsed
	* @param final Whether this is an initial- or final-state branching
	*/
	string addSplitting(string in ,bool final);

	/**
	* Delete a splitting
	* @param in string to be parsed
	* @param final Whether this is an initial- or final-state branching
	*/
	string deleteSplitting(string in ,bool final);

	private:

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

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

	private:

	/**
	* Switches to control the radiation
	*/
	//@{
	/**
	* Is inqitial-state radiation on/off
	*/
	bool _isr_Mode;

	/**
	* Is final-state radiation on/off
	*/
	bool _fsr_Mode;
	//@}

	/**
	* List of the branchings and the appropriate Sudakovs for forward branchings
	*/
	BranchingList _fbranchings;

	/**
	* Lists of the branchings and the appropriate Sudakovs for backward branchings.
	*/
	BranchingList _bbranchings;
	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	/** @endcond */

	}

	#endif /* HERWIG_SplittingGenerator_H */

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 5:23 PM (14 h, 48 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023620
Default Alt Text: (245 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions