No OneTemporary
Actions

Size

115 KB

Subscribers

None

View Options

	diff --git a/Hadronization/ClusterFinder.cc b/Hadronization/ClusterFinder.cc
	--- a/Hadronization/ClusterFinder.cc
	+++ b/Hadronization/ClusterFinder.cc
	@@ -1,366 +1,359 @@
	// -- C++ --
	//
	// ClusterFinder.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ClusterFinder class.
	//

	#include "ClusterFinder.h"
	#include <ThePEG/Interface/ClassDocumentation.h>
	#include <ThePEG/PDT/StandardMatchers.h>
	#include <ThePEG/PDT/EnumParticles.h>
	#include <ThePEG/Repository/EventGenerator.h>
	#include <ThePEG/EventRecord/Collision.h>
	#include "CheckId.h"
	#include "Herwig++/Utilities/EnumParticles.h"
	#include "Cluster.h"
	#include <ThePEG/Utilities/DescribeClass.h>

	using namespace Herwig;

	DescribeNoPIOClass<ClusterFinder,Interfaced>
	describeClusterFinder("Herwig::ClusterFinder","");

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

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

	void ClusterFinder::Init() {

	static ClassDocumentation<ClusterFinder> documentation
	("This class is responsible of finding clusters.");

	}


	-ClusterVector ClusterFinder::formClusters(const PVector & partons)
	- {
	+ClusterVector ClusterFinder::formClusters(const PVector & partons) {

	set<tPPtr> examinedSet; // colour particles already included in a cluster
	map<tColinePtr, pair<tPPtr,tPPtr> > quarkQuark; // quark quark
	map<tColinePtr, pair<tPPtr,tPPtr> > aQuarkQuark; // anti quark anti quark
	ParticleSet inputParticles(partons.begin(),partons.end());

	ClusterVector clusters;

	// Loop over all current particles.
	for(PVector::const_iterator pit=partons.begin();pit!=partons.end();++pit){
	// Skip to the next particle if it is not coloured or already examined.
	assert(*pit);
	assert((*pit)->dataPtr());
	if(!(**pit).data().coloured()
	\|\| examinedSet.find(*pit) != examinedSet.end()) {
	continue;
	}
	// We assume that a cluster is made of, at most, 3 constituents although
	// in most cases the number will be 2; however, for baryon violating decays
	// (for example in Susy model without R parity conservation) we can have 3
	// constituents. In the latter case, a quark (antiquark) do not have an
	// anticolour (colour) partner as usual, but its colour line either stems
	// from a colour source, or ends in a colour sink. In the case of double
	// baryon violating decays, but with overall baryon conservation
	// ( for instance:
	// tilde_u_R -> dbar_1 + dbar_2
	// tilde_u_R_star -> d1 + d2
	// where tilde_u_R and tilde_u_R_star are colour connected )
	// a special treatment is needed, because first we have to process all
	// partons in the current step, and then for each left pair of quarks which
	// stem from a colour source we have to find the corresponding pair of
	// anti-quarks which ends in a colour sink and is connected with the
	// above colour source. These special pairs are kept into the maps:
	// spec/CluHadConfig.hialQuarkQuarkMap and specialAntiQuarkAntiQuarkMap.

	tParticleVector connected(3);
	int iElement = 0;
	connected[iElement++] = *pit;
	bool specialCase = false;

	if((*pit)->hasColour()) {
	tPPtr partner =
	(*pit)->colourLine()->getColouredParticle(partons.begin(),
	partons.end(),
	true);

	if(partner) {
	connected[iElement++]= partner;
	}
	// colour source : baryon-violating process
	else {
	if((*pit)->colourLine()->sourceNeighbours() != tColinePair()) {
	tColinePair sourcePair = (*pit)->colourLine()->sourceNeighbours();
	tColinePtr intCL = tColinePtr();
	for(int i = 0; i < 2; ++i) {
	tColinePtr pLine = i==0 ? sourcePair.first : sourcePair.second;
	int saveNumElements = iElement;
	for(tPVector::const_iterator cit = pLine->coloured().begin();
	cit != pLine->coloured().end(); ++cit ) {
	ParticleSet::const_iterator cjt = inputParticles.find(*cit);
	if(cjt!=inputParticles.end()) connected[iElement++]= (*cit);
	}
	if(iElement == saveNumElements) intCL = pLine;
	}
	if(intCL && iElement == 2) {
	specialCase = true;
	pair<tPPtr,tPPtr> qp=pair<tPPtr,tPPtr>(connected[0],connected[1]);
	quarkQuark.insert(pair<tColinePtr,pair<tPPtr,tPPtr> >(intCL,qp));
	}
	else if(iElement != 3) {
	throw Exception() << "Colour connections fail in the hadronization for "
	<< **pit << "in ClusterFinder::formClusters"
	<< " for a coloured particle."
	<< " Failed to find particles from a source"
	<< Exception::runerror;
	}
	}
	else {
	throw Exception() << "Colour connections fail in the hadronization for "
	<< **pit << "in ClusterFinder::formClusters for"
	<< " a coloured particle"
	<< Exception::runerror;
	}
	}
	}

	if((*pit)->hasAntiColour()) {
	tPPtr partner =
	(*pit)->antiColourLine()->getColouredParticle(partons.begin(),
	partons.end(),
	false);
	if(partner) {
	connected[iElement++]=partner;
	}
	// colour sink : baryon-violating process
	else {
	if((*pit)->antiColourLine()->sinkNeighbours() != tColinePair()) {
	tColinePair sinkPair = (*pit)->antiColourLine()->sinkNeighbours();
	tColinePtr intCL = tColinePtr();
	for(int i = 0; i < 2; ++i) {
	tColinePtr pLine = i==0 ? sinkPair.first : sinkPair.second;
	int saveNumElements = iElement;
	for(tPVector::const_iterator cit = pLine->antiColoured().begin();
	cit != pLine->antiColoured().end(); ++cit ) {
	ParticleSet::const_iterator cjt = inputParticles.find(*cit);
	if(cjt!=inputParticles.end()) connected[iElement++]= (*cit);
	}
	if(iElement == saveNumElements) intCL = pLine;
	}
	if(intCL && iElement == 2) {
	specialCase = true;
	pair<tPPtr,tPPtr> aqp=pair<tPPtr,tPPtr>(connected[0],connected[1]);
	aQuarkQuark.insert(pair<tColinePtr,pair<tPPtr,tPPtr> >(intCL,aqp));
	}
	else if( iElement !=3) {
	throw Exception() << "Colour connections fail in the hadronization for "
	<< **pit << "in ClusterFinder::formClusters for"
	<< " an anti-coloured particle."
	<< " Failed to find particles from a sink"
	<< Exception::runerror;
	}
	}
	else {
	throw Exception() << "Colour connections fail in the hadronization for "
	<< **pit << "in ClusterFinder::formClusters for"
	<< " an anti-coloured particle"
	<< Exception::runerror;
	}
	}
	}

	if(!specialCase) {
	// Tag the components of the found cluster as already examined.
	for (int i=0; i<iElement; ++i) examinedSet.insert(connected[i]);
	// Create the cluster object with the colour connected particles
	ClusterPtr cluPtr = new_ptr(Cluster(connected[0],connected[1],
	connected[2]));
	// add to the step
	connected[0]->addChild(cluPtr);
	connected[1]->addChild(cluPtr);
	if(connected[2]) connected[2]->addChild(cluPtr);
	clusters.push_back(cluPtr);
	// Check if any of the components is a beam remnant, and if this
	// is the case then inform the cluster.
	// this will only work for baryon collisions
	for (int i=0; i<iElement; ++i) {
	bool fromRemnant = false;
	tPPtr parent=connected[i];
	while(parent) {
	if(parent->id()==ParticleID::Remnant) {
	fromRemnant = true;
	break;
	}
	parent = parent->parents().empty() ? tPPtr() : parent->parents()[0];
	}
	if(fromRemnant&&DiquarkMatcher::Check(connected[i]->id()))
	cluPtr->isBeamCluster(connected[i]);
	}
	}

	}

	// Treat now the special cases, if any. The idea is to find for each pair
	// of quarks coming from a common colour source the corresponding pair of
	// antiquarks coming from a common colour sink, connected to the above
	// colour source via the same colour line. Then, randomly couple one of
	// the two quarks with one of the two antiquarks, and do the same with the
	// quark and antiquark left.
	for(map<tColinePtr, pair<tPPtr,tPPtr> >::const_iterator
	cit = quarkQuark.begin(); cit != quarkQuark.end(); ++cit ) {
	tColinePtr coline = cit->first;
	pair<tPPtr,tPPtr> quarkPair = cit->second;
	if(aQuarkQuark.find( coline ) != aQuarkQuark.end()) {
	pair<tPPtr,tPPtr> antiQuarkPair = aQuarkQuark.find(coline)->second;
	ClusterPtr cluPtr1, cluPtr2;
	if ( UseRandom::rndbool() ) {
	cluPtr1 = new_ptr(Cluster(quarkPair.first , antiQuarkPair.first));
	cluPtr2 = new_ptr(Cluster(quarkPair.second , antiQuarkPair.second));
	quarkPair.first->addChild(cluPtr1);
	antiQuarkPair.first->addChild(cluPtr1);
	quarkPair.second->addChild(cluPtr2);
	antiQuarkPair.second->addChild(cluPtr2);
	} else {
	cluPtr1 = new_ptr(Cluster(quarkPair.first , antiQuarkPair.second));
	cluPtr2 = new_ptr(Cluster(quarkPair.second , antiQuarkPair.first));
	quarkPair.second->addChild(cluPtr2);
	antiQuarkPair.first->addChild(cluPtr2);
	quarkPair.first->addChild(cluPtr1);
	antiQuarkPair.second->addChild(cluPtr1);
	}
	clusters.push_back(cluPtr1);
	clusters.push_back(cluPtr2);
	}
	else {
	throw Exception() << "ClusterFinder::formClusters : "
	<< "***Skip event: unable to match pairs in "
	<< "Baryon-violating processes***"
	<< Exception::eventerror;
	}
	}
	return clusters;
	}


	void ClusterFinder::reduceToTwoComponents(ClusterVector & clusters) {

	// In order to preserve all of the information, we do not modify the
	// directly the 3-component clusters, but instead we define new clusters,
	// which are related to the original ones by a child-parent relationship,
	// by considering two randomly chosen components as a diquark (or anti-diquark).
	// These new clusters are first added to the vector vecNewRedefinedCluPtr,
	// and at the end, when all input clusters have been examined, the elements of
	// this vector will be copied in collecCluPtr (the reason is that it is not
	// allowed to modify a STL container while iterating over it).
	vector<tClusterPtr> redefinedClusters;
	for(ClusterVector::iterator cluIter = clusters.begin() ;
	cluIter != clusters.end() ; ++cluIter) {
	tParticleVector vec;

	if ( ! (*cluIter)->isAvailable()
	\|\| (*cluIter)->numComponents() != 3 ) continue;

	tPPtr other;
	int iCharge1(0);
	for(int i = 0; i<(*cluIter)->numComponents(); i++) {
	tPPtr part = (*cluIter)->particle(i);
	if(!DiquarkMatcher::Check(*(part->dataPtr())))
	vec.push_back(part);
	else
	other = part;
	iCharge1 += part->dataPtr()->iCharge();
	}
	if(vec.size()<2) {
	throw Exception() << "Could not make a diquark for a baryonic cluster decay from "
	<< (*cluIter)->particle(0)->PDGName() << " "
	<< (*cluIter)->particle(1)->PDGName() << " "
	<< (*cluIter)->particle(2)->PDGName() << " "
	<< " in ClusterFinder::reduceToTwoComponents()."
	<< Exception::eventerror;
	}

	// Randomly selects two components to be considered as a (anti)diquark
	// and place them as the first and second element of vec.
	int choice = vec.size()==2 ? 0 : UseRandom::rnd3(1.0, 1.0, 1.0);
	switch (choice) {
	case 0:
	break;
	case 1:
	swap(vec[2],vec[0]);
	break;
	case 2:
	swap(vec[2],vec[1]);
	break;
	}
	tcPDPtr temp1 = vec[0]->dataPtr();
	tcPDPtr temp2 = vec[1]->dataPtr();
	if(!other) other = vec[2];

	tcPDPtr dataDiquark = CheckId::makeDiquark(temp1,temp2);

	if(!dataDiquark)
	throw Exception() << "Could not make a diquark from"
	<< temp1->PDGName() << " and "
	<< temp2->PDGName()
	<< " in ClusterFinder::reduceToTwoComponents()"
	<< Exception::eventerror;


	// Create the new cluster (with two components) and assign to it the same
	// momentum and position of the original (with three components) one.
	// Furthermore, assign to the diquark component a momentum given by the
	// sum of the two original components from which has been formed; for the
	// position, we are assuming, very simply, that the diquark position is
	// the average positions of the two original components.
	// Notice that the mass (5-th component of the 5-momentum) of the diquark
	// is set by hand to the constituent mass of the diquark (which is equal
	// to the sum of the constituent masses of the two quarks which form the
	// diquark) because the sum of 5-component vectors do add only the "normal"
	// 4-components, not the 5-th one. After that, the 5-momentum of the diquark
	// is in an inconsistent state, because the mass (5-th component) is not
	// equal to the invariant mass obtained from the 4-momemtum. This is not
	// unique to this kind of component (all perturbative components are in
	// a similar situation), but it is not harmful.

	+ // construct the diquark
	PPtr diquark = dataDiquark->produceParticle();
	vec[0]->addChild(diquark);
	vec[1]->addChild(diquark);
	+ diquark->set5Momentum(Lorentz5Momentum(vec[0]->momentum() + vec[1]->momentum(),
	+ dataDiquark->constituentMass()));
	+ diquark->setVertex(0.5*(vec[0]->vertex() + vec[1]->vertex()));
	+ // make the new cluster
	ClusterPtr nclus = new_ptr(Cluster(other,diquark));
	-
	//vec[0]->addChild(nclus);
	//diquark->addChild(nclus);
	- (*cluIter)->addChild(nclus);

	- nclus->set5Momentum((*cluIter)->momentum());
	- nclus->setVertex((*cluIter)->vertex());
	- for(int i = 0; i<nclus->numComponents(); i++) {
	- if(nclus->particle(i)->id() == dataDiquark->id()) {
	- nclus->particle(i)->set5Momentum(Lorentz5Momentum(vec[0]->momentum()
	- + vec[1]->momentum(), dataDiquark->constituentMass()));
	- nclus->particle(i)->setVertex(0.5*(vec[0]->vertex()
	- + vec[1]->vertex()));
	- }
	- }
	// Set the parent/children relationship between the original cluster
	// (the one with three components) with the new one (the one with two components)
	// and add the latter to the vector of new redefined clusters.
	- //(*cluIter)->addChild(nclus);
	+ (*cluIter)->addChild(nclus);
	+
	redefinedClusters.push_back(nclus);
	}

	// Add to collecCluPtr all of the redefined new clusters (indeed the
	// pointers to them are added) contained in vecNewRedefinedCluPtr.
	/// \todo why do we keep the original of the redefined clusters?
	for (tClusterVector::const_iterator it = redefinedClusters.begin();
	it != redefinedClusters.end(); ++it) {
	clusters.push_back(*it);
	}

	}
	diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc
	--- a/Hadronization/ClusterHadronizationHandler.cc
	+++ b/Hadronization/ClusterHadronizationHandler.cc
	@@ -1,328 +1,291 @@
	// -- C++ --
	//
	// ClusterHadronizationHandler.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ClusterHadronizationHandler class.
	//

	#include "ClusterHadronizationHandler.h"
	#include <ThePEG/Interface/ClassDocumentation.h>
	#include <ThePEG/Persistency/PersistentOStream.h>
	#include <ThePEG/Persistency/PersistentIStream.h>
	#include <ThePEG/Interface/Parameter.h>
	#include <ThePEG/Interface/Reference.h>
	#include <ThePEG/Handlers/EventHandler.h>
	#include <ThePEG/Handlers/Hint.h>
	#include <ThePEG/PDT/ParticleData.h>
	#include <ThePEG/EventRecord/Particle.h>
	#include <ThePEG/EventRecord/Step.h>
	#include <ThePEG/PDT/PDT.h>
	#include <ThePEG/PDT/EnumParticles.h>
	#include <ThePEG/Utilities/Throw.h>
	#include "Herwig++/Utilities/EnumParticles.h"
	#include "CluHadConfig.h"
	#include "Cluster.h"
	#include <ThePEG/Utilities/DescribeClass.h>

	using namespace Herwig;

	ClusterHadronizationHandler * ClusterHadronizationHandler::currentHandler_ = 0;

	DescribeClass<ClusterHadronizationHandler,HadronizationHandler>
	describeClusterHadronizationHandler("Herwig::ClusterHadronizationHandler","");

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

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

	void ClusterHadronizationHandler::persistentOutput(PersistentOStream & os)
	const {
	os << _partonSplitter
	<< _clusterFinder
	<< _colourReconnector
	<< _clusterFissioner
	<< _lightClusterDecayer
	<< _clusterDecayer
	<< ounit(_minVirtuality2,GeV2)
	<< ounit(_maxDisplacement,mm)
	<< _underlyingEventHandler;
	}


	void ClusterHadronizationHandler::persistentInput(PersistentIStream & is, int) {
	is >> _partonSplitter
	>> _clusterFinder
	>> _colourReconnector
	>> _clusterFissioner
	>> _lightClusterDecayer
	>> _clusterDecayer
	>> iunit(_minVirtuality2,GeV2)
	>> iunit(_maxDisplacement,mm)
	>> _underlyingEventHandler;
	}


	void ClusterHadronizationHandler::Init() {

	static ClassDocumentation<ClusterHadronizationHandler> documentation
	("This is the main handler class for the Cluster Hadronization",
	"The hadronization was performed using the cluster model of \\cite{Webber:1983if}.",
	"%\\cite{Webber:1983if}\n"
	"\\bibitem{Webber:1983if}\n"
	" B.~R.~Webber,\n"
	" ``A QCD Model For Jet Fragmentation Including Soft Gluon Interference,''\n"
	" Nucl.\\ Phys.\\ B {\\bf 238}, 492 (1984).\n"
	" %%CITATION = NUPHA,B238,492;%%\n"
	// main manual
	);

	static Reference<ClusterHadronizationHandler,PartonSplitter>
	interfacePartonSplitter("PartonSplitter",
	"A reference to the PartonSplitter object",
	&Herwig::ClusterHadronizationHandler::_partonSplitter,
	false, false, true, false);

	static Reference<ClusterHadronizationHandler,ClusterFinder>
	interfaceClusterFinder("ClusterFinder",
	"A reference to the ClusterFinder object",
	&Herwig::ClusterHadronizationHandler::_clusterFinder,
	false, false, true, false);

	static Reference<ClusterHadronizationHandler,ColourReconnector>
	interfaceColourReconnector("ColourReconnector",
	"A reference to the ColourReconnector object",
	&Herwig::ClusterHadronizationHandler::_colourReconnector,
	false, false, true, false);

	static Reference<ClusterHadronizationHandler,ClusterFissioner>
	interfaceClusterFissioner("ClusterFissioner",
	"A reference to the ClusterFissioner object",
	&Herwig::ClusterHadronizationHandler::_clusterFissioner,
	false, false, true, false);

	static Reference<ClusterHadronizationHandler,LightClusterDecayer>
	interfaceLightClusterDecayer("LightClusterDecayer",
	"A reference to the LightClusterDecayer object",
	&Herwig::ClusterHadronizationHandler::_lightClusterDecayer,
	false, false, true, false);

	static Reference<ClusterHadronizationHandler,ClusterDecayer>
	interfaceClusterDecayer("ClusterDecayer",
	"A reference to the ClusterDecayer object",
	&Herwig::ClusterHadronizationHandler::_clusterDecayer,
	false, false, true, false);

	static Parameter<ClusterHadronizationHandler,Energy2> interfaceMinVirtuality2
	("MinVirtuality2",
	"Minimum virtuality^2 of partons to use in calculating distances (unit [GeV2]).",
	&ClusterHadronizationHandler::_minVirtuality2, GeV2, 0.1GeV2, ZERO, 10.0GeV2,false,false,false);

	static Parameter<ClusterHadronizationHandler,Length> interfaceMaxDisplacement
	("MaxDisplacement",
	"Maximum displacement that is allowed for a particle (unit [millimeter]).",
	&ClusterHadronizationHandler::_maxDisplacement, mm, 1.0e-10*mm,
	0.0mm, 1.0e-9mm,false,false,false);

	static Reference<ClusterHadronizationHandler,StepHandler> interfaceUnderlyingEventHandler
	("UnderlyingEventHandler",
	"Pointer to the handler for the Underlying Event. "
	"Set to NULL to disable.",
	&ClusterHadronizationHandler::_underlyingEventHandler, false, false, true, true, false);
	}

	namespace {
	void extractChildren(tPPtr p, set<PPtr> & all) {
	if (p->children().empty()) return;

	for (PVector::const_iterator child = p->children().begin();
	child != p->children().end(); ++child) {
	all.insert(*child);
	extractChildren(*child, all);
	}
	}
	}

	void ClusterHadronizationHandler::
	handle(EventHandler & ch, const tPVector & tagged,
	const Hint &) {
	useMe();
	currentHandler_ = this;
	PVector currentlist(tagged.begin(),tagged.end());
	// set the scale for coloured particles to just above the gluon mass squared
	// if less than this so they are classed as perturbative
	Energy2 Q02 = 1.01*sqr(getParticleData(ParticleID::g)->constituentMass());
	for(unsigned int ix=0;ix<currentlist.size();++ix) {
	if(currentlist[ix]->scale()<Q02) currentlist[ix]->scale(Q02);
	}
	+
	// split the gluons
	_partonSplitter->split(currentlist);
	// form the clusters
	ClusterVector clusters =
	_clusterFinder->formClusters(currentlist);

	_clusterFinder->reduceToTwoComponents(clusters);

	// perform colour reconnection if needed and then
	// decay the clusters into one hadron
	bool lightOK = false;
	short tried = 0;
	const ClusterVector savedclusters = clusters;
	tPVector finalHadrons; // only needed for partonic decayer
	while (!lightOK && tried++ < 10) {

	// no colour reconnection with baryon-number-violating (BV) clusters
	ClusterVector CRclusters, BVclusters;
	CRclusters.reserve( clusters.size() );
	BVclusters.reserve( clusters.size() );
	for (size_t ic = 0; ic < clusters.size(); ++ic) {
	ClusterPtr cl = clusters.at(ic);
	bool hasClusterParent = false;
	for (unsigned int ix=0; ix < cl->parents().size(); ++ix) {
	if (cl->parents()[ix]->id() == ParticleID::Cluster) {
	hasClusterParent = true;
	break;
	}
	}
	if (cl->numComponents() > 2 \|\| hasClusterParent) BVclusters.push_back(cl);
	else CRclusters.push_back(cl);
	}

	// colour reconnection
	_colourReconnector->rearrange(CRclusters);

	// tag new clusters as children of the partons to hadronize
	_setChildren(CRclusters);

	// recombine vectors of (possibly) reconnected and BV clusters
	clusters.clear();
	clusters.insert( clusters.end(), CRclusters.begin(), CRclusters.end() );
	clusters.insert( clusters.end(), BVclusters.begin(), BVclusters.end() );

	// fission of heavy clusters
	// NB: during cluster fission, light hadrons might be produced straight away
	finalHadrons = _clusterFissioner->fission(clusters,isSoftUnderlyingEventON());

	lightOK = _lightClusterDecayer->decay(clusters,finalHadrons);

	// if the decay of the light clusters was not successful, undo the cluster
	// fission and decay steps and revert to the original state of the event
	// record
	if (!lightOK) {
	clusters = savedclusters;
	for_each(clusters.begin(),
	clusters.end(),
	mem_fun(&Particle::undecay));
	}
	}
	if (!lightOK) {
	// currentHandler_ = 0;
	throw Exception("CluHad::handle(): tried LightClusterDecayer 10 times!",
	Exception::eventerror);
	}


	// decay the remaining clusters
	_clusterDecayer->decay(clusters,finalHadrons);

	// *****************************************
	// *****************************************
	// *****************************************

	StepPtr pstep = newStep();
	set<PPtr> allDecendants;
	for (tPVector::const_iterator it = tagged.begin();
	it != tagged.end(); ++it) {
	extractChildren(*it, allDecendants);
	}

	for(set<PPtr>::const_iterator it = allDecendants.begin();
	it != allDecendants.end(); ++it) {
	// this is a workaround because the set sometimes
	// re-orders parents after their children
	if ((*it)->children().empty())
	pstep->addDecayProduct(*it);
	else {
	pstep->addDecayProduct(*it);
	pstep->addIntermediate(*it);
	}
	}

	// *****************************************
	// *****************************************
	// *****************************************

	// soft underlying event if needed
	if (isSoftUnderlyingEventON()) {
	assert(_underlyingEventHandler);
	ch.performStep(_underlyingEventHandler,Hint::Default());
	}
	- // calculate positions
	- // extract all particles from the event
	- tEventPtr event=ch.currentEvent();
	- vector<tPPtr> particles;
	- particles.reserve(256);
	- event->select(back_inserter(particles), ThePEG::AllSelector());
	- for(vector<tPPtr>::const_iterator pit=particles.begin();
	- pit!=particles.end(); ++pit) {
	- if ((**pit).parents().empty())
	- continue;
	- tPPtr parent = (**pit).parents()[0];
	- // fudged fix for the shower's technical vertices:
	- // make lifelength for 1-1 vertices zero
	- // \todo sort out shower inserted vertices properly
	- if ( (**pit).id() == parent->id()
	- && (**pit).parents().size() == 1
	- && parent->children().size() == 1 ) {
	- parent->setLifeLength(LorentzDistance());
	- // ??? (**pit).setVertex( parent->vertex() );
	- }
	- // fix the vertex position for particles from clusters
	- bool inClusters = false;
	- bool newVertex = false;
	- // iterate up the ancestry to find the origin of the parent cluster
	- while( !parent->parents().empty() ) {
	- bool parentIsCluster = parent->id() == ParticleID::Cluster;
	- if ( !inClusters && parentIsCluster ) {
	- inClusters = true;
	- }
	- else if ( inClusters && !parentIsCluster ) {
	- newVertex = true;
	- break;
	- }
	- parent = parent->parents()[0];
	- }
	- // parent is now the ancestor of the cluster chain
	- if ( newVertex ) (**pit).setVertex( parent->vertex() );
	- }
	}

	void ClusterHadronizationHandler::_setChildren(ClusterVector clusters) const {

	// erase existing information about the partons' children
	tPVector partons;
	for (ClusterVector::const_iterator cl = clusters.begin();
	cl != clusters.end(); cl++) {
	partons.push_back( (*cl)->colParticle() );
	partons.push_back( (*cl)->antiColParticle() );
	}
	for_each(partons.begin(), partons.end(), mem_fun(&Particle::undecay));

	// give new parents to the clusters: their constituents
	for (ClusterVector::iterator cl = clusters.begin();
	cl != clusters.end(); cl++) {
	(cl)->colParticle()->addChild(cl);
	(cl)->antiColParticle()->addChild(cl);
	}

	}

	diff --git a/Hadronization/ColourReconnector.cc b/Hadronization/ColourReconnector.cc
	--- a/Hadronization/ColourReconnector.cc
	+++ b/Hadronization/ColourReconnector.cc
	@@ -1,407 +1,426 @@
	// -- C++ --
	//
	// ColourReconnector.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ColourReconnector class.
	//

	#include "ColourReconnector.h"
	#include "Cluster.h"
	#include "Herwig++/Utilities/Maths.h"
	-
	#include <ThePEG/Interface/Switch.h>
	#include "ThePEG/Interface/Parameter.h"
	#include <ThePEG/Persistency/PersistentOStream.h>
	#include <ThePEG/Persistency/PersistentIStream.h>
	#include <ThePEG/Repository/UseRandom.h>
	#include <algorithm>
	#include <ThePEG/Utilities/DescribeClass.h>
	+#include <ThePEG/Repository/EventGenerator.h>


	using namespace Herwig;

	typedef ClusterVector::iterator CluVecIt;

	DescribeClass<ColourReconnector,Interfaced>
	describeColourReconnector("Herwig::ColourReconnector","");

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

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

	void ColourReconnector::rearrange(ClusterVector & clusters) {
	if (_clreco == 0) return;

	// need at least two clusters
	if (clusters.size() < 2) return;

	// do the colour reconnection
	switch (_algorithm) {
	case 0: _doRecoPlain(clusters);
	break;
	case 1: _doRecoStatistical(clusters);
	break;
	}

	return;
	}


	Energy2 ColourReconnector::_clusterMassSum(const PVector & q,
	const PVector & aq) const {
	const size_t nclusters = q.size();
	assert (aq.size() == nclusters);
	Energy2 sum = ZERO;
	for (size_t i = 0; i < nclusters; i++)
	sum += ( q[i]->momentum() + aq[i]->momentum() ).m2();
	return sum;
	}


	bool ColourReconnector::_containsColour8(const ClusterVector & cv,
	const vector<size_t> & P) const {
	assert (P.size() == cv.size());
	for (size_t i = 0; i < cv.size(); i++) {
	tcPPtr p = cv[i]->colParticle();
	tcPPtr q = cv[P[i]]->antiColParticle();
	if (isColour8(p, q)) return true;
	}
	return false;
	}


	void ColourReconnector::_doRecoStatistical(ClusterVector & cv) const {

	const size_t nclusters = cv.size();

	// initially, enumerate (anti)quarks as given in the cluster vector
	ParticleVector q, aq;
	for (size_t i = 0; i < nclusters; i++) {
	q.push_back( cv[i]->colParticle() );
	aq.push_back( cv[i]->antiColParticle() );
	}

	// annealing scheme
	Energy2 t, delta;
	Energy2 lambda = _clusterMassSum(q,aq);
	const unsigned _ntries = _triesPerStepFactor * nclusters;

	// find appropriate starting temperature by measuring the largest lambda
	// difference in some dry-run random rearrangements
	{
	vector<Energy2> typical;
	for (int i = 0; i < 10; i++) {
	const pair <int,int> toswap = _shuffle(q,aq,5);
	ParticleVector newaq = aq;
	swap (newaq[toswap.first], newaq[toswap.second]);
	Energy2 newlambda = _clusterMassSum(q,newaq);
	typical.push_back( abs(newlambda - lambda) );
	}
	t = _initTemp * Math::median(typical);
	}

	// anneal in up to _annealingSteps temperature steps
	for (unsigned step = 0; step < _annealingSteps; step++) {

	// For this temperature step, try to reconnect _ntries times. Stop the
	// algorithm if no successful reconnection happens.
	unsigned nSuccess = 0;
	for (unsigned it = 0; it < _ntries; it++) {

	// make a random rearrangement
	const unsigned maxtries = 10;
	const pair <int,int> toswap = _shuffle(q,aq,maxtries);
	const int i = toswap.first;
	const int j = toswap.second;

	// stop here if we cannot find any allowed reconfiguration
	if (i == -1) break;

	// create a new antiquark vector with the two partons swapped
	ParticleVector newaq = aq;
	swap (newaq[i], newaq[j]);

	// Check if lambda would decrease. If yes, accept the reconnection. If no,
	// accept it only with a probability given by the current Boltzmann
	// factor. In the latter case we set p = 0 if the temperature is close to
	// 0, to avoid division by 0.
	Energy2 newlambda = _clusterMassSum(q,newaq);
	delta = newlambda - lambda;
	double prob = 1.0;
	if (delta > ZERO) prob = ( abs(t) < 1e-8*MeV2 ) ? 0.0 : exp(-delta/t);
	if (UseRandom::rnd() < prob) {
	lambda = newlambda;
	swap (newaq, aq);
	nSuccess++;
	}
	}
	if (nSuccess == 0) break;

	// reduce temperature
	t *= _annealingFactor;
	}

	// construct the new cluster vector
	ClusterVector newclusters;
	for (size_t i = 0; i < nclusters; i++) {
	ClusterPtr cl = new_ptr( Cluster( q[i], aq[i] ) );
	newclusters.push_back(cl);
	}
	swap(newclusters,cv);
	return;
	}


	void ColourReconnector::_doRecoPlain(ClusterVector & cv) const {

	ClusterVector newcv = cv;

	// try to avoid systematic errors by randomising the reconnection order
	long (*p_irnd)(long) = UseRandom::irnd;
	random_shuffle( newcv.begin(), newcv.end(), p_irnd );

	// iterate over all clusters
	for (CluVecIt cit = newcv.begin(); cit != newcv.end(); cit++) {

	// find the cluster which, if reconnected with *cit, would result in the
	// smallest sum of cluster masses
	// NB this method returns *cit if no reconnection partner can be found
	CluVecIt candidate = _findRecoPartner(cit, newcv);

	// skip this cluster if no possible reshuffling partner can be found
	if (candidate == cit) continue;

	// accept the reconnection with probability _preco.
	if (UseRandom::rnd() < _preco) {

	pair <ClusterPtr,ClusterPtr> reconnected = _reconnect(cit, candidate);

	// Replace the clusters in the ClusterVector. The order of the
	// colour-triplet partons in the cluster vector is retained here.

	// replace *cit by reconnected.first
	*cit = reconnected.first;

	// replace candidate by reconnected.second
	*candidate = reconnected.second;
	}
	}

	swap(cv,newcv);
	return;
	}


	CluVecIt ColourReconnector::_findRecoPartner(CluVecIt cl,
	ClusterVector & cv) const {

	CluVecIt candidate = cl;
	Energy minMass = 1*TeV;
	for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {

	// don't even look at original cluster
	if(cit==cl) continue;

	// don't allow colour octet clusters
	if ( isColour8( (*cl)->colParticle(),
	(*cit)->antiColParticle() ) \|\|
	isColour8( (*cit)->colParticle(),
	(*cl)->antiColParticle() ) ) {
	continue;
	}

	// stop it putting beam remnants together
	if((cl)->isBeamCluster() && (cit)->isBeamCluster()) continue;

	+ // stop it putting far apart clusters together
	+ if(((cl).vertex()-(cit).vertex()).m()>_maxDistance) continue;
	+
	// momenta of the old clusters
	Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() +
	(*cl)->antiColParticle()->momentum();
	Lorentz5Momentum p2 = (*cit)->colParticle()->momentum() +
	(*cit)->antiColParticle()->momentum();

	// momenta of the new clusters
	Lorentz5Momentum p3 = (*cl)->colParticle()->momentum() +
	(*cit)->antiColParticle()->momentum();
	Lorentz5Momentum p4 = (*cit)->colParticle()->momentum() +
	(*cl)->antiColParticle()->momentum();

	Energy oldMass = abs( p1.m() ) + abs( p2.m() );
	Energy newMass = abs( p3.m() ) + abs( p4.m() );

	if ( newMass < oldMass && newMass < minMass ) {
	minMass = newMass;
	candidate = cit;
	}
	}
	return candidate;
	}


	pair <ClusterPtr,ClusterPtr>
	ColourReconnector::_reconnect(ClusterPtr c1, ClusterPtr c2) const {

	// choose the other possibility to form two clusters from the given
	// constituents
	assert(c1->numComponents()==2);
	assert(c2->numComponents()==2);
	int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1);
	for(unsigned int ix=0;ix<2;++ix) {
	if (c1->particle(ix)->hasColour(false)) c1_col = ix;
	else if(c1->particle(ix)->hasColour(true )) c1_anti = ix;
	if (c2->particle(ix)->hasColour(false)) c2_col = ix;
	else if(c2->particle(ix)->hasColour(true )) c2_anti = ix;
	}
	assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0);

	ClusterPtr newCluster1
	= new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) );
	+
	+ newCluster1->setVertex(0.5*( c1->colParticle()->vertex() +
	+ c2->antiColParticle()->vertex() ));
	+
	if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true);
	if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true);

	ClusterPtr newCluster2
	= new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) );
	+
	+ newCluster2->setVertex(0.5*( c2->colParticle()->vertex() +
	+ c1->antiColParticle()->vertex() ));
	+
	if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true);
	if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true);

	return pair <ClusterPtr,ClusterPtr> (newCluster1, newCluster2);
	}


	pair <int,int> ColourReconnector::_shuffle
	(const PVector & q, const PVector & aq, unsigned maxtries) const {

	const size_t nclusters = q.size();
	assert (nclusters > 1);
	assert (aq.size() == nclusters);

	int i, j;
	unsigned tries = 0;
	bool octet;

	do {
	// find two different random integers in the range [0, nclusters)
	i = UseRandom::irnd( nclusters );
	do { j = UseRandom::irnd( nclusters ); } while (i == j);

	// check if one of the two potential clusters would be a colour octet state
	octet = isColour8( q[i], aq[j] ) \|\| isColour8( q[j], aq[i] ) ;
	tries++;
	} while (octet && tries < maxtries);

	if (octet) i = j = -1;
	return make_pair(i,j);
	}


	bool ColourReconnector::isColour8(cPPtr p, cPPtr q) {
	bool octet = false;

	// make sure we have a triplet and an anti-triplet
	if ( ( p->hasColour() && q->hasAntiColour() ) \|\|
	( p->hasAntiColour() && q->hasColour() ) ) {
	if ( !p->parents().empty() && !q->parents().empty() ) {
	// true if p and q are originated from a colour octet
	octet = ( p->parents()[0] == q->parents()[0] ) &&
	( p->parents()[0]->data().iColour() == PDT::Colour8 );
	}
	}

	return octet;
	}


	void ColourReconnector::persistentOutput(PersistentOStream & os) const {
	os << _clreco << _preco << _algorithm << _initTemp << _annealingFactor
	- << _annealingSteps << _triesPerStepFactor;
	+ << _annealingSteps << _triesPerStepFactor << ounit(_maxDistance,femtometer);
	}

	void ColourReconnector::persistentInput(PersistentIStream & is, int) {
	is >> _clreco >> _preco >> _algorithm >> _initTemp >> _annealingFactor
	- >> _annealingSteps >> _triesPerStepFactor;
	+ >> _annealingSteps >> _triesPerStepFactor >> iunit(_maxDistance,femtometer);
	}


	void ColourReconnector::Init() {

	static ClassDocumentation<ColourReconnector> documentation
	("This class is responsible of the colour reconnection.");


	static Switch<ColourReconnector,int> interfaceColourReconnection
	("ColourReconnection",
	"Colour reconnections",
	&ColourReconnector::_clreco, 0, true, false);
	static SwitchOption interfaceColourReconnectionOff
	(interfaceColourReconnection,
	"No",
	"Colour reconnections off",
	0);
	static SwitchOption interfaceColourReconnectionOn
	(interfaceColourReconnection,
	"Yes",
	"Colour reconnections on",
	1);


	static Parameter<ColourReconnector,double> interfaceMtrpAnnealingFactor
	("AnnealingFactor",
	"The annealing factor is the ratio of the temperatures in two successive "
	"temperature steps.",
	&ColourReconnector::_annealingFactor, 0.9, 0.0, 1.0,
	false, false, Interface::limited);

	static Parameter<ColourReconnector,unsigned> interfaceMtrpAnnealingSteps
	("AnnealingSteps",
	"Number of temperature steps in the statistical annealing algorithm",
	&ColourReconnector::_annealingSteps, 50, 1, 10000,
	false, false, Interface::limited);

	static Parameter<ColourReconnector,double> interfaceMtrpTriesPerStepFactor
	("TriesPerStepFactor",
	"The number of reconnection tries per temperature steps is the number of "
	"clusters times this factor.",
	&ColourReconnector::_triesPerStepFactor, 5.0, 0.0, 100.0,
	false, false, Interface::limited);


	static Parameter<ColourReconnector,double> interfaceMtrpInitialTemp
	("InitialTemperature",
	"Factor used to determine the initial temperature from the median of the "
	"energy change in a few random rearrangements.",
	&ColourReconnector::_initTemp, 0.1, 0.00001, 100.0,
	false, false, Interface::limited);


	static Parameter<ColourReconnector,double> interfaceRecoProb
	("ReconnectionProbability",
	"Probability that a found reconnection possibility is actually accepted",
	&ColourReconnector::_preco, 0.5, 0.0, 1.0,
	false, false, Interface::limited);


	static Switch<ColourReconnector,int> interfaceAlgorithm
	("Algorithm",
	"Specifies the colour reconnection algorithm",
	&ColourReconnector::_algorithm, 0, true, false);
	static SwitchOption interfaceAlgorithmPlain
	(interfaceAlgorithm,
	"Plain",
	"Plain colour reconnection as in Herwig++ 2.5.0",
	0);
	static SwitchOption interfaceAlgorithmStatistical
	(interfaceAlgorithm,
	"Statistical",
	"Statistical colour reconnection using simulated annealing",
	1);

	+
	+ static Parameter<ColourReconnector,Length> interfaceMaxDistance
	+ ("MaxDistance",
	+ "Maximum distance between the clusters at which to consider rearrangement"
	+ " to avoid colour reconneections of displaced vertices",
	+ &ColourReconnector::_maxDistance, femtometer, 1000.femtometer, 0.0femtometer, 1e100*femtometer,
	+ false, false, Interface::limited);
	+
	}
	diff --git a/Hadronization/ColourReconnector.h b/Hadronization/ColourReconnector.h
	--- a/Hadronization/ColourReconnector.h
	+++ b/Hadronization/ColourReconnector.h
	@@ -1,235 +1,241 @@
	// -- C++ --
	//
	// ColourReconnector.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_ColourReconnector_H
	#define HERWIG_ColourReconnector_H

	#include <ThePEG/Interface/Interfaced.h>
	#include "CluHadConfig.h"
	#include "ColourReconnector.fh"

	namespace Herwig {


	using namespace ThePEG;

	/** \ingroup Hadronization
	* \class ColourReconnector
	* \brief Class for changing colour reconnections of partons.
	* \author Alberto Ribon, Christian Roehr
	*
	* This class does the nonperturbative colour rearrangement, after the
	* nonperturbative gluon splitting and the "normal" cluster formation.
	* It uses the list of particles in the event record, and the collections of
	* "usual" clusters which is passed to the main method. If the colour
	* reconnection is actually accepted, then the previous collections of "usual"
	* clusters is first deleted and then the new one is created.
	*
	* * @see \ref ColourReconnectorInterfaces "The interfaces"
	* defined for ColourReconnector.
	*/
	class ColourReconnector: public Interfaced {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Default constructor.
	*/
	ColourReconnector() :
	_algorithm(0),
	_annealingFactor(0.9),
	_annealingSteps(50),
	_clreco(0),
	_initTemp(0.1),
	_preco(0.5),
	- _triesPerStepFactor(5.0)
	+ _triesPerStepFactor(5.0),
	+ _maxDistance(1000.*femtometer)
	{}
	//@}

	/**
	* Does the colour rearrangement, starting out from the list of particles in
	* the event record and the collection of "usual" clusters passed as
	* arguments. If the actual rearrangement is accepted, the initial collection of
	* clusters is overridden by the old ones.
	*/
	void rearrange(ClusterVector & clusters);


	private:

	/** PRIVATE MEMBER FUNCTIONS */

	/**
	* @brief Calculates the sum of the squared cluster masses.
	* @arguments q, aq vectors containing the quarks and antiquarks respectively
	* @return Sum of cluster squared masses M^2_{q[i],aq[i]}.
	*/
	Energy2 _clusterMassSum(const PVector & q, const PVector & aq) const;


	/**
	* @brief Examines whether the cluster vector (under the given permutation of
	* the antiquarks) contains colour-octet clusters
	* @param cv Cluster vector
	* @param P Permutation, a vector of permutated indices from 0 to
	* cv.size()-1
	*/
	bool _containsColour8(const ClusterVector & cv, const vector<size_t> & P) const;

	/**
	* @brief A Metropolis-type algorithm which finds a local minimum in the
	* total sum of cluster masses
	* @arguments cv cluster vector
	*/
	void _doRecoStatistical(ClusterVector & cv) const;

	/**
	* @brief Plain colour reconnection as used in Herwig++ 2.5.0
	* @arguments cv cluster vector
	*/
	void _doRecoPlain(ClusterVector & cv) const;

	/**
	* @brief Finds the cluster in cv which, if reconnected with the given
	* cluster cl, would result in the smallest sum of cluster masses.
	* If no reconnection partner can be found, a pointer to the
	* original Cluster cl is returned.
	* @arguments cv cluster vector
	* cl cluster iterator (must be from cv) which wants to have a reconnection partner
	* @return iterator to the found cluster, or the original cluster pointer if
	* no mass-reducing combination can be found
	*/
	ClusterVector::iterator _findRecoPartner(ClusterVector::iterator cl,
	ClusterVector & cv) const;

	/**
	* @brief Reconnects the constituents of the given clusters to the (only)
	* other possible cluster combination.
	* @return pair of pointers to the two new clusters
	*/
	pair <ClusterPtr,ClusterPtr> _reconnect(ClusterPtr c1, ClusterPtr c2) const;

	/**
	* @brief At random, swap two antiquarks, if not excluded by the
	* constraint that there must not be any colour-octet clusters.
	* @arguments q, aq vectors containing the quarks and antiquarks respectively
	* maxtries maximal number of tries to find a non-colour-octet
	* reconfiguration
	* @return Pair of ints indicating the indices of the antiquarks to be
	* swapped. Returns (-1,-1) if no valid reconfiguration could be
	* found after maxtries trials
	*/
	pair <int,int>
	_shuffle(const PVector & q, const PVector & aq, unsigned maxtries = 10) const;


	/** DATA MEMBERS */

	/**
	* Specifies the colour reconnection algorithm to be used.
	*/
	int _algorithm;

	/**
	* The annealing factor is the ratio of two successive temperature steps:
	* T_n = _annealingFactor * T_(n-1)
	*/
	double _annealingFactor;

	/**
	* Number of temperature steps in the statistical annealing algorithm
	*/
	unsigned _annealingSteps;

	/**
	* Do we do colour reconnections?
	*/
	int _clreco;

	/**
	* Factor used to determine the initial temperature according to
	* InitialTemperature = _initTemp * median {energy changes in a few random
	* rearrangements}
	*/
	double _initTemp;

	/**
	* Probability that a found reconnection possibility is actually accepted.
	*/
	double _preco;

	/**
	* The number of tries per temperature steps is the number of clusters times
	* this factor.
	*/
	double _triesPerStepFactor;

	/**
	+ * Maximium distance for reconnections
	+ */
	+ Length _maxDistance;
	+
	+ /**
	* @return true, if the two partons are splitting products of the same
	* gluon
	*/
	static bool isColour8(cPPtr p, cPPtr q);


	public:

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

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

	/**
	* Standard Init function used to initialize the interfaces.
	*/
	static void Init();

	protected:

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

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


	private:

	/**
	* Private and non-existent assignment operator.
	*/
	ColourReconnector & operator=(const ColourReconnector &);


	};


	}

	#endif /* HERWIG_ColourReconnector_H */
	diff --git a/Hadronization/LightClusterDecayer.cc b/Hadronization/LightClusterDecayer.cc
	--- a/Hadronization/LightClusterDecayer.cc
	+++ b/Hadronization/LightClusterDecayer.cc
	@@ -1,417 +1,417 @@
	// -- C++ --
	//
	// LightClusterDecayer.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the LightClusterDecayer class.
	//

	#include "LightClusterDecayer.h"
	#include <ThePEG/Interface/ClassDocumentation.h>
	#include <ThePEG/Interface/Parameter.h>
	#include <ThePEG/Interface/Reference.h>
	#include <ThePEG/Persistency/PersistentOStream.h>
	#include <ThePEG/Persistency/PersistentIStream.h>
	#include <ThePEG/PDT/EnumParticles.h>
	#include <ThePEG/Repository/EventGenerator.h>
	#include "Cluster.h"
	#include "CheckId.h"
	#include "Herwig++/Utilities/Kinematics.h"
	#include <ThePEG/Utilities/DescribeClass.h>

	using namespace Herwig;

	DescribeClass<LightClusterDecayer,Interfaced>
	describeLightClusterDecayer("Herwig::LightClusterDecayer","");

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

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


	void LightClusterDecayer::persistentOutput(PersistentOStream & os) const {
	os << _hadronSelector << _limBottom << _limCharm << _limExotic;
	}

	void LightClusterDecayer::persistentInput(PersistentIStream & is, int) {
	is >> _hadronSelector >> _limBottom >> _limCharm >> _limExotic ;
	}

	void LightClusterDecayer::Init() {

	static ClassDocumentation<LightClusterDecayer> documentation
	("There is the class responsible for the one-hadron decay of light clusters");

	static Reference<LightClusterDecayer,HadronSelector>
	interfaceHadronSelector("HadronSelector",
	"A reference to the HadronSelector object",
	&Herwig::LightClusterDecayer::_hadronSelector,
	false, false, true, false);

	static Parameter<LightClusterDecayer,double>
	interfaceSingleHadronLimitBottom ("SingleHadronLimitBottom","threshold for one-hadron decay of b-cluster",
	&LightClusterDecayer::_limBottom, 0, 0.0, 0.0, 100.0,false,false,false);
	static Parameter<LightClusterDecayer,double>
	interfaceSingleHadronLimitCharm ("SingleHadronLimitCharm","threshold for one-hadron decay of c-cluster",
	&LightClusterDecayer::_limCharm, 0, 0.0, 0.0, 100.0,false,false,false);
	static Parameter<LightClusterDecayer,double>
	interfaceSingleHadronLimitExotic ("SingleHadronLimitExotic","threshold for one-hadron decay of exotic cluster",
	&LightClusterDecayer::_limExotic, 0, 0.0, 0.0, 100.0,false,false,false);

	}


	bool LightClusterDecayer::decay(ClusterVector & clusters, tPVector & finalhadrons) {

	// Loop over all clusters, and for those that were not heavy enough
	// to undergo to fission, check if they are below the threshold
	// for normal two-hadron decays. If this is the case, then the cluster
	// should be decayed into a single hadron: this can happen only if
	// it is possible to reshuffle momenta between the cluster and
	// another one; in the rare occasions in which such exchange of momenta
	// is not possible (because all of the clusters are too light) then
	// the event is skipped.
	// Notice that, differently from what happens in Fortran Herwig,
	// light (that is below the threshold for the production of the lightest
	// pair of hadrons with the proper flavours) fission products, produced
	// by the fission of heavy clusters in class ClusterFissioner
	// have been already "decayed" into single hadron (the lightest one
	// with proper flavour) by the same latter class, without requiring
	// any reshuffling. Therefore the light clusters that are treated in
	// this LightClusterDecayer class are produced directly
	// (originally) by the class ClusterFinder.

	// To preserve all of the information, the cluster partner with which
	// the light cluster (that decays into a single hadron) exchanges
	// momentum in the reshuffling procedure is redefined and inserted
	// in the vector vecNewRedefinedCluPtr. Only at the end, when all
	// light clusters have been examined, the elements this vector will be
	// copied in collecCluPtr (the reason is that it is not allowed to
	// modify a STL container while iterating over it. At the same time,
	// this ensures that a cluster can be redefined only once, which seems
	// sensible although not strictly necessary).
	// Notice that the cluster reshuffling partner is normally redefined
	// and inserted in the vector vecNewRedefinedCluPtr, but not always:
	// in the case it is also light, then it is also decayed immediately
	// into a single hadron, without redefining it (the reason being that,
	// otherwise, the would-be redefined cluster could have undefined
	// components).
	vector<tClusterPtr> redefinedClusters;

	for (ClusterVector::const_iterator it = clusters.begin();
	it != clusters.end(); ++it) {

	// Skip the clusters that are not available or that are
	// heavy, intermediate, clusters that have undergone to fission,
	if ( ! (it)->isAvailable() \|\| ! (it)->isReadyToDecay() ) continue;

	// We need to require (at least at the moment, maybe in the future we
	// could change it) that the cluster has exactly two components,
	// because otherwise we don't know how to deal with the kinematics.
	// If this is not the case, then send a warning because it is not suppose
	// to happen, and then do nothing with (ignore) such cluster.
	if ( (*it)->numComponents() != 2 ) {
	generator()->logWarning( Exception("LightClusterDecayer::decay "
	"***Still cluster with not exactly"
	" 2 components*** ",
	Exception::warning) );
	continue;
	}

	// Extract the particle pointer of the two components of the cluster.

	tPPtr ptrQ1 = (*it)->particle(0);
	tPPtr ptrQ2 = (*it)->particle(1);

	tcPDPtr par1 = ptrQ1->dataPtr();
	tcPDPtr par2 = ptrQ2->dataPtr();

	// Determine the sum of the nominal masses of the two lightest hadrons
	// with the right flavour numbers as the cluster under consideration.
	// Notice that we don't need real masses (drawn by a Breit-Wigner
	// distribution) because the lightest pair of hadrons does not involve
	// any broad resonance.
	Energy threshold = _hadronSelector->massLightestHadronPair(par1,par2);
	// Special: it allows one-hadron decays also above threshold.
	if (CheckId::isExotic(par1,par2))
	threshold = (1.0 + UseRandom::rnd()_limExotic);
	else if (CheckId::hasBottom(par1,par2))
	threshold = (1.0 + UseRandom::rnd()_limBottom);
	else if (CheckId::hasCharm(par1,par2))
	threshold = (1.0 + UseRandom::rnd()_limCharm);


	// only do one hadron decay is mass less than the threshold
	if((*it)->mass()>=threshold) continue;

	tcPDPtr hadron= _hadronSelector->lightestHadron(par1,par2);

	// We assume that the candidate reshuffling cluster partner,
	// with whom the light cluster can exchange momenta,
	// is chosen as the closest in space-time between the available
	// clusters. Notice that an alternative, sensible approach
	// could be to consider instead the "closeness" in the colour
	// structure...
	// Notice that nor a light cluster (which decays into a single hadron)
	// neither its cluster reshuffling partner (which either has a
	// redefined cluster or also decays into a single hadron) can be
	// a reshuffling partner of another light cluster.
	// This because we are requiring that the considered candidate cluster
	// reshuffling partner has the status "isAvailable && isReadyToDecay" true;
	// furthermore, the new redefined clusters are not added to the collection
	// of cluster before the end of the entire reshuffling procedure, avoiding
	// in this way that the redefined cluster of a cluster reshuffling partner
	// is used again later. Needless to say, this is just an assumption,
	// although reasonable, but nothing more than that!

	// Build a multimap of available reshuffling cluster partners,
	// with key given by the module of the invariant space-time distance
	// w.r.t. the light cluster, so that this new collection is automatically
	// ordered in increasing distance values.
	// We use a multimap, rather than a map, just for precaution against not properly
	// defined cluster positions which could produce all identical (null) distances.
	multimap<Length,tClusterPtr> candidates;
	for ( ClusterVector::iterator jt = clusters.begin();
	jt != clusters.end(); ++jt ) {
	if ((jt)->isAvailable() && (jt)->isReadyToDecay() && jt != it) {
	Length distance = abs (((it)->vertex() - (jt)->vertex()).m());
	candidates.insert(pair<Length,tClusterPtr>(distance,*jt));
	}
	}

	// Loop sequentially the multimap.
	multimap<Length,tClusterPtr>::const_iterator mmapIt = candidates.begin();
	bool found = false;
	while (!found && mmapIt != candidates.end()) {
	found = reshuffling(hadron, it, (mmapIt).second, redefinedClusters, finalhadrons);
	if (!found) ++mmapIt;
	}

	if (!found) return partonicReshuffle(hadron,*it,finalhadrons);
	} // end loop over collecCluPtr

	// Add to collecCluPtr all of the redefined new clusters (indeed the
	// pointers to them are added) contained in vecNewRedefinedCluPtr.
	for (tClusterVector::const_iterator it = redefinedClusters.begin();
	it != redefinedClusters.end(); ++it) {
	clusters.push_back(*it);
	}
	return true;
	}


	bool LightClusterDecayer::reshuffling(const tcPDPtr pdata1,
	tClusterPtr cluPtr1,
	tClusterPtr cluPtr2,
	tClusterVector & redefinedClusters,
	tPVector & finalhadrons)
	{
	// don't reshuffle with beam clusters
	if(cluPtr2->isBeamCluster()) return false;

	// This method does the reshuffling of momenta between the cluster "1",
	// that must decay into a single hadron (with id equal to idhad1), and
	// the candidate cluster "2". It returns true if the reshuffling succeed,
	// false otherwise.

	PPtr ptrhad1 = pdata1->produceParticle();
	if ( ! ptrhad1 ) {
	generator()->logWarning( Exception("LightClusterDecayer::reshuffling"
	"*Cannot create a particle with specified id*",
	Exception::warning) );
	return false;
	}
	Energy mhad1 = ptrhad1->mass();

	// Let's call "3" and "4" the two constituents of the second cluster
	tPPtr part3 = cluPtr2->particle(0);
	tPPtr part4 = cluPtr2->particle(1);

	// Check if the system of the two clusters can kinematically be replaced by
	// an hadron of mass mhad1 (which is the lightest single hadron with the
	// same flavour numbers as the first cluster) and the second cluster.
	// If not, then try to replace the second cluster with the lightest hadron
	// with the same flavour numbers; if it still fails, then give up!
	Lorentz5Momentum pSystem = cluPtr1->momentum() + cluPtr2->momentum();
	pSystem.rescaleMass(); // set the mass as the invariant of the quadri-vector
	Energy mSystem = pSystem.mass();
	Energy mclu2 = cluPtr2->mass();
	bool singleHadron = false;
	Energy mLHP2 = _hadronSelector->massLightestHadronPair(part3->dataPtr(),part4->dataPtr());
	Energy mLH2 = _hadronSelector->massLightestHadron(part3->dataPtr(),part4->dataPtr());

	if(mSystem > mhad1 + mclu2 && mclu2 > mLHP2) { singleHadron = false; }
	else if(mSystem > mhad1 + mLH2) { singleHadron = true; mclu2 = mLH2; }
	else return false;
	// Let's call from now on "Sys" the system of the two clusters, and
	// had1 (of mass mhad1) the lightest hadron in which the first
	// cluster decays, and clu2 (of mass mclu2) either the second
	// cluster or the lightest hadron in which it decays (depending
	// which one is kinematically allowed, see above).
	// The idea behind the reshuffling is to replace the system of the
	// two clusters by the system of the hadron had1 and (cluster or hadron) clu2,
	// but leaving the overall system unchanged. Furthermore, the motion
	// of had1 and clu2 in the Sys frame is assumed to be parallel to, respectively,
	// those of the original cluster1 and cluster2 in the same Sys frame.

	// Calculate the unit three-vector, in the frame "Sys" along which the
	// two initial clusters move.
	Lorentz5Momentum u( cluPtr1->momentum() );
	u.boost( - pSystem.boostVector() ); // boost from LAB to Sys

	// Calculate the momenta of had1 and clu2 in the Sys frame first,
	// and then boost back in the LAB frame.
	Lorentz5Momentum phad1, pclu2;

	if (pSystem.m() < mhad1 + mclu2 ) {
	throw Exception() << "Impossible Kinematics in LightClusterDecayer::reshuffling()"
	<< Exception::eventerror;
	}

	Kinematics::twoBodyDecay(pSystem, mhad1, mclu2, u.vect().unit(), phad1, pclu2);

	- ptrhad1->set5Momentum( phad1 ); // set momentum of first hadron.
	- ptrhad1->setLabVertex(cluPtr1->vertex()); // set hadron vertex position to the
	- // parent cluster position.
	+ ptrhad1->set5Momentum( phad1 ); // set momentum of first hadron.
	+ ptrhad1->setVertex(cluPtr1->vertex()); // set hadron vertex position to the
	+ // parent cluster position.
	cluPtr1->addChild(ptrhad1);
	finalhadrons.push_back(ptrhad1);
	cluPtr1->flagAsReshuffled();
	cluPtr2->flagAsReshuffled();


	if(singleHadron) {

	// In the case that also the cluster reshuffling partner is light
	// it is decayed into a single hadron, without creating the
	// redefined cluster (this choice is justified in order to avoid
	// clusters that could have undefined components).

	PPtr ptrhad2 = _hadronSelector->lightestHadron(part3->dataPtr(),part4->dataPtr())
	->produceParticle();
	ptrhad2->set5Momentum( pclu2 );
	- ptrhad2->setLabVertex( cluPtr2->vertex() ); // set hadron vertex position to the
	- // parent cluster position.
	+ ptrhad2->setVertex( cluPtr2->vertex() ); // set hadron vertex position to the
	+ // parent cluster position.
	cluPtr2->addChild(ptrhad2);
	finalhadrons.push_back(ptrhad2);
	} else {

	// Create the new cluster which is the redefinitions of the cluster
	// partner (cluster "2") used in the reshuffling procedure of the
	// light cluster (cluster "1").
	// The rationale of this is to preserve completely all of the information.
	ClusterPtr cluPtr2new = ClusterPtr();
	if(part3 && part4) cluPtr2new = new_ptr(Cluster(part3,part4));

	cluPtr2new->set5Momentum( pclu2 );
	cluPtr2new->setVertex( cluPtr2->vertex() );
	cluPtr2->addChild( cluPtr2new );
	redefinedClusters.push_back( cluPtr2new );

	// Set consistently the momenta of the two components of the second cluster
	// after the reshuffling. To do that we first calculate the momenta of the
	// constituents in the initial cluster rest frame; then we boost them back
	// in the lab but using this time the new cluster rest frame. Finally we store
	// these information in the new cluster. Notice that we do not set
	// consistently also the momenta of the (eventual) particles pointed by the
	// two components: that's because we do not need to do so, being the momentum
	// an explicit private member of the class Component (which is set equal
	// to the momentum of the eventual particle pointed only in the constructor,
	// but then later should not necessary be the same), and furthermore it allows
	// us not to loose any information, in the sense that we can always, later on,
	// to find the original momenta of the two components before the reshuffling.
	Lorentz5Momentum p3 = part3->momentum(); //p3new->momentum();
	p3.boost( - (cluPtr2->momentum()).boostVector() ); // from LAB to clu2 (old) frame
	p3.boost( pclu2.boostVector() ); // from clu2 (new) to LAB frame
	Lorentz5Momentum p4 = part4->momentum(); //p4new->momentum();
	p4.boost( - (cluPtr2->momentum()).boostVector() ); // from LAB to clu2 (old) frame
	p4.boost( pclu2.boostVector() ); // from clu2 (new) to LAB frame
	cluPtr2new->particle(0)->set5Momentum(p3);
	cluPtr2new->particle(1)->set5Momentum(p4);
	} // end of if (singleHadron)
	return true;
	}

	bool LightClusterDecayer::partonicReshuffle(const tcPDPtr had,
	const PPtr cluster,
	tPVector & finalhadrons) {
	tPPtr meson(cluster);
	if(!meson->parents().empty()) meson=meson->parents()[0];
	if(!meson->parents().empty()) meson=meson->parents()[0];
	// check b/c hadron decay
	int ptype(abs(meson->id())%10000);
	bool heavy = (ptype/1000 == 5 \|\| ptype/1000 ==4 );
	heavy \|= (ptype/100 == 5 \|\| ptype/100 ==4 );
	heavy \|= (ptype/10 == 5 \|\| ptype/10 ==4 );
	if(!heavy) return false;
	// find the leptons
	tPVector leptons;
	for(unsigned int ix=0;ix<meson->children().size();++ix) {
	if(!(meson->children()[ix]->dataPtr()->coloured())) {
	leptons.push_back(meson->children()[ix]);
	}
	}
	if(leptons.size()==1) {
	tPPtr w=leptons[0];
	leptons.pop_back();
	for(unsigned int ix=0;ix<w->children().size();++ix) {
	if(!w->children()[ix]->dataPtr()->coloured()) {
	leptons.push_back(w->children()[ix]);
	}
	}
	}
	if(leptons.size()!=2) return false;
	// get momentum of leptonic system and the its minimum possible mass
	Energy mmin(ZERO);
	Lorentz5Momentum pw;
	for(unsigned int ix=0;ix<leptons.size();++ix) {
	pw+=leptons[ix]->momentum();
	mmin+=leptons[ix]->mass();
	}
	pw.rescaleMass();
	// check we can do the reshuffling
	PPtr ptrhad = had->produceParticle();
	// total momentum fo the system
	Lorentz5Momentum pSystem = pw + cluster->momentum();
	pSystem.rescaleMass();
	// normal case get additional energy by rescaling momentum in rest frame of
	// system
	if(pSystem.mass()>ptrhad->mass()+pw.mass()&&pw.mass()>mmin) {
	// Calculate the unit three-vector, in the frame "Sys" along which the
	// two initial clusters move.
	Lorentz5Momentum u(cluster->momentum());
	u.boost( - pSystem.boostVector() );
	// Calculate the momenta of had1 and clu2 in the Sys frame first,
	// and then boost back in the LAB frame.
	Lorentz5Momentum phad1, pclu2;
	Kinematics::twoBodyDecay(pSystem, ptrhad->mass(), pw.mass(),
	u.vect().unit(), phad1, pclu2);
	// set momentum of first hadron.
	ptrhad->set5Momentum( phad1 );
	// set hadron vertex position to the parent cluster position.
	ptrhad->setLabVertex(cluster->vertex());
	// add hadron
	cluster->addChild(ptrhad);
	finalhadrons.push_back(ptrhad);
	// reshuffle the leptons
	// boost the leptons to the rest frame of the system
	Boost boost1(-pw.boostVector());
	Boost boost2( pclu2.boostVector());
	for(unsigned int ix=0;ix<leptons.size();++ix) {
	leptons[ix]->deepBoost(boost1);
	leptons[ix]->deepBoost(boost2);
	}
	return true;
	}
	else {
	return false;
	}
	}
	diff --git a/Hadronization/PartonSplitter.cc b/Hadronization/PartonSplitter.cc
	--- a/Hadronization/PartonSplitter.cc
	+++ b/Hadronization/PartonSplitter.cc
	@@ -1,132 +1,138 @@
	// -- C++ --
	//
	// PartonSplitter.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the PartonSplitter class.
	//

	#include "PartonSplitter.h"
	#include <ThePEG/Interface/ClassDocumentation.h>
	#include <ThePEG/Interface/Reference.h>
	#include <ThePEG/Persistency/PersistentOStream.h>
	#include <ThePEG/Persistency/PersistentIStream.h>
	#include <ThePEG/PDT/EnumParticles.h>
	#include <ThePEG/EventRecord/Step.h>
	#include <ThePEG/Repository/EventGenerator.h>
	#include <ThePEG/Repository/CurrentGenerator.h>
	#include "Herwig++/Utilities/Kinematics.h"
	#include <ThePEG/Utilities/DescribeClass.h>

	using namespace Herwig;

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

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

	void PartonSplitter::persistentOutput(PersistentOStream & os) const {
	os << _quarkSelector;
	}

	void PartonSplitter::persistentInput(PersistentIStream & is, int) {
	is >> _quarkSelector;
	}

	DescribeClass<PartonSplitter,Interfaced>
	describePartonSplitter("Herwig::PartonSplitter","");

	void PartonSplitter::Init() {

	static ClassDocumentation<PartonSplitter> documentation
	("This class is reponsible of the nonperturbative splitting of partons");

	}

	void PartonSplitter::split(PVector & tagged) {
	PVector newtag;
	Energy2 Q02 = 0.99*sqr(getParticleData(ParticleID::g)->constituentMass());
	// Loop over all of the particles in the event.
	for(PVector::const_iterator pit = tagged.begin(); pit!=tagged.end(); ++pit) {
	// only considering gluons so add other particles to list of particles
	if( (**pit).data().id() != ParticleID::g ) {
	newtag.push_back(*pit);
	continue;
	}
	// should not have been called for massless or space-like gluons
	if((*pit).momentum().m2() <= 0.0sqr(MeV) ) {
	throw Exception()
	<< "Spacelike or massless gluon m2= " << (**pit).momentum().m2()/GeV2
	<< "GeV2 in PartonSplitter::split()"
	<< Exception::eventerror;
	}
	// time like gluon gets split
	PPtr ptrQ = PPtr();
	PPtr ptrQbar = PPtr();
	splitTimeLikeGluon(*pit,ptrQ,ptrQbar);
	ptrQ->scale(Q02);
	ptrQbar->scale(Q02);

	(*pit)->colourLine()->addColoured(ptrQ);
	(*pit)->addChild(ptrQ);
	newtag.push_back(ptrQ);

	(*pit)->antiColourLine()->addAntiColoured(ptrQbar);
	(*pit)->addChild(ptrQbar);
	newtag.push_back(ptrQbar);
	+
	+ // assume same position as gluon
	+ ptrQ ->setVertex((**pit).decayVertex());
	+ ptrQ ->setLifeLength(Lorentz5Distance());
	+ ptrQbar->setVertex((**pit).decayVertex());
	+ ptrQbar->setLifeLength(Lorentz5Distance());
	}
	swap(tagged,newtag);
	}

	void PartonSplitter::splitTimeLikeGluon(tcPPtr ptrGluon,
	PPtr & ptrQ,
	PPtr & ptrQbar){
	// select the quark flavour
	tPDPtr quark = _quarkSelector.select(UseRandom::rnd());
	// Solve the kinematics of the two body decay G --> Q + Qbar
	Lorentz5Momentum momentumQ;
	Lorentz5Momentum momentumQbar;
	double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 );
	using Constants::pi;
	double phiStar = UseRandom::rnd( -pi , pi );
	Energy constituentQmass = quark->constituentMass();

	if (ptrGluon->momentum().m() < 2.0*constituentQmass) {
	throw Exception() << "Impossible Kinematics in PartonSplitter::splitTimeLikeGluon()"
	<< Exception::eventerror;
	}
	Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass,
	constituentQmass, cosThetaStar, phiStar, momentumQ,
	momentumQbar );
	// Create quark and anti-quark particles of the chosen flavour
	// and set they 5-momentum (the mass is the constituent one).
	ptrQ = new_ptr(Particle(quark ));
	ptrQbar = new_ptr(Particle(quark->CC()));
	ptrQ ->set5Momentum( momentumQ );
	ptrQbar ->set5Momentum( momentumQbar );
	}

	void PartonSplitter::doinit() {
	Interfaced::doinit();
	// calculate the probabilties for the gluon to branch into each quark type
	// based on the available phase-space, as in fortran.
	Energy mg=getParticleData(ParticleID::g)->constituentMass();
	for( int ix=1; ix<6; ++ix ) {
	PDPtr quark = getParticleData(ix);
	Energy pcm = Kinematics::pstarTwoBodyDecay(mg,quark->constituentMass(),
	quark->constituentMass());
	if(pcm>ZERO) _quarkSelector.insert(pcm/GeV,quark);
	}
	if(_quarkSelector.empty())
	throw InitException() << "At least one quark must have constituent mass less "
	<< "then the constituent mass of the gluon in "
	<< "PartonSplitter::doinit()" << Exception::runerror;
	}
	diff --git a/Shower/Base/ShowerTree.cc b/Shower/Base/ShowerTree.cc
	--- a/Shower/Base/ShowerTree.cc
	+++ b/Shower/Base/ShowerTree.cc
	@@ -1,1125 +1,1177 @@
	// -- C++ --
	//
	// ShowerTree.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#include "ShowerProgenitor.h"
	#include "ThePEG/EventRecord/MultiColour.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::isolateLine(vector<PPair>::const_iterator cit,
	vector<PPair> & particles,
	tcColinePtr oldline,
	tColinePtr newline) {
	// loop over particles
	for(vector<PPair>::const_iterator cjt=particles.begin();
	cjt!=particles.end();++cjt) {
	if(cjt==cit) continue;
	// if particle has colour line
	if((*cjt).second->colourLine()) {
	// if only one check if current line and reset
	if(int((*cjt).second->colourInfo()->colourLines().size())==1) {
	if((*cjt).second->colourLine()==oldline)
	newline->addColoured((*cjt).first);
	}
	// if more than one check if each line current line and reset
	else {
	Ptr<MultiColour>::pointer colour1 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>
	((*cjt).second->colourInfo());
	Ptr<MultiColour>::pointer colour2 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>
	((*cjt).first ->colourInfo());
	for(unsigned int ix=0;ix<colour1->colourLines().size();++ix) {
	if(colour1->colourLines()[ix]==oldline)
	colour2->colourLine(newline,int(ix)+1);
	}
	}
	}
	// if particle has anticolour line
	if((*cjt).second->antiColourLine()) {
	// if only one check if current line and reset
	if(int((*cjt).second->colourInfo()->antiColourLines().size())==1) {
	if((*cjt).second->antiColourLine()==oldline)
	newline->addColoured((*cjt).first,true);
	}
	// if more than one check if each line current line and reset
	else {
	Ptr<MultiColour>::pointer colour1 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>
	((*cjt).second->colourInfo());
	Ptr<MultiColour>::pointer colour2 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>
	((*cjt).first ->colourInfo());
	for(unsigned int ix=0;ix<colour1->antiColourLines().size();++ix) {
	if(colour1->antiColourLines()[ix]==oldline)
	colour2->antiColourLine(newline, int(ix)+1);
	}
	}
	}
	}
	}

	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;
	// make the colour connections of the copies
	for(cit=particles.begin();cit!=particles.end();++cit) {
	if((*cit).first->colourInfo()) {
	if((*cit).first->dataPtr()->iColour() == PDT::Colour6 \|\|
	(*cit).first->dataPtr()->iColour() == PDT::Colour6bar)
	(*cit).first->colourInfo(new_ptr(MultiColour()));
	else
	(*cit).first->colourInfo(new_ptr(ColourBase()));
	}
	}
	map<tcColinePtr,tColinePtr> cmap;
	// make the colour connections of the copies
	// loop over the particles
	for(cit=particles.begin();cit!=particles.end();++cit) {
	// if particle has at least one colour line
	if((*cit).second->colourLine()) {
	// one and only one line
	if(int((*cit).second->colourInfo()->colourLines().size())==1) {
	// if not already change
	if(!(*cit).first->colourLine()) {
	// make new line
	tcColinePtr oldline=(*cit).second->colourLine();
	ColinePtr newline=ColourLine::create((*cit).first);
	cmap[oldline]=newline;
	isolateLine(cit,particles,oldline,newline);
	}
	}
	// more than one line
	else {
	Ptr<MultiColour>::pointer colour1 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>
	((*cit).second->colourInfo());
	vector<tcColinePtr> lines1 = colour1->colourLines();
	Ptr<MultiColour>::pointer colour2 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>
	((*cit).first->colourInfo());
	vector<tcColinePtr> lines2 = colour2->colourLines();
	// loop over lines
	for(unsigned int ix=0;ix<lines1.size();++ix) {
	if( (lines2.size()>ix && !lines2[ix]) \|\|
	lines2.size()<=ix) {
	tcColinePtr oldline = lines1[ix];
	ColinePtr newline = new_ptr(ColourLine());
	cmap[oldline]=newline;
	colour2->colourLine(newline, int(ix)+1);
	isolateLine(cit,particles,oldline,newline);
	}
	}
	}
	}
	// if anticolour line
	if((*cit).second->antiColourLine()) {
	// one and only one line
	if(int((*cit).second->colourInfo()->antiColourLines().size())==1) {
	// if not already change
	if(!(*cit).first->antiColourLine()) {
	// make new line
	tcColinePtr oldline=(*cit).second->antiColourLine();
	ColinePtr newline=ColourLine::create((*cit).first, true);
	cmap[oldline]=newline;
	isolateLine(cit,particles,oldline,newline);
	}
	}
	// more than one line
	else {
	Ptr<MultiColour>::pointer colour1 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>
	((*cit).second->colourInfo());
	vector<tcColinePtr> lines1 = colour1->antiColourLines();
	Ptr<MultiColour>::pointer colour2 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>
	((*cit).first->colourInfo());
	vector<tcColinePtr> lines2 = colour2->antiColourLines();
	// loop over lines
	for(unsigned int ix=0;ix<lines1.size();++ix) {
	if( (lines2.size()>ix && !lines2[ix]) \|\|
	lines2.size()<=ix) {
	tcColinePtr oldline = lines1[ix];
	ColinePtr newline = new_ptr(ColourLine());
	cmap[oldline]=newline;
	colour2->antiColourLine(newline, int(ix)+1);
	isolateLine(cit,particles,oldline,newline);
	}
	}
	}
	}
	}

	// 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<tcColinePtr,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::mapColour(PPtr original,
	PPtr copy) {
	// has colour line
	if(copy->colourLine()) {
	// one and only one
	if(copy->colourInfo()->colourLines().size()==1) {
	_colour.insert(make_pair(copy->colourLine(),
	original->colourLine()));
	}
	// more than one
	else {
	Ptr<MultiColour>::pointer colour1 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>(copy->colourInfo());
	vector<tcColinePtr> lines1 = colour1->colourLines();
	Ptr<MultiColour>::pointer colour2 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>(original->colourInfo());
	vector<tcColinePtr> lines2 = colour2->colourLines();
	for(unsigned int ix=0;ix<lines1.size();++ix)
	_colour.insert(make_pair(const_ptr_cast<ColinePtr>(lines1[ix]),
	const_ptr_cast<ColinePtr>(lines2[ix])));
	}
	}
	// has anticolour line
	if(copy->antiColourLine()) {
	// one and only one
	if(copy->colourInfo()->antiColourLines().size()==1) {
	_colour.insert(make_pair(copy->antiColourLine(),
	original->antiColourLine()));
	}
	// more than one
	else {
	Ptr<MultiColour>::pointer colour1 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>(copy->colourInfo());
	vector<tcColinePtr> lines1 = colour1->antiColourLines();

	Ptr<MultiColour>::pointer colour2 =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>(original->colourInfo());
	vector<tcColinePtr> lines2 = colour2->antiColourLines();

	for(unsigned int ix=0;ix<lines1.size();++ix)
	_colour.insert(make_pair(const_ptr_cast<ColinePtr>(lines1[ix]),
	const_ptr_cast<ColinePtr>(lines2[ix])));
	}
	}
	}

	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;
	mapColour(cit->first->original(),cit->first->copy());
	}
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
	for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) {
	if(!cjt->first->perturbative()) continue;
	mapColour(cjt->first->original(),cjt->first->copy());
	}
	// initial-state radiation
	if(ISR) {
	for(cit=incomingLines().begin();cit!=incomingLines().end();++cit) {
	ShowerParticlePtr init=(*cit).first->progenitor();
	assert(init->thePEGBase());
	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);
	+ copy->setLifeLength(Lorentz5Distance());
	+ copy->setVertex(LorentzPoint());
	// 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);
	}
	}
	+ init->setLifeLength(Lorentz5Distance());
	+ init->setVertex(LorentzPoint());
	}
	}
	else {
	for(cit=incomingLines().begin();cit!=incomingLines().end();++cit) {
	ShowerParticlePtr init=(*cit).first->progenitor();
	assert(init->thePEGBase());
	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);
	+ init->setLifeLength(Lorentz5Distance());
	+ init->setVertex(LorentzPoint());
	+ original->setLifeLength(Lorentz5Distance());
	+ original->setVertex(LorentzPoint());
	}
	}
	// final-state radiation
	for(cjt=outgoingLines().begin();cjt!=outgoingLines().end();++cjt) {
	ShowerParticlePtr init=(*cjt).first->progenitor();
	assert(init->thePEGBase());
	+ // ZERO the distance of original
	+ (*cjt).first->original()->setLifeLength(Lorentz5Distance());
	+ (*cjt).first->original()->setVertex(LorentzPoint());
	// 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);
	}
	+ // ZERO the distance of copy ??? \todo change if add space-time
	+ (*cjt).first->copy()->setLifeLength(Lorentz5Distance());
	+ (*cjt).first->copy()->setVertex(LorentzPoint());
	}
	+ // copy so travels no distance
	+ init->setLifeLength(Lorentz5Distance());
	+ init->setVertex(init->parents()[0]->decayVertex());
	+ // sort out the colour
	updateColour(init);
	// insert shower products
	addFinalStateShower(init,pstep);
	}
	_colour.clear();
	}

	void ShowerTree::addFinalStateShower(PPtr p, StepPtr s) {
	- if(p->children().empty()) return;
	+ if(p->children().empty()) {
	+ p->setLifeLength(Lorentz5Distance());
	+ return;
	+ }
	+ // \todo the space time-distance should be set properly here !!!!
	+ else {
	+ p->setLifeLength(Lorentz5Distance());
	+ }
	ParticleVector::const_iterator child;
	for(child=p->children().begin(); child != p->children().end(); ++child) {
	updateColour(*child);
	s->addDecayProduct(*child);
	+ (**child).setVertex(p->decayVertex());
	addFinalStateShower(*child,s);
	}
	}

	void ShowerTree::updateColour(PPtr particle) {
	// if attached to a colour line
	if(particle->colourLine()) {
	// one and only one
	if(particle->colourInfo()->colourLines().size()==1) {
	bool reset=false;
	// if colour line from hard process reconnect
	ColinePtr c1=particle->colourLine();
	if(_colour.find(c1)!=_colour.end()) {
	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);
	}
	}
	else {
	Ptr<MultiColour>::pointer colour =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>(particle->colourInfo());
	vector<tcColinePtr> lines = colour->colourLines();
	for(unsigned int ix=0;ix<lines.size();++ix) {
	ColinePtr c1 = const_ptr_cast<ColinePtr>(lines[ix]);
	if(_colour.find(c1)!=_colour.end()) {
	colour->colourLine(_colour[c1],int(ix)+1);
	c1->removeColoured(particle);
	}
	}
	}
	}
	// if attached to an anticolour line
	if(particle->antiColourLine()) {
	bool reset=false;
	// one and only one
	if(particle->colourInfo()->antiColourLines().size()==1) {
	// if anti colour line from hard process reconnect
	ColinePtr c1=particle->antiColourLine();
	if(_colour.find(c1)!=_colour.end()) {
	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);
	}
	}
	else {
	Ptr<MultiColour>::pointer colour =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>(particle->colourInfo());
	vector<tcColinePtr> lines = colour->antiColourLines();
	for(unsigned int ix=0;ix<lines.size();++ix) {
	ColinePtr c1 = const_ptr_cast<ColinePtr>(lines[ix]);
	if(_colour.find(c1)!=_colour.end()) {
	colour->antiColourLine(_colour[c1],int(ix)+1);
	c1->removeColoured(particle,true);
	}
	}
	}
	}
	}

	void ShowerTree::addInitialStateShower(PPtr p, PPtr hadron,
	StepPtr s, bool addchildren) {
	+ p->setLifeLength(Lorentz5Distance());
	+ p->setVertex(LorentzPoint());
	// 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);
	+ (**child).setLifeLength(Lorentz5Distance());
	+ (**child).setVertex(p->vertex());
	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();
	mapColour(final,copy);
	+ // now this is the ONE instance of the particle which should have a life length
	+ // \todo change if space-time picture added
	+ // set the lifelength, need this so that still in right direction after
	+ // any ISR recoils
	+ Length ctau = copy->lifeTime();
	+ Lorentz5Distance lifeLength(ctau,final->momentum().vect()*(ctau/final->mass()));
	+ final->setLifeLength(lifeLength);
	// 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);
	+ // just a copy doesn't travel
	+ init->setLifeLength(Lorentz5Distance());
	+ init->setVertex(final->decayVertex());
	// insert shower products
	addFinalStateShower(init,pstep);
	// sort out colour
	final=_incomingLines.begin()->second;
	_colour.clear();
	mapColour(final,copy);
	}
	// 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);
	+ // just a copy does move
	+ copy->setLifeLength(Lorentz5Distance());
	+ copy->setVertex(final->decayVertex());
	// final-state radiation
	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
	for(cit=outgoingLines().begin();cit!=outgoingLines().end();++cit) {
	ShowerParticlePtr init=cit->first->progenitor();
	+ // ZERO the distance
	+ init->setLifeLength(Lorentz5Distance());
	if(!init->thePEGBase())
	throw Exception() << "Final-state particle must have a ThePEGBase"
	<< " in ShowerTree::insertDecay()"
	<< Exception::runerror;
	// if not from matrix element correction
	if(cit->first->perturbative()) {
	// add the child
	updateColour(cit->first->copy());
	PPtr orig=cit->first->original();
	+ orig->setLifeLength(Lorentz5Distance());
	+ orig->setVertex(copy->decayVertex());
	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);
	}
	+ // ZERO the distances as just copies
	+ cit->first->copy()->setLifeLength(Lorentz5Distance());
	+ init->setLifeLength(Lorentz5Distance());
	+ cit->first->copy()->setVertex(copy->decayVertex());
	+ init->setVertex(copy->decayVertex());
	// 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;
	}

	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, 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(_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(_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->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();
	}

	void ShowerTree::fixColour(tShowerParticlePtr part) {
	if(!part->colourInfo()->colourLines().empty()) {
	if(part->colourInfo()->colourLines().size()==1) {
	ColinePtr line=part->colourLine();
	if(line) {
	line->removeColoured(part);
	line->addColoured(part);
	}
	}
	else {
	Ptr<MultiColour>::pointer colour =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>(part->colourInfo());
	vector<tcColinePtr> lines = colour->colourLines();
	for(unsigned int ix=0;ix<lines.size();++ix) {
	ColinePtr line = const_ptr_cast<ColinePtr>(lines[ix]);
	if(line) {
	line->removeColoured(part);
	line->addColoured(part);
	}
	}
	}
	}
	if(!part->colourInfo()->antiColourLines().empty()) {
	if(part->colourInfo()->antiColourLines().size()==1) {
	ColinePtr line=part->antiColourLine();
	if(line) {
	line->removeAntiColoured(part);
	line->addAntiColoured(part);
	}
	}
	else {
	Ptr<MultiColour>::pointer colour =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>(part->colourInfo());
	vector<tcColinePtr> lines = colour->antiColourLines();
	for(unsigned int ix=0;ix<lines.size();++ix) {
	ColinePtr line = const_ptr_cast<ColinePtr>(lines[ix]);
	if(line) {
	line->removeAntiColoured(part);
	line->addAntiColoured(part);
	}
	}
	}
	}
	}

	vector<ShowerParticlePtr> ShowerTree::extractProgenitorParticles() {
	vector<ShowerParticlePtr> particles;
	// incoming particles
	for(map<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator
	cit=incomingLines().begin(); cit!=incomingLines().end();++cit)
	particles.push_back(cit->first->progenitor());
	// outgoing particles
	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
	cjt=outgoingLines().begin(); cjt!=outgoingLines().end();++cjt)
	particles.push_back(cjt->first->progenitor());
	return particles;
	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 12:58 PM (1 d, 6 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4022860
Default Alt Text: (115 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions