No OneTemporary
Actions

Size

183 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,3225 +1,3227 @@
	// -- 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 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();
	}
	+ progenitor()->progenitor()->scales().EW = progenitor()->progenitor()->mass();
	+ progenitor()->progenitor()->scales().EW_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/ShowerParticle.h b/Shower/Base/ShowerParticle.h
	--- a/Shower/Base/ShowerParticle.h
	+++ b/Shower/Base/ShowerParticle.h
	@@ -1,541 +1,543 @@
	// -- 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
	+ << " EW=" << es.EW / GeV
	<< " QED_noAO=" << es.QED_noAO / GeV
	<< " QCD_c_noAO=" << es.QCD_c_noAO / GeV
	<< " QCD_ac_noAO=" << es.QCD_ac_noAO / GeV
	+ << " EW_noAO=" << es.EW_noAO / GeV
	<< '\n';
	return os;
	}

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	/** @endcond */

	}

	#endif /* HERWIG_ShowerParticle_H */
	diff --git a/Shower/Base/SudakovFormFactor.cc b/Shower/Base/SudakovFormFactor.cc
	--- a/Shower/Base/SudakovFormFactor.cc
	+++ b/Shower/Base/SudakovFormFactor.cc
	@@ -1,326 +1,326 @@
	// -- C++ --
	//
	// SudakovFormFactor.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 SudakovFormFactor class.
	//

	#include "SudakovFormFactor.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ShowerKinematics.h"
	#include "ShowerParticle.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "Herwig/Shower/ShowerHandler.h"

	using namespace Herwig;

	DescribeAbstractClass<SudakovFormFactor,Interfaced>
	describeSudakovFormFactor ("Herwig::SudakovFormFactor","");

	void SudakovFormFactor::persistentOutput(PersistentOStream & os) const {
	os << splittingFn_ << alpha_ << pdfmax_ << particles_ << pdffactor_
	<< a_ << b_ << ounit(c_,GeV) << ounit(kinCutoffScale_,GeV) << cutOffOption_
	<< ounit(vgcut_,GeV) << ounit(vqcut_,GeV)
	<< ounit(pTmin_,GeV) << ounit(pT2min_,GeV2)
	<< theFactorizationScaleFactor << theRenormalizationScaleFactor;
	}

	void SudakovFormFactor::persistentInput(PersistentIStream & is, int) {
	is >> splittingFn_ >> alpha_ >> pdfmax_ >> particles_ >> pdffactor_
	>> a_ >> b_ >> iunit(c_,GeV) >> iunit(kinCutoffScale_,GeV) >> cutOffOption_
	>> iunit(vgcut_,GeV) >> iunit(vqcut_,GeV)
	>> iunit(pTmin_,GeV) >> iunit(pT2min_,GeV2)
	>> theFactorizationScaleFactor >> theRenormalizationScaleFactor;
	}

	void SudakovFormFactor::Init() {

	static ClassDocumentation<SudakovFormFactor> documentation
	("The SudakovFormFactor class is the base class for the implementation of Sudakov"
	" form factors in Herwig");

	static Reference<SudakovFormFactor,SplittingFunction>
	interfaceSplittingFunction("SplittingFunction",
	"A reference to the SplittingFunction object",
	&Herwig::SudakovFormFactor::splittingFn_,
	false, false, true, false);

	static Reference<SudakovFormFactor,ShowerAlpha>
	interfaceAlpha("Alpha",
	"A reference to the Alpha object",
	&Herwig::SudakovFormFactor::alpha_,
	false, false, true, false);

	static Parameter<SudakovFormFactor,double> interfacePDFmax
	("PDFmax",
	"Maximum value of PDF weight. ",
	&SudakovFormFactor::pdfmax_, 35.0, 1.0, 100000.0,
	false, false, Interface::limited);

	static Switch<SudakovFormFactor,unsigned int> interfacePDFFactor
	("PDFFactor",
	"Include additional factors in the overestimate for the PDFs",
	&SudakovFormFactor::pdffactor_, 0, false, false);
	static SwitchOption interfacePDFFactorOff
	(interfacePDFFactor,
	"Off",
	"Don't include any factors",
	0);
	static SwitchOption interfacePDFFactorOverZ
	(interfacePDFFactor,
	"OverZ",
	"Include an additional factor of 1/z",
	1);
	static SwitchOption interfacePDFFactorOverOneMinusZ
	(interfacePDFFactor,
	"OverOneMinusZ",
	"Include an additional factor of 1/(1-z)",
	2);
	static SwitchOption interfacePDFFactorOverZOneMinusZ
	(interfacePDFFactor,
	"OverZOneMinusZ",
	"Include an additional factor of 1/z/(1-z)",
	3);


	static Switch<SudakovFormFactor,unsigned int> interfaceCutOffOption
	("CutOffOption",
	"The type of cut-off to use to end the shower",
	&SudakovFormFactor::cutOffOption_, 0, false, false);
	static SwitchOption interfaceCutOffOptionDefault
	(interfaceCutOffOption,
	"Default",
	"Use the standard Herwig cut-off on virtualities with the minimum"
	" virtuality depending on the mass of the branching particle",
	0);
	static SwitchOption interfaceCutOffOptionFORTRAN
	(interfaceCutOffOption,
	"FORTRAN",
	"Use a FORTRAN-like cut-off on virtualities",
	1);
	static SwitchOption interfaceCutOffOptionpT
	(interfaceCutOffOption,
	"pT",
	"Use a cut on the minimum allowed pT",
	2);

	static Parameter<SudakovFormFactor,double> interfaceaParameter
	("aParameter",
	"The a parameter for the kinematic cut-off",
	&SudakovFormFactor::a_, 0.3, -10.0, 10.0,
	false, false, Interface::limited);

	static Parameter<SudakovFormFactor,double> interfacebParameter
	("bParameter",
	"The b parameter for the kinematic cut-off",
	&SudakovFormFactor::b_, 2.3, -10.0, 10.0,
	false, false, Interface::limited);

	static Parameter<SudakovFormFactor,Energy> interfacecParameter
	("cParameter",
	"The c parameter for the kinematic cut-off",
	&SudakovFormFactor::c_, GeV, 0.3GeV, 0.1GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<SudakovFormFactor,Energy>
	interfaceKinScale ("cutoffKinScale",
	"kinematic cutoff scale for the parton shower phase"
	" space (unit [GeV])",
	&SudakovFormFactor::kinCutoffScale_, GeV,
	2.3GeV, 0.001GeV, 10.0*GeV,false,false,false);

	static Parameter<SudakovFormFactor,Energy> interfaceGluonVirtualityCut
	("GluonVirtualityCut",
	"For the FORTRAN cut-off option the minimum virtuality of the gluon",
	&SudakovFormFactor::vgcut_, GeV, 0.85GeV, 0.1GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<SudakovFormFactor,Energy> interfaceQuarkVirtualityCut
	("QuarkVirtualityCut",
	"For the FORTRAN cut-off option the minimum virtuality added to"
	" the mass for particles other than the gluon",
	&SudakovFormFactor::vqcut_, GeV, 0.85GeV, 0.1GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<SudakovFormFactor,Energy> interfacepTmin
	("pTmin",
	"The minimum pT if using a cut-off on the pT",
	&SudakovFormFactor::pTmin_, GeV, 1.0GeV, ZERO, 10.0GeV,
	false, false, Interface::limited);
	}

	bool SudakovFormFactor::alphaSVeto(Energy2 pt2) const {
	pt2 *= sqr(renormalizationScaleFactor());
	return UseRandom::rnd() > ThePEG::Math::powi(alpha_->ratio(pt2),
	splittingFn_->interactionOrder());
	}

	bool SudakovFormFactor::
	PDFVeto(const Energy2 t, const double x,
	const tcPDPtr parton0, const tcPDPtr parton1,
	Ptr<BeamParticleData>::transient_const_pointer beam) const {
	assert(pdf_);

	Energy2 theScale = t * sqr(factorizationScaleFactor());
	if (theScale < sqr(freeze_)) theScale = sqr(freeze_);

	double newpdf(0.0), oldpdf(0.0);
	//different treatment of MPI ISR is done via CascadeHandler::resetPDFs()
	newpdf=pdf_->xfx(beam,parton0,theScale,x/z());
	oldpdf=pdf_->xfx(beam,parton1,theScale,x);


	if(newpdf<=0.) return true;
	if(oldpdf<=0.) return false;
	double ratio = newpdf/oldpdf;

	double maxpdf(pdfmax_);
	switch (pdffactor_) {
	case 1:
	maxpdf /= z();
	break;
	case 2:
	maxpdf /= 1.-z();
	break;
	case 3:
	maxpdf /= (z()*(1.-z()));
	break;
	}

	// ratio / PDFMax must be a probability <= 1.0
	if (ratio > maxpdf) {
	generator()->log() << "PDFVeto warning: Ratio > " << name()
	<< ":PDFmax (by a factor of "
	<< ratio/maxpdf <<") for "
	<< parton0->PDGName() << " to "
	<< parton1->PDGName() << "\n";
	}
	return ratio < UseRandom::rnd()*maxpdf;
	}

	void SudakovFormFactor::addSplitting(const IdList & in) {
	bool add=true;
	for(unsigned int ix=0;ix<particles_.size();++ix) {
	if(particles_[ix].size()==in.size()) {
	bool match=true;
	for(unsigned int iy=0;iy<in.size();++iy) {
	if(particles_[ix][iy]!=in[iy]) {
	match=false;
	break;
	}
	}
	if(match) {
	add=false;
	break;
	}
	}
	}
	if(add) particles_.push_back(in);
	}

	void SudakovFormFactor::removeSplitting(const IdList & in) {
	for(vector<IdList>::iterator it=particles_.begin();
	it!=particles_.end();++it) {
	if(it->size()==in.size()) {
	bool match=true;
	for(unsigned int iy=0;iy<in.size();++iy) {
	if((*it)[iy]!=in[iy]) {
	match=false;
	break;
	}
	}
	if(match) {
	vector<IdList>::iterator itemp=it;
	--itemp;
	particles_.erase(it);
	it = itemp;
	}
	}
	}
	}

	Energy2 SudakovFormFactor::guesst(Energy2 t1,unsigned int iopt,
	- const IdList &ids,
	- double enhance,bool ident) const {
	+ const IdList &ids,
	+ double enhance,bool ident) const {
	unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_;
	double c =
	1./((splittingFn_->integOverP(zlimits_.second,ids,pdfopt) -
	splittingFn_->integOverP(zlimits_.first ,ids,pdfopt))*
	alpha_->overestimateValue()/Constants::twopi*enhance);
	assert(iopt<=2);
	if(iopt==1) {
	c/=pdfmax_;
	//symmetry of FS gluon splitting
	if(ident) c*=0.5;
	}
	else if(iopt==2) c*=-1.;
	if(splittingFn_->interactionOrder()==1) {
	double r = UseRandom::rnd();
	if(iopt!=2 \|\| c*log(r)<log(Constants::MaxEnergy2/t1)) {
	return t1*pow(r,c);
	}
	else
	return Constants::MaxEnergy2;
	}
	else {
	assert(false && "Units are dubious here.");
	int nm(splittingFn()->interactionOrder()-1);
	c/=Math::powi(alpha_->overestimateValue()/Constants::twopi,nm);
	return t1 / pow (1. - nmclog(UseRandom::rnd())
	* Math::powi(t1*UnitRemoval::InvE2,nm)
	,1./double(nm));
	}
	}

	double SudakovFormFactor::guessz (unsigned int iopt, const IdList &ids) const {
	unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_;
	double lower = splittingFn_->integOverP(zlimits_.first,ids,pdfopt);
	return splittingFn_->invIntegOverP
	(lower + UseRandom::rnd()*(splittingFn_->integOverP(zlimits_.second,ids,pdfopt) -
	lower),ids,pdfopt);
	}

	void SudakovFormFactor::doinit() {
	Interfaced::doinit();
	pT2min_ = cutOffOption()==2 ? sqr(pTmin_) : ZERO;
	}

	const vector<Energy> & SudakovFormFactor::virtualMasses(const IdList & ids) {
	static vector<Energy> output;
	output.clear();
	if(cutOffOption() == 0) {
	for(unsigned int ix=0;ix<ids.size();++ix)
	output.push_back(getParticleData(ids[ix])->mass());
	Energy kinCutoff=
	kinematicCutOff(kinScale(),*std::max_element(output.begin(),output.end()));
	for(unsigned int ix=0;ix<output.size();++ix)
	output[ix]=max(kinCutoff,output[ix]);
	}
	else if(cutOffOption() == 1) {
	for(unsigned int ix=0;ix<ids.size();++ix) {
	output.push_back(getParticleData(ids[ix])->mass());
	output.back() += ids[ix]==ParticleID::g ? vgCut() : vqCut();
	}
	}
	else if(cutOffOption() == 2) {
	for(unsigned int ix=0;ix<ids.size();++ix)
	output.push_back(getParticleData(ids[ix])->mass());
	}
	else {
	throw Exception() << "Unknown option for the cut-off"
	<< " in SudakovFormFactor::virtualMasses()"
	<< Exception::runerror;
	}
	return output;
	}
	diff --git a/Shower/SplittingFunctions/SplittingFunction.cc b/Shower/SplittingFunctions/SplittingFunction.cc
	--- a/Shower/SplittingFunctions/SplittingFunction.cc
	+++ b/Shower/SplittingFunctions/SplittingFunction.cc
	@@ -1,1002 +1,1028 @@
	// -- 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;
	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;
	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().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().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) {
	// 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 if(partnerType==ShowerPartnerType::EW) {
	if(abs(spacelike->id())!=ParticleID::Wplus &&
	spacelike->id() !=ParticleID::Z0 ) {
	// QCD scales
	parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c );
	parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO );
	parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac );
	parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
	timelike->scales().QCD_c = ZERO;
	timelike->scales().QCD_c_noAO = ZERO;
	timelike->scales().QCD_ac = ZERO;
	timelike->scales().QCD_ac_noAO = ZERO;
	// QED scales
	if(timelike->id()==ParticleID::Z0) {
	parent ->scales().QED = min(scale,spacelike->scales().QED );
	parent ->scales().QED_noAO = min(scale,spacelike->scales().QED_noAO );
	timelike->scales().QED = ZERO;
	timelike->scales().QED_noAO = ZERO;
	}
	else {
	parent ->scales().QED = scale;
	parent ->scales().QED_noAO = scale;
	timelike->scales().QED = AOScale;
	timelike->scales().QED_noAO = scale;
	}
	// EW scales
	parent ->scales().EW = scale;
	parent ->scales().EW_noAO = scale;
	timelike->scales().EW = AOScale;
	timelike->scales().EW_noAO = scale;
	}
	else assert(false);
	}
	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;
	+ timelike->scales().EW = ZERO;
	+ timelike->scales().EW_noAO = ZERO;
	// spacelike
	spacelike->scales().QED = scale;
	spacelike->scales().QED_noAO = scale;
	+ spacelike->scales().EW = max(scale,parent->scales().EW );
	+ spacelike->scales().EW_noAO = max(scale,parent->scales().EW_noAO );
	}
	// QCD
	else if(partnerType==ShowerPartnerType::QCDColourLine \|\|
	partnerType==ShowerPartnerType::QCDAntiColourLine) {
	// timelike
	timelike->scales().QED = ZERO;
	timelike->scales().QED_noAO = ZERO;
	timelike->scales().QCD_c = AOScale;
	timelike->scales().QCD_c_noAO = scale;
	timelike->scales().QCD_ac = AOScale;
	timelike->scales().QCD_ac_noAO = scale;
	+ timelike->scales().EW = ZERO;
	+ timelike->scales().EW_noAO = ZERO;
	// spacelike
	spacelike->scales().QED = max(scale,parent->scales().QED );
	spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO );
	+ spacelike->scales().EW = max(scale,parent->scales().EW );
	+ spacelike->scales().EW_noAO = max(scale,parent->scales().EW_noAO );
	}
	- else {
	+ else if(partnerType==ShowerPartnerType::EW) {
	+ // EW
	+ timelike->scales().EW = AOScale;
	+ timelike->scales().EW_noAO = scale;
	+ spacelike->scales().EW = max(scale,parent->scales().EW );
	+ spacelike->scales().EW_noAO = max(scale,parent->scales().EW_noAO );
	+ // QCD
	+ timelike->scales().QCD_c = ZERO;
	+ timelike->scales().QCD_c_noAO = ZERO;
	+ timelike->scales().QCD_ac = ZERO;
	+ timelike->scales().QCD_ac_noAO = ZERO;
	+ timelike->scales().EW = ZERO;
	+ timelike->scales().EW_noAO = ZERO;
	+ // QED
	+ timelike->scales().QED = ZERO;
	+ timelike->scales().QED_noAO = ZERO;
	+ spacelike->scales().QED = max(scale,parent->scales().QED );
	+ spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO );
	+ }
	+ else
	assert(false);
	- }
	spacelike->scales().QCD_c = max(scale,parent->scales().QCD_c );
	spacelike->scales().QCD_c_noAO = max(scale,parent->scales().QCD_c_noAO );
	spacelike->scales().QCD_ac = max(scale,parent->scales().QCD_ac );
	spacelike->scales().QCD_ac_noAO = max(scale,parent->scales().QCD_ac_noAO);
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 3:24 PM (1 d, 18 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3804942
Default Alt Text: (183 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions