No OneTemporary
Actions

Size

179 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,1502 +1,1503 @@
	// -- C++ --
	//
	// Evolver.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2007 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 "ShowerKinematics.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Handlers/EventHandler.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ShowerTree.h"
	#include "ShowerProgenitor.h"
	#include "KinematicsReconstructor.h"
	#include "PartnerFinder.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "Herwig++/Shower/ShowerHandler.h"

	using namespace Herwig;

	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 << _limitEmissions
	<< ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
	<< _vetoes << _hardonly << _trunc_Mode << _hardEmissionMode;
	}

	void Evolver::persistentInput(PersistentIStream & is, int) {
	is >> _model >> _splittingGenerator >> _maxtry
	>> _meCorrMode >> _hardVetoMode >> _hardVetoRead >> _limitEmissions
	>> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
	>> _vetoes >> _hardonly >> _trunc_Mode >> _hardEmissionMode;
	}

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

	void Evolver::Init() {

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

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

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

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

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

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

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

	static 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 Switch<Evolver,bool> interfaceHardOnly
	("HardOnly",
	"Only generate the emission supplied by the hardest emission"
	" generator, for testing only.",
	&Evolver::_hardonly, false, false, false);
	static SwitchOption interfaceHardOnlyNo
	(interfaceHardOnly,
	"No",
	"Generate full shower",
	false);
	static SwitchOption interfaceHardOnlyYes
	(interfaceHardOnly,
	"Yes",
	"Only the hardest emission",
	true);

	static Switch<Evolver,bool> interfaceTruncMode
	("TruncatedShower", "Include the truncated shower?",
	&Evolver::_trunc_Mode, 1, false, false);
	static SwitchOption interfaceTruncMode0
	(interfaceTruncMode,"No","Truncated Shower is OFF", 0);
	static SwitchOption interfaceTruncMode1
	(interfaceTruncMode,"Yes","Truncated Shower is ON", 1);

	static Switch<Evolver,unsigned int> interfaceHardEmissionMode
	("HardEmissionMode",
	"Whether to use ME corrections or POWHEG for the hardest emission",
	&Evolver::_hardEmissionMode, 0, false, false);
	static SwitchOption interfaceHardEmissionModeMECorrection
	(interfaceHardEmissionMode,
	"MECorrection",
	"Old fashioned ME correction",
	0);
	static SwitchOption interfaceHardEmissionModePOWHEG
	(interfaceHardEmissionMode,
	"POWHEG",
	"Powheg style hard emission",
	1);

	}

	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(ShowerTreePtr hard,
	vector<ShowerProgenitorPtr> p) {
	// find out if hard partonic subprocess.
	bool isPartonic(false);
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit = _currenttree->incomingLines().begin();
	Lorentz5Momentum pcm;
	for(; cit!=hard->incomingLines().end(); ++cit) {
	pcm += cit->first->progenitor()->momentum();
	isPartonic \|= cit->first->progenitor()->coloured();
	}
	// find maximum pt from hard process, the maximum pt from all outgoing
	// coloured lines (this is simpler and more general than
	// 2stu/(s^2+t^2+u^2)). Maximum scale for scattering processes will
	// be transverse mass.
	Energy ptmax = -1.0*GeV;
	// general case calculate the scale
	if (!hardVetoXComb()) {
	// scattering process
	if(hard->isHard()) {
	// coloured incoming particles
	if (isPartonic) {
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt = hard->outgoingLines().begin();
	for(; cjt!=hard->outgoingLines().end(); ++cjt) {
	if (cjt->first->progenitor()->coloured())
	ptmax = max(ptmax,cjt->first->progenitor()->momentum().mt());
	}
	}
	if (ptmax < ZERO) ptmax = pcm.m();
	}
	// decay, incoming() is the decaying particle.
	else {
	ptmax = hard->incomingLines().begin()->first
	->progenitor()->momentum().mass();
	}
	}
	// hepeup.SCALUP is written into the lastXComb by the
	// LesHouchesReader itself - use this by user's choice.
	// Can be more general than this.
	else {
	ptmax = sqrt( ShowerHandler::currentHandler()
	->lastXCombPtr()->lastScale() );
	}
	// set maxHardPt for all progenitors. For partonic processes this
	// is now the max pt in the FS, for non-partonic processes or
	// processes with no coloured FS the invariant mass of the IS
	vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
	for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(ptmax);
	}

	void Evolver::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) {
	_hardme = HwMEBasePtr();
	// extract the matrix element
	tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(xcomb);
	if(lastXC) {
	_hardme = dynamic_ptr_cast<HwMEBasePtr>(lastXC->matrixElement());
	}
	_decayme = HwDecayerBasePtr();
	// set the current tree
	currentTree(hard);
	// zero number of emissions
	_nis = _nfs = 0;
	// extract particles to shower
	vector<ShowerProgenitorPtr> particlesToShower=setupShower(true);
	// setup the maximum scales for the shower, given by the hard process
	if (hardVetoOn()) setupMaximumScales(currentTree(), particlesToShower);
	// generate the intrinsic p_T once and for all
	generateIntrinsicpT(particlesToShower);
	// main shower loop
	unsigned int ntry(0);
	do {
	// clear results of last attempt if needed
	if(ntry!=0) {
	currentTree()->clear();
	setEvolutionPartners(true,ShowerInteraction::QCD);
	_nis = _nfs = 0;
	}
	// generate the shower
	// pick random starting point
	unsigned int istart=UseRandom::irnd(particlesToShower.size());
	unsigned int istop = particlesToShower.size();
	// loop over particles with random starting point
	for(unsigned int ix=istart;ix<=istop;++ix) {
	if(ix==particlesToShower.size()) {
	if(istart!=0) {
	istop = istart-1;
	ix=0;
	}
	else break;
	}
	// set the progenitor
	_progenitor=particlesToShower[ix];
	// initial-state
	if(!_progenitor->progenitor()->isFinalState()) {
	if(!isISRadiationON()) continue;
	// 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,
	ShowerInteraction::QCD));
	}
	// final-state
	else {
	if(!isFSRadiationON()) continue;
	// perform shower
	progenitor()->hasEmitted(startTimeLikeShower(ShowerInteraction::QCD));
	}
	}
	}
	while(!showerModel()->kinematicsReconstructor()->
	reconstructHardJets(hard,intrinsicpT())&&
	maximumTries()>++ntry);
	_hardme=HwMEBasePtr();
	_hardtree=HardTreePtr();
	if(_maxtry==ntry) throw ShowerHandler::ShowerTriesVeto(ntry);

	// the tree has now showered
	_currenttree->hasShowered(true);
	}

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

	bool Evolver::timeLikeShower(tShowerParticlePtr particle,
	ShowerInteraction::Type type) {
	// don't do anything if not needed
	if(_limitEmissions == 1 \|\| _limitEmissions == 3 \|\|
	( _limitEmissions == 2 && _nfs != 0) ) return false;
	// generate the emission
	Branching fb;
	while (true) {
	fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
	// no emission return
	if(!fb.kinematics) return false;
	// if emission OK break
	if(!timeLikeVetoed(fb,particle)) break;
	// otherwise reset scale and continue - SO IS involved in veto algorithm
	particle->setEvolutionScale(fb.kinematics->scale());
	}
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->setShowerKinematics(fb.kinematics);
	// Assign the splitting function to the emitting particle.
	// For the time being we are considering only 1->2 branching
	// Create the ShowerParticle objects for the two children of
	// the emitting particle; set the parent/child relationship
	// if same as definition create particles, otherwise create cc
	tcPDPtr pdata[2];
	for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
	if(particle->id()!=fb.ids[0]) {
	for(unsigned int ix=0;ix<2;++ix) {
	tPDPtr cc(pdata[ix]->CC());
	if(cc) pdata[ix]=cc;
	}
	}
	ShowerParticleVector theChildren;
	theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
	theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
	// update the children
	particle->showerKinematics()->updateChildren(particle, theChildren,true);
	// update the history if needed
	if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
	_currenttree->updateFinalStateShowerProduct(_progenitor,
	particle,theChildren);
	_currenttree->addFinalStateBranching(particle,theChildren);
	// update number of emissions
	++_nfs;
	if(_limitEmissions!=0) return true;
	// shower the first particle
	timeLikeShower(theChildren[0],type);
	// shower the second particle
	timeLikeShower(theChildren[1],type);
	// branching has happened
	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 \|\| _limitEmissions == 3 \|\|
	( _limitEmissions == 1 && _nis != 0 ) ) return false;
	Branching bb;
	// generate branching
	while (true) {
	bb=_splittingGenerator->chooseBackwardBranching(*particle,beam,
	_initialenhance,
	_beam,type,
	pdf,freeze);
	// return if no emission
	if(!bb.kinematics) return false;
	// if not vetoed break
	if(!spaceLikeVetoed(bb,particle)) break;
	// otherwise reset scale and continue
	particle->setEvolutionScale(bb.kinematics->scale());
	}
	// assign the splitting function and shower kinematics
	particle->setShowerKinematics(bb.kinematics);
	// For the time being we are considering only 1->2 branching
	// particles as in Sudakov form factor
	tcPDPtr part[2]={getParticleData(bb.ids[0]),
	getParticleData(bb.ids[2])};
	if(particle->id()!=bb.ids[1]) {
	if(part[0]->CC()) part[0]=part[0]->CC();
	if(part[1]->CC()) part[1]=part[1]->CC();
	}
	// Now create the actual particles, make the otherChild a final state
	// particle, while the newParent is not
	ShowerParticlePtr newParent=new_ptr(ShowerParticle(part[0],false));
	ShowerParticlePtr otherChild = new_ptr(ShowerParticle(part[1],true,true));
	ShowerParticleVector theChildren;
	theChildren.push_back(particle);
	theChildren.push_back(otherChild);
	//this updates the evolution scale
	particle->showerKinematics()->updateParent(newParent, theChildren,true);
	// update the history if needed
	_currenttree->updateInitialStateShowerProduct(_progenitor,newParent);
	_currenttree->addInitialStateBranching(particle,newParent,otherChild);
	// for the reconstruction of kinematics, parent/child
	// relationships are according to the branching process:
	// now continue the shower
	++_nis;
	bool emitted = _limitEmissions==0 ?
	spaceLikeShower(newParent,beam,type) : false;
	// now reconstruct the momentum
	if(!emitted) {
	if(_intrinsic.find(_progenitor)==_intrinsic.end()) {
	bb.kinematics->updateLast(newParent,ZERO,ZERO);
	}
	else {
	pair<Energy,double> kt=_intrinsic[_progenitor];
	bb.kinematics->updateLast(newParent,
	kt.first*cos(kt.second),
	kt.first*sin(kt.second));
	}
	}
	particle->showerKinematics()->updateChildren(newParent, theChildren,true);
	if(_limitEmissions!=0) return true;
	// perform the shower of the final-state particle
	timeLikeShower(otherChild,type);
	// return the emitted
	return true;
	}

	void Evolver::showerDecay(ShowerTreePtr decay) {
	_decayme = HwDecayerBasePtr();
	_hardme = HwMEBasePtr();
	// find the decayer
	// try the normal way if possible
	tDMPtr dm = decay->incomingLines().begin()->first->original() ->decayMode();
	if(!dm) dm = decay->incomingLines().begin()->first->copy() ->decayMode();
	if(!dm) dm = decay->incomingLines().begin()->first->progenitor()->decayMode();
	// otherwise make a string and look it up
	if(!dm) {
	string tag = decay->incomingLines().begin()->first->original()->dataPtr()->name()
	+ "->";
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) {
	if(it!=decay->outgoingLines().begin()) tag += ",";
	tag += it->first->original()->dataPtr()->name();
	}
	tag += ";";
	dm = generator()->findDecayMode(tag);
	}
	if(dm) _decayme = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
	// set the ShowerTree to be showered
	currentTree(decay);
	+ decay->applyTransforms();
	// extract particles to be shower, set scales and
	// perform hard matrix element correction
	vector<ShowerProgenitorPtr> particlesToShower=setupShower(false);
	setupMaximumScales(currentTree(), particlesToShower);
	// compute the minimum mass of the final-state
	Energy minmass(ZERO);
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(particlesToShower[ix]->progenitor()->isFinalState())
	minmass+=particlesToShower[ix]->progenitor()->mass();
	}
	// main showering loop
	unsigned int ntry(0);
	do {
	// clear results of last attempt
	if(ntry!=0) {
	currentTree()->clear();
	setEvolutionPartners(false,ShowerInteraction::QCD);
	}
	unsigned int istart=UseRandom::irnd(particlesToShower.size());
	unsigned int istop = particlesToShower.size();
	// loop over particles with random starting point
	for(unsigned int ix=istart;ix<=istop;++ix) {
	if(ix==particlesToShower.size()) {
	if(istart!=0) {
	istop = istart-1;
	ix=0;
	}
	else break;
	}
	// extract the progenitor
	progenitor(particlesToShower[ix]);
	// final-state radiation
	if(progenitor()->progenitor()->isFinalState()) {
	if(!isFSRadiationON()) continue;
	// perform shower
	progenitor()->hasEmitted(startTimeLikeShower(ShowerInteraction::QCD));
	}
	// initial-state radiation
	else {
	if(!isISRadiationON()) continue;
	// perform shower
	// set the scales correctly. The current scale is the maximum scale for
	// emission not the starting scale
	Energy maxscale=progenitor()->progenitor()->evolutionScale();
	Energy startScale=progenitor()->progenitor()->mass();
	progenitor()->progenitor()->setEvolutionScale(startScale);
	// perform the shower
	progenitor()->hasEmitted(startSpaceLikeDecayShower(maxscale,minmass,
	ShowerInteraction::QCD));
	}
	}
	}
	while(!showerModel()->kinematicsReconstructor()->reconstructDecayJets(decay)&&
	maximumTries()>++ntry);
	_decayme = HwDecayerBasePtr();
	if(maximumTries()==ntry)
	throw Exception() << "Failed to generate the shower after "
	<< ntry << " attempts in Evolver::showerDecay()"
	<< Exception::eventerror;

	// tree has now showered
	_currenttree->hasShowered(true);
	}

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

	vector<ShowerProgenitorPtr> Evolver::setupShower(bool hard) {
	// generate POWHEG hard emission if needed
	if(_hardEmissionMode==1) hardestEmission(hard);
	// set the initial colour partners
	setEvolutionPartners(hard,ShowerInteraction::QCD);
	// get the particles to be showered
	map<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator cit;
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
	// generate hard me if needed
	if(_hardEmissionMode==0) hardMatrixElementCorrection(hard);
	// get the particles to be showered
	vector<ShowerProgenitorPtr> particlesToShower;
	// incoming particles
	for(cit=currentTree()->incomingLines().begin();
	cit!=currentTree()->incomingLines().end();++cit)
	particlesToShower.push_back((*cit).first);
	assert((particlesToShower.size()==1&&!hard)\|\|
	(particlesToShower.size()==2&&hard));
	// outgoing particles
	for(cjt=currentTree()->outgoingLines().begin();
	cjt!=currentTree()->outgoingLines().end();++cjt)
	particlesToShower.push_back(((*cjt).first));
	// remake the colour partners if needed
	if(_hardEmissionMode==0 && _currenttree->hardMatrixElementCorrection()) {
	setEvolutionPartners(hard,ShowerInteraction::QCD);
	_currenttree->resetShowerProducts();
	}
	return particlesToShower;
	}

	void Evolver::setEvolutionPartners(bool hard,ShowerInteraction::Type ) {
	map<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator cit;
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
	vector<ShowerParticlePtr> particles;
	// match the particles in the ShowerTree and hardTree
	if(hardTree() && !hardTree()->connect(currentTree()))
	throw Exception() << "Can't match trees in "
	<< "Evolver::setEvolutionPartners()"
	<< Exception::eventerror;
	// sort out the colour partners
	for(cit=currentTree()->incomingLines().begin();
	cit!=currentTree()->incomingLines().end();++cit)
	particles.push_back(cit->first->progenitor());
	assert((particles.size()==1&&!hard)\|\|(particles.size()==2&&hard));
	// outgoing particles
	for(cjt=currentTree()->outgoingLines().begin();
	cjt!=currentTree()->outgoingLines().end();++cjt)
	particles.push_back(cjt->first->progenitor());
	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]->setPartner(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,ShowerInteraction::QCD,!_hardtree);
	}

	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() ) {
	return truncatedTimeLikeShower(progenitor()->progenitor(), mit->second ,type);
	}
	}
	return hardOnly() ? false :
	timeLikeShower(progenitor()->progenitor() ,type) ;
	}

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

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

	bool Evolver::timeLikeVetoed(const Branching & fb,
	ShowerParticlePtr particle) {
	// 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()) return true;
	// general vetos
	if (fb.kinematics && !_vetoes.empty()) {
	bool vetoed=false;
	for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
	v != _vetoes.end(); ++v) {
	bool test = (**v).vetoTimeLike(_progenitor,particle,fb);
	switch((**v).vetoType()) {
	case ShowerVeto::Emission:
	vetoed \|= test;
	break;
	case ShowerVeto::Shower:
	if(test) throw VetoShower();
	break;
	case ShowerVeto::Event:
	if(test) throw Veto();
	break;
	}
	}
	if(vetoed) return true;
	}
	return false;
	}

	bool Evolver::spaceLikeVetoed(const Branching & bb,ShowerParticlePtr particle) {
	// 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()) return true;

	if (!_vetoes.empty()) {
	bool vetoed=false;
	for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
	v != _vetoes.end(); ++v) {
	bool test = (**v).vetoSpaceLike(_progenitor,particle,bb);
	switch ((**v).vetoType()) {
	case ShowerVeto::Emission:
	vetoed \|= test;
	break;
	case ShowerVeto::Shower:
	if(test) throw VetoShower();
	break;
	case ShowerVeto::Event:
	if(test) throw Veto();
	break;
	}
	}
	if (vetoed) return true;
	}
	return false;
	}

	bool Evolver::spaceLikeDecayVetoed( const Branching & fb,
	ShowerParticlePtr particle ) {
	// 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()) return true;
	// general vetos
	if (!_vetoes.empty()) {
	bool vetoed=false;
	for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
	v != _vetoes.end(); ++v) {
	bool test = (**v).vetoSpaceLike(_progenitor,particle,fb);
	switch((**v).vetoType()) {
	case ShowerVeto::Emission:
	vetoed \|= test;
	break;
	case ShowerVeto::Shower:
	if(test) throw VetoShower();
	break;
	case ShowerVeto::Event:
	if(test) throw Veto();
	break;
	}
	if (vetoed) return true;
	}
	}
	return false;
	}

	void Evolver::hardestEmission(bool hard) {
	if( ( _hardme && _hardme->hasPOWHEGCorrection()) \|\|
	(_decayme && _decayme->hasPOWHEGCorrection())) {
	if(_hardme)
	_hardtree = _hardme->generateHardest( currentTree() );
	else
	_hardtree = _decayme->generateHardest( currentTree() );
	if(!_hardtree) return;
	// join up the two tree
	connectTrees(currentTree(),_hardtree,hard);
	}
	else {
	_hardtree = ShowerHandler::currentHandler()->generateCKKW(currentTree());
	}
	}

	bool Evolver::truncatedTimeLikeShower(tShowerParticlePtr particle,
	HardBranchingPtr branch,
	ShowerInteraction::Type type) {
	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->setEvolutionScale(fb.kinematics->scale());
	continue;
	}
	double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
	// only if same interaction for forced branching
	// and evolution
	if(type==branch->sudakov()->interactionType()) {
	if(zsplit < 0.5 \|\| // hardest line veto
	fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
	particle->setEvolutionScale(fb.kinematics->scale());
	continue;
	}
	}
	// pt veto
	if(fb.kinematics->pT() > progenitor()->maximumpT()) {
	particle->setEvolutionScale(fb.kinematics->scale());
	continue;
	}
	// should do base class vetos as well
	if(timeLikeVetoed(fb,particle)) {
	particle->setEvolutionScale(fb.kinematics->scale());
	continue;
	}
	break;
	}
	// if no branching set decay matrix and return
	if(!fb.kinematics) {
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin=
	branch->sudakov()->createFinalStateBranching(branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT());
	particle->setEvolutionScale(branch->scale() );
	showerKin->initialize( *particle,PPtr() );
	IdList idlist(3);
	idlist[0] = particle->id();
	idlist[1] = branch->children()[0]->branchingParticle()->id();
	idlist[2] = branch->children()[1]->branchingParticle()->id();
	fb = Branching( showerKin, idlist, branch->sudakov() );
	// Assign the shower kinematics to the emitting particle.
	particle->setShowerKinematics( fb.kinematics );
	// Assign the splitting function to the emitting particle.
	// For the time being we are considering only 1->2 branching
	// Create the ShowerParticle objects for the two children of
	// the emitting particle; set the parent/child relationship
	// if same as definition create particles, otherwise create cc
	ShowerParticleVector theChildren;
	theChildren.push_back(new_ptr(ShowerParticle(branch->children()[0]->
	branchingParticle()->dataPtr(),true)));
	theChildren.push_back(new_ptr(ShowerParticle(branch->children()[1]->
	branchingParticle()->dataPtr(),true)));
	particle->showerKinematics()->
	updateChildren(particle, theChildren,type==branch->sudakov()->interactionType());
	// update the history if needed
	if(particle==currentTree()->getFinalStateShowerProduct(progenitor()))
	currentTree()->updateFinalStateShowerProduct(progenitor(),
	particle,theChildren);
	currentTree()->addFinalStateBranching(particle,theChildren);
	// shower the first particle
	if( branch->children()[0]->children().empty() ) {
	if( ! hardOnly() )
	timeLikeShower(theChildren[0],type);
	}
	else {
	truncatedTimeLikeShower( theChildren[0],branch->children()[0],type);
	}
	// shower the second particle
	if( branch->children()[1]->children().empty() ) {
	if( ! hardOnly() )
	timeLikeShower( theChildren[1] , type);
	}
	else {
	truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type);
	}
	return true;
	}
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->setShowerKinematics(fb.kinematics);
	// Assign the splitting function to the emitting particle.
	// For the time being we are considering only 1->2 branching
	// Create the ShowerParticle objects for the two children of
	// the emitting particle; set the parent/child relationship
	// if same as definition create particles, otherwise create cc
	ShowerParticleVector theChildren;
	theChildren.push_back( new_ptr( ShowerParticle( pdata[0], true ) ) );
	theChildren.push_back( new_ptr( ShowerParticle( pdata[1], true ) ) );
	particle->showerKinematics()->
	updateChildren( particle, theChildren , true);
	// update the history if needed
	if( particle == currentTree()->getFinalStateShowerProduct( progenitor() ) )
	currentTree()->updateFinalStateShowerProduct( progenitor(),
	particle, theChildren );
	currentTree()->addFinalStateBranching( particle, theChildren );
	// shower the first particle
	if( iout == 1 ) truncatedTimeLikeShower( theChildren[0], branch , type );
	else timeLikeShower( theChildren[0] , type);
	// shower the second particle
	if( iout == 2 ) truncatedTimeLikeShower( theChildren[1], branch , type );
	else timeLikeShower( theChildren[1] , type);
	// branching has happened
	return true;
	}

	bool Evolver::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
	HardBranchingPtr branch,
	ShowerInteraction::Type type) {
	tcPDFPtr pdf;
	if(ShowerHandler::currentHandler()->firstPDF().particle() == beamParticle())
	pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
	if(ShowerHandler::currentHandler()->secondPDF().particle() == beamParticle())
	pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
	Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
	Branching bb;
	// generate branching
	tcPDPtr part[2];
	while (true) {
	if( !isTruncatedShowerON() \|\| hardOnly() ) break;
	bb = splittingGenerator()->chooseBackwardBranching( *particle,
	beam, 1., beamParticle(),
	type , pdf,freeze);
	if( !bb.kinematics \|\| bb.kinematics->scale() < branch->scale() ) {
	bb = Branching();
	break;
	}
	// particles as in Sudakov form factor
	part[0] = getParticleData( bb.ids[0] );
	part[1] = getParticleData( bb.ids[2] );

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

	bool Evolver::
	truncatedSpaceLikeDecayShower(tShowerParticlePtr particle, Energy maxscale,
	Energy minmass, HardBranchingPtr branch,
	ShowerInteraction::Type type) {
	Branching fb;
	unsigned int iout=0;
	tcPDPtr pdata[2];
	while (true) {
	// no truncated shower break
	if(!isTruncatedShowerON()\|\|hardOnly()) break;
	fb=splittingGenerator()->chooseDecayBranching(*particle,maxscale,minmass,1.,type);
	// return if no radiation
	if(!fb.kinematics) break;
	// check haven't evolved too far
	if(fb.kinematics->scale() < branch->scale()) {
	fb=Branching();
	break;
	}
	// get the particle data objects
	for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
	if(particle->id()!=fb.ids[0]) {
	for(unsigned int ix=0;ix<2;++ix) {
	tPDPtr cc(pdata[ix]->CC());
	if(cc) pdata[ix]=cc;
	}
	}
	// find the truncated line
	iout=0;
	if(pdata[0]->id()!=pdata[1]->id()) {
	if(pdata[0]->id()==particle->id()) iout=1;
	else if (pdata[1]->id()==particle->id()) iout=2;
	}
	else if(pdata[0]->id()==particle->id()) {
	if(fb.kinematics->z()>0.5) iout=1;
	else iout=2;
	}
	// apply the vetos for the truncated shower
	// no flavour changing branchings
	if(iout==0) {
	particle->setEvolutionScale(fb.kinematics->scale());
	continue;
	}
	double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
	if(type==branch->sudakov()->interactionType()) {
	if(zsplit < 0.5 \|\| // hardest line veto
	fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
	particle->setEvolutionScale(fb.kinematics->scale());
	continue;
	}
	}
	// pt veto
	if(fb.kinematics->pT() > progenitor()->maximumpT()) {
	particle->setEvolutionScale(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->setEvolutionScale(fb.kinematics->scale());
	}
	// if no branching set decay matrix and return
	if(!fb.kinematics) {
	// construct the kinematics for the hard emission
	ShoKinPtr showerKin=
	branch->sudakov()->createDecayBranching(branch->scale(),
	branch->children()[0]->z(),
	branch->phi(),
	branch->children()[0]->pT());
	particle->setEvolutionScale(branch->scale() );
	showerKin->initialize( *particle,PPtr() );
	IdList idlist(3);
	idlist[0] = particle->id();
	idlist[1] = branch->children()[0]->branchingParticle()->id();
	idlist[2] = branch->children()[1]->branchingParticle()->id();
	fb = Branching( showerKin, idlist, branch->sudakov() );
	// Assign the shower kinematics to the emitting particle.
	particle->setShowerKinematics( fb.kinematics );
	// Assign the splitting function to the emitting particle.
	// For the time being we are considering only 1->2 branching
	// Create the ShowerParticle objects for the two children of
	// the emitting particle; set the parent/child relationship
	// if same as definition create particles, otherwise create cc
	ShowerParticleVector theChildren;
	theChildren.push_back(new_ptr(ShowerParticle(branch->children()[0]->
	branchingParticle()->dataPtr(),true)));
	theChildren.push_back(new_ptr(ShowerParticle(branch->children()[1]->
	branchingParticle()->dataPtr(),true)));
	particle->showerKinematics()->
	updateChildren(particle, theChildren,
	type==branch->sudakov()->interactionType());
	if(theChildren[0]->id()==particle->id()) {
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[0]);
	currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
	// shower the space-like particle
	if( branch->children()[0]->children().empty() ) {
	if( ! hardOnly() ) spaceLikeDecayShower(theChildren[0],maxscale,minmass,type);
	}
	else {
	truncatedSpaceLikeDecayShower( theChildren[0],maxscale,minmass,
	branch->children()[0],type);
	}
	// shower the second particle
	if( branch->children()[1]->children().empty() ) {
	if( ! hardOnly() ) timeLikeShower( theChildren[1] , type);
	}
	else {
	truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type);
	}
	}
	else {
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[1]);
	currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
	// shower the space-like particle
	if( branch->children()[1]->children().empty() ) {
	if( ! hardOnly() ) spaceLikeDecayShower(theChildren[1],maxscale,minmass,type);
	}
	else {
	truncatedSpaceLikeDecayShower( theChildren[1],maxscale,minmass,
	branch->children()[1],type);
	}
	// shower the second particle
	if( branch->children()[0]->children().empty() ) {
	if( ! hardOnly() ) timeLikeShower( theChildren[0] , type);
	}
	else {
	truncatedTimeLikeShower( theChildren[0],branch->children()[0] ,type);
	}
	}
	return true;
	}
	// has emitted
	// Assign the shower kinematics to the emitting particle.
	particle->setShowerKinematics(fb.kinematics);
	// For the time being we are considering only 1->2 branching
	// Create the ShowerParticle objects for the two children of
	// the emitting particle; set the parent/child relationship
	// if same as definition create particles, otherwise create cc
	ShowerParticleVector theChildren;
	theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
	theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
	particle->showerKinematics()->updateChildren(particle, theChildren,true);
	// In the case of splittings which involves coloured particles,
	// set properly the colour flow of the branching.
	// update the history if needed
	currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[0]);
	currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
	// shower the first particle
	truncatedSpaceLikeDecayShower(theChildren[0],maxscale,minmass,branch,type);
	// shower the second particle
	timeLikeShower(theChildren[1],type);
	// branching has happened
	return true;
	}

	void Evolver::connectTrees(ShowerTreePtr showerTree, HardTreePtr hardTree, bool hard )const {
	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) {
	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()) {
	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::generateHardest()"
	<< 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(),true);
	// inverse reconstruction
	if(hard)
	showerModel()->kinematicsReconstructor()->
	deconstructHardJets(hardTree,ShowerHandler::currentHandler()->evolver(),
	hardTree->interaction());
	else
	showerModel()->kinematicsReconstructor()->
	deconstructDecayJets(hardTree,ShowerHandler::currentHandler()->evolver(),
	hardTree->interaction());
	// now reset the momenta of the showering particles
	vector<ShowerProgenitorPtr> particlesToShower;
	for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit=showerTree->incomingLines().begin();
	cit!=showerTree->incomingLines().end();++cit )
	particlesToShower.push_back(cit->first);
	// extract the showering particles
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cit=showerTree->outgoingLines().begin();
	cit!=showerTree->outgoingLines().end();++cit )
	particlesToShower.push_back(cit->first);
	// match them
	vector<bool> matched(particlesToShower.size(),false);
	for(set<HardBranchingPtr>::const_iterator cit=hardTree->branchings().begin();
	cit!=hardTree->branchings().end();++cit) {
	Energy2 dmin( 1e30*GeV2 );
	int iloc(-1);
	for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
	if(matched[ix]) continue;
	if( (**cit).branchingParticle()->id() != particlesToShower[ix]->progenitor()->id() ) continue;
	if( (**cit).branchingParticle()->isFinalState() !=
	particlesToShower[ix]->progenitor()->isFinalState() ) continue;
	Energy2 dtest =
	sqr( particlesToShower[ix]->progenitor()->momentum().x() - (**cit).showerMomentum().x() ) +
	sqr( particlesToShower[ix]->progenitor()->momentum().y() - (**cit).showerMomentum().y() ) +
	sqr( particlesToShower[ix]->progenitor()->momentum().z() - (**cit).showerMomentum().z() ) +
	sqr( particlesToShower[ix]->progenitor()->momentum().t() - (**cit).showerMomentum().t() );
	// add mass difference for identical particles (e.g. Z0 Z0 production)
	dtest += 1e10*sqr(particlesToShower[ix]->progenitor()->momentum().m()-
	(**cit).showerMomentum().m());
	if( dtest < dmin ) {
	iloc = ix;
	dmin = dtest;
	}
	}
	if(iloc<0) throw Exception() << "Failed to match shower and hard trees in Evolver::hardestEmission"
	<< Exception::eventerror;
	particlesToShower[iloc]->progenitor()->set5Momentum((**cit).showerMomentum());
	matched[iloc] = true;
	}
	// correction boosts for daughter trees
	for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	tit = showerTree->treelinks().begin();
	tit != showerTree->treelinks().end();++tit) {
	ShowerTreePtr decayTree = tit->first;
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
	cit = decayTree->incomingLines().begin();
	// reset the momentum of the decay particle
	Lorentz5Momentum oldMomentum = cit->first->progenitor()->momentum();
	Lorentz5Momentum newMomentum = tit->second.second->momentum();
	LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
	boost.boost (-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
	- decayTree->transform(boost);
	+ decayTree->transform(boost,true);
	}
	}
	diff --git a/Shower/Base/ShowerTree.cc b/Shower/Base/ShowerTree.cc
	--- a/Shower/Base/ShowerTree.cc
	+++ b/Shower/Base/ShowerTree.cc
	@@ -1,874 +1,914 @@
	// -- C++ --
	//
	// ShowerTree.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2007 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.
	//
	#include "ShowerProgenitor.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ShowerTree.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "ThePEG/Handlers/EventHandler.h"
	#include "ThePEG/Handlers/XComb.h"
	#include "KinematicsReconstructor.h"
	#include <cassert>
	#include "ThePEG/Repository/CurrentGenerator.h"

	using namespace Herwig;
	using namespace ThePEG;

	set<long> ShowerTree::_decayInShower = set<long>();

	namespace {
	void findBeam(tPPtr & beam, PPtr incoming) {
	while(!beam->children().empty()) {
	bool found=false;
	for(unsigned int ix=0;ix<beam->children().size();++ix) {
	if(beam->children()[ix]==incoming) {
	found = true;
	break;
	}
	}
	if(found) break;
	beam = beam->children()[0];
	}
	}
	}

	// constructor from hard process
	ShowerTree::ShowerTree(const PPair incoming, const ParticleVector & out,
	ShowerDecayMap& decay)
	: _hardMECorrection(false), _wasHard(true),
	_parent(), _hasShowered(false) {
	tPPair beam = CurrentGenerator::current().currentEvent()->incoming();
	findBeam(beam.first ,incoming.first );
	findBeam(beam.second,incoming.second);
	_incoming = incoming;
	double x1(_incoming.first ->momentum().rho()/beam.first ->momentum().rho());
	double x2(_incoming.second->momentum().rho()/beam.second->momentum().rho());
	// must have two incoming particles
	assert(_incoming.first && _incoming.second);
	// set the parent tree
	_parent=ShowerTreePtr();
	// temporary vectors to contain all the particles before insertion into
	// the data structure
	vector<PPtr> original,copy;
	vector<ShowerParticlePtr> shower;
	// create copies of ThePEG particles for the incoming particles
	original.push_back(_incoming.first);
	copy.push_back(new_ptr(Particle(*_incoming.first)));
	original.push_back(_incoming.second);
	copy.push_back(new_ptr(Particle(*_incoming.second)));
	// and same for outgoing
	map<PPtr,ShowerTreePtr> trees;
	for (ParticleVector::const_iterator it = out.begin();
	it != out.end(); ++it) {
	// if decayed or should be decayed in shower make the tree
	PPtr orig = *it;
	if(!orig->children().empty()\|\|
	(decaysInShower(orig->id())&&!orig->dataPtr()->stable())) {
	ShowerTreePtr newtree=new_ptr(ShowerTree(orig,decay));
	newtree->setParents();
	trees.insert(make_pair(orig,newtree));
	Energy width=orig->dataPtr()->generateWidth(orig->mass());
	decay.insert(make_pair(width,newtree));
	}
	original.push_back(orig);
	copy.push_back(new_ptr(Particle(*orig)));
	}
	// colour isolate the hard process
	colourIsolate(original,copy);
	// now create the Shower particles
	// create ShowerParticles for the incoming particles
	assert(original.size() == copy.size());
	for(unsigned int ix=0;ix<original.size();++ix) {
	ShowerParticlePtr temp=new_ptr(ShowerParticle(*copy[ix],1,ix>=2));
	fixColour(temp);
	// incoming
	if(ix<2) {
	temp->x(ix==0 ? x1 : x2);
	_incomingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[ix],
	copy[ix],temp)),
	temp));
	_backward.insert(temp);
	}
	// outgoing
	else {
	_outgoingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[ix],
	copy[ix],temp)),
	temp));
	_forward.insert(temp);
	}
	}
	// set up the map of daughter trees
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mit;
	for(mit=_outgoingLines.begin();mit!=_outgoingLines.end();++mit) {
	map<PPtr,ShowerTreePtr>::const_iterator tit=trees.find(mit->first->original());
	if(tit!=trees.end())
	_treelinks.insert(make_pair(tit->second,
	make_pair(mit->first,mit->first->progenitor())));
	}
	}

	ShowerTree::ShowerTree(PPtr in,
	ShowerDecayMap& decay)
	: _hardMECorrection(false), _wasHard(false), _hasShowered(false) {
	// there must be an incoming particle
	assert(in);
	// temporary vectors to contain all the particles before insertion into
	// the data structure
	vector<PPtr> original,copy;
	// insert place holder for incoming particle
	original.push_back(in);
	copy.push_back(PPtr());
	// we need to deal with the decay products if decayed
	map<PPtr,ShowerTreePtr> trees;
	if(!in->children().empty()) {
	ParticleVector children=in->children();
	for(unsigned int ix=0;ix<children.size();++ix) {
	// if decayed or should be decayed in shower make the tree
	PPtr orig=children[ix];
	in->abandonChild(orig);
	if(!orig->children().empty()\|\|
	(decaysInShower(orig->id())&&!orig->dataPtr()->stable())) {
	ShowerTreePtr newtree=new_ptr(ShowerTree(orig,decay));
	trees.insert(make_pair(orig,newtree));
	Energy width=orig->dataPtr()->generateWidth(orig->mass());
	decay.insert(make_pair(width,newtree));
	newtree->setParents();
	newtree->_parent=this;
	}
	original.push_back(orig);
	copy.push_back(new_ptr(Particle(*orig)));
	}
	}
	// create the incoming particle
	copy[0] = new_ptr(Particle(*in));
	// isolate the colour
	colourIsolate(original,copy);
	// create the parent
	ShowerParticlePtr sparent(new_ptr(ShowerParticle(*copy[0],2,false)));
	fixColour(sparent);
	_incomingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[0],copy[0],sparent))
	,sparent));
	// return if not decayed
	if(original.size()==1) return;
	// create the children
	assert(copy.size() == original.size());
	for (unsigned int ix=1;ix<original.size();++ix) {
	ShowerParticlePtr stemp= new_ptr(ShowerParticle(*copy[ix],2,true));
	fixColour(stemp);
	_outgoingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[ix],copy[ix],
	stemp)),
	stemp));
	_forward.insert(stemp);
	}
	// set up the map of daughter trees
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mit;
	for(mit=_outgoingLines.begin();mit!=_outgoingLines.end();++mit) {
	map<PPtr,ShowerTreePtr>::const_iterator tit=trees.find(mit->first->original());
	if(tit!=trees.end())
	_treelinks.insert(make_pair(tit->second,
	make_pair(mit->first,mit->first->progenitor())));
	}
	}

	void ShowerTree::updateFinalStateShowerProduct(ShowerProgenitorPtr progenitor,
	ShowerParticlePtr parent,
	const ShowerParticleVector & children) {
	assert(children.size()==2);
	bool matches[2];
	for(unsigned int ix=0;ix<2;++ix) {
	matches[ix] = children[ix]->id()==progenitor->id();
	}
	ShowerParticlePtr newpart;
	if(matches[0]&&matches[1]) {
	if(parent->showerKinematics()->z()>0.5) newpart=children[0];
	else newpart=children[1];
	}
	else if(matches[0]) newpart=children[0];
	else if(matches[1]) newpart=children[1];
	_outgoingLines[progenitor]=newpart;
	}

	void ShowerTree::updateInitialStateShowerProduct(ShowerProgenitorPtr progenitor,
	ShowerParticlePtr newParent) {
	_incomingLines[progenitor]=newParent;
	}

	void ShowerTree::colourIsolate(const vector<PPtr> & original,
	const vector<PPtr> & copy) {
	// vectors must have same size
	assert(original.size()==copy.size());
	// create a temporary map with all the particles to make looping easier
	vector<PPair> particles;
	particles.reserve(original.size());
	for(unsigned int ix=0;ix<original.size();++ix)
	particles.push_back(make_pair(copy[ix],original[ix]));
	// reset the colour of the copies
	vector<PPair>::const_iterator cit,cjt;
	for(cit=particles.begin();cit!=particles.end();++cit)
	if((cit).first->colourInfo()) (cit).first->colourInfo(new_ptr(ColourBase()));
	map<tColinePtr,tColinePtr> cmap;
	// make the colour connections of the copies
	for(cit=particles.begin();cit!=particles.end();++cit) {
	ColinePtr c1,newline;
	// if particle has a colour line
	if((cit).second->colourLine()&&!(cit).first->colourLine()) {
	c1=(*cit).second->colourLine();
	newline=ColourLine::create((*cit).first);
	cmap[c1]=newline;
	for(cjt=particles.begin();cjt!=particles.end();++cjt) {
	if(cjt==cit) continue;
	if((*cjt).second->colourLine()==c1)
	newline->addColoured((*cjt).first);
	else if((*cjt).second->antiColourLine()==c1)
	newline->addColoured((*cjt).first,true);
	}
	}
	// if anticolour line
	if((cit).second->antiColourLine()&&!(cit).first->antiColourLine()) {
	c1=(*cit).second->antiColourLine();
	newline=ColourLine::create((*cit).first,true);
	cmap[c1]=newline;
	for(cjt=particles.begin();cjt!=particles.end();++cjt) {
	if(cjt==cit) continue;
	if((*cjt).second->colourLine()==c1)
	newline->addColoured((*cjt).first);
	else if((*cjt).second->antiColourLine()==c1)
	newline->addColoured((*cjt).first,true);
	}
	}
	}
	// sort out sinks and sources
	for(cit=particles.begin();cit!=particles.end();++cit) {
	tColinePtr cline[2];
	tColinePair cpair;
	for(unsigned int ix=0;ix<4;++ix) {
	cline[0] = ix<2 ? cit->second->colourLine() : cit->second->antiColourLine();
	cline[1] = ix<2 ? cit->first ->colourLine() : cit->first ->antiColourLine();
	if(cline[0]) {
	switch (ix) {
	case 0: case 2:
	cpair = cline[0]->sinkNeighbours();
	break;
	case 1: case 3:
	cpair = cline[0]->sourceNeighbours();
	break;
	};
	}
	else {
	cpair = make_pair(tColinePtr(),tColinePtr());
	}
	if(cline[0]&&cpair.first) {
	map<tColinePtr,tColinePtr>::const_iterator
	mit[2] = {cmap.find(cpair.first),cmap.find(cpair.second)};
	if(mit[0]!=cmap.end()&&mit[1]!=cmap.end()) {
	if(ix==0\|\|ix==2) {
	cline[1]->setSinkNeighbours(mit[0]->second,mit[1]->second);
	}
	else {
	cline[1]->setSourceNeighbours(mit[0]->second,mit[1]->second);
	}
	}
	}
	}
	}
	}

	void ShowerTree::insertHard(StepPtr pstep, bool ISR, bool) {
	assert(_incomingLines.size()==2);
	_colour.clear();
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
	// construct the map of colour lines for hard process
	for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) {
	if(!cit->first->perturbative()) continue;
	if((*cit).first->copy()->colourLine())
	_colour.insert(make_pair((*cit).first->copy()->colourLine(),
	(*cit).first->original()->colourLine()));
	if((*cit).first->copy()->antiColourLine())
	_colour.insert(make_pair((*cit).first->copy()->antiColourLine(),
	(*cit).first->original()->antiColourLine()));
	}
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
	for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) {
	if(!cjt->first->perturbative()) continue;
	if((*cjt).first->copy()->colourLine())
	_colour.insert(make_pair((*cjt).first->copy()->colourLine(),
	(*cjt).first->original()->colourLine()));
	if((*cjt).first->copy()->antiColourLine())
	_colour.insert(make_pair((*cjt).first->copy()->antiColourLine(),
	(*cjt).first->original()->antiColourLine()));
	}
	// initial-state radiation
	if(ISR) {
	for(cit=incomingLines().begin();cit!=incomingLines().end();++cit) {
	ShowerParticlePtr init=(*cit).first->progenitor();
	assert(init->getThePEGBase());
	PPtr original = (*cit).first->original();
	if(original->parents().empty()) continue;
	PPtr hadron= original->parents()[0];
	assert(!original->children().empty());
	PPtr copy=cit->first->copy();
	ParticleVector intermediates=original->children();
	for(unsigned int ix=0;ix<intermediates.size();++ix) {
	init->abandonChild(intermediates[ix]);
	copy->abandonChild(intermediates[ix]);
	}
	// if not from a matrix element correction
	if(cit->first->perturbative()) {
	// break mother/daugther relations
	init->addChild(original);
	hadron->abandonChild(original);
	// if particle showers add shower
	if(cit->first->hasEmitted()) {
	addInitialStateShower(init,hadron,pstep,false);
	}
	// no showering for this particle
	else {
	updateColour(init);
	hadron->addChild(init);
	pstep->addIntermediate(init);
	}
	}
	// from matrix element correction
	else {
	// break mother/daugther relations
	hadron->abandonChild(original);
	copy->addChild(original);
	updateColour(copy);
	init->addChild(copy);
	pstep->addIntermediate(copy);
	// if particle showers add shower
	if(cit->first->hasEmitted()) {
	addInitialStateShower(init,hadron,pstep,false);
	}
	// no showering for this particle
	else {
	updateColour(init);
	hadron->addChild(init);
	pstep->addIntermediate(init);
	}
	}
	}
	}
	else {
	for(cit=incomingLines().begin();cit!=incomingLines().end();++cit) {
	ShowerParticlePtr init=(*cit).first->progenitor();
	assert(init->getThePEGBase());
	PPtr original = (*cit).first->original();
	if(original->parents().empty()) continue;
	PPtr hadron= original->parents()[0];
	assert(!original->children().empty());
	PPtr copy=cit->first->copy();
	ParticleVector intermediates=original->children();
	for(unsigned int ix=0;ix<intermediates.size();++ix) {
	init->abandonChild(intermediates[ix]);
	copy->abandonChild(intermediates[ix]);
	}
	// break mother/daugther relations
	init->addChild(original);
	hadron->abandonChild(original);
	// no showering for this particle
	updateColour(init);
	hadron->addChild(init);
	pstep->addIntermediate(init);
	}
	}
	// final-state radiation
	for(cjt=outgoingLines().begin();cjt!=outgoingLines().end();++cjt) {
	ShowerParticlePtr init=(*cjt).first->progenitor();
	assert(init->getThePEGBase());
	// if not from a matrix element correction
	if(cjt->first->perturbative()) {
	// register the shower particle as a
	// copy of the one from the hard process
	tParticleVector parents=init->parents();
	for(unsigned int ix=0;ix<parents.size();++ix)
	parents[ix]->abandonChild(init);
	(*cjt).first->original()->addChild(init);
	pstep->addDecayProduct(init);
	}
	// from a matrix element correction
	else {
	if(cjt->first->original()==_incoming.first\|\|
	cjt->first->original()==_incoming.second) {
	updateColour((*cjt).first->copy());
	(*cjt).first->original()->parents()[0]->
	addChild((*cjt).first->copy());
	pstep->addDecayProduct((*cjt).first->copy());
	(*cjt).first->copy()->addChild(init);
	pstep->addDecayProduct(init);
	}
	else {
	updateColour((*cjt).first->copy());
	(cjt).first->original()->addChild((cjt).first->copy());
	pstep->addDecayProduct((*cjt).first->copy());
	(*cjt).first->copy()->addChild(init);
	pstep->addDecayProduct(init);
	}
	}
	updateColour(init);
	// insert shower products
	addFinalStateShower(init,pstep);
	}
	_colour.clear();
	}

	void ShowerTree::addFinalStateShower(PPtr p, StepPtr s) {
	if(p->children().empty()) return;
	ParticleVector::const_iterator child;
	for(child=p->children().begin(); child != p->children().end(); ++child) {
	updateColour(*child);
	s->addDecayProduct(*child);
	addFinalStateShower(*child,s);
	}
	}

	void ShowerTree::updateColour(PPtr particle) {
	// if attached to a colour line
	if(particle->colourLine()) {
	bool reset=false;
	// if colour line from hard process reconnect
	if(_colour.find(particle->colourLine())!=_colour.end()) {
	ColinePtr c1=particle->colourLine();
	c1->removeColoured(particle);
	_colour[c1]->addColoured(particle);
	reset=true;
	}
	// ensure properly connected to the line
	if(!reset) {
	ColinePtr c1=particle->colourLine();
	c1->removeColoured(particle);
	c1->addColoured(particle);
	}
	}
	// if attached to an anticolour line
	if(particle->antiColourLine()) {
	bool reset=false;
	// if anti colour line from hard process reconnect
	if(_colour.find(particle->antiColourLine())!=_colour.end()) {
	ColinePtr c1=particle->antiColourLine();
	c1->removeColoured(particle,true);
	_colour[c1]->addColoured(particle,true);
	reset=true;
	}
	if(!reset) {
	ColinePtr c1=particle->antiColourLine();
	c1->removeColoured(particle,true);
	c1->addColoured(particle,true);
	}
	}
	}

	void ShowerTree::addInitialStateShower(PPtr p, PPtr hadron,
	StepPtr s, bool addchildren) {
	// Each parton here should only have one parent
	if(!p->parents().empty()) {
	if(p->parents().size()!=1)
	throw Exception() << "Particle must only have one parent in ShowerTree"
	<< "::addInitialStateShower" << Exception::runerror;
	addInitialStateShower(p->parents()[0],hadron,s);
	}
	else {
	hadron->addChild(p);
	s->addIntermediate(p);
	}
	updateColour(p);
	ParticleVector::const_iterator child;
	// if not adding children return
	if(!addchildren) return;
	// add children
	for(child = p->children().begin(); child != p->children().end(); ++child) {
	// if a final-state particle update the colour
	ShowerParticlePtr schild =
	dynamic_ptr_cast<ShowerParticlePtr>(*child);
	if(schild && schild->isFinalState()) updateColour(*child);
	// if there are grandchildren of p
	if(!(*child)->children().empty()) {
	// Add child as intermediate
	s->addIntermediate(*child);
	// If child is shower particle and final-state, add children
	if(schild && schild->isFinalState()) addFinalStateShower(schild,s);
	}
	else
	s->addDecayProduct(*child);
	}
	}

	void ShowerTree::decay(ShowerDecayMap & decay) {
	// must be one incoming particle
	assert(_incomingLines.size()==1);
	// if already decayed return
	if(!_outgoingLines.empty()) return;
	// otherwise decay it
	// now we need to replace the particle with a new copy after the shower
	// find particle after the shower
	ShowerParticlePtr newparent=_parent->_treelinks[this].second;
	// now make the new progenitor
	vector<PPtr> original,copy;
	original.push_back(newparent);
	copy.push_back(new_ptr(Particle(*newparent)));
	// reisolate the colour
	colourIsolate(original,copy);
	// make the new progenitor
	ShowerParticlePtr stemp=new_ptr(ShowerParticle(*copy[0],2,false));
	fixColour(stemp);
	ShowerProgenitorPtr newprog=new_ptr(ShowerProgenitor(original[0],copy[0],stemp));
	_incomingLines.clear();
	_incomingLines.insert(make_pair(newprog,stemp));
	// now we need to decay the copy
	PPtr parent=copy[0];
	unsigned int ntry = 0;
	while (true) {
	// exit if fails
	if (++ntry>=200)
	throw Exception() << "Failed to perform decay in ShowerTree::decay()"
	<< " after " << 200
	<< " attempts for " << parent->PDGName()
	<< Exception::eventerror;
	// select decay mode
	tDMPtr dm(parent->data().selectMode(*parent));
	if(!dm)
	throw Exception() << "Failed to select decay mode in ShowerTree::decay()"
	<< "for " << newparent->PDGName()
	<< Exception::eventerror;
	if(!dm->decayer())
	throw Exception() << "No Decayer for selected decay mode "
	<< " in ShowerTree::decay()"
	<< Exception::runerror;
	// start of try block
	try {
	ParticleVector children = dm->decayer()->decay(dm, parent);
	// if no children have another go
	if(children.empty()) continue;
	// set up parent
	parent->decayMode(dm);
	// add children
	for (unsigned int i = 0, N = children.size(); i < N; ++i ) {
	children[i]->setLabVertex(parent->labDecayVertex());
	parent->addChild(children[i]);
	parent->scale(ZERO);
	}
	// if succeeded break out of loop
	break;
	}
	catch(KinematicsReconstructionVeto) {}
	}
	// insert the trees from the children
	ParticleVector children=parent->children();
	map<PPtr,ShowerTreePtr> trees;
	for(unsigned int ix=0;ix<children.size();++ix) {
	PPtr orig=children[ix];
	parent->abandonChild(orig);
	// if particle has children or decays in shower
	if(!orig->children().empty()\|\|
	(decaysInShower(orig->id())&&!orig->dataPtr()->stable())) {
	ShowerTreePtr newtree=new_ptr(ShowerTree(orig,decay));
	trees.insert(make_pair(orig,newtree));
	Energy width=orig->dataPtr()->generateWidth(orig->mass());
	decay.insert(make_pair(width,newtree));
	}
	// now create the shower progenitors
	PPtr ncopy=new_ptr(Particle(*orig));
	//copy[0]->addChild(ncopy);
	ShowerParticlePtr nshow=new_ptr(ShowerParticle(*ncopy,2,true));
	fixColour(nshow);
	ShowerProgenitorPtr prog=new_ptr(ShowerProgenitor(children[ix],
	ncopy,nshow));
	_outgoingLines.insert(make_pair(prog,nshow));
	}
	// set up the map of daughter trees
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mit;
	for(mit=_outgoingLines.begin();mit!=_outgoingLines.end();++mit) {
	map<PPtr,ShowerTreePtr>::const_iterator tit=trees.find(mit->first->original());
	if(tit!=trees.end()) {
	_treelinks.insert(make_pair(tit->second,
	make_pair(mit->first,
	mit->first->progenitor())));
	tit->second->_parent=this;
	}
	}
	}

	void ShowerTree::insertDecay(StepPtr pstep,bool ISR, bool) {
	assert(_incomingLines.size()==1);
	_colour.clear();
	// find final particle from previous tree
	PPtr final;
	if(_parent&&!_parent->_treelinks.empty())
	final = _parent->_treelinks[this].second;
	else
	final=_incomingLines.begin()->first->original();
	// construct the map of colour lines
	PPtr copy=_incomingLines.begin()->first->copy();
	if(copy->colourLine())
	_colour.insert(make_pair(copy->colourLine(),final->colourLine()));
	if(copy->antiColourLine())
	_colour.insert(make_pair(copy->antiColourLine(),final->antiColourLine()));
	// initial-state radiation
	if(ISR&&!_incomingLines.begin()->first->progenitor()->children().empty()) {
	ShowerParticlePtr init=_incomingLines.begin()->first->progenitor();
	updateColour(init);
	final->addChild(init);
	pstep->addDecayProduct(init);
	// insert shower products
	addFinalStateShower(init,pstep);
	// sort out colour
	final=_incomingLines.begin()->second;
	_colour.clear();
	if(copy->colourLine())
	_colour.insert(make_pair(copy->colourLine(),final->colourLine()));
	if(copy->antiColourLine())
	_colour.insert(make_pair(copy->antiColourLine(),final->antiColourLine()));
	}
	// get the decaying particles
	// make the copy
	tColinePair cline=make_pair(copy->colourLine(),copy->antiColourLine());
	updateColour(copy);
	// sort out sinks and sources if needed
	if(cline.first) {
	if(cline.first->sourceNeighbours().first) {
	copy->colourLine()->setSourceNeighbours(cline.first->sourceNeighbours().first,
	cline.first->sourceNeighbours().second);
	}
	else if (cline.first->sinkNeighbours().first) {
	copy->colourLine()->setSinkNeighbours(cline.first->sinkNeighbours().first,
	cline.first->sinkNeighbours().second);
	}
	}
	if(cline.second) {
	if(cline.second->sourceNeighbours().first) {
	copy->antiColourLine()->setSourceNeighbours(cline.second->sourceNeighbours().first,
	cline.second->sourceNeighbours().second);
	}
	else if (cline.second->sinkNeighbours().first) {
	copy->antiColourLine()->setSinkNeighbours(cline.second->sinkNeighbours().first,
	cline.second->sinkNeighbours().second);
	}
	}
	// copy of the one from the hard process
	tParticleVector dpar=copy->parents();
	for(unsigned int ix=0;ix<dpar.size();++ix) dpar[ix]->abandonChild(copy);
	final->addChild(copy);
	pstep->addDecayProduct(copy);
	// final-state radiation
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
	for(cit=outgoingLines().begin();cit!=outgoingLines().end();++cit) {
	ShowerParticlePtr init=cit->first->progenitor();
	if(!init->getThePEGBase())
	throw Exception() << "Final-state particle must have a ThePEGBase"
	<< " in ShowerTree::fillEventRecord()"
	<< Exception::runerror;
	// if not from matrix element correction
	if(cit->first->perturbative()) {
	// add the child
	updateColour(cit->first->copy());
	PPtr orig=cit->first->original();
	copy->addChild(orig);
	pstep->addDecayProduct(orig);
	orig->addChild(cit->first->copy());
	pstep->addDecayProduct(cit->first->copy());
	// register the shower particle as a
	// copy of the one from the hard process
	tParticleVector parents=init->parents();
	for(unsigned int ix=0;ix<parents.size();++ix)
	{parents[ix]->abandonChild(init);}
	(*cit).first->copy()->addChild(init);
	pstep->addDecayProduct(init);
	updateColour(init);
	}
	// from a matrix element correction
	else {
	if(copy->children().end()==
	find(copy->children().begin(),copy->children().end(),
	cit->first->original())) {
	updateColour(cit->first->original());
	copy->addChild(cit->first->original());
	pstep->addDecayProduct(cit->first->original());
	}
	updateColour(cit->first->copy());
	cit->first->original()->addChild(cit->first->copy());
	pstep->addDecayProduct(cit->first->copy());
	// register the shower particle as a
	// copy of the one from the hard process
	tParticleVector parents=init->parents();
	for(unsigned int ix=0;ix<parents.size();++ix)
	{parents[ix]->abandonChild(init);}
	(*cit).first->copy()->addChild(init);
	pstep->addDecayProduct(init);
	updateColour(init);
	}
	// insert shower products
	addFinalStateShower(init,pstep);
	}
	_colour.clear();
	}

	void ShowerTree::clear() {
	// reset the has showered flag
	_hasShowered=false;
	// clear the colour map
	_colour.clear();
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
	map<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator cjt;
	// abandon the children of the outgoing particles
	for(cit=_outgoingLines.begin();cit!=_outgoingLines.end();++cit) {
	ShowerParticlePtr orig=cit->first->progenitor();
	orig->set5Momentum(cit->first->copy()->momentum());
	ParticleVector children=orig->children();
	for(unsigned int ix=0;ix<children.size();++ix) orig->abandonChild(children[ix]);
	_outgoingLines[cit->first]=orig;
	cit->first->hasEmitted(false);
	}
	// forward products
	_forward.clear();
	for(cit=_outgoingLines.begin();cit!=_outgoingLines.end();++cit)
	_forward.insert(cit->first->progenitor());
	// if a decay
	if(!_wasHard) {
	ShowerParticlePtr orig=_incomingLines.begin()->first->progenitor();
	orig->set5Momentum(_incomingLines.begin()->first->copy()->momentum());
	ParticleVector children=orig->children();
	for(unsigned int ix=0;ix<children.size();++ix) orig->abandonChild(children[ix]);
	}
	// if a hard process
	else {
	for(cjt=_incomingLines.begin();cjt!=_incomingLines.end();++cjt) {
	tPPtr parent = cjt->first->original()->parents().empty() ?
	tPPtr() : cjt->first->original()->parents()[0];
	ShowerParticlePtr temp=
	new_ptr(ShowerParticle(*cjt->first->copy(),
	cjt->first->progenitor()->perturbative(),
	cjt->first->progenitor()->isFinalState()));
	fixColour(temp);
	temp->x(cjt->first->progenitor()->x());
	cjt->first->hasEmitted(false);
	if(!(cjt->first->progenitor()==cjt->second)&&cjt->second&&parent)
	parent->abandonChild(cjt->second);
	cjt->first->progenitor(temp);
	_incomingLines[cjt->first]=temp;
	}
	}
	// reset the particles at the end of the shower
	_backward.clear();
	// if hard process backward products
	if(_wasHard)
	for(cjt=_incomingLines.begin();cjt!=_incomingLines.end();++cjt)
	_backward.insert(cjt->first->progenitor());
	+ clearTransforms();
	}

	void ShowerTree::resetShowerProducts() {
	map<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator cit;
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
	_backward.clear();
	_forward.clear();
	for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit)
	_backward.insert(cit->second);
	for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt)
	_forward.insert(cjt->second);
	}

	void ShowerTree::updateAfterShower(ShowerDecayMap & decay) {
	// update the links
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mit;
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::iterator tit;
	for(tit=_treelinks.begin();tit!=_treelinks.end();++tit) {
	if(tit->second.first) {
	mit=_outgoingLines.find(tit->second.first);
	if(mit!=_outgoingLines.end()) tit->second.second=mit->second;
	}
	}
	// get the particles coming from those in the hard process
	set<tShowerParticlePtr> hard;
	for(mit=_outgoingLines.begin();mit!=_outgoingLines.end();++mit)
	hard.insert(mit->second);
	// find the shower particles which should be decayed in the
	// shower but didn't come from the hard process
	set<tShowerParticlePtr>::const_iterator cit;
	for(cit=_forward.begin();cit!=_forward.end();++cit) {
	if(decaysInShower((**cit).id())&&
	hard.find(*cit)==hard.end()) {
	ShowerTreePtr newtree=new_ptr(ShowerTree(*cit,decay));
	newtree->setParents();
	newtree->_parent=this;
	Energy width=(cit).dataPtr()->generateWidth((cit).mass());
	decay.insert(make_pair(width,newtree));
	_treelinks.insert(make_pair(newtree,
	make_pair(tShowerProgenitorPtr(),*cit)));
	}
	}
	}

	void ShowerTree::addFinalStateBranching(ShowerParticlePtr parent,
	const ShowerParticleVector & children) {
	assert(children.size()==2);
	_forward.erase(parent);
	for(unsigned int ix=0; ix<children.size(); ++ix) {
	_forward.insert(children[ix]);
	}
	}

	void ShowerTree::addInitialStateBranching(ShowerParticlePtr oldParent,
	ShowerParticlePtr newParent,
	ShowerParticlePtr otherChild) {
	_backward.erase(oldParent);
	_backward.insert(newParent);
	_forward.insert(otherChild);
	}


	void ShowerTree::setParents() {
	// set the parent tree of the children
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
	for(tit=_treelinks.begin();tit!=_treelinks.end();++tit)
	tit->first->_parent=this;
	}

	void ShowerTree::fixColour(tShowerParticlePtr part) {
	ColinePtr line=part->colourLine();
	if(line) {
	line->removeColoured(part);
	line->addColoured(part);
	}
	line=part->antiColourLine();
	if(line) {
	line->removeAntiColoured(part);
	line->addAntiColoured(part);
	}
	}

	vector<ShowerProgenitorPtr> ShowerTree::extractProgenitors() {
	// extract the particles from the ShowerTree
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator mit;
	vector<ShowerProgenitorPtr> ShowerHardJets;
	for(mit=incomingLines().begin();mit!=incomingLines().end();++mit)
	ShowerHardJets.push_back((*mit).first);
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mjt;
	for(mjt=outgoingLines().begin();mjt!=outgoingLines().end();++mjt)
	ShowerHardJets.push_back((*mjt).first);
	return ShowerHardJets;
	}

	-void ShowerTree::transform(const LorentzRotation & boost) {
	+void ShowerTree::transform(const LorentzRotation & boost, bool applyNow) {
	+ if(applyNow) {
	+ // now boost all the particles
	+ map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
	+ // incoming
	+ for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) {
	+ cit->first->progenitor()->deepTransform(boost);
	+ cit->first->copy()->deepTransform(boost);
	+ }
	+ // outgoing
	+ map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
	+ for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) {
	+ cjt->first->progenitor()->deepTransform(boost);
	+ cjt->first->copy()->deepTransform(boost);
	+ }
	+ }
	+ else {
	+ Lorentz5Momentum ptemp1 = _incomingLines.begin()->first->progenitor()->momentum();
	+ Lorentz5Momentum ptemp2 = ptemp1;
	+ ptemp1 *= _transforms;
	+ ptemp1 *= boost;
	+ _transforms.transform(boost);
	+ ptemp2 *= _transforms;
	+ }
	+ // child trees
	+ for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	+ tit=_treelinks.begin();tit!=_treelinks.end();++tit)
	+ tit->first->transform(boost,applyNow);
	+}
	+
	+void ShowerTree::applyTransforms() {
	// now boost all the particles
	map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
	// incoming
	for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) {
	- cit->first->progenitor()->deepTransform(boost);
	- cit->first->copy()->deepTransform(boost);
	+ cit->first->progenitor()->deepTransform(_transforms);
	+ cit->first->copy()->deepTransform(_transforms);
	}
	// outgoing
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
	for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) {
	- cjt->first->progenitor()->deepTransform(boost);
	- cjt->first->copy()->deepTransform(boost);
	+ cjt->first->progenitor()->deepTransform(_transforms);
	+ cjt->first->copy()->deepTransform(_transforms);
	}
	// child trees
	for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	tit=_treelinks.begin();tit!=_treelinks.end();++tit)
	- tit->first->transform(boost);
	+ tit->first->applyTransforms();
	+ _transforms = LorentzRotation();
	}
	+
	+void ShowerTree::clearTransforms() {
	+ _transforms = LorentzRotation();
	+ // child trees
	+ for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
	+ tit=_treelinks.begin();tit!=_treelinks.end();++tit)
	+ tit->first->clearTransforms();
	+}
	diff --git a/Shower/Base/ShowerTree.h b/Shower/Base/ShowerTree.h
	--- a/Shower/Base/ShowerTree.h
	+++ b/Shower/Base/ShowerTree.h
	@@ -1,396 +1,416 @@
	// -- C++ --
	//
	// ShowerTree.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2007 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_ShowerTree_H
	#define HERWIG_ShowerTree_H

	#include "ThePEG/Config/ThePEG.h"
	#include "Herwig++/Shower/ShowerHandler.fh"
	#include "Herwig++/Shower/ShowerConfig.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "ShowerProgenitor.h"
	#include "ThePEG/EventRecord/Step.h"
	#include <cassert>
	#include "ShowerTree.fh"

	namespace Herwig {

	/**
	* Typedef for map of ShowerTrees for decays
	*/
	typedef multimap<Energy,ShowerTreePtr,std::greater<Energy> > ShowerDecayMap;

	using namespace ThePEG;

	/** \ingroup Shower
	*
	* The ShowerTree class stores the basic information needed for
	* each hard interaction, either a scattering process or decay, which
	* needs to be showered.
	*
	*/
	class ShowerTree : public Base {

	/**
	* ShowerHandler is a friend
	*/
	friend class ShowerHandler;

	public:

	/**
	* Constructors and Destructors
	*/
	//@{
	/**
	* Constructor for a scattering process
	* @param incoming The incoming particles
	* @param out The outgoing particles
	* @param decay Map into which the trees for any unstable particles are inserted
	*/
	ShowerTree(const PPair incoming, const ParticleVector & out,
	ShowerDecayMap & decay);

	/**
	* Constructor for a decay
	* @param in The decaying particle
	* @param decay Map into which the trees for any unstable particles are inserted
	*/
	ShowerTree(PPtr in, ShowerDecayMap & decay);
	//@}

	public:

	/**
	* Insert the tree into the event record
	* @param pstep The step into which the particles should be inserted
	* @param ISR Whether or not ISR is switched on
	* @param FSR Whether or not FSR is switched on
	*/
	void fillEventRecord(StepPtr pstep,bool ISR,bool FSR) {
	if(_wasHard)
	insertHard(pstep,ISR,FSR);
	else
	insertDecay(pstep,ISR,FSR);
	}


	/**
	* Set the parent tree to this tree for trees which come from this one.
	* This needs to be run after the constructor.
	*/
	void setParents();

	/**
	* Perform the decay for a tree starting with an unstable particle
	* @param decay The map of widths and ShowerTrees for the decays so that
	* any unstable decay products can be added.
	*/
	void decay(ShowerDecayMap & decay);

	/**
	* Access methods for the type of interaction
	*/
	//@{
	/**
	* Whether or not this is a scattering process
	*/
	bool isHard() const { return _wasHard; }


	/**
	* Whether or not this is a decay.
	*/
	bool isDecay() const { return !_wasHard; }
	//@}

	/**
	* Flags relating to the application of the hard matrix element correction
	*/
	//@{
	/**
	* Was the hard matrix element correction applied
	*/
	bool hardMatrixElementCorrection() const { return _hardMECorrection; }

	/**
	* Set whether or not the hard matrix element correction was applied
	*/
	void hardMatrixElementCorrection(bool in) { _hardMECorrection=in; }
	//@}

	/**
	* Get the incoming shower particles
	*/
	map<ShowerProgenitorPtr,ShowerParticlePtr> & incomingLines() {
	return _incomingLines;
	}

	/**
	* Get the outgoing shower particles
	*/
	map<ShowerProgenitorPtr,tShowerParticlePtr> & outgoingLines() {
	return _outgoingLines;
	}

	/**
	* Update the shower product for a final-state particle
	*/
	void updateFinalStateShowerProduct(ShowerProgenitorPtr progenitor,
	ShowerParticlePtr parent,
	const ShowerParticleVector & children);

	/**
	* Update the shower product for an initial-state particle
	*/
	void updateInitialStateShowerProduct(ShowerProgenitorPtr progenitor,
	ShowerParticlePtr newParent);

	/**
	* Get the current final shower product for a final-state particle
	*/
	tShowerParticlePtr getFinalStateShowerProduct(ShowerProgenitorPtr progenitor) {
	return _outgoingLines.find(progenitor)==_outgoingLines.end()
	? tShowerParticlePtr() : _outgoingLines[progenitor];
	}

	/**
	* Add a final-state branching. This method removes the parent of the branching
	* from the list of particles at the end of the shower and inserts the children
	* @param parent The parent for the branching
	* @param children The outgoing particles in the branching
	*/
	void addFinalStateBranching(ShowerParticlePtr parent,
	const ShowerParticleVector & children);

	/**
	* Add an initial-state branching. This method removes the oldParent of the
	* branching and inserts the result of the backward evolution and the
	* outgoing particle into the relevant lists.
	* @param oldParent The particle being backward evolved
	* @param newParent The initial-state particle resulting from the backward evolution
	* @param otherChild The final-state particle produced in the evolution.
	*/
	void addInitialStateBranching(ShowerParticlePtr oldParent,
	ShowerParticlePtr newParent,
	ShowerParticlePtr otherChild);

	/**
	* Member called at the end of the shower of a tree to perform a number of
	* updates.
	* @param decay The map of widths and ShowerTrees for the decays so that
	* any unstable decay products can be added.
	*/
	void updateAfterShower(ShowerDecayMap & decay);

	/**
	* Access and set the flag for whether this tree has been showered
	*/
	//@{
	/**
	* Access the flag
	*/
	bool hasShowered() const { return _hasShowered; }

	/**
	* Set the flag
	*/
	void hasShowered(bool in) { _hasShowered=in; }
	//@}

	/**
	* Access the parent tree
	*/
	ShowerTreePtr parent() const { return _parent; }

	/**
	* Clear all the shower products so that the event can be reshowered
	* if the first attempt fail
	*/
	void clear();

	/**
	* Reset the particles resulting from the shower to those which started
	* the shower
	*/
	void resetShowerProducts();

	/**
	* Extract the progenitors for the reconstruction
	*/
	vector<ShowerProgenitorPtr> extractProgenitors();

	/**
	* Access to the outgoing particles
	*/
	const set<tShowerParticlePtr> & forwardParticles() const { return _forward; }

	/**
	* Map of particles in this Tree which are the initial particles in other
	* trees
	*/
	const map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> > &
	treelinks() const {return _treelinks;}

	/**
	* Update the link between shower particle and tree
	*/
	void updateLink(tShowerTreePtr tree,
	pair<tShowerProgenitorPtr,tShowerParticlePtr> in) {
	_treelinks[tree] = in;
	}

	/**
	* Transform the tree
	*/
	- void transform(const LorentzRotation & rot);
	+ void transform(const LorentzRotation & rot, bool applyNow);
	+
	+ /**
	+ * Apply any postphoned transformations
	+ */
	+ void applyTransforms();
	+
	+ /**
	+ * Clear any postphoned transformations
	+ */
	+ void clearTransforms();
	+
	+ /**
	+ * Transform which needs to be applied
	+ */
	+ const LorentzRotation & transform() {return _transforms;}

	protected:

	/**
	* Functions to add the shower to the event record.
	*/
	//@{
	/**
	* Insert a hard process
	* @param pstep The step into which the particles should be inserted
	* @param ISR Whether or not ISR is switched on
	* @param FSR Whether or not FSR is switched on
	*/
	void insertHard(StepPtr pstep,bool ISR,bool FSR);

	/**
	* Insert a decay process
	* @param pstep The step into which the particles should be inserted
	* @param ISR Whether or not ISR is switched on
	* @param FSR Whether or not FSR is switched on
	*/
	void insertDecay(StepPtr pstep,bool ISR,bool FSR);

	/**
	* Recursively add the final-state shower from the particle to the event record.
	* @param particle The final-state particle
	* @param step The step
	*/
	void addFinalStateShower(PPtr particle, StepPtr step);

	/**
	* Add the initial-state shwoer from the particle to the step
	* @param particle The final-state particle
	* @param hadron The incoming hadron
	* @param step The step
	* @param addchildren Add the children of the particle
	*/
	void addInitialStateShower(PPtr particle, PPtr hadron,
	StepPtr step, bool addchildren=true);

	/**
	* Update the colour information of a particle prior to insertion into the
	* event record.
	* @param particle The particle for which the colour is updated.
	*/
	void updateColour(PPtr particle);
	//@}

	/**
	* Isolate the colour of the process from the rest of the event.
	* Called in the constructor
	* @param original The original particles
	* @param copy The colour isolated copies
	*/
	void colourIsolate(const vector<PPtr> & original, const vector<PPtr> & copy);

	/**
	* After the creatation of a ShowerParticle make sure it is properly attached
	* to its ColourLine
	* @param part The particle
	*/
	void fixColour(tShowerParticlePtr part);

	private:

	/**
	* Incoming partons for the hard process
	*/
	PPair _incoming;

	/**
	* The incoming ShowerParticles connected to the interaction
	* as the index of a map with the particle the shower backward evolves
	* them to as the value
	*/
	map<ShowerProgenitorPtr,ShowerParticlePtr> _incomingLines;

	/**
	* The outgoing ShowerParticles connected to the interaction
	* as the index of a map with the particle the shower
	* evolves them to as the value
	*/
	map<ShowerProgenitorPtr,tShowerParticlePtr> _outgoingLines;

	/**
	* The outgoing ShowerParticles at the end of the final-state shower
	*/
	set<tShowerParticlePtr> _forward;

	/**
	* The incoming ShowerParticles at the end of the initial-state shower
	*/
	set<tShowerParticlePtr> _backward;

	/**
	* Was the hard matrix element correction applied
	*/
	bool _hardMECorrection;

	/**
	* Was this a scattering process or a decay
	*/
	bool _wasHard;

	/**
	* Map of colour lines used to reset colours when inserted into the event
	*/
	map<ColinePtr,ColinePtr> _colour;

	/**
	* Map of particles in this Tree which are the initial particles in other
	* trees
	*/
	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> > _treelinks;

	/**
	* The parent tree
	*/
	tShowerTreePtr _parent;

	/**
	* Has this tree showered
	*/
	bool _hasShowered;

	+ /**
	+ * The transforms which still need to be applied
	+ */
	+ LorentzRotation _transforms;
	+
	private:

	/**
	* Copy of decay in shower from ShowerHandler
	*/
	static set<long> _decayInShower;

	/**
	* Check if a particle decays in the shower
	* @param id The PDG code for the particle
	*/
	static bool decaysInShower(long id) {
	return ( _decayInShower.find( abs(id) ) !=
	_decayInShower.end() );
	}

	};
	}

	#endif
	diff --git a/Shower/Default/QTildeReconstructor.cc b/Shower/Default/QTildeReconstructor.cc
	--- a/Shower/Default/QTildeReconstructor.cc
	+++ b/Shower/Default/QTildeReconstructor.cc
	@@ -1,2181 +1,2182 @@
	// -- C++ --
	//
	// QTildeReconstructor.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2007 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 <cassert>

	using namespace Herwig;

	namespace {

	enum SystemType { UNDEFINED=-1, II, IF, F ,I };

	struct ColourSingletSystem {

	ColourSingletSystem() : type(UNDEFINED) {}

	ColourSingletSystem(SystemType intype,ShowerProgenitorPtr inpart)
	: type(intype),jets(1,inpart) {}

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

	/**
	* The jets in the system
	*/
	vector<ShowerProgenitorPtr> jets;
	};

	struct ColourSingletShower {

	ColourSingletShower() : type(UNDEFINED) {}

	ColourSingletShower(SystemType intype,HardBranchingPtr inpart)
	: type(intype),jets(1,inpart) {}

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

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

	};

	}

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

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

	ClassDescription<QTildeReconstructor> QTildeReconstructor::initQTildeReconstructor;
	// Definition of the static class description member.

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

	}

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

	bool QTildeReconstructor::
	reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent,
	unsigned int iopt) const {
	assert(particleJetParent);
	bool emitted=true;
	// if this is not a fixed point in the reconstruction
	if( !particleJetParent->isReconstructionFixedPoint() ) {
	// 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),iopt);
	}
	// 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,iopt,
	_progenitor->mass());
	}
	else {
	jetGrandParent->showerKinematics()->reconstructLast(particleJetParent,iopt);
	}
	}
	}
	// 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();
	if(dm>dum.e()) throw KinematicsReconstructionVeto();
	dum.setMass(dm);
	dum.rescaleRho();
	particleJetParent->set5Momentum(dum);
	}
	else {
	emitted=false;
	}
	}
	}
	// recursion has reached an endpoint once, ie we can reconstruct the
	// kinematics from the children.
	if( !(particleJetParent->isReconstructionFixedPoint()) )
	particleJetParent->showerKinematics()
	->reconstructParent( particleJetParent, particleJetParent->children() );
	return emitted;
	}

	bool QTildeReconstructor::
	reconstructHardJets(ShowerTreePtr hard,
	const map<tShowerProgenitorPtr,
	pair<Energy,double> > & intrinsic) const {
	_currentTree = hard;
	_intrinsic=intrinsic;
	// extract the particles from the ShowerTree
	vector<ShowerProgenitorPtr> ShowerHardJets=hard->extractProgenitors();
	try {
	// old recon method, using new member functions
	if(_reconopt==0) {
	reconstructGeneralSystem(ShowerHardJets);
	}
	// reconstruction based on coloured systems
	else {
	// identify the colour singlet systems
	vector<ColourSingletSystem> systems;
	vector<bool> done(ShowerHardJets.size(),false);
	for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
	// if not treated create new system
	if(done[ix]) continue;
	systems.push_back(ColourSingletSystem(UNDEFINED,ShowerHardJets[ix]));
	done[ix] = true;
	if(!ShowerHardJets[ix]->progenitor()->coloured()) continue;
	// now find the colour connected particles
	vector<unsigned int> iloc(1,ix);
	do {
	vector<unsigned int> temp=findPartners(iloc.back(),ShowerHardJets);
	iloc.pop_back();
	for(unsigned int iy=0;iy<temp.size();++iy) {
	if(!done[temp[iy]]) {
	done[temp[iy]] = true;
	iloc.push_back(temp[iy]);
	systems.back().jets.push_back(ShowerHardJets[temp[iy]]);
	}
	}
	}
	while(!iloc.empty());
	}
	// catagorize the systems
	unsigned int nnun(0),nnii(0),nnif(0),nnf(0),nni(0);
	for(unsigned int ix=0;ix<systems.size();++ix) {
	unsigned int ni(0),nf(0);
	for(unsigned int iy=0;iy<systems[ix].jets.size();++iy) {
	if(systems[ix].jets[iy]->progenitor()->isFinalState()) ++nf;
	else ++ni;
	}
	// type
	// initial-initial
	if(ni==2&&nf==0) {
	systems[ix].type = II;
	++nnii;
	}
	// initial only
	else if(ni==1&&nf==0) {
	systems[ix].type = I;
	++nni;
	}
	// initial-final
	else if(ni==1&&nf>0) {
	systems[ix].type = IF;
	++nnif;
	}
	// final only
	else if(ni==0&&nf>0) {
	systems[ix].type = F;
	++nnf;
	}
	// otherwise unknown
	else {
	systems[ix].type = UNDEFINED;
	++nnun;
	}
	}
	// now decide what to do
	// initial-initial connection and final-state colour singlet systems
	LorentzRotation toRest,fromRest;
	bool applyBoost(false);
	bool 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);
	}
	}
	// 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) {
	// only FS needed
	}
	// 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);
	}
	}
	}
	catch(KinematicsReconstructionVeto) {
	_progenitor=tShowerParticlePtr();
	_intrinsic.clear();
	_currentTree = tShowerTreePtr();
	return false;
	}
	_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()) {
	_currentTree = tShowerTreePtr();
	return false;
	}
	}
	_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();
	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,0);
	// 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) const {
	_currentTree = decay;
	try {
	// extract the particles from the ShowerTree
	vector<ShowerProgenitorPtr> ShowerHardJets=decay->extractProgenitors();
	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;
	atLeastOnce \|= reconstructTimeLikeJet(tempJetKin.parent,0);
	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)) {
	_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 {
	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) {
	_currentTree = tShowerTreePtr();
	return false;
	}
	_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,1);
	// 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);
	Lorentz5Momentum ptest;
	for(unsigned int ix=0;ix<jetKinematics.size();++ix) {
	pmag.push_back(jetKinematics[ix].p.vect().mag2());
	total+=jetKinematics[ix].q.mass();
	ptest+=jetKinematics[ix].p;
	}
	// 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;
	root.resize(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);
	}
	// do the initial-state reconstruction
	LorentzRotation toRest,fromRest;
	bool applyBoost(false);
	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);
	}
	// inverse of reconstruction based on coloured systems
	else {
	// identify the colour singlet systems
	vector<ColourSingletShower> systems;
	set<HardBranchingPtr> done;
	for(set<HardBranchingPtr>::const_iterator it=tree->branchings().begin();
	it!=tree->branchings().end();++it) {
	// if not treated create new system
	if(done.find(*it)!=done.end()) continue;
	done.insert(*it);
	systems.push_back(ColourSingletShower(UNDEFINED,*it));
	if(!(**it).branchingParticle()->coloured()) continue;
	// now find the colour connected particles
	findPartners(*it,done,tree->branchings(),systems.back().jets);
	}
	// catagorize the systems
	unsigned int nnun(0),nnii(0),nnif(0),nnf(0),nni(0);
	for(unsigned int ix=0;ix<systems.size();++ix) {
	unsigned int ni(0),nf(0);
	for(unsigned int iy=0;iy<systems[ix].jets.size();++iy) {
	if(systems[ix].jets[iy]->status()==HardBranching::Outgoing) ++nf;
	else ++ni;
	}
	// type
	// initial-initial
	if(ni==2&&nf==0) {
	systems[ix].type = II;
	++nnii;
	}
	// initial only
	else if(ni==1&&nf==0) {
	systems[ix].type = I;
	++nni;
	}
	// initial-final
	else if(ni==1&&nf>0) {
	systems[ix].type = IF;
	++nnif;
	}
	// final only
	else if(ni==0&&nf>0) {
	systems[ix].type = F;
	++nnf;
	}
	// otherwise unknown
	else {
	systems[ix].type = UNDEFINED;
	++nnun;
	}
	}
	// 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);
	}
	}
	// 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();
	}
	// 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);
	}
	}
	else {
	return deconstructGeneralSystem(tree,evolver,type);
	}
	return true;
	}
	}

	vector<unsigned int> QTildeReconstructor::
	findPartners(unsigned int iloc ,
	vector<ShowerProgenitorPtr> jets) const {
	vector<unsigned int> output;
	for(unsigned int iy=0;iy<jets.size();++iy) {
	if(!jets[iy]->progenitor()->data().coloured()\|\|iy==iloc) continue;
	bool isPartner = false;
	// both in either initial or final state
	if(jets[iloc]->progenitor()->isFinalState()!=jets[iy]->progenitor()->isFinalState()) {
	if(jets[iloc]->progenitor()->colourLine() &&
	jets[iloc]->progenitor()->colourLine() == jets[iy]->progenitor()->colourLine())
	isPartner = true;
	if(jets[iloc]->progenitor()->antiColourLine() &&
	jets[iloc]->progenitor()->antiColourLine() == jets[iy]->progenitor()->antiColourLine())
	isPartner = true;
	}
	else {
	if(jets[iloc]->progenitor()->colourLine() &&
	jets[iloc]->progenitor()->colourLine() == jets[iy]->progenitor()->antiColourLine())
	isPartner = true;
	if(jets[iloc]->progenitor()->antiColourLine() &&
	jets[iloc]->progenitor()->antiColourLine() == jets[iy]->progenitor()->colourLine())
	isPartner = true;
	}
	// special for sources/sinks
	if(jets[iloc]->progenitor()->colourLine()) {
	if(jets[iloc]->progenitor()->colourLine()->sourceNeighbours().first) {
	tColinePair lines = jets[iloc]->progenitor()->colourLine()->sourceNeighbours();
	if(lines.first == jets[iy]->progenitor()-> colourLine() \|\|
	lines.first == jets[iy]->progenitor()-> colourLine() \|\|
	lines.second == jets[iy]->progenitor()->antiColourLine() \|\|
	lines.second == jets[iy]->progenitor()->antiColourLine())
	isPartner = true;
	}
	if(jets[iloc]->progenitor()->colourLine()->sinkNeighbours().first) {
	tColinePair lines = jets[iloc]->progenitor()->colourLine()->sinkNeighbours();
	if(lines.first == jets[iy]->progenitor()-> colourLine() \|\|
	lines.first == jets[iy]->progenitor()-> colourLine() \|\|
	lines.second == jets[iy]->progenitor()->antiColourLine() \|\|
	lines.second == jets[iy]->progenitor()->antiColourLine())
	isPartner = true;
	}
	}
	if(jets[iloc]->progenitor()->antiColourLine()) {
	if(jets[iloc]->progenitor()->antiColourLine()->sourceNeighbours().first) {
	tColinePair lines = jets[iloc]->progenitor()->antiColourLine()->sourceNeighbours();
	if(lines.first == jets[iy]->progenitor()-> colourLine() \|\|
	lines.first == jets[iy]->progenitor()-> colourLine() \|\|
	lines.second == jets[iy]->progenitor()->antiColourLine() \|\|
	lines.second == jets[iy]->progenitor()->antiColourLine())
	isPartner = true;
	}
	if(jets[iloc]->progenitor()->antiColourLine()->sinkNeighbours().first) {
	tColinePair lines = jets[iloc]->progenitor()->antiColourLine()->sinkNeighbours();
	if(lines.first == jets[iy]->progenitor()-> colourLine() \|\|
	lines.first == jets[iy]->progenitor()-> colourLine() \|\|
	lines.second == jets[iy]->progenitor()->antiColourLine() \|\|
	lines.second == jets[iy]->progenitor()->antiColourLine())
	isPartner = true;
	}
	}
	if(isPartner) output.push_back(iy);
	}
	return output;
	}

	void QTildeReconstructor::
	reconstructInitialFinalSystem(vector<ShowerProgenitorPtr> jets) const {
	Lorentz5Momentum pin[2],pout[2];
	bool atLeastOnce(false);
	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();
	atLeastOnce \|= reconstructTimeLikeJet(jets[ix]->progenitor(),0);
	}
	// initial-state parton
	else {
	pin[0] +=jets[ix]->progenitor()->momentum();
	atLeastOnce \|= reconstructSpaceLikeJet(jets[ix]->progenitor());
	assert(!jets[ix]->original()->parents().empty());
	}
	}
	// add intrinsic pt if needed
	atLeastOnce \|= 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];
	Lorentz5Momentum pb = pin[0];
	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());
	Lorentz5Momentum ptemp=rot*pb;
	Boost trans = -1./ptemp.e()*ptemp.vect();
	trans.setZ(0.);
	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];
	double a[2],b[2];
	a[0] = n2*qbp/n1n2;
	b[0] = n1*qbp/n1n2;
	Lorentz5Momentum qperp = qbp-a[0]n1-b[0]n2;
	a[1] = 0.5*(qcp.m2()-qperp.m2())/n1n2;
	b[1] = 1.;
	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();
	double kb = 0.5(-B+sqrt(sqr(B)-4.A*C))/A;
	double kc = (a[0]*kb-0.5)/a[1];
	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();
	LorentzRotation transb=rotinvsolveBoostZ(pnew[0],qbp)rot;
	LorentzRotation transc=rotinvsolveBoost(pnew[1],qcp)rot;
	for(unsigned int ix=0;ix<jets.size();++ix) {
	if(jets[ix]->progenitor()->isFinalState())
	deepTransform(jets[ix]->progenitor(),transc);
	else {
	tPPtr parent;
	boostChain(jets[ix]->progenitor(),transb,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()) 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);
	Boost beta = -betam(k/kp)oldp.vect();
	// 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;
	if ( ax.mag2()/GeV2/MeV2 > 1e-16 ) {
	R.rotate( delta, unitVector(ax) ).boost( beta );
	}
	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 );
	}
	else {
	R.boost( beta );
	}
	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());
	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 ( 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 Energy2 eps2 = 1e-8*GeV2;
	static const Energy eps = 1e-4 *GeV;
	LorentzRotation R;
	double beta;
	Energy2 den = (p.t()q.t()-p.z()q.z());
	Energy2 num = -(p.z()q.t()-q.z()p.t());
	if(abs(den)<eps2\|\|abs(num)<eps2) {
	if(abs(p.t()-abs(p.z()))<eps&&abs(q.t()-abs(q.z()))<eps) {
	double ratio = sqr(q.t()/p.t());
	beta = -(1.-ratio)/(1.+ratio);
	}
	else {
	beta=0.;
	}
	}
	else {
	beta = num/den;

	}
	R.boostZ(beta);
	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 {
	// special for case of individual particle
	if(jets.size()==1) {
	LorentzRotation trans(toRest);
	trans.transform(fromRest);
	deepTransform(jets[0]->progenitor(),trans);
	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();
	}
	// check if in CMF frame
	Boost beta_cm = pcm.findBoostToCM();
	bool gottaBoost = (beta_cm.mag() > 1e-12);
	// 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(gottaBoost) {
	tempJetKin.parent->boost(beta_cm);
	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==tempJetKin.parent)
	- tit->first->transform(LorentzRotation(beta_cm));
	+ tit->first->transform(LorentzRotation(beta_cm),false);
	}
	}
	tempJetKin.p = (*cit)->progenitor()->momentum();
	_progenitor=tempJetKin.parent;
	radiated \|= reconstructTimeLikeJet((*cit)->progenitor(),0);
	tempJetKin.q = (*cit)->progenitor()->momentum();
	jetKinematics.push_back(tempJetKin);
	}
	// 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 = LorentzRotation();
	if(radiated) Trafo = solveBoost(k, it->q, it->p);
	if(gottaBoost) Trafo.boost(-beta_cm);
	if(applyBoost) {
	Trafo.transform( toRest);
	Trafo.transform(fromRest);
	}
	if(radiated \|\| gottaBoost \|\| applyBoost) deepTransform(it->parent,Trafo);
	}
	}

	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();
	}
	// check if intrinsic pt to be added
	radiated \|= !_intrinsic.empty();
	// if no radiation return
	if(!radiated) return;
	// initial state shuffling
	applyBoost=false;
	vector<Lorentz5Momentum> p, pq, p_in;
	for(unsigned int ix=0;ix<jets.size();++ix) {
	// at momentum to vector
	p_in.push_back(jets[ix]->progenitor()->momentum());
	// reconstruct the jet
	radiated \|= reconstructSpaceLikeJet(jets[ix]->progenitor());
	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);
	}
	// 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();
	Energy MDY = (p_in[0] + p_in[1]).m();
	Energy2 S = (pq[0]+pq[1]).m2();
	// if not need don't apply boosts
	if(!(radiated && p.size() == 2 && pq.size() == 2)) return;
	applyBoost=true;
	// 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
	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, k2
	rad = kp(b[0]+kpb[1])/(kpa[0]+a[1])(x1/x2);
	if(rad <= 0.) throw KinematicsReconstructionVeto();
	double k1 = sqrt(rad);
	double k2 = kp/k1;
	double beta[2] =
	{getBeta((a[0]+b[0]), (a[0]-b[0]), (k1a[0]+b[0]/k1), (k1a[0]-b[0]/k1)),
	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];
	}
	// 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(betaboost),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();
	// do one boost
	toRest = LorentzRotation(pcm.findBoostToCM());
	if(_initialBoost==0) {
	fromRest = LorentzRotation(newcmf.boostVector());
	}
	else if(_initialBoost==1) {
	// first apply longitudinal boost
	double beta = newcmf.z()/sqrt(newcmf.m2()+sqr(newcmf.z()));
	fromRest=LorentzRotation(Boost(0.,0.,beta));
	// then transverse one
	Energy pT = sqrt(sqr(newcmf.x())+sqr(newcmf.y()));
	beta = pT/newcmf.t();
	fromRest.boost(Boost(betanewcmf.x()/pT,betanewcmf.y()/pT,0.));
	}
	else
	assert(false);
	}

	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);
	// decompose the momenta
	double alpha[2],beta[2];
	Energy2 p12=pq[0]*pq[1];
	Lorentz5Momentum pt[2];
	for(unsigned int ix=0;ix<2;++ix) {
	alpha[ix] = pin[ix]*pq[1]/p12;
	beta [ix] = pin[ix]*pq[0]/p12;
	pt[ix] = pin[ix]-alpha[ix]pq[0]-beta[ix]pq[1];
	}
	// parton level centre-of-mass
	Lorentz5Momentum pcm=pin[0]+pin[1];
	pcm.rescaleMass();
	double rap=pcm.rapidity();
	// hadron level cmf
	Energy2 s = (pq[0] +pq[1] ).m2();
	// calculate the x values
	double x[2]={sqrt(pcm.mass2()/sexp(2.rap)),pcm.mass2()/s/x[0]};
	if(pq[0].z()<ZERO) swap(x[0],x[1]);
	double k1=alpha[0]/x[0],k2=beta[1]/x[1];
	double alphanew[2]={alpha[0]/k1,alpha[1]*k2};
	double betanew [2]={beta [0]*k1,beta [1]/k2};
	double boost[2];
	for(unsigned int ix=0;ix<2;++ix) {
	boost[ix] = getBeta(alpha [ix]+beta [ix], alpha[ix] -beta [ix],
	alphanew[ix]+betanew[ix], alphanew[ix]-betanew[ix]);
	if (pq[0].z() > ZERO) beta[ix]*=-1.;
	}
	// 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 {
	if(jets.size()==1) {
	LorentzRotation R(toRest);
	R.transform(fromRest);
	// \todo What does this do? tree->showerRot( R );
	jets[0]->original(R*jets[0]->branchingParticle()->momentum());
	jets[0]->showerMomentum(R*jets[0]->branchingParticle()->momentum());
	// 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;
	for(cit=jets.begin();cit!=jets.end();++cit) {
	pout.push_back((*cit)->branchingParticle()->momentum());
	// 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
	if((*cit)->branchingParticle()->children().size()==0 \|\|
	(!(*cit)->branchingParticle()->dataPtr()->coloured() &&
	!(*cit)->branchingParticle()->dataPtr()->stable()) )
	mon.push_back(pout.back().mass());
	else
	mon.push_back((*cit)->branchingParticle()->dataPtr()->mass());
	}
	// boost all the momenta to the rest frame of the decaying particle
	Lorentz5Momentum pin;
	for(unsigned int ix=0;ix<pout.size();++ix) {
	pout[ix].transform(toRest);
	pin += pout[ix];
	}
	pin.rescaleMass();
	// rescaling factor
	double lambda=inverseRescalingFactor(pout,mon,pin.mass());
	// now calculate the p reference vectors
	for(cit=jets.begin();cit!=jets.end();++cit){
	Lorentz5Momentum pvect = (*cit)->branchingParticle()->momentum();
	pvect.transform(toRest);
	pvect /= lambda;
	if((*cit)->branchingParticle()->children().size()==0 \|\|
	(!(*cit)->branchingParticle()->dataPtr()->coloured() &&
	!(*cit)->branchingParticle()->dataPtr()->stable()) )
	pvect.setMass((*cit)->branchingParticle()->momentum().mass());
	else
	pvect.setMass((*cit)->branchingParticle()->dataPtr()->mass());
	pvect.rescaleEnergy();
	pvect.transform(fromRest);
	(*cit)->pVector(pvect);
	(*cit)->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);
	}
	}

	void QTildeReconstructor::
	reconstructGeneralSystem(vector<ShowerProgenitorPtr> & ShowerHardJets) const {
	// general recon, all initial-state in one system and final-state
	// in another
	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);
	reconstructInitialInitialSystem(applyBoost,toRest,fromRest,in.jets);
	// reconstruct the final-state systems
	reconstructFinalStateSystem(applyBoost,toRest,fromRest,out.jets);
	}

	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::
	findPartners(HardBranchingPtr branch,set<HardBranchingPtr> & done,
	const set<HardBranchingPtr> & branchings,
	vector<HardBranchingPtr> & jets) const {
	tShowerParticlePtr part=branch->branchingParticle();
	for(set<HardBranchingPtr>::const_iterator cit=branchings.begin();
	cit!=branchings.end();++cit) {
	if(done.find(cit)!=done.end()\|\|!(*cit).branchingParticle()->coloured())
	continue;
	bool isPartner = false;
	// one initial and one final
	if(branch->status()!=(**cit).status()) {
	if(part->colourLine() &&
	part->colourLine() == (**cit).branchingParticle()->colourLine())
	isPartner = true;
	if(part->antiColourLine() &&
	part->antiColourLine() == (**cit).branchingParticle()->antiColourLine())
	isPartner = true;
	}
	// both in either initial or final state
	else {
	if(part->colourLine() &&
	part->colourLine() == (**cit).branchingParticle()->antiColourLine())
	isPartner = true;
	if(part->antiColourLine() &&
	part->antiColourLine() == (**cit).branchingParticle()->colourLine())
	isPartner = true;
	}
	if(isPartner) {
	jets.push_back(*cit);
	done.insert(*cit);
	findPartners(*cit,done,branchings,jets);
	}
	}
	}

	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()->getThePEGBase() ?
	jets[ix]->branchingParticle()->getThePEGBase()->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 {
	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);
	if ( particle->next() ) deepTransform(particle->next(),r,match,original);
	if(!match) return;
	if(!particle->children().empty()) return;
	// 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);
	+ tit->first->transform(r*rot,false);
	}
	}
	}

File Metadata

Mime Type: text/x-diff
Expires: Thu, Apr 24, 6:37 AM (1 d, 19 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4864948
Default Alt Text: (179 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions