No OneTemporary
Actions

Size

167 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/DrellYanBase.cc b/MatrixElement/DrellYanBase.cc
	--- a/MatrixElement/DrellYanBase.cc
	+++ b/MatrixElement/DrellYanBase.cc
	@@ -1,1047 +1,1048 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the DrellYanBase class.
	//

	#include "DrellYanBase.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/EnumParticles.h"

	using namespace Herwig;

	DrellYanBase::DrellYanBase()
	: _channelwgtA(0.12), _channelwgtB(2.00), _nover(0), _maxwgt(0.),
	_power(2.0),_preqqbar(6.5),_preqg(4.0),_pregqbar(4.0),
	_min_pt(2.*GeV) {}

	void DrellYanBase::persistentOutput(PersistentOStream & os) const {
	os << _channelwgtA << _channelwgtB << _channelweights << _alpha
	<< _power << _preqqbar << _preqg << _pregqbar << ounit( _min_pt, GeV )
	<< _prefactor;
	}

	void DrellYanBase::persistentInput(PersistentIStream & is, int) {
	is >> _channelwgtA >> _channelwgtB >> _channelweights >> _alpha
	>> _power >> _preqqbar >> _preqg >> _pregqbar >> iunit( _min_pt, GeV )
	>> _prefactor;
	}

	AbstractClassDescription<DrellYanBase> DrellYanBase::initDrellYanBase;
	// Definition of the static class description member.

	void DrellYanBase::Init() {

	static ClassDocumentation<DrellYanBase> documentation
	("The DrellYanBase class provides a base class for the"
	" corrections to Drell-Yan type processes");

	static Parameter<DrellYanBase,double> interfaceQQbarChannelWeight
	("QQbarChannelWeight",
	"The relative weights of the q qbar abd q g channels for selection."
	" This is a technical parameter for the phase-space generation and "
	"should not affect the results only the efficiency and fraction"
	" of events with weight > 1.",
	&DrellYanBase::_channelwgtA, 0.12, 0.01, 100.,
	false, false, Interface::limited);

	static Parameter<DrellYanBase,double> interfaceQGChannelWeight
	("QGChannelWeight",
	"The relative weights of the qg abd qbar g channels for selection."
	" This is a technical parameter for the phase-space generation and "
	"should not affect the results only the efficiency and fraction",
	&DrellYanBase::_channelwgtB, 2., 0.01, 100.,
	false, false, Interface::limited);

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

	static Parameter<DrellYanBase,double> interfacePower
	("Power",
	"The power for the sampling of the matrix elements",
	&DrellYanBase::_power, 2.0, 1.0, 10.0,
	false, false, Interface::limited);

	static Parameter<DrellYanBase,double> interfacePrefactorqqbar
	("Prefactorqqbar",
	"The prefactor for the sampling of the q qbar channel",
	&DrellYanBase::_preqqbar, 5.0, 0.0, 1000.0,
	false, false, Interface::limited);

	static Parameter<DrellYanBase,double> interfacePrefactorqg
	("Prefactorqg",
	"The prefactor for the sampling of the q g channel",
	&DrellYanBase::_preqg, 3.0, 0.0, 1000.0,
	false, false, Interface::limited);

	static Parameter<DrellYanBase,double> interfacePrefactorgqbar
	("Prefactorgqbar",
	"The prefactor for the sampling of the g qbar channel",
	&DrellYanBase::_pregqbar, 3.0, 0.0, 1000.0,
	false, false, Interface::limited);

	static Parameter<DrellYanBase, Energy> interfacePtMin
	("minPt",
	"The pt cut on hardest emision generation"
	"2(1-Beta)exp(-sqr(intrinsicpT/RMS))/sqr(RMS)",
	&DrellYanBase::_min_pt, GeV, 2.GeV, ZERO, 100000.0GeV,
	false, false, Interface::limited);
	}

	void DrellYanBase::doinit() {
	HwMEBase::doinit();
	_channelweights.push_back(_channelwgtA/(1.+_channelwgtA));
	_channelweights.push_back(_channelweights[0]+1./(1.+_channelwgtA)/(1+_channelwgtB));
	_channelweights.push_back(1.0);
	_prefactor.push_back(_preqqbar);
	_prefactor.push_back(_preqg);
	_prefactor.push_back(_pregqbar);
	}

	void DrellYanBase::dofinish() {
	HwMEBase::dofinish();
	if(_nover==0) return;
	generator()->log() << "DrellYanBase when applying the hard correction "
	<< _nover << " weights larger than one were generated of which"
	<< " the largest was " << _maxwgt << "\n";
	}

	void DrellYanBase::applyHardMatrixElementCorrection(ShowerTreePtr tree) {
	// get the quark,antiquark and the gauge boson
	// get the quarks
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
	ShowerParticleVector incoming;
	vector<tcBeamPtr> beams;
	for(cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
	incoming.push_back(cit->first->progenitor());
	beams.push_back(cit->first->beam());
	}
	// ensure that the quark is first
	if(incoming[0]->id()<incoming[1]->id()) {
	swap(incoming[0],incoming[1]);
	swap(beams[0],beams[1]);
	}
	Lorentz5Momentum pboson;
	for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt= tree->outgoingLines().begin();
	cjt != tree->outgoingLines().end();++cjt ) {
	pboson += cjt->first->original()->momentum();
	}
	pboson.rescaleMass();
	// calculate the momenta
	unsigned int iemit,itype;
	vector<Lorentz5Momentum> pnew;
	LorentzRotation trans;
	pair<double,double> xnew;
	// if not accepted return
	if(!applyHard(incoming,beams,pboson,iemit,itype,pnew,trans,xnew)) return;
	// if applying ME correction create the new particles
	if(itype==0) {
	// get the momenta of the new particles
	Lorentz5Momentum pquark(pnew[0]),panti(pnew[1]),pgluon(pnew[2]);
	if(iemit==2) swap(pquark,panti);
	// ensure gluon can be put on shell
	Lorentz5Momentum ptest(pgluon);
	if(ptest.boost(-(pquark+panti).boostVector()).e() <
	getParticleData(ParticleID::g)->constituentMass()) return;
	// create the new gluon
	PPtr newg= getParticleData(ParticleID::g)->produceParticle(pgluon);
	PPtr newq,newa;
	ColinePtr col;
	// make the new particles
	if(iemit==1) {
	col=incoming[0]->colourLine();
	newq = getParticleData(incoming[0]->id())->produceParticle(pquark);
	newa = new_ptr(Particle(*incoming[1]));
	col->removeAntiColoured(newa);
	newa->set5Momentum(panti);
	}
	else {
	col=incoming[1]->antiColourLine();
	newa = getParticleData(incoming[1]->id())->produceParticle(panti);
	newq = new_ptr(Particle(*incoming[0]));
	col->removeColoured(newq);
	newq->set5Momentum(pquark);
	}
	// set the colour lines
	ColinePtr newline=new_ptr(ColourLine());
	if(iemit==1) {
	newline->addColoured(newq);
	newline->addColoured(newg);
	col->addAntiColoured(newg);
	col->addAntiColoured(newa);
	}
	else {
	newline->addAntiColoured(newa);
	newline->addAntiColoured(newg);
	col->addColoured(newg);
	col->addColoured(newq);
	}
	// change the existing quark and antiquark
	PPtr orig;
	for(cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
	if(cit->first->progenitor()->id()==newq->id()) {
	// remove old particles from colour line
	col->removeColoured(cit->first->copy());
	col->removeColoured(cit->first->progenitor());
	// insert new particles
	cit->first->copy(newq);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newq,1,false)));
	sp->x(xnew.first);
	cit->first->progenitor(sp);
	tree->incomingLines()[cit->first]=sp;
	cit->first->perturbative(iemit!=1);
	if(iemit==1) orig=cit->first->original();
	}
	else {
	// remove old particles from colour line
	col->removeAntiColoured(cit->first->copy());
	col->removeColoured(cit->first->progenitor());
	// insert new particles
	cit->first->copy(newa);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newa,1,false)));
	sp->x(xnew.second);
	cit->first->progenitor(sp);
	tree->incomingLines()[cit->first]=sp;
	cit->first->perturbative(iemit==1);
	if(iemit==2) orig=cit->first->original();
	}
	}
	// fix the momentum of the gauge boson
	Boost boostv=pboson.findBoostToCM();
	trans *=LorentzRotation(boostv);
	for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt= tree->outgoingLines().begin();
	cjt != tree->outgoingLines().end();++cjt ) {
	cjt->first->progenitor()->transform(trans);
	cjt->first->copy()->transform(trans);
	}
	tree->hardMatrixElementCorrection(true);
	// add the gluon
	ShowerParticlePtr sg=new_ptr(ShowerParticle(*newg,1,true));
	ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(orig,newg,sg));
	gluon->perturbative(false);
	tree->outgoingLines().insert(make_pair(gluon,sg));
	}
	else if(itype==1) {
	Lorentz5Momentum pin(pnew[0]),pout(pnew[1]),pgluon(pnew[2]);
	if(iemit==2) swap(pin,pout);
	// ensure outgoing quark can be put on-shell
	Lorentz5Momentum ptest(pout);
	if(ptest.boost(-(pin+pgluon).boostVector()).e() <
	incoming[1]->dataPtr()->constituentMass()) return;
	// create the new gluon
	PPtr newg = getParticleData(ParticleID::g)->produceParticle(pgluon);
	// create the new outgoing quark
	PPtr newout= getParticleData(-incoming[1]->id())->produceParticle(pout);
	// create the new incoming quark
	PPtr newin = new_ptr(Particle(*incoming[0]));
	newin->set5Momentum(pin);
	// colour info
	ColinePtr col=incoming[0]->colourLine();
	col->removeColoured(newin);
	ColinePtr newline=new_ptr(ColourLine());
	newline->addColoured(newout);
	newline->addColoured(newg);
	col->addAntiColoured(newg);
	col->addColoured(newin);
	// change the existing incoming partons
	PPtr orig;
	for(cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
	if(cit->first->progenitor()->id()==newin->id()) {
	// remove old particles from colour line
	col->removeColoured(cit->first->copy());
	col->removeColoured(cit->first->progenitor());
	// insert new particles
	cit->first->copy(newin);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newin,1,false)));
	sp->x(xnew.first);
	cit->first->progenitor(sp);
	tree->incomingLines()[cit->first]=sp;
	cit->first->perturbative(true);
	}
	else {
	// remove old particles from colour line
	col->removeAntiColoured(cit->first->copy());
	col->removeColoured(cit->first->progenitor());
	// insert new particles
	cit->first->copy(newg);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newg,1,false)));
	sp->x(xnew.second);
	cit->first->progenitor(sp);
	tree->incomingLines()[cit->first]=sp;
	cit->first->perturbative(false);
	orig=cit->first->original();
	}
	}
	// fix the momentum of the gauge boson
	Boost boostv=pboson.findBoostToCM();
	trans *=LorentzRotation(boostv);
	for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt= tree->outgoingLines().begin();
	cjt != tree->outgoingLines().end();++cjt ) {
	cjt->first->progenitor()->transform(trans);
	cjt->first->copy()->transform(trans);
	}
	tree->hardMatrixElementCorrection(true);
	// add the outgoing quark
	ShowerParticlePtr sout=new_ptr(ShowerParticle(*newout,1,true));
	ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(orig,newout,sout));
	out->perturbative(false);
	tree->outgoingLines().insert(make_pair(out,sout));
	}
	else if(itype==2) {
	Lorentz5Momentum pin(pnew[0]),pout(pnew[1]),pgluon(pnew[2]);
	if(iemit==2) swap(pin,pout);
	// ensure outgoing antiquark can be put on-shell
	Lorentz5Momentum ptest(pout);
	if(ptest.boost(-(pin+pgluon).boostVector()).e() <
	incoming[0]->dataPtr()->constituentMass()) return;
	// create the new gluon
	PPtr newg = getParticleData(ParticleID::g)->produceParticle(pgluon);
	// create the new outgoing antiquark
	PPtr newout= getParticleData(-incoming[0]->id())->produceParticle(pout);
	// create the new incoming antiquark
	PPtr newin = new_ptr(Particle(*incoming[1]));
	newin->set5Momentum(pin);
	// colour info
	ColinePtr col=incoming[0]->colourLine();
	col->removeAntiColoured(newin);
	ColinePtr newline=new_ptr(ColourLine());
	newline->addAntiColoured(newout);
	newline->addAntiColoured(newg);
	col->addColoured(newg);
	col->addAntiColoured(newin);
	// change the existing incoming partons
	PPtr orig;
	for(cit=tree->incomingLines().begin();cit!=tree->incomingLines().end();++cit) {
	if(cit->first->progenitor()->id()==newin->id()) {
	// remove old particles from colour line
	col->removeAntiColoured(cit->first->copy());
	col->removeColoured(cit->first->progenitor());
	// insert new particles
	cit->first->copy(newin);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newin,1,false)));
	sp->x(xnew.second);
	cit->first->progenitor(sp);
	tree->incomingLines()[cit->first]=sp;
	cit->first->perturbative(true);
	}
	else {
	// remove old particles from colour line
	col->removeColoured(cit->first->copy());
	col->removeColoured(cit->first->progenitor());
	// insert new particles
	cit->first->copy(newg);
	ShowerParticlePtr sp(new_ptr(ShowerParticle(*newg,1,false)));
	sp->x(xnew.first);
	cit->first->progenitor(sp);
	tree->incomingLines()[cit->first]=sp;
	cit->first->perturbative(false);
	orig=cit->first->original();
	}
	}
	// fix the momentum of the gauge boson
	Boost boostv=pboson.findBoostToCM();
	trans *=LorentzRotation(boostv);
	for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt= tree->outgoingLines().begin();
	cjt != tree->outgoingLines().end();++cjt ) {
	cjt->first->progenitor()->transform(trans);
	cjt->first->copy()->transform(trans);
	}
	tree->hardMatrixElementCorrection(true);
	// add the outgoing antiquark
	ShowerParticlePtr sout=new_ptr(ShowerParticle(*newout,1,true));
	ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(orig,newout,sout));
	out->perturbative(false);
	tree->outgoingLines().insert(make_pair(out,sout));
	}
	}

	bool DrellYanBase::applyHard(ShowerParticleVector quarks,
	vector<tcBeamPtr> beams, Lorentz5Momentum boson,
	unsigned int & iemit,unsigned int & itype,
	vector<Lorentz5Momentum> & pnew,
	LorentzRotation & trans,
	pair<double,double> & xout) {
	// check that quark along +z and qbar along -z
	bool quarkplus=quarks[0]->momentum().z()>quarks[1]->momentum().z();
	// calculate the limits on s
	Energy mb = sqrt(mb2_);
	Energy2 smin=mb2_;
	Energy2 s=
	(generator()->currentEvent()->incoming().first ->momentum()+
	generator()->currentEvent()->incoming().second->momentum()).m2();
	Energy2 smax(s);
	// calculate the rapidity of the boson
	double yB=0.5*log((boson.e()+boson.z())/
	(boson.e()-boson.z()));
	if(!quarkplus) yB=-yB;
	// if no phase-space return
	if(smax<smin) return false;
	// get the evolution scales (this needs improving)
	double kappa[2]={1.,1.};
	// get the momentum fractions for the leading order process
	// and the values of the PDF's
	double x[2],fx[2];
	tcPDFPtr pdf[2];
	for(unsigned int ix=0;ix<quarks.size();++ix) {
	x[ix]=quarks[ix]->x();
	assert(beams[ix]);
	pdf[ix]=beams[ix]->pdf();
	assert(pdf[ix]);
	fx[ix]=pdf[ix]->xfx(beams[ix],quarks[ix]->dataPtr(),mb2_,x[ix]);
	}
	// select the type of process and generate the kinematics
	double rn(UseRandom::rnd());
	Energy2 shat(ZERO),uhat(ZERO),that(ZERO);
	double weight(0.),xnew[2]={1.,1.};
	// generate the value of s according to 1/s^2
	shat = smaxsmin/(smin+UseRandom::rnd()(smax-smin));
	Energy2 jacobian = sqr(shat)*(1./smin-1./smax);
	double sbar=shat/mb2_;
	// calculate limits on that
	Energy2 tmax=mb2_kappa[0](1.-sbar)/(kappa[0]+sbar);
	Energy2 tmin=shat*(1.-sbar)/(kappa[1]+sbar);
	// calculate the limits on uhat
	Energy2 umax(mb2_-shat-tmin),umin(mb2_-shat-tmax);
	// check inside phase space
	if(tmax<tmin\|\|umax<umin) return false;
	// q qbar -> g V
	if(rn<_channelweights[0]) {
	// generate t and u according to 1/t+1/u
	// generate in 1/t
	if(UseRandom::rndbool(0.5)) {
	that=tmax*pow(tmin/tmax,UseRandom::rnd());
	uhat=mb2_-shat-that;
	jacobian *=log(tmin/tmax);
	}
	// generate in 1/u
	else {
	uhat=umax*pow(umin/umax,UseRandom::rnd());
	that=mb2_-shat-uhat;
	jacobian *=log(umin/umax);
	}
	Energy4 jacobian2 = jacobian * 2.uhatthat/(shat-mb2_);
	// new scale (this is mt^2=pt^2+mb^2)
	Energy2 scale(uhat*that/shat+mb2_);
	// the PDF's with the emitted gluon
	double fxnew[2];
	xnew[0]=exp(yB)/sqrt(s)sqrt(shat(mb2_-uhat)/(mb2_-that));
	xnew[1]=shat/(s*xnew[0]);
	for(unsigned int ix=0;ix<2;++ix)
	{fxnew[ix]=pdf[ix]->xfx(beams[ix],quarks[ix]->dataPtr(),scale,xnew[ix]);}
	// jacobian and me parts of the weight
	weight=jacobian2(sqr(mb2_-that)+sqr(mb2_-uhat))/(sqr(shat)that*uhat);
	// pdf part of the weight
	weight =fxnew[0]fxnew[1]x[0]x[1]/(fx[0]fx[1]xnew[0]*xnew[1]);
	// finally coupling, colour factor and different channel pieces
	weight = 2./3./Constants::pi/_channelweights[0]_alpha->value(scale);
	// select the emiting particle
	iemit=1;
	if(UseRandom::rnd()<sqr(mb2_-uhat)/(sqr(mb2_-uhat)+sqr(mb2_-that))) iemit=2;
	itype=0;
	}
	// incoming gluon
	else {
	// generate t
	if(rn>_channelweights[1]) {
	swap(tmax,tmin);
	tmax=mb2_-shat-tmax;
	tmin=mb2_-shat-tmin;
	}
	that=tmax*pow(tmin/tmax,UseRandom::rnd());
	uhat=mb2_-shat-that;
	Energy4 jacobian2 = jacobian * that*log(tmax/tmin);
	// new scale (this is mt^2=pt^2+mb^2)
	Energy2 scale(uhat*that/shat+mb2_);
	// g qbar -> qbar V
	double fxnew[2];
	if(rn<_channelweights[1]) {
	itype=2;
	xnew[0]=exp(yB)/sqrt(s)sqrt(shat(mb2_-uhat)/(mb2_-that));
	xnew[1]=shat/(s*xnew[0]);
	fxnew[0]=pdf[0]->xfx(beams[0],getParticleData(ParticleID::g),scale,xnew[0]);
	fxnew[1]=pdf[1]->xfx(beams[1],quarks[1]->dataPtr(),scale,xnew[1]);
	jacobian2/=(_channelweights[1]-_channelweights[0]);
	}
	// q g -> q V
	else {
	itype=1;
	xnew[0]=exp(yB)/sqrt(s)sqrt(shat(mb2_-that)/(mb2_-uhat));
	xnew[1]=shat/(s*xnew[0]);
	fxnew[0]=pdf[0]->xfx(beams[0],quarks[0]->dataPtr(),scale,xnew[0]);
	fxnew[1]=pdf[1]->xfx(beams[1],getParticleData(ParticleID::g),scale,xnew[1]);
	jacobian2/=(_channelweights[2]-_channelweights[1]);
	}
	// jacobian and me parts of the weight
	weight=-jacobian2(sqr(mb2_-that)+sqr(mb2_-shat))/(sqr(shat)shat*that);
	// pdf part of the weight
	weight =fxnew[0]fxnew[1]x[0]x[1]/(fx[0]fx[1]xnew[0]*xnew[1]);
	// finally coupling, colour factor and different channel pieces
	weight = 0.25/Constants::pi_alpha->value(scale);
	// select the emiting particle
	iemit=1;
	if(UseRandom::rnd()<sqr(mb2_-that)/(sqr(mb2_-that)+sqr(mb2_-shat))) iemit=2;
	}
	// if me correction should be applied
	if(weight>1.) {
	++_nover;
	_maxwgt=max(_maxwgt,weight);
	weight=1.;
	}
	if(UseRandom::rnd()>weight) return false;
	// construct the momenta in the rest frame of the boson
	Lorentz5Momentum pb(ZERO,ZERO,ZERO,mb,mb),pspect,pg,pemit;
	double cos3 = 0.0;
	if(itype==0)
	{
	pg = Lorentz5Momentum(ZERO,ZERO,ZERO,0.5*(shat-mb2_)/mb,ZERO);
	Energy2 tp(that),up(uhat);
	double zsign(-1.);
	if(iemit==2)
	{
	tp=uhat;
	up=that;
	zsign=1.;
	}
	pspect = Lorentz5Momentum(ZERO,ZERO,zsign0.5(mb2_-tp)/mb,
	0.5*(mb2_-tp)/mb,ZERO);
	Energy eemit=0.5*(mb2_-up)/mb;
	cos3 = 0.5/pspect.z()/pg.e()*(sqr(pspect.e())+sqr(pg.e())-sqr(eemit));
	}
	else
	{
	pg=Lorentz5Momentum(ZERO,ZERO,ZERO,0.5*(mb2_-uhat)/mb,ZERO);
	double zsign(1.);
	if(iemit==1)
	{
	if(itype==1) zsign=-1.;
	pspect=Lorentz5Momentum(ZERO,ZERO,0.5zsign(shat-mb2_)/mb,
	0.5*(shat-mb2_)/mb);
	Energy eemit=0.5*(mb2_-that)/mb;
	cos3 = 0.5/pspect.z()/pg.e()*(sqr(pspect.e())+sqr(pg.e())-sqr(eemit));
	}
	else
	{
	if(itype==2) zsign=-1.;
	pspect=Lorentz5Momentum(ZERO,ZERO,0.5zsign(mb2_-that)/mb,
	0.5*(mb2_-that)/mb);
	Energy eemit=0.5*(shat-mb2_)/mb;
	cos3 = 0.5/pspect.z()/pg.e()*(-sqr(pspect.e())-sqr(pg.e())+sqr(eemit));
	}
	}

	// rotate the gluon
	double sin3(sqrt(1.-sqr(cos3)));
	double phi(Constants::twopi*UseRandom::rnd());
	pg.setX(pg.e()sin3cos(phi));
	pg.setY(pg.e()sin3sin(phi));
	pg.setZ(pg.e()*cos3);
	if(itype==0) {
	pemit=pb+pg-pspect;
	}
	else {
	if(iemit==1) pemit=pb+pspect-pg;
	else pemit=pspect+pg-pb;
	}
	pemit.rescaleMass();
	// find the new CMF
	Lorentz5Momentum pp[2];
	if(itype==0) {
	if(iemit==1) {
	pp[0]=pemit;
	pp[1]=pspect;
	}
	else {
	pp[0]=pspect;
	pp[1]=pemit;
	}
	}
	else if(itype==1) {
	pp[1]=pg;
	if(iemit==1) pp[0]=pemit;
	else pp[0]=pspect;
	}
	else {
	pp[0]=pg;
	if(iemit==1) pp[1]=pemit;
	else pp[1]=pspect;
	}
	Lorentz5Momentum pz= quarkplus ? pp[0] : pp[1];
	pp[0]/=xnew[0];
	pp[1]/=xnew[1];
	Lorentz5Momentum plab(pp[0]+pp[1]);
	plab.rescaleMass();
	// construct the boost to rest frame of plab
	trans=LorentzRotation(plab.findBoostToCM());
	pz.transform(trans);
	// rotate so emitting particle along z axis
	// rotate so in x-z plane
	trans.rotateZ(-atan2(pz.y(),pz.x()));
	// rotate so along
	trans.rotateY(-acos(pz.z()/pz.vect().mag()));
	// undo azimuthal rotation
	trans.rotateZ(atan2(pz.y(),pz.x()));
	// perform the transforms
	pb .transform(trans);
	pspect.transform(trans);
	pg .transform(trans);
	pemit .transform(trans);
	// momenta to be returned
	pnew.push_back(pemit);
	pnew.push_back(pspect);
	pnew.push_back(pg);
	pnew.push_back(pb);
	xout.first=xnew[0];
	xout.second=xnew[1];
	return true;
	}

	bool DrellYanBase::softMatrixElementVeto(ShowerProgenitorPtr initial,
	ShowerParticlePtr parent,Branching br) {
	if(parent->isFinalState()) return false;
	// check if me correction should be applied
	long id[2]={initial->id(),parent->id()};
	if(id[0]!=id[1]\|\|id[1]==ParticleID::g) return false;
	// get the pT
	Energy pT=br.kinematics->pT();
	// check if hardest so far
	if(pT<initial->highestpT()) return false;
	// compute the invariants
	double kappa(sqr(br.kinematics->scale())/mb2_),z(br.kinematics->z());
	Energy2 shat(mb2_/z(1.+(1.-z)kappa)),that(-(1.-z)kappamb2_),uhat(-(1.-z)*shat);
	// check which type of process
	// g qbar
	double wgt(1.);
	if(id[0]>0&&br.ids[0]==ParticleID::g)
	wgt=mb2_/(shat+uhat)*(sqr(mb2_-that)+sqr(mb2_-shat))/(sqr(shat+uhat)+sqr(uhat));
	else if(id[0]>0&&br.ids[0]==id[0])
	wgt=mb2_/(shat+uhat)*(sqr(mb2_-uhat)+sqr(mb2_-that))/(sqr(shat)+sqr(shat+uhat));
	else if(id[0]<0&&br.ids[0]==ParticleID::g)
	wgt=mb2_/(shat+uhat)*(sqr(mb2_-that)+sqr(mb2_-shat))/(sqr(shat+uhat)+sqr(uhat));
	else if(id[0]<0&&br.ids[0]==-id[0])
	wgt=mb2_/(shat+uhat)*(sqr(mb2_-uhat)+sqr(mb2_-that))/(sqr(shat)+sqr(shat+uhat));
	else return false;
	if(wgt<.0\|\|wgt>1.) generator()->log() << "Soft ME correction weight too large or "
	<< "negative in DrellYanBase::"
	<< "softMatrixElementVeto()soft weight "
	<< " sbar = " << shat/mb2_
	<< " tbar = " << that/mb2_
	<< "weight = " << wgt << "\n";
	// if not vetoed
	if(UseRandom::rndbool(wgt)) {
	initial->highestpT(pT);
	return false;
	}
	// otherwise
	- parent->evolutionScale(br.kinematics->scale());
	+ parent->vetoEmission(parent->id()>0 ? ShowerPartnerType::QCDColourLine :
	+ ShowerPartnerType::QCDColourLine,br.kinematics->scale());
	return true;
	}

	HardTreePtr DrellYanBase::generateHardest(ShowerTreePtr tree,
	vector<ShowerInteraction::Type> inter) {
	bool found = false;
	// check if generating QCD radiation
	for(unsigned int ix=0;ix<inter.size();++ix) {
	found \|= inter[ix]==ShowerInteraction::QCD;
	}
	if(!found) return HardTreePtr();
	useMe();
	// get the particles to be showered
	_beams.clear();
	_partons.clear();
	// find the incoming particles
	ShowerParticleVector incoming;
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
	_quarkplus = true;
	vector<ShowerProgenitorPtr> particlesToShower;
	//progenitor particles are produced in z direction.
	for( cit = tree->incomingLines().begin(); cit != tree->incomingLines().end(); ++cit ) {
	incoming.push_back( cit->first->progenitor() );
	_beams.push_back( cit->first->beam() );
	_partons.push_back( cit->first->progenitor()->dataPtr() );
	// check that quark is along +ve z direction
	if(cit->first->progenitor()->id() > 0 &&
	cit->first->progenitor()->momentum().z() < ZERO )
	_quarkplus = false;
	particlesToShower.push_back( cit->first );
	}
	Lorentz5Momentum pboson;
	for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt= tree->outgoingLines().begin();
	cjt != tree->outgoingLines().end();++cjt ) {
	pboson += cjt->first->original()->momentum();
	}
	pboson.rescaleMass();
	// calculate the rapidity of the boson
	_yb = 0.5 * log((pboson.e()+pboson.z())/(pboson.e()-pboson.z()));
	_yb *= _quarkplus ? 1. : -1.;
	_mass = pboson.mass();
	// we are assuming quark first, swap order to ensure this
	// if antiquark first
	if(_partons[0]->id()<_partons[1]->id()) {
	swap(_partons[0],_partons[1]);
	swap(_beams[0],_beams[1]);
	}
	vector<Lorentz5Momentum> pnew;
	int emission_type(-1);
	// generate the hard emission and return if no emission
	if(!getEvent(pnew,emission_type)) {
	// for(unsigned int ix=0;ix<particlesToShower.size();++ix)
	// particlesToShower[ix]->maximumpT(_min_pt,ShowerInteraction::QCD);
	return HardTreePtr();
	}
	// construct the HardTree object needed to perform the showers
	ShowerParticleVector newparticles;
	// make the particles for the HardTree
	tcPDPtr gluon=getParticleData(ParticleID::g);
	// create the partons
	int iemit;
	// q qbar -> g V
	if(emission_type==0) {
	newparticles.push_back(new_ptr(ShowerParticle(_partons[0] ,false)));
	newparticles.push_back(new_ptr(ShowerParticle(_partons[1] ,false)));
	newparticles.push_back(new_ptr(ShowerParticle(gluon , true)));
	iemit = pnew[0].z()/pnew[2].rapidity()>ZERO ? 0 : 1;
	}
	// q g -> q V
	else if(emission_type==1) {
	iemit=1;
	newparticles.push_back(new_ptr(ShowerParticle(_partons[0] ,false)));
	newparticles.push_back(new_ptr(ShowerParticle(gluon ,false)));
	newparticles.push_back(new_ptr(ShowerParticle(_partons[1]->CC(), true)));
	}
	// g qbar -> qbar V
	else {
	iemit=0;
	newparticles.push_back(new_ptr(ShowerParticle(gluon ,false)));
	newparticles.push_back(new_ptr(ShowerParticle(_partons[1] ,false)));
	newparticles.push_back(new_ptr(ShowerParticle(_partons[0]->CC(), true)));
	}
	// set the momenta
	for(unsigned int ix=0;ix<3;++ix) newparticles[ix]->set5Momentum(pnew[ix]);
	// create the off-shell particle
	Lorentz5Momentum poff=pnew[iemit]-pnew[2];
	poff.rescaleMass();
	newparticles.push_back(new_ptr(ShowerParticle(_partons[iemit],false)));
	newparticles.back()->set5Momentum(poff);
	// compute the boost for the bosons
	LorentzRotation boost(pboson.findBoostToCM());
	boost.boost(pnew[3].boostVector());
	for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt= tree->outgoingLines().begin();
	cjt != tree->outgoingLines().end();++cjt ) {
	newparticles.push_back(new_ptr(ShowerParticle(cjt->first->original()->dataPtr(),
	true)));
	newparticles.back()->set5Momentum(boost*cjt->first->original()->momentum());
	}
	vector<HardBranchingPtr> inBranch,hardBranch;
	// create the branchings for the incoming particles
	inBranch.push_back(new_ptr(HardBranching(newparticles[0],SudakovPtr(),
	HardBranchingPtr(),HardBranching::Incoming)));
	inBranch.push_back(new_ptr(HardBranching(newparticles[1],SudakovPtr(),
	HardBranchingPtr(),HardBranching::Incoming)));
	// intermediate IS particle
	hardBranch.push_back(new_ptr(HardBranching(newparticles[3],SudakovPtr(),
	inBranch[iemit],HardBranching::Incoming)));
	inBranch[iemit]->addChild(hardBranch.back());
	// create the branching for the emitted jet
	inBranch[iemit]->addChild(new_ptr(HardBranching(newparticles[2],SudakovPtr(),
	inBranch[iemit],HardBranching::Outgoing)));
	// set the colour partners
	hardBranch.back()->colourPartner(inBranch[iemit==0 ? 1 : 0]);
	inBranch[iemit==0 ? 1 : 0]->colourPartner(hardBranch.back());
	// add other particle
	hardBranch.push_back(inBranch[iemit==0 ? 1 : 0]);
	// outgoing particles
	for(unsigned int ix=4;ix<newparticles.size();++ix) {
	hardBranch.push_back(new_ptr(HardBranching(newparticles[ix],SudakovPtr(),
	HardBranchingPtr(),HardBranching::Outgoing)));
	}
	// make the tree
	HardTreePtr hardtree=new_ptr(HardTree(hardBranch,inBranch,ShowerInteraction::QCD));
	// connect the ShowerParticles with the branchings
	// and set the maximum pt for the radiation
	set<HardBranchingPtr> hard=hardtree->branchings();
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if( _pt < _min_pt ) particlesToShower[ix]->maximumpT(_min_pt,ShowerInteraction::QCD);
	else particlesToShower[ix]->maximumpT(_pt,ShowerInteraction::QCD);
	for(set<HardBranchingPtr>::const_iterator mit=hard.begin();
	mit!=hard.end();++mit) {
	if(particlesToShower[ix]->progenitor()->id()==(*mit)->branchingParticle()->id()&&
	(( particlesToShower[ix]->progenitor()->isFinalState()&&
	(**mit).status()==HardBranching::Outgoing)\|\|
	(!particlesToShower[ix]->progenitor()->isFinalState()&&
	(**mit).status()==HardBranching::Incoming))) {
	hardtree->connect(particlesToShower[ix]->progenitor(),*mit);
	if((**mit).status()==HardBranching::Incoming) {
	(*mit)->beam(particlesToShower[ix]->original()->parents()[0]);
	}
	HardBranchingPtr parent=(*mit)->parent();
	while(parent) {
	parent->beam(particlesToShower[ix]->original()->parents()[0]);
	parent=parent->parent();
	};
	}
	}
	}
	ColinePtr newline=new_ptr(ColourLine());
	for(set<HardBranchingPtr>::const_iterator cit=hardtree->branchings().begin();
	cit!=hardtree->branchings().end();++cit) {
	if((**cit).branchingParticle()->dataPtr()->iColour()==PDT::Colour3)
	newline->addColoured((**cit).branchingParticle());
	else if((**cit).branchingParticle()->dataPtr()->iColour()==PDT::Colour3bar)
	newline->addAntiColoured((**cit).branchingParticle());
	}
	// return the tree
	return hardtree;
	}

	double DrellYanBase::getResult(int emis_type, Energy pt, double yj) {
	Energy2 s=sqr(generator()->maximumCMEnergy());
	Energy2 m2(sqr(_mass));
	Energy2 scale = m2+sqr(pt);
	Energy et=sqrt(scale);
	// longitudinal real correction fractions
	double x = ptexp( yj)/sqrt(s)+etexp( _yb)/sqrt(s);
	double y = ptexp(-yj)/sqrt(s)+etexp(-_yb)/sqrt(s);
	// reject if outside region
	if(x<0.\|\|x>1.\|\|y<0.\|\|y>1.\|\|x*y<m2/s) return 0.;
	// longitudinal born fractions
	double x1 = _mass*exp( _yb)/sqrt(s);
	double y1 = _mass*exp(-_yb)/sqrt(s);
	// mandelstam variables
	Energy2 th = -sqrt(s)xpt*exp(-yj);
	Energy2 uh = -sqrt(s)ypt*exp( yj);
	Energy2 sh = m2-th-uh;
	double res;
	// pdf part of the cross section
	double pdf[4];
	pdf[0]=_beams[0]->pdf()->xfx(_beams[0],_partons[0],m2,x1);
	pdf[1]=_beams[1]->pdf()->xfx(_beams[1],_partons[1],m2,y1);
	//qqbar2Zg
	using Constants::pi;
	if(emis_type==0) {
	pdf[2]=_beams[0]->pdf()->xfx(_beams[0],_partons[0],scale,x);
	pdf[3]=_beams[1]->pdf()->xfx(_beams[1],_partons[1],scale,y);
	res = 4./3./pi(sqr(th-m2)+sqr(uh-m2))pt/(shuhth)*GeV;
	}
	//qg2Zq
	else if(emis_type==1) {
	pdf[2]=_beams[0]->pdf()->xfx(_beams[0],_partons[0],scale,x);
	pdf[3]=_beams[1]->pdf()->xfx(_beams[1],getParticleData(ParticleID::g),scale,y);
	res = -1./2./pi(sqr(uh-m2)+sqr(sh-m2))pt/(shshuh)*GeV;
	}
	//qbarg2Zqbar
	else {
	pdf[2]=_beams[0]->pdf()->xfx(_beams[0],getParticleData(ParticleID::g),scale,x);
	pdf[3]=_beams[1]->pdf()->xfx(_beams[1],_partons[1],scale,y);
	res =- 1./2./pi(sqr(th-m2)+sqr(sh-m2))pt/(shshth)*GeV;
	}
	//deals with pdf zero issue at large x
	if(pdf[0]<=0.\|\|pdf[1]<=0.\|\|pdf[2]<=0.\|\|pdf[3]<=0.) {
	res=0.;
	}
	else {
	res=pdf[2]pdf[3]/pdf[0]/pdf[1]*m2/sh;
	}
	res*=_alpha->ratio(sqr(pt));
	return res;
	}

	bool DrellYanBase::getEvent(vector<Lorentz5Momentum> & pnew,
	int & emis_type){
	// pt cut-off
	// Energy minp = 0.1*GeV;
	// maximum pt (half of centre-of-mass energy)
	Energy maxp = 0.5*generator()->maximumCMEnergy();
	// set pt of emission to zero
	_pt=ZERO;
	//Working Variables
	Energy pt;
	double yj;
	// limits on the rapidity of the jet
	double minyj = -8.0,maxyj = 8.0;
	bool reject;
	double wgt;
	emis_type=-1;
	for(int j=0;j<3;j++) {
	pt=maxp;
	double a = _alpha->overestimateValue()_prefactor[j](maxyj-minyj)/(_power-1.);
	do {
	// generate next pt
	pt=GeV/pow(pow(GeV/pt,_power-1)-log(UseRandom::rnd())/a,1./(_power-1.));
	// generate rapidity of the jet
	yj=UseRandom::rnd()*(maxyj-minyj)+ minyj;
	// calculate rejection weight
	wgt=getResult(j,pt,yj);
	wgt/= _prefactor[j]*pow(GeV/pt,_power);
	reject = UseRandom::rnd()>wgt;
	//no emission event if p goes past p min - basically set to outside
	//of the histogram bounds (hopefully hist object just ignores it)
	if(pt<_min_pt){
	pt=ZERO;
	reject = false;
	}
	if(wgt>1.0) {
	ostringstream s;
	s << "DrellYanBase::getEvent weight for channel " << j
	<< "is " << wgt << " which is greater than 1";
	generator()->logWarning( Exception(s.str(), Exception::warning) );
	}
	}
	while(reject);
	// set pt of emission etc
	if(pt>_pt){
	emis_type = j;
	_pt=pt;
	_yj=yj;
	}
	}
	//was this an (overall) no emission event?
	if(_pt<_min_pt){
	_pt=ZERO;
	emis_type = 3;
	}
	if(emis_type==3) return false;
	// generate the momenta of the particles
	// hadron-hadron cmf
	Energy2 s=sqr(generator()->maximumCMEnergy());
	// transverse energy
	Energy2 m2(sqr(_mass));
	Energy et=sqrt(m2+sqr(_pt));
	// first calculate all the kinematic variables
	// longitudinal real correction fractions
	double x = _ptexp( _yj)/sqrt(s)+etexp( _yb)/sqrt(s);
	double y = _ptexp(-_yj)/sqrt(s)+etexp(-_yb)/sqrt(s);
	// that and uhat
	Energy2 th = -sqrt(s)x_pt*exp(-_yj);
	Energy2 uh = -sqrt(s)y_pt*exp( _yj);
	Energy2 sh = xys;
	if(emis_type==1) swap(th,uh);
	// decide which was the emitting particle
	unsigned int iemit=1;
	// from q qbar
	if(emis_type==0) {
	if(UseRandom::rnd()<sqr(m2-uh)/(sqr(m2-uh)+sqr(m2-th))) iemit=2;
	}
	else {
	if(UseRandom::rnd()<sqr(m2-th)/(sqr(m2-th)+sqr(m2-sh))) iemit=2;
	}
	// reconstruct the momenta in the rest frame of the gauge boson
	Lorentz5Momentum pb(ZERO,ZERO,ZERO,_mass,_mass),pspect,pg,pemit;
	double cos3;
	if(emis_type==0) {
	pg=Lorentz5Momentum(ZERO,ZERO,ZERO,0.5*(sh-m2)/_mass,ZERO);
	Energy2 tp(th),up(uh);
	double zsign(-1.);
	if(iemit==2) {
	swap(tp,up);
	zsign=1;
	}
	pspect = Lorentz5Momentum(ZERO,ZERO
	,zsign0.5(m2-tp)/_mass,0.5*(m2-tp)/_mass,
	ZERO);
	Energy eemit=0.5*(m2-up)/_mass;
	cos3 = 0.5/pspect.z()/pg.e()*(sqr(pspect.e())+sqr(pg.e())-sqr(eemit));
	}
	else {
	pg=Lorentz5Momentum(ZERO,ZERO,ZERO,0.5*(m2-uh)/_mass,ZERO);
	double zsign(1.);
	if(iemit==1) {
	if(emis_type==1) zsign=-1.;
	pspect=Lorentz5Momentum(ZERO,ZERO,0.5zsign(sh-m2)/_mass,0.5*(sh-m2)/_mass);
	Energy eemit=0.5*(m2-th)/_mass;
	cos3 = 0.5/pspect.z()/pg.e()*(sqr(pspect.e())+sqr(pg.e())-sqr(eemit));
	}
	else {
	if(emis_type==2) zsign=-1.;
	pspect=Lorentz5Momentum(ZERO,ZERO,0.5zsign(m2-th)/_mass,0.5*(m2-th)/_mass);
	Energy eemit=0.5*(sh-m2)/_mass;
	cos3 = 0.5/pspect.z()/pg.e()*(-sqr(pspect.e())-sqr(pg.e())+sqr(eemit));
	}
	}
	// rotate the gluon
	double sin3(sqrt(1.-sqr(cos3))),phi(Constants::twopi*UseRandom::rnd());
	pg.setX(pg.e()sin3cos(phi));
	pg.setY(pg.e()sin3sin(phi));
	pg.setZ(pg.e()*cos3);
	if(emis_type==0) {
	pemit=pb+pg-pspect;
	}
	else {
	if(iemit==1) pemit=pb+pspect-pg;
	else pemit=pspect+pg-pb;
	}
	pemit .setMass(ZERO);
	pg .setMass(ZERO);
	pspect.setMass(ZERO);
	// find the new CMF
	Lorentz5Momentum pp[2];
	if(emis_type==0) {
	pp[0]=pemit;
	pp[1]=pspect;
	if(iemit==2) swap(pp[0],pp[1]);
	}
	else if(emis_type==1) {
	pp[1]=pg;
	if(iemit==1) pp[0]=pemit;
	else pp[0]=pspect;
	}
	else {
	pp[0]=pg;
	if(iemit==1) pp[1]=pemit;
	else pp[1]=pspect;
	}
	Lorentz5Momentum pz= _quarkplus ? pp[0] : pp[1];
	pp[0]/=x;
	pp[1]/=y;
	Lorentz5Momentum plab(pp[0]+pp[1]);
	plab.rescaleMass();
	// construct the boost to rest frame of plab
	LorentzRotation trans=LorentzRotation(plab.findBoostToCM());
	pz.transform(trans);
	// rotate so emitting particle along z axis
	// rotate so in x-z plane
	trans.rotateZ(-atan2(pz.y(),pz.x()));
	// rotate so along
	trans.rotateY(-acos(pz.z()/pz.vect().mag()));
	// undo azimuthal rotation
	trans.rotateZ(atan2(pz.y(),pz.x()));
	// perform the transforms
	pb .transform(trans);
	pspect.transform(trans);
	pg .transform(trans);
	pemit .transform(trans);
	// copy the momenta for the new particles
	pnew.resize(4);
	if(emis_type==0) {
	pnew[0]=pemit;
	pnew[1]=pspect;
	if(iemit==2) swap(pnew[0],pnew[1]);
	pnew[2]=pg;
	}
	else if(emis_type==1) {
	pnew[0]=pemit;
	pnew[2]=pspect;
	if(iemit==2) swap(pnew[0],pnew[2]);
	pnew[1]=pg;
	}
	else if(emis_type==2) {
	pnew[1]=pspect;
	pnew[2]=pemit;
	if(iemit==1) swap(pnew[1],pnew[2]);
	pnew[0]=pg;
	}
	pnew[3]=pb;
	return true;
	}
	diff --git a/Shower/Base/Evolver.cc b/Shower/Base/Evolver.cc
	--- a/Shower/Base/Evolver.cc
	+++ b/Shower/Base/Evolver.cc
	@@ -1,1990 +1,1990 @@
	// -- 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->evolutionScale(bb.kinematics->scale());
	+ 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,
	Energy maxscale,
	Energy minmass,ShowerInteraction::Type type) {
	Branching fb;
	while (true) {
	fb=_splittingGenerator->chooseDecayBranching(*particle,maxscale,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->evolutionScale(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],maxscale,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(Energy maxscale,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(),maxscale,
	minimumMass, branch ,type);
	}
	}
	return hardOnly() ? false :
	spaceLikeDecayShower(progenitor()->progenitor(),maxscale,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)
	_hardtree = _hardme->generateHardest( currentTree(),interactions_ );
	else
	_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()));
	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->evolutionScale(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->evolutionScale(bb.kinematics->scale() );
	continue;
	}
	// and those from the base class
	if(spaceLikeVetoed(bb,particle,type)) {
	particle->evolutionScale(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, bb.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, Energy maxscale,
	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,maxscale,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->evolutionScale(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->evolutionScale(fb.kinematics->scale());
	continue;
	}
	}
	// pt veto
	if(fb.kinematics->pT() > progenitor()->maximumpT(type)) {
	particle->evolutionScale(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->evolutionScale(fb.kinematics->scale());
	}
	// if no branching set decay matrix and return
	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->scale() );
	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],maxscale,minmass,type);
	}
	else {
	truncatedSpaceLikeDecayShower( theChildren[0],maxscale,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],maxscale,minmass,type);
	}
	else {
	truncatedSpaceLikeDecayShower( theChildren[1],maxscale,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],maxscale,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);
	}
	}
	}

	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()) {
	if(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 {
	// perform shower
	// set the scales correctly. The current scale is the maximum scale for
	// emission not the starting scale
	Energy maxscale=progenitor()->progenitor()->evolutionScale();
	Energy startScale=progenitor()->progenitor()->mass();
	progenitor()->progenitor()->evolutionScale(startScale);
	// perform the shower
	progenitor()->hasEmitted(startSpaceLikeDecayShower(maxscale,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/ShowerParticle.cc b/Shower/Base/ShowerParticle.cc
	--- a/Shower/Base/ShowerParticle.cc
	+++ b/Shower/Base/ShowerParticle.cc
	@@ -1,42 +1,43 @@
	// -- C++ --
	//
	// ShowerParticle.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 ShowerParticle class.
	//

	#include "ShowerParticle.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace Herwig;

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

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

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

	void ShowerParticle::vetoEmission(ShowerPartnerType::Type type, Energy scale) {
	for(map<ShowerPartnerType::Type ,pair<Energy,Energy> >::iterator
	it=evolutionScales_.begin();it!=evolutionScales_.end();++it) {
	- if(it->first==type \|\| (it->first!=ShowerPartnerType::QED &&
	- type !=ShowerPartnerType::QED)) {
	+ if(it->first==type) {
	it->second = make_pair(scale,scale);
	}
	- else
	- assert(false);
	+ else {
	+ if(it->second.first >scale) it->second.first = scale;
	+ if(it->second.second>scale) it->second.second = scale;
	+ }
	}
	}
	diff --git a/Shower/Default/FS_QTildeShowerKinematics1to2.cc b/Shower/Default/FS_QTildeShowerKinematics1to2.cc
	--- a/Shower/Default/FS_QTildeShowerKinematics1to2.cc
	+++ b/Shower/Default/FS_QTildeShowerKinematics1to2.cc
	@@ -1,311 +1,314 @@
	// -- C++ --
	//
	// FS_QTildeShowerKinematics1to2.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 FS_QTildeShowerKinematics1to2 class.
	//

	#include "FS_QTildeShowerKinematics1to2.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "Herwig++/Shower/SplittingFunctions/SplittingFunction.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	+#include "ThePEG/Utilities/Debug.h"

	using namespace Herwig;

	void FS_QTildeShowerKinematics1to2::
	updateChildren(const tShowerParticlePtr parent,
	const ShowerParticleVector & children,
	ShowerPartnerType::Type partnerType,
	bool angularOrder) const {
	assert(children.size()==2);
	// resize the parameter vectors
	if(parent->showerVariables().empty()) {
	parent->showerVariables().resize(3);
	parent->showerParameters().resize(2);
	parent->showerParameters()[0]=1.;
	}
	children[0]->showerVariables() .resize(3);
	children[0]->showerParameters().resize(2);
	children[1]->showerVariables() .resize(3);
	children[1]->showerParameters().resize(2);
	// identify emittor and spectator
	tShowerParticlePtr emittor = children[0];
	tShowerParticlePtr emitted = children[1];
	double zChild[2] = {z(),1.-z()};
	bool bosonSplitting(false);
	// special for g -> gg, particle highest z is emittor
	if(emittor->id()==emitted->id()&&emittor->id()==parent->id()) {
	if(zChild[1]>zChild[0]) {
	swap(zChild[0],zChild[1]);
	swap(emitted,emittor);
	}
	}
	// otherwise if particle ID same
	else if(emitted->id()==parent->id()) {
	swap(zChild[0],zChild[1]);
	swap(emitted,emittor);
	}
	// no real emittor/eemitted
	else if(emittor->id()!=parent->id()) {
	bosonSplitting = true;
	}
	// scales for angular ordering
	Energy AOScale[2];
	if(angularOrder) {
	for(unsigned int ix=0;ix<2;++ix) AOScale [ix] = zChild[ix]*scale();
	}
	else {
	for(unsigned int ix=0;ix<2;++ix) AOScale [ix] = scale();
	}
	// if not these cases doesn't matter
	// now the various scales
	if(partnerType==ShowerPartnerType::QED) {
	// normal case
	if(!bosonSplitting) {
	for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	it = parent->evolutionScales().begin();
	it!=parent->evolutionScales().end();++it) {
	emittor->evolutionScale(it->first,make_pair(min(AOScale[0],it->second.first ),
	min(scale() ,it->second.second)));
	if(it->first==ShowerPartnerType::QED) {
	emitted->evolutionScale(it->first,make_pair(min(AOScale[1],it->second.first ),
	min(scale() ,it->second.second)));
	}
	}
	}
	// gamma -> f fbar
	else {
	// set QED scales as normal
	for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	it = parent->evolutionScales().begin();
	it!=parent->evolutionScales().end();++it) {
	assert(it->first==ShowerPartnerType::QED);
	emittor->evolutionScale(it->first,make_pair(AOScale[0],scale()));
	emitted->evolutionScale(it->first,make_pair(AOScale[1],scale()));
	}
	// and any QCD scales needed
	PDT::Colour emittorColour = emittor->dataPtr()->iColour();
	if(emittorColour==PDT::Colour3\|\|emittorColour==PDT::Colour8\|\|emittorColour==PDT::Colour6) {
	emittor->evolutionScale(ShowerPartnerType::QCDColourLine,
	make_pair(AOScale[0],scale()));
	}
	if(emittorColour==PDT::Colour3bar\|\|emittorColour==PDT::Colour8\|\|emittorColour==PDT::Colour6bar) {
	emittor->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
	make_pair(AOScale[0],scale()));
	}
	PDT::Colour emittedColour = emitted->dataPtr()->iColour();
	if(emittedColour==PDT::Colour3\|\|emittedColour==PDT::Colour8\|\|emittedColour==PDT::Colour6) {
	emitted->evolutionScale(ShowerPartnerType::QCDColourLine,
	make_pair(AOScale[1],scale()));
	}
	if(emittedColour==PDT::Colour3bar\|\|emittedColour==PDT::Colour8\|\|emittedColour==PDT::Colour6bar) {
	emitted->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
	make_pair(AOScale[1],scale()));
	}
	}
	}
	else {
	// normal case
	if(!bosonSplitting) {
	// scales for the emittor
	for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	it = parent->evolutionScales().begin();
	it!=parent->evolutionScales().end();++it) {
	emittor->evolutionScale(it->first,make_pair(min(AOScale[0],it->second.first ),
	min(scale() ,it->second.second)));
	}
	}
	else {
	// QCD scales needed
	PDT::Colour emittorColour = emittor->dataPtr()->iColour();
	if(emittorColour==PDT::Colour3\|\|emittorColour==PDT::Colour8\|\|emittorColour==PDT::Colour6) {
	emittor->evolutionScale(ShowerPartnerType::QCDColourLine,
	make_pair(AOScale[0],scale()));
	}
	if(emittorColour==PDT::Colour3bar\|\|emittorColour==PDT::Colour8\|\|emittorColour==PDT::Colour6bar) {
	emittor->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
	make_pair(AOScale[0],scale()));
	}
	// QED scales
	if(emittor->dataPtr()->charged()) {
	emittor->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale[0],scale()));
	}
	}
	PDT::Colour emittedColour = emitted->dataPtr()->iColour();
	if(emittedColour==PDT::Colour3\|\|emittedColour==PDT::Colour8\|\|emittedColour==PDT::Colour6) {
	emitted->evolutionScale(ShowerPartnerType::QCDColourLine,
	make_pair(AOScale[1],scale()));
	}
	if(emittedColour==PDT::Colour3bar\|\|emittedColour==PDT::Colour8\|\|emittedColour==PDT::Colour6bar) {
	emitted->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
	make_pair(AOScale[1],scale()));
	}
	if(emitted->dataPtr()->charged()) {
	emitted->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale[1],scale()));
	}
	}
	+ // debugging printout ifd needed
	+ if(Debug::level >= 10 ) printScales(parent,children[0],children[1]);
	// determine alphas of children according to interpretation of z
	children[0]->showerParameters()[0]= z() *parent->showerParameters()[0];
	children[1]->showerParameters()[0]= (1.-z())*parent->showerParameters()[0];
	// set the values
	children[0]->showerVariables()[0]=
	pT()cos(phi()) + z() parent->showerVariables()[0];
	children[0]->showerVariables()[1]=
	pT()sin(phi()) + z() parent->showerVariables()[1];
	children[1]->showerVariables()[0]=
	- pT()cos(phi()) + (1.-z())parent->showerVariables()[0];
	children[1]->showerVariables()[1]=
	- pT()sin(phi()) + (1.-z())parent->showerVariables()[1];
	for(unsigned int ix=0;ix<2;++ix)
	children[ix]->showerVariables()[2]=
	sqrt(sqr(children[ix]->showerVariables()[0])+
	sqr(children[ix]->showerVariables()[1]));
	// set up the colour connections
	splittingFn()->colourConnection(parent,children[0],children[1],partnerType,false);
	// make the products children of the parent
	parent->addChild(children[0]);
	parent->addChild(children[1]);
	}

	void FS_QTildeShowerKinematics1to2::
	reconstructParent(const tShowerParticlePtr parent,
	const ParticleVector & children ) const {
	if(children.size() != 2)
	throw Exception() << "FS_QTildeShowerKinematics1to2::updateParent() "
	<< "Warning! too many children!"
	<< Exception::eventerror;
	ShowerParticlePtr c1 = dynamic_ptr_cast<ShowerParticlePtr>(children[0]);
	ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(children[1]);
	parent->showerParameters()[1]=
	c1->showerParameters()[1] + c2->showerParameters()[1];
	parent->set5Momentum( c1->momentum() + c2->momentum() );
	}

	void FS_QTildeShowerKinematics1to2::reconstructLast(const tShowerParticlePtr theLast,
	unsigned int iopt,
	Energy mass) const {
	// set beta component and consequently all missing data from that,
	// using the nominal (i.e. PDT) mass.
	Energy theMass = mass > ZERO ? mass : theLast->data().constituentMass();
	theLast->showerParameters()[1]=
	(sqr(theMass) + sqr(theLast->showerVariables()[2])
	- sqr( theLast->showerParameters()[0] )*pVector().m2())
	/ ( 2.theLast->showerParameters()[0]p_dot_n() );
	// set that new momentum
	theLast->set5Momentum(sudakov2Momentum( theLast->showerParameters()[0],
	theLast->showerParameters()[1],
	theLast->showerVariables()[0],
	theLast->showerVariables()[1],iopt));
	}

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

	void FS_QTildeShowerKinematics1to2::updateParent(const tShowerParticlePtr parent,
	const ShowerParticleVector & children,
	ShowerPartnerType::Type,
	bool) const {
	IdList ids(3);
	ids[0] = parent->id();
	ids[1] = children[0]->id();
	ids[2] = children[1]->id();
	vector<Energy> virtualMasses = SudakovFormFactor()->virtualMasses(ids);
	if(children[0]->children().empty()) children[0]->virtualMass(virtualMasses[1]);
	if(children[1]->children().empty()) children[1]->virtualMass(virtualMasses[2]);
	// compute the new pT of the branching
	Energy2 pt2=sqr(z()(1.-z()))sqr(scale())
	- sqr(children[0]->virtualMass())*(1.-z())
	- sqr(children[1]->virtualMass())* z() ;
	if(ids[0]!=ParticleID::g) pt2 += z()(1.-z())sqr(virtualMasses[0]);
	Energy2 q2 =
	sqr(children[0]->virtualMass())/z() +
	sqr(children[1]->virtualMass())/(1.-z()) +
	pt2/z()/(1.-z());
	if(pt2<ZERO) {
	parent->virtualMass(ZERO);
	}
	else {
	parent->virtualMass(sqrt(q2));
	pT(sqrt(pt2));
	}
	}

	void FS_QTildeShowerKinematics1to2::
	resetChildren(const tShowerParticlePtr parent,
	const ShowerParticleVector & children) const {
	// set the values
	children[0]->showerVariables()[0]=
	pT()cos(phi()) + z() parent->showerVariables()[0];
	children[0]->showerVariables()[1]=
	pT()sin(phi()) + z() parent->showerVariables()[1];
	children[1]->showerVariables()[0]=
	- pT()cos(phi()) + (1.-z())parent->showerVariables()[0];
	children[1]->showerVariables()[1]=
	- pT()sin(phi()) + (1.-z())parent->showerVariables()[1];
	for(unsigned int ix=0;ix<2;++ix)
	children[ix]->showerVariables()[2]=
	sqrt(sqr(children[ix]->showerVariables()[0])+
	sqr(children[ix]->showerVariables()[1]));
	for(unsigned int ix=0;ix<children.size();++ix) {
	if(children[ix]->children().empty()) continue;
	ShowerParticleVector newChildren;
	for(unsigned int iy=0;iy<children[ix]->children().size();++iy)
	newChildren.push_back(dynamic_ptr_cast<ShowerParticlePtr>
	(children[ix]->children()[iy]));
	children[ix]->showerKinematics()->resetChildren(children[ix],newChildren);
	}
	}
	diff --git a/Shower/Default/IS_QTildeShowerKinematics1to2.cc b/Shower/Default/IS_QTildeShowerKinematics1to2.cc
	--- a/Shower/Default/IS_QTildeShowerKinematics1to2.cc
	+++ b/Shower/Default/IS_QTildeShowerKinematics1to2.cc
	@@ -1,168 +1,264 @@
	// -- C++ --
	//
	// IS_QTildeShowerKinematics1to2.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 IS_QTildeShowerKinematics1to2 class.
	//

	#include "IS_QTildeShowerKinematics1to2.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	+#include "ThePEG/Utilities/Debug.h"
	#include <cassert>

	using namespace Herwig;

	void IS_QTildeShowerKinematics1to2::
	updateChildren( const tShowerParticlePtr theParent,
	const ShowerParticleVector & theChildren,
	ShowerPartnerType::Type,
	bool ) const {
	theChildren[1]->showerVariables().resize(3);
	theChildren[1]->showerParameters().resize(2);
	theChildren[1]->showerParameters()[0]=
	(1.-z())*theParent->showerParameters()[0];
	double cphi = cos(phi()),sphi = sin(phi());
	theChildren[1]->showerVariables()[0] =
	(1.-z())theParent->showerVariables()[0] - cphipT();
	theChildren[1]->showerVariables()[1]=
	(1.-z())theParent->showerVariables()[1] - sphipT();
	theChildren[1]->showerVariables()[2]=
	sqrt(sqr(theChildren[1]->showerVariables()[0])+
	sqr(theChildren[1]->showerVariables()[1]));
	// space-like child
	theChildren[0]->showerParameters()[0] =
	theParent->showerParameters()[0] - theChildren[1]->showerParameters()[0];
	theChildren[0]->showerParameters()[1] =
	theParent->showerParameters()[1] - theChildren[1]->showerParameters()[1];
	theChildren[0]->showerVariables()[0]=
	theParent->showerVariables()[0] - theChildren[1]->showerVariables()[0];
	theChildren[0]->showerVariables()[1]=
	theParent->showerVariables()[1] - theChildren[1]->showerVariables()[1];
	}


	void IS_QTildeShowerKinematics1to2::
	-updateParent(const tShowerParticlePtr theParent,
	- const ShowerParticleVector & theChildren,
	+updateParent(const tShowerParticlePtr parent,
	+ const ShowerParticleVector & children,
	ShowerPartnerType::Type partnerType,
	bool angularOrder) const {
	- // no z for angular ordering in backward branchings
	- theParent->evolutionScale(scale());
	- if(angularOrder)
	- theChildren[1]->evolutionScale((1.-z())*scale());
	- else
	- theChildren[1]->evolutionScale(scale());
	+ // 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,0);
	+ 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,0);
	+ 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,0);
	+ 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 ,0);
	+ 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,0);
	+ 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,0);
	+ 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 ,0);
	+ 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)));
	+ }
	+ }
	+ // debugging printout if needed
	+ if(Debug::level >= 10 ) printScales(parent,children[0],children[1]);
	// set proper colour connections
	- splittingFn()->colourConnection(theParent,theChildren[0],theChildren[1],
	+ splittingFn()->colourConnection(parent,children[0],children[1],
	partnerType,true);
	// set proper parent/child relationships
	- theParent->addChild(theChildren[0]);
	- theParent->addChild(theChildren[1]);
	- theParent->x(theChildren[0]->x()/z());
	+ parent->addChild(children[0]);
	+ parent->addChild(children[1]);
	+ parent->x(children[0]->x()/z());
	// create the storage for the shower variables
	- theParent->showerVariables().resize(3);
	- theParent->showerParameters().resize(2);
	- if(theChildren[0]->showerVariables().empty()) {
	- theChildren[0]->showerVariables().resize(3);
	- theChildren[0]->showerParameters().resize(2);
	+ parent->showerVariables().resize(3);
	+ parent->showerParameters().resize(2);
	+ if(children[0]->showerVariables().empty()) {
	+ children[0]->showerVariables().resize(3);
	+ children[0]->showerParameters().resize(2);
	}
	}

	void IS_QTildeShowerKinematics1to2::
	reconstructParent(const tShowerParticlePtr theParent,
	const ParticleVector & theChildren ) const {
	ShowerParticlePtr c1 = dynamic_ptr_cast<ShowerParticlePtr>(theChildren[0]);
	ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(theChildren[1]);
	// get shower variables from 1st child in order to keep notation
	// parent->(c1, c2) clean even though the splitting was initiated
	// from c1. The name updateParent is still referring to the
	// timelike branching though.
	// on-shell child
	c2->showerParameters()[1]=
	0.5*(sqr(c2->data().constituentMass())+sqr(c2->showerVariables()[2]))
	/(c2->showerParameters()[0]*p_dot_n());
	c2->set5Momentum(sudakov2Momentum(c2->showerParameters()[0],
	c2->showerParameters()[1],
	c2->showerVariables()[0],
	c2->showerVariables()[1],0));
	// spacelike child
	Lorentz5Momentum pc1(theParent->momentum() - c2->momentum());
	pc1.rescaleMass();
	c1->set5Momentum(pc1);
	}

	void IS_QTildeShowerKinematics1to2::
	updateLast( const tShowerParticlePtr theLast,Energy px,Energy py) const {
	if(theLast->isFinalState()) return;
	Energy2 pt2=sqr(px)+sqr(py);
	theLast->showerParameters()[0]=theLast->x();
	theLast->showerParameters()[1]=0.5*pt2/
	theLast->showerParameters()[0]/p_dot_n();
	theLast->showerVariables().resize(3);
	theLast->showerParameters().resize(2);
	for(unsigned int ix=0;ix<3;++ix) theLast->showerVariables()[ix]=ZERO;
	// momentum
	Lorentz5Momentum ntemp=Lorentz5Momentum(ZERO,-pVector().vect());
	double beta = 0.5*pt2/
	theLast->showerParameters()[0]/(pVector()*ntemp);
	Lorentz5Momentum plast =
	Lorentz5Momentum(pVector().z()>ZERO ? px : -px ,py,ZERO,ZERO)
	+theLast->x()pVector()+betantemp;
	plast.rescaleMass();
	theLast->set5Momentum(plast);
	}

	void IS_QTildeShowerKinematics1to2::initialize(ShowerParticle & particle, PPtr parent) {
	// For the time being we are considering only 1->2 branching
	Lorentz5Momentum p, n, pthis, pcm;
	assert(particle.perturbative()!=2);
	if(particle.perturbative()==1) {
	// find the partner and its momentum
	ShowerParticlePtr partner=particle.partner();
	if(partner->isFinalState()) {
	Lorentz5Momentum pa = -particle.momentum()+partner->momentum();
	Lorentz5Momentum pb = particle.momentum();
	Energy scale=parent->momentum().t();
	Lorentz5Momentum pbasis(ZERO,parent->momentum().vect().unit()*scale);
	Axis axis(pa.vect().unit());
	LorentzRotation rot;
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	if(axis.perp2()>1e-20) {
	rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	rot.rotateX(Constants::pi);
	}
	if(abs(1.-pa.e()/pa.vect().mag())>1e-6) rot.boostZ( pa.e()/pa.vect().mag());
	pb *= rot;


	if(pb.perp2()/GeV2>1e-20) {
	Boost trans = -1./pb.e()*pb.vect();
	trans.setZ(0.);
	rot.boost(trans);
	}
	pbasis *=rot;
	rot.invert();
	n = rot*Lorentz5Momentum(ZERO,-pbasis.vect());
	p = rot*Lorentz5Momentum(ZERO, pbasis.vect());
	}
	else {
	pcm = parent->momentum();
	p = Lorentz5Momentum(ZERO, pcm.vect());
	n = Lorentz5Momentum(ZERO, -pcm.vect());
	}
	}
	else {
	p = dynamic_ptr_cast<ShowerParticlePtr>(particle.children()[0])
	->showerKinematics()->getBasis()[0];
	n = dynamic_ptr_cast<ShowerParticlePtr>(particle.children()[0])
	->showerKinematics()->getBasis()[1];
	}
	setBasis(p,n);
	}
	diff --git a/Shower/Default/QTildeShowerKinematics1to2.cc b/Shower/Default/QTildeShowerKinematics1to2.cc
	--- a/Shower/Default/QTildeShowerKinematics1to2.cc
	+++ b/Shower/Default/QTildeShowerKinematics1to2.cc
	@@ -1,102 +1,132 @@
	// -- C++ --
	//
	// QTildeShowerKinematics1to2.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 QTildeShowerKinematics1to2 class.
	//

	#include "QTildeShowerKinematics1to2.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	+#include "Herwig++/Shower/Base/ShowerParticle.h"
	+#include "ThePEG/Repository/CurrentGenerator.h"

	using namespace Herwig;

	vector<Lorentz5Momentum> QTildeShowerKinematics1to2::getBasis() const {
	vector<Lorentz5Momentum> dum;
	dum.push_back( _pVector );
	dum.push_back( _nVector );
	return dum;
	}

	void QTildeShowerKinematics1to2::setBasis(const Lorentz5Momentum &p,
	const Lorentz5Momentum & n) {
	_pVector=p;
	_nVector=n;
	}

	Lorentz5Momentum QTildeShowerKinematics1to2::
	sudakov2Momentum(double alpha, double beta, Energy px,
	Energy py,unsigned int iopt) const {
	if(isnan(beta)\|\|isinf(beta))
	throw Exception() << "beta infinite in "
	<< "QTildeShowerKinematics1to2::sudakov2Momentum()"
	<< Exception::eventerror;
	Lorentz5Momentum dq;
	if(iopt==0) {
	const Boost beta_bb = -(_pVector + _nVector).boostVector();
	Lorentz5Momentum p_bb = _pVector;
	Lorentz5Momentum n_bb = _nVector;
	p_bb.boost( beta_bb );
	n_bb.boost( beta_bb );
	// set first in b2b frame along z-axis (assuming that p and n are
	// b2b as checked above)
	dq=Lorentz5Momentum(ZERO, ZERO, (alpha - beta)*p_bb.vect().mag(),
	alphap_bb.t() + betan_bb.t());
	// add transverse components
	dq.setX(px);
	dq.setY(py);
	// rotate to have z-axis parallel to p
	// this rotation changed by PR to a different rotation with the same effect
	// but different azimuthal angle to make implementing spin correlations easier
	// dq.rotateUz( unitVector(p_bb.vect()) );
	Axis axis(p_bb.vect().unit());
	if(axis.perp2()>0.) {
	LorentzRotation rot;
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	dq.transform(rot);
	}
	else if(axis.z()<0.) {
	dq.setZ(-dq.z());
	}
	// boost back
	dq.boost( -beta_bb );
	dq.rescaleMass();
	// return the momentum
	}
	else {
	const Boost beta_bb = -pVector().boostVector();
	Lorentz5Momentum p_bb = pVector();
	Lorentz5Momentum n_bb = nVector();
	p_bb.boost( beta_bb );
	n_bb.boost( beta_bb );
	// set first in b2b frame along z-axis (assuming that p and n are
	// b2b as checked above)
	dq=Lorentz5Momentum (ZERO, ZERO, 0.5betapVector().mass(),
	alphapVector().mass() + 0.5beta*pVector().mass());
	// add transverse components
	dq.setX(px);
	dq.setY(py);
	// changed to be same as other case
	// dq.rotateUz( unitVector(n_bb.vect()) );
	Axis axis(n_bb.vect().unit());
	if(axis.perp2()>0.) {
	LorentzRotation rot;
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	dq.transform(rot);
	}
	else if(axis.z()<0.) {
	dq.setZ(-dq.z());
	}
	// boost back
	dq.boost( -beta_bb );
	dq.rescaleMass();
	}
	return dq;
	}
	+
	+void QTildeShowerKinematics1to2::printScales(tShowerParticlePtr parent,
	+ tShowerParticlePtr child1,
	+ tShowerParticlePtr child2) const {
	+ CurrentGenerator::log() << parent << "\n" << child1 << " "
	+ << *child2 << "\n";
	+ CurrentGenerator::log() << "testing parent\n";
	+ for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	+ it = parent->evolutionScales().begin();
	+ it!=parent->evolutionScales().end();++it) {
	+ CurrentGenerator::log() << it->first << " " << it->second.first/GeV << " "
	+ << it->second.second/GeV << "\n";
	+ }
	+ CurrentGenerator::log() << "testing child[0]\n";
	+ for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	+ it = child1->evolutionScales().begin();
	+ it!=child1->evolutionScales().end();++it) {
	+ CurrentGenerator::log() << it->first << " " << it->second.first/GeV << " "
	+ << it->second.second/GeV << "\n";
	+ }
	+ CurrentGenerator::log() << "testing child[1]\n";
	+ for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
	+ it = child2->evolutionScales().begin();
	+ it!=child2->evolutionScales().end();++it) {
	+ CurrentGenerator::log() << it->first << " " << it->second.first/GeV << " "
	+ << it->second.second/GeV << "\n";
	+ }
	+}
	diff --git a/Shower/Default/QTildeShowerKinematics1to2.h b/Shower/Default/QTildeShowerKinematics1to2.h
	--- a/Shower/Default/QTildeShowerKinematics1to2.h
	+++ b/Shower/Default/QTildeShowerKinematics1to2.h
	@@ -1,108 +1,114 @@
	// -- C++ --
	//
	// QTildeShowerKinematics1to2.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_QTildeShowerKinematics1to2_H
	#define HERWIG_QTildeShowerKinematics1to2_H
	//
	// This is the declaration of the QTildeShowerKinematics1to2 class.
	//

	#include "Herwig++/Shower/Base/ShowerKinematics.h"
	#include "ThePEG/Vectors/Lorentz5Vector.h"
	#include "QTildeShowerKinematics1to2.fh"

	namespace Herwig {

	using namespace ThePEG;

	/** \ingroup Shower
	*
	* This abstract class describes the common features for initial and final
	* state radiation kinematics for \f$1\to2\f$ branchings and for
	* the choice of \f$\tilde{q}\f$ as evolution variable.
	*
	* @see ShowerKinematics
	* @see IS_QTildeShowerKinematics1to2
	* @see FS_QTildeShowerKinematics1to2
	* @see KinematicsReconstructor
	*/
	class QTildeShowerKinematics1to2: public ShowerKinematics {

	public:

	/**
	* Implementation of the virtual function returning a set of basis vectors, specific to
	* the type of evolution. This function will be used by the
	* ForwardShowerEvolver in order to access \f$p\f$
	* and \f$n\f$.
	*/
	virtual vector<Lorentz5Momentum> getBasis() const;

	/**
	* Access to the \f$p\f$ vector used to describe the kinematics.
	*/
	const Lorentz5Momentum & pVector() const {return _pVector;}

	/**
	* Access to the \f$n\f$ vector used to describe the kinematics.
	*/
	const Lorentz5Momentum & nVector() const {return _nVector;}

	/**
	* Dot product of thew basis vectors
	*/
	Energy2 p_dot_n() const {return _pVector*_nVector;}

	/**
	* Converts a Sudakov parametrization of a momentum w.r.t. the given
	* basis \f$p\f$ and \f$n\f$ into a 5 momentum.
	* @param alpha The \f$\alpha\f$ parameter of the Sudakov parameterisation
	* @param beta The \f$\beta\f$ parameter of the Sudakov parameterisation
	* @param px The \f$x\f$-component of the transverse momentum in the Sudakov
	* parameterisation
	* @param py The \f$x\f$-component of the transverse momentum in the Sudakov
	* parameterisation
	* @param iopt The option for the momentum reconstruction
	* - 0 is in the rest frame of the pair of reference vectors
	* - 1 is in the rest frame of the p vector
	*/
	Lorentz5Momentum sudakov2Momentum(double alpha, double beta, Energy px, Energy py,
	unsigned int iopt) const;

	protected:

	/**
	* Set the basis vectors
	*/
	void setBasis(const Lorentz5Momentum &p, const Lorentz5Momentum & n);

	+ /**
	+ * Print out scales after branching for debugging
	+ */
	+ void printScales(tShowerParticlePtr parent,tShowerParticlePtr child1,
	+ tShowerParticlePtr child2) const;
	+
	private:

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

	private:

	/**
	* The \f$p\f$ reference vector
	*/
	Lorentz5Momentum _pVector;

	/**
	* The \f$n\f$ reference vector
	*/
	Lorentz5Momentum _nVector;

	};

	}

	#endif /* HERWIG_QTildeShowerKinematics1to2_H */
	diff --git a/Shower/SplittingFunctions/SplittingGenerator.cc b/Shower/SplittingFunctions/SplittingGenerator.cc
	--- a/Shower/SplittingFunctions/SplittingGenerator.cc
	+++ b/Shower/SplittingFunctions/SplittingGenerator.cc
	@@ -1,407 +1,470 @@
	// -- 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;
	newKin = cit->second.first->
	generateNextTimeBranching(particle.evolutionScale(angularOrdered,type),
	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
	+ // 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->
	generateNextTimeBranching(particle.evolutionScale(angularOrdered,type),
	cit->second.second,particle.id()!=cit->first,
	0.5*enhance);
	ShoKinPtr newKin2 = cit->second.first->
	generateNextTimeBranching(particle.evolutionScale(angularOrdered,
	ShowerPartnerType::QCDAntiColourLine),
	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
	else {
	Energy startingScale =
	max(particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDColourLine),
	particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDAntiColourLine));
	newKin= cit->second.first->
	generateNextTimeBranching(startingScale,
	cit->second.second,particle.id()!=cit->first,
	enhance);
	- type = UseRandom::rndbool() ?
	+ type = cit->second.second[0]<0 ?
	ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
	}
	}
	+ // everything else
	else {
	type = particle.hasColour() ?
	ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
	newKin= cit->second.first->
	generateNextTimeBranching(particle.evolutionScale(angularOrdered,type),
	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, Energy stoppingScale,
	Energy minmass, double enhance,
	ShowerInteraction::Type type) const {
	Energy newQ = Constants::MaxEnergy;
	ShoKinPtr kinematics;
	SudakovPtr sudakov;
	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) {
	if(type!=cit->second.first->interactionType()) continue;
	ShoKinPtr newKin;
	if(particle.evolutionScale() < stoppingScale)
	newKin = cit->second.first->
	generateNextDecayBranching(particle.evolutionScale(),
	stoppingScale,minmass,
	cit->second.second,
	particle.id()!=cit->first,enhance);
	if(!newKin) continue;
	if(newKin->scale() < newQ && newKin->scale() > particle.evolutionScale()) {
	newQ = newKin->scale();
	ids = cit->second.second;
	kinematics=newKin;
	sudakov=cit->second.first;
	}
	}
	// 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
	assert(false);
	return Branching(kinematics, ids,sudakov,ShowerPartnerType::Undefined);
	}

	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(),sudakov,
	- ShowerPartnerType::Undefined);
	- // select the branching
	+ 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 ) {
	- if(type!=cit->second.first->interactionType()) continue;
	+ // 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);
	- ShoKinPtr newKin=cit->second.first->
	- generateNextSpaceBranching(particle.evolutionScale(),
	- cit->second.second, particle.x(),
	- particle.id()!=cit->first,enhance,beam);
	+ // 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);
	+ // 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
	- assert(false);
	- return Branching(kinematics, ids,sudakov,ShowerPartnerType::Undefined);
	+ 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;
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:31 PM (1 d, 11 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805842
Default Alt Text: (167 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions