No OneTemporary
Actions

Size

254 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,3155 +1,3162 @@
	// -- 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
	<< _massVetoOption << isMCatNLOSEvent << isMCatNLOHEvent
	<< isPowhegSEvent << isPowhegHEvent
	<< theFactorizationScaleFactor << theRenormalizationScaleFactor
	<< interaction_<< interactions_.size();
	for(unsigned int ix=0;ix<interactions_.size();++ix)
	os << oenum(interactions_[ix]);
	}

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

	void Evolver::doinit() {
	Interfaced::doinit();
	// interactions may have been changed through a setup file so we
	// clear it up here
	interactions_.clear();
	if(interaction_==0) {
	interactions_.push_back(ShowerInteraction::QCD);
	interactions_.push_back(ShowerInteraction::QED);
	}
	else if(interaction_==1) {
	interactions_.push_back(ShowerInteraction::QCD);
	}
	else if(interaction_==2) {
	interactions_.push_back(ShowerInteraction::QED);
	interactions_.push_back(ShowerInteraction::QCD);
	}
	else if(interaction_==3) {
	interactions_.push_back(ShowerInteraction::QED);
	}
	else if(interaction_==4) {
	interactions_.push_back(ShowerInteraction::Both);
	}
	}

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

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

	static Switch<Evolver,unsigned int> interfaceMassVetoOption
	("MassVetoOption",
	"Option for the handling of the mass vetos",
	&Evolver::_massVetoOption, 1, false, false);
	static SwitchOption interfaceMassVetoOptionReset
	(interfaceMassVetoOption,
	"Reset",
	"Try another branching without resetting the starting scale",
	0);
	static SwitchOption interfaceMassVetoOptionInclude
	(interfaceMassVetoOption,
	"Include",
	"Include the veto in the scale generation via the veto algorithm",
	1);


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

	}

	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->lastCentralScale()));
	}
	}
	// 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->lastCentralScale() );
	}
	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::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);
	}
	}
	}

	Branching Evolver::selectTimeLikeBranching(tShowerParticlePtr particle,
	ShowerInteraction::Type type) {
	Branching fb;
	while (true) {
	fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
	// no emission return
	if(!fb.kinematics) return fb;
	// if emission OK break
	if(!timeLikeVetoed(fb,particle)) break;
	// otherwise reset scale and continue - SO IS involved in veto algorithm
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	}
	return fb;
	}

	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;
	}
	ShowerParticleVector children;
	int ntry=0;
	while (ntry<50) {
	++ntry;
	// generate the emission
	if(!fb.kinematics)
	fb = selectTimeLikeBranching(particle,type);
	// no emission, return
	if(!fb.kinematics) {
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return false;
	}
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	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);
	// update number of emissions
	++_nfs;
	if(_limitEmissions!=0) {
	if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
	if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}
	// select branchings for children
	Branching fc[2] = {selectTimeLikeBranching(children[0],type),
	selectTimeLikeBranching(children[1],type)};
	// old recon option
	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();
	}
	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());
	if(_massVetoOption==1) particle->vetoEmission(fb.type,fb.kinematics->scale());
	fb = Branching();
	continue;
	}
	}
	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) {
	// branching has happened
	particle->showerKinematics()->
	updateParent(particle, children,fb.type);
	// 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();
	}
	else {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<children.size();++ix)
	particle->abandonChild(children[ix]);
	children.clear();
	if(_massVetoOption==1) particle->vetoEmission(fb.type,fb.kinematics->scale());
	if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
	fb = Branching();
	continue;
	}
	}
	break;
	};
	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);
	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;
	while (true) {
	fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass,
	_initialenhance,type);
	// return if no radiation
	if(!fb.kinematics) return false;
	// if not vetoed break
	if(!spaceLikeDecayVetoed(fb,particle)) break;
	// otherwise reset scale and continue
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	}
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// Create the ShowerParticle objects for the two children of
	// the emitting particle; set the parent/child relationship
	// if same as definition create particles, otherwise create cc
	tcPDPtr pdata[2];
	for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
	if(particle->id()!=fb.ids[0]) {
	for(unsigned int ix=0;ix<2;++ix) {
	tPDPtr cc(pdata[ix]->CC());
	if(cc) pdata[ix]=cc;
	}
	}
	ShowerParticleVector theChildren;
	for(unsigned int ix=0;ix<2;++ix) {
	theChildren.push_back(new_ptr(ShowerParticle(pdata[ix],true)));
	if(theChildren[ix]->id()==_progenitor->id()&&!pdata[ix]->stable())
	theChildren[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
	else
	theChildren[ix]->set5Momentum(Lorentz5Momentum(pdata[ix]->mass()));
	}
	// some code moved to updateChildren
	particle->showerKinematics()->
	updateChildren(particle, theChildren, fb.type);
	// In the case of splittings which involves coloured particles,
	// set properly the colour flow of the branching.
	// update the history if needed
	_currenttree->updateInitialStateShowerProduct(_progenitor,theChildren[0]);
	_currenttree->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
	// shower the first particle
	spaceLikeDecayShower(theChildren[0],maxScales,minmass,type);
	// shower the second particle
	timeLikeShower(theChildren[1],type,Branching(),true);
	updateHistory(theChildren[1]);
	// 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 = interactions_[0];
	if(_hardtree&&inter!=ShowerInteraction::Both) {
	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,interactions_[0],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;
	}
	}
	}
	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) {
	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,true);
	if(output) updateHistory(progenitor()->progenitor());
	return output;
	}
	}
	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) {
	if(hardTree()) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit =hardTree()->particles().end(),
	mit = hardTree()->particles().find(progenitor()->progenitor());
	if( mit != eit && mit->second->parent() ) {
	return truncatedSpaceLikeShower( progenitor()->progenitor(),
	parent, mit->second->parent(), type );
	}
	}
	return hardOnly() ? false :
	spaceLikeShower(progenitor()->progenitor(),parent,type);
	}

	bool Evolver::
	startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
	Energy minimumMass,ShowerInteraction::Type type) {
	if(hardTree()) {
	map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
	eit =hardTree()->particles().end(),
	mit = hardTree()->particles().find(progenitor()->progenitor());
	if( mit != eit && mit->second->parent() ) {
	HardBranchingPtr branch=mit->second;
	while(branch->parent()) branch=branch->parent();
	return truncatedSpaceLikeDecayShower(progenitor()->progenitor(),maxScales,
	minimumMass, branch ,type);
	}
	}
	return hardOnly() ? false :
	spaceLikeDecayShower(progenitor()->progenitor(),maxScales,minimumMass,type);
	}

	bool Evolver::timeLikeVetoed(const Branching & fb,
	ShowerParticlePtr particle) {
	// work out type of interaction
	ShowerInteraction::Type type = fb.type==ShowerPartnerType::QED ?
	ShowerInteraction::QED : ShowerInteraction::QCD;
	// 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()->profileScales() &&
	ShowerHandler::currentHandler()->firstInteraction() ) {
	double weight =
	ShowerHandler::currentHandler()->profileScales()->
	hardScaleProfile(_progenitor->maxHardPt(),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 = bb.type==ShowerPartnerType::QED ?
	ShowerInteraction::QED : ShowerInteraction::QCD;
	// 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()->profileScales() &&
	ShowerHandler::currentHandler()->firstInteraction() ) {
	double weight =
	ShowerHandler::currentHandler()->profileScales()->
	hardScaleProfile(_progenitor->maxHardPt(),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 = fb.type==ShowerPartnerType::QED ?
	ShowerInteraction::QED : ShowerInteraction::QCD;
	// 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);
	if(interaction_==4) {
	vector<ShowerInteraction::Type> inter(2);
	inter[0] = ShowerInteraction::QCD;
	inter[1] = ShowerInteraction::QED;
	_hardtree = _hardme->generateHardest( currentTree(),inter );
	}
	else {
	_hardtree = _hardme->generateHardest( currentTree(),interactions_ );
	}
	}
	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()->lastCentralScale()),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,bool first) {
	int ntry=0;
	do {
	++ntry;
	Branching fb;
	unsigned int iout=0;
	tcPDPtr pdata[2];
	while (true) {
	// no truncated shower break
	if(!isTruncatedShowerON()\|\|hardOnly()) break;
	// generate emission
	fb=splittingGenerator()->chooseForwardBranching(*particle,1.,type);
	// no emission break
	if(!fb.kinematics) break;
	// check haven't evolved too far
	if(fb.kinematics->scale() < branch->scale()) {
	fb=Branching();
	break;
	}
	// get the particle data objects
	for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
	if(particle->id()!=fb.ids[0]) {
	for(unsigned int ix=0;ix<2;++ix) {
	tPDPtr cc(pdata[ix]->CC());
	if(cc) pdata[ix]=cc;
	}
	}
	// find the truncated line
	iout=0;
	if(pdata[0]->id()!=pdata[1]->id()) {
	if(pdata[0]->id()==particle->id()) iout=1;
	else if (pdata[1]->id()==particle->id()) iout=2;
	}
	else if(pdata[0]->id()==particle->id()) {
	if(fb.kinematics->z()>0.5) iout=1;
	else iout=2;
	}
	// apply the vetos for the truncated shower
	// no flavour changing branchings
	if(iout==0) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
	// only if same interaction for forced branching
	ShowerInteraction::Type type2 = fb.type==ShowerPartnerType::QED ?
	ShowerInteraction::QED : ShowerInteraction::QCD;
	// 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;
	}
	// should do base class vetos as well
	if(timeLikeVetoed(fb,particle)) {
	particle->vetoEmission(fb.type,fb.kinematics->scale());
	continue;
	}
	break;
	}
	// if no branching force truncated emission
	if(!fb.kinematics) {
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin=
	branch->sudakov()->createFinalStateBranching(branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT());
	showerKin->initialize( *particle,PPtr() );
	IdList idlist(3);
	idlist[0] = particle->id();
	idlist[1] = branch->children()[0]->branchingParticle()->id();
	idlist[2] = branch->children()[1]->branchingParticle()->id();
	fb = Branching( showerKin, idlist, branch->sudakov(),branch->type() );
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics( fb.kinematics );
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// Assign the splitting function to the emitting particle.
	// For the time being we are considering only 1->2 branching
	// Create the ShowerParticle objects for the two children of
	// the emitting particle; set the parent/child relationship
	// if same as definition create particles, otherwise create cc
	ShowerParticleVector theChildren;
	for(unsigned int ix=0;ix<2;++ix) {
	theChildren.push_back(new_ptr(ShowerParticle(branch->children()[ix]->
	branchingParticle()->dataPtr(),true)));
	if(theChildren[ix]->id()==_progenitor->id()&&!theChildren[ix]->dataPtr()->stable())
	theChildren[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
	else
	theChildren[ix]->set5Momentum(Lorentz5Momentum(theChildren[ix]->dataPtr()->mass()));
	}
	particle->showerKinematics()->
	updateChildren(particle, theChildren,fb.type);
	for(unsigned int ix=0;ix<2;++ix) {
	theChildren[ix]->scales().QED = min(theChildren[ix]->scales().QED ,particle->scales().QED );
	theChildren[ix]->scales().QED_noAO = min(theChildren[ix]->scales().QED_noAO ,particle->scales().QED_noAO );
	theChildren[ix]->scales().QCD_c = min(theChildren[ix]->scales().QCD_c ,particle->scales().QCD_c );
	theChildren[ix]->scales().QCD_c_noAO = min(theChildren[ix]->scales().QCD_c_noAO ,particle->scales().QCD_c_noAO );
	theChildren[ix]->scales().QCD_ac = min(theChildren[ix]->scales().QCD_ac ,particle->scales().QCD_ac );
	theChildren[ix]->scales().QCD_ac_noAO = min(theChildren[ix]->scales().QCD_ac_noAO,particle->scales().QCD_ac_noAO);
	}
	// shower the first particle
	if( branch->children()[0]->children().empty() ) {
	if( ! hardOnly() )
	timeLikeShower(theChildren[0],type,Branching(),false);
	}
	else {
	truncatedTimeLikeShower( theChildren[0],branch->children()[0],type,false);
	}
	// shower the second particle
	if( branch->children()[1]->children().empty() ) {
	if( ! hardOnly() )
	timeLikeShower( theChildren[1] , type,Branching(),false);
	}
	else {
	truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type,false);
	}
	// that's if for old approach
	if(_reconOpt==0) return true;
	// branching has happened
	particle->showerKinematics()->updateParent(particle, theChildren,fb.type);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<theChildren.size();++ix)
	particle->abandonChild(theChildren[ix]);
	theChildren.clear();
	continue;
	}
	else {
	if(first&&!theChildren.empty())
	particle->showerKinematics()->resetChildren(particle,theChildren);
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}
	}
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->showerKinematics(fb.kinematics);
	if(fb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(fb.kinematics->pT());
	// Assign the splitting function to the emitting particle.
	// For the time being we are considering only 1->2 branching
	// Create the ShowerParticle objects for the two children of
	// the emitting particle; set the parent/child relationship
	// if same as definition create particles, otherwise create cc
	ShowerParticleVector theChildren;
	for(unsigned int ix=0;ix<2;++ix) {
	theChildren.push_back( new_ptr( ShowerParticle( pdata[ix], true ) ) );
	if(theChildren[ix]->id()==_progenitor->id()&&!pdata[ix]->stable())
	theChildren[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
	else
	theChildren[ix]->set5Momentum(Lorentz5Momentum(pdata[ix]->mass()));
	}
	particle->showerKinematics()->
	updateChildren( particle, theChildren , fb.type);
	// shower the first particle
	if( iout == 1 ) truncatedTimeLikeShower( theChildren[0], branch , type ,false);
	else timeLikeShower( theChildren[0] , type,Branching(),false);
	// shower the second particle
	if( iout == 2 ) truncatedTimeLikeShower( theChildren[1], branch , type ,false);
	else timeLikeShower( theChildren[1] , type,Branching(),false);
	// that's if for old approach
	if(_reconOpt==0) return true;
	// branching has happened
	particle->showerKinematics()->updateParent(particle, theChildren,fb.type);
	// clean up the vetoed emission
	if(particle->virtualMass()==ZERO) {
	particle->showerKinematics(ShoKinPtr());
	for(unsigned int ix=0;ix<theChildren.size();++ix)
	particle->abandonChild(theChildren[ix]);
	theChildren.clear();
	}
	else {
	if(first&&!theChildren.empty())
	particle->showerKinematics()->resetChildren(particle,theChildren);
	if(particle->spinInfo()) particle->spinInfo()->develop();
	return true;
	}
	}
	while(ntry<50);
	return false;
	}

	bool Evolver::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
	HardBranchingPtr branch,
	ShowerInteraction::Type type) {
	tcPDFPtr pdf;
	if(ShowerHandler::currentHandler()->firstPDF().particle() == beamParticle())
	pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
	if(ShowerHandler::currentHandler()->secondPDF().particle() == beamParticle())
	pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
	Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
	Branching bb;
	// parameters of the force branching
	double z(0.);
	HardBranchingPtr timelike;
	for( unsigned int ix = 0; ix < branch->children().size(); ++ix ) {
	if( branch->children()[ix]->status() ==HardBranching::Outgoing) {
	timelike = branch->children()[ix];
	}
	if( branch->children()[ix]->status() ==HardBranching::Incoming )
	z = branch->children()[ix]->z();
	}
	// generate truncated branching
	tcPDPtr part[2];
	if(z>=0.&&z<=1.) {
	while (true) {
	if( !isTruncatedShowerON() \|\| hardOnly() ) break;
	bb = splittingGenerator()->chooseBackwardBranching( *particle,
	beam, 1., beamParticle(),
	type , pdf,freeze);
	if( !bb.kinematics \|\| bb.kinematics->scale() < branch->scale() ) {
	bb = Branching();
	break;
	}
	// particles as in Sudakov form factor
	part[0] = getParticleData( bb.ids[0] );
	part[1] = getParticleData( bb.ids[2] );

	//is emitter anti-particle
	if( particle->id() != bb.ids[1]) {
	if( part[0]->CC() ) part[0] = part[0]->CC();
	if( part[1]->CC() ) part[1] = part[1]->CC();
	}
	double zsplit = bb.kinematics->z();
	// apply the vetos for the truncated shower
	// if doesn't carry most of momentum
	ShowerInteraction::Type type2 = bb.type==ShowerPartnerType::QED ?
	ShowerInteraction::QED : ShowerInteraction::QCD;
	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() );
	kinematics->initialize( *particle, beam );
	// assign the splitting function and shower kinematics
	particle->showerKinematics( kinematics );
	if(kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent =
	new_ptr( ShowerParticle( branch->branchingParticle()->dataPtr(), false ) );
	ShowerParticlePtr otherChild =
	new_ptr( ShowerParticle( timelike->branchingParticle()->dataPtr(),
	true, true ) );
	ShowerParticleVector theChildren;
	theChildren.push_back( particle );
	theChildren.push_back( otherChild );
	particle->showerKinematics()->
	updateParent( newParent, theChildren, branch->type());
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
	currentTree()->addInitialStateBranching( particle, newParent, otherChild );
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	bool emitted=false;
	if(!hardOnly()) {
	if( branch->parent() ) {
	emitted = truncatedSpaceLikeShower( newParent, beam, branch->parent() , type);
	}
	else {
	emitted = spaceLikeShower( newParent, beam , type);
	}
	}
	if( !emitted ) {
	if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
	kinematics->updateLast( newParent, ZERO, ZERO );
	}
	else {
	pair<Energy,double> kt = intrinsicpT()[progenitor()];
	kinematics->updateLast( newParent,
	kt.first*cos( kt.second ),
	kt.first*sin( kt.second ) );
	}
	}
	particle->showerKinematics()->
	updateChildren( newParent, theChildren,bb.type);
	if(hardOnly()) return true;
	// perform the shower of the final-state particle
	if( timelike->children().empty() ) {
	timeLikeShower( otherChild , type,Branching(),true);
	}
	else {
	truncatedTimeLikeShower( otherChild, timelike , type,true);
	}
	updateHistory(otherChild);
	// return the emitted
	return true;
	}
	// assign the splitting function and shower kinematics
	particle->showerKinematics( bb.kinematics );
	if(bb.kinematics->pT()>progenitor()->highestpT())
	progenitor()->highestpT(bb.kinematics->pT());
	// For the time being we are considering only 1->2 branching
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent = new_ptr( ShowerParticle( part[0], false ) );
	ShowerParticlePtr otherChild = new_ptr( ShowerParticle( part[1], true, true ) );
	ShowerParticleVector theChildren;
	theChildren.push_back( particle );
	theChildren.push_back( otherChild );
	particle->showerKinematics()->
	updateParent( newParent, theChildren, bb.type);
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
	currentTree()->addInitialStateBranching( particle, newParent, otherChild );
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	bool emitted = truncatedSpaceLikeShower( newParent, beam, branch,type);
	// now reconstruct the momentum
	if( !emitted ) {
	if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
	bb.kinematics->updateLast( newParent, ZERO, ZERO );
	}
	else {
	pair<Energy,double> kt = intrinsicpT()[ progenitor() ];
	bb.kinematics->updateLast( newParent,
	kt.first*cos( kt.second ),
	kt.first*sin( kt.second ) );
	}
	}
	particle->showerKinematics()->
	updateChildren( newParent, theChildren, bb.type);
	// perform the shower of the final-state particle
	timeLikeShower( otherChild , type,Branching(),true);
	updateHistory(otherChild);
	// return the emitted
	return true;
	}

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

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

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

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

	}
	else
	assert(false);
	}
	}
	}

	void Evolver::constructSpaceLikeLine(tShowerParticlePtr particle,
	HardBranchingPtr & first,
	HardBranchingPtr & last,
	SudakovPtr sud,PPtr beam) {
	if(!particle) return;
	if(!particle->parents().empty()) {
	tShowerParticlePtr parent =
	dynamic_ptr_cast<ShowerParticlePtr>(particle->parents()[0]);
	SudakovPtr newSud=particle->showerKinematics()->SudakovFormFactor();
	constructSpaceLikeLine(parent,first,last,newSud,beam);
	}
	HardBranchingPtr newBranch =
	new_ptr(HardBranching(particle,sud,last,HardBranching::Incoming));
	newBranch->beam(beam);
	if(!first) {
	first=newBranch;
	last =newBranch;
	return;
	}
	last->addChild(newBranch);
	tShowerParticlePtr timeChild =
	dynamic_ptr_cast<ShowerParticlePtr>(particle->parents()[0]->children()[1]);
	HardBranchingPtr timeBranch;
	if(!timeChild->children().empty()) {
	timeBranch =
	new_ptr(HardBranching(timeChild,
	timeChild->showerKinematics()->SudakovFormFactor(),
	last,HardBranching::Outgoing));
	constructTimeLikeLine(timeBranch,timeChild);
	}
	else {
	timeBranch =
	new_ptr(HardBranching(timeChild,SudakovPtr(),last,HardBranching::Outgoing));
	}
	last->addChild(timeBranch);
	// sort out the type
	if(last->branchingParticle() ->id() == ParticleID::gamma \|\|
	newBranch->branchingParticle() ->id() == ParticleID::gamma \|\|
	timeBranch->branchingParticle()->id() == ParticleID::gamma) {
	last->type(ShowerPartnerType::QED);
	}
	else if(last->branchingParticle()->id()==newBranch->branchingParticle()->id()) {
	if(last->branchingParticle()->id()==ParticleID::g) {
	if(last->branchingParticle()->colourLine()==
	newBranch->branchingParticle()->colourLine()) {
	last->type(ShowerPartnerType::QCDAntiColourLine);
	}
	else {
	last->type(ShowerPartnerType::QCDColourLine);
	}
	}
	else if(last->branchingParticle()->hasColour()) {
	last->type(ShowerPartnerType::QCDColourLine);
	}
	else if(last->branchingParticle()->hasAntiColour()) {
	last->type(ShowerPartnerType::QCDAntiColourLine);
	}
	else
	assert(false);
	}
	else if(newBranch->branchingParticle()->id()==ParticleID::g) {
	if(last->branchingParticle()->hasColour()) {
	last->type(ShowerPartnerType::QCDAntiColourLine);
	}
	else if(last->branchingParticle()->hasAntiColour()) {
	last->type(ShowerPartnerType::QCDColourLine);
	}
	else
	assert(false);
	}
	else if(newBranch->branchingParticle()->hasColour()) {
	last->type(ShowerPartnerType::QCDColourLine);
	}
	else if(newBranch->branchingParticle()->hasAntiColour()) {
	last->type(ShowerPartnerType::QCDAntiColourLine);
	}
	else {
	assert(false);
	}
	last=newBranch;
	}

	void Evolver::connectTrees(ShowerTreePtr showerTree,
	HardTreePtr hardTree, bool hard ) {
	ShowerParticleVector particles;
	// find the Sudakovs
	for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
	cit!=hardTree->branchings().end();++cit) {
	// Sudakovs for ISR
	if((cit).parent()&&(cit).status()==HardBranching::Incoming) {
	++_nis;
	IdList br(3);
	br[0] = (**cit).parent()->branchingParticle()->id();
	br[1] = (**cit). branchingParticle()->id();
	br[2] = (*cit).parent()->children()[0]==cit ?
	(**cit).parent()->children()[1]->branchingParticle()->id() :
	(**cit).parent()->children()[0]->branchingParticle()->id();
	BranchingList branchings = splittingGenerator()->initialStateBranchings();
	if(br[1]<0&&br[0]==br[1]) {
	br[0] = abs(br[0]);
	br[1] = abs(br[1]);
	}
	else if(br[1]<0) {
	br[1] = -br[1];
	br[2] = -br[2];
	}
	long index = abs(br[1]);
	SudakovPtr sudakov;
	for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
	cjt != branchings.upper_bound(index); ++cjt ) {
	IdList ids = cjt->second.second;
	if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
	sudakov=cjt->second.first;
	break;
	}
	}
	if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
	<< "Evolver::connectTrees() for ISR"
	<< Exception::runerror;
	(**cit).parent()->sudakov(sudakov);
	}
	// Sudakovs for FSR
	else if(!(**cit).children().empty()) {
	++_nfs;
	IdList br(3);
	br[0] = (**cit) .branchingParticle()->id();
	br[1] = (**cit).children()[0]->branchingParticle()->id();
	br[2] = (**cit).children()[1]->branchingParticle()->id();
	BranchingList branchings = splittingGenerator()->finalStateBranchings();
	if(br[0]<0) {
	br[0] = abs(br[0]);
	br[1] = abs(br[1]);
	br[2] = abs(br[2]);
	}
	long index = br[0];
	SudakovPtr sudakov;
	for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
	cjt != branchings.upper_bound(index); ++cjt ) {
	IdList ids = cjt->second.second;
	if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
	sudakov=cjt->second.first;
	break;
	}
	}
	if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
	<< "Evolver::connectTrees()"
	<< Exception::runerror;
	(**cit).sudakov(sudakov);
	}
	}
	// calculate the evolution scale
	for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
	cit!=hardTree->branchings().end();++cit) {
	particles.push_back((*cit)->branchingParticle());
	}
	showerModel()->partnerFinder()->
	setInitialEvolutionScales(particles,!hard,hardTree->interaction(),
	!hardTree->partnersSet());
	hardTree->partnersSet(true);
	// inverse reconstruction
	if(hard) {
	showerModel()->kinematicsReconstructor()->
	deconstructHardJets(hardTree,ShowerHandler::currentHandler()->evolver(),
	hardTree->interaction());
	}
	else
	showerModel()->kinematicsReconstructor()->
	deconstructDecayJets(hardTree,ShowerHandler::currentHandler()->evolver(),
	hardTree->interaction());
	// now reset the momenta of the showering particles
	vector<ShowerProgenitorPtr> particlesToShower;
	for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit=showerTree->incomingLines().begin();
	cit!=showerTree->incomingLines().end();++cit )
	particlesToShower.push_back(cit->first);
	// extract the showering particles
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cit=showerTree->outgoingLines().begin();
	cit!=showerTree->outgoingLines().end();++cit )
	particlesToShower.push_back(cit->first);
	// match them
	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) {
	// order of the interactions
	bool showerOrder(true);
	// 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);
	// 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;
	}
	}
	// check if interactions in right order
	if(hardTree() && interaction_!=4 &&
	hardTree()->interaction()!=interactions_[0]) {
	assert(interactions_.size()==2);
	showerOrder = false;
	swap(interactions_[0],interactions_[1]);
	}
	// loop over possible interactions
	for(unsigned int inter=0;inter<interactions_.size();++inter) {
	// set up for second pass if required
	if(inter!=0) {
	// zero intrinsic pt so only added first time round
	intrinsicpT().clear();
	// construct the tree and throw veto if not possible
	if(!(hard ?
	constructHardTree (particlesToShower,interactions_[inter]) :
	constructDecayTree(particlesToShower,interactions_[inter])))
	throw InteractionVeto();
	}
	+ // create random particle vector (only need to do once)
	+ vector<ShowerProgenitorPtr> tmp;
	+ unsigned int nColouredIncoming = 0;
	+ while(particlesToShower.size()>0){
	+ unsigned int xx=UseRandom::irnd(particlesToShower.size());
	+ tmp.push_back(particlesToShower[xx]);
	+ particlesToShower.erase(particlesToShower.begin()+xx);
	+ }
	+ particlesToShower=tmp;
	+ for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	+ if(!particlesToShower[ix]->progenitor()->isFinalState() &&
	+ particlesToShower[ix]->progenitor()->coloured()) ++nColouredIncoming;
	+ }
	+ bool switchRecon = hard && nColouredIncoming !=1;
	// main shower loop
	unsigned int ntry(0);
	bool reconstructed = false;
	do {
	// clear results of last attempt if needed
	if(ntry!=0) {
	currentTree()->clear();
	setEvolutionPartners(hard,interactions_[inter],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());
	}
	}
	}
	- // create random particle vector
	- vector<ShowerProgenitorPtr> tmp;
	- while(particlesToShower.size()>0){
	- unsigned int xx=UseRandom::irnd(particlesToShower.size());
	- tmp.push_back(particlesToShower[xx]);
	- particlesToShower.erase(particlesToShower.begin()+xx);
	- }
	- particlesToShower=tmp;
	// 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(interactions_[inter]));
	}
	// initial-state radiation
	else {
	if(!isISRadiationON()) continue;
	// hard process
	if(hard) {
	// get the PDF
	setBeamParticle(_progenitor->beam());
	assert(beamParticle());
	// perform the shower
	// set the beam particle
	tPPtr beamparticle=progenitor()->original();
	if(!beamparticle->parents().empty())
	beamparticle=beamparticle->parents()[0];
	// generate the shower
	progenitor()->hasEmitted(startSpaceLikeShower(beamparticle,
	interactions_[inter]));
	}
	// decay
	else {
	// skip colour and electrically neutral particles
	if(!progenitor()->progenitor()->dataPtr()->coloured() &&
	!progenitor()->progenitor()->dataPtr()->charged()) {
	progenitor()->hasEmitted(false);
	continue;
	}
	// perform shower
	// set the scales correctly. The current scale is the maximum scale for
	// emission not the starting scale
	ShowerParticle::EvolutionScales maxScales(progenitor()->progenitor()->scales());
	progenitor()->progenitor()->scales() = ShowerParticle::EvolutionScales();
	if(progenitor()->progenitor()->dataPtr()->charged()) {
	progenitor()->progenitor()->scales().QED = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QED_noAO = progenitor()->progenitor()->mass();
	}
	if(progenitor()->progenitor()->hasColour()) {
	progenitor()->progenitor()->scales().QCD_c = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QCD_c_noAO = progenitor()->progenitor()->mass();
	}
	if(progenitor()->progenitor()->hasAntiColour()) {
	progenitor()->progenitor()->scales().QCD_ac = progenitor()->progenitor()->mass();
	progenitor()->progenitor()->scales().QCD_ac_noAO = progenitor()->progenitor()->mass();
	}
	// perform the shower
	progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass,
	interactions_[inter]));
	}
	}
	}
	// do the kinematic reconstruction, checking if it worked
	reconstructed = hard ?
	showerModel()->kinematicsReconstructor()->
	- reconstructHardJets (currentTree(),intrinsicpT(),interactions_[inter]) :
	+ reconstructHardJets (currentTree(),intrinsicpT(),interactions_[inter],
	+ switchRecon && ntry>maximumTries()/2) :
	showerModel()->kinematicsReconstructor()->
	reconstructDecayJets(currentTree(),interactions_[inter]);
	}
	while(!reconstructed&&maximumTries()>++ntry);
	// check if failed to generate the shower
	if(ntry==maximumTries()) {
	if(hard)
	throw ShowerHandler::ShowerTriesVeto(ntry);
	else
	throw Exception() << "Failed to generate the shower after "
	<< ntry << " attempts in Evolver::showerDecay()"
	<< Exception::eventerror;
	}
	}
	// tree has now showered
	_currenttree->hasShowered(true);
	if(!showerOrder) swap(interactions_[0],interactions_[1]);
	hardTree(HardTreePtr());
	}

	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);
	}
	diff --git a/Shower/Base/KinematicsReconstructor.h b/Shower/Base/KinematicsReconstructor.h
	--- a/Shower/Base/KinematicsReconstructor.h
	+++ b/Shower/Base/KinematicsReconstructor.h
	@@ -1,134 +1,135 @@
	// -- C++ --
	//
	// KinematicsReconstructor.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_KinematicsReconstructor_H
	#define HERWIG_KinematicsReconstructor_H
	//
	// This is the declaration of the KinematicsReconstructor class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "ShowerParticle.h"
	#include "ShowerProgenitor.h"
	#include "ShowerTree.h"
	#include "HardTree.h"
	#include "Evolver.fh"
	#include "KinematicsReconstructor.fh"
	#include <cassert>

	namespace Herwig {

	using namespace ThePEG;

	/**\ingroup Shower
	* Exception class
	* used to communicate failure of kinematics
	* reconstruction.
	*/
	struct KinematicsReconstructionVeto {};


	/** \ingroup Shower
	*
	* This class is responsible for the kinematical reconstruction
	* after each showering step, and also for the necessary Lorentz boosts
	* in order to preserve energy-momentum conservation in the overall collision,
	* and also the invariant mass and the rapidity of the hard subprocess system.
	* In the case of multi-step showering, there will be not unnecessary
	* kinematical reconstructions.
	*
	* Notice:
	* - although we often use the term "jet" in either methods or variables names,
	* or in comments, which could appear applicable only for QCD showering,
	* there is indeed no "dynamics" represented in this class: only kinematics
	* is involved, as the name of this class remainds. Therefore it can be used
	* for any kind of showers (QCD-,QED-,EWK-,... bremsstrahlung).
	*
	* @see \ref KinematicsReconstructorInterfaces "The interfaces"
	* defined for KinematicsReconstructor.
	*/
	class KinematicsReconstructor: public Interfaced {

	public:

	/**
	* Methods to reconstruct the kinematics of a scattering or decay process
	*/
	//@{
	/**
	* Given the ShowerTree for the shower from a hard process
	* the method does the reconstruction of the jets,
	* including the appropriate boosts (kinematics reshufflings)
	* needed to conserve the total energy-momentum of the collision
	* and preserving the invariant mass and the rapidity of the
	* hard subprocess system.
	*/
	virtual bool reconstructHardJets(ShowerTreePtr hard,
	const map<tShowerProgenitorPtr,
	pair<Energy,double> > & pt,
	- ShowerInteraction::Type type) const=0;
	+ ShowerInteraction::Type type,
	+ bool switchRecon) const=0;

	/**
	* Given the ShowerTree for a decay shower
	* the method does the reconstruction of the jets,
	* including the appropriate boosts (kinematics reshufflings)
	* needed to conserve the total energy-momentum of the collision
	* and preserving the invariant mass and the rapidity of the
	* hard subprocess system.
	*/
	virtual bool reconstructDecayJets(ShowerTreePtr decay,
	ShowerInteraction::Type type) const=0;
	//@}

	/**
	* Methods to invert the reconstruction of the shower for
	* a scattering or decay process and calculate
	* the variables used to generate the
	* shower given the particles produced.
	* This is needed for the CKKW and POWHEG approaches
	*/
	//@{
	/**
	* Given the particles, with a history which we wish to interpret
	* as a shower reconstruct the variables used to generate the
	* shower for a decay process
	*/
	virtual bool deconstructDecayJets(HardTreePtr decay,cEvolverPtr,
	ShowerInteraction::Type) const=0;

	/**
	* Given the particles, with a history which we wish to interpret
	* as a shower reconstruct the variables used to generate the shower
	* for a hard process
	*/
	virtual bool deconstructHardJets(HardTreePtr hard,cEvolverPtr,
	ShowerInteraction::Type) const=0;
	//@}

	public:

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

	private:

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

	}

	#endif /* HERWIG_KinematicsReconstructor_H */
	diff --git a/Shower/Default/QTildeReconstructor.cc b/Shower/Default/QTildeReconstructor.cc
	--- a/Shower/Default/QTildeReconstructor.cc
	+++ b/Shower/Default/QTildeReconstructor.cc
	@@ -1,2944 +1,2945 @@
	// -- C++ --
	//
	// QTildeReconstructor.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 QTildeReconstructor class.
	//

	#include "QTildeReconstructor.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/EventRecord/Event.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "Herwig/Shower/Base/Evolver.h"
	#include "Herwig/Shower/Base/PartnerFinder.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig/Shower/SplittingFunctions/SplittingFunction.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/EventRecord/ColourLine.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include <cassert>

	using namespace Herwig;

	DescribeClass<QTildeReconstructor,KinematicsReconstructor>
	describeQTildeReconstructor("Herwig::QTildeReconstructor", "HwShower.so");

	namespace {

	/**
	* Struct to order the jets in off-shellness
	*/
	struct JetOrdering {

	bool operator() (const JetKinStruct & j1, const JetKinStruct & j2) {
	Energy diff1 = j1.q.m()-j1.p.m();
	Energy diff2 = j2.q.m()-j2.p.m();
	if(diff1!=diff2) {
	return diff1>diff2;
	}
	else if( j1.q.e() != j2.q.e() )
	return j1.q.e()>j2.q.e();
	else
	return j1.parent->uniqueId>j2.parent->uniqueId;
	}
	};

	}

	void QTildeReconstructor::persistentOutput(PersistentOStream & os) const {
	os << _reconopt << _initialBoost << ounit(_minQ,GeV) << _noRescale
	<< _noRescaleVector << _finalStateReconOption
	<< _initialStateReconOption;
	}

	void QTildeReconstructor::persistentInput(PersistentIStream & is, int) {
	is >> _reconopt >> _initialBoost >> iunit(_minQ,GeV) >> _noRescale
	>> _noRescaleVector >> _finalStateReconOption
	>> _initialStateReconOption;
	}

	void QTildeReconstructor::Init() {

	static ClassDocumentation<QTildeReconstructor> documentation
	( "This class is responsible for the kinematics reconstruction of the showering,",
	" including the kinematics reshuffling necessary to compensate for the recoil"
	"of the emissions." );

	static Switch<QTildeReconstructor,unsigned int> interfaceReconstructionOption
	("ReconstructionOption",
	"Option for the kinematics reconstruction",
	&QTildeReconstructor::_reconopt, 0, false, false);
	static SwitchOption interfaceReconstructionOptionGeneral
	(interfaceReconstructionOption,
	"General",
	"Use the general solution which ignores the colour structure for all processes",
	0);
	static SwitchOption interfaceReconstructionOptionColour
	(interfaceReconstructionOption,
	"Colour",
	"Use the colour structure of the process to determine the reconstruction procedure.",
	1);
	static SwitchOption interfaceReconstructionOptionColour2
	(interfaceReconstructionOption,
	"Colour2",
	"Make the most use possible of the colour structure of the process to determine the reconstruction procedure. "
	"Start with FF, then IF then II colour connections",
	2);
	static SwitchOption interfaceReconstructionOptionColour3
	(interfaceReconstructionOption,
	"Colour3",
	"Make the most use possible of the colour structure of the process to determine the reconstruction procedure. "
	"Do the colour connections in order of the pT's emitted in the shower starting with the hardest."
	" The colour partner is fully reconstructed at the same time.",
	3);
	static SwitchOption interfaceReconstructionOptionColour4
	(interfaceReconstructionOption,
	"Colour4",
	"Make the most use possible of the colour structure of the process to determine the reconstruction procedure. "
	"Do the colour connections in order of the pT's emitted in the shower starting with the hardest, while leaving"
	" the colour partner on mass-shell",
	4);

	static Parameter<QTildeReconstructor,Energy> interfaceMinimumQ2
	("MinimumQ2",
	"The minimum Q2 for the reconstruction of initial-final systems",
	&QTildeReconstructor::_minQ, GeV, 0.001GeV, 1e-6GeV, 10.0*GeV,
	false, false, Interface::limited);

	static RefVector<QTildeReconstructor,ParticleData> interfaceNoRescale
	("NoRescale",
	"Particles which shouldn't be rescaled to be on shell by the shower",
	&QTildeReconstructor::_noRescaleVector, -1, false, false, true, false, false);

	static Switch<QTildeReconstructor,unsigned int> interfaceInitialInitialBoostOption
	("InitialInitialBoostOption",
	"Option for how the boost from the system before ISR to that after ISR is applied.",
	&QTildeReconstructor::_initialBoost, 0, false, false);
	static SwitchOption interfaceInitialInitialBoostOptionOneBoost
	(interfaceInitialInitialBoostOption,
	"OneBoost",
	"Apply one boost from old CMS to new CMS",
	0);
	static SwitchOption interfaceInitialInitialBoostOptionLongTransBoost
	(interfaceInitialInitialBoostOption,
	"LongTransBoost",
	"First apply a longitudinal and then a transverse boost",
	1);

	static Switch<QTildeReconstructor,unsigned int> interfaceFinalStateReconOption
	("FinalStateReconOption",
	"Option for how to reconstruct the momenta of the final-state system",
	&QTildeReconstructor::_finalStateReconOption, 0, false, false);
	static SwitchOption interfaceFinalStateReconOptionDefault
	(interfaceFinalStateReconOption,
	"Default",
	"All the momenta are rescaled in the rest frame",
	0);
	static SwitchOption interfaceFinalStateReconOptionMostOffShell
	(interfaceFinalStateReconOption,
	"MostOffShell",
	"All particles put on the new-mass shell and then the most off-shell and"
	" recoiling system are rescaled to ensure 4-momentum is conserved.",
	1);
	static SwitchOption interfaceFinalStateReconOptionRecursive
	(interfaceFinalStateReconOption,
	"Recursive",
	"Recursively put on shell by putting the most off-shell particle which"
	" hasn't been rescaled on-shell by rescaling the particles and the recoiling system. ",
	2);
	static SwitchOption interfaceFinalStateReconOptionRestMostOffShell
	(interfaceFinalStateReconOption,
	"RestMostOffShell",
	"The most off-shell is put on shell by rescaling it and the recoiling system,"
	" the recoiling system is then put on-shell in its rest frame.",
	3);
	static SwitchOption interfaceFinalStateReconOptionRestRecursive
	(interfaceFinalStateReconOption,
	"RestRecursive",
	"As 3 but recursive treated the currently most-off shell,"
	" only makes a difference if more than 3 partons.",
	4);

	static Switch<QTildeReconstructor,unsigned int> interfaceInitialStateReconOption
	("InitialStateReconOption",
	"Option for the reconstruction of initial state radiation",
	&QTildeReconstructor::_initialStateReconOption, 0, false, false);
	static SwitchOption interfaceInitialStateReconOptionRapidity
	(interfaceInitialStateReconOption,
	"Rapidity",
	"Preserve shat and rapidity",
	0);
	static SwitchOption interfaceInitialStateReconOptionLongitudinal
	(interfaceInitialStateReconOption,
	"Longitudinal",
	"Preserve longitudinal momentum",
	1);
	static SwitchOption interfaceInitialStateReconOptionSofterFraction
	(interfaceInitialStateReconOption,
	"SofterFraction",
	"Preserve the momentum fraction of the parton which has emitted softer.",
	2);

	}

	void QTildeReconstructor::doinit() {
	KinematicsReconstructor::doinit();
	_noRescale = set<cPDPtr>(_noRescaleVector.begin(),_noRescaleVector.end());
	}

	bool QTildeReconstructor::
	reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent) const {
	assert(particleJetParent);
	bool emitted=true;
	// if this is not a fixed point in the reconstruction
	if( !particleJetParent->children().empty() ) {
	// if not a reconstruction fixpoint, dig deeper for all children:
	for ( ParticleVector::const_iterator cit =
	particleJetParent->children().begin();
	cit != particleJetParent->children().end(); ++cit )
	reconstructTimeLikeJet(dynamic_ptr_cast<ShowerParticlePtr>(*cit));
	}
	// it is a reconstruction fixpoint, ie kinematical data has to be available
	else {
	// check if the parent was part of the shower
	ShowerParticlePtr jetGrandParent;
	if(!particleJetParent->parents().empty())
	jetGrandParent= dynamic_ptr_cast<ShowerParticlePtr>
	(particleJetParent->parents()[0]);
	// update if so
	if (jetGrandParent) {
	if (jetGrandParent->showerKinematics()) {
	if(particleJetParent->id()==_progenitor->id()&&
	!_progenitor->data().stable()) {
	jetGrandParent->showerKinematics()->reconstructLast(particleJetParent,
	_progenitor->mass());
	}
	else {
	jetGrandParent->showerKinematics()->reconstructLast(particleJetParent);
	}
	}
	}
	// otherwise
	else {
	Energy dm = particleJetParent->data().constituentMass();
	if (abs(dm-particleJetParent->momentum().m())>0.001*MeV
	&&particleJetParent->dataPtr()->stable()
	&&particleJetParent->id()!=ParticleID::gamma
	&&_noRescale.find(particleJetParent->dataPtr())==_noRescale.end()) {
	Lorentz5Momentum dum = particleJetParent->momentum();
	dum.setMass(dm);
	dum.rescaleEnergy();
	particleJetParent->set5Momentum(dum);
	}
	else {
	emitted=false;
	}
	}
	}
	// recursion has reached an endpoint once, ie we can reconstruct the
	// kinematics from the children.
	if( !particleJetParent->children().empty() )
	particleJetParent->showerKinematics()
	->reconstructParent( particleJetParent, particleJetParent->children() );
	return emitted;
	}

	bool QTildeReconstructor::
	reconstructHardJets(ShowerTreePtr hard,
	const map<tShowerProgenitorPtr,
	pair<Energy,double> > & intrinsic,
	- ShowerInteraction::Type type) const {
	+ ShowerInteraction::Type type,
	+ bool switchRecon) const {
	_currentTree = hard;
	_intrinsic=intrinsic;
	// extract the particles from the ShowerTree
	vector<ShowerProgenitorPtr> ShowerHardJets=hard->extractProgenitors();
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	_boosts[ShowerHardJets[ix]->progenitor()] = vector<LorentzRotation>();
	}
	for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	tit = _currentTree->treelinks().begin();
	tit != _currentTree->treelinks().end();++tit) {
	_treeBoosts[tit->first] = vector<LorentzRotation>();
	}
	try {
	// old recon method, using new member functions
	- if(_reconopt == 0 ) {
	+ if(_reconopt == 0 \|\| switchRecon ) {
	reconstructGeneralSystem(ShowerHardJets);
	}
	// reconstruction based on coloured systems
	else if( _reconopt == 1) {
	reconstructColourSinglets(ShowerHardJets,type);
	}
	// reconstruction of FF, then IF, then II
	else if( _reconopt == 2) {
	reconstructFinalFirst(ShowerHardJets);
	}
	// reconstruction based on coloured systems
	else if( _reconopt == 3 \|\| _reconopt == 4) {
	reconstructColourPartner(ShowerHardJets);
	}
	else
	assert(false);
	}
	catch(KinematicsReconstructionVeto) {
	_progenitor=tShowerParticlePtr();
	_intrinsic.clear();
	for(map<tPPtr,vector<LorentzRotation> >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) {
	for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
	LorentzRotation rot = rit->inverse();
	bit->first->transform(rot);
	}
	}
	_boosts.clear();
	for(map<tShowerTreePtr,vector<LorentzRotation> >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) {
	for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
	LorentzRotation rot = rit->inverse();
	bit->first->transform(rot,false);
	}
	}
	_currentTree = tShowerTreePtr();
	_treeBoosts.clear();
	return false;
	}
	catch (Exception & ex) {
	_progenitor=tShowerParticlePtr();
	_intrinsic.clear();
	_currentTree = tShowerTreePtr();
	_boosts.clear();
	_treeBoosts.clear();
	throw ex;
	}
	_progenitor=tShowerParticlePtr();
	_intrinsic.clear();
	// ensure x<1
	for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit=hard->incomingLines().begin();cit!=hard->incomingLines().end();++cit) {
	tPPtr parent = cit->first->progenitor();
	while (!parent->parents().empty()) {
	parent = parent->parents()[0];
	}
	tPPtr hadron;
	if ( cit->first->original()->parents().empty() ) {
	hadron = cit->first->original();
	}
	else {
	hadron = cit->first->original()->parents()[0];
	}
	if( ! (hadron->id() == parent->id() && hadron->children().size() <= 1)
	&& parent->momentum().rho() > hadron->momentum().rho()) {
	_progenitor=tShowerParticlePtr();
	_intrinsic.clear();
	for(map<tPPtr,vector<LorentzRotation> >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) {
	for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
	LorentzRotation rot = rit->inverse();
	bit->first->transform(rot);
	}
	}
	_boosts.clear();
	for(map<tShowerTreePtr,vector<LorentzRotation> >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) {
	for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
	LorentzRotation rot = rit->inverse();
	bit->first->transform(rot,false);
	}
	}
	_currentTree = tShowerTreePtr();
	_treeBoosts.clear();
	return false;
	}
	}
	_boosts.clear();
	_treeBoosts.clear();
	_currentTree = tShowerTreePtr();
	return true;
	}

	double
	QTildeReconstructor::solveKfactor(const Energy & root_s,
	const JetKinVect & jets) const {
	Energy2 s = sqr(root_s);
	// must be at least two jets
	if ( jets.size() < 2) throw KinematicsReconstructionVeto();
	// sum of jet masses must be less than roots
	if(momConsEq( 0.0, root_s, jets )>ZERO) throw KinematicsReconstructionVeto();
	// if two jets simple solution
	if ( jets.size() == 2 ) {
	static const Energy2 eps = 1.0e-4 * MeV2;
	if ( sqr(jets[0].p.x()+jets[1].p.x()) < eps &&
	sqr(jets[0].p.y()+jets[1].p.y()) < eps &&
	sqr(jets[0].p.z()+jets[1].p.z()) < eps ) {
	Energy test = (jets[0].p+jets[1].p).vect().mag();
	if(test > 1.0e-4 * MeV) throw KinematicsReconstructionVeto();
	if ( jets[0].p.vect().mag2() < eps ) throw KinematicsReconstructionVeto();
	Energy2 m1sq(jets[0].q.m2()),m2sq(jets[1].q.m2());
	return sqrt( ( sqr(s - m1sq - m2sq) - 4.m1sqm2sq )
	/(4.sjets[0].p.vect().mag2()) );
	}
	else throw KinematicsReconstructionVeto();
	}
	// i.e. jets.size() > 2, numerically
	// check convergence, if it's a problem maybe use Newton iteration?
	else {
	double k1 = 0.,k2 = 1.,k = 0.;
	if ( momConsEq( k1, root_s, jets ) < ZERO ) {
	while ( momConsEq( k2, root_s, jets ) < ZERO ) {
	k1 = k2;
	k2 *= 2;
	}
	while ( fabs( (k1 - k2)/(k1 + k2) ) > 1.e-10 ) {
	if( momConsEq( k2, root_s, jets ) == ZERO ) {
	return k2;
	} else {
	k = (k1+k2)/2.;
	if ( momConsEq( k, root_s, jets ) > ZERO ) {
	k2 = k;
	} else {
	k1 = k;
	}
	}
	}
	return k1;
	} else throw KinematicsReconstructionVeto();
	}
	throw KinematicsReconstructionVeto();
	}

	bool QTildeReconstructor::
	reconstructSpaceLikeJet( const tShowerParticlePtr p) const {
	bool emitted = true;
	tShowerParticlePtr child;
	tShowerParticlePtr parent;
	if(!p->parents().empty())
	parent = dynamic_ptr_cast<ShowerParticlePtr>(p->parents()[0]);
	if(parent) {
	emitted=true;
	reconstructSpaceLikeJet(parent);
	}
	// if branching reconstruct time-like child
	if(p->children().size()==2)
	child = dynamic_ptr_cast<ShowerParticlePtr>(p->children()[1]);
	if(p->perturbative()==0 && child) {
	dynamic_ptr_cast<ShowerParticlePtr>(p->children()[0])->
	showerKinematics()->reconstructParent(p,p->children());
	if(!child->children().empty()) {
	_progenitor=child;
	reconstructTimeLikeJet(child);
	// calculate the momentum of the particle
	Lorentz5Momentum pnew=p->momentum()-child->momentum();
	pnew.rescaleMass();
	p->children()[0]->set5Momentum(pnew);
	}
	}
	return emitted;
	}

	Boost QTildeReconstructor::
	solveBoostBeta( const double k, const Lorentz5Momentum & newq,
	const Lorentz5Momentum & oldp ) {
	// try something different, purely numerical first:
	// a) boost to rest frame of newq, b) boost with kp/E
	Energy q = newq.vect().mag();
	Energy2 qs = sqr(q);
	Energy2 Q2 = newq.m2();
	Energy kp = k*(oldp.vect().mag());
	Energy2 kps = sqr(kp);

	// usually we take the minus sign, since this boost will be smaller.
	// we only require \|k \vec p\| = \|\vec q'\| which leaves the sign of
	// the boost open but the 'minus' solution gives a smaller boost
	// parameter, i.e. the result should be closest to the previous
	// result. this is to be changed if we would get many momentum
	// conservation violations at the end of the shower from a hard
	// process.
	double betam = (qsqrt(qs + Q2) - kpsqrt(kps + Q2))/(kps + qs + Q2);
	// move directly to 'return'
	Boost beta = -betam(k/kp)oldp.vect();
	// note that (k/kp)*oldp.vect() = oldp.vect()/oldp.vect().mag() but cheaper.
	// leave this out if it's running properly!
	if ( betam >= 0 ) return beta;
	else return Boost(0., 0., 0.);
	}

	bool QTildeReconstructor::
	reconstructDecayJets(ShowerTreePtr decay,
	ShowerInteraction::Type) const {
	_currentTree = decay;
	// extract the particles from the ShowerTree
	vector<ShowerProgenitorPtr> ShowerHardJets=decay->extractProgenitors();
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	_boosts[ShowerHardJets[ix]->progenitor()] = vector<LorentzRotation>();
	}
	for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	tit = _currentTree->treelinks().begin();
	tit != _currentTree->treelinks().end();++tit) {
	_treeBoosts[tit->first] = vector<LorentzRotation>();
	}
	try {
	bool radiated[2]={false,false};
	// find the decaying particle and check if particles radiated
	ShowerProgenitorPtr initial;
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	// only consider initial-state jets
	if(ShowerHardJets[ix]->progenitor()->isFinalState()) {
	radiated[1] \|=ShowerHardJets[ix]->hasEmitted();
	}
	else {
	initial=ShowerHardJets[ix];
	radiated[0]\|=ShowerHardJets[ix]->hasEmitted();
	}
	}
	// find boost to the rest frame if needed
	Boost boosttorest=-initial->progenitor()->momentum().boostVector();
	double gammarest =
	initial->progenitor()->momentum().e()/
	initial->progenitor()->momentum().mass();
	// check if need to boost to rest frame
	bool gottaBoost = (boosttorest.mag() > 1e-12);
	// if initial state radiation reconstruct the jet and set up the basis vectors
	Lorentz5Momentum pjet;
	Lorentz5Momentum nvect;
	// find the partner
	ShowerParticlePtr partner = initial->progenitor()->partner();
	Lorentz5Momentum ppartner[2];
	if(partner) ppartner[0]=partner->momentum();
	// get the n reference vector
	if(partner) {
	if(initial->progenitor()->showerKinematics()) {
	nvect = initial->progenitor()->showerKinematics()->getBasis()[1];
	}
	else {
	Lorentz5Momentum ppartner=initial->progenitor()->partner()->momentum();
	if(gottaBoost) ppartner.boost(boosttorest,gammarest);
	nvect = Lorentz5Momentum( ZERO,0.5initial->progenitor()->mass()
	ppartner.vect().unit());
	nvect.boost(-boosttorest,gammarest);
	}
	}
	// if ISR
	if(radiated[0]) {
	// reconstruct the decay jet
	reconstructDecayJet(initial->progenitor());
	// momentum of decaying particle after ISR
	pjet=initial->progenitor()->momentum()
	-decay->incomingLines().begin()->second->momentum();
	pjet.rescaleMass();
	}
	// boost initial state jet and basis vector if needed
	if(gottaBoost) {
	pjet.boost(boosttorest,gammarest);
	nvect.boost(boosttorest,gammarest);
	ppartner[0].boost(boosttorest,gammarest);
	}
	// loop over the final-state particles and do the reconstruction
	JetKinVect possiblepartners;
	JetKinVect jetKinematics;
	bool atLeastOnce = radiated[0];
	LorentzRotation restboost(boosttorest,gammarest);
	Energy inmass(ZERO);
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	// only consider final-state jets
	if(!ShowerHardJets[ix]->progenitor()->isFinalState()) {
	inmass=ShowerHardJets[ix]->progenitor()->mass();
	continue;
	}
	// do the reconstruction
	JetKinStruct tempJetKin;
	tempJetKin.parent = ShowerHardJets[ix]->progenitor();
	if(ShowerHardJets.size()==2) {
	Lorentz5Momentum dum=ShowerHardJets[ix]->progenitor()->momentum();
	dum.setMass(inmass);
	dum.rescaleRho();
	tempJetKin.parent->set5Momentum(dum);
	}
	tempJetKin.p = ShowerHardJets[ix]->progenitor()->momentum();
	if(gottaBoost) tempJetKin.p.boost(boosttorest,gammarest);
	_progenitor=tempJetKin.parent;
	if(ShowerHardJets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) {
	atLeastOnce \|= reconstructTimeLikeJet(tempJetKin.parent);
	ShowerHardJets[ix]->reconstructed(ShowerProgenitor::done);
	}
	if(gottaBoost) deepTransform(tempJetKin.parent,restboost);
	tempJetKin.q = ShowerHardJets[ix]->progenitor()->momentum();
	jetKinematics.push_back(tempJetKin);
	}
	if(partner) ppartner[1]=partner->momentum();
	// calculate the rescaling parameters
	double k1,k2;
	Lorentz5Momentum qt;
	if(!solveDecayKFactor(initial->progenitor()->mass(),nvect,pjet,
	jetKinematics,partner,ppartner,k1,k2,qt)) {
	for(map<tPPtr,vector<LorentzRotation> >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) {
	for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
	LorentzRotation rot = rit->inverse();
	bit->first->transform(rot);
	}
	}
	_boosts.clear();
	for(map<tShowerTreePtr,vector<LorentzRotation> >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) {
	for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
	LorentzRotation rot = rit->inverse();
	bit->first->transform(rot,false);
	}
	}
	_treeBoosts.clear();
	_currentTree = tShowerTreePtr();
	return false;
	}
	// apply boosts and rescalings to final-state jets
	for(JetKinVect::iterator it = jetKinematics.begin();
	it != jetKinematics.end(); ++it) {
	LorentzRotation Trafo = LorentzRotation();
	if(it->parent!=partner) {
	// boost for rescaling
	if(atLeastOnce) {
	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==it->parent)
	break;
	}
	if(it->parent->children().empty()&&!it->parent->spinInfo() &&
	tit==_currentTree->treelinks().end()) {
	Lorentz5Momentum pnew(k2*it->p.vect(),
	sqrt(sqr(k2*it->p.vect().mag())+it->q.mass2()),
	it->q.mass());
	it->parent->set5Momentum(pnew);
	}
	else {
	// rescaling boost can't ever work in this case
	if(k2<0. && it->q.mass()==ZERO)
	throw KinematicsReconstructionVeto();
	Trafo = solveBoost(k2, it->q, it->p);
	}
	}
	if(gottaBoost) Trafo.boost(-boosttorest,gammarest);
	if(atLeastOnce \|\| gottaBoost) deepTransform(it->parent,Trafo);
	}
	else {
	Lorentz5Momentum pnew=ppartner[0];
	pnew *=k1;
	pnew-=qt;
	pnew.setMass(ppartner[1].mass());
	pnew.rescaleEnergy();
	LorentzRotation Trafo=solveBoost(1.,ppartner[1],pnew);
	if(gottaBoost) Trafo.boost(-boosttorest,gammarest);
	deepTransform(partner,Trafo);
	}
	}
	}
	catch(KinematicsReconstructionVeto) {
	for(map<tPPtr,vector<LorentzRotation> >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) {
	for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
	LorentzRotation rot = rit->inverse();
	bit->first->transform(rot);
	}
	}
	_boosts.clear();
	for(map<tShowerTreePtr,vector<LorentzRotation> >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) {
	for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
	LorentzRotation rot = rit->inverse();
	bit->first->transform(rot,false);
	}
	}
	_treeBoosts.clear();
	_currentTree = tShowerTreePtr();
	return false;
	}
	catch (Exception & ex) {
	_currentTree = tShowerTreePtr();
	_boosts.clear();
	_treeBoosts.clear();
	throw ex;
	}
	_boosts.clear();
	_treeBoosts.clear();
	_currentTree = tShowerTreePtr();
	return true;
	}

	bool QTildeReconstructor::
	reconstructDecayJet( const tShowerParticlePtr p) const {
	if(p->children().empty()) return false;
	tShowerParticlePtr child;
	// if branching reconstruct time-like child
	child = dynamic_ptr_cast<ShowerParticlePtr>(p->children()[1]);
	if(child) {
	_progenitor=child;
	reconstructTimeLikeJet(child);
	// calculate the momentum of the particle
	Lorentz5Momentum pnew=p->momentum()-child->momentum();
	pnew.rescaleMass();
	p->children()[0]->set5Momentum(pnew);
	child=dynamic_ptr_cast<ShowerParticlePtr>(p->children()[0]);
	reconstructDecayJet(child);
	return true;
	}
	return false;
	}

	bool QTildeReconstructor::
	solveDecayKFactor(Energy mb,
	const Lorentz5Momentum & n,
	const Lorentz5Momentum & pjet,
	const JetKinVect & jetKinematics,
	ShowerParticlePtr partner,
	Lorentz5Momentum ppartner[2],
	double & k1, double & k2,
	Lorentz5Momentum & qt) const {
	Energy2 pjn = partner ? pjet.vect()*n.vect() : ZERO;
	Energy2 pcn = partner ? ppartner[0].vect()n.vect() : 1.MeV2;
	Energy2 nmag = n.vect().mag2();
	Lorentz5Momentum pn = partner ? (pjn/nmag)*n : Lorentz5Momentum();
	qt=pjet-pn; qt.setE(ZERO);
	Energy2 pt2=qt.vect().mag2();
	Energy Ejet = pjet.e();
	// magnitudes of the momenta for fast access
	vector<Energy2> pmag;
	Energy total(Ejet);
	for(unsigned int ix=0;ix<jetKinematics.size();++ix) {
	pmag.push_back(jetKinematics[ix].p.vect().mag2());
	total+=jetKinematics[ix].q.mass();
	}
	// return if no possible solution
	if(total>mb) return false;
	Energy2 pcmag=ppartner[0].vect().mag2();
	// used newton-raphson to get the rescaling
	static const Energy eps=1e-8*GeV;
	long double d1(1.),d2(1.);
	Energy roots, ea, ec, ds;
	unsigned int ix=0;
	do {
	++ix;
	d2 = d1 + pjn/pcn;
	roots = Ejet;
	ds = ZERO;
	for(unsigned int iy=0;iy<jetKinematics.size();++iy) {
	if(jetKinematics[iy].parent==partner) continue;
	ea = sqrt(sqr(d2)*pmag[iy]+jetKinematics[iy].q.mass2());
	roots += ea;
	ds += d2/ea*pmag[iy];
	}
	if(partner) {
	ec = sqrt(sqr(d1)*pcmag + pt2 + ppartner[1].mass2());
	roots += ec;
	ds += d1/ec*pcmag;
	}
	d1 += (mb-roots)/ds;
	d2 = d1 + pjn/pcn;
	}
	while(abs(mb-roots)>eps && ix<100);
	k1=d1;
	k2=d2;
	// return true if N-R succeed, otherwise false
	return ix<100;
	}

	bool QTildeReconstructor::
	deconstructDecayJets(HardTreePtr decay, cEvolverPtr,
	ShowerInteraction::Type) const {
	// extract the momenta of the particles
	vector<Lorentz5Momentum> pin;
	vector<Lorentz5Momentum> pout;
	// on-shell masses of the decay products
	vector<Energy> mon;
	Energy mbar(-GeV);
	// the hard branchings of the particles
	set<HardBranchingPtr>::iterator cit;
	set<HardBranchingPtr> branchings=decay->branchings();
	// properties of the incoming particle
	bool ISR = false;
	HardBranchingPtr initial;
	Lorentz5Momentum qisr;
	// find the incoming particle, both before and after
	// any ISR
	for(cit=branchings.begin();cit!=branchings.end();++cit){
	if((*cit)->status()==HardBranching::Incoming\|\|
	(*cit)->status()==HardBranching::Decay) {
	// search back up isr if needed
	HardBranchingPtr branch = *cit;
	while(branch->parent()) branch=branch->parent();
	initial=branch;
	// momentum or original parent
	pin.push_back(branch->branchingParticle()->momentum());
	// ISR?
	ISR = !branch->branchingParticle()->children().empty();
	// ISR momentum
	qisr = pin.back()-(**cit).branchingParticle()->momentum();
	qisr.rescaleMass();
	}
	}
	assert(pin.size()==1);
	// compute boost to rest frame
	Boost boostv=-pin[0].boostVector();
	// partner for ISR
	ShowerParticlePtr partner;
	Lorentz5Momentum ppartner;
	if(initial->branchingParticle()->partner()) {
	partner=initial->branchingParticle()->partner();
	ppartner=partner->momentum();
	}
	// momentum of the decay products
	for(cit=branchings.begin();cit!=branchings.end();++cit) {
	if((*cit)->status()!=HardBranching::Outgoing) continue;
	// find the mass of the particle
	// including special treatment for off-shell resonances
	// to preserve off-shell mass
	Energy mass;
	if(!(**cit).branchingParticle()->dataPtr()->stable()) {
	HardBranchingPtr branch=*cit;
	while(!branch->children().empty()) {
	for(unsigned int ix=0;ix<branch->children().size();++ix) {
	if(branch->children()[ix]->branchingParticle()->id()==
	(**cit).branchingParticle()->id()) {
	branch = branch->children()[ix];
	continue;
	}
	}
	};
	mass = branch->branchingParticle()->mass();
	}
	else {
	mass = (**cit).branchingParticle()->dataPtr()->mass();
	}
	// if not evolution partner of decaying particle
	if((*cit)->branchingParticle()!=partner) {
	pout.push_back((*cit)->branchingParticle()->momentum());
	mon.push_back(mass);
	}
	// evolution partner of decaying particle
	else {
	mbar = mass;
	}
	}
	// boost all the momenta to the rest frame of the decaying particle
	for(unsigned int ix=0;ix<pout.size();++ix) pout[ix].boost(boostv);
	if(initial->branchingParticle()->partner()) {
	ppartner.boost(boostv);
	qisr.boost(boostv);
	}
	// compute the rescaling factors
	double k1,k2;
	if(!ISR) {
	if(partner) {
	pout.push_back(ppartner);
	mon.push_back(mbar);
	}
	k1=k2=inverseRescalingFactor(pout,mon,pin[0].mass());
	if(partner) {
	pout.pop_back();
	mon.pop_back();
	}
	}
	else {
	if(!inverseDecayRescalingFactor(pout,mon,pin[0].mass(),
	ppartner,mbar,k1,k2)) return false;
	}
	// now calculate the p reference vectors
	unsigned int ifinal=0;
	for(cit=branchings.begin();cit!=branchings.end();++cit) {
	if((**cit).status()!=HardBranching::Outgoing) continue;
	// for partners other than colour partner of decaying particle
	if((*cit)->branchingParticle()!=partner) {
	Lorentz5Momentum pvect = (*cit)->branchingParticle()->momentum();
	pvect.boost(boostv);
	pvect /= k1;
	pvect.setMass(mon[ifinal]);
	++ifinal;
	pvect.rescaleEnergy();
	pvect.boost(-boostv);
	(*cit)->pVector(pvect);
	(*cit)->showerMomentum(pvect);
	}
	// for colour partner of decaying particle
	else {
	Lorentz5Momentum pvect = (*cit)->branchingParticle()->momentum();
	pvect.boost(boostv);
	Lorentz5Momentum qtotal;
	for(unsigned int ix=0;ix<pout.size();++ix) qtotal+=pout[ix];
	Lorentz5Momentum qperp =
	qisr-(qisr.vect()qtotal.vect())/(qtotal.vect().mag2())qtotal;
	pvect +=qperp;
	pvect /=k2;
	pvect.setMass(mbar);
	pvect.rescaleEnergy();
	pvect.boost(-boostv);
	(*cit)->pVector(pvect);
	(*cit)->showerMomentum(pvect);
	}
	}
	// // find the evolution partners
	// ShowerParticleVector particles;
	// particles.push_back((**decay->incoming().begin()).branchingParticle());
	// for(cit=branchings.begin();cit!=branchings.end();++cit) {
	// if((**cit).status()==HardBranching::Outgoing)
	// particles.push_back((*cit)->branchingParticle());
	// }
	// // partners should
	// evolver->showerModel()->partnerFinder()
	// ->setInitialEvolutionScales(particles,true,type,false);
	// For initial-state if needed
	if(initial) {
	tShowerParticlePtr newPartner=initial->branchingParticle()->partner();
	if(newPartner) {
	tHardBranchingPtr branch;
	for( set<HardBranchingPtr>::iterator clt = branchings.begin();
	clt != branchings.end(); ++clt ) {
	if((**clt).branchingParticle()==newPartner) {
	initial->colourPartner(*clt);
	branch=*clt;
	break;
	}
	}
	Lorentz5Momentum pvect = initial->branchingParticle()->momentum();
	initial->pVector(pvect);
	Lorentz5Momentum ptemp = branch->pVector();
	ptemp.boost(boostv);
	Lorentz5Momentum nvect = Lorentz5Momentum( ZERO,
	0.5initial->branchingParticle()->mass()
	ptemp.vect().unit());
	nvect.boost(-boostv);
	initial->nVector(nvect);
	}
	}
	// calculate the reference vectors, then for outgoing particles
	for(cit=branchings.begin();cit!=branchings.end();++cit){
	if((**cit).status()!=HardBranching::Outgoing) continue;
	// find the partner branchings
	tShowerParticlePtr newPartner=(*cit)->branchingParticle()->partner();
	if(!newPartner) continue;
	tHardBranchingPtr branch;
	for( set<HardBranchingPtr>::iterator clt = branchings.begin();
	clt != branchings.end(); ++clt ) {
	if(cit==clt) continue;
	if((**clt).branchingParticle()==newPartner) {
	(*cit).colourPartner(clt);
	branch=*clt;
	break;
	}
	}
	if((**decay->incoming().begin()).branchingParticle()==newPartner) {
	(*cit).colourPartner(decay->incoming().begin());
	branch = *decay->incoming().begin();
	}
	// final-state colour partner
	if(branch->status()==HardBranching::Outgoing) {
	Boost boost=((*cit)->pVector()+branch->pVector()).findBoostToCM();
	Lorentz5Momentum pcm = branch->pVector();
	pcm.boost(boost);
	Lorentz5Momentum nvect = Lorentz5Momentum(ZERO,pcm.vect());
	nvect.boost( -boost);
	(*cit)->nVector(nvect);
	}
	// initial-state colour partner
	else {
	Boost boost=branch->pVector().findBoostToCM();
	Lorentz5Momentum pcm = (*cit)->pVector();
	pcm.boost(boost);
	Lorentz5Momentum nvect = Lorentz5Momentum( ZERO, -pcm.vect());
	nvect.boost( -boost);
	(*cit)->nVector(nvect);
	}
	}
	// now compute the new momenta
	// and calculate the shower variables
	for(cit=branchings.begin();cit!=branchings.end();++cit) {
	if((**cit).status()!=HardBranching::Outgoing) continue;
	LorentzRotation B=LorentzRotation(-boostv);
	LorentzRotation A=LorentzRotation(boostv),R;
	if((*cit)->branchingParticle()==partner) {
	Lorentz5Momentum qnew;
	Energy2 dot=(cit)->pVector()(*cit)->nVector();
	double beta = 0.5((cit)->branchingParticle()->momentum().m2()
	-sqr((*cit)->pVector().mass()))/dot;
	qnew=(cit)->pVector()+beta(*cit)->nVector();
	qnew.rescaleMass();
	// compute the boost
	R=BsolveBoost(Aqnew,A(cit)->branchingParticle()->momentum())*A;
	}
	else {
	Lorentz5Momentum qnew;
	if((*cit)->branchingParticle()->partner()) {
	Energy2 dot=(cit)->pVector()(*cit)->nVector();
	double beta = 0.5((cit)->branchingParticle()->momentum().m2()
	-sqr((*cit)->pVector().mass()))/dot;
	qnew=(cit)->pVector()+beta(*cit)->nVector();
	qnew.rescaleMass();
	}
	else {
	qnew = (*cit)->pVector();
	}
	// compute the boost
	R=BsolveBoost(Aqnew,A(cit)->branchingParticle()->momentum())*A;
	}
	// reconstruct the momenta
	(*cit)->setMomenta(R,1.0,Lorentz5Momentum());
	}
	if(initial) {
	initial->setMomenta(LorentzRotation(),1.0,Lorentz5Momentum());
	}
	return true;
	}

	double QTildeReconstructor::
	inverseRescalingFactor(vector<Lorentz5Momentum> pout,
	vector<Energy> mon, Energy roots) const {
	double lambda=1.;
	if(pout.size()==2) {
	double mu_q1(pout[0].m()/roots), mu_q2(pout[1].m()/roots);
	double mu_p1(mon[0]/roots) , mu_p2(mon[1]/roots);
	lambda =
	((1.+mu_q1+mu_q2)(1.-mu_q1-mu_q2)(mu_q1-1.-mu_q2)*(mu_q2-1.-mu_q1))/
	((1.+mu_p1+mu_p2)(1.-mu_p1-mu_p2)(mu_p1-1.-mu_p2)*(mu_p2-1.-mu_p1));
	if(lambda<0.)
	throw Exception() << "Rescaling factor is imaginary in QTildeReconstructor::"
	<< "inverseRescalingFactor lambda^2= " << lambda
	<< Exception::eventerror;
	lambda = sqrt(lambda);
	}
	else {
	unsigned int ntry=0;
	// compute magnitudes once for speed
	vector<Energy2> pmag;
	for(unsigned int ix=0;ix<pout.size();++ix) {
	pmag.push_back(pout[ix].vect().mag2());
	}
	// Newton-Raphson for the rescaling
	vector<Energy> root(pout.size());
	do {
	// compute new energies
	Energy sum(ZERO);
	for(unsigned int ix=0;ix<pout.size();++ix) {
	root[ix] = sqrt(pmag[ix]/sqr(lambda)+sqr(mon[ix]));
	sum+=root[ix];
	}
	// if accuracy reached exit
	if(abs(sum/roots-1.)<1e-10) break;
	// use Newton-Raphson to compute new guess for lambda
	Energy numer(ZERO),denom(ZERO);
	for(unsigned int ix=0;ix<pout.size();++ix) {
	numer +=root[ix];
	denom +=pmag[ix]/root[ix];
	}
	numer-=roots;
	double fact = 1.+sqr(lambda)*numer/denom;
	if(fact<0.) fact=0.5;
	lambda *=fact;
	++ntry;
	}
	while(ntry<100);
	}
	if(isnan(lambda))
	throw Exception() << "Rescaling factor is nan in QTildeReconstructor::"
	<< "inverseRescalingFactor "
	<< Exception::eventerror;
	return lambda;
	}

	bool QTildeReconstructor::
	deconstructGeneralSystem(HardTreePtr tree,
	cEvolverPtr evolver,
	ShowerInteraction::Type type) const {
	// extract incoming and outgoing particles
	ColourSingletShower in,out;
	for(set<HardBranchingPtr>::const_iterator it=tree->branchings().begin();
	it!=tree->branchings().end();++it) {
	if((*it).status()==HardBranching::Incoming) in .jets.push_back(it);
	else out.jets.push_back(*it);
	}
	LorentzRotation toRest,fromRest;
	bool applyBoost(false);
	// do the initial-state reconstruction
	deconstructInitialInitialSystem(applyBoost,toRest,fromRest,
	tree,in.jets,type);
	// do the final-state reconstruction
	deconstructFinalStateSystem(toRest,fromRest,tree,
	out.jets,evolver,type);
	// only at this point that we can be sure all the reference vectors
	// are correct
	for(set<HardBranchingPtr>::const_iterator it=tree->branchings().begin();
	it!=tree->branchings().end();++it) {
	if((**it).status()==HardBranching::Incoming) continue;
	if((**it).branchingParticle()->coloured())
	(**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false);
	}
	for(set<HardBranchingPtr>::const_iterator it=tree->incoming().begin();
	it!=tree->incoming().end();++it) {
	(**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false);
	}
	return true;
	}

	bool QTildeReconstructor::deconstructHardJets(HardTreePtr tree,
	cEvolverPtr evolver,
	ShowerInteraction::Type type) const {
	// inverse of old recon method
	if(_reconopt == 0) {
	return deconstructGeneralSystem(tree,evolver,type);
	}
	else if(_reconopt == 1) {
	return deconstructColourSinglets(tree,evolver,type);
	}
	else if(_reconopt == 2) {
	throw Exception() << "Inverse reconstruction is not currently supported for ReconstructionOption Colour2 "
	<< "in QTildeReconstructor::deconstructHardJets(). Please use one of the other options\n"
	<< Exception::runerror;
	}
	else if(_reconopt == 3 \|\| _reconopt == 4 ) {
	return deconstructColourPartner(tree,evolver,type);
	}
	else
	assert(false);
	}

	bool QTildeReconstructor::
	deconstructColourSinglets(HardTreePtr tree,cEvolverPtr evolver,
	ShowerInteraction::Type type) const {
	// identify the colour singlet systems
	unsigned int nnun(0),nnii(0),nnif(0),nnf(0),nni(0);
	vector<ColourSingletShower>
	systems(identifySystems(tree->branchings(),nnun,nnii,nnif,nnf,nni));
	// now decide what to do
	LorentzRotation toRest,fromRest;
	bool applyBoost(false);
	bool general(false);
	// initial-initial connection and final-state colour singlet systems
	// Drell-Yan type
	if(nnun==0&&nnii==1&&nnif==0&&nnf>0&&nni==0) {
	// reconstruct initial-initial system
	for(unsigned int ix=0;ix<systems.size();++ix) {
	if(systems[ix].type==II)
	deconstructInitialInitialSystem(applyBoost,toRest,fromRest,tree,
	systems[ix].jets,type);
	}
	if(type==ShowerInteraction::QED\|\|type==ShowerInteraction::Both) {
	combineFinalState(systems);
	general=false;
	}
	}
	// DIS and VBF type
	else if(nnun==0&&nnii==0&&((nnif==1&&nnf>0&&nni==1)\|\|
	(nnif==2&& nni==0))) {
	for(unsigned int ix=0;ix<systems.size();++ix) {
	if(systems[ix].type==IF)
	deconstructInitialFinalSystem(tree,systems[ix].jets,evolver,type);
	}
	}
	// e+e- type
	else if(nnun==0&&nnii==0&&nnif==0&&nnf>0&&nni==2) {
	// only FS needed
	// but need to boost to rest frame if QED ISR
	Lorentz5Momentum ptotal;
	for(unsigned int ix=0;ix<systems.size();++ix) {
	if(systems[ix].type==I)
	ptotal += systems[ix].jets[0]->branchingParticle()->momentum();
	}
	toRest = LorentzRotation(ptotal.findBoostToCM());
	fromRest = toRest;
	fromRest.invert();
	if(type==ShowerInteraction::QED\|\|type==ShowerInteraction::Both) {
	combineFinalState(systems);
	general=false;
	}
	}
	// general type
	else {
	general = true;
	}
	// final-state systems except for general recon
	if(!general) {
	for(unsigned int ix=0;ix<systems.size();++ix) {
	if(systems[ix].type==F)
	deconstructFinalStateSystem(toRest,fromRest,tree,
	systems[ix].jets,evolver,type);
	}
	// only at this point that we can be sure all the reference vectors
	// are correct
	for(set<HardBranchingPtr>::const_iterator it=tree->branchings().begin();
	it!=tree->branchings().end();++it) {
	if((**it).status()==HardBranching::Incoming) continue;
	if((**it).branchingParticle()->coloured())
	(**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false);
	}
	for(set<HardBranchingPtr>::const_iterator it=tree->incoming().begin();
	it!=tree->incoming().end();++it) {
	(**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false);
	}
	return true;
	}
	else {
	return deconstructGeneralSystem(tree,evolver,type);
	}
	return true;
	}

	bool QTildeReconstructor::
	deconstructColourPartner(HardTreePtr tree,cEvolverPtr evolver,
	ShowerInteraction::Type type) const {
	Lorentz5Momentum ptotal;
	HardBranchingPtr emitter;
	ColourSingletShower incomingShower,outgoingShower;
	for(set<HardBranchingPtr>::const_iterator it=tree->branchings().begin();
	it!=tree->branchings().end();++it) {
	if((**it).status()==HardBranching::Incoming) {
	incomingShower.jets.push_back(*it);
	ptotal += (*it)->branchingParticle()->momentum();
	// check for emitting particle
	if((**it).parent() ) {
	if(!emitter)
	emitter = *it;
	else
	throw Exception() << "Only one emitting particle allowed in "
	<< "QTildeReconstructor::deconstructColourPartner()"
	<< Exception::runerror;
	}
	}
	else if ((**it).status()==HardBranching::Outgoing) {
	outgoingShower.jets.push_back(*it);
	// check for emitting particle
	if(!(**it).children().empty() ) {
	if(!emitter)
	emitter = *it;
	else
	throw Exception() << "Only one emitting particle allowed in "
	<< "QTildeReconstructor::deconstructColourPartner()"
	<< Exception::runerror;
	}
	}
	}
	assert(emitter);
	assert(emitter->colourPartner());
	ColourSingletShower system;
	system.jets.push_back(emitter);
	system.jets.push_back(emitter->colourPartner());
	LorentzRotation toRest,fromRest;
	bool applyBoost(false);
	// identify the colour singlet system
	if(emitter->status() == HardBranching::Outgoing &&
	emitter->colourPartner()->status() == HardBranching::Outgoing ) {
	system.type=F;
	// need to boost to rest frame if QED ISR
	if( !incomingShower.jets[0]->branchingParticle()->coloured() &&
	!incomingShower.jets[1]->branchingParticle()->coloured() ) {
	Boost boost = ptotal.findBoostToCM();
	toRest = LorentzRotation( boost);
	fromRest = LorentzRotation(-boost);
	}
	else
	findInitialBoost(ptotal,ptotal,toRest,fromRest);
	deconstructFinalStateSystem(toRest,fromRest,tree,
	system.jets,evolver,type);
	}
	else if (emitter->status() == HardBranching::Incoming &&
	emitter->colourPartner()->status() == HardBranching::Incoming) {
	system.type=II;
	deconstructInitialInitialSystem(applyBoost,toRest,fromRest,tree,system.jets,type);
	// make sure the recoil gets applied
	deconstructFinalStateSystem(toRest,fromRest,tree,
	outgoingShower.jets,evolver,type);
	}
	else if ((emitter->status() == HardBranching::Outgoing &&
	emitter->colourPartner()->status() == HardBranching::Incoming ) \|\|
	(emitter->status() == HardBranching::Incoming &&
	emitter->colourPartner()->status() == HardBranching::Outgoing)) {
	system.type=IF;
	// enusre incoming first
	if(system.jets[0]->status() == HardBranching::Outgoing)
	swap(system.jets[0],system.jets[1]);
	deconstructInitialFinalSystem(tree,system.jets,evolver,type);
	}
	else {
	throw Exception() << "Unknown type of system in "
	<< "QTildeReconstructor::deconstructColourPartner()"
	<< Exception::runerror;
	}
	// only at this point that we can be sure all the reference vectors
	// are correct
	for(set<HardBranchingPtr>::const_iterator it=tree->branchings().begin();
	it!=tree->branchings().end();++it) {
	if((**it).status()==HardBranching::Incoming) continue;
	if((**it).branchingParticle()->coloured())
	(**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false);
	}
	for(set<HardBranchingPtr>::const_iterator it=tree->incoming().begin();
	it!=tree->incoming().end();++it) {
	(**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false);
	}
	for(set<HardBranchingPtr>::const_iterator it=tree->branchings().begin();
	it!=tree->branchings().end();++it) {
	if((**it).status()!=HardBranching::Incoming) continue;
	if(it==system.jets[0] \|\| it==system.jets[1]) continue;
	if((**it).branchingParticle()->momentum().z()>ZERO) {
	(it).z((it).branchingParticle()->momentum().plus()/(**it).beam()->momentum().plus());
	}
	else {
	(it).z((it).branchingParticle()->momentum().minus()/(**it).beam()->momentum().minus());
	}
	}
	return true;
	}

	void QTildeReconstructor::
	reconstructInitialFinalSystem(vector<ShowerProgenitorPtr> jets) const {
	Lorentz5Momentum pin[2],pout[2],pbeam;
	for(unsigned int ix=0;ix<jets.size();++ix) {
	// final-state parton
	if(jets[ix]->progenitor()->isFinalState()) {
	pout[0] +=jets[ix]->progenitor()->momentum();
	_progenitor = jets[ix]->progenitor();
	if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) {
	reconstructTimeLikeJet(jets[ix]->progenitor());
	jets[ix]->reconstructed(ShowerProgenitor::done);
	}
	}
	// initial-state parton
	else {
	pin[0] +=jets[ix]->progenitor()->momentum();
	if(jets[ix]->progenitor()->showerKinematics()) {
	pbeam = jets[ix]->progenitor()->showerKinematics()->getBasis()[0];
	}
	else {
	if ( jets[ix]->original()->parents().empty() ) {
	pbeam = jets[ix]->progenitor()->momentum();
	}
	else {
	pbeam = jets[ix]->original()->parents()[0]->momentum();
	}
	}
	if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) {
	reconstructSpaceLikeJet(jets[ix]->progenitor());
	jets[ix]->reconstructed(ShowerProgenitor::done);
	}
	assert(!jets[ix]->original()->parents().empty());
	}
	}
	// add intrinsic pt if needed
	addIntrinsicPt(jets);
	// momenta after showering
	for(unsigned int ix=0;ix<jets.size();++ix) {
	if(jets[ix]->progenitor()->isFinalState())
	pout[1] += jets[ix]->progenitor()->momentum();
	else
	pin[1] += jets[ix]->progenitor()->momentum();
	}
	// work out the boost to the Breit frame
	Lorentz5Momentum pa = pout[0]-pin[0];
	Axis axis(pa.vect().unit());
	LorentzRotation rot;
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	if ( sinth > 1.e-9 )
	rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	rot.rotateX(Constants::pi);
	rot.boostZ( pa.e()/pa.vect().mag());
	Lorentz5Momentum ptemp=rot*pbeam;
	Boost trans = -1./ptemp.e()*ptemp.vect();
	trans.setZ(0.);
	if ( trans.mag2() - 1. >= 0. ) throw KinematicsReconstructionVeto();
	rot.boost(trans);
	pa *=rot;
	// project and calculate rescaling
	// reference vectors
	Lorentz5Momentum n1(ZERO,ZERO,-pa.z(),-pa.z());
	Lorentz5Momentum n2(ZERO,ZERO, pa.z(),-pa.z());
	Energy2 n1n2 = n1*n2;
	// decompose the momenta
	Lorentz5Momentum qbp=rotpin[1],qcp=rotpout[1];
	qbp.rescaleMass();
	qcp.rescaleMass();
	double a[2],b[2];
	a[0] = n2*qbp/n1n2;
	b[0] = n1*qbp/n1n2;
	Lorentz5Momentum qperp = qbp-a[0]n1-b[0]n2;
	b[1] = 0.5;
	a[1] = 0.5*(qcp.m2()-qperp.m2())/n1n2/b[1];
	double kb;
	if(a[0]!=0.) {
	double A(0.5a[0]),B(b[0]a[0]-a[1]b[1]-0.25),C(-0.5b[0]);
	if(sqr(B)-4.AC<0.) throw KinematicsReconstructionVeto();
	kb = 0.5(-B+sqrt(sqr(B)-4.A*C))/A;
	}
	else {
	kb = 0.5b[0]/(b[0]a[0]-a[1]*b[1]-0.25);
	}
	// changed to improve stability
	if(kb==0.) throw KinematicsReconstructionVeto();
	if ( a[1]>b[1] && abs(a[1]) < 1e-12 )
	throw KinematicsReconstructionVeto();
	if ( a[1]<=b[1] && abs(0.5+b[0]/kb) < 1e-12 )
	throw KinematicsReconstructionVeto();
	double kc = (a[1]>b[1]) ? (a[0]*kb-0.5)/a[1] : b[1]/(0.5+b[0]/kb);
	if(kc==0.) throw KinematicsReconstructionVeto();
	Lorentz5Momentum pnew[2] = { a[0]kbn1+b[0]/kb*n2+qperp,
	a[1]kcn1+b[1]/kc*n2+qperp};
	LorentzRotation rotinv=rot.inverse();
	for(unsigned int ix=0;ix<jets.size();++ix) {
	if(jets[ix]->progenitor()->isFinalState()) {
	deepTransform(jets[ix]->progenitor(),rot);
	deepTransform(jets[ix]->progenitor(),solveBoost(pnew[1],qcp));
	Energy delta = jets[ix]->progenitor()->momentum().m()-jets[ix]->progenitor()->momentum().mass();
	if ( abs(delta) > MeV ) throw KinematicsReconstructionVeto();
	deepTransform(jets[ix]->progenitor(),rotinv);
	}
	else {
	tPPtr parent;
	boostChain(jets[ix]->progenitor(),rot,parent);
	boostChain(jets[ix]->progenitor(),solveBoostZ(pnew[0],qbp),parent);
	// check the first boost worked, and if not apply small correction to
	// fix energy/momentum conservation
	// this is a kludge but it reduces momentum non-conservation dramatically
	Lorentz5Momentum pdiff = pnew[0]-jets[ix]->progenitor()->momentum();
	Energy2 delta = sqr(pdiff.x())+sqr(pdiff.y())+sqr(pdiff.z())+sqr(pdiff.t());
	unsigned int ntry=0;
	while(delta>1e-6*GeV2 && ntry<5 ) {
	ntry +=1;
	boostChain(jets[ix]->progenitor(),solveBoostZ(pnew[0],jets[ix]->progenitor()->momentum()),parent);
	pdiff = pnew[0]-jets[ix]->progenitor()->momentum();
	delta = sqr(pdiff.x())+sqr(pdiff.y())+sqr(pdiff.z())+sqr(pdiff.t());
	}
	// apply test in breit-frame
	Lorentz5Momentum ptest1 = parent->momentum();
	Lorentz5Momentum ptest2 = rot*pbeam;
	if(ptest1.z()/ptest2.z()<0. \|\| ptest1.z()/ptest2.z()>1.)
	throw KinematicsReconstructionVeto();
	boostChain(jets[ix]->progenitor(),rotinv,parent);
	}
	}
	}

	bool QTildeReconstructor::addIntrinsicPt(vector<ShowerProgenitorPtr> jets) const {
	bool added=false;
	// add the intrinsic pt if needed
	for(unsigned int ix=0;ix<jets.size();++ix) {
	// only for initial-state particles which haven't radiated
	if(jets[ix]->progenitor()->isFinalState()\|\|
	jets[ix]->hasEmitted()\|\|
	jets[ix]->reconstructed()==ShowerProgenitor::dontReconstruct) continue;
	if(_intrinsic.find(jets[ix])==_intrinsic.end()) continue;
	pair<Energy,double> pt=_intrinsic[jets[ix]];
	Energy etemp = jets[ix]->original()->parents()[0]->momentum().z();
	Lorentz5Momentum
	p_basis(ZERO, ZERO, etemp, abs(etemp)),
	n_basis(ZERO, ZERO,-etemp, abs(etemp));
	double alpha = jets[ix]->progenitor()->x();
	double beta = 0.5*(sqr(jets[ix]->progenitor()->data().mass())+
	sqr(pt.first))/alpha/(p_basis*n_basis);
	Lorentz5Momentum pnew=alphap_basis+betan_basis;
	pnew.setX(pt.first*cos(pt.second));
	pnew.setY(pt.first*sin(pt.second));
	pnew.rescaleMass();
	jets[ix]->progenitor()->set5Momentum(pnew);
	added = true;
	}
	return added;
	}

	LorentzRotation QTildeReconstructor::
	solveBoost(const double k, const Lorentz5Momentum & newq,
	const Lorentz5Momentum & oldp ) const {
	Energy q = newq.vect().mag();
	Energy2 qs = sqr(q);
	Energy2 Q2 = newq.mass2();
	Energy kp = k*(oldp.vect().mag());
	Energy2 kps = sqr(kp);
	double betam = (qnewq.e() - kpsqrt(kps + Q2))/(kps + qs + Q2);
	if ( abs(betam) - 1. >= 0. ) throw KinematicsReconstructionVeto();
	Boost beta = -betam(k/kp)oldp.vect();
	double gamma = 0.;
	if(Q2/sqr(oldp.e())>1e-4) {
	if(betam<0.5) {
	gamma = 1./sqrt(1.-sqr(betam));
	}
	else {
	gamma = ( kps+ qs + Q2)/
	sqrt(2.kpsqs + kpsQ2 + qsQ2 + sqr(Q2) + 2.qnewq.e()kpsqrt(kps + Q2));
	}
	}
	else {
	if(k>0) {
	gamma = 4.kpsqs/sqr(kps +qs) + 2.sqr(kps-qs)Q2/pow<3,1>(kps +qs)
	- 0.25( sqr(kps) + 14.kpsqs + sqr(qs))sqr(kps-qs)/(pow<4,1>(kps +qs)kpsqs)*sqr(Q2);
	}
	else {
	gamma = 0.25sqr(Q2)/(kpsqs)(1. - 0.5(kps+qs)/(kpsqs)Q2);
	}
	gamma = 1./sqrt(gamma);
	}
	// note that (k/kp)*oldp.vect() = oldp.vect()/oldp.vect().mag() but cheaper.
	ThreeVector<Energy2> ax = newq.vect().cross( oldp.vect() );
	double delta = newq.vect().angle( oldp.vect() );
	LorentzRotation R;
	using Constants::pi;
	Energy2 scale1 = sqr(newq.x())+ sqr(newq.y())+sqr(newq.z());
	Energy2 scale2 = sqr(oldp.x())+ sqr(oldp.y())+sqr(oldp.z());
	if ( ax.mag2()/scale1/scale2 > 1e-28 ) {
	R.rotate( delta, unitVector(ax) ).boost( beta , gamma );
	}
	else if(abs(delta-pi)/pi < 0.001) {
	double phi=2.piUseRandom::rnd();
	Axis axis(cos(phi),sin(phi),0.);
	axis.rotateUz(newq.vect().unit());
	R.rotate(delta,axis).boost( beta , gamma );
	}
	else {
	R.boost( beta , gamma );
	}
	return R;
	}

	LorentzRotation QTildeReconstructor::solveBoost(const Lorentz5Momentum & q,
	const Lorentz5Momentum & p ) const {
	Energy modp = p.vect().mag();
	Energy modq = q.vect().mag();
	double betam = (p.e()modp-q.e()modq)/(sqr(modq)+sqr(modp)+p.mass2());
	if ( abs(betam)-1. >= 0. ) throw KinematicsReconstructionVeto();
	Boost beta = -betam*q.vect().unit();
	ThreeVector<Energy2> ax = p.vect().cross( q.vect() );
	double delta = p.vect().angle( q.vect() );
	LorentzRotation R;
	using Constants::pi;
	if ( beta.mag2() - 1. >= 0. ) throw KinematicsReconstructionVeto();
	if ( ax.mag2()/GeV2/MeV2 > 1e-16 ) {
	R.rotate( delta, unitVector(ax) ).boost( beta );
	}
	else {
	R.boost( beta );
	}
	return R;
	}

	LorentzRotation QTildeReconstructor::solveBoostZ(const Lorentz5Momentum & q,
	const Lorentz5Momentum & p ) const {
	static const double eps = 1e-6;
	LorentzRotation R;
	double beta;
	Energy2 mt2 = p.mass()<ZERO ? -sqr(p.mass())+sqr(p.x())+sqr(p.y()) : sqr(p.mass())+sqr(p.x())+sqr(p.y()) ;
	double ratio = mt2/(sqr(p.t())+sqr(q.t()));
	if(abs(ratio)>eps) {
	double erat = (q.t()+q.z())/(p.t()+p.z());
	Energy2 den = mt2*(erat+1./erat);
	Energy2 num = (q.z()-p.z())(q.t()+p.t()) + (p.z()+q.z())(p.t()-q.t());
	beta = num/den;
	if ( abs(beta) - 1. >= 0. ) throw KinematicsReconstructionVeto();
	R.boostZ(beta);
	}
	else {
	double er = sqr(p.t()/q.t());
	double x = ratio+0.125(er+10.+1./er)sqr(ratio);
	beta = -(p.t()-q.t())(p.t()+q.t())/(sqr(p.t())+sqr(q.t()))(1.+x);
	double gamma = (4.sqr(p.t()q.t()) +sqr(p.t()-q.t())sqr(p.t()+q.t())
	(-2.*x+sqr(x)))/sqr(sqr(p.t())+sqr(q.t()));
	if ( abs(beta) - 1. >= 0. ) throw KinematicsReconstructionVeto();
	gamma = 1./sqrt(gamma);
	R.boost(0.,0.,beta,gamma);
	}
	Lorentz5Momentum ptest = R*p;
	if(ptest.z()/q.z() < 0. \|\| ptest.t()/q.t() < 0. ) {
	throw KinematicsReconstructionVeto();
	}
	return R;
	}

	void QTildeReconstructor::
	reconstructFinalStateSystem(bool applyBoost,
	const LorentzRotation & toRest,
	const LorentzRotation & fromRest,
	vector<ShowerProgenitorPtr> jets) const {
	LorentzRotation trans = applyBoost? toRest : LorentzRotation();
	// special for case of individual particle
	if(jets.size()==1) {
	deepTransform(jets[0]->progenitor(),trans);
	deepTransform(jets[0]->progenitor(),fromRest);
	return;
	}
	bool radiated(false);
	// find the hard process centre-of-mass energy
	Lorentz5Momentum pcm;
	// check if radiated and calculate total momentum
	for(unsigned int ix=0;ix<jets.size();++ix) {
	radiated \|=jets[ix]->hasEmitted();
	pcm += jets[ix]->progenitor()->momentum();
	}
	if(applyBoost) pcm *= trans;
	// check if in CMF frame
	Boost beta_cm = pcm.findBoostToCM();
	bool gottaBoost(false);
	if(beta_cm.mag() > 1e-12) {
	gottaBoost = true;
	trans.boost(beta_cm);
	}
	// collection of pointers to initial hard particle and jet momenta
	// for final boosts
	JetKinVect jetKinematics;
	vector<ShowerProgenitorPtr>::const_iterator cit;
	for(cit = jets.begin(); cit != jets.end(); cit++) {
	JetKinStruct tempJetKin;
	tempJetKin.parent = (*cit)->progenitor();
	if(applyBoost \|\| gottaBoost) {
	deepTransform(tempJetKin.parent,trans);
	}
	tempJetKin.p = (*cit)->progenitor()->momentum();
	_progenitor=tempJetKin.parent;
	if((**cit).reconstructed()==ShowerProgenitor::notReconstructed) {
	radiated \|= reconstructTimeLikeJet((*cit)->progenitor());
	(**cit).reconstructed(ShowerProgenitor::done);
	}
	else {
	radiated \|= !(*cit)->progenitor()->children().empty();
	}
	tempJetKin.q = (*cit)->progenitor()->momentum();
	jetKinematics.push_back(tempJetKin);
	}
	// default option rescale everything with the same factor
	if( _finalStateReconOption == 0 \|\| jetKinematics.size() <= 2 ) {
	// find the rescaling factor
	double k = 0.0;
	if(radiated) {
	k = solveKfactor(pcm.m(), jetKinematics);
	// perform the rescaling and boosts
	for(JetKinVect::iterator it = jetKinematics.begin();
	it != jetKinematics.end(); ++it) {
	LorentzRotation Trafo = solveBoost(k, it->q, it->p);
	deepTransform(it->parent,Trafo);
	}
	}
	}
	// different treatment of most off-shell
	else if ( _finalStateReconOption <= 4 ) {
	// sort the jets by virtuality
	std::sort(jetKinematics.begin(),jetKinematics.end(),JetOrdering());
	// Bryan's procedures from FORTRAN
	if( _finalStateReconOption <=2 ) {
	// loop over the off-shell partons, _finalStateReconOption==1 only first ==2 all
	JetKinVect::const_iterator jend = _finalStateReconOption==1 ? jetKinematics.begin()+1 : jetKinematics.end();
	for(JetKinVect::const_iterator jit=jetKinematics.begin(); jit!=jend;++jit) {
	// calculate the 4-momentum of the recoiling system
	Lorentz5Momentum psum;
	bool done = true;
	for(JetKinVect::const_iterator it=jetKinematics.begin();it!=jetKinematics.end();++it) {
	if(it==jit) {
	done = false;
	continue;
	}
	// first option put on-shell and sum 4-momenta
	if( _finalStateReconOption == 1 ) {
	LorentzRotation Trafo = solveBoost(1., it->q, it->p);
	deepTransform(it->parent,Trafo);
	psum += it->parent->momentum();
	}
	// second option, sum momenta
	else {
	// already rescaled
	if(done) psum += it->parent->momentum();
	// still needs to be rescaled
	else psum += it->p;
	}
	}
	// set the mass
	psum.rescaleMass();
	// calculate the 3-momentum rescaling factor
	Energy2 s(pcm.m2());
	Energy2 m1sq(jit->q.m2()),m2sq(psum.m2());
	Energy4 num = sqr(s - m1sq - m2sq) - 4.m1sqm2sq;
	if(num<ZERO) throw KinematicsReconstructionVeto();
	double k = sqrt( num / (4.sjit->p.vect().mag2()) );
	// boost the off-shell parton
	LorentzRotation B1 = solveBoost(k, jit->q, jit->p);
	deepTransform(jit->parent,B1);
	// boost everything else to rescale
	LorentzRotation B2 = solveBoost(k, psum, psum);
	for(JetKinVect::iterator it=jetKinematics.begin();it!=jetKinematics.end();++it) {
	if(it==jit) continue;
	deepTransform(it->parent,B2);
	it->p *= B2;
	it->q *= B2;
	}
	}
	}
	// Peter's C++ procedures
	else {
	reconstructFinalFinalOffShell(jetKinematics,pcm.m2(), _finalStateReconOption == 4);
	}
	}
	else
	assert(false);
	// apply the final boosts
	if(gottaBoost \|\| applyBoost) {
	LorentzRotation finalBoosts;
	if(gottaBoost) finalBoosts.boost(-beta_cm);
	if(applyBoost) finalBoosts.transform(fromRest);
	for(JetKinVect::iterator it = jetKinematics.begin();
	it != jetKinematics.end(); ++it) {
	deepTransform(it->parent,finalBoosts);
	}
	}
	}

	void QTildeReconstructor::
	reconstructInitialInitialSystem(bool & applyBoost,
	LorentzRotation & toRest,
	LorentzRotation & fromRest,
	vector<ShowerProgenitorPtr> jets) const {
	bool radiated = false;
	Lorentz5Momentum pcm;
	// check whether particles radiated and calculate total momentum
	for( unsigned int ix = 0; ix < jets.size(); ++ix ) {
	radiated \|= jets[ix]->hasEmitted();
	pcm += jets[ix]->progenitor()->momentum();
	if(jets[ix]->original()->parents().empty()) return;
	}
	pcm.rescaleMass();
	// check if intrinsic pt to be added
	radiated \|= !_intrinsic.empty();
	// if no radiation return
	if(!radiated) {
	for(unsigned int ix=0;ix<jets.size();++ix) {
	if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed)
	jets[ix]->reconstructed(ShowerProgenitor::done);
	}
	return;
	}
	// initial state shuffling
	applyBoost=false;
	vector<Lorentz5Momentum> p, pq, p_in;
	vector<Energy> pts;
	for(unsigned int ix=0;ix<jets.size();++ix) {
	// add momentum to vector
	p_in.push_back(jets[ix]->progenitor()->momentum());
	// reconstruct the jet
	if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) {
	radiated \|= reconstructSpaceLikeJet(jets[ix]->progenitor());
	jets[ix]->reconstructed(ShowerProgenitor::done);
	}
	assert(!jets[ix]->original()->parents().empty());
	Energy etemp = jets[ix]->original()->parents()[0]->momentum().z();
	Lorentz5Momentum ptemp = Lorentz5Momentum(ZERO, ZERO, etemp, abs(etemp));
	pq.push_back(ptemp);
	pts.push_back(jets[ix]->highestpT());
	}
	// add the intrinsic pt if needed
	radiated \|=addIntrinsicPt(jets);
	for(unsigned int ix=0;ix<jets.size();++ix) {
	p.push_back(jets[ix]->progenitor()->momentum());
	}
	double x1 = p_in[0].z()/pq[0].z();
	double x2 = p_in[1].z()/pq[1].z();
	vector<double> beta=initialStateRescaling(x1,x2,p_in[0]+p_in[1],p,pq,pts);
	// if not need don't apply boosts
	if(!(radiated && p.size() == 2 && pq.size() == 2)) return;
	applyBoost=true;
	// apply the boosts
	Lorentz5Momentum newcmf;
	for(unsigned int ix=0;ix<jets.size();++ix) {
	tPPtr toBoost = jets[ix]->progenitor();
	Boost betaboost(0, 0, beta[ix]);
	tPPtr parent;
	boostChain(toBoost, LorentzRotation(0.,0.,beta[ix]),parent);
	if(parent->momentum().e()/pq[ix].e()>1.\|\|
	parent->momentum().z()/pq[ix].z()>1.) throw KinematicsReconstructionVeto();
	newcmf+=toBoost->momentum();
	}
	if(newcmf.m()<ZERO\|\|newcmf.e()<ZERO) throw KinematicsReconstructionVeto();
	findInitialBoost(pcm,newcmf,toRest,fromRest);
	}

	void QTildeReconstructor::
	deconstructInitialInitialSystem(bool & applyBoost,
	LorentzRotation & toRest,
	LorentzRotation & fromRest,
	HardTreePtr tree,
	vector<HardBranchingPtr> jets,
	ShowerInteraction::Type) const {
	// get the momenta of the particles
	vector<Lorentz5Momentum> pin;
	vector<Lorentz5Momentum> pq;
	vector<HardBranchingPtr>::iterator cit;
	for(cit=jets.begin();cit!=jets.end();++cit) {
	pin.push_back((*cit)->branchingParticle()->momentum());
	Energy etemp = (*cit)->beam()->momentum().z();
	pq.push_back(Lorentz5Momentum(ZERO, ZERO,etemp, abs(etemp)));
	}
	bool order = (*tree->incoming().begin())->beam()->momentum().z()/pq[0].z()<0.;
	assert(pin.size()==2);
	// calculate the x values
	double x[2];
	Lorentz5Momentum pcm=pin[0]+pin[1];
	assert(pcm.mass2()>ZERO);
	pcm.rescaleMass();
	vector<double> boost = inverseInitialStateRescaling(x[0],x[1],pcm,pin,pq);
	// apply the boost the the particles
	// first incoming particle
	if(order) swap(pq[0],pq[1]);
	// now apply the boosts
	Boost betaboost(0.,0.,-boost[0]);
	LorentzRotation R;
	R.boost(betaboost);
	set<HardBranchingPtr>::const_iterator cjt=tree->incoming().begin();
	(*cjt)->pVector(pq[0]);
	(*cjt)->nVector(pq[1]);
	(*cjt)->setMomenta(R,1.,Lorentz5Momentum());
	// second incoming particle
	betaboost = Boost(0.,0.,-boost[1]);
	R=LorentzRotation(betaboost);
	++cjt;
	(*cjt)->pVector(pq[1]);
	(*cjt)->nVector(pq[0]);
	(*cjt)->setMomenta(R,1.,Lorentz5Momentum());
	jets[0]->showerMomentum(x[0]*jets[0]->pVector());
	jets[1]->showerMomentum(x[1]*jets[1]->pVector());
	// and calculate the boosts
	applyBoost=true;
	// do one boost
	if(_initialBoost==0) {
	toRest = LorentzRotation(-pcm.boostVector());
	}
	else if(_initialBoost==1) {
	// first the transverse boost
	Energy pT = sqrt(sqr(pcm.x())+sqr(pcm.y()));
	double beta = -pT/pcm.t();
	toRest=LorentzRotation(Boost(betapcm.x()/pT,betapcm.y()/pT,0.));
	// the longitudinal
	beta = pcm.z()/sqrt(pcm.m2()+sqr(pcm.z()));
	toRest.boost(Boost(0.,0.,-beta));
	}
	else
	assert(false);
	fromRest = LorentzRotation((jets[0]->showerMomentum()+
	jets[1]->showerMomentum()).boostVector());
	}

	void QTildeReconstructor::
	deconstructFinalStateSystem(const LorentzRotation & toRest,
	const LorentzRotation & fromRest,
	HardTreePtr tree, vector<HardBranchingPtr> jets,
	cEvolverPtr evolver,
	ShowerInteraction::Type type) const {
	LorentzRotation trans = toRest;
	if(jets.size()==1) {
	Lorentz5Momentum pnew = toRest*(jets[0]->branchingParticle()->momentum());
	pnew *= fromRest;
	jets[0]-> original(pnew);
	jets[0]->showerMomentum(pnew);
	// find the colour partners
	ShowerParticleVector particles;
	vector<Lorentz5Momentum> ptemp;
	set<HardBranchingPtr>::const_iterator cjt;
	for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) {
	ptemp.push_back((**cjt).branchingParticle()->momentum());
	(cjt).branchingParticle()->set5Momentum((cjt).showerMomentum());
	particles.push_back((**cjt).branchingParticle());
	}
	evolver->showerModel()->partnerFinder()
	->setInitialEvolutionScales(particles,false,type,false);
	// calculate the reference vectors
	unsigned int iloc(0);
	set<HardBranchingPtr>::iterator clt;
	for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) {
	// reset the momentum
	(**cjt).branchingParticle()->set5Momentum(ptemp[iloc]);
	++iloc;
	// sort out the partners
	tShowerParticlePtr partner =
	(*cjt)->branchingParticle()->partner();
	if(!partner) continue;
	for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) {
	if((**clt).branchingParticle()==partner) {
	(*cjt).colourPartner(clt);
	break;
	}
	}
	tHardBranchingPtr branch;
	for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) {
	if(clt==cjt) continue;
	if((*clt)->branchingParticle()==partner) {
	branch=*clt;
	break;
	}
	}
	}
	return;
	}
	vector<HardBranchingPtr>::iterator cit;
	vector<Lorentz5Momentum> pout;
	vector<Energy> mon;
	Lorentz5Momentum pin;
	for(cit=jets.begin();cit!=jets.end();++cit) {
	pout.push_back((*cit)->branchingParticle()->momentum());
	mon.push_back(findMass(*cit));
	pin+=pout.back();
	}
	// boost all the momenta to the rest frame of the decaying particle
	pin.rescaleMass();
	pin *=trans;
	Boost beta_cm = pin.findBoostToCM();
	bool gottaBoost(false);
	if(beta_cm.mag() > 1e-12) {
	gottaBoost = true;
	trans.boost(beta_cm);
	pin.boost(beta_cm);
	}
	for(unsigned int ix=0;ix<pout.size();++ix) {
	pout[ix].transform(trans);
	}
	// rescaling factor
	double lambda=inverseRescalingFactor(pout,mon,pin.mass());
	if (lambda< 1.e-10) throw KinematicsReconstructionVeto();
	// now calculate the p reference vectors
	for(unsigned int ix=0;ix<jets.size();++ix) {
	Lorentz5Momentum pvect = jets[ix]->branchingParticle()->momentum();
	pvect.transform(trans);
	pvect /= lambda;
	pvect.setMass(mon[ix]);
	pvect.rescaleEnergy();
	if(gottaBoost) pvect.boost(-beta_cm);
	pvect.transform(fromRest);
	jets[ix]->pVector(pvect);
	jets[ix]->showerMomentum(pvect);
	}
	// find the colour partners
	ShowerParticleVector particles;
	vector<Lorentz5Momentum> ptemp;
	set<HardBranchingPtr>::const_iterator cjt;
	for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) {
	ptemp.push_back((**cjt).branchingParticle()->momentum());
	(cjt).branchingParticle()->set5Momentum((cjt).showerMomentum());
	particles.push_back((**cjt).branchingParticle());
	}
	evolver->showerModel()->partnerFinder()
	->setInitialEvolutionScales(particles,false,type,false);
	// calculate the reference vectors
	unsigned int iloc(0);
	set<HardBranchingPtr>::iterator clt;
	for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) {
	// reset the momentum
	(**cjt).branchingParticle()->set5Momentum(ptemp[iloc]);
	++iloc;
	}
	for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) {
	// sort out the partners
	tShowerParticlePtr partner =
	(*cjt)->branchingParticle()->partner();
	if(!partner) continue;
	for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) {
	if((**clt).branchingParticle()==partner) {
	(*cjt).colourPartner(clt);
	break;
	}
	}
	tHardBranchingPtr branch;
	for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) {
	if(clt==cjt) continue;
	if((*clt)->branchingParticle()==partner) {
	branch=*clt;
	break;
	}
	}
	// compute the reference vectors
	// both incoming, should all ready be done
	if((**cjt).status()==HardBranching::Incoming &&
	(**clt).status()==HardBranching::Incoming) {
	continue;
	}
	// both outgoing
	else if((**cjt).status()!=HardBranching::Incoming&&
	branch->status()==HardBranching::Outgoing) {
	Boost boost=((*cjt)->pVector()+branch->pVector()).findBoostToCM();
	Lorentz5Momentum pcm = branch->pVector();
	pcm.boost(boost);
	Lorentz5Momentum nvect = Lorentz5Momentum(ZERO,pcm.vect());
	nvect.boost( -boost);
	(**cjt).nVector(nvect);
	}
	else if((**cjt).status()==HardBranching::Incoming) {
	Lorentz5Momentum pa = -(**cjt).showerMomentum()+branch->showerMomentum();
	Lorentz5Momentum pb = (**cjt).showerMomentum();
	Axis axis(pa.vect().unit());
	LorentzRotation rot;
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	rot.rotateX(Constants::pi);
	rot.boostZ( pa.e()/pa.vect().mag());
	pb*=rot;
	Boost trans = -1./pb.e()*pb.vect();
	trans.setZ(0.);
	rot.boost(trans);
	Energy scale=(**cjt).beam()->momentum().e();
	Lorentz5Momentum pbasis(ZERO,(*cjt).beam()->momentum().vect().unit()scale);
	Lorentz5Momentum pcm = rot*pbasis;
	rot.invert();
	(*cjt).nVector(rotLorentz5Momentum(ZERO,-pcm.vect()));
	tHardBranchingPtr branch2 = *cjt;;
	while (branch2->parent()) {
	branch2=branch2->parent();
	branch2->nVector(rot*Lorentz5Momentum(ZERO,-pcm.vect()));
	}
	}
	else if(branch->status()==HardBranching::Incoming) {
	(**cjt).nVector(Lorentz5Momentum(ZERO,branch->showerMomentum().vect()));
	}
	}
	// now compute the new momenta
	for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) {
	if(!(*cjt)->branchingParticle()->isFinalState()) continue;
	Lorentz5Momentum qnew;
	if((*cjt)->branchingParticle()->partner()) {
	Energy2 dot=(cjt)->pVector()(*cjt)->nVector();
	double beta = 0.5((cjt)->branchingParticle()->momentum().m2()
	-sqr((*cjt)->pVector().mass()))/dot;
	qnew=(cjt)->pVector()+beta(*cjt)->nVector();
	qnew.rescaleMass();
	}
	else {
	qnew = (*cjt)->pVector();
	}
	// qnew is the unshuffled momentum in the rest frame of the p basis vectors,
	// for the simple case Z->q qbar g this was checked against analytic formulae.
	// compute the boost
	LorentzRotation R=solveBoost(qnew,
	toRest(cjt)->branchingParticle()->momentum())*toRest;
	(*cjt)->setMomenta(R,1.0,Lorentz5Momentum());
	}
	}

	Energy QTildeReconstructor::momConsEq(const double & k,
	const Energy & root_s,
	const JetKinVect & jets) const {
	static const Energy2 eps=1e-8*GeV2;
	Energy dum = ZERO;
	for(JetKinVect::const_iterator it = jets.begin(); it != jets.end(); ++it) {
	Energy2 dum2 = (it->q).m2() + sqr(k)*(it->p).vect().mag2();
	if(dum2 < ZERO) {
	if(dum2 < -eps) throw KinematicsReconstructionVeto();
	dum2 = ZERO;
	}
	dum += sqrt(dum2);
	}
	return dum - root_s;
	}

	void QTildeReconstructor::boostChain(tPPtr p, const LorentzRotation &bv,
	tPPtr & parent) const {
	if(!p->parents().empty()) boostChain(p->parents()[0], bv,parent);
	else parent=p;
	p->transform(bv);
	if(p->children().size()==2) {
	if(dynamic_ptr_cast<ShowerParticlePtr>(p->children()[1]))
	deepTransform(p->children()[1],bv);
	}
	}

	namespace {

	bool sortJets(ShowerProgenitorPtr j1, ShowerProgenitorPtr j2) {
	return j1->highestpT()>j2->highestpT();
	}

	}

	void QTildeReconstructor::
	reconstructGeneralSystem(vector<ShowerProgenitorPtr> & ShowerHardJets) const {
	// find initial- and final-state systems
	ColourSingletSystem in,out;
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	if(ShowerHardJets[ix]->progenitor()->isFinalState())
	out.jets.push_back(ShowerHardJets[ix]);
	else
	in.jets.push_back(ShowerHardJets[ix]);
	}
	// reconstruct initial-initial system
	LorentzRotation toRest,fromRest;
	bool applyBoost(false);
	// reconstruct initial-initial system
	reconstructInitialInitialSystem(applyBoost,toRest,fromRest,in.jets);
	// reconstruct the final-state systems
	reconstructFinalStateSystem(applyBoost,toRest,fromRest,out.jets);
	}


	void QTildeReconstructor::
	reconstructFinalFirst(vector<ShowerProgenitorPtr> & ShowerHardJets) const {
	static const Energy2 minQ2 = 1e-4*GeV2;
	map<ShowerProgenitorPtr,bool> used;
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	used[ShowerHardJets[ix]] = false;
	} // first to the final-state reconstruction of any systems which need it
	set<ShowerProgenitorPtr> outgoing;
	// first find any particles with final state partners
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	if(ShowerHardJets[ix]->progenitor()->isFinalState()&&
	ShowerHardJets[ix]->progenitor()->partner()&&
	ShowerHardJets[ix]->progenitor()->partner()->isFinalState()) outgoing.insert(ShowerHardJets[ix]);
	}
	// then find the colour partners
	if(!outgoing.empty()) {
	set<ShowerProgenitorPtr> partners;
	for(set<ShowerProgenitorPtr>::const_iterator it=outgoing.begin();it!=outgoing.end();++it) {
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	if((**it).progenitor()->partner()==ShowerHardJets[ix]->progenitor()) {
	partners.insert(ShowerHardJets[ix]);
	break;
	}
	}
	}
	outgoing.insert(partners.begin(),partners.end());
	}
	// do the final-state reconstruction if needed
	if(!outgoing.empty()) {
	assert(outgoing.size()!=1);
	LorentzRotation toRest,fromRest;
	vector<ShowerProgenitorPtr> outgoingJets(outgoing.begin(),outgoing.end());
	reconstructFinalStateSystem(false,toRest,fromRest,outgoingJets);
	}
	// Now do any initial-final systems which are needed
	vector<ColourSingletSystem> IFSystems;
	// find the systems N.B. can have duplicates
	// find initial-state with FS partners or FS with IS partners
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	if(!ShowerHardJets[ix]->progenitor()->isFinalState()&&
	ShowerHardJets[ix]->progenitor()->partner()&&
	ShowerHardJets[ix]->progenitor()->partner()->isFinalState()) {
	IFSystems.push_back(ColourSingletSystem(IF,ShowerHardJets[ix]));
	}
	else if(ShowerHardJets[ix]->progenitor()->isFinalState()&&
	ShowerHardJets[ix]->progenitor()->partner()&&
	!ShowerHardJets[ix]->progenitor()->partner()->isFinalState()) {
	IFSystems.push_back(ColourSingletSystem(IF,ShowerHardJets[ix]));
	}
	}
	// then add the partners
	for(unsigned int is=0;is<IFSystems.size();++is) {
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	if(IFSystems[is].jets[0]->progenitor()->partner()==ShowerHardJets[ix]->progenitor()) {
	IFSystems[is].jets.push_back(ShowerHardJets[ix]);
	}
	}
	// ensure incoming first
	if(IFSystems[is].jets[0]->progenitor()->isFinalState())
	swap(IFSystems[is].jets[0],IFSystems[is].jets[1]);
	}
	if(!IFSystems.empty()) {
	unsigned int istart = UseRandom::irnd(IFSystems.size());
	unsigned int istop=IFSystems.size();
	for(unsigned int is=istart;is<=istop;++is) {
	if(is==IFSystems.size()) {
	if(istart!=0) {
	istop = istart-1;
	is=0;
	}
	else break;
	}
	// skip duplicates
	if(used[IFSystems[is].jets[0]] &&
	used[IFSystems[is].jets[1]] ) continue;
	if(IFSystems[is].jets[0]->original()&&IFSystems[is].jets[0]->original()->parents().empty()) continue;
	Lorentz5Momentum psum;
	for(unsigned int ix=0;ix<IFSystems[is].jets.size();++ix) {
	if(IFSystems[is].jets[ix]->progenitor()->isFinalState())
	psum += IFSystems[is].jets[ix]->progenitor()->momentum();
	else
	psum -= IFSystems[is].jets[ix]->progenitor()->momentum();
	}
	if(-psum.m2()>minQ2) {
	reconstructInitialFinalSystem(IFSystems[is].jets);
	for(unsigned int ix=0;ix<IFSystems[is].jets.size();++ix) {
	used[IFSystems[is].jets[ix]] = true;
	}
	}
	}
	}
	// now we finally need to handle the initial state system
	ColourSingletSystem in,out;
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	if(ShowerHardJets[ix]->progenitor()->isFinalState())
	out.jets.push_back(ShowerHardJets[ix]);
	else
	in.jets.push_back(ShowerHardJets[ix]);
	}
	// reconstruct initial-initial system
	bool doRecon = false;
	for(unsigned int ix=0;ix<in.jets.size();++ix) {
	if(!used[in.jets[ix]]) {
	doRecon = true;
	break;
	}
	}
	LorentzRotation toRest,fromRest;
	bool applyBoost(false);
	if(doRecon) {
	reconstructInitialInitialSystem(applyBoost,toRest,fromRest,in.jets);
	}
	// reconstruct the final-state systems
	if(!doRecon) {
	for(unsigned int ix=0;ix<out.jets.size();++ix) {
	if(!used[out.jets[ix]]) {
	doRecon = true;
	break;
	}
	}
	}
	if(doRecon) {
	reconstructFinalStateSystem(applyBoost,toRest,fromRest,out.jets);
	}
	}

	void QTildeReconstructor::
	reconstructColourPartner(vector<ShowerProgenitorPtr> & ShowerHardJets) const {
	static const Energy2 minQ2 = 1e-4*GeV2;
	// sort the vector by hardness of emission
	std::sort(ShowerHardJets.begin(),ShowerHardJets.end(),sortJets);
	// map between particles and progenitors for easy lookup
	map<ShowerParticlePtr,ShowerProgenitorPtr> progenitorMap;
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	progenitorMap[ShowerHardJets[ix]->progenitor()] = ShowerHardJets[ix];
	}
	// check that the IF systems can be reconstructed
	bool canReconstruct = true;
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	tShowerParticlePtr progenitor = ShowerHardJets[ix]->progenitor();
	tShowerParticlePtr partner = progenitor->partner();
	if(!partner) continue;
	else if((progenitor->isFinalState() &&
	!partner->isFinalState()) \|\|
	(!progenitor->isFinalState() &&
	partner->isFinalState()) ) {
	vector<ShowerProgenitorPtr> jets(2);
	jets[0] = ShowerHardJets[ix];
	jets[1] = progenitorMap[partner];
	Lorentz5Momentum psum;
	for(unsigned int iy=0;iy<jets.size();++iy) {
	if(jets[iy]->progenitor()->isFinalState())
	psum += jets[iy]->progenitor()->momentum();
	else
	psum -= jets[iy]->progenitor()->momentum();
	}
	if(-psum.m2()<minQ2) {
	canReconstruct = false;
	break;
	}
	}
	}
	if(!canReconstruct) {
	reconstructGeneralSystem(ShowerHardJets);
	return;
	}
	map<ShowerProgenitorPtr,bool> used;
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	used[ShowerHardJets[ix]] = false;
	}
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	// skip jets which have already been handled
	if(ShowerHardJets[ix]->reconstructed()==ShowerProgenitor::done) continue;
	// already reconstructed
	if(used[ShowerHardJets[ix]]) continue;
	// no partner continue
	tShowerParticlePtr progenitor = ShowerHardJets[ix]->progenitor();
	tShowerParticlePtr partner = progenitor->partner();
	if(!partner) {
	// check if there's a daughter tree which also needs boosting
	Lorentz5Momentum porig = progenitor->momentum();
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
	for(tit = _currentTree->treelinks().begin();
	tit != _currentTree->treelinks().end();++tit) {
	// if there is, boost it
	if(tit->second.first && tit->second.second==progenitor) {
	Lorentz5Momentum pnew = tit->first->incomingLines().begin()
	->first->progenitor()->momentum();
	pnew *= tit->first->transform();
	Lorentz5Momentum pdiff = porig-pnew;
	Energy2 test = sqr(pdiff.x()) + sqr(pdiff.y()) +
	sqr(pdiff.z()) + sqr(pdiff.t());
	LorentzRotation rot;
	if(test>1e-6*GeV2) rot = solveBoost(porig,pnew);
	tit->first->transform(rot,false);
	_treeBoosts[tit->first].push_back(rot);
	}
	}
	ShowerHardJets[ix]->reconstructed(ShowerProgenitor::done);
	continue;
	}
	// do the reconstruction
	// final-final
	if(progenitor->isFinalState() &&
	partner->isFinalState() ) {
	LorentzRotation toRest,fromRest;
	vector<ShowerProgenitorPtr> jets(2);
	jets[0] = ShowerHardJets[ix];
	jets[1] = progenitorMap[partner];
	if(_reconopt==4 && jets[1]->reconstructed()==ShowerProgenitor::notReconstructed)
	jets[1]->reconstructed(ShowerProgenitor::dontReconstruct);
	reconstructFinalStateSystem(false,toRest,fromRest,jets);
	if(_reconopt==4 && jets[1]->reconstructed()==ShowerProgenitor::dontReconstruct)
	jets[1]->reconstructed(ShowerProgenitor::notReconstructed);
	used[jets[0]] = true;
	if(_reconopt==3) used[jets[1]] = true;
	}
	// initial-final
	else if((progenitor->isFinalState() &&
	!partner->isFinalState()) \|\|
	(!progenitor->isFinalState() &&
	partner->isFinalState()) ) {
	vector<ShowerProgenitorPtr> jets(2);
	jets[0] = ShowerHardJets[ix];
	jets[1] = progenitorMap[partner];
	if(jets[0]->progenitor()->isFinalState()) swap(jets[0],jets[1]);
	if(jets[0]->original()&&jets[0]->original()->parents().empty()) continue;
	Lorentz5Momentum psum;
	for(unsigned int iy=0;iy<jets.size();++iy) {
	if(jets[iy]->progenitor()->isFinalState())
	psum += jets[iy]->progenitor()->momentum();
	else
	psum -= jets[iy]->progenitor()->momentum();
	}
	if(_reconopt==4 && progenitorMap[partner]->reconstructed()==ShowerProgenitor::notReconstructed)
	progenitorMap[partner]->reconstructed(ShowerProgenitor::dontReconstruct);
	reconstructInitialFinalSystem(jets);
	if(_reconopt==4 && progenitorMap[partner]->reconstructed()==ShowerProgenitor::dontReconstruct)
	progenitorMap[partner]->reconstructed(ShowerProgenitor::notReconstructed);
	used[ShowerHardJets[ix]] = true;
	if(_reconopt==3) used[progenitorMap[partner]] = true;
	}
	// initial-initial
	else if(!progenitor->isFinalState() &&
	!partner->isFinalState() ) {
	ColourSingletSystem in,out;
	in.jets.push_back(ShowerHardJets[ix]);
	in.jets.push_back(progenitorMap[partner]);
	for(unsigned int iy=0;iy<ShowerHardJets.size();++iy) {
	if(ShowerHardJets[iy]->progenitor()->isFinalState())
	out.jets.push_back(ShowerHardJets[iy]);
	}
	LorentzRotation toRest,fromRest;
	bool applyBoost(false);
	if(_reconopt==4 && in.jets[1]->reconstructed()==ShowerProgenitor::notReconstructed)
	in.jets[1]->reconstructed(ShowerProgenitor::dontReconstruct);
	reconstructInitialInitialSystem(applyBoost,toRest,fromRest,in.jets);
	if(_reconopt==4 && in.jets[1]->reconstructed()==ShowerProgenitor::dontReconstruct)
	in.jets[1]->reconstructed(ShowerProgenitor::notReconstructed);
	used[in.jets[0]] = true;
	if(_reconopt==3) used[in.jets[1]] = true;
	for(unsigned int iy=0;iy<out.jets.size();++iy) {
	if(out.jets[iy]->reconstructed()==ShowerProgenitor::notReconstructed)
	out.jets[iy]->reconstructed(ShowerProgenitor::dontReconstruct);
	}
	// reconstruct the final-state systems
	LorentzRotation finalBoosts;
	finalBoosts.transform( toRest);
	finalBoosts.transform(fromRest);
	for(unsigned int iy=0;iy<out.jets.size();++iy) {
	deepTransform(out.jets[iy]->progenitor(),finalBoosts);
	}
	for(unsigned int iy=0;iy<out.jets.size();++iy) {
	if(out.jets[iy]->reconstructed()==ShowerProgenitor::dontReconstruct)
	out.jets[iy]->reconstructed(ShowerProgenitor::notReconstructed);
	}
	}
	}
	}

	bool QTildeReconstructor::
	inverseDecayRescalingFactor(vector<Lorentz5Momentum> pout,
	vector<Energy> mon,Energy roots,
	Lorentz5Momentum ppartner, Energy mbar,
	double & k1, double & k2) const {
	ThreeVector<Energy> qtotal;
	vector<Energy2> pmag;
	for(unsigned int ix=0;ix<pout.size();++ix) {
	pmag.push_back(pout[ix].vect().mag2());
	qtotal+=pout[ix].vect();
	}
	Energy2 dot1 = qtotal*ppartner.vect();
	Energy2 qmag2=qtotal.mag2();
	double a = -dot1/qmag2;
	static const Energy eps=1e-10*GeV;
	unsigned int itry(0);
	Energy numer(ZERO),denom(ZERO);
	k1=1.;
	do {
	++itry;
	numer=denom=0.*GeV;
	double k12=sqr(k1);
	for(unsigned int ix=0;ix<pout.size();++ix) {
	Energy en = sqrt(pmag[ix]/k12+sqr(mon[ix]));
	numer += en;
	denom += pmag[ix]/en;
	}
	Energy en = sqrt(qmag2/k12+sqr(mbar));
	numer += en-roots;
	denom += qmag2/en;
	k1 += numer/denomk12k1;
	if(abs(k1)>1e10) return false;
	}
	while (abs(numer)>eps&&itry<100);
	k1 = abs(k1);
	k2 = a*k1;
	return itry<100;
	}

	void QTildeReconstructor::
	deconstructInitialFinalSystem(HardTreePtr tree,vector<HardBranchingPtr> jets,
	cEvolverPtr evolver,
	ShowerInteraction::Type type) const {
	HardBranchingPtr incoming;
	Lorentz5Momentum pin[2],pout[2],pbeam;
	HardBranchingPtr initial;
	Energy mc(ZERO);
	for(unsigned int ix=0;ix<jets.size();++ix) {
	// final-state parton
	if(jets[ix]->status()==HardBranching::Outgoing) {
	pout[0] += jets[ix]->branchingParticle()->momentum();
	mc = jets[ix]->branchingParticle()->thePEGBase() ?
	jets[ix]->branchingParticle()->thePEGBase()->mass() :
	jets[ix]->branchingParticle()->dataPtr()->mass();
	}
	// initial-state parton
	else {
	pin[0] += jets[ix]->branchingParticle()->momentum();
	initial = jets[ix];
	pbeam = jets[ix]->beam()->momentum();
	Energy scale=pbeam.t();
	pbeam = Lorentz5Momentum(ZERO,pbeam.vect().unit()*scale);
	incoming = jets[ix];
	while(incoming->parent()) incoming = incoming->parent();
	}
	}
	if(jets.size()>2) {
	pout[0].rescaleMass();
	mc = pout[0].mass();
	}
	// work out the boost to the Breit frame
	Lorentz5Momentum pa = pout[0]-pin[0];
	Axis axis(pa.vect().unit());
	LorentzRotation rot;
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	if(axis.perp2()>0.) {
	rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	rot.rotateX(Constants::pi);
	rot.boostZ( pa.e()/pa.vect().mag());
	}
	// transverse part
	Lorentz5Momentum paxis=rot*pbeam;
	Boost trans = -1./paxis.e()*paxis.vect();
	trans.setZ(0.);
	rot.boost(trans);
	pa *= rot;
	// reference vectors
	Lorentz5Momentum n1(ZERO,ZERO,-pa.z(),-pa.z());
	Lorentz5Momentum n2(ZERO,ZERO, pa.z(),-pa.z());
	Energy2 n1n2 = n1*n2;
	// decompose the momenta
	Lorentz5Momentum qbp=rotpin[0],qcp= rotpout[0];
	double a[2],b[2];
	a[0] = n2*qbp/n1n2;
	b[0] = n1*qbp/n1n2;
	a[1] = n2*qcp/n1n2;
	b[1] = n1*qcp/n1n2;
	Lorentz5Momentum qperp = qbp-a[0]n1-b[0]n2;
	// before reshuffling
	Energy Q = abs(pa.z());
	double c = sqr(mc/Q);
	Lorentz5Momentum pb(ZERO,ZERO,0.5Q(1.+c),0.5Q(1.+c));
	Lorentz5Momentum pc(ZERO,ZERO,0.5Q(c-1.),0.5Q(1.+c));
	double anew[2],bnew[2];
	anew[0] = pb*n2/n1n2;
	bnew[0] = 0.5*(qbp.m2()-qperp.m2())/n1n2/anew[0];
	bnew[1] = pc*n1/n1n2;
	anew[1] = 0.5*qcp.m2()/bnew[1]/n1n2;
	Lorentz5Momentum qnewb = (anew[0]n1+bnew[0]n2+qperp);
	Lorentz5Momentum qnewc = (anew[1]n1+bnew[1]n2);
	// initial-state boost
	LorentzRotation rotinv=rot.inverse();
	LorentzRotation transb=rotinvsolveBoostZ(qnewb,qbp)rot;
	// final-state boost
	LorentzRotation transc=rotinvsolveBoost(qnewc,qcp)rot;
	// this will need changing for more than one outgoing particle
	// set the pvectors
	for(unsigned int ix=0;ix<jets.size();++ix) {
	if(jets[ix]->status()==HardBranching::Incoming) {
	jets[ix]->pVector(pbeam);
	jets[ix]->showerMomentum(rotinv*pb);
	incoming->pVector(jets[ix]->pVector());
	}
	else {
	jets[ix]->pVector(rotinv*pc);
	jets[ix]->showerMomentum(jets[ix]->pVector());
	}
	}
	// find the colour partners
	ShowerParticleVector particles;
	vector<Lorentz5Momentum> ptemp;
	set<HardBranchingPtr>::const_iterator cjt;
	for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) {
	ptemp.push_back((**cjt).branchingParticle()->momentum());
	(cjt).branchingParticle()->set5Momentum((cjt).showerMomentum());
	particles.push_back((**cjt).branchingParticle());
	}
	evolver->showerModel()->partnerFinder()
	->setInitialEvolutionScales(particles,false,type,false);
	unsigned int iloc(0);
	for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) {
	// reset the momentum
	(**cjt).branchingParticle()->set5Momentum(ptemp[iloc]);
	++iloc;
	}
	for(vector<HardBranchingPtr>::const_iterator cjt=jets.begin();
	cjt!=jets.end();++cjt) {
	// sort out the partners
	tShowerParticlePtr partner =
	(*cjt)->branchingParticle()->partner();
	if(!partner) continue;
	tHardBranchingPtr branch;
	for(set<HardBranchingPtr>::const_iterator
	clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) {
	if((**clt).branchingParticle()==partner) {
	(*cjt).colourPartner(clt);
	branch=*clt;
	break;
	}
	}
	// compute the reference vectors
	// both incoming, should all ready be done
	if((**cjt).status()==HardBranching::Incoming &&
	branch->status()==HardBranching::Incoming) {
	Energy etemp = (*cjt)->beam()->momentum().z();
	Lorentz5Momentum nvect(ZERO, ZERO,-etemp, abs(etemp));
	tHardBranchingPtr branch2 = *cjt;
	(**cjt).nVector(nvect);
	while (branch2->parent()) {
	branch2=branch2->parent();
	branch2->nVector(nvect);
	}
	}
	// both outgoing
	else if((**cjt).status()==HardBranching::Outgoing&&
	branch->status()==HardBranching::Outgoing) {
	Boost boost=((*cjt)->pVector()+branch->pVector()).findBoostToCM();
	Lorentz5Momentum pcm = branch->pVector();
	pcm.boost(boost);
	Lorentz5Momentum nvect = Lorentz5Momentum(ZERO,pcm.vect());
	nvect.boost( -boost);
	(**cjt).nVector(nvect);
	}
	else if((**cjt).status()==HardBranching::Incoming) {
	Lorentz5Momentum pa = -(**cjt).showerMomentum()+branch->showerMomentum();
	Lorentz5Momentum pb = (**cjt).showerMomentum();
	Axis axis(pa.vect().unit());
	LorentzRotation rot;
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	if(axis.perp2()>1e-20) {
	rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	rot.rotateX(Constants::pi);
	}
	if(abs(1.-pa.e()/pa.vect().mag())>1e-6) rot.boostZ( pa.e()/pa.vect().mag());
	pb*=rot;
	Boost trans = -1./pb.e()*pb.vect();
	trans.setZ(0.);
	rot.boost(trans);
	Energy scale=(**cjt).beam()->momentum().t();
	Lorentz5Momentum pbasis(ZERO,(*cjt).beam()->momentum().vect().unit()scale);
	Lorentz5Momentum pcm = rot*pbasis;
	rot.invert();
	Lorentz5Momentum nvect = rot*Lorentz5Momentum(ZERO,-pcm.vect());
	(**cjt).nVector(nvect);
	tHardBranchingPtr branch2 = *cjt;
	while (branch2->parent()) {
	branch2=branch2->parent();
	branch2->nVector(nvect);
	}
	}
	else if(branch->status()==HardBranching::Incoming) {
	Lorentz5Momentum nvect=Lorentz5Momentum(ZERO,branch->showerMomentum().vect());
	(**cjt).nVector(nvect);
	}
	}
	// now compute the new momenta
	for(vector<HardBranchingPtr>::const_iterator cjt=jets.begin();
	cjt!=jets.end();++cjt) {
	if((**cjt).status()==HardBranching::Outgoing) {
	(**cjt).setMomenta(transc,1.,Lorentz5Momentum());
	}
	}
	incoming->setMomenta(transb,1.,Lorentz5Momentum());
	}

	void QTildeReconstructor::deepTransform(PPtr particle,
	const LorentzRotation & r,
	bool match,
	PPtr original) const {
	if(_boosts.find(particle)!=_boosts.end()) {
	_boosts[particle].push_back(r);
	}
	Lorentz5Momentum porig = particle->momentum();
	if(!original) original = particle;
	for ( int i = 0, N = particle->children().size(); i < N; ++i ) {
	deepTransform(particle->children()[i],r,
	particle->children()[i]->id()==original->id()&&match,original);
	}
	particle->transform(r);
	// transform the p and n vectors
	ShowerParticlePtr sparticle = dynamic_ptr_cast<ShowerParticlePtr>(particle);
	if(sparticle && sparticle->showerKinematics()) {
	sparticle->showerKinematics()->transform(r);
	}
	if ( particle->next() ) deepTransform(particle->next(),r,match,original);
	if(!match) return;
	if(!particle->children().empty()) return;
	// force the mass shell
	if(particle->dataPtr()->stable()) {
	Lorentz5Momentum ptemp = particle->momentum();
	ptemp.rescaleEnergy();
	particle->set5Momentum(ptemp);
	}
	// check if there's a daughter tree which also needs boosting
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
	for(tit = _currentTree->treelinks().begin();
	tit != _currentTree->treelinks().end();++tit) {
	// if there is, boost it
	if(tit->second.first && tit->second.second==original) {
	Lorentz5Momentum pnew = tit->first->incomingLines().begin()
	->first->progenitor()->momentum();
	pnew *= tit->first->transform();
	Lorentz5Momentum pdiff = porig-pnew;
	Energy2 test = sqr(pdiff.x()) + sqr(pdiff.y()) +
	sqr(pdiff.z()) + sqr(pdiff.t());
	LorentzRotation rot;
	if(test>1e-6*GeV2) rot = solveBoost(porig,pnew);
	tit->first->transform(r*rot,false);
	_treeBoosts[tit->first].push_back(r*rot);
	}
	}
	}

	void QTildeReconstructor::reconstructFinalFinalOffShell(JetKinVect orderedJets,
	Energy2 s,
	bool recursive) const {
	JetKinVect::iterator jit;
	jit = orderedJets.begin(); ++jit;
	// 4-momentum of recoiling system
	Lorentz5Momentum psum;
	for( ; jit!=orderedJets.end(); ++jit) psum += jit->p;
	psum.rescaleMass();
	// calculate the 3-momentum rescaling factor
	Energy2 m1sq(orderedJets.begin()->q.m2()),m2sq(psum.m2());
	Energy4 num = sqr(s - m1sq - m2sq) - 4.m1sqm2sq;
	if(num<ZERO) throw KinematicsReconstructionVeto();
	double k = sqrt( num / (4.sorderedJets.begin()->p.vect().mag2()) );
	// boost the most off-shell
	LorentzRotation B1 = solveBoost(k, orderedJets.begin()->q, orderedJets.begin()->p);
	deepTransform(orderedJets.begin()->parent,B1);
	// boost everything else
	// first to rescale
	LorentzRotation B2 = solveBoost(k, psum, psum);
	// and then to rest frame of new system
	Lorentz5Momentum pnew = B2*psum;
	pnew.rescaleMass();
	B2.transform(pnew.findBoostToCM());
	// apply transform (calling routine ensures at least 3 elements)
	jit = orderedJets.begin(); ++jit;
	for(;jit!=orderedJets.end();++jit) {
	deepTransform(jit->parent,B2);
	jit->p *= B2;
	jit->q *= B2;
	}
	JetKinVect newJets(orderedJets.begin()+1,orderedJets.end());
	// final reconstruction
	if(newJets.size()==2 \|\| !recursive ) {
	// rescaling factor
	double k = solveKfactor(psum.m(), newJets);
	// rescale jets in the new CMF
	for(JetKinVect::iterator it = newJets.begin(); it != newJets.end(); ++it) {
	LorentzRotation Trafo = solveBoost(k, it->q, it->p);
	deepTransform(it->parent,Trafo);
	}
	}
	// recursive
	else {
	std::sort(newJets.begin(),newJets.end(),JetOrdering());
	reconstructFinalFinalOffShell(newJets,psum.m2(),recursive);
	}
	// finally boost back from new CMF
	LorentzRotation back(-pnew.findBoostToCM());
	for(JetKinVect::iterator it = newJets.begin(); it != newJets.end(); ++it) {
	deepTransform(it->parent,back);
	}
	}

	Energy QTildeReconstructor::findMass(HardBranchingPtr branch) const {
	// KH - 230909 - If the particle has no children then it will
	// not have showered and so it should be "on-shell" so we can
	// get it's mass from it's momentum. This means that the
	// inverseRescalingFactor doesn't give any nans or do things
	// it shouldn't if it gets e.g. two Z bosons generated with
	// off-shell masses. This is for sure not the best solution.
	// PR 1/1/10 modification to previous soln
	// PR 28/8/14 change to procedure and factorize into a function
	if(branch->children().empty()) {
	return branch->branchingParticle()->mass();
	}
	else if(!branch->children().empty() &&
	!branch->branchingParticle()->dataPtr()->stable() ) {
	for(unsigned int ix=0;ix<branch->children().size();++ix) {
	if(branch->branchingParticle()->id()==
	branch->children()[ix]->branchingParticle()->id())
	return findMass(branch->children()[ix]);
	}
	}
	return branch->branchingParticle()->dataPtr()->mass();
	}

	vector<double>
	QTildeReconstructor::inverseInitialStateRescaling(double & x1, double & x2,
	const Lorentz5Momentum & pold,
	const vector<Lorentz5Momentum> & p,
	const vector<Lorentz5Momentum> & pq) const {
	// hadronic CMS
	Energy2 s = (pq[0] +pq[1] ).m2();
	// partonic CMS
	Energy MDY = pold.m();
	// find alpha, beta and pt
	Energy2 p12=pq[0]*pq[1];
	double a[2],b[2];
	Lorentz5Momentum pt[2];
	for(unsigned int ix=0;ix<2;++ix) {
	a[ix] = p[ix]*pq[1]/p12;
	b [ix] = p[ix]*pq[0]/p12;
	pt[ix] = p[ix]-a[ix]pq[0]-b[ix]pq[1];
	}
	// compute kappa
	// we always want to preserve the mass of the system
	double k1(1.),k2(1.);
	if(_initialStateReconOption==0) {
	double rap=pold.rapidity();
	x1 = MDY/sqrt(sexp(2.rap));
	x2 = sqr(MDY)/s/x1;
	if(pq[0].z()>ZERO) swap(x1,x2);
	k1=a[0]/x1;
	k2=b[1]/x2;
	}
	// longitudinal momentum
	else if(_initialStateReconOption==1) {
	double A = 1.;
	double C = -sqr(MDY)/s;
	double B = pq[0].z()>ZERO ? 2.pold.z()/sqrt(s) : -2.pold.z()/sqrt(s);
	if(abs(B)>1e-10) {
	double discrim = 1.-4.AC/sqr(B);
	if(discrim < 0.) throw KinematicsReconstructionVeto();
	x1 = B>0. ? 0.5B/A(1.+sqrt(discrim)) : 0.5B/A(1.-sqrt(discrim));
	}
	else {
	x1 = -C/A;
	if( x1 <= 0.) throw KinematicsReconstructionVeto();
	x1 = sqrt(x1);
	}
	x2 = sqr(MDY)/s/x1;
	k1=a[0]/x1;
	k2=b[1]/x2;
	}
	// preserve mass and don't scale the softer system
	// to reproduce the dipole kinematics
	else if(_initialStateReconOption==2) {
	// in this case kp = k1 or k2 depending on who's the harder guy
	k1 = a[0]b[1]s/sqr(MDY);
	if ( pt[0].perp2() < pt[1].perp2() ) swap(k1,k2);
	x1 = a[0]/k1;
	x2 = b[1]/k2;
	}
	else
	assert(false);
	// decompose the momenta
	double anew[2] = {a[0]/k1,a[1]*k2};
	double bnew[2] = {b[0]*k1,b[1]/k2};
	vector<double> boost(2);
	for(unsigned int ix=0;ix<2;++ix) {
	boost[ix] = getBeta(a [ix]+b [ix], a[ix] -b [ix],
	anew[ix]+bnew[ix], anew[ix]-bnew[ix]);
	if (pq[0].z() > ZERO) boost[ix] *= -1.;
	}
	return boost;
	}

	vector<double>
	QTildeReconstructor::initialStateRescaling(double x1, double x2,
	const Lorentz5Momentum & pold,
	const vector<Lorentz5Momentum> & p,
	const vector<Lorentz5Momentum> & pq,
	const vector<Energy>& highestpts) const {
	Energy2 S = (pq[0]+pq[1]).m2();
	// find alphas and betas in terms of desired basis
	Energy2 p12 = pq[0]*pq[1];
	double a[2] = {p[0]pq[1]/p12,p[1]pq[1]/p12};
	double b[2] = {p[0]pq[0]/p12,p[1]pq[0]/p12};
	Lorentz5Momentum p1p = p[0] - a[0]pq[0] - b[0]pq[1];
	Lorentz5Momentum p2p = p[1] - a[1]pq[0] - b[1]pq[1];
	// compute kappa
	// we always want to preserve the mass of the system
	Energy MDY = pold.m();
	Energy2 A = a[0]b[1]S;
	Energy2 B = Energy2(sqr(MDY)) - (a[0]b[0]+a[1]b[1])*S - (p1p+p2p).m2();
	Energy2 C = a[1]b[0]S;
	double rad = 1.-4.AC/sqr(B);
	if(rad < 0.) throw KinematicsReconstructionVeto();
	double kp = B/(2.A)(1.+sqrt(rad));
	// now compute k1
	// conserve rapidity
	double k1(0.);
	double k2(0.);
	if(_initialStateReconOption==0) {
	rad = kp(b[0]+kpb[1])/(kp*a[0]+a[1]);
	rad = pq[0].z()<ZERO ? exp(-2.pold.rapidity()) : exp(2.*pold.rapidity());
	if(rad <= 0.) throw KinematicsReconstructionVeto();
	k1 = sqrt(rad);
	k2 = kp/k1;
	}
	// conserve longitudinal momentum
	else if(_initialStateReconOption==1) {
	double a2 = (a[0]+a[1]/kp);
	double b2 = -x2+x1;
	double c2 = -(b[1]*kp+b[0]);
	if(abs(b2)>1e-10) {
	double discrim = 1.-4.a2c2/sqr(b2);
	if(discrim < 0.) throw KinematicsReconstructionVeto();
	k1 = b2>0. ? 0.5b2/a2(1.+sqrt(discrim)) : 0.5b2/a2(1.-sqrt(discrim));
	}
	else {
	k1 = -c2/a2;
	if( k1 <= 0.) throw KinematicsReconstructionVeto();
	k1 = sqrt(k1);
	}
	k2 = kp/k1;
	}
	// preserve mass and don't scale the softer system
	// to reproduce the dipole kinematics
	else if(_initialStateReconOption==2) {
	// in this case kp = k1 or k2 depending on who's the harder guy
	k1 = kp; k2 = 1.;
	if ( highestpts[0] < highestpts[1] )
	swap(k1,k2);
	}
	else
	assert(false);
	// calculate the boosts
	vector<double> beta(2);
	beta[0] = getBeta((a[0]+b[0]), (a[0]-b[0]), (k1a[0]+b[0]/k1), (k1a[0]-b[0]/k1));
	beta[1] = getBeta((a[1]+b[1]), (a[1]-b[1]), (a[1]/k2+k2b[1]), (a[1]/k2-k2b[1]));
	if (pq[0].z() > ZERO) {
	beta[0] = -beta[0];
	beta[1] = -beta[1];
	}
	return beta;
	}

	void QTildeReconstructor::
	reconstructColourSinglets(vector<ShowerProgenitorPtr> & ShowerHardJets,
	ShowerInteraction::Type type) const {
	// identify and catagorize the colour singlet systems
	unsigned int nnun(0),nnii(0),nnif(0),nnf(0),nni(0);
	vector<ColourSingletSystem>
	systems(identifySystems(set<ShowerProgenitorPtr>(ShowerHardJets.begin(),ShowerHardJets.end()),
	nnun,nnii,nnif,nnf,nni));
	// now decide what to do
	// initial-initial connection and final-state colour singlet systems
	LorentzRotation toRest,fromRest;
	bool applyBoost(false),general(false);
	// Drell-Yan type
	if(nnun==0&&nnii==1&&nnif==0&&nnf>0&&nni==0) {
	// reconstruct initial-initial system
	for(unsigned int ix=0;ix<systems.size();++ix) {
	if(systems[ix].type==II)
	reconstructInitialInitialSystem(applyBoost,toRest,fromRest,
	systems[ix].jets);
	}
	if(type==ShowerInteraction::QED\|\|
	type==ShowerInteraction::Both) {
	combineFinalState(systems);
	general=false;
	}
	}
	// DIS and VBF type
	else if(nnun==0&&nnii==0&&((nnif==1&&nnf>0&&nni==1)\|\|
	(nnif==2&& nni==0))) {
	// check these systems can be reconstructed
	for(unsigned int ix=0;ix<systems.size();++ix) {
	// compute q^2
	if(systems[ix].type!=IF) continue;
	Lorentz5Momentum q;
	for(unsigned int iy=0;iy<systems[ix].jets.size();++iy) {
	if(systems[ix].jets[iy]->progenitor()->isFinalState())
	q += systems[ix].jets[iy]->progenitor()->momentum();
	else
	q -= systems[ix].jets[iy]->progenitor()->momentum();
	}
	q.rescaleMass();
	// check above cut
	if(abs(q.m())>=_minQ) continue;
	if(nnif==1&&nni==1) {
	throw KinematicsReconstructionVeto();
	}
	else {
	general = true;
	break;
	}
	}
	if(!general) {
	for(unsigned int ix=0;ix<systems.size();++ix) {
	if(systems[ix].type==IF)
	reconstructInitialFinalSystem(systems[ix].jets);
	}
	}
	}
	// e+e- type
	else if(nnun==0&&nnii==0&&nnif==0&&nnf>0&&nni==2) {
	general = type==ShowerInteraction::QED \|\| type == ShowerInteraction::Both;
	}
	// general type
	else {
	general = true;
	}
	// final-state systems except for general recon
	if(!general) {
	for(unsigned int ix=0;ix<systems.size();++ix) {
	if(systems[ix].type==F)
	reconstructFinalStateSystem(applyBoost,toRest,fromRest,
	systems[ix].jets);
	}
	}
	else {
	reconstructGeneralSystem(ShowerHardJets);
	}
	}

	void QTildeReconstructor::findInitialBoost(const Lorentz5Momentum & pold,
	const Lorentz5Momentum & pnew,
	LorentzRotation & toRest,
	LorentzRotation & fromRest) const {
	// do one boost
	if(_initialBoost==0) {
	toRest = LorentzRotation(pold.findBoostToCM());
	fromRest = LorentzRotation(pnew.boostVector());
	}
	else if(_initialBoost==1) {
	// boost to rest frame
	// first transverse
	toRest = Boost(-pold.x()/pold.t(),-pold.y()/pold.t(),0.);
	// then longitudinal
	double beta = pold.z()/sqrt(pold.m2()+sqr(pold.z()));
	toRest.boost((Boost(0.,0.,-beta)));
	// boost from rest frame
	// first apply longitudinal boost
	beta = pnew.z()/sqrt(pnew.m2()+sqr(pnew.z()));
	fromRest=LorentzRotation(Boost(0.,0.,beta));
	// then transverse one
	fromRest.boost(Boost(pnew.x()/pnew.t(),
	pnew.y()/pnew.t(),0.));
	}
	else
	assert(false);
	}
	diff --git a/Shower/Default/QTildeReconstructor.h b/Shower/Default/QTildeReconstructor.h
	--- a/Shower/Default/QTildeReconstructor.h
	+++ b/Shower/Default/QTildeReconstructor.h
	@@ -1,645 +1,646 @@
	// -- C++ --
	//
	// QTildeReconstructor.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_QTildeReconstructor_H
	#define HERWIG_QTildeReconstructor_H
	//
	// This is the declaration of the QTildeReconstructor class.
	//

	#include "Herwig/Shower/Base/KinematicsReconstructor.h"

	namespace Herwig {

	using namespace ThePEG;

	/** \ingroup Shower
	* A simple struct to store the information we need on the
	* showering
	*/
	struct JetKinStruct {

	/**
	* Parent particle of the jet
	*/
	tShowerParticlePtr parent;

	/**
	* Momentum of the particle before reconstruction
	*/
	Lorentz5Momentum p;

	/**
	* Momentum of the particle after reconstruction
	*/
	Lorentz5Momentum q;
	};

	/**
	* typedef for a vector of JetKinStruct
	*/
	typedef vector<JetKinStruct> JetKinVect;

	/**
	* Enum to identify types of colour singlet systems
	*/
	enum SystemType { UNDEFINED=-1, II, IF, F ,I };

	/**
	* Struct to store colour singlets
	*/
	template<typename Value> struct ColourSinglet {

	typedef vector<ColourSinglet<Value> > VecType;

	ColourSinglet() : type(UNDEFINED) {}

	ColourSinglet(SystemType intype,Value inpart)
	: type(intype),jets(1,inpart) {}

	/**
	* The type of system
	*/
	SystemType type;

	/**
	* The jets in the system
	*/
	vector<Value> jets;

	};

	/**
	* Struct to store a colour singlet system of particles
	*/
	typedef ColourSinglet<ShowerProgenitorPtr> ColourSingletSystem;

	/**
	* Struct to store a colour singlet shower
	*/
	typedef ColourSinglet<HardBranchingPtr> ColourSingletShower;

	/** \ingroup Shower
	*
	* This class is responsible for the kinematical reconstruction
	* after each showering step, and also for the necessary Lorentz boosts
	* in order to preserve energy-momentum conservation in the overall collision,
	* and also the invariant mass and the rapidity of the hard subprocess system.
	* In the case of multi-step showering, there will be not unnecessary
	* kinematical reconstructions.
	*
	* There is also the option of taking a set of momenta for the particles
	* and inverting the reconstruction to give the evolution variables for the
	* shower.
	*
	* Notice:
	* - although we often use the term "jet" in either methods or variables names,
	* or in comments, which could appear applicable only for QCD showering,
	* there is indeed no "dynamics" represented in this class: only kinematics
	* is involved, as the name of this class remainds. Therefore it can be used
	* for any kind of showers (QCD-,QED-,EWK-,... bremsstrahlung).
	*
	* @see ShowerParticle
	* @see ShowerKinematics
	* @see \ref QTildeReconstructorInterfaces "The interfaces"
	* defined for QTildeReconstructor.
	*/
	class QTildeReconstructor: public KinematicsReconstructor {

	public:

	/**
	* Default constructor
	*/
	QTildeReconstructor() : _reconopt(0), _initialBoost(0),
	_finalStateReconOption(0),
	_initialStateReconOption(0), _minQ(MeV) {};

	/**
	* Methods to reconstruct the kinematics of a scattering or decay process
	*/
	//@{
	/**
	* Given in input a vector of the particles which initiated the showers
	* the method does the reconstruction of such jets,
	* including the appropriate boosts (kinematics reshufflings)
	* needed to conserve the total energy-momentum of the collision
	* and preserving the invariant mass and the rapidity of the
	* hard subprocess system.
	*/
	virtual bool reconstructHardJets(ShowerTreePtr hard,
	const map<tShowerProgenitorPtr,
	pair<Energy,double> > & pt,
	- ShowerInteraction::Type type) const;
	+ ShowerInteraction::Type type,
	+ bool switchRecon) const;

	/**
	* Given in input a vector of the particles which initiated the showers
	* the method does the reconstruction of such jets,
	* including the appropriate boosts (kinematics reshufflings)
	* needed to conserve the total energy-momentum of the collision
	* and preserving the invariant mass and the rapidity of the
	* hard subprocess system.
	*/
	virtual bool reconstructDecayJets(ShowerTreePtr decay,
	ShowerInteraction::Type type) const;
	//@}

	/**
	* Methods to invert the reconstruction of the shower for
	* a scattering or decay process and calculate
	* the variables used to generate the
	* shower given the particles produced.
	* This is needed for the CKKW and POWHEG approaches
	*/
	//@{
	/**
	* Given the particles, with a history which we wish to interpret
	* as a shower reconstruct the variables used to generate the
	* shower
	*/
	virtual bool deconstructDecayJets(HardTreePtr, cEvolverPtr,
	ShowerInteraction::Type) const;

	/**
	* Given the particles, with a history which we wish to interpret
	* as a shower reconstruct the variables used to generate the shower
	* for a hard process
	*/
	virtual bool deconstructHardJets(HardTreePtr, cEvolverPtr,
	ShowerInteraction::Type) const;
	//@}

	public:

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

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

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

	protected:

	/**
	* Methods to reconstruct the kinematics of individual jets
	*/
	//@{
	/**
	* Given the particle (ShowerParticle object) that
	* originates a forward (time-like) jet, this method reconstructs the kinematics
	* of the jet. That is, by starting from the final grand-children (which
	* originates directly or indirectly from particleJetParent,
	* and which don't have children), and moving "backwards" (in a physical
	* time picture), towards the particleJetParent, the
	* ShowerKinematics objects associated with the various particles,
	* which have been created during the showering, are now completed.
	* In particular, at the end, we get the mass of the jet, which is the
	* main information we want.
	* This methods returns false if there was no radiation or rescaling required
	*/
	virtual bool reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent) const;

	/**
	* Exactly similar to the previous one, but for a space-like jet.
	* Also in this case we start from the final grand-children (which
	* are childless) of the particle which originates the jet, but in
	* this case we proceed "forward" (in the physical time picture)
	* towards the particleJetParent.
	* This methods returns false if there was no radiation or rescaling required
	*/
	bool reconstructSpaceLikeJet(const tShowerParticlePtr particleJetParent) const;

	/**
	* Exactly similar to the previous one, but for a decay jet
	* This methods returns false if there was no radiation or rescaling required
	*/
	bool reconstructDecayJet(const tShowerParticlePtr particleJetParent) const;
	//@}

	/**
	* Methods to perform the reconstruction of various types of colour
	* singlet systems
	*/
	//@{
	/**
	* Perform the reconstruction of a system with one incoming and at least one
	* outgoing particle
	*/
	void reconstructInitialFinalSystem(vector<ShowerProgenitorPtr>) const;

	/**
	* Perform the reconstruction of a system with only final-state
	* particles
	*/
	void reconstructFinalStateSystem(bool applyBoost,
	const LorentzRotation & toRest,
	const LorentzRotation & fromRest,
	vector<ShowerProgenitorPtr>) const;

	/**
	* Reconstruction of a general coloured system
	*/
	void reconstructGeneralSystem(vector<ShowerProgenitorPtr> & ShowerHardJets) const;

	/**
	* Reconstruction of a general coloured system doing
	* final-final, then initial-final and then initial-initial
	*/
	void reconstructFinalFirst(vector<ShowerProgenitorPtr> & ShowerHardJets) const;

	/**
	* Reconstruction of a general coloured system doing
	* colour parners
	*/
	void reconstructColourPartner(vector<ShowerProgenitorPtr> & ShowerHardJets) const;

	/**
	* Reconstruction based on colour singlet systems
	*/
	void reconstructColourSinglets(vector<ShowerProgenitorPtr> & ShowerHardJets,
	ShowerInteraction::Type type) const;

	/**
	* Perform the reconstruction of a system with only final-state
	* particles
	*/
	void reconstructInitialInitialSystem(bool & applyBoost,
	LorentzRotation & toRest,
	LorentzRotation & fromRest,
	vector<ShowerProgenitorPtr>) const;
	//@}

	/**
	* Methods to perform the inverse reconstruction of various types of
	* colour singlet systems
	*/
	//@{
	/**
	* Perform the inverse reconstruction of a system with only final-state
	* particles
	*/
	void deconstructFinalStateSystem(const LorentzRotation & toRest,
	const LorentzRotation & fromRest,
	HardTreePtr,
	vector<HardBranchingPtr>,
	cEvolverPtr,
	ShowerInteraction::Type) const;

	/**
	* Perform the inverse reconstruction of a system with only initial-state
	* particles
	*/
	void deconstructInitialInitialSystem(bool & applyBoost,
	LorentzRotation & toRest,
	LorentzRotation & fromRest,
	HardTreePtr,
	vector<HardBranchingPtr>,
	ShowerInteraction::Type ) const;

	/**
	* Perform the inverse reconstruction of a system with only initial-state
	* particles
	*/
	void deconstructInitialFinalSystem(HardTreePtr,
	vector<HardBranchingPtr>,
	cEvolverPtr,
	ShowerInteraction::Type ) const;

	bool deconstructGeneralSystem(HardTreePtr, cEvolverPtr,
	ShowerInteraction::Type) const;

	bool deconstructColourSinglets(HardTreePtr, cEvolverPtr,
	ShowerInteraction::Type) const;

	bool deconstructColourPartner(HardTreePtr, cEvolverPtr,
	ShowerInteraction::Type) const;
	//@}

	/**
	* Recursively treat the most off-shell paricle seperately
	* for final-final reconstruction
	*/
	void reconstructFinalFinalOffShell(JetKinVect orderedJets, Energy2 s,
	bool recursive) const;

	/**
	* Various methods for the Lorentz transforms needed to do the
	* rescalings
	*/
	//@{
	/**
	* Compute the boost to get from the the old momentum to the new
	*/
	LorentzRotation solveBoost(const double k,
	const Lorentz5Momentum & newq,
	const Lorentz5Momentum & oldp) const;

	/**
	* Compute the boost to get from the the old momentum to the new
	*/
	LorentzRotation solveBoost(const Lorentz5Momentum & newq,
	const Lorentz5Momentum & oldq) const;

	/**
	* Compute the boost to get from the the old momentum to the new
	*/
	LorentzRotation solveBoostZ(const Lorentz5Momentum & newq,
	const Lorentz5Momentum & oldq) const;

	/**
	* Recursively boost the initial-state shower
	* @param p The particle
	* @param bv The boost
	* @param parent The parent of the chain
	*/
	void boostChain(tPPtr p, const LorentzRotation & bv, tPPtr & parent) const;

	/**
	* Given a 5-momentum and a scale factor, the method returns the
	* Lorentz boost that transforms the 3-vector vec{momentum} --->
	* k*vec{momentum}. The method returns the null boost in the case no
	* solution exists. This will only work in the case where the
	* outgoing jet-momenta are parallel to the momenta of the particles
	* leaving the hard subprocess.
	*/
	Boost solveBoostBeta( const double k, const Lorentz5Momentum & newq,
	const Lorentz5Momentum & oldp);

	/**
	* Compute boost parameter along z axis to get (Ep, any perp, qp)
	* from (E, same perp, q).
	*/
	double getBeta(const double E, const double q,
	const double Ep, const double qp) const
	{return (qE-qpEp)/(sqr(qp)+sqr(E));}
	//@}

	/**
	* Methods to calculate the various scaling factors
	*/
	//@{
	/**
	* Given a vector of 5-momenta of jets, where the 3-momenta are the initial
	* ones before showering and the masses are reconstructed after the showering,
	* this method returns the overall scaling factor for the 3-momenta of the
	* vector of particles, vec{P}_i -> k * vec{P}_i, such to preserve energy-
	* momentum conservation, i.e. after the rescaling the center of mass 5-momentum
	* is equal to the one specified in input, cmMomentum.
	* The method returns 0 if such factor cannot be found.
	* @param root_s Centre-of-mass energy
	* @param jets The jets
	*/
	double solveKfactor( const Energy & root_s, const JetKinVect & jets ) const;

	/**
	* Calculate the rescaling factors for the jets in a particle decay where
	* there was initial-state radiation
	* @param mb The mass of the decaying particle
	* @param n The reference vector for the initial state radiation
	* @param pjet The momentum of the initial-state jet
	* @param jetKinematics The JetKinStruct objects for the jets
	* @param partner The colour partner
	* @param ppartner The momentum of the colour partner of the decaying particle
	* before and after radiation
	* @param k1 The rescaling parameter for the partner
	* @param k2 The rescaling parameter for the outgoing singlet
	* @param qt The transverse momentum vector
	*/
	bool solveDecayKFactor(Energy mb,
	const Lorentz5Momentum & n,
	const Lorentz5Momentum & pjet,
	const JetKinVect & jetKinematics,
	ShowerParticlePtr partner,
	Lorentz5Momentum ppartner[2],
	double & k1,
	double & k2,
	Lorentz5Momentum & qt) const;

	/**
	* Compute the momentum rescaling factor needed to invert the shower
	* @param pout The momenta of the outgoing particles
	* @param mon The on-shell masses
	* @param roots The mass of the decaying particle
	*/
	double inverseRescalingFactor(vector<Lorentz5Momentum> pout,
	vector<Energy> mon,Energy roots) const;

	/**
	* Compute the momentum rescaling factor needed to invert the shower
	* @param pout The momenta of the outgoing particles
	* @param mon The on-shell masses
	* @param roots The mass of the decaying particle
	* @param ppartner The momentum of the colour partner
	* @param mbar The mass of the decaying particle
	* @param k1 The first scaling factor
	* @param k2 The second scaling factor
	*/
	bool inverseDecayRescalingFactor(vector<Lorentz5Momentum> pout,
	vector<Energy> mon,Energy roots,
	Lorentz5Momentum ppartner, Energy mbar,
	double & k1, double & k2) const;

	/**
	* Check the rescaling conserves momentum
	* @param k The rescaling
	* @param root_s The centre-of-mass energy
	* @param jets The jets
	*/
	Energy momConsEq(const double & k, const Energy & root_s,
	const JetKinVect & jets) const;


	void findInitialBoost(const Lorentz5Momentum & pold, const Lorentz5Momentum & pnew,
	LorentzRotation & toRest, LorentzRotation & fromRest) const;
	//@}

	/**
	* Find the colour partners of a particle to identify the colour singlet
	* systems for the reconstruction.
	*/
	template<typename Value> void findPartners(Value branch,set<Value> & done,
	const set<Value> & branchings,
	vector<Value> & jets) const;

	/**
	* Add the intrinsic \f$p_T\f$ to the system if needed
	*/
	bool addIntrinsicPt(vector<ShowerProgenitorPtr>) const;

	/**
	* Apply a transform to the particle and any child, including child ShowerTree
	* objects
	* @param particle The particle
	* @param r The Lorentz transformation
	* @param match Whether or not to look at children etc
	* @param original The original particle
	*/
	void deepTransform(PPtr particle,const LorentzRotation & r,
	bool match=true,PPtr original=PPtr()) const;

	/**
	* Find the mass of a particle in the hard branching
	*/
	Energy findMass(HardBranchingPtr) const;

	/**
	* Calculate the initial-state rescaling factors
	*/
	vector<double> initialStateRescaling(double x1, double x2,
	const Lorentz5Momentum & pold,
	const vector<Lorentz5Momentum> & p,
	const vector<Lorentz5Momentum> & pq,
	const vector<Energy>& highespts) const;

	/**
	* Calculate the inverse of the initial-state rescaling factor
	*/
	vector<double> inverseInitialStateRescaling(double & x1, double & x2,
	const Lorentz5Momentum & pold,
	const vector<Lorentz5Momentum> & p,
	const vector<Lorentz5Momentum> & pq) const;

	/**
	* Find the colour singlet systems
	*/
	template<typename Value >
	typename ColourSinglet<Value>::VecType identifySystems(set<Value> jets,
	unsigned int & nnun,unsigned int & nnii,
	unsigned int & nnif,unsigned int & nnf,
	unsigned int & nni) const;

	/**
	* Combine final-state colour systems
	*/
	template<typename Value>
	void combineFinalState(vector<ColourSinglet<Value> > & systems) const;

	protected:

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

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

	protected:

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

	private:

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

	private:

	/**
	* Option for handling the reconstruction
	*/
	unsigned int _reconopt;

	/**
	* Option for the boost for initial-initial reconstruction
	*/
	unsigned int _initialBoost;

	/**
	* Option for the reconstruction of final state systems
	*/
	unsigned int _finalStateReconOption;

	/**
	* Option for the initial state reconstruction
	*/
	unsigned int _initialStateReconOption;

	/**
	* Minimum invariant mass for initial-final dipoles to allow the
	* reconstruction
	*/
	Energy _minQ;

	/**
	* The progenitor of the jet currently being reconstructed
	*/
	mutable tShowerParticlePtr _progenitor;

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

	/**
	* Current ShowerTree
	*/
	mutable tShowerTreePtr _currentTree;

	/**
	* Particles which shouldn't have their masses rescaled as
	* vector for the interface
	*/
	PDVector _noRescaleVector;

	/**
	* Particles which shouldn't have their masses rescaled as
	* set for quick access
	*/
	set<cPDPtr> _noRescale;

	/**
	* Storage of the boosts applied to enable resetting after failure
	*/
	mutable map<tPPtr,vector<LorentzRotation> > _boosts;

	/**
	* Storage of the boosts applied to enable resetting after failure
	*/
	mutable map<tShowerTreePtr,vector<LorentzRotation> > _treeBoosts;
	};

	}

	#include "QTildeReconstructor.tcc"
	#endif /* HERWIG_QTildeReconstructor_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:50 PM (1 d, 12 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805241
Default Alt Text: (254 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions