No OneTemporary
Actions

Size

210 KB

Subscribers

None

View Options

	diff --git a/Shower/Base/Evolver.cc b/Shower/Base/Evolver.cc
	--- a/Shower/Base/Evolver.cc
	+++ b/Shower/Base/Evolver.cc
	@@ -1,3220 +1,3225 @@
	// -- 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"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ShowerVertex.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
	#include "ThePEG/Handlers/StandardXComb.h"

	using namespace Herwig;

	namespace {

	/**
	* A struct to order the particles in the same way as in the DecayMode's
	*/
	struct ParticleOrdering {
	/**
	* Operator for the ordering
	* @param p1 The first ParticleData object
	* @param p2 The second ParticleData object
	*/
	bool operator() (tcPDPtr p1, tcPDPtr p2) {
	return abs(p1->id()) > abs(p2->id()) \|\|
	( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) \|\|
	( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
	}
	};
	typedef multiset<tcPDPtr,ParticleOrdering> OrderedParticles;

	/**
	* Cached lookup of decay modes.
	* Generator::findDecayMode() is not efficient.
	*/
	tDMPtr findDecayMode(const string & tag) {
	static map<string,DMPtr> cache;
	map<string,DMPtr>::const_iterator pos = cache.find(tag);

	if ( pos != cache.end() )
	return pos->second;

	tDMPtr dm = CurrentGenerator::current().findDecayMode(tag);
	cache[tag] = dm;
	return dm;
	}
	}

	DescribeClass<Evolver,Interfaced>
	describeEvolver ("Herwig::Evolver","HwShower.so");

	bool Evolver::_hardEmissionModeWarn = true;
	bool Evolver::_missingTruncWarn = true;

	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 << _spinOpt << _softOpt << _hardPOWHEG
	<< ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
	<< _vetoes << _trunc_Mode << _hardEmissionMode << _reconOpt
	<< isMCatNLOSEvent << isMCatNLOHEvent
	<< isPowhegSEvent << isPowhegHEvent
	<< theFactorizationScaleFactor << theRenormalizationScaleFactor << ounit(muPt,GeV)
	<< oenum(interaction_) << _maxTryFSR << _maxFailFSR << _fracFSR;
	}

	void Evolver::persistentInput(PersistentIStream & is, int) {
	is >> _model >> _splittingGenerator >> _maxtry
	>> _meCorrMode >> _hardVetoMode >> _hardVetoRead >> _hardVetoReadOption
	>> _limitEmissions >> _spinOpt >> _softOpt >> _hardPOWHEG
	>> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
	>> _vetoes >> _trunc_Mode >> _hardEmissionMode >> _reconOpt
	>> isMCatNLOSEvent >> isMCatNLOHEvent
	>> isPowhegSEvent >> isPowhegHEvent
	>> theFactorizationScaleFactor >> theRenormalizationScaleFactor >> iunit(muPt,GeV)
	>> ienum(interaction_) >> _maxTryFSR >> _maxFailFSR >> _fracFSR;
	}

	void Evolver::doinit() {
	Interfaced::doinit();
	// calculate max no of FSR vetos
	_maxFailFSR = max(int(_maxFailFSR), int(_fracFSR*double(generator()->N())));
	}

	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,int> interfaceHardEmissionMode
	("HardEmissionMode",
	"Whether to use ME corrections or POWHEG for the hardest emission",
	&Evolver::_hardEmissionMode, 0, false, false);
	static SwitchOption interfaceHardEmissionModeDecayMECorrection
	(interfaceHardEmissionMode,
	"DecayMECorrection",
	"Old fashioned ME correction for decays only",
	-1);
	static SwitchOption interfaceHardEmissionModeMECorrection
	(interfaceHardEmissionMode,
	"MECorrection",
	"Old fashioned ME correction",
	0);
	static SwitchOption interfaceHardEmissionModePOWHEG
	(interfaceHardEmissionMode,
	"POWHEG",
	"Powheg style hard emission using internal matrix elements",
	1);
	static SwitchOption interfaceHardEmissionModeMatchboxPOWHEG
	(interfaceHardEmissionMode,
	"MatchboxPOWHEG",
	"Powheg style emission for the hard process using Matchbox",
	2);
	static SwitchOption interfaceHardEmissionModeFullPOWHEG
	(interfaceHardEmissionMode,
	"FullPOWHEG",
	"Powheg style emission for the hard process using Matchbox "
	"and decays using internal matrix elements",
	3);

	static Switch<Evolver,ShowerInteraction::Type > interfaceInteraction
	("Interaction",
	"The interactions to be used in the shower",
	&Evolver::interaction_, ShowerInteraction::QEDQCD, false, false);
	static SwitchOption interfaceInteractionQCD
	(interfaceInteraction,
	"QCDOnly",
	"Only QCD",
	ShowerInteraction::QCD);
	static SwitchOption interfaceInteractionQEDOnly
	(interfaceInteraction,
	"QEDOnly",
	"Only QED",
	ShowerInteraction::QED);
	- static SwitchOption interfaceInteractionBothAtOnce
	+ static SwitchOption interfaceInteractionQEDQCD
	(interfaceInteraction,
	"QEDQCD",
	"QED and QCD",
	ShowerInteraction::QEDQCD);
	+ static SwitchOption interfaceInteractionALL
	+ (interfaceInteraction,
	+ "ALL",
	+ "QED, QCD and EW",
	+ ShowerInteraction::ALL);

	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 veto the emission of the parent,"
	" no veto in generation of emissions from children",
	1);
	static SwitchOption interfaceReconstructionOptionOffShell2
	(interfaceReconstructionOption,
	"OffShell2",
	"Use the off-shell masses in the calculation veto the emissions from the children."
	" no veto in generation of emissions from children",
	2);
	static SwitchOption interfaceReconstructionOptionOffShell3
	(interfaceReconstructionOption,
	"OffShell3",
	"Use the off-shell masses in the calculation veto the emissions from the children."
	" veto in generation of emissions from children using cut-off for second parton",
	3);

	static Switch<Evolver,unsigned int> interfaceSpinCorrelations
	("SpinCorrelations",
	"Treatment of spin correlations in the parton shower",
	&Evolver::_spinOpt, 1, false, false);
	static SwitchOption interfaceSpinCorrelationsOff
	(interfaceSpinCorrelations,
	"No",
	"No spin correlations",
	0);
	static SwitchOption interfaceSpinCorrelationsSpin
	(interfaceSpinCorrelations,
	"Yes",
	"Include the azimuthal spin correlations only",
	1);

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

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

	static Parameter<Evolver,unsigned int> interfaceMaxTryFSR
	("MaxTryFSR",
	"The maximum number of attempted FSR emissions in"
	" the generation of the FSR",
	&Evolver::_maxTryFSR, 100000, 10, 100000000,
	false, false, Interface::limited);

	static Parameter<Evolver,unsigned int> interfaceMaxFailFSR
	("MaxFailFSR",
	"Maximum number of failures generating the FSR",
	&Evolver::_maxFailFSR, 100, 1, 100000000,
	false, false, Interface::limited);


	static Parameter<Evolver,double> interfaceFSRFailureFraction
	("FSRFailureFraction",
	"Maximum fraction of events allowed to fail due to too many FSR emissions",
	&Evolver::_fracFSR, 0.001, 1e-10, 1,
	false, false, Interface::limited);
	}

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

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

	void Evolver::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) {


	isMCatNLOSEvent = false;
	isMCatNLOHEvent = false;
	isPowhegSEvent = false;
	isPowhegHEvent = false;

	Ptr<SubtractedME>::tptr subme;
	Ptr<MatchboxMEBase>::tptr me;
	Ptr<SubtractionDipole>::tptr dipme;

	Ptr<StandardXComb>::ptr sxc = dynamic_ptr_cast<Ptr<StandardXComb>::ptr>(xcomb);

	if ( sxc ) {
	subme = dynamic_ptr_cast<Ptr<SubtractedME>::tptr>(sxc->matrixElement());
	me = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(sxc->matrixElement());
	dipme = dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(sxc->matrixElement());
	}

	if ( subme ) {
	if ( subme->showerApproximation() ) {
	theShowerApproximation = subme->showerApproximation();
	// separate MCatNLO and POWHEG-type corrections
	if ( !subme->showerApproximation()->needsSplittingGenerator() ) {
	if ( subme->realShowerSubtraction() )
	isMCatNLOHEvent = true;
	else if ( subme->virtualShowerSubtraction() )
	isMCatNLOSEvent = true;
	}
	else {
	if ( subme->realShowerSubtraction() )
	isPowhegHEvent = true;
	else if ( subme->virtualShowerSubtraction() \|\| subme->loopSimSubtraction() )
	isPowhegSEvent = true;
	}
	}
	} else if ( me ) {
	if ( me->factory()->showerApproximation() ) {
	theShowerApproximation = me->factory()->showerApproximation();
	if ( !me->factory()->showerApproximation()->needsSplittingGenerator() )
	isMCatNLOSEvent = true;
	else
	isPowhegSEvent = true;
	}
	}

	string error = "Inconsistent hard emission set-up in Evolver::showerHardProcess(). ";
	if ( ( isMCatNLOSEvent \|\| isMCatNLOHEvent ) ){
	if (_hardEmissionMode > 1)
	throw Exception() << error
	<< "Cannot generate POWHEG matching with MC@NLO shower "
	<< "approximation. Add 'set Evolver:HardEmissionMode 0' to input file."
	<< Exception::runerror;
	if ( ShowerHandler::currentHandler()->canHandleMatchboxTrunc())
	throw Exception() << error
	<< "Cannot use truncated qtilde shower with MC@NLO shower "
	<< "approximation. Set LHCGenerator:EventHandler"
	<< ":CascadeHandler to '/Herwig/Shower/ShowerHandler' or "
	<< "'/Herwig/DipoleShower/DipoleShowerHandler'."
	<< Exception::runerror;
	}
	else if ( ((isPowhegSEvent \|\| isPowhegHEvent) \|\| dipme) &&
	_hardEmissionMode < 2){
	if ( ShowerHandler::currentHandler()->canHandleMatchboxTrunc())
	throw Exception() << error
	<< "Unmatched events requested for POWHEG shower "
	<< "approximation. Set Evolver:HardEmissionMode to "
	<< "'MatchboxPOWHEG' or 'FullPOWHEG'."
	<< Exception::runerror;
	else if (_hardEmissionModeWarn){
	_hardEmissionModeWarn = false;
	_hardEmissionMode+=2;
	throw Exception() << error
	<< "Unmatched events requested for POWHEG shower "
	<< "approximation. Changing Evolver:HardEmissionMode from "
	<< _hardEmissionMode-2 << " to " << _hardEmissionMode
	<< Exception::warning;
	}
	}

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

	_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,xcomb);
	// 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(hard)) 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(hard))
	_hardme->applyHardMatrixElementCorrection(_currenttree);
	}
	}
	else {
	if(_decayme&&_decayme->hasMECorrection()) {
	_decayme->initializeMECorrection(_currenttree,
	_initialenhance,_finalenhance);
	if(hardMEC(hard))
	_decayme->applyHardMatrixElementCorrection(_currenttree);
	}
	}
	}

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

	bool Evolver::timeLikeShower(tShowerParticlePtr particle,
	ShowerInteraction::Type type,
	Branching fb, bool first) {
	// don't do anything if not needed
	if(_limitEmissions == 1 \|\| hardOnly() \|\|
	( _limitEmissions == 2 && _nfs != 0) \|\|
	( _limitEmissions == 4 && _nfs + _nis != 0) ) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	// too many tries
	if(_nFSR>=_maxTryFSR) {
	++_nFailedFSR;
	// too many failed events
	if(_nFailedFSR>=_maxFailFSR)
	throw Exception() << "Too many events have failed due to too many shower emissions, in\n"
	<< "Evolver::timeLikeShower(). Terminating run\n"
	<< Exception::runerror;
	throw Exception() << "Too many attempted emissions in Evolver::timeLikeShower()\n"
	<< Exception::eventerror;
	}
	// generate the emission
	ShowerParticleVector children;
	int ntry=0;
	// generate the emission
	if(!fb.kinematics)
	fb = selectTimeLikeBranching(particle,type,HardBranchingPtr());
	// no emission, return
	if(!fb.kinematics) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	Branching fc[2];
	bool setupChildren = true;
	while (ntry<50) {
	fc[0] = Branching();
	fc[1] = Branching();
	++ntry;
	assert(fb.kinematics);
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	if(setupChildren) {
	++_nFSR;
	particle->showerKinematics(fb.kinematics);
	// generate phi
	fb.kinematics->phi(fb.sudakov->generatePhiForward(*particle,fb.ids,fb.kinematics));
	// check highest pT
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the children
	children = createTimeLikeChildren(particle,fb.ids);
	// update the children
	particle->showerKinematics()->
	updateChildren(particle, children,fb.type,_reconOpt>=3);
	// update number of emissions
	++_nfs;
	if(_limitEmissions!=0) {
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}
	setupChildren = false;
	}
	// select branchings for children
	fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
	fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
	// old default
	if(_reconOpt==0) {
	// shower the first particle
	if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	break;
	}
	// Herwig default
	else if(_reconOpt==1) {
	// shower the first particle
	if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<children.size();++ix)
	particle->abandonChild(children[ix]);
	children.clear();
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	// generate the new emission
	fb = selectTimeLikeBranching(particle,type,HardBranchingPtr());
	// no emission, return
	if(!fb.kinematics) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	setupChildren = true;
	continue;
	}
	else
	break;
	}
	// veto children
	else if(_reconOpt>=2) {
	// cut-off masses for the branching
	const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
	// compute the masses of the children
	Energy masses[3];
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].kinematics) {
	const vector<Energy> & vm = fc[ix].sudakov->virtualMasses(fc[ix].ids);
	Energy2 q2 =
	fc[ix].kinematics->z()(1.-fc[ix].kinematics->z())sqr(fc[ix].kinematics->scale());
	if(fc[ix].ids[0]!=ParticleID::g) q2 += sqr(vm[0]);
	masses[ix+1] = sqrt(q2);
	}
	else {
	masses[ix+1] = virtualMasses[ix+1];
	}
	}
	masses[0] = fb.ids[0]!=ParticleID::g ? virtualMasses[0] : ZERO;
	double z = fb.kinematics->z();
	Energy2 pt2 = z(1.-z)(z(1.-z)sqr(fb.kinematics->scale()) + sqr(masses[0]))
	- sqr(masses[1])(1.-z) - sqr(masses[2])z;
	if(pt2>=ZERO) {
	break;
	}
	else {
	// reset the scales for the children
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].kinematics)
	children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
	else
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
	children[ix]->virtualMass(ZERO);
	}
	}
	}
	};
	if(_reconOpt>=2) {
	// shower the first particle
	if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	}
	if(first&&!children.empty())
	particle->showerKinematics()->resetChildren(particle,children);
	if(particle->spinInfo()) particle->spinInfo()->develop();
	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 ) ) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	Branching bb;
	// generate branching
	while (true) {
	bb=_splittingGenerator->chooseBackwardBranching(*particle,beam,
	_initialenhance,
	_beam,type,
	pdf,freeze);
	// return if no emission
	if(!bb.kinematics) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	// if not vetoed break
	if(!spaceLikeVetoed(bb,particle)) break;
	// otherwise reset scale and continue
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	}
	// assign the splitting function and shower kinematics
	particle->showerKinematics(bb.kinematics);
	if(bb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(bb.kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// particles as in Sudakov form factor
	tcPDPtr part[2]={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);
	// update the history if needed
	_currenttree->updateInitialStateShowerProduct(_progenitor,newParent);
	_currenttree->addInitialStateBranching(particle,newParent,otherChild);
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	++_nis;
	bool emitted = _limitEmissions==0 ?
	spaceLikeShower(newParent,beam,type) : false;
	if(newParent->spinInfo()) newParent->spinInfo()->develop();
	// now reconstruct the momentum
	if(!emitted) {
	if(_intrinsic.find(_progenitor)==_intrinsic.end()) {
	bb.kinematics->updateLast(newParent,ZERO,ZERO);
	}
	else {
	pair<Energy,double> kt=_intrinsic[_progenitor];
	bb.kinematics->updateLast(newParent,
	kt.first*cos(kt.second),
	kt.first*sin(kt.second));
	}
	}
	particle->showerKinematics()->
	updateChildren(newParent, theChildren,bb.type,false);
	if(_limitEmissions!=0) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}
	// perform the shower of the final-state particle
	timeLikeShower(otherChild,type,Branching(),true);
	updateHistory(otherChild);
	if(theChildren[1]->spinInfo()) theChildren[1]->spinInfo()->develop();
	// return the emitted
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}

	void 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()
	+ "->";
	OrderedParticles outgoing;
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) {
	if(abs(decay->incomingLines().begin()->first->original()->id()) == ParticleID::t &&
	abs(it->first->original()->id())==ParticleID::Wplus &&
	decay->treelinks().size() == 1) {
	ShowerTreePtr Wtree = decay->treelinks().begin()->first;
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	it2=Wtree->outgoingLines().begin();it2!=Wtree->outgoingLines().end();++it2) {
	outgoing.insert(it2->first->original()->dataPtr());
	}
	}
	else {
	outgoing.insert(it->first->original()->dataPtr());
	}
	}
	for(OrderedParticles::const_iterator it=outgoing.begin(); it!=outgoing.end();++it) {
	if(it!=outgoing.begin()) tag += ",";
	tag +=(**it).name();
	}
	tag += ";";
	dm = findDecayMode(tag);
	}
	if(dm) _decayme = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
	// set the ShowerTree to be showered
	currentTree(decay);
	decay->applyTransforms();
	hardTree(HardTreePtr());
	unsigned int interactionTry=0;
	do {
	try {
	// generate the showering
	doShowering(false,XCPtr());
	// if no vetos
	// force calculation of spin correlations
	SpinPtr spInfo = decay->incomingLines().begin()->first->progenitor()->spinInfo();
	if(spInfo) {
	if(!spInfo->developed()) spInfo->needsUpdate();
	spInfo->develop();
	}
	// and then return
	return;
	}
	catch (InteractionVeto) {
	currentTree()->clear();
	++interactionTry;
	}
	}
	while(interactionTry<=5);
	throw Exception() << "Too many tries for QED shower in Evolver::showerDecay()"
	<< Exception::eventerror;
	}

	bool Evolver::spaceLikeDecayShower(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass,ShowerInteraction::Type type,
	Branching fb) {
	// too many tries
	if(_nFSR>=_maxTryFSR) {
	++_nFailedFSR;
	// too many failed events
	if(_nFailedFSR>=_maxFailFSR)
	throw Exception() << "Too many events have failed due to too many shower emissions, in\n"
	<< "Evolver::timeLikeShower(). Terminating run\n"
	<< Exception::runerror;
	throw Exception() << "Too many attempted emissions in Evolver::timeLikeShower()\n"
	<< Exception::eventerror;
	}
	// generate the emission
	ShowerParticleVector children;
	int ntry=0;
	// generate the emission
	if(!fb.kinematics)
	fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,
	HardBranchingPtr());
	// no emission, return
	if(!fb.kinematics) return false;
	Branching fc[2];
	bool setupChildren = true;
	while (ntry<50) {
	if(particle->virtualMass()==ZERO)
	particle->virtualMass(_progenitor->progenitor()->mass());
	fc[0] = Branching();
	fc[1] = Branching();
	++ntry;
	assert(fb.kinematics);
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	if(setupChildren) {
	++_nFSR;
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the ShowerParticle objects for the two children
	children = createTimeLikeChildren(particle,fb.ids);
	// updateChildren the children
	particle->showerKinematics()->
	updateChildren(particle, children, fb.type,_reconOpt>=3);
	setupChildren = false;
	}
	// select branchings for children
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,
	type,HardBranchingPtr());
	fc[1] = selectTimeLikeBranching (children[1],type,HardBranchingPtr());
	// old default
	if(_reconOpt==0) {
	// shower the first particle
	_currenttree->updateInitialStateShowerProduct(_progenitor,children[0]);
	_currenttree->addInitialStateBranching(particle,children[0],children[1]);
	if(fc[0].kinematics) spaceLikeDecayShower(children[0],maxScales,minmass,type,Branching());
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],true);
	updateHistory(children[1]);
	// branching has happened
	break;
	}
	// Herwig default
	else if(_reconOpt==1) {
	// shower the first particle
	_currenttree->updateInitialStateShowerProduct(_progenitor,children[0]);
	_currenttree->addInitialStateBranching(particle,children[0],children[1]);
	if(fc[0].kinematics) spaceLikeDecayShower(children[0],maxScales,minmass,type,Branching());
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],true);
	updateHistory(children[1]);
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<children.size();++ix)
	particle->abandonChild(children[ix]);
	children.clear();
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	// generate the new emission
	fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,
	HardBranchingPtr());
	// no emission, return
	if(!fb.kinematics) {
	return false;
	}
	setupChildren = true;
	continue;
	}
	else
	break;
	}
	else if(_reconOpt>=2) {
	// cut-off masses for the branching
	const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
	// compute the masses of the children
	Energy masses[3];
	// space-like children
	masses[1] = children[0]->virtualMass();
	// time-like child
	if(fc[1].kinematics) {
	const vector<Energy> & vm = fc[1].sudakov->virtualMasses(fc[1].ids);
	Energy2 q2 =
	fc[1].kinematics->z()(1.-fc[1].kinematics->z())sqr(fc[1].kinematics->scale());
	if(fc[1].ids[0]!=ParticleID::g) q2 += sqr(vm[0]);
	masses[2] = sqrt(q2);
	}
	else {
	masses[2] = virtualMasses[2];
	}
	masses[0]=particle->virtualMass();
	double z = fb.kinematics->z();
	Energy2 pt2 = (1.-z)(zsqr(masses[0])-sqr(masses[1])-z/(1.-z)*sqr(masses[2]));
	if(pt2>=ZERO) {
	break;
	}
	else {
	// reset the scales for the children
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].kinematics)
	children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
	else {
	if(ix==0)
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,Constants::MaxEnergy);
	else
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
	}
	}
	children[0]->virtualMass(_progenitor->progenitor()->mass());
	children[1]->virtualMass(ZERO);
	}
	}
	};
	if(_reconOpt>=2) {
	// 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,children[0]);
	_currenttree->addInitialStateBranching(particle,children[0],children[1]);
	// shower the first particle
	if(fc[0].kinematics) spaceLikeDecayShower(children[0],maxScales,minmass,type,Branching());
	// shower the second particle
	if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],true);
	updateHistory(children[1]);
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	}
	// branching has happened
	return true;
	}

	vector<ShowerProgenitorPtr> Evolver::setupShower(bool hard) {
	// generate POWHEG hard emission if needed
	if(_hardEmissionMode>0) hardestEmission(hard);
	ShowerInteraction::Type inter = interaction_;
	if(_hardtree&&inter!=ShowerInteraction::QEDQCD && inter!=ShowerInteraction::ALL) {
	inter = _hardtree->interaction();
	}
	// set the initial colour partners
	setEvolutionPartners(hard,inter,false);
	// generate hard me if needed
	if(_hardEmissionMode==0 \|\|
	(!hard && _hardEmissionMode==-1)) hardMatrixElementCorrection(hard);
	// get the particles to be showered
	vector<ShowerProgenitorPtr> particlesToShower =
	currentTree()->extractProgenitors();
	// remake the colour partners if needed
	if(_currenttree->hardMatrixElementCorrection()) {
	setEvolutionPartners(hard,interaction_,true);
	_currenttree->resetShowerProducts();
	}
	// return the answer
	return particlesToShower;
	}

	void Evolver::setEvolutionPartners(bool hard,ShowerInteraction::Type type,
	bool clear) {
	// 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();
	// clear the partners if needed
	if(clear) {
	for(unsigned int ix=0;ix<particles.size();++ix) {
	particles[ix]->partner(ShowerParticlePtr());
	particles[ix]->clearPartners();
	}
	}
	// sort out the colour partners
	if(hardTree()) {
	// find the partner
	for(unsigned int ix=0;ix<particles.size();++ix) {
	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);
	if(hardTree() && _hardPOWHEG) {
	bool tooHard=false;
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit=hardTree()->particles().end();
	for(unsigned int ix=0;ix<particles.size();++ix) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	mit = hardTree()->particles().find(particles[ix]);
	Energy hardScale(ZERO);
	ShowerPartnerType::Type type(ShowerPartnerType::Undefined);
	// final-state
	if(particles[ix]->isFinalState()) {
	if(mit!= eit && !mit->second->children().empty()) {
	hardScale = mit->second->scale();
	type = mit->second->type();
	}
	}
	// initial-state
	else {
	if(mit!= eit && mit->second->parent()) {
	hardScale = mit->second->parent()->scale();
	type = mit->second->parent()->type();
	}
	}
	if(type!=ShowerPartnerType::Undefined) {
	if(type==ShowerPartnerType::QED) {
	tooHard \|= particles[ix]->scales().QED_noAO<hardScale;
	}
	else if(type==ShowerPartnerType::QCDColourLine) {
	tooHard \|= particles[ix]->scales().QCD_c_noAO<hardScale;
	}
	else if(type==ShowerPartnerType::QCDAntiColourLine) {
	tooHard \|= particles[ix]->scales().QCD_ac_noAO<hardScale;
	}
	else if(type==ShowerPartnerType::EW) {
	tooHard \|= particles[ix]->scales().EW<hardScale;
	}
	}
	}
	if(tooHard) convertHardTree(hard,type);
	}
	}

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

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

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

	bool Evolver::timeLikeVetoed(const Branching & fb,
	ShowerParticlePtr particle) {
	// work out type of interaction
	ShowerInteraction::Type type = convertInteraction(fb.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;
	}
	if ( ShowerHandler::currentHandler()->firstInteraction() &&
	ShowerHandler::currentHandler()->profileScales() ) {
	double weight =
	ShowerHandler::currentHandler()->profileScales()->
	hardScaleProfile(_progenitor->hardScale(),fb.kinematics->pT());
	if ( UseRandom::rnd() > weight )
	return true;
	}
	return false;
	}

	bool Evolver::spaceLikeVetoed(const Branching & bb,
	ShowerParticlePtr particle) {
	// work out type of interaction
	ShowerInteraction::Type type = convertInteraction(bb.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;
	}
	if ( ShowerHandler::currentHandler()->firstInteraction() &&
	ShowerHandler::currentHandler()->profileScales() ) {
	double weight =
	ShowerHandler::currentHandler()->profileScales()->
	hardScaleProfile(_progenitor->hardScale(),bb.kinematics->pT());
	if ( UseRandom::rnd() > weight )
	return true;
	}
	return false;
	}

	bool Evolver::spaceLikeDecayVetoed( const Branching & fb,
	ShowerParticlePtr particle) {
	// work out type of interaction
	ShowerInteraction::Type type = convertInteraction(fb.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) {
	HardTreePtr ISRTree;
	if( ( _hardme && _hardme->hasPOWHEGCorrection()!=0 && _hardEmissionMode< 2) \|\|
	( _decayme && _decayme->hasPOWHEGCorrection()!=0 && _hardEmissionMode!=2) ) {
	if(_hardme) {
	assert(hard);
	_hardtree = _hardme->generateHardest( currentTree(),interaction_);
	}
	else {
	assert(!hard);
	_hardtree = _decayme->generateHardest( currentTree() );
	}
	// store initial state POWHEG radiation
	if(_hardtree && _hardme && _hardme->hasPOWHEGCorrection()==1)
	ISRTree=_hardtree;
	}

	else if (_hardEmissionMode>1 && hard) {
	// Get minimum pT cutoff used in shower approximation
	Energy maxpt = 1.*GeV;
	int colouredIn = 0;
	int colouredOut = 0;
	for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
	= currentTree()->outgoingLines().begin();
	it != currentTree()->outgoingLines().end(); ++it ) {
	if( it->second->coloured() ) colouredOut+=1;
	}
	for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
	= currentTree()->incomingLines().begin();
	it != currentTree()->incomingLines().end(); ++it ) {
	if( ! it->second->coloured() ) colouredIn+=1;
	}

	if ( theShowerApproximation ){
	if ( theShowerApproximation->ffPtCut() == theShowerApproximation->fiPtCut() &&
	theShowerApproximation->ffPtCut() == theShowerApproximation->iiPtCut() )
	maxpt = theShowerApproximation->ffPtCut();
	else if ( colouredIn == 2 && colouredOut == 0 )
	maxpt = theShowerApproximation->iiPtCut();
	else if ( colouredIn == 0 && colouredOut > 1 )
	maxpt = theShowerApproximation->ffPtCut();
	else if ( colouredIn == 2 && colouredOut == 1 )
	maxpt = min(theShowerApproximation->iiPtCut(), theShowerApproximation->fiPtCut());
	else if ( colouredIn == 1 && colouredOut > 1 )
	maxpt = min(theShowerApproximation->ffPtCut(), theShowerApproximation->fiPtCut());
	else
	maxpt = min(min(theShowerApproximation->iiPtCut(), theShowerApproximation->fiPtCut()),
	theShowerApproximation->ffPtCut());
	}

	// Generate hardtree from born and real emission subprocesses
	_hardtree = ShowerHandler::currentHandler()->generateCKKW(currentTree());

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


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

	// Set maxpt to pT of emission when showering POWHEG real-emission subprocesses
	if (!isPowhegSEvent && !isPowhegHEvent){
	vector<int> outGluon;
	vector<int> outQuark;
	map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it;
	for( it = currentTree()->outgoingLines().begin();
	it != currentTree()->outgoingLines().end(); ++it ) {
	if ( abs(it->second->id())< 6) outQuark.push_back(it->second->id());
	if ( it->second->id()==21 ) outGluon.push_back(it->second->id());
	}
	if (outGluon.size() + outQuark.size() == 1){
	for( it = currentTree()->outgoingLines().begin();
	it != currentTree()->outgoingLines().end(); ++it ) {
	if ( abs(it->second->id())< 6 \|\| it->second->id()==21 )
	maxpt = it->second->momentum().perp();
	}
	}
	else if (outGluon.size() + outQuark.size() > 1){
	// assume qqbar pair from a Z/gamma
	if (outGluon.size()==1 && outQuark.size() == 2 && outQuark[0]==-outQuark[1]){
	for( it = currentTree()->outgoingLines().begin();
	it != currentTree()->outgoingLines().end(); ++it ) {
	if ( it->second->id()==21 )
	maxpt = it->second->momentum().perp();
	}
	}
	// otherwise take the lowest pT avoiding born DY events
	else {
	maxpt = generator()->maximumCMEnergy();
	for( it = currentTree()->outgoingLines().begin();
	it != currentTree()->outgoingLines().end(); ++it ) {
	if ( abs(it->second->id())< 6 \|\| it->second->id()==21 )
	maxpt = min(maxpt,it->second->momentum().perp());
	}
	}
	}
	}

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

	// if hard me doesn't have a FSR powheg
	// correction use decay powheg correction
	if (_hardme && _hardme->hasPOWHEGCorrection()<2) {
	// check for intermediate colour singlet resonance
	const ParticleVector inter = _hardme->subProcess()->intermediates();
	if (inter.size()!=1 \|\|
	inter[0]->momentum().m2()/GeV2 < 0 \|\|
	inter[0]->dataPtr()->iColour()!=PDT::Colour0){
	if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
	return;
	}

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

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

	// get the decayer
	HwDecayerBasePtr decayer;
	if(dm) decayer = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
	// check if decayer has a FSR POWHEG correction
	if (!decayer \|\| decayer->hasPOWHEGCorrection()<2){
	if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
	return;
	}

	// generate the hardest emission
	ShowerDecayMap decay;
	PPtr in = new_ptr(*inter[0]);
	ShowerTreePtr decayTree = new_ptr(ShowerTree(in, decay));
	HardTreePtr FSRTree = decayer->generateHardest(decayTree);
	if (!FSRTree) {
	if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
	return;
	}

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

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

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

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

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

	bool Evolver::truncatedTimeLikeShower(tShowerParticlePtr particle,
	HardBranchingPtr branch,
	ShowerInteraction::Type type,
	Branching fb, bool first) {
	// select a branching if we don't have one
	if(!fb.kinematics)
	fb = selectTimeLikeBranching(particle,type,branch);
	// must be an emission, the forced one it not a truncated one
	assert(fb.kinematics);
	ShowerParticleVector children;
	int ntry=0;
	Branching fc[2];
	bool setupChildren = true;
	while (ntry<50) {
	if(!fc[0].hard) fc[0] = Branching();
	if(!fc[1].hard) fc[1] = Branching();
	++ntry;
	// Assign the shower kinematics to the emitting particle.
	if(setupChildren) {
	++_nFSR;
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// if not hard generate phi
	if(!fb.hard)
	fb.kinematics->phi(fb.sudakov->generatePhiForward(*particle,fb.ids,fb.kinematics));
	// create the children
	children = createTimeLikeChildren(particle,fb.ids);
	// update the children
	particle->showerKinematics()->
	updateChildren(particle, children,fb.type,_reconOpt>=3);
	setupChildren = false;
	}
	// select branchings for children
	if(!fc[0].kinematics) {
	// select branching for first particle
	if(!fb.hard && fb.iout ==1 )
	fc[0] = selectTimeLikeBranching(children[0],type,branch);
	else if(fb.hard && !branch->children()[0]->children().empty() )
	fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]);
	else
	fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
	}
	// select branching for the second particle
	if(!fc[1].kinematics) {
	// select branching for first particle
	if(!fb.hard && fb.iout ==2 )
	fc[1] = selectTimeLikeBranching(children[1],type,branch);
	else if(fb.hard && !branch->children()[1]->children().empty() )
	fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]);
	else
	fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
	}
	// old default
	if(_reconOpt==0 \|\| (_reconOpt==1 && fb.hard) ) {
	// shower the first particle
	if(fc[0].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 1)
	truncatedTimeLikeShower(children[0],branch,type,fc[0],false);
	// hard emission and subsquent hard emissions
	else if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 2)
	truncatedTimeLikeShower(children[1],branch,type,fc[1],false);
	// hard emission and subsquent hard emissions
	else if(fb.hard && !branch->children()[1]->children().empty() )
	truncatedTimeLikeShower(children[1],branch->children()[1],type,fc[1],false);
	else
	timeLikeShower(children[1],type,fc[1],false);
	}
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	break;
	}
	// H7 default
	else if(_reconOpt==1) {
	// shower the first particle
	if(fc[0].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 1)
	truncatedTimeLikeShower(children[0],branch,type,fc[0],false);
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 2)
	truncatedTimeLikeShower(children[1],branch,type,fc[1],false);
	else
	timeLikeShower(children[1],type,fc[1],false);
	}
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<children.size();++ix)
	particle->abandonChild(children[ix]);
	children.clear();
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	// generate the new emission
	fb = selectTimeLikeBranching(particle,type,branch);
	// must be at least hard emission
	assert(fb.kinematics);
	setupChildren = true;
	continue;
	}
	else
	break;
	}
	else if(_reconOpt>=2) {
	// cut-off masses for the branching
	const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
	// compute the masses of the children
	Energy masses[3];
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].kinematics) {
	const vector<Energy> & vm = fc[ix].sudakov->virtualMasses(fc[ix].ids);
	Energy2 q2 =
	fc[ix].kinematics->z()(1.-fc[ix].kinematics->z())sqr(fc[ix].kinematics->scale());
	if(fc[ix].ids[0]!=ParticleID::g) q2 += sqr(vm[0]);
	masses[ix+1] = sqrt(q2);
	}
	else {
	masses[ix+1] = virtualMasses[ix+1];
	}
	}
	masses[0] = fb.ids[0]!=ParticleID::g ? virtualMasses[0] : ZERO;
	double z = fb.kinematics->z();
	Energy2 pt2 = z(1.-z)(z(1.-z)sqr(fb.kinematics->scale()) + sqr(masses[0]))
	- sqr(masses[1])(1.-z) - sqr(masses[2])z;
	if(pt2>=ZERO) {
	break;
	}
	// if only the hard emission have to accept it
	else if ((fc[0].hard && !fc[1].kinematics) \|\|
	(fc[1].hard && !fc[0].kinematics) ) {
	break;
	}
	else {
	// reset the scales for the children
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].hard) continue;
	if(fc[ix].kinematics && ! fc[ix].hard )
	children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
	else
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
	children[ix]->virtualMass(ZERO);
	}
	}
	}
	};
	if(_reconOpt>=2) {
	// shower the first particle
	if(fc[0].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 1)
	truncatedTimeLikeShower(children[0],branch,type,fc[0],false);
	// hard emission and subsquent hard emissions
	else if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	// shower the second particle
	if(fc[1].kinematics) {
	// the parent has truncated emission and following line
	if(!fb.hard && fb.iout == 2)
	truncatedTimeLikeShower(children[1],branch,type,fc[1],false);
	// hard emission and subsquent hard emissions
	else if(fb.hard && !branch->children()[1]->children().empty() )
	truncatedTimeLikeShower(children[1],branch->children()[1],type,fc[1],false);
	else
	timeLikeShower(children[1],type,fc[1],false);
	}
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	// branching has happened
	particle->showerKinematics()->updateParent(particle, children,fb.type);
	}
	if(first&&!children.empty())
	particle->showerKinematics()->resetChildren(particle,children);
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}

	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
	ShowerInteraction::Type type2 = convertInteraction(bb.type);
	if(type2==branch->sudakov()->interactionType() &&
	zsplit < 0.5) {
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	continue;
	}
	// others
	if( part[0]->id() != particle->id() \|\| // if particle changes type
	bb.kinematics->pT() > progenitor()->maximumpT(type2) \|\| // pt veto
	bb.kinematics->scale() < branch->scale()) { // angular ordering veto
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	continue;
	}
	// and those from the base class
	if(spaceLikeVetoed(bb,particle)) {
	particle->vetoEmission(bb.type,bb.kinematics->scale());
	continue;
	}
	break;
	}
	}
	if( !bb.kinematics ) {
	//do the hard emission
	ShoKinPtr kinematics =
	branch->sudakov()->createInitialStateBranching( branch->scale(), z, branch->phi(),
	branch->children()[0]->pT() );
	// assign the splitting function and shower kinematics
	particle->showerKinematics( kinematics );
	if(kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent =
	new_ptr( ShowerParticle( branch->branchingParticle()->dataPtr(), false ) );
	ShowerParticlePtr otherChild =
	new_ptr( ShowerParticle( timelike->branchingParticle()->dataPtr(),
	true, true ) );
	ShowerParticleVector theChildren;
	theChildren.push_back( particle );
	theChildren.push_back( otherChild );
	particle->showerKinematics()->
	updateParent( newParent, theChildren, branch->type());
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
	currentTree()->addInitialStateBranching( particle, newParent, otherChild );
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	bool emitted=false;
	if(!hardOnly()) {
	if( branch->parent() ) {
	emitted = truncatedSpaceLikeShower( newParent, beam, branch->parent() , type);
	}
	else {
	emitted = spaceLikeShower( newParent, beam , type);
	}
	}
	if( !emitted ) {
	if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
	kinematics->updateLast( newParent, ZERO, ZERO );
	}
	else {
	pair<Energy,double> kt = intrinsicpT()[progenitor()];
	kinematics->updateLast( newParent,
	kt.first*cos( kt.second ),
	kt.first*sin( kt.second ) );
	}
	}
	particle->showerKinematics()->
	updateChildren( newParent, theChildren,bb.type,false);
	if(hardOnly()) return true;
	// perform the shower of the final-state particle
	if( timelike->children().empty() ) {
	timeLikeShower( otherChild , type,Branching(),true);
	}
	else {
	truncatedTimeLikeShower( otherChild, timelike , type,Branching(), true);
	}
	updateHistory(otherChild);
	// return the emitted
	return true;
	}
	// assign the splitting function and shower kinematics
	particle->showerKinematics( bb.kinematics );
	if(bb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(bb.kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent = new_ptr( ShowerParticle( part[0], false ) );
	ShowerParticlePtr otherChild = new_ptr( ShowerParticle( part[1], true, true ) );
	ShowerParticleVector theChildren;
	theChildren.push_back( particle );
	theChildren.push_back( otherChild );
	particle->showerKinematics()->
	updateParent( newParent, theChildren, bb.type);
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
	currentTree()->addInitialStateBranching( particle, newParent, otherChild );
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	bool emitted = truncatedSpaceLikeShower( newParent, beam, branch,type);
	// now reconstruct the momentum
	if( !emitted ) {
	if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
	bb.kinematics->updateLast( newParent, ZERO, ZERO );
	}
	else {
	pair<Energy,double> kt = intrinsicpT()[ progenitor() ];
	bb.kinematics->updateLast( newParent,
	kt.first*cos( kt.second ),
	kt.first*sin( kt.second ) );
	}
	}
	particle->showerKinematics()->
	updateChildren( newParent, theChildren, bb.type,false);
	// perform the shower of the final-state particle
	timeLikeShower( otherChild , type,Branching(),true);
	updateHistory(otherChild);
	// return the emitted
	return true;
	}

	bool Evolver::
	truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass, HardBranchingPtr branch,
	ShowerInteraction::Type type, Branching fb) {
	// select a branching if we don't have one
	if(!fb.kinematics)
	fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,branch);
	// must be an emission, the forced one it not a truncated one
	assert(fb.kinematics);
	ShowerParticleVector children;
	int ntry=0;
	Branching fc[2];
	bool setupChildren = true;
	while (ntry<50) {
	if(!fc[0].hard) fc[0] = Branching();
	if(!fc[1].hard) fc[1] = Branching();
	++ntry;
	if(setupChildren) {
	++_nFSR;
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// create the ShowerParticle objects for the two children
	children = createTimeLikeChildren(particle,fb.ids);
	// updateChildren the children
	particle->showerKinematics()->
	updateChildren(particle, children, fb.type,_reconOpt>=3);
	setupChildren = false;
	}
	// select branchings for children
	if(!fc[0].kinematics) {
	if(children[0]->id()==particle->id()) {
	// select branching for first particle
	if(!fb.hard)
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,branch);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,
	branch->children()[0]);
	else
	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,
	HardBranchingPtr());
	}
	else {
	// select branching for first particle
	if(fb.hard && !branch->children()[0]->children().empty() )
	fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]);
	else
	fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
	}
	}
	// select branching for the second particle
	if(!fc[1].kinematics) {
	if(children[1]->id()==particle->id()) {
	// select branching for first particle
	if(!fb.hard)
	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,branch);
	else if(fb.hard && ! branch->children()[1]->children().empty() )
	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,
	branch->children()[1]);
	else
	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,
	HardBranchingPtr());
	}
	else {
	if(fb.hard && !branch->children()[1]->children().empty() )
	fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]);
	else
	fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
	}
	}
	// old default
	if(_reconOpt==0 \|\| (_reconOpt==1 && fb.hard) ) {
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),children[0]);
	currentTree()->addInitialStateBranching(particle,children[0],children[1]);
	// shower the first particle
	if(fc[0].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[0]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[0]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[0]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	}
	// shower the second particle
	if(fc[1].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[1]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[1]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[1]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[1],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[1],false);
	}
	}
	updateHistory(children[1]);
	// branching has happened
	break;
	}
	// H7 default
	else if(_reconOpt==1) {
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),children[0]);
	currentTree()->addInitialStateBranching(particle,children[0],children[1]);
	// shower the first particle
	if(fc[0].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[0]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[0]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[0]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	}
	// shower the second particle
	if(fc[1].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[1]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[1]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[1]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[1],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[1],false);
	}
	}
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<children.size();++ix)
	particle->abandonChild(children[ix]);
	children.clear();
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	// generate the new emission
	fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,branch);
	// must be at least hard emission
	assert(fb.kinematics);
	setupChildren = true;
	continue;
	}
	else {
	updateHistory(children[1]);
	break;
	}
	}
	else if(_reconOpt>=2) {
	// cut-off masses for the branching
	const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
	// compute the masses of the children
	Energy masses[3];
	// space-like children
	masses[1] = children[0]->virtualMass();
	// time-like child
	if(fc[1].kinematics) {
	const vector<Energy> & vm = fc[1].sudakov->virtualMasses(fc[1].ids);
	Energy2 q2 =
	fc[1].kinematics->z()(1.-fc[1].kinematics->z())sqr(fc[1].kinematics->scale());
	if(fc[1].ids[0]!=ParticleID::g) q2 += sqr(vm[0]);
	masses[2] = sqrt(q2);
	}
	else {
	masses[2] = virtualMasses[2];
	}
	masses[0]=particle->virtualMass();
	double z = fb.kinematics->z();
	Energy2 pt2 = (1.-z)(zsqr(masses[0])-sqr(masses[1])-z/(1.-z)*sqr(masses[2]));
	if(pt2>=ZERO) {
	break;
	}
	else {
	// reset the scales for the children
	for(unsigned int ix=0;ix<2;++ix) {
	if(fc[ix].kinematics)
	children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
	else {
	if(ix==0)
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,Constants::MaxEnergy);
	else
	children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
	}
	}
	children[0]->virtualMass(_progenitor->progenitor()->mass());
	children[1]->virtualMass(ZERO);
	}
	}
	};
	if(_reconOpt>=2) {
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),children[0]);
	currentTree()->addInitialStateBranching(particle,children[0],children[1]);
	// shower the first particle
	if(fc[0].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[0]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[0]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[0]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[0],false);
	}
	}
	// shower the second particle
	if(fc[1].kinematics) {
	if(children[0]->id()==particle->id()) {
	if(!fb.hard)
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch,type,fc[1]);
	else if(fb.hard && ! branch->children()[0]->children().empty() )
	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
	branch->children()[0],type,fc[1]);
	else
	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[1]);
	}
	else {
	if(fb.hard && !branch->children()[0]->children().empty() )
	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[1],false);
	// normal shower
	else
	timeLikeShower(children[0],type,fc[1],false);
	}
	}
	updateHistory(children[1]);
	}
	return true;
	}

	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
	map<ShowerProgenitorPtr,HardBranchingPtr> partners;
	for(set<HardBranchingPtr>::const_iterator bit=hardTree->branchings().begin();
	bit!=hardTree->branchings().end();++bit) {
	Energy2 dmin( 1e30*GeV2 );
	ShowerProgenitorPtr partner;
	for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
	pit!=particlesToShower.end();++pit) {
	if(partners.find(*pit)!=partners.end()) continue;
	if( (bit).branchingParticle()->id() != (pit).progenitor()->id() ) continue;
	if( (**bit).branchingParticle()->isFinalState() !=
	(**pit).progenitor()->isFinalState() ) continue;
	if( (**pit).progenitor()->isFinalState() ) {
	Energy2 dtest =
	sqr( (pit).progenitor()->momentum().x() - (bit).showerMomentum().x() ) +
	sqr( (pit).progenitor()->momentum().y() - (bit).showerMomentum().y() ) +
	sqr( (pit).progenitor()->momentum().z() - (bit).showerMomentum().z() ) +
	sqr( (pit).progenitor()->momentum().t() - (bit).showerMomentum().t() );
	// add mass difference for identical particles (e.g. Z0 Z0 production)
	dtest += 1e10sqr((pit).progenitor()->momentum().m()-(*bit).showerMomentum().m());
	if( dtest < dmin ) {
	partner = *pit;
	dmin = dtest;
	}
	}
	else {
	// ensure directions are right
	if((pit).progenitor()->momentum().z()/(bit).showerMomentum().z()>ZERO) {
	partner = *pit;
	break;
	}
	}
	}
	if(!partner) throw Exception() << "Failed to match shower and hard trees in Evolver::hardestEmission"
	<< Exception::eventerror;
	partners[partner] = *bit;
	}
	for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
	pit!=particlesToShower.end();++pit) {
	HardBranchingPtr partner = partners[*pit];
	if((**pit).progenitor()->dataPtr()->stable()) {
	(**pit).progenitor()->set5Momentum(partner->showerMomentum());
	(**pit).copy()->set5Momentum(partner->showerMomentum());
	}
	else {
	Lorentz5Momentum oldMomentum = (**pit).progenitor()->momentum();
	Lorentz5Momentum newMomentum = partner->showerMomentum();
	LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
	(**pit).progenitor()->transform(boost);
	(**pit).copy() ->transform(boost);
	boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
	(**pit).progenitor()->transform(boost);
	(**pit).copy() ->transform(boost);
	}
	}
	// correction boosts for daughter trees
	for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	tit = showerTree->treelinks().begin();
	tit != showerTree->treelinks().end();++tit) {
	ShowerTreePtr decayTree = tit->first;
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit = decayTree->incomingLines().begin();
	// reset the momentum of the decay particle
	Lorentz5Momentum oldMomentum = cit->first->progenitor()->momentum();
	Lorentz5Momentum newMomentum = tit->second.second->momentum();
	LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
	decayTree->transform(boost,true);
	boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
	decayTree->transform(boost,true);
	}
	}

	void Evolver::doShowering(bool hard,XCPtr xcomb) {
	// zero number of emissions
	_nis = _nfs = 0;
	// if MC@NLO H event and limited emissions
	// indicate both final and initial state emission
	if ( isMCatNLOHEvent && _limitEmissions != 0 ) {
	_nis = _nfs = 1;
	}
	// extract particles to shower
	vector<ShowerProgenitorPtr> particlesToShower(setupShower(hard));
	// setup the maximum scales for the shower
	if (hardVetoOn()) setupMaximumScales(particlesToShower,xcomb);
	// set the hard scales for the profiles
	setupHardScales(particlesToShower,xcomb);
	// specific stuff for hard processes and decays
	Energy minmass(ZERO), mIn(ZERO);
	// hard process generate the intrinsic p_T once and for all
	if(hard) {
	generateIntrinsicpT(particlesToShower);
	}
	// decay compute the minimum mass of the final-state
	else {
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(particlesToShower[ix]->progenitor()->isFinalState()) {
	if(particlesToShower[ix]->progenitor()->dataPtr()->stable())
	minmass += particlesToShower[ix]->progenitor()->dataPtr()->constituentMass();
	else
	minmass += particlesToShower[ix]->progenitor()->mass();
	}
	else {
	mIn = particlesToShower[ix]->progenitor()->mass();
	}
	}
	// throw exception if decay can't happen
	if ( minmass > mIn ) {
	throw Exception() << "Evolver.cc: Mass of decaying particle is "
	<< "below constituent masses of decay products."
	<< Exception::eventerror;
	}
	}
	// create random particle vector
	vector<ShowerProgenitorPtr> tmp;
	unsigned int nColouredIncoming = 0;
	while(particlesToShower.size()>0){
	unsigned int xx=UseRandom::irnd(particlesToShower.size());
	tmp.push_back(particlesToShower[xx]);
	particlesToShower.erase(particlesToShower.begin()+xx);
	}
	particlesToShower=tmp;
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(!particlesToShower[ix]->progenitor()->isFinalState() &&
	particlesToShower[ix]->progenitor()->coloured()) ++nColouredIncoming;
	}
	bool switchRecon = hard && nColouredIncoming !=1;
	// main shower loop
	unsigned int ntry(0);
	bool reconstructed = false;
	do {
	// clear results of last attempt if needed
	if(ntry!=0) {
	currentTree()->clear();
	setEvolutionPartners(hard,interaction_,true);
	_nis = _nfs = 0;
	// if MC@NLO H event and limited emissions
	// indicate both final and initial state emission
	if ( isMCatNLOHEvent && _limitEmissions != 0 ) {
	_nis = _nfs = 1;
	}
	for(unsigned int ix=0; ix<particlesToShower.size();++ix) {
	SpinPtr spin = particlesToShower[ix]->progenitor()->spinInfo();
	if(spin && spin->decayVertex() &&
	dynamic_ptr_cast<tcSVertexPtr>(spin->decayVertex())) {
	spin->decayVertex(VertexPtr());
	}
	}
	}
	// loop over particles
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	// extract the progenitor
	progenitor(particlesToShower[ix]);
	// final-state radiation
	if(progenitor()->progenitor()->isFinalState()) {
	if(!isFSRadiationON()) continue;
	// perform shower
	progenitor()->hasEmitted(startTimeLikeShower(interaction_));
	}
	// 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,
	interaction_));
	}
	// decay
	else {
	// skip colour and electrically neutral particles
	if(!progenitor()->progenitor()->dataPtr()->coloured() &&
	!progenitor()->progenitor()->dataPtr()->charged()) {
	progenitor()->hasEmitted(false);
	continue;
	}
	// perform shower
	// set the scales correctly. The current scale is the maximum scale for
	// emission not the starting scale
	ShowerParticle::EvolutionScales maxScales(progenitor()->progenitor()->scales());
	progenitor()->progenitor()->scales() = ShowerParticle::EvolutionScales();
	if(progenitor()->progenitor()->dataPtr()->charged()) {
	progenitor()->progenitor()->scales().QED = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QED_noAO = progenitor()->progenitor()->mass();
	}
	if(progenitor()->progenitor()->hasColour()) {
	progenitor()->progenitor()->scales().QCD_c = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QCD_c_noAO = progenitor()->progenitor()->mass();
	}
	if(progenitor()->progenitor()->hasAntiColour()) {
	progenitor()->progenitor()->scales().QCD_ac = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QCD_ac_noAO = progenitor()->progenitor()->mass();
	}
	// perform the shower
	progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass,
	interaction_));
	}
	}
	}
	// do the kinematic reconstruction, checking if it worked
	reconstructed = hard ?
	showerModel()->kinematicsReconstructor()->
	reconstructHardJets (currentTree(),intrinsicpT(),interaction_,
	switchRecon && ntry>maximumTries()/2) :
	showerModel()->kinematicsReconstructor()->
	reconstructDecayJets(currentTree(),interaction_);
	}
	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);
	hardTree(HardTreePtr());
	}

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

	Branching Evolver::selectTimeLikeBranching(tShowerParticlePtr particle,
	ShowerInteraction::Type type,
	HardBranchingPtr branch) {
	Branching fb;
	unsigned int iout=0;
	tcPDPtr pdata[2];
	while (true) {
	// break if doing truncated shower and no truncated shower needed
	if(branch && (!isTruncatedShowerON()\|\|hardOnly())) break;
	fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
	// no emission break
	if(!fb.kinematics) break;
	// special for truncated shower
	if(branch) {
	// check haven't evolved too far
	if(fb.kinematics->scale() < branch->scale()) {
	fb=Branching();
	break;
	}
	// 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
	ShowerInteraction::Type type2 = convertInteraction(fb.type);
	// and evolution
	if(type2==branch->sudakov()->interactionType()) {
	if(zsplit < 0.5 \|\| // hardest line veto
	fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// pt veto
	if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// standard vetos for all emissions
	if(timeLikeVetoed(fb,particle)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	continue;
	}
	break;
	}
	// normal case
	if(!branch) {
	if(fb.kinematics) fb.hard = false;
	return fb;
	}
	// truncated emission
	if(fb.kinematics) {
	fb.hard = false;
	fb.iout = iout;
	return fb;
	}
	// otherwise need to return the hard emission
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin=
	branch->sudakov()->createFinalStateBranching(branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT());
	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() );
	fb.hard = true;
	fb.iout=0;
	// return it
	return fb;
	}

	Branching Evolver::selectSpaceLikeDecayBranching(tShowerParticlePtr particle,
	const ShowerParticle::EvolutionScales & maxScales,
	Energy minmass,ShowerInteraction::Type type,
	HardBranchingPtr branch) {
	Branching fb;
	unsigned int iout=0;
	tcPDPtr pdata[2];
	while (true) {
	// break if doing truncated shower and no truncated shower needed
	if(branch && (!isTruncatedShowerON()\|\|hardOnly())) break;
	// select branching
	fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass,
	_initialenhance,type);
	// return if no radiation
	if(!fb.kinematics) break;
	// special for truncated shower
	if(branch) {
	// check haven't evolved too far
	if(fb.kinematics->scale() < branch->scale()) {
	fb=Branching();
	break;
	}
	// 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;
	}
	ShowerInteraction::Type type2 = convertInteraction(fb.type);
	double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
	if(type2==branch->sudakov()->interactionType()) {
	if(zsplit < 0.5 \|\| // hardest line veto
	fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// pt veto
	if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	}
	// if not vetoed break
	if(spaceLikeDecayVetoed(fb,particle)) {
	// otherwise reset scale and continue
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	break;
	}
	// normal case
	if(!branch) {
	if(fb.kinematics) fb.hard = false;
	return fb;
	}
	// truncated emission
	if(fb.kinematics) {
	fb.hard = false;
	fb.iout = iout;
	return fb;
	}
	// otherwise need to return the hard emission
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin=
	branch->sudakov()->createDecayBranching(branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT());
	IdList idlist(3);
	idlist[0] = particle->id();
	idlist[1] = branch->children()[0]->branchingParticle()->id();
	idlist[2] = branch->children()[1]->branchingParticle()->id();
	// create the branching
	fb = Branching( showerKin, idlist, branch->sudakov(),ShowerPartnerType::QCDColourLine );
	fb.hard=true;
	fb.iout=0;
	// return it
	return fb;
	}
	diff --git a/Shower/Base/PartnerFinder.cc b/Shower/Base/PartnerFinder.cc
	--- a/Shower/Base/PartnerFinder.cc
	+++ b/Shower/Base/PartnerFinder.cc
	@@ -1,528 +1,623 @@
	// -- C++ --
	//
	// PartnerFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the PartnerFinder class.
	//

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

	using namespace Herwig;

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

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

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

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

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

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

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

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

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

	void PartnerFinder::Init() {

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

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

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

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

	}

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

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

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

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

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

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

	-vector< pair<double, tShowerParticlePtr> >
	+vector< pair<double, tShowerParticlePtr> >
	PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
	const ShowerParticleVector &particles,
	const bool isDecayCase) {
	vector< pair<double, tShowerParticlePtr> > partners;
	ShowerParticleVector::const_iterator cjt;
	for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
	if(!(cjt)->data().charged() \|\| particle == cjt) continue;
	double charge = double(particle->data().iCharge()(cjt)->data().iCharge());
	if( FS(particle) != FS(cjt) ) charge =-1.;
	if( QEDPartner_ != 0 && !isDecayCase ) {
	// only include II and FF as requested
	if( QEDPartner_ == 1 && FS(particle) != FS(*cjt) )
	continue;
	// only include IF is requested
	else if(QEDPartner_ == 2 && FS(particle) == FS(*cjt) )
	continue;
	}
	// only keep positive dipoles
	if(charge<0.) partners.push_back(make_pair(-charge,*cjt));
	}
	return partners;
	}
	+
	+void PartnerFinder::setInitialEWEvolutionScales(const ShowerParticleVector &particles,
	+ const bool isDecayCase,
	+ const bool setPartners) {
	+ // loop over all the particles
	+ for(ShowerParticleVector::const_iterator cit = particles.begin();
	+ cit != particles.end(); ++cit) {
	+ // if not weakly interacting continue
	+ if( !weaklyInteracting( (**cit).dataPtr()))
	+ continue;
	+ // find the potential partners
	+ vector<pair<double,tShowerParticlePtr> > partners = findEWPartners(*cit,particles,isDecayCase);
	+ if(partners.empty()) {
	+ throw Exception() << "Failed to find partner in "
	+ << "PartnerFinder::setEWInitialEvolutionScales"
	+ << (**cit) << Exception::eventerror;
	+ }
	+ // calculate the probabilities
	+ double prob(0.);
	+ for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
	+ // normalise
	+ for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
	+ // set the partner if required
	+ int position(-1);
	+ // use QCD partner if set
	+ if(!setPartners&&(*cit)->partner()) {
	+ for(unsigned int ix=0;ix<partners.size();++ix) {
	+ if((*cit)->partner()==partners[ix].second) {
	+ position = ix;
	+ break;
	+ }
	+ }
	+ }
	+ // set the partner
	+ if(setPartners\|\|!(*cit)->partner()\|\|position<0) {
	+ prob = UseRandom::rnd();
	+ for(unsigned int ix=0;ix<partners.size();++ix) {
	+ if(partners[ix].first>prob) {
	+ position = ix;
	+ break;
	+ }
	+ prob -= partners[ix].first;
	+ }
	+ if(position>=0&&(setPartners\|\|!(*cit)->partner())) {
	+ (*cit)->partner(partners[position].second);
	+ }
	+ }
	+ // must have a partner
	+ if(position<0) throw Exception() << "Failed to find partner in "
	+ << "PartnerFinder::setEWInitialEvolutionScales"
	+ << (**cit) << Exception::eventerror;
	+ // Calculate the evolution scales for all possible pairs of of particles
	+ vector<pair<Energy,Energy> > scales;
	+ for(unsigned int ix=0;ix< partners.size();++ix) {
	+ scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
	+ isDecayCase));
	+ }
	+ // store all the possible partners
	+ for(unsigned int ix=0;ix<partners.size();++ix) {
	+ (**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
	+ partners[ix].first,
	+ ShowerPartnerType::EW,
	+ scales[ix].first));
	+ }
	+ // set scales
	+ (**cit).scales().EW = scales[position].first;
	+ (**cit).scales().EW_noAO = scales[position].first;
	+ }
	+}
	+
	+vector< pair<double, tShowerParticlePtr> >
	+PartnerFinder::findEWPartners(tShowerParticlePtr particle,
	+ const ShowerParticleVector &particles,
	+ const bool ) {
	+ vector< pair<double, tShowerParticlePtr> > partners;
	+ ShowerParticleVector::const_iterator cjt;
	+ for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
	+ if(!weaklyInteracting((*cjt)->dataPtr()) \|\|
	+ particle == *cjt) continue;
	+ double charge = 1.;
	+ // only keep positive dipoles
	+ if(charge>0.) partners.push_back(make_pair(charge,*cjt));
	+ }
	+ return partners;
	+}
	diff --git a/Shower/Base/PartnerFinder.h b/Shower/Base/PartnerFinder.h
	--- a/Shower/Base/PartnerFinder.h
	+++ b/Shower/Base/PartnerFinder.h
	@@ -1,209 +1,236 @@
	// -- C++ --
	//
	// PartnerFinder.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_PartnerFinder_H
	#define HERWIG_PartnerFinder_H
	//
	// This is the declaration of the PartnerFinder class.
	//
	#include "Herwig/Shower/ShowerConfig.h"
	+#include "ThePEG/PDT/ParticleData.h"
	#include "ThePEG/Interface/Interfaced.h"
	#include "Evolver.fh"
	#include "PartnerFinder.fh"

	namespace Herwig {

	using namespace ThePEG;

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

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

	public:

	/**
	* The default constructor.
	*/
	PartnerFinder() : partnerMethod_(0), QEDPartner_(0), scaleChoice_(0) {}

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

	public:

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

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

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

	protected:

	/**
	* Members to set the scales for different interactions
	*/
	//@{
	/**
	* Set initial scales for a QCD interaction
	*/
	virtual void setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
	const bool isDecayCase,
	const bool setPartners=true);

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

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

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

	/**
	+ * Find the EW partners
	+ * @param particle The particle to find the partners for
	+ * @param particles The full set of particles to search
	+ */
	+ vector< pair<double, tShowerParticlePtr> >
	+ findEWPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles,
	+ const bool isDecayCase);
	+
	+ /**
	* Given a pair of particles, supposedly partners w.r.t. an interaction,
	* this method returns their initial evolution scales as a pair.
	* If something wrong happens, it returns the null (ZERO,ZERO) pair.
	* This method is used by the above setXXXInitialEvolutionScales
	* methods.
	* These methods must be overiden in inheriting classes
	*/
	//@{
	/**
	* General method to calculate the initial evolution scales
	*/
	virtual pair<Energy,Energy> calculateInitialEvolutionScales(const ShowerPPair &,
	const bool isDecayCase);

	/**
	* Calculate the initial evolution scales for two final-state particles
	*/
	virtual pair<Energy,Energy> calculateFinalFinalScales(const ShowerPPair &)=0;

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

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

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

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

	private:

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

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

	/**
	* Choice of the scale
	*/
	int scaleChoice_;
	};

	}

	#endif /* HERWIG_PartnerFinder_H */
	diff --git a/Shower/Base/ShowerParticle.h b/Shower/Base/ShowerParticle.h
	--- a/Shower/Base/ShowerParticle.h
	+++ b/Shower/Base/ShowerParticle.h
	@@ -1,536 +1,541 @@
	// -- C++ --
	//
	// ShowerParticle.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_ShowerParticle_H
	#define HERWIG_ShowerParticle_H
	//
	// This is the declaration of the ShowerParticle class.
	//

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

	namespace Herwig {

	using namespace ThePEG;

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

	public:

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

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

	/**
	* The partner
	*/
	tShowerParticlePtr partner;

	/**
	* Weight
	*/
	double weight;

	/**
	* Type
	*/
	ShowerPartnerType::Type type;

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

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

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

	/**
	* QED scale
	*/
	Energy QED;

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

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

	/**
	* QED scale
	*/
	Energy QED_noAO;

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

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

	/**
	* EW scale
	*/
	Energy EW;

	/**
	+ * EW scale
	+ */
	+ Energy EW_noAO;
	+
	+ /**
	* Maximum allowed virtuality of the particle
	*/
	Energy2 Max_Q2;
	};


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

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

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

	public:

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

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

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

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

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

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

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

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

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

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

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


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

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


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

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

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

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

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

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

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

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

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

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

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

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

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


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


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

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

	public:

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

	protected:

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


	protected:

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

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

	private:

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

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

	private:

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

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

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

	/**
	* Dimensionless parameters
	*/
	Parameters _parameters;

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

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

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

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

	/**
	* Virtual mass
	*/
	Energy _vMass;

	/**
	* Partners
	*/
	tShowerParticlePtr _partner;

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

	/**
	* Progenitor
	*/
	tShowerParticlePtr _progenitor;

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

	};

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

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	/** @endcond */

	}

	#endif /* HERWIG_ShowerParticle_H */
	diff --git a/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc b/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc
	--- a/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc
	+++ b/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc
	@@ -1,216 +1,214 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the HalfHalfOneEWSplitFn class.
	//

	#include "HalfHalfOneEWSplitFn.h"
	#include "ThePEG/StandardModel/StandardModelBase.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"

	using namespace Herwig;

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

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

	void HalfHalfOneEWSplitFn::persistentOutput(PersistentOStream & os) const {
	os << gZ_ << gWL_;
	}

	void HalfHalfOneEWSplitFn::persistentInput(PersistentIStream & is, int) {
	is >> gZ_ >> gWL_;
	}

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

	void HalfHalfOneEWSplitFn::Init() {

	static ClassDocumentation<HalfHalfOneEWSplitFn> documentation
	("The HalfHalfOneEWSplitFn class implements the splitting q->qWand q->qZ");

	}

	void HalfHalfOneEWSplitFn::doinit() {
	- cerr << "testing in do init\n";
	SplittingFunction::doinit();
	- cerr << "testing in do init\n";
	tcSMPtr sm = generator()->standardModel();
	double sw2 = sm->sin2ThetaW();
	// left-handled W coupling
	gWL_ = 1./sqrt(2.*sw2);
	// Z couplings
	double fact = 0.25/sqrt(sw2*(1.-sw2));
	for(int ix=1;ix<4;++ix) {
	gZ_[2ix-1] = make_pair(fact(sm->vd() + sm->ad()),
	fact*(sm->vd() - sm->ad() ));
	gZ_[2ix ] = make_pair(fact(sm->vu() + sm->au() ),
	fact*(sm->vu() - sm->au() ));
	gZ_[2ix+9 ] = make_pair(fact(sm->ve() + sm->ae() ),
	fact*(sm->ve() - sm->ae() ));
	gZ_[2ix+10] = make_pair(fact(sm->vnu() + sm->anu()),
	fact*(sm->vnu() - sm->anu()));
	}
	}

	void HalfHalfOneEWSplitFn::getCouplings(double & gL, double & gR, const IdList & ids) const {
	if(ids[2]==ParticleID::Z0) {
	map<long,pair<double,double> >::const_iterator it = gZ_.find(abs(ids[0]));
	assert(it!=gZ_.end());
	gL = it->second.first ;
	gR = it->second.second;
	}
	else if(abs(ids[2])==ParticleID::Wplus) {
	gL = gWL_;
	}
	else
	assert(false);
	}

	double HalfHalfOneEWSplitFn::P(const double z, const Energy2 t,
	const IdList &ids, const bool mass) const {
	double gL(0.),gR(0.);
	getCouplings(gL,gR,ids);
	double val = (1. + sqr(z))/(1.-z);
	if(mass) {
	Energy m = getParticleData(ids[2])->mass();
	val -= sqr(m)/t;
	}
	val *= (sqr(gL)+sqr(gR));
	return colourFactor(ids)*val;
	}


	double HalfHalfOneEWSplitFn::overestimateP(const double z,
	const IdList & ids) const {
	double gL(0.),gR(0.);
	getCouplings(gL,gR,ids);
	return 2.sqr(max(gL,gR))colourFactor(ids)/(1.-z);
	}

	double HalfHalfOneEWSplitFn::ratioP(const double z, const Energy2 t,
	const IdList & ids, const bool mass) const {
	double gL(0.),gR(0.);
	getCouplings(gL,gR,ids);
	double val = 1. + sqr(z);
	if(mass) {
	Energy m = getParticleData(ids[2])->mass();
	val -= (1.-z)*sqr(m)/t;
	}
	val = 0.5(sqr(gL)+sqr(gR))/sqr(max(gL,gR));
	return val;
	}

	double HalfHalfOneEWSplitFn::integOverP(const double z,
	const IdList & ids,
	unsigned int PDFfactor) const {
	double gL(0.),gR(0.);
	getCouplings(gL,gR,ids);
	double pre = colourFactor(ids)*sqr(max(gL,gR));
	switch (PDFfactor) {
	case 0:
	return -2.preMath::log1m(z);
	case 1:
	return 2.prelog(z/(1.-z));
	case 2:
	return 2.*pre/(1.-z);
	case 3:
	default:
	throw Exception() << "HalfHalfOneEWSplitFn::integOverP() invalid PDFfactor = "
	<< PDFfactor << Exception::runerror;
	}
	}

	double HalfHalfOneEWSplitFn::invIntegOverP(const double r, const IdList & ids,
	unsigned int PDFfactor) const {
	double gL(0.),gR(0.);
	getCouplings(gL,gR,ids);
	double pre = colourFactor(ids)*sqr(max(gL,gR));
	switch (PDFfactor) {
	case 0:
	return 1. - exp(- 0.5*r/pre);
	case 1:
	return 1./(1.-exp(-0.5*r/pre));
	case 2:
	return 1.-2.*pre/r;
	case 3:
	default:
	throw Exception() << "HalfHalfOneEWSplitFn::invIntegOverP() invalid PDFfactor = "
	<< PDFfactor << Exception::runerror;
	}
	}

	bool HalfHalfOneEWSplitFn::accept(const IdList &ids) const {
	if(ids.size()!=3) return false;
	if(ids[2]==ParticleID::Z0) {
	if(ids[0]==ids[1] &&
	((ids[0]>=1 && ids[0]<=6) \|\| (ids[0]>=11&&ids[0]<=16) )) return true;
	}
	else if(abs(ids[2])==ParticleID::Wplus) {
	if(!((ids[0]>=1 && ids[0]<=6) \|\| (ids[0]>=11&&ids[0]<=16) )) return false;
	if(!((ids[1]>=1 && ids[1]<=6) \|\| (ids[1]>=11&&ids[1]<=16) )) return false;
	if(ids[0]+1!=ids[1] && ids[0]-1!=ids[1]) return false;
	int out = getParticleData(ids[1])->iCharge()+getParticleData(ids[2])->iCharge();
	if(getParticleData(ids[0])->iCharge()==out) return true;
	}
	return false;
	}

	vector<pair<int, Complex> >
	HalfHalfOneEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
	const RhoDMatrix &) {
	// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
	// and rest = splitting function for Tr(rho)=1 as required by defn
	return vector<pair<int, Complex> >(1,make_pair(0,1.));
	}

	vector<pair<int, Complex> >
	HalfHalfOneEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & ,
	const RhoDMatrix &) {
	// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
	// and rest = splitting function for Tr(rho)=1 as required by defn
	return vector<pair<int, Complex> >(1,make_pair(0,1.));
	}

	DecayMEPtr HalfHalfOneEWSplitFn::matrixElement(const double z, const Energy2 t,
	const IdList & ids, const double phi,
	bool) {
	// calculate the kernal
	DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin1)));
	Energy m = getParticleData(ids[2])->mass();
	double gL(0.),gR(0.);
	getCouplings(gL,gR,ids);
	double mt = m/sqrt(t);
	double root = sqrt(1.-sqr(m)/t/(1-z));
	double romz = sqrt(1.-z);
	double rz = sqrt(z);
	double r2 = sqrt(2.);
	Complex phase = exp(Complex(0.,1.)*phi);
	Complex cphase = conj(phase);
	(kernal)(0,0,0) = -phaseroot*gL/romz;
	(kernal)(1,1,2) = cphaseroot*gR/romz;
	(kernal)(0,0,2) = cphasezrootgL/romz;
	(kernal)(1,1,0) = -phasezrootgR/romz;
	// long terms
	(kernal)(1,1,1) =-gRmtr2rz/(1-z);
	(kernal)(0,0,1) =-gLmtr2rz/(1-z);
	// +- -+ terms zero due quark mass
	for(unsigned int ix=0;ix<3;++ix) {
	(*kernal)(1,0,ix) = 0.;
	(*kernal)(0,1,ix) = 0.;
	}
	// return the answer
	return kernal;
	}
	diff --git a/Shower/SplittingFunctions/SplittingFunction.cc b/Shower/SplittingFunctions/SplittingFunction.cc
	--- a/Shower/SplittingFunctions/SplittingFunction.cc
	+++ b/Shower/SplittingFunctions/SplittingFunction.cc
	@@ -1,940 +1,968 @@
	// -- C++ --
	//
	// SplittingFunction.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the SplittingFunction class.
	//

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

	using namespace Herwig;

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

	void SplittingFunction::Init() {

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

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

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

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

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

	}

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

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

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

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

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

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

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

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

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

	namespace {

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

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

	}

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

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

	void SplittingFunction::evaluateDecayScales(ShowerPartnerType::Type partnerType,
	Energy scale, double z,
	tShowerParticlePtr parent,
	tShowerParticlePtr spacelike,
	tShowerParticlePtr timelike) {
	assert(parent->id()==spacelike->id());
	// angular-ordered scale for 2nd child
	Energy AOScale = (1.-z)*scale;
	// QED
	if(partnerType==ShowerPartnerType::QED) {
	// timelike
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	timelike->scales().QCD_c = ZERO;
	timelike->scales().QCD_c_noAO = ZERO;
	timelike->scales().QCD_ac = ZERO;
	timelike->scales().QCD_ac_noAO = ZERO;
	// spacelike
	spacelike->scales().QED = scale;
	spacelike->scales().QED_noAO = scale;
	}
	// QCD
	else if(partnerType==ShowerPartnerType::QCDColourLine \|\|
	partnerType==ShowerPartnerType::QCDAntiColourLine) {
	// timelike
	timelike->scales().QED = ZERO;
	timelike->scales().QED_noAO = ZERO;
	timelike->scales().QCD_c = AOScale;
	timelike->scales().QCD_c_noAO = scale;
	timelike->scales().QCD_ac = AOScale;
	timelike->scales().QCD_ac_noAO = scale;
	// spacelike
	spacelike->scales().QED = max(scale,parent->scales().QED );
	spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO );
	}
	else
	assert(false);
	spacelike->scales().QCD_c = max(scale,parent->scales().QCD_c );
	spacelike->scales().QCD_c_noAO = max(scale,parent->scales().QCD_c_noAO );
	spacelike->scales().QCD_ac = max(scale,parent->scales().QCD_ac );
	spacelike->scales().QCD_ac_noAO = max(scale,parent->scales().QCD_ac_noAO);
	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 5:52 PM (10 h, 6 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4015963
Default Alt Text: (210 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions