No OneTemporary
Actions

Size

239 KB

Subscribers

None

View Options

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

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

	using namespace Herwig;

	DescribeClass<ClusterFissioner,Interfaced>
	describeClusterFissioner("Herwig::ClusterFissioner","Herwig.so");

	ClusterFissioner::ClusterFissioner() :
	_clMaxLight(3.35*GeV),
	_clMaxBottom(3.35*GeV),
	_clMaxCharm(3.35*GeV),
	_clMaxExotic(3.35*GeV),
	_clPowLight(2.0),
	_clPowBottom(2.0),
	_clPowCharm(2.0),
	_clPowExotic(2.0),
	_pSplitLight(1.0),
	_pSplitBottom(1.0),
	_pSplitCharm(1.0),
	_pSplitExotic(1.0),
	_fissionPwtUquark(1),
	_fissionPwtDquark(1),
	_fissionPwtSquark(0.5),
	_fissionCluster(0),
	_btClM(1.0*GeV),
	_iopRem(1),
	_kappa(1.0e15*GeV/meter),
	_enhanceSProb(0),
	_m0Fission(2.*GeV),
	_massMeasure(0)
	{}

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

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

	void ClusterFissioner::persistentOutput(PersistentOStream & os) const {
	os << _hadronsSelector << ounit(_clMaxLight,GeV)
	<< ounit(_clMaxBottom,GeV) << ounit(_clMaxCharm,GeV)
	<< ounit(_clMaxExotic,GeV) << _clPowLight << _clPowBottom
	<< _clPowCharm << _clPowExotic << _pSplitLight
	<< _pSplitBottom << _pSplitCharm << _pSplitExotic
	<< _fissionCluster << _fissionPwtUquark << _fissionPwtDquark << _fissionPwtSquark
	<< ounit(_btClM,GeV)
	<< _iopRem << ounit(_kappa, GeV/meter)
	<< _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure;
	}

	void ClusterFissioner::persistentInput(PersistentIStream & is, int) {
	is >> _hadronsSelector >> iunit(_clMaxLight,GeV)
	>> iunit(_clMaxBottom,GeV) >> iunit(_clMaxCharm,GeV)
	>> iunit(_clMaxExotic,GeV) >> _clPowLight >> _clPowBottom
	>> _clPowCharm >> _clPowExotic >> _pSplitLight
	>> _pSplitBottom >> _pSplitCharm >> _pSplitExotic
	>> _fissionCluster >> _fissionPwtUquark >> _fissionPwtDquark >> _fissionPwtSquark
	>> iunit(_btClM,GeV) >> _iopRem
	>> iunit(_kappa, GeV/meter)
	>> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure;
	}

	void ClusterFissioner::Init() {

	static ClassDocumentation<ClusterFissioner> documentation
	("Class responsibles for chopping up the clusters");

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

	// ClMax for light, Bottom, Charm and exotic (e.g. Susy) quarks
	static Parameter<ClusterFissioner,Energy>
	interfaceClMaxLight ("ClMaxLight","cluster max mass for light quarks (unit [GeV])",
	&ClusterFissioner::_clMaxLight, GeV, 3.35GeV, ZERO, 10.0GeV,
	false,false,false);
	static Parameter<ClusterFissioner,Energy>
	interfaceClMaxBottom ("ClMaxBottom","cluster max mass for b quarks (unit [GeV])",
	&ClusterFissioner::_clMaxBottom, GeV, 3.35GeV, ZERO, 10.0GeV,
	false,false,false);
	static Parameter<ClusterFissioner,Energy>
	interfaceClMaxCharm ("ClMaxCharm","cluster max mass for c quarks (unit [GeV])",
	&ClusterFissioner::_clMaxCharm, GeV, 3.35GeV, ZERO, 10.0GeV,
	false,false,false);
	static Parameter<ClusterFissioner,Energy>
	interfaceClMaxExotic ("ClMaxExotic","cluster max mass for exotic quarks (unit [GeV])",
	&ClusterFissioner::_clMaxExotic, GeV, 3.35GeV, ZERO, 10.0GeV,
	false,false,false);

	// ClPow for light, Bottom, Charm and exotic (e.g. Susy) quarks
	static Parameter<ClusterFissioner,double>
	interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks",
	&ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false);
	static Parameter<ClusterFissioner,double>
	interfaceClPowBottom ("ClPowBottom","cluster mass exponent for b quarks",
	&ClusterFissioner::_clPowBottom, 0, 2.0, 0.0, 10.0,false,false,false);
	static Parameter<ClusterFissioner,double>
	interfaceClPowCharm ("ClPowCharm","cluster mass exponent for c quarks",
	&ClusterFissioner::_clPowCharm, 0, 2.0, 0.0, 10.0,false,false,false);
	static Parameter<ClusterFissioner,double>
	interfaceClPowExotic ("ClPowExotic","cluster mass exponent for exotic quarks",
	&ClusterFissioner::_clPowExotic, 0, 2.0, 0.0, 10.0,false,false,false);

	// PSplit for light, Bottom, Charm and exotic (e.g. Susy) quarks
	static Parameter<ClusterFissioner,double>
	interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks",
	&ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false);
	static Parameter<ClusterFissioner,double>
	interfacePSplitBottom ("PSplitBottom","cluster mass splitting param for b quarks",
	&ClusterFissioner::_pSplitBottom, 0, 1.0, 0.0, 10.0,false,false,false);
	static Parameter<ClusterFissioner,double>
	interfacePSplitCharm ("PSplitCharm","cluster mass splitting param for c quarks",
	&ClusterFissioner::_pSplitCharm, 0, 1.0, 0.0, 10.0,false,false,false);
	static Parameter<ClusterFissioner,double>
	interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks",
	&ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false);


	static Switch<ClusterFissioner,int> interfaceFission
	("Fission",
	"Option for different Fission options",
	&ClusterFissioner::_fissionCluster, 1, false, false);
	static SwitchOption interfaceFissionDefault
	(interfaceFission,
	"default",
	"Normal cluster fission which depends on the hadron selector class.",
	0);
	static SwitchOption interfaceFissionNew
	(interfaceFission,
	"new",
	"Alternative cluster fission which does not depend on the hadron selector class",
	1);


	static Parameter<ClusterFissioner,double> interfaceFissionPwtUquark
	("FissionPwtUquark",
	"Weight for fission in U quarks",
	&ClusterFissioner::_fissionPwtUquark, 1, 0.0, 1.0,
	false, false, Interface::limited);
	static Parameter<ClusterFissioner,double> interfaceFissionPwtDquark
	("FissionPwtDquark",
	"Weight for fission in D quarks",
	&ClusterFissioner::_fissionPwtDquark, 1, 0.0, 1.0,
	false, false, Interface::limited);
	static Parameter<ClusterFissioner,double> interfaceFissionPwtSquark
	("FissionPwtSquark",
	"Weight for fission in S quarks",
	&ClusterFissioner::_fissionPwtSquark, 0.5, 0.0, 1.0,
	false, false, Interface::limited);


	static Switch<ClusterFissioner,int> interfaceRemnantOption
	("RemnantOption",
	"Option for the treatment of remnant clusters",
	&ClusterFissioner::_iopRem, 1, false, false);
	static SwitchOption interfaceRemnantOptionSoft
	(interfaceRemnantOption,
	"Soft",
	"Both clusters produced in the fission of the beam cluster"
	" are treated as soft clusters.",
	0);
	static SwitchOption interfaceRemnantOptionHard
	(interfaceRemnantOption,
	"Hard",
	"Only the cluster containing the remnant is treated as a soft cluster.",
	1);
	static SwitchOption interfaceRemnantOptionVeryHard
	(interfaceRemnantOption,
	"VeryHard",
	"Even remnant clusters are treated as hard, i.e. all clusters the same",
	2);

	static Parameter<ClusterFissioner,Energy> interfaceBTCLM
	("SoftClusterFactor",
	"Parameter for the mass spectrum of remnant clusters",
	&ClusterFissioner::_btClM, GeV, 1.GeV, 0.1GeV, 10.0*GeV,
	false, false, Interface::limited);


	static Parameter<ClusterFissioner,Tension> interfaceStringTension
	("StringTension",
	"String tension used in vertex displacement calculation",
	&ClusterFissioner::_kappa, GeV/meter,
	1.0e15*GeV/meter, ZERO, ZERO,
	false, false, Interface::lowerlim);

	static Switch<ClusterFissioner,int> interfaceEnhanceSProb
	("EnhanceSProb",
	"Option for enhancing strangeness",
	&ClusterFissioner::_enhanceSProb, 0, false, false);
	static SwitchOption interfaceEnhanceSProbNo
	(interfaceEnhanceSProb,
	"No",
	"No strangeness enhancement.",
	0);
	static SwitchOption interfaceEnhanceSProbScaled
	(interfaceEnhanceSProb,
	"Scaled",
	"Scaled strangeness enhancement",
	1);
	static SwitchOption interfaceEnhanceSProbExponential
	(interfaceEnhanceSProb,
	"Exponential",
	"Exponential strangeness enhancement",
	2);

	static Switch<ClusterFissioner,int> interfaceMassMeasure
	("MassMeasure",
	"Option to use different mass measures",
	&ClusterFissioner::_massMeasure,0,false,false);
	static SwitchOption interfaceMassMeasureMass
	(interfaceMassMeasure,
	"Mass",
	"Mass Measure",
	0);
	static SwitchOption interfaceMassMeasureLambda
	(interfaceMassMeasure,
	"Lambda",
	"Lambda Measure",
	1);

	static Parameter<ClusterFissioner,Energy> interfaceFissionMassScale
	("FissionMassScale",
	"Cluster fission mass scale",
	&ClusterFissioner::_m0Fission, GeV, 2.0GeV, 0.1GeV, 50.*GeV,
	false, false, Interface::limited);

	}

	tPVector ClusterFissioner::fission(ClusterVector & clusters, bool softUEisOn) {
	// return if no clusters
	if (clusters.empty()) return tPVector();

	/*****************
	* Loop over the (input) collection of cluster pointers, and store in
	* the vector splitClusters all the clusters that need to be split
	* (these are beam clusters, if soft underlying event is off, and
	* heavy non-beam clusters).
	********************/

	stack<ClusterPtr> splitClusters;
	for(ClusterVector::iterator it = clusters.begin() ;
	it != clusters.end() ; ++it) {
	/**************
	* Skip 3-component clusters that have been redefined (as 2-component
	* clusters) or not available clusters. The latter check is indeed
	* redundant now, but it is used for possible future extensions in which,
	* for some reasons, some of the clusters found by ClusterFinder are tagged
	* straight away as not available.
	**************/
	if((it)->isRedefined() \|\| !(it)->isAvailable()) continue;
	// if the cluster is a beam cluster add it to the vector of clusters
	// to be split or if it is heavy
	if((it)->isBeamCluster() \|\| isHeavy(it)) splitClusters.push(*it);
	}
	tPVector finalhadrons;
	cut(splitClusters, clusters, finalhadrons, softUEisOn);
	return finalhadrons;
	}

	void ClusterFissioner::cut(stack<ClusterPtr> & clusterStack,
	ClusterVector &clusters, tPVector & finalhadrons,
	bool softUEisOn) {
	/**************************************************
	* This method does the splitting of the cluster pointed by cluPtr
	* and "recursively" by all of its cluster children, if heavy. All of these
	* new children clusters are added (indeed the pointers to them) to the
	* collection of cluster pointers collecCluPtr. The method works as follows.
	* Initially the vector vecCluPtr contains just the input pointer to the
	* cluster to be split. Then it will be filled "recursively" by all
	* of the cluster's children that are heavy enough to require, in their turn,
	* to be split. In each loop, the last element of the vector vecCluPtr is
	* considered (only once because it is then removed from the vector).
	* This approach is conceptually recursive, but avoid the overhead of
	* a concrete recursive function. Furthermore it requires minimal changes
	* in the case that the fission of an heavy cluster could produce more
	* than two cluster children as assumed now.
	*
	* Draw the masses: for normal, non-beam clusters a power-like mass dist
	* is used, whereas for beam clusters a fast-decreasing exponential mass
	* dist is used instead (to avoid many iterative splitting which could
	* produce an unphysical large transverse energy from a supposed soft beam
	* remnant process).
	****************************************/
	// Here we recursively loop over clusters in the stack and cut them
	while (!clusterStack.empty()) {
	// take the last element of the vector
	ClusterPtr iCluster = clusterStack.top(); clusterStack.pop();
	// split it
	cutType ct = iCluster->numComponents() == 2 ?
	cutTwo(iCluster, finalhadrons, softUEisOn) :
	cutThree(iCluster, finalhadrons, softUEisOn);

	// There are cases when we don't want to split, even if it fails mass test
	if(!ct.first.first \|\| !ct.second.first) {
	// if an unsplit beam cluster leave if for the underlying event
	if(iCluster->isBeamCluster() && softUEisOn)
	iCluster->isAvailable(false);
	continue;
	}
	// check if clusters
	ClusterPtr one = dynamic_ptr_cast<ClusterPtr>(ct.first.first);
	ClusterPtr two = dynamic_ptr_cast<ClusterPtr>(ct.second.first);
	// is a beam cluster must be split into two clusters
	if(iCluster->isBeamCluster() && (!one\|\|!two) && softUEisOn) {
	iCluster->isAvailable(false);
	continue;
	}

	// There should always be a intermediate quark(s) from the splitting
	assert(ct.first.second && ct.second.second);
	/// \todo sort out motherless quark pairs here. Watch out for 'quark in final state' errors
	iCluster->addChild(ct.first.first);
	// iCluster->addChild(ct.first.second);
	// ct.first.second->addChild(ct.first.first);

	iCluster->addChild(ct.second.first);
	// iCluster->addChild(ct.second.second);
	// ct.second.second->addChild(ct.second.first);

	// Sometimes the clusters decay C -> H + C' rather then C -> C' + C''
	if(one) {
	clusters.push_back(one);
	if(one->isBeamCluster() && softUEisOn)
	one->isAvailable(false);
	if(isHeavy(one) && one->isAvailable())
	clusterStack.push(one);
	}
	if(two) {
	clusters.push_back(two);
	if(two->isBeamCluster() && softUEisOn)
	two->isAvailable(false);
	if(isHeavy(two) && two->isAvailable())
	clusterStack.push(two);
	}
	}
	}

	namespace {

	/**
	* Check if can't make a hadron from the partons
	*/
	bool cantMakeHadron(tcPPtr p1, tcPPtr p2) {
	return ! CheckId::canBeHadron(p1->dataPtr(), p2->dataPtr());
	}

	/**
	* Check if can't make a diquark from the partons
	*/
	bool cantMakeDiQuark(tcPPtr p1, tcPPtr p2) {
	long id1 = p1->id(), id2 = p2->id();
	return ! (QuarkMatcher::Check(id1) && QuarkMatcher::Check(id2) && id1*id2>0);
	}
	}

	ClusterFissioner::cutType
	ClusterFissioner::cutTwo(ClusterPtr & cluster, tPVector & finalhadrons,
	bool softUEisOn) {
	// need to make sure only 2-cpt clusters get here
	assert(cluster->numComponents() == 2);
	tPPtr ptrQ1 = cluster->particle(0);
	tPPtr ptrQ2 = cluster->particle(1);
	Energy Mc = cluster->mass();
	assert(ptrQ1);
	assert(ptrQ2);

	// And check if those particles are from a beam remnant
	bool rem1 = cluster->isBeamRemnant(0);
	bool rem2 = cluster->isBeamRemnant(1);
	// workout which distribution to use
	bool soft1(false),soft2(false);
	switch (_iopRem) {
	case 0:
	soft1 = rem1 \|\| rem2;
	soft2 = rem2 \|\| rem1;
	break;
	case 1:
	soft1 = rem1;
	soft2 = rem2;
	break;
	}
	// Initialization for the exponential ("soft") mass distribution.
	static const int max_loop = 1000;
	int counter = 0;
	Energy Mc1 = ZERO, Mc2 = ZERO,m1=ZERO,m2=ZERO,m=ZERO;
	tcPDPtr toHadron1, toHadron2;
	PPtr newPtr1 = PPtr ();
	PPtr newPtr2 = PPtr ();
	bool succeeded = false;
	+ Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2;
	+ // pClu1(Mc1), pClu2(Mc2), pQ1(m1), pQone(m), pQtwo(m), pQ2(m2);
	do
	{
	succeeded = false;
	++counter;
	- if (_enhanceSProb == 0){
	- drawNewFlavour(newPtr1,newPtr2);
	- }
	- else {
	- Energy2 mass2 = clustermass(cluster);
	- drawNewFlavourEnhanced(newPtr1,newPtr2,mass2);
	- }
	+
	+ // get a flavour for the qqbar pair
	+ drawNewFlavour(newPtr1,newPtr2,cluster);
	+
	// check for right ordering
	assert (ptrQ2);
	assert (newPtr2);
	assert (ptrQ2->dataPtr());
	assert (newPtr2->dataPtr());
	if(cantMakeHadron(ptrQ1, newPtr1) \|\| cantMakeHadron(ptrQ2, newPtr2)) {
	swap(newPtr1, newPtr2);
	// check again
	if(cantMakeHadron(ptrQ1, newPtr1) \|\| cantMakeHadron(ptrQ2, newPtr2)) {
	throw Exception()
	<< "ClusterFissioner cannot split the cluster ("
	<< ptrQ1->PDGName() << ' ' << ptrQ2->PDGName()
	<< ") into hadrons.\n" << Exception::runerror;
	}
	}
	// Check that new clusters can produce particles and there is enough
	// phase space to choose the drawn flavour
	m1 = ptrQ1->data().constituentMass();
	m2 = ptrQ2->data().constituentMass();
	m = newPtr1->data().constituentMass();
	// Do not split in the case there is no phase space available
	if(Mc < m1+m + m2+m) continue;
	- // power for splitting
	- double exp1=_pSplitLight;
	- double exp2=_pSplitLight;

	- if (CheckId::isExotic(ptrQ1->dataPtr())) exp1 = _pSplitExotic;
	- else if(CheckId::hasBottom(ptrQ1->dataPtr()))exp1 = _pSplitBottom;
	- else if(CheckId::hasCharm(ptrQ1->dataPtr())) exp1 = _pSplitCharm;
	+ pQ1.setMass(m1);
	+ pQone.setMass(m);
	+ pQtwo.setMass(m);
	+ pQ2.setMass(m2);

	- if (CheckId::isExotic(ptrQ2->dataPtr())) exp2 = _pSplitExotic;
	- else if(CheckId::hasBottom(ptrQ2->dataPtr())) exp2 = _pSplitBottom;
	- else if(CheckId::hasCharm(ptrQ2->dataPtr())) exp2 = _pSplitCharm;
	+ pair<Energy,Energy> res = drawNewMasses(Mc, soft1, soft2, pClu1, pClu2,
	+ ptrQ1, pQ1, newPtr1, pQone,
	+ newPtr2, pQtwo, ptrQ2, pQ2);

	- // If, during the drawing of candidate masses, too many attempts fail
	- // (because the phase space available is tiny)
	- /// \todo run separate loop here?
	- Mc1 = drawChildMass(Mc,m1,m2,m,exp1,soft1);
	- Mc2 = drawChildMass(Mc,m2,m1,m,exp2,soft2);
	+ Mc1 = res.first; Mc2 = res.second;
	+
	if(Mc1 < m1+m \|\| Mc2 < m+m2 \|\| Mc1+Mc2 > Mc) continue;
	/**************************
	* New (not present in Fortran Herwig):
	* check whether the fragment masses Mc1 and Mc2 are above the
	* threshold for the production of the lightest pair of hadrons with the
	* right flavours. If not, then set by hand the mass to the lightest
	* single hadron with the right flavours, in order to solve correctly
	* the kinematics, and (later in this method) create directly such hadron
	* and add it to the children hadrons of the cluster that undergoes the
	* fission (i.e. the one pointed by iCluPtr). Notice that in this special
	* case, the heavy cluster that undergoes the fission has one single
	* cluster child and one single hadron child. We prefer this approach,
	* rather than to create a light cluster, with the mass set equal to
	* the lightest hadron, and let then the class LightClusterDecayer to do
	* the job to decay it to that single hadron, for two reasons:
	* First, because the sum of the masses of the two constituents can be,
	* in this case, greater than the mass of that hadron, hence it would
	* be impossible to solve the kinematics for such two components, and
	* therefore we would have a cluster whose components are undefined.
	* Second, the algorithm is faster, because it avoids the reshuffling
	* procedure that would be necessary if we used LightClusterDecayer
	* to decay the light cluster to the lightest hadron.
	****************************/
	+
	+ // override chosen masses if needed
	toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1);
	- if(toHadron1) Mc1 = toHadron1->mass();
	+ if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); }
	toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2);
	- if(toHadron2) Mc2 = toHadron2->mass();
	+ if(toHadron2) { Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); }
	// if a beam cluster not allowed to decay to hadrons
	if(cluster->isBeamCluster() && (toHadron1\|\|toHadron2) && softUEisOn)
	continue;
	// Check if the decay kinematics is still possible: if not then
	// force the one-hadron decay for the other cluster as well.
	if(Mc1 + Mc2 > Mc) {
	if(!toHadron1) {
	toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2);
	- if(toHadron1) Mc1 = toHadron1->mass();
	+ if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); }
	}
	else if(!toHadron2) {
	toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1);
	- if(toHadron2) Mc2 = toHadron2->mass();
	+ if(toHadron2) { Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); }
	}
	}
	succeeded = (Mc >= Mc1+Mc2);
	}
	while (!succeeded && counter < max_loop);

	if(counter >= max_loop) {
	static const PPtr null = PPtr();
	return cutType(PPair(null,null),PPair(null,null));
	}

	// Determined the (5-components) momenta (all in the LAB frame)
	Lorentz5Momentum pClu = cluster->momentum(); // known
	Lorentz5Momentum p0Q1 = ptrQ1->momentum(); // known (mom Q1 before fission)
	- Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2; //unknown
	- pClu1.setMass(Mc1);
	- pClu2.setMass(Mc2);
	- pQ1.setMass(m1);
	- pQ2.setMass(m2);
	- pQone.setMass(m);
	- pQtwo.setMass(m);
	calculateKinematics(pClu,p0Q1,toHadron1,toHadron2,
	- pClu1,pClu2,pQ1,pQone,pQtwo,pQ2); // out
	+ pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);

	/******************
	* The previous methods have determined the kinematics and positions
	* of C -> C1 + C2.
	* In the case that one of the two product is light, that means either
	* decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta
	* of the components of that light product have not been determined,
	* and a (light) cluster will not be created: the heavy father cluster
	* decays, in this case, into a single (not-light) cluster and a
	* single hadron. In the other, "normal", cases the father cluster
	* decays into two clusters, each of which has well defined components.
	* Notice that, in the case of components which point to particles, the
	* momenta of the components is properly set to the new values, whereas
	* we do not change the momenta of the pointed particles, because we
	* want to keep all of the information (that is the new momentum of a
	* component after the splitting, which is contained in the _momentum
	* member of the Component class, and the (old) momentum of that component
	* before the splitting, which is contained in the momentum of the
	* pointed particle). Please not make confusion of this only apparent
	* inconsistency!
	********************/
	LorentzPoint posC,pos1,pos2;
	posC = cluster->vertex();
	calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2);
	cutType rval;
	if(toHadron1) {
	rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1);
	finalhadrons.push_back(rval.first.first);
	}
	else {
	rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1);
	}
	if(toHadron2) {
	rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2);
	finalhadrons.push_back(rval.second.first);
	}
	else {
	rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2);
	}
	return rval;
	}


	ClusterFissioner::cutType
	ClusterFissioner::cutThree(ClusterPtr & cluster, tPVector & finalhadrons,
	bool softUEisOn) {
	// need to make sure only 3-cpt clusters get here
	assert(cluster->numComponents() == 3);
	// extract quarks
	tPPtr ptrQ[3] = {cluster->particle(0),cluster->particle(1),cluster->particle(2)};
	assert( ptrQ[0] && ptrQ[1] && ptrQ[2] );
	// find maximum mass pair
	Energy mmax(ZERO);
	Lorentz5Momentum pDiQuark;
	int iq1(-1),iq2(-1);
	Lorentz5Momentum psum;
	for(int q1=0;q1<3;++q1) {
	psum+= ptrQ[q1]->momentum();
	for(int q2=q1+1;q2<3;++q2) {
	Lorentz5Momentum ptest = ptrQ[q1]->momentum()+ptrQ[q2]->momentum();
	ptest.rescaleMass();
	Energy mass = ptest.m();
	if(mass>mmax) {
	mmax = mass;
	pDiQuark = ptest;
	iq1 = q1;
	iq2 = q2;
	}
	}
	}
	// and the spectators
	int iother(-1);
	for(int ix=0;ix<3;++ix) if(ix!=iq1&&ix!=iq2) iother=ix;
	assert(iq1>=0&&iq2>=0&&iother>=0);

	// And check if those particles are from a beam remnant
	bool rem1 = cluster->isBeamRemnant(iq1);
	bool rem2 = cluster->isBeamRemnant(iq2);
	// workout which distribution to use
	bool soft1(false),soft2(false);
	switch (_iopRem) {
	case 0:
	soft1 = rem1 \|\| rem2;
	soft2 = rem2 \|\| rem1;
	break;
	case 1:
	soft1 = rem1;
	soft2 = rem2;
	break;
	}
	// Initialization for the exponential ("soft") mass distribution.
	static const int max_loop = 1000;
	int counter = 0;
	Energy Mc1 = ZERO, Mc2 = ZERO, m1=ZERO, m2=ZERO, m=ZERO;
	tcPDPtr toHadron;
	bool toDiQuark(false);
	PPtr newPtr1 = PPtr(),newPtr2 = PPtr();
	PDPtr diquark;
	bool succeeded = false;
	+ Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2;
	do {
	succeeded = false;
	++counter;
	- if (_enhanceSProb == 0) {
	- drawNewFlavour(newPtr1,newPtr2);
	- }
	- else {
	- Energy2 mass2 = clustermass(cluster);
	- drawNewFlavourEnhanced(newPtr1,newPtr2, mass2);
	- }
	+
	+ // get a flavour for the qqbar pair
	+ drawNewFlavour(newPtr1,newPtr2,cluster);
	+
	// randomly pick which will be (anti)diquark and which a mesonic cluster
	if(UseRandom::rndbool()) {
	swap(iq1,iq2);
	swap(rem1,rem2);
	}
	// check first order
	if(cantMakeHadron(ptrQ[iq1], newPtr1) \|\| cantMakeDiQuark(ptrQ[iq2], newPtr2)) {
	swap(newPtr1,newPtr2);
	}
	// check again
	if(cantMakeHadron(ptrQ[iq1], newPtr1) \|\| cantMakeDiQuark(ptrQ[iq2], newPtr2)) {
	throw Exception()
	<< "ClusterFissioner cannot split the cluster ("
	<< ptrQ[iq1]->PDGName() << ' ' << ptrQ[iq2]->PDGName()
	<< ") into a hadron and diquark.\n" << Exception::runerror;
	}
	// Check that new clusters can produce particles and there is enough
	// phase space to choose the drawn flavour
	m1 = ptrQ[iq1]->data().constituentMass();
	m2 = ptrQ[iq2]->data().constituentMass();
	m = newPtr1->data().constituentMass();
	// Do not split in the case there is no phase space available
	if(mmax < m1+m + m2+m) continue;
	- // power for splitting
	- double exp1(_pSplitLight),exp2(_pSplitLight);
	- if (CheckId::isExotic (ptrQ[iq1]->dataPtr())) exp1 = _pSplitExotic;
	- else if(CheckId::hasBottom(ptrQ[iq1]->dataPtr())) exp1 = _pSplitBottom;
	- else if(CheckId::hasCharm (ptrQ[iq1]->dataPtr())) exp1 = _pSplitCharm;

	- if (CheckId::isExotic (ptrQ[iq2]->dataPtr())) exp2 = _pSplitExotic;
	- else if(CheckId::hasBottom(ptrQ[iq2]->dataPtr())) exp2 = _pSplitBottom;
	- else if(CheckId::hasCharm (ptrQ[iq2]->dataPtr())) exp2 = _pSplitCharm;
	+ pQ1.setMass(m1);
	+ pQone.setMass(m);
	+ pQtwo.setMass(m);
	+ pQ2.setMass(m2);

	- // If, during the drawing of candidate masses, too many attempts fail
	- // (because the phase space available is tiny)
	- /// \todo run separate loop here?
	- Mc1 = drawChildMass(mmax,m1,m2,m,exp1,soft1);
	- Mc2 = drawChildMass(mmax,m2,m1,m,exp2,soft2);
	+ pair<Energy,Energy> res = drawNewMasses(mmax, soft1, soft2, pClu1, pClu2,
	+ ptrQ[iq1], pQ1, newPtr1, pQone,
	+ newPtr2, pQtwo, ptrQ[iq1], pQ2);
	+
	+ Mc1 = res.first; Mc2 = res.second;
	+
	if(Mc1 < m1+m \|\| Mc2 < m+m2 \|\| Mc1+Mc2 > mmax) continue;
	+
	// check if need to force meson clster to hadron
	toHadron = _hadronsSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),Mc1);
	- if(toHadron) Mc1 = toHadron->mass();
	+ if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); }
	// check if need to force diquark cluster to be on-shell
	toDiQuark = false;
	diquark = CheckId::makeDiquark(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr());
	if(Mc2 < diquark->constituentMass()) {
	- Mc2 = diquark->constituentMass();
	+ Mc2 = diquark->constituentMass(); pClu2.setMass(Mc2);
	toDiQuark = true;
	}
	// if a beam cluster not allowed to decay to hadrons
	if(cluster->isBeamCluster() && toHadron && softUEisOn)
	continue;
	// Check if the decay kinematics is still possible: if not then
	// force the one-hadron decay for the other cluster as well.
	if(Mc1 + Mc2 > mmax) {
	if(!toHadron) {
	toHadron = _hadronsSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2);
	- if(toHadron) Mc1 = toHadron->mass();
	+ if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); }
	}
	else if(!toDiQuark) {
	- Mc2 = _hadronsSelector->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr());
	+ Mc2 = _hadronsSelector->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); pClu2.setMass(Mc2);
	toDiQuark = true;
	}
	}
	succeeded = (mmax >= Mc1+Mc2);
	}
	while (!succeeded && counter < max_loop);
	// check no of tries
	if(counter >= max_loop) return cutType();

	// Determine the (5-components) momenta (all in the LAB frame)
	Lorentz5Momentum p0Q1 = ptrQ[iq1]->momentum();
	- // to be determined
	- Lorentz5Momentum pClu1(Mc1), pClu2(Mc2), pQ1(m1), pQone(m), pQtwo(m), pQ2(m2);
	calculateKinematics(pDiQuark,p0Q1,toHadron,toDiQuark,
	pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);
	// positions of the new clusters
	LorentzPoint pos1,pos2;
	Lorentz5Momentum pBaryon = pClu2+ptrQ[iother]->momentum();
	calculatePositions(cluster->momentum(), cluster->vertex(), pClu1, pBaryon, pos1, pos2);
	// first the mesonic cluster/meson
	cutType rval;
	if(toHadron) {
	rval.first = produceHadron(toHadron, newPtr1, pClu1, pos1);
	finalhadrons.push_back(rval.first.first);
	}
	else {
	rval.first = produceCluster(ptrQ[iq1], newPtr1, pClu1, pos1, pQ1, pQone, rem1);
	}
	if(toDiQuark) {
	rem2 \|= cluster->isBeamRemnant(iother);
	PPtr newDiQuark = diquark->produceParticle(pClu2);
	rval.second = produceCluster(newDiQuark, ptrQ[iother], pBaryon, pos2, pClu2,
	ptrQ[iother]->momentum(), rem2);
	}
	else {
	rval.second = produceCluster(ptrQ[iq2], newPtr2, pBaryon, pos2, pQ2, pQtwo, rem2,
	ptrQ[iother],cluster->isBeamRemnant(iother));
	}
	cluster->isAvailable(false);
	return rval;
	}

	ClusterFissioner::PPair
	ClusterFissioner::produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a,
	const LorentzPoint &b) const {
	PPair rval;
	if(hadron->coloured()) {
	rval.first = (_hadronsSelector->lightestHadron(hadron,newPtr->dataPtr()))->produceParticle();
	}
	else
	rval.first = hadron->produceParticle();
	rval.second = newPtr;
	rval.first->set5Momentum(a);
	rval.first->setVertex(b);
	return rval;
	}

	ClusterFissioner::PPair ClusterFissioner::produceCluster(tPPtr ptrQ, tPPtr newPtr,
	const Lorentz5Momentum & a,
	const LorentzPoint & b,
	const Lorentz5Momentum & c,
	const Lorentz5Momentum & d,
	bool isRem,
	tPPtr spect, bool remSpect) const {
	PPair rval;
	rval.second = newPtr;
	ClusterPtr cluster = !spect ? new_ptr(Cluster(ptrQ,rval.second)) : new_ptr(Cluster(ptrQ,rval.second,spect));
	rval.first = cluster;
	cluster->set5Momentum(a);
	cluster->setVertex(b);
	assert(cluster->particle(0)->id() == ptrQ->id());
	cluster->particle(0)->set5Momentum(c);
	cluster->particle(1)->set5Momentum(d);
	cluster->setBeamRemnant(0,isRem);
	if(remSpect) cluster->setBeamRemnant(2,remSpect);
	return rval;
	}

	void ClusterFissioner::drawNewFlavour(PPtr& newPtrPos,PPtr& newPtrNeg) const {

	// Flavour is assumed to be only u, d, s, with weights
	// (which are not normalized probabilities) given
	// by the same weights as used in HadronsSelector for
	// the decay of clusters into two hadrons.


	double prob_d;
	double prob_u;
	double prob_s;
	switch(_fissionCluster){
	case 0:
	prob_d = _hadronsSelector->pwtDquark();
	prob_u = _hadronsSelector->pwtUquark();
	prob_s = _hadronsSelector->pwtSquark();
	break;
	case 1:
	prob_d = _fissionPwtDquark;
	prob_u = _fissionPwtUquark;
	prob_s = _fissionPwtSquark;
	break;
	default :
	assert(false);
	}
	int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
	long idNew = 0;
	switch (choice) {
	case 0: idNew = ThePEG::ParticleID::u; break;
	case 1: idNew = ThePEG::ParticleID::d; break;
	case 2: idNew = ThePEG::ParticleID::s; break;
	}
	newPtrPos = getParticle(idNew);
	newPtrNeg = getParticle(-idNew);
	assert (newPtrPos);
	assert(newPtrNeg);
	assert (newPtrPos->dataPtr());
	assert(newPtrNeg->dataPtr());

	}

	void ClusterFissioner::drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg,
	Energy2 mass2) const {

	// Flavour is assumed to be only u, d, s, with weights
	// (which are not normalized probabilities) given
	// by the same weights as used in HadronsSelector for
	// the decay of clusters into two hadrons.

	double prob_d;
	double prob_u;
	double prob_s = 0.;
	double scale = abs(double(sqr(_m0Fission)/mass2));
	// Choose which splitting weights you wish to use
	switch(_fissionCluster){
	// 0: ClusterFissioner and ClusterDecayer use the same weights
	case 0:
	prob_d = _hadronsSelector->pwtDquark();
	prob_u = _hadronsSelector->pwtUquark();
	/* Strangeness enhancement:
	Case 1: probability scaling
	Case 2: Exponential scaling
	*/
	if (_enhanceSProb == 1)
	prob_s = (_maxScale < scale) ? 0. : pow(_hadronsSelector->pwtSquark(),scale);
	else if (_enhanceSProb == 2)
	prob_s = (_maxScale < scale) ? 0. : exp(-scale);
	break;
	/* 1: ClusterFissioner uses its own unique set of weights,
	i.e. decoupled from ClusterDecayer */
	case 1:
	prob_d = _fissionPwtDquark;
	prob_u = _fissionPwtUquark;
	if (_enhanceSProb == 1)
	prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwtSquark,scale);
	else if (_enhanceSProb == 2)
	prob_s = (_maxScale < scale) ? 0. : exp(-scale);
	break;
	}

	int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
	long idNew = 0;
	switch (choice) {
	case 0: idNew = ThePEG::ParticleID::u; break;
	case 1: idNew = ThePEG::ParticleID::d; break;
	case 2: idNew = ThePEG::ParticleID::s; break;
	}
	newPtrPos = getParticle(idNew);
	newPtrNeg = getParticle(-idNew);
	assert (newPtrPos);
	assert(newPtrNeg);
	assert (newPtrPos->dataPtr());
	assert(newPtrNeg->dataPtr());

	}


	-Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster){
	+Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster) const{
	Lorentz5Momentum pIn = cluster->momentum();
	Energy2 endpointmass2 = sqr(cluster->particle(0)->mass() +
	cluster->particle(1)->mass());
	Energy2 singletm2 = pIn.m2();
	// Return either the cluster mass, or the lambda measure
	return (_massMeasure == 0) ? singletm2 : singletm2 - endpointmass2;
	}


	Energy ClusterFissioner::drawChildMass(const Energy M, const Energy m1,
	const Energy m2, const Energy m,
	const double expt, const bool soft) const {

	/***************************
	* This method, given in input the cluster mass Mclu of an heavy cluster C,
	* made of consituents of masses m1 and m2, draws the masses Mclu1 and Mclu2
	* of, respectively, the children cluster C1, made of constituent masses m1
	* and m, and cluster C2, of mass Mclu2 and made of constituent masses m2
	* and m. The mass is extracted from one of the two following mass
	* distributions:
	* --- power-like ("normal" distribution)
	* d(Prob) / d(M^exponent) = const
	* where the exponent can be different from the two children C1 (exp1)
	* and C2 (exponent2).
	* --- exponential ("soft" distribution)
	* d(Prob) / d(M^2) = exp(-b*M)
	* where b = 2.0 / average.
	* Such distributions are limited below by the masses of
	* the constituents quarks, and above from the mass of decaying cluster C.
	* The choice of which of the two mass distributions to use for each of the
	* two cluster children is dictated by iRemnant (see below).
	* If the number of attempts to extract a pair of mass values that are
	* kinematically acceptable is above some fixed number (max_loop, see below)
	* the method gives up and returns false; otherwise, when it succeeds, it
	* returns true.
	*
	* These distributions have been modified from HERWIG:
	* Before these were:
	* Mclu1 = m1 + (Mclu - m1 - m2)*pow( rnd(), 1.0/exponent1 );
	* The new one coded here is a more efficient version, same density
	* but taking into account 'in phase space from' beforehand
	***************************/
	// hard cluster
	if(!soft) {
	return pow(UseRandom::rnd(pow((M-m1-m2-m)*UnitRemoval::InvE, expt),
	pow(m*UnitRemoval::InvE, expt)), 1./expt
	)*UnitRemoval::E + m1;
	}
	// Otherwise it uses a soft mass distribution
	else {
	static const InvEnergy b = 2.0 / _btClM;
	Energy max = M-m1-m2-2.0*m;
	double rmin = b*max;
	rmin = ( rmin < 50 ) ? exp(-rmin) : 0.;
	double r1;
	do {
	r1 = UseRandom::rnd(rmin, 1.0) * UseRandom::rnd(rmin, 1.0);
	}
	while (r1 < rmin);
	return m1 + m - log(r1)/b;
	}
	}


	void ClusterFissioner::calculateKinematics(const Lorentz5Momentum & pClu,
	const Lorentz5Momentum & p0Q1,
	const bool toHadron1,
	const bool toHadron2,
	Lorentz5Momentum & pClu1,
	Lorentz5Momentum & pClu2,
	Lorentz5Momentum & pQ1,
	Lorentz5Momentum & pQbar,
	Lorentz5Momentum & pQ,
	Lorentz5Momentum & pQ2bar) const {

	/******************
	* This method solves the kinematics of the two body cluster decay:
	* C (Q1 Q2bar) ---> C1 (Q1 Qbar) + C2 (Q Q2bar)
	* In input we receive the momentum of C, pClu, and the momentum
	* of the quark Q1 (constituent of C), p0Q1, both in the LAB frame.
	* Furthermore, two boolean variables inform whether the two fission
	* products (C1, C2) decay immediately into a single hadron (in which
	* case the cluster itself is identify with that hadron) and we do
	* not have to solve the kinematics of the components (Q1,Qbar) for
	* C1 and (Q,Q2bar) for C2.
	* The output is given by the following momenta (all 5-components,
	* and all in the LAB frame):
	* pClu1 , pClu2 respectively of C1 , C2
	* pQ1 , pQbar respectively of Q1 , Qbar in C1
	* pQ , pQ2bar respectively of Q , Q2 in C2
	* The assumption, suggested from the string model, is that, in C frame,
	* C1 and its constituents Q1 and Qbar are collinear, and collinear to
	* the direction of Q1 in C (that is before cluster decay); similarly,
	* (always in the C frame) C2 and its constituents Q and Q2bar are
	* collinear (and therefore anti-collinear with C1,Q1,Qbar).
	* The solution is then obtained by using Lorentz boosts, as follows.
	* The kinematics of C1 and C2 is solved in their parent C frame,
	* and then boosted back in the LAB. The kinematics of Q1 and Qbar
	* is solved in their parent C1 frame and then boosted back in the LAB;
	* similarly, the kinematics of Q and Q2bar is solved in their parent
	* C2 frame and then boosted back in the LAB. In each of the three
	* "two-body decay"-like cases, we use the fact that the direction
	* of the motion of the decay products is known in the rest frame of
	* their parent. This is obvious for the first case in which the
	* parent rest frame is C; but it is also true in the other two cases
	* where the rest frames are C1 and C2. This is because C1 and C2
	* are boosted w.r.t. C in the same direction where their components,
	* respectively (Q1,Qbar) and (Q,Q2bar) move in C1 and C2 rest frame
	* respectively.
	* Of course, although the notation used assumed that C = (Q1 Q2bar)
	* where Q1 is a quark and Q2bar an antiquark, indeed everything remain
	* unchanged also in all following cases:
	* Q1 quark, Q2bar antiquark; --> Q quark;
	* Q1 antiquark , Q2bar quark; --> Q antiquark;
	* Q1 quark, Q2bar diquark; --> Q quark
	* Q1 antiquark, Q2bar anti-diquark; --> Q antiquark
	* Q1 diquark, Q2bar quark --> Q antiquark
	* Q1 anti-diquark, Q2bar antiquark; --> Q quark
	**************************/

	// Calculate the unit three-vector, in the C frame, along which
	// all of the constituents and children clusters move.
	Lorentz5Momentum u(p0Q1);
	u.boost( -pClu.boostVector() ); // boost from LAB to C
	// the unit three-vector is then u.vect().unit()

	// Calculate the momenta of C1 and C2 in the (parent) C frame first,
	// where the direction of C1 is u.vect().unit(), and then boost back in the
	// LAB frame.

	if (pClu.m() < pClu1.mass() + pClu2.mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (A)"
	<< Exception::eventerror;
	}
	Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(),
	u.vect().unit(), pClu1, pClu2);

	// In the case that cluster1 does not decay immediately into a single hadron,
	// calculate the momenta of Q1 (as constituent of C1) and Qbar in the
	// (parent) C1 frame first, where the direction of Q1 is u.vect().unit(),
	// and then boost back in the LAB frame.
	if(!toHadron1) {
	if (pClu1.m() < pQ1.mass() + pQbar.mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (B)"
	<< Exception::eventerror;
	}
	Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(),
	u.vect().unit(), pQ1, pQbar);
	}

	// In the case that cluster2 does not decay immediately into a single hadron,
	// Calculate the momenta of Q and Q2bar (as constituent of C2) in the
	// (parent) C2 frame first, where the direction of Q is u.vect().unit(),
	// and then boost back in the LAB frame.
	if(!toHadron2) {
	if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (C)"
	<< Exception::eventerror;
	}
	Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(),
	u.vect().unit(), pQ, pQ2bar);
	}
	}


	void ClusterFissioner::calculatePositions(const Lorentz5Momentum & pClu,
	const LorentzPoint & positionClu,
	const Lorentz5Momentum & pClu1,
	const Lorentz5Momentum & pClu2,
	LorentzPoint & positionClu1,
	LorentzPoint & positionClu2) const {
	// Determine positions of cluster children.
	// See Marc Smith's thesis, page 127, formulas (4.122) and (4.123).
	Energy Mclu = pClu.m();
	Energy Mclu1 = pClu1.m();
	Energy Mclu2 = pClu2.m();

	// Calculate the unit three-vector, in the C frame, along which
	// children clusters move.
	Lorentz5Momentum u(pClu1);
	u.boost( -pClu.boostVector() ); // boost from LAB to C frame
	// the unit three-vector is then u.vect().unit()

	Energy pstarChild = Kinematics::pstarTwoBodyDecay(Mclu,Mclu1,Mclu2);

	// First, determine the relative positions of the children clusters
	// in the parent cluster reference frame.
	Length x1 = ( 0.25Mclu + 0.5( pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
	Length t1 = Mclu/_kappa - x1;
	LorentzDistance distanceClu1( x1 * u.vect().unit(), t1 );

	Length x2 = (-0.25Mclu + 0.5(-pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
	Length t2 = Mclu/_kappa + x2;
	LorentzDistance distanceClu2( x2 * u.vect().unit(), t2 );

	// Then, transform such relative positions from the parent cluster
	// reference frame to the Lab frame.
	distanceClu1.boost( pClu.boostVector() );
	distanceClu2.boost( pClu.boostVector() );

	// Finally, determine the absolute positions in the Lab frame.
	positionClu1 = positionClu + distanceClu1;
	positionClu2 = positionClu + distanceClu2;

	}

	bool ClusterFissioner::isHeavy(tcClusterPtr clu) {
	// default
	double clpow = _clPowLight;
	Energy clmax = _clMaxLight;
	// particle data for constituents
	tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()};
	for(int ix=0;ix<min(clu->numComponents(),3);++ix) {
	cptr[ix]=clu->particle(ix)->dataPtr();
	}
	// different parameters for exotic, bottom and charm clusters
	if(CheckId::isExotic(cptr[0],cptr[1],cptr[1])) {
	clpow = _clPowExotic;
	clmax = _clMaxExotic;
	}
	else if(CheckId::hasBottom(cptr[0],cptr[1],cptr[1])) {
	clpow = _clPowBottom;
	clmax = _clMaxBottom;
	}
	else if(CheckId::hasCharm(cptr[0],cptr[1],cptr[1])) {
	clpow = _clPowCharm;
	clmax = _clMaxCharm;
	}
	bool aboveCutoff = (
	pow(clu->mass()*UnitRemoval::InvE , clpow)
	>
	pow(clmax*UnitRemoval::InvE, clpow)
	+ pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow)
	);
	// required test for SUSY clusters, since aboveCutoff alone
	// cannot guarantee (Mc > m1 + m2 + 2*m) in cut()
	static const Energy minmass
	= getParticleData(ParticleID::d)->constituentMass();
	bool canSplitMinimally
	= clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass;
	return aboveCutoff && canSplitMinimally;
	}
	diff --git a/Hadronization/ClusterFissioner.h b/Hadronization/ClusterFissioner.h
	--- a/Hadronization/ClusterFissioner.h
	+++ b/Hadronization/ClusterFissioner.h
	@@ -1,439 +1,490 @@
	// -- C++ --
	//
	// ClusterFissioner.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_ClusterFissioner_H
	#define HERWIG_ClusterFissioner_H

	#include <ThePEG/Interface/Interfaced.h>
	#include "CluHadConfig.h"
	#include "HadronSelector.h"
	#include "ClusterFissioner.fh"
	+#include "CheckId.h"

	namespace Herwig {
	using namespace ThePEG;

	//class Cluster; // forward declaration

	/** \ingroup Hadronization
	* \class ClusterFissioner
	* \brief This class handles clusters which are too heavy.
	* \author Philip Stephens
	* \author Alberto Ribon
	* \author Stefan Gieseke
	*
	* This class does the job of chopping up either heavy clusters or beam
	* clusters in two lighter ones. The procedure is repeated recursively until
	* all of the cluster children have masses below some threshold values.
	*
	* For the beam remnant clusters, at the moment what is done is the following.
	* In the case that the soft underlying event is switched on, the
	* beam remnant clusters are tagged as not available,
	* therefore they will not be treated at all during the hadronization.
	* In the case instead that the soft underlying event is switched off,
	* then the beam remnant clusters are treated exactly as "normal" clusters,
	* with the only exception of the mass spectrum used to generate the
	* cluster children masses. For non-beam clusters, the masses of the cluster
	* children are draw from a power-like mass distribution; for beam clusters,
	* according to the value of the flag _IOpRem, either both
	* children masses are draw from a fast-decreasing exponential mass
	* distribution (case _IOpRem == 0, or, indendently by
	* _IOpRem, in the special case that the beam cluster contains two
	* beam remnants), or one mass from the exponential distribution (corresponding
	* of the cluster child with the beam remnant) and the other with the usual
	* power-like distribution (case _IOpRem == 1, which is the
	* default one, as in Herwig 6.3).
	*
	* The reason behind the use of a fast-decreasing exponential distribution
	* is that to avoid a large transverse energy from the many sequential
	* fissions that would otherwise occur due to the typical large cluster
	* mass of beam clusters. Using instead an exponential distribution
	* the masses of the two cluster children will be very small (order of
	* GeV).
	*
	* The rationale behind the implementation of the splitting of clusters
	* has been to preserve all of the information about such splitting
	* process. More explicitly a ThePEG::Step class is passed in and the
	* new clusters are added to the step as the decay products of the
	* heavy cluster. This approach has the twofold
	* advantage to provide all of the information that could be needed
	* (expecially in future developments), without any information loss,
	* and furthermore it allows a better debugging.
	*
	* @see HadronSelector
	* @see \ref ClusterFissionerInterfaces "The interfaces"
	* defined for ClusterFissioner.
	*/
	class ClusterFissioner: public Interfaced {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Default constructor.
	*/
	ClusterFissioner();

	//@}

	/** Splits the clusters which are too heavy.
	*
	* Split either heavy clusters or beam clusters recursively until all
	* children have mass below some threshold. Heavy clusters are those that
	* satisfy the condition
	* \f[ M^P > C^P + S^P \f]
	* where \f$ M \f$ is the clusters mass, \f$ P \f$ is the parameter
	* ClPow, \f$ C \f$ is the parameter ClMax and \f$ S \f$ is the
	* sum of the clusters constituent partons.
	* For beam clusters, they are split only if the soft underlying event
	* is switched off, otherwise these clusters will be tagged as unavailable
	* and they will not be treated by the hadronization altogether.
	* In the case beam clusters will be split, the procedure is exactly
	* the same as for normal non-beam clusters, with the only exception
	* of the mass spectrum from which to draw the masses of the two
	* cluster children (see method drawChildrenMasses for details).
	*/
	tPVector fission(ClusterVector & clusters, bool softUEisOn);

	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.
	*/
	ClusterFissioner & operator=(const ClusterFissioner &) = delete;

	/**
	* This method directs the splitting of the heavy clusters
	*
	* This method does the splitting of the clusters and all of its cluster
	* children, if heavy. All of these new children clusters are added to the
	* collection of clusters. The method works as follows.
	* Initially the vector contains just the stack of input pointers to the
	* clusters to be split. Then it will be filled recursively by all
	* of the cluster's children that are heavy enough to require
	* to be split. In each loop, the last element of the vector is
	* considered (only once because it is then removed from the vector).
	*
	* \todo is the following still true?
	* For normal, non-beam clusters, a power-like mass distribution
	* is used, whereas for beam clusters a fast-decreasing exponential mass
	* distribution is used instead. This avoids many iterative splitting which
	* could produce an unphysical large transverse energy from a supposed
	* soft beam remnant process.
	*/
	void cut(stack<ClusterPtr> &,
	ClusterVector&, tPVector & finalhadrons, bool softUEisOn);

	public:

	/**
	* Definition for easy passing of two particles.
	*/
	typedef pair<PPtr,PPtr> PPair;

	/**
	* Definition for use in the cut function.
	*/
	typedef pair<PPair,PPair> cutType;

	/**
	* Splits the input cluster.
	*
	* Split the input cluster (which can be either an heavy non-beam
	* cluster or a beam cluster). The result is two pairs of particles. The
	* first element of each pair is new cluster/hadron, while the second
	* element of each pair is the particle drawn from the vacuum to create
	* the new cluster/hadron.
	* Notice that this method treats also beam clusters by using a different
	* mass spectrum used to generate the cluster child masses (see method
	* drawChildMass).
	*/
	//@{
	/**
	* Split two-component cluster
	*/
	virtual cutType cutTwo(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);

	/**
	* Split three-component cluster
	*/
	virtual cutType cutThree(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);
	//@}
	public:

	/**
	* Produces a hadron and returns the flavour drawn from the vacuum.
	*
	* This routine produces a new hadron. It
	* also sets the momentum and vertex to the values given.
	*/
	PPair produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a,
	const LorentzPoint &b) const;
	+
	+protected:
	+
	+ /**
	+ * Function that returns either the cluster mass or the lambda measure
	+ */
	+ Energy2 clustermass(const ClusterPtr & cluster) const;
	+
	+ /**
	+ * Draw a new flavour for the given cluster; currently defaults to
	+ * the default model
	+ */
	+ virtual void drawNewFlavour(PPtr& newPtr1, PPtr& newPtr2, const ClusterPtr & cluster) const {
	+ if (_enhanceSProb == 0){
	+ drawNewFlavour(newPtr1,newPtr2);
	+ }
	+ else {
	+ drawNewFlavourEnhanced(newPtr1,newPtr2,clustermass(cluster));
	+ }
	+ }
	+
	+ /**
	+ * Calculate the masses and possibly kinematics of the cluster
	+ * fission at hand; if claculateKineamtics is perfomring non-trivial
	+ * steps kinematics claulcated here will be overriden. Currentl;y resorts to the default
	+ */
	+ virtual pair<Energy,Energy> drawNewMasses(Energy Mc, bool soft1, bool soft2,
	+ Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
	+ tPPtr ptrQ1, Lorentz5Momentum& pQ1,
	+ tPPtr, Lorentz5Momentum& pQone,
	+ tPPtr, Lorentz5Momentum& pQtwo,
	+ tPPtr ptrQ2, Lorentz5Momentum& pQ2) const {
	+
	+ pair<Energy,Energy> result;
	+
	+ double exp1=_pSplitLight;
	+ double exp2=_pSplitLight;
	+
	+ if (CheckId::isExotic(ptrQ1->dataPtr())) exp1 = _pSplitExotic;
	+ else if(CheckId::hasBottom(ptrQ1->dataPtr()))exp1 = _pSplitBottom;
	+ else if(CheckId::hasCharm(ptrQ1->dataPtr())) exp1 = _pSplitCharm;
	+
	+ if (CheckId::isExotic(ptrQ2->dataPtr())) exp2 = _pSplitExotic;
	+ else if(CheckId::hasBottom(ptrQ2->dataPtr())) exp2 = _pSplitBottom;
	+ else if(CheckId::hasCharm(ptrQ2->dataPtr())) exp2 = _pSplitCharm;
	+
	+ result.first = drawChildMass(Mc,pQ1.mass(),pQ2.mass(),pQone.mass(),exp1,soft1);
	+ result.second = drawChildMass(Mc,pQ2.mass(),pQ1.mass(),pQtwo.mass(),exp2,soft2);
	+
	+ pClu1.setMass(result.first);
	+ pClu2.setMass(result.second);
	+
	+ return result;
	+
	+ }
	+
	+ /**
	+ * Calculate the final kinematics of a heavy cluster decay C->C1 +
	+ * C2, if not already performed by drawNewMasses
	+ */
	+ virtual void calculateKinematics(const Lorentz5Momentum &pClu,
	+ const Lorentz5Momentum &p0Q1,
	+ const bool toHadron1, const bool toHadron2,
	+ Lorentz5Momentum &pClu1, Lorentz5Momentum &pClu2,
	+ Lorentz5Momentum &pQ1, Lorentz5Momentum &pQb,
	+ Lorentz5Momentum &pQ2, Lorentz5Momentum &pQ2b) const;
	+
	protected:

	/**
	* Produces a cluster from the flavours passed in.
	*
	* This routine produces a new cluster with the flavours given by ptrQ and newPtr.
	* The new 5 momentum is a and the parent momentum are c and d. C is for the
	* ptrQ and d is for the new particle newPtr. rem specifies whether the existing
	* particle is a beam remnant or not.
	*/
	PPair produceCluster(tPPtr ptrQ, tPPtr newPtr, const Lorentz5Momentum &a,
	const LorentzPoint &b, const Lorentz5Momentum &c,
	const Lorentz5Momentum &d, const bool rem,
	tPPtr spect=tPPtr(), bool remSpect=false) const;

	/**
	* Returns the new quark-antiquark pair
	* needed for fission of a heavy cluster. Equal probabilities
	* are assumed for producing u, d, or s pairs.
	*/
	- void drawNewFlavour(PPtr& newPtrPos,PPtr& newPtrNeg) const;
	+ void drawNewFlavour(PPtr& newPtrPos, PPtr& newPtrNeg) const;

	/**
	* Returns the new quark-antiquark pair
	* needed for fission of a heavy cluster. Equal probabilities
	* are assumed for producing u, d, or s pairs.
	* Extra argument is used when performing strangeness enhancement
	*/
	void drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, Energy2 mass2) const;


	/**
	* Produces the mass of a child cluster.
	*
	* Draw the masses \f$M'\f$ of the the cluster child produced
	* by the fission of an heavy cluster (of mass M). m1, m2 are the masses
	* of the constituents of the cluster; m is the mass of the quark extract
	* from the vacuum (together with its antiparticle). The algorithm produces
	* the mass of the cluster formed with consituent m1.
	* Two mass distributions can be used for the child cluster mass:
	* -# power-like mass distribution ("normal" mass) with power exp
	* \f[ M' = {\rm rnd}((M-m_1-m_2-m)^P, m^p)^{1/P} + m_1 \f]
	* where \f$ P \f$ is a parameter of the model and \f$ \rm{rnd} \f$ is
	* the function:
	* \f[ \rm{rnd}(a,b) = (1-r)a + r b \f]
	* and here \f$ r \f$ is a random number [0,1].
	* -# fast-decreasing exponential mass distribution ("soft" mass) with
	* rmin. rmin is given by
	* \f[ r_{\rm min} = \exp(-b (M - m_1 - m_2 - 2 m)) \f]
	* where \f$ b \f$ is a parameter of the model. The generated mass is
	* given by
	* \f[ M' = m_1 + m - \frac{\log\left(
	* {\rm rnd}(r_{\rm min}, 1-r_{\rm min})\right)}{b} \f].
	*
	* The choice of which mass distribution should be used for each of the two
	* cluster children is dictated by the parameter soft.
	*/
	Energy drawChildMass(const Energy M, const Energy m1, const Energy m2,
	const Energy m, const double exp, const bool soft) const;

	/**
	- * Determines the kinematics of a heavy cluster decay C->C1 + C2
	- */
	- void calculateKinematics(const Lorentz5Momentum &pClu,
	- const Lorentz5Momentum &p0Q1,
	- const bool toHadron1, const bool toHadron2,
	- Lorentz5Momentum &pClu1, Lorentz5Momentum &pClu2,
	- Lorentz5Momentum &pQ1, Lorentz5Momentum &pQb,
	- Lorentz5Momentum &pQ2, Lorentz5Momentum &pQ2b) const;
	-
	- /**
	* Determine the positions of the two children clusters.
	*
	* This routine generates the momentum of the decay products. It also
	* generates the momentum in the lab frame of the partons drawn out of
	* the vacuum.
	*/
	void calculatePositions(const Lorentz5Momentum &pClu,
	const LorentzPoint & positionClu,
	const Lorentz5Momentum & pClu1,
	const Lorentz5Momentum & pClu2,
	LorentzPoint & positionClu1,
	LorentzPoint & positionClu2 ) const;

	-
	-
	- /**
	- * Function that returns either the cluster mass or the lambda measure
	- */
	- Energy2 clustermass(const ClusterPtr & cluster);
	-
	protected:

	/** @name Access members for child classes. */
	//@{
	/**
	* Access to the hadron selector
	*/
	HadronSelectorPtr hadronsSelector() const {return _hadronsSelector;}

	/**
	* Access to soft-cluster parameter
	*/
	Energy btClM() const {return _btClM;}

	/**
	* Cluster splitting paramater for light quarks
	*/
	double pSplitLight() const {return _pSplitLight;}

	/**
	* Cluster splitting paramater for bottom quarks
	*/
	double pSplitBottom() const {return _pSplitBottom;}

	/**
	* Cluster splitting paramater for charm quarks
	*/
	double pSplitCharm() const {return _pSplitCharm;}

	/**
	* Cluster splitting paramater for exotic particles
	*/
	double pSplitExotic() const {return _pSplitExotic;}
	//@}

	private:

	/**
	* Check if a cluster is heavy enough to split again
	*/
	bool isHeavy(tcClusterPtr );

	/**
	* A pointer to a Herwig::HadronSelector object for generating hadrons.
	*/
	HadronSelectorPtr _hadronsSelector;

	/**
	* @name The Cluster max mass,dependant on which quarks are involved, used to determine when
	* fission will occur.
	*/
	//@{
	Energy _clMaxLight;
	Energy _clMaxBottom;
	Energy _clMaxCharm;
	Energy _clMaxExotic;
	//@}
	/**
	* @name The power used to determine when cluster fission will occur.
	*/
	//@{
	double _clPowLight;
	double _clPowBottom;
	double _clPowCharm;
	double _clPowExotic;
	//@}
	/**
	* @name The power, dependant on whic quarks are involved, used in the cluster mass generation.
	*/
	//@{
	double _pSplitLight;
	double _pSplitBottom;
	double _pSplitCharm;
	double _pSplitExotic;


	// weights for alternaive cluster fission
	double _fissionPwtUquark;
	double _fissionPwtDquark;
	double _fissionPwtSquark;

	/**
	* Flag used to determine between normal cluster fission and alternative cluster fission
	*/
	int _fissionCluster;

	//@}
	/**
	* Parameter used (2/b) for the beam cluster mass generation.
	* Currently hard coded value.
	*/
	Energy _btClM;

	/**
	* Flag used to determine what distributions to use for the cluster masses.
	*/
	int _iopRem;

	/**
	* The string constant
	*/
	Tension _kappa;

	/**
	* Flag that switches between no strangeness enhancement, scaling enhancement,
	* and exponential enhancement (in numerical order)
	*/
	int _enhanceSProb;

	/**
	* Parameter that governs the strangeness enhancement scaling
	*/
	Energy _m0Fission;

	/**
	* Flag that switches between mass measures used in strangeness enhancement:
	* cluster mass, or the lambda measure - ( m_{clu}^2 - (m_q + m_{qbar})^2 )
	*/
	int _massMeasure;

	/**
	* Constant variable which stops the scale from being to large, and not worth
	* calculating
	*/
	const double _maxScale = 20.;


	};

	}

	#endif /* HERWIG_ClusterFissioner_H */
	diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc
	--- a/Hadronization/ClusterHadronizationHandler.cc
	+++ b/Hadronization/ClusterHadronizationHandler.cc
	@@ -1,320 +1,463 @@
	// -- C++ --
	//
	// ClusterHadronizationHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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/Switch.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","Herwig.so");

	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
	+ << reshuffle_ << reshuffleMode_ << gluonMassGenerator_
	<< ounit(_minVirtuality2,GeV2) << ounit(_maxDisplacement,mm)
	<< _underlyingEventHandler << _reduceToTwoComponents;
	}


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


	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 Reference<ClusterHadronizationHandler,GluonMassGenerator> interfaceGluonMassGenerator
	+ ("GluonMassGenerator",
	+ "Set a reference to a gluon mass generator.",
	+ &ClusterHadronizationHandler::gluonMassGenerator_, false, false, true, true, false);
	+
	+ static Switch<ClusterHadronizationHandler,bool> interfaceReshuffle
	+ ("Reshuffle",
	+ "Perform reshuffling if constituent masses have not yet been included by the shower",
	+ &ClusterHadronizationHandler::reshuffle_, false, false, false);
	+ static SwitchOption interfaceReshuffleYes
	+ (interfaceReshuffle,
	+ "Global",
	+ "Do reshuffle.",
	+ true);
	+ static SwitchOption interfaceReshuffleNo
	+ (interfaceReshuffle,
	+ "No",
	+ "Do not reshuffle.",
	+ false);
	+
	+ static Switch<ClusterHadronizationHandler,int> interfaceReshuffleMode
	+ ("ReshuffleMode",
	+ "Which mode is used for the reshuffling to constituent masses",
	+ &ClusterHadronizationHandler::reshuffleMode_, 0, false, false);
	+ static SwitchOption interfaceReshuffleModeGlobal
	+ (interfaceReshuffleMode,
	+ "Global",
	+ "Global reshuffling on all final state partons",
	+ 0);
	+ static SwitchOption interfaceReshuffleModeColourConnected
	+ (interfaceReshuffleMode,
	+ "ColourConnected",
	+ "Separate reshuffling for colour connected partons",
	+ 1);
	+
	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);

	static Switch<ClusterHadronizationHandler,bool> interfaceReduceToTwoComponents
	("ReduceToTwoComponents",
	"Whether or not to reduce three component baryon-number violating clusters to two components before cluster splitting or leave"
	" this till after the cluster splitting",
	&ClusterHadronizationHandler::_reduceToTwoComponents, true, false, false);
	static SwitchOption interfaceReduceToTwoComponentsYes
	(interfaceReduceToTwoComponents,
	"BeforeSplitting",
	"Reduce to two components",
	true);
	static SwitchOption interfaceReduceToTwoComponentsNo
	(interfaceReduceToTwoComponents,
	"AfterSplitting",
	"Treat as three components",
	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());
	+
	+ PVector theList(tagged.begin(),tagged.end());
	+
	+ if ( reshuffle_ ) {
	+
	+ vector<PVector> reshufflelists;
	+
	+ if (reshuffleMode_==0){ // global reshuffling
	+ reshufflelists.push_back(theList);
	+ }
	+ else if (reshuffleMode_==1){// colour connected reshuffling
	+ splitIntoColourSinglets(theList, reshufflelists);
	+ }
	+
	+ for (auto currentlist : reshufflelists){
	+ // get available energy and energy needed for constituent mass shells
	+ LorentzMomentum totalQ;
	+ Energy needQ = ZERO;
	+ size_t nGluons = 0; // number of gluons for which a mass need be generated
	+ for ( auto p : currentlist ) {
	+ totalQ += p->momentum();
	+ if ( p->id() == ParticleID::g && gluonMassGenerator() ) {
	+ ++nGluons;
	+ continue;
	+ }
	+ needQ += p->dataPtr()->constituentMass();
	+ }
	+ Energy Q = totalQ.m();
	+ if ( needQ > Q )
	+ throw Exception() << "cannot reshuffle to constituent mass shells" << Exception::eventerror;
	+
	+ // generate gluon masses if needed
	+ list<Energy> gluonMasses;
	+ if ( nGluons && gluonMassGenerator() )
	+ gluonMasses = gluonMassGenerator()->generateMany(nGluons,Q-needQ);
	+
	+ // set masses for inidividual particles
	+ vector<Energy> masses;
	+ for ( auto p : currentlist ) {
	+ if ( p->id() == ParticleID::g && gluonMassGenerator() ) {
	+ list<Energy>::const_iterator it = gluonMasses.begin();
	+ advance(it,UseRandom::irnd(gluonMasses.size()));
	+ masses.push_back(*it);
	+ gluonMasses.erase(it);
	+ }
	+ else {
	+ masses.push_back(p->dataPtr()->constituentMass());
	+ }
	+ }
	+
	+ // reshuffle to new masses
	+ reshuffle(currentlist,masses);
	+
	+ }
	+
	+ }
	+
	// 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);
	+ for(unsigned int ix=0;ix<theList.size();++ix) {
	+ if(theList[ix]->scale()<Q02) theList[ix]->scale(Q02);
	}

	// split the gluons
	- _partonSplitter->split(currentlist);
	+ _partonSplitter->split(theList);

	// form the clusters
	ClusterVector clusters =
	- _clusterFinder->formClusters(currentlist);
	+ _clusterFinder->formClusters(theList);
	// reduce BV clusters to two components now if needed
	if(_reduceToTwoComponents)
	_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);


	// forms diquarks
	_clusterFinder->reduceToTwoComponents(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());


	// if clusters not previously reduced to two components do it now
	if(!_reduceToTwoComponents)
	_clusterFinder->reduceToTwoComponents(clusters);

	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(),
	std::mem_fn(&Particle::undecay));
	}
	}
	if (!lightOK) {
	throw Exception( "CluHad::handle(): tried LightClusterDecayer 10 times!",
	Exception::eventerror);
	}

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

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

	bool finalStateCluster=false;
	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()){
	// If there is a cluster in the final state throw an event error
	if((*it)->id()==81) {
	finalStateCluster=true;
	}
	pstep->addDecayProduct(*it);
	}
	else {
	pstep->addDecayProduct(*it);
	pstep->addIntermediate(*it);
	}
	}

	// For very small center of mass energies it might happen that baryonic clusters cannot decay into hadrons
	if (finalStateCluster){
	throw Exception( "CluHad::Handle(): Cluster in the final state",
	Exception::eventerror);
	}
	// *****************************************
	// *****************************************
	// *****************************************

	// soft underlying event if needed
	if (isSoftUnderlyingEventON()) {
	assert(_underlyingEventHandler);
	ch.performStep(_underlyingEventHandler,Hint::Default());
	}
	}


	// Sets parent child relationship of all clusters with two components
	// Relationships for clusters with more than two components are set elsewhere in the Colour Reconnector
	void ClusterHadronizationHandler::_setChildren(const ClusterVector & clusters) const {
	// erase existing information about the partons' children
	tPVector partons;
	for ( const auto & cl : clusters ) {
	if ( cl->numComponents() > 2 ) continue;
	partons.push_back( cl->colParticle() );
	partons.push_back( cl->antiColParticle() );
	}
	// erase all previous information about parent child relationship
	for_each(partons.begin(), partons.end(), std::mem_fn(&Particle::undecay));

	// give new parents to the clusters: their constituents
	for ( const auto & cl : clusters ) {
	if ( cl->numComponents() > 2 ) continue;
	cl->colParticle()->addChild(cl);
	cl->antiColParticle()->addChild(cl);
	}
	}
	+
	+void ClusterHadronizationHandler::splitIntoColourSinglets(PVector copylist,
	+ vector<PVector>& reshufflelists){
	+
	+ PVector currentlist;
	+ bool gluonloop;
	+ PPtr firstparticle, temp;
	+ reshufflelists.clear();
	+
	+ while (copylist.size()>0){
	+ gluonloop=false;
	+ currentlist.clear();
	+
	+ firstparticle=copylist.back();
	+ copylist.pop_back();
	+
	+ if (!firstparticle->coloured()){
	+ continue; //non-coloured particles are not included
	+ }
	+
	+ currentlist.push_back(firstparticle);
	+
	+ //go up the anitColourLine and check if we are in a gluon loop
	+ temp=firstparticle;
	+ while( temp->hasAntiColour()){
	+ temp = temp->antiColourLine()->endParticle();
	+ if(temp==firstparticle){
	+ gluonloop=true;
	+ break;
	+ }
	+ else{
	+ currentlist.push_back(temp);
	+ copylist.erase(remove(copylist.begin(),copylist.end(), temp), copylist.end());
	+ }
	+ }
	+
	+ //if not a gluon loop, go up the ColourLine
	+ if(!gluonloop){
	+ temp=firstparticle;
	+ while( temp->hasColour()){
	+ temp=temp->colourLine()->startParticle();
	+ currentlist.push_back(temp);
	+ copylist.erase(remove(copylist.begin(),copylist.end(), temp), copylist.end());
	+ }
	+ }
	+
	+ reshufflelists.push_back(currentlist);
	+ }
	+
	+}
	diff --git a/Hadronization/ClusterHadronizationHandler.h b/Hadronization/ClusterHadronizationHandler.h
	--- a/Hadronization/ClusterHadronizationHandler.h
	+++ b/Hadronization/ClusterHadronizationHandler.h
	@@ -1,214 +1,247 @@
	// -- C++ --
	//
	// ClusterHadronizationHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_ClusterHadronizationHandler_H
	#define HERWIG_ClusterHadronizationHandler_H

	#include <ThePEG/Handlers/HadronizationHandler.h>
	#include "PartonSplitter.h"
	#include "ClusterFinder.h"
	#include "ColourReconnector.h"
	#include "ClusterFissioner.h"
	#include "LightClusterDecayer.h"
	#include "ClusterDecayer.h"
	#include "ClusterHadronizationHandler.fh"
	+#include "Herwig/Utilities/Reshuffler.h"
	+#include "GluonMassGenerator.h"

	namespace Herwig {
	using namespace ThePEG;


	/** \ingroup Hadronization
	* \class ClusterHadronizationHandler
	* \brief Class that controls the cluster hadronization algorithm.
	* \author Philip Stephens // cerr << *ch.currentEvent() << '\n';
	cerr << finalHadrons.size() << '\n';

	cerr << "Finished hadronizing \n";

	* \author Alberto Ribon
	*
	* This class is the main driver of the cluster hadronization: it is
	* responsible for the proper handling of all other specific collaborating
	* classes PartonSplitter, ClusterFinder, ColourReconnector, ClusterFissioner,
	* LightClusterDecayer, ClusterDecayer;
	* and for the storing of the produced particles in the Event record.
	*
	* @see PartonSplitter
	* @see ClusterFinder
	* @see ColourReconnector
	* @see ClusterFissioner
	* @see LightClusterDecayer
	* @see ClusterDecayer
	* @see Cluster
	* @see \ref ClusterHadronizationHandlerInterfaces "The interfaces"
	* defined for ClusterHadronizationHandler.
	*/
	-class ClusterHadronizationHandler: public HadronizationHandler {
	+class ClusterHadronizationHandler:
	+ public HadronizationHandler, public Reshuffler {

	public:

	/**
	* The main method which manages the all cluster hadronization.
	*
	* This routine directs "traffic". It determines which function is called
	* and on which particles/clusters. This function also handles the
	* situation of vetos on the hadronization.
	*/
	virtual void handle(EventHandler & ch, const tPVector & tagged,
	const Hint & hint);

	/**
	* It returns minimum virtuality^2 of partons to use in calculating
	* distances. It is used both in the Showering and Hadronization.
	*/
	Energy2 minVirtuality2() const
	{ return _minVirtuality2; }

	/**
	* It returns the maximum displacement that is allowed for a particle
	* (used to determine the position of a cluster with two components).
	*/
	Length maxDisplacement() const
	{ return _maxDisplacement; }

	/**
	* It returns true/false according if the soft underlying model
	* is switched on/off.
	*/
	bool isSoftUnderlyingEventON() const
	{ return _underlyingEventHandler; }

	/**
	* pointer to "this", the current HadronizationHandler.
	*/
	static const ClusterHadronizationHandler * currentHandler() {
	if(!currentHandler_){
	cerr<< " \nCreating new ClusterHadronizationHandler without input from infiles.";
	cerr<< " \nWhen using for example the string model ";
	cerr<< " hadronic decays are still treated by the Cluster model\n";
	currentHandler_=new ClusterHadronizationHandler();;
	}
	return currentHandler_;
	}

	+ /**
	+ * A pointer to a gluon mass generator for the reshuffling
	+ */
	+ Ptr<GluonMassGenerator>::tptr gluonMassGenerator() const {
	+ return gluonMassGenerator_;
	+ }
	+
	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.
	*/
	ClusterHadronizationHandler & operator=(const ClusterHadronizationHandler &) = delete;

	/**
	* This is a pointer to a Herwig::PartonSplitter object.
	*/
	PartonSplitterPtr _partonSplitter;

	/**
	* This is a pointer to a Herwig::ClusterFinder object.
	*/
	ClusterFinderPtr _clusterFinder;

	/**
	* This is a pointer to a Herwig::ColourReconnector object.
	*/
	ColourReconnectorPtr _colourReconnector;

	/**
	* This is a pointer to a Herwig::ClusterFissioner object.
	*/
	ClusterFissionerPtr _clusterFissioner;

	/**
	* This is a pointer to a Herwig::LightClusterDecayer object.
	*/
	LightClusterDecayerPtr _lightClusterDecayer;

	/**
	* This is a pointer to a Herwig::ClusterDecayer object.
	*/
	ClusterDecayerPtr _clusterDecayer;

	/**
	+ * Perform reshuffling to constituent masses.
	+ */
	+ bool reshuffle_ = false;
	+
	+ /**
	+ * Which type of reshuffling (global (default) or colour connected) is used
	+ */
	+ int reshuffleMode_ = 0;
	+
	+ /**
	+ * A pointer to a gluon mass generator for the reshuffling
	+ */
	+ Ptr<GluonMassGenerator>::ptr gluonMassGenerator_;
	+
	+ /**
	* The minimum virtuality^2 of partons to use in calculating
	* distances.
	*/
	Energy2 _minVirtuality2 = 0.1_GeV2;

	/**
	* The maximum displacement that is allowed for a particle
	* (used to determine the position of a cluster with two components).
	*/
	Length _maxDisplacement = 1.0e-10_mm;

	/**
	* The pointer to the Underlying Event handler.
	*/
	StepHdlPtr _underlyingEventHandler;

	/**
	* How to handle baryon-number clusters
	*/
	bool _reduceToTwoComponents = true;

	/**
	* Tag the constituents of the clusters as their parents
	*/
	void _setChildren(const ClusterVector & clusters) const;
	+
	+
	+ /**
	+ * Split the list of partons into colour connected sub-lists before reshuffling
	+ */
	+ void splitIntoColourSinglets(PVector thelist,
	+ vector<PVector>& reshufflelists);
	+

	/**
	* pointer to "this", the current HadronizationHandler.
	*/
	static ClusterHadronizationHandler * currentHandler_;


	};


	}

	#endif /* HERWIG_ClusterHadronizationHandler_H */
	diff --git a/Hadronization/GluonMassGenerator.cc b/Hadronization/GluonMassGenerator.cc
	new file mode 100644
	--- /dev/null
	+++ b/Hadronization/GluonMassGenerator.cc
	@@ -0,0 +1,63 @@
	+// -- C++ --
	+//
	+// GluonMassGenerator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	+// Copyright (C) 2002-2019 The Herwig Collaboration
	+//
	+// Herwig is licenced under version 3 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 GluonMassGenerator class.
	+//
	+
	+#include "GluonMassGenerator.h"
	+#include "ThePEG/Interface/ClassDocumentation.h"
	+#include "ThePEG/EventRecord/Particle.h"
	+#include "ThePEG/Repository/UseRandom.h"
	+#include "ThePEG/Repository/EventGenerator.h"
	+#include "ThePEG/Utilities/DescribeClass.h"
	+#include "ThePEG/StandardModel/StandardModelBase.h"
	+#include "ClusterHadronizationHandler.h"
	+
	+#include "ThePEG/Persistency/PersistentOStream.h"
	+#include "ThePEG/Persistency/PersistentIStream.h"
	+
	+using namespace Herwig;
	+
	+GluonMassGenerator::GluonMassGenerator() {}
	+
	+GluonMassGenerator::~GluonMassGenerator() {}
	+
	+IBPtr GluonMassGenerator::clone() const {
	+ return new_ptr(*this);
	+}
	+
	+IBPtr GluonMassGenerator::fullclone() const {
	+ return new_ptr(*this);
	+}
	+
	+// If needed, insert default implementations of virtual function defined
	+// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
	+
	+
	+void GluonMassGenerator::persistentOutput(PersistentOStream &) const {}
	+
	+void GluonMassGenerator::persistentInput(PersistentIStream &, int) {}
	+
	+
	+// * Attention * The following static variable is needed for the type
	+// description system in ThePEG. Please check that the template arguments
	+// are correct (the class and its base class), and that the constructor
	+// arguments are correct (the class name and the name of the dynamically
	+// loadable library where the class implementation can be found).
	+DescribeClass<GluonMassGenerator,HandlerBase>
	+ describeHerwigGluonMassGenerator("Herwig::GluonMassGenerator", "");
	+
	+void GluonMassGenerator::Init() {
	+
	+ static ClassDocumentation<GluonMassGenerator> documentation
	+ ("Dynamic gluon mass generation");
	+
	+}
	+
	diff --git a/Hadronization/GluonMassGenerator.h b/Hadronization/GluonMassGenerator.h
	new file mode 100644
	--- /dev/null
	+++ b/Hadronization/GluonMassGenerator.h
	@@ -0,0 +1,171 @@
	+// -- C++ --
	+//
	+// GluonMassGenerator.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	+// Copyright (C) 2002-2019 The Herwig Collaboration
	+//
	+// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	+//
	+#ifndef Herwig_GluonMassGenerator_H
	+#define Herwig_GluonMassGenerator_H
	+//
	+// This is the declaration of the GluonMassGenerator class.
	+//
	+
	+#include "ThePEG/Handlers/HandlerBase.h"
	+#include "ThePEG/EventRecord/Particle.h"
	+
	+namespace Herwig {
	+
	+using namespace ThePEG;
	+
	+/**
	+ * \ingroup Hadronization
	+ * \brief Dynamic gluon mass generator; the default returns a constant mass.
	+ *
	+ * @see \ref GluonMassGeneratorInterfaces "The interfaces"
	+ * defined for GluonMassGenerator.
	+ */
	+class GluonMassGenerator: public HandlerBase {
	+
	+public:
	+
	+ /** @name Standard constructors and destructors. */
	+ //@{
	+ /**
	+ * The default constructor.
	+ */
	+ GluonMassGenerator();
	+
	+ /**
	+ * The destructor.
	+ */
	+ virtual ~GluonMassGenerator();
	+ //@}
	+
	+public:
	+
	+ /**
	+ * Generate a single gluon mass with possible reference to a hard
	+ * scale Q and up to a maximum value
	+ */
	+ virtual Energy generate(Energy, Energy) const {
	+ return generate();
	+ }
	+
	+ /**
	+ * Generate a single gluon mass with possible reference to a hard
	+ * scale Q
	+ */
	+ virtual Energy generate(Energy) const {
	+ return generate();
	+ }
	+
	+ /**
	+ * Generate a single gluon mass without further constraints
	+ */
	+ virtual Energy generate() const {
	+ return getParticleData(ThePEG::ParticleID::g)->constituentMass();
	+ }
	+
	+ /**
	+ * Generate a list of n gluon masses, with a maximum available energy
	+ */
	+ list<Energy> generateMany(size_t n, Energy QMax) const {
	+ list<Energy> res;
	+ Energy m0, mu, md, ms, mg, mgmax, summg;
	+
	+ mu=getParticleData(ThePEG::ParticleID::u)->constituentMass();
	+ md=getParticleData(ThePEG::ParticleID::d)->constituentMass();
	+ ms=getParticleData(ThePEG::ParticleID::s)->constituentMass();
	+
	+ m0=md;
	+ if(mu<m0){m0=mu;}
	+ if(ms<m0){m0=ms;}
	+
	+ if( QMax<2.0m0n ){
	+ throw Exception() << "cannot reshuffle to constituent mass shells" << Exception::eventerror;
	+ }
	+
	+ bool repeat=true;
	+
	+ while( repeat ){
	+ repeat=false;
	+ summg = 0.0*GeV;
	+ res.clear();
	+ for( size_t k = 0; k < n; ++k ){
	+ mg = generate();
	+ res.push_back(mg);
	+ summg += mg;
	+ if( summg > QMax - 2.0m0(n-k-1) ){
	+ repeat=true;
	+ break;
	+ }
	+ }
	+ }
	+
	+ return res;
	+
	+ }
	+
	+public:
	+
	+ /** @name Functions used by the persistent I/O system. */
	+ //@{
	+ /**
	+ * Function used to write out object persistently.
	+ * @param os the persistent output stream written to.
	+ */
	+ void persistentOutput(PersistentOStream & os) const;
	+
	+ /**
	+ * Function used to read in object persistently.
	+ * @param is the persistent input stream read from.
	+ * @param version the version number of the object when written.
	+ */
	+ void persistentInput(PersistentIStream & is, int version);
	+ //@}
	+
	+ /**
	+ * The standard Init function used to initialize the interfaces.
	+ * Called exactly once for each class by the class description system
	+ * before the main function starts or
	+ * when this class is dynamically loaded.
	+ */
	+ static void Init();
	+
	+protected:
	+
	+ /** @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;
	+ //@}
	+
	+
	+// If needed, insert declarations of virtual function defined in the
	+// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
	+
	+
	+private:
	+
	+ /**
	+ * The assignment operator is private and must never be called.
	+ * In fact, it should not even be implemented.
	+ */
	+ GluonMassGenerator & operator=(const GluonMassGenerator &);
	+
	+};
	+
	+}
	+
	+#endif /* Herwig_GluonMassGenerator_H */
	diff --git a/Hadronization/Makefile.am b/Hadronization/Makefile.am
	--- a/Hadronization/Makefile.am
	+++ b/Hadronization/Makefile.am
	@@ -1,17 +1,18 @@
	noinst_LTLIBRARIES = libHwHadronization.la
	libHwHadronization_la_SOURCES = \
	CheckId.cc CheckId.h \
	CluHadConfig.h \
	Cluster.h Cluster.cc Cluster.fh \
	ClusterDecayer.cc ClusterDecayer.h ClusterDecayer.fh \
	ClusterFinder.cc ClusterFinder.h ClusterFinder.fh \
	ClusterFissioner.cc ClusterFissioner.h ClusterFissioner.fh \
	ClusterHadronizationHandler.cc ClusterHadronizationHandler.h \
	ClusterHadronizationHandler.fh \
	ColourReconnector.cc ColourReconnector.h ColourReconnector.fh\
	+GluonMassGenerator.h GluonMassGenerator.cc \
	HadronSelector.cc HadronSelector.h HadronSelector.fh\
	Hw64Selector.cc Hw64Selector.h Hw64Selector.fh\
	HwppSelector.cc HwppSelector.h HwppSelector.fh\
	LightClusterDecayer.cc LightClusterDecayer.h LightClusterDecayer.fh \
	PartonSplitter.cc PartonSplitter.h PartonSplitter.fh \
	SpinHadronizer.h SpinHadronizer.cc
	diff --git a/Shower/Dipole/Utility/ConstituentReshuffler.cc b/Shower/Dipole/Utility/ConstituentReshuffler.cc
	--- a/Shower/Dipole/Utility/ConstituentReshuffler.cc
	+++ b/Shower/Dipole/Utility/ConstituentReshuffler.cc
	@@ -1,663 +1,622 @@
	// -- C++ --
	//
	// ConstituentReshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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 ConstituentReshuffler class.
	//

	#include <config.h>
	#include "ConstituentReshuffler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"

	#include <limits>

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "DipolePartonSplitter.h"

	#include "Herwig/Utilities/GSLBisection.h"

	#include "Herwig/Shower/Dipole/DipoleShowerHandler.h"

	#include "Herwig/Shower/ShowerHandler.h"

	using namespace Herwig;

	ConstituentReshuffler::ConstituentReshuffler()
	: HandlerBase() {}

	ConstituentReshuffler::~ConstituentReshuffler() {}

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

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

	-double ConstituentReshuffler::ReshuffleEquation::aUnit() {
	- return 1.;
	-}
	-
	-double ConstituentReshuffler::ReshuffleEquation::vUnit() {
	- return 1.;
	-}
	-
	-double ConstituentReshuffler::DecayReshuffleEquation::aUnit() {
	- return 1.;
	-}
	-
	-double ConstituentReshuffler::DecayReshuffleEquation::vUnit() {
	- return 1.;
	-}
	-
	-
	-double ConstituentReshuffler::ReshuffleEquation::operator() (double xi) const {
	-
	- double r = - w/GeV;
	-
	- for (PList::iterator p = p_begin; p != p_end; ++p) {
	- r += sqrt(sqr((**p).dataPtr()->constituentMass()) +
	- xixi(sqr((p).momentum().t())-sqr((p).dataPtr()->mass()))) / GeV;
	- }
	- return r;
	-
	-}
	-
	-double ConstituentReshuffler::DecayReshuffleEquation::operator() (double xi) const {
	- double r = - w/GeV;
	-
	- for (PList::iterator pIt = p_begin; pIt != p_end; ++pIt) {
	- r += sqrt(sqr((**pIt).dataPtr()->constituentMass()) +
	- xixi(sqr((pIt).momentum().t())-sqr((pIt).dataPtr()->mass()))) / GeV;
	- }
	-
	- for (PList::iterator rIt = r_begin; rIt != r_end; ++rIt) {
	- r += sqrt(sqr((**rIt).momentum().m()) +
	- xixi(sqr((rIt).momentum().t())-sqr((rIt).momentum().m()))) / GeV;
	- }
	-
	- return r;
	-
	-}
	-
	void ConstituentReshuffler::reshuffle(PList& out,
	PPair& in,
	PList& intermediates,
	const bool decay,
	PList& decayPartons,
	PList& decayRecoilers) {

	assert(ShowerHandler::currentHandler()->retConstituentMasses());

	if ( !decay ) {

	if (out.size() == 0)
	return;

	if (out.size() == 1) {

	PPtr recoiler;
	PPtr parton = out.front();

	if (DipolePartonSplitter::colourConnected(parton,in.first) &&
	DipolePartonSplitter::colourConnected(parton,in.second)) {
	if (UseRandom::rnd() < .5)
	recoiler = in.first;
	else
	recoiler = in.second;
	} else if (DipolePartonSplitter::colourConnected(parton,in.first)) {
	recoiler = in.first;
	} else if (DipolePartonSplitter::colourConnected(parton,in.second)) {
	recoiler = in.second;
	} else assert(false);

	assert(abs(recoiler->momentum().vect().perp2()/GeV2) < 1e-6);

	double sign = recoiler->momentum().z() < 0.*GeV ? -1. : 1.;

	Energy2 qperp2 = parton->momentum().perp2();

	if (qperp2/GeV2 < Constants::epsilon) {
	// no emission off a 2 -> singlet process which
	// needed a single forced splitting: should never happen (?)
	assert(false);
	throw Veto();
	}

	Energy2 m2 = sqr(parton->dataPtr()->constituentMass());

	Energy abs_q = parton->momentum().vect().mag();
	Energy qz = parton->momentum().z();
	Energy abs_pz = recoiler->momentum().t();
	assert(abs_pz > 0.*GeV);

	Energy xi_pz = sign(2.qperp2abs_pz + m2(abs_q + signqz))/(2.qperp2);
	Energy x_qz = (2.qperp2qz + m2(qz+signabs_q))/(2.*qperp2);

	Lorentz5Momentum recoiler_momentum
	(0.GeV,0.GeV,xi_pz,xi_pz < 0.*GeV ? - xi_pz : xi_pz);

	recoiler_momentum.rescaleMass();

	Lorentz5Momentum parton_momentum
	(parton->momentum().x(),parton->momentum().y(),x_qz,sqrt(m2+qperp2+x_qz*x_qz));

	parton_momentum.rescaleMass();

	PPtr n_parton = new_ptr(Particle(parton->dataPtr()));
	n_parton->set5Momentum(parton_momentum);

	DipolePartonSplitter::change(parton,n_parton,false);

	out.pop_front();
	intermediates.push_back(parton);
	out.push_back(n_parton);

	PPtr n_recoiler = new_ptr(Particle(recoiler->dataPtr()));
	n_recoiler->set5Momentum(recoiler_momentum);

	DipolePartonSplitter::change(recoiler,n_recoiler,true);

	intermediates.push_back(recoiler);

	if (recoiler == in.first) {
	in.first = n_recoiler;
	}

	if (recoiler == in.second) {
	in.second = n_recoiler;
	}

	return;

	}

	}

	Energy zero (0.*GeV);
	Lorentz5Momentum Q (zero,zero,zero,zero);

	for (PList::iterator p = out.begin();
	p != out.end(); ++p) {
	Q += (**p).momentum();
	}

	Boost beta = Q.findBoostToCM();

	list<Lorentz5Momentum> mbackup;

	bool need_boost = (beta.mag2() > Constants::epsilon);

	if (need_boost) {

	for (PList::iterator p = out.begin();
	p != out.end(); ++p) {
	Lorentz5Momentum mom = (**p).momentum();
	mbackup.push_back(mom);
	(**p).set5Momentum(mom.boost(beta));
	}

	}

	double xi;

	// Only partons
	if ( decayRecoilers.size()==0 ) {
	- ReshuffleEquation solve (Q.m(),out.begin(),out.end());
	+ list<Energy> masses;
	+ for ( auto p : out )
	+ masses.push_back(p->dataPtr()->constituentMass());
	+ ReshuffleEquation<PList::iterator,list<Energy>::const_iterator>
	+ solve (Q.m(),out.begin(),out.end(),
	+ masses.begin(),masses.end());

	GSLBisection solver(1e-10,1e-8,10000);

	try {
	xi = solver.value(solve,0.0,1.1);
	} catch (GSLBisection::GSLerror) {
	throw DipoleShowerHandler::RedoShower();
	} catch (GSLBisection::IntervalError) {
	throw DipoleShowerHandler::RedoShower();
	}
	}

	// Partons and decaying recoilers
	else {
	DecayReshuffleEquation solve (Q.m(),decayPartons.begin(),decayPartons.end(),decayRecoilers.begin(),decayRecoilers.end());

	GSLBisection solver(1e-10,1e-8,10000);

	try {
	xi = solver.value(solve,0.0,1.1);
	} catch (GSLBisection::GSLerror) {
	throw DipoleShowerHandler::RedoShower();
	} catch (GSLBisection::IntervalError) {
	throw DipoleShowerHandler::RedoShower();
	}
	}


	PList reshuffled;

	list<Lorentz5Momentum>::const_iterator backup_it;
	if (need_boost)
	backup_it = mbackup.begin();


	// Reshuffling of non-decaying partons only
	if ( decayRecoilers.size()==0 ) {

	for (PList::iterator p = out.begin();
	p != out.end(); ++p) {

	PPtr rp = new_ptr(Particle((**p).dataPtr()));

	DipolePartonSplitter::change(*p,rp,false);

	Lorentz5Momentum rm;

	rm = Lorentz5Momentum (xi(*p).momentum().x(),
	xi(*p).momentum().y(),
	xi(*p).momentum().z(),
	sqrt(sqr((**p).dataPtr()->constituentMass()) +
	xixi(sqr((p).momentum().t())-sqr((p).dataPtr()->mass()))));

	rm.rescaleMass();

	if (need_boost) {
	(*p).set5Momentum(backup_it);
	++backup_it;
	rm.boost(-beta);
	}

	rp->set5Momentum(rm);

	intermediates.push_back(*p);
	reshuffled.push_back(rp);


	}
	}

	// For the case of a decay process with non-partonic recoilers
	else {
	assert ( decay );

	for (PList::iterator p = out.begin();
	p != out.end(); ++p) {

	// Flag to update spinInfo
	bool updateSpin = false;

	PPtr rp = new_ptr(Particle((**p).dataPtr()));

	DipolePartonSplitter::change(*p,rp,false);

	Lorentz5Momentum rm;

	// If the particle is a parton and not a recoiler
	if ( find( decayRecoilers.begin(), decayRecoilers.end(), *p ) == decayRecoilers.end() ) {
	rm = Lorentz5Momentum (xi(*p).momentum().x(),
	xi(*p).momentum().y(),
	xi(*p).momentum().z(),
	sqrt(sqr((**p).dataPtr()->constituentMass()) +
	xixi(sqr((p).momentum().t())-sqr((p).dataPtr()->mass()))));
	}


	// Otherwise the parton is a recoiler
	// and its invariant mass must be preserved
	else {
	if ( (*p)-> spinInfo() )
	updateSpin = true;
	rm = Lorentz5Momentum (xi(*p).momentum().x(),
	xi(*p).momentum().y(),
	xi(*p).momentum().z(),
	sqrt(sqr((**p).momentum().m()) +
	xixi(sqr((p).momentum().t())-sqr((p).momentum().m()))));
	}

	rm.rescaleMass();

	if (need_boost) {
	(*p).set5Momentum(backup_it);
	++backup_it;
	rm.boost(-beta);
	}

	rp->set5Momentum(rm);

	// Update SpinInfo if required
	if ( updateSpin )
	updateSpinInfo(*p, rp);

	intermediates.push_back(*p);
	reshuffled.push_back(rp);


	}
	}

	out.clear();
	out.splice(out.end(),reshuffled);

	}


	void ConstituentReshuffler::hardProcDecayReshuffle(PList& decaying,
	PList& eventOutgoing,
	PList& eventHard,
	PPair& eventIncoming,
	PList& eventIntermediates) {

	// Note, when this function is called, the particle pointers
	// in theDecays/decaying are those prior to the showering.
	// Here we find the newest pointers in the outgoing.
	// The update of the PPtrs in theDecays is done in DipoleShowerHandler::constituentReshuffle()
	// as this needs to be done if ConstituentReshuffling is switched off.


	//Make sure the shower should return constituent masses:
	assert(ShowerHandler::currentHandler()->retConstituentMasses());

	// Find the outgoing decaying particles
	PList recoilers;
	for ( PList::iterator decIt = decaying.begin(); decIt != decaying.end(); ++decIt) {

	// First find the particles in the intermediates
	PList::iterator pos = find(eventIntermediates.begin(),eventIntermediates.end(), *decIt);

	// Colourless particle or coloured particle that did not radiate.
	if(pos==eventIntermediates.end()) {

	// Check that this is not a particle from a subsequent decay.
	// e.g. the W from a top decay from an LHE file.
	if ( find( eventHard.begin(), eventHard.end(), *decIt ) == eventHard.end() &&
	find( eventOutgoing.begin(), eventOutgoing.end(), *decIt ) == eventOutgoing.end() )
	continue;

	else
	recoilers.push_back( *decIt );

	}

	// Coloured decaying particle that radiated
	else {
	PPtr unstable = *pos;
	while(!unstable->children().empty()) {
	unstable = unstable->children()[0];
	}
	assert( find( eventOutgoing.begin(),eventOutgoing.end(), unstable ) != eventOutgoing.end() );

	recoilers.push_back( unstable );
	}
	}

	// Make a list of partons
	PList partons;
	for ( PList::iterator outPos = eventOutgoing.begin(); outPos != eventOutgoing.end(); ++outPos ) {
	if ( find (recoilers.begin(), recoilers.end(), *outPos ) == recoilers.end() ) {
	partons.push_back( *outPos );
	}
	}


	// If no outgoing partons, do nothing
	if ( partons.size() == 0 ){
	return;
	}

	// Otherwise reshuffling needs to be done.

	// If there is only one parton, attempt to reshuffle with
	// the incoming to be consistent with the reshuffle for a
	// hard process with no decays.
	else if ( partons.size() == 1 &&
	( DipolePartonSplitter::colourConnected(partons.front(),eventIncoming.first) \|\|
	DipolePartonSplitter::colourConnected(partons.front(),eventIncoming.second) ) ) {

	// Erase the parton from the event outgoing
	eventOutgoing.erase( find( eventOutgoing.begin(), eventOutgoing.end(), partons.front() ) );
	// Perform the reshuffle, this update the intermediates and the incoming
	reshuffle(partons, eventIncoming, eventIntermediates);
	// Update the outgoing
	eventOutgoing.push_back(partons.front());
	return;
	}

	// If reshuffling amongst the incoming is not possible
	// or if we have multiple outgoing partons.
	else {

	// Create a complete list of the outgoing from the process
	PList out;
	// Make an empty list for storing the new intermediates
	PList intermediates;
	// Empty incoming particles pair
	PPair in;

	// A single parton which cannot be reshuffled
	// with the incoming.
	if ( partons.size() == 1 ) {

	// Populate the out for the reshuffling
	out.insert(out.end(),partons.begin(),partons.end());
	out.insert(out.end(),recoilers.begin(),recoilers.end());
	assert( out.size() > 1 );

	// Perform the reshuffle with the temporary particle lists
	reshuffle(out, in, intermediates, true, partons, recoilers);
	}

	// If there is more than one parton, reshuffle only
	// amongst the partons
	else {
	assert(partons.size() > 1);

	// Populate the out for the reshuffling
	out.insert(out.end(),partons.begin(),partons.end());
	assert( out.size() > 1 );

	// Perform the reshuffle with the temporary particle lists
	reshuffle(out, in, intermediates, true);
	}

	// Update the dipole event record
	updateEvent(intermediates, eventIntermediates, out, eventOutgoing, eventHard );
	return;
	}
	}


	void ConstituentReshuffler::decayReshuffle(PerturbativeProcessPtr& decayProc,
	PList& eventOutgoing,
	PList& eventHard,
	PList& eventIntermediates ) {

	// Separate particles into those to be assigned constituent masses
	// i.e. non-decaying coloured partons
	// and those which must only absorb recoil
	// i.e. non-coloured and decaying particles
	PList partons;
	PList recoilers;

	//Make sure the shower should return constituent masses:
	assert(ShowerHandler::currentHandler()->retConstituentMasses());


	// Populate the particle lists from the outgoing of the decay process
	for( unsigned int ix = 0; ix<decayProc->outgoing().size(); ++ix) {

	// Identify recoilers
	if ( !decayProc->outgoing()[ix].first->coloured() \|\|
	ShowerHandler::currentHandler()->decaysInShower(decayProc->outgoing()[ix].first->id() ) )
	recoilers.push_back(decayProc->outgoing()[ix].first);

	else
	partons.push_back(decayProc->outgoing()[ix].first);
	}

	// If there are no outgoing partons, then no reshuffling
	// needs to be done
	if ( partons.size() == 0 )
	return;

	// Reshuffling needs to be done:
	else {

	// Create a complete list of the outgoing from the process
	PList out;
	// Make an empty list for storing the new intermediates
	PList intermediates;
	// Empty incoming particles pair
	PPair in;


	// SW - 15/06/2018, 31/01/2019 - Always include 'recoilers' in
	// reshuffling, regardless of the number of partons to be put on their
	// constituent mass shell. This is because reshuffling between 2 partons
	// frequently leads to a redoShower exception. This treatment is
	// consistent with the AO shower

	// Populate the out for the reshuffling
	out.insert(out.end(),partons.begin(),partons.end());
	out.insert(out.end(),recoilers.begin(),recoilers.end());
	assert( out.size() > 1 );

	// Perform the reshuffle with the temporary particle lists
	reshuffle(out, in, intermediates, true, partons, recoilers);

	// Update the dipole event record and the decay process
	updateEvent(intermediates, eventIntermediates, out, eventOutgoing, eventHard, decayProc );
	return;
	}
	}


	void ConstituentReshuffler::updateEvent( PList& intermediates,
	PList& eventIntermediates,
	#ifndef NDEBUG
	PList& out,
	#else
	PList&,
	#endif
	PList& eventOutgoing,
	PList& eventHard,
	PerturbativeProcessPtr decayProc ) {

	// Loop over the new intermediates following the reshuffling
	for (PList::iterator p = intermediates.begin();
	p != intermediates.end(); ++p) {

	// Update the event record intermediates
	eventIntermediates.push_back(*p);

	// Identify the reshuffled particle
	assert( (*p)->children().size()==1 );
	PPtr reshuffled = (*p)->children()[0];
	assert( find(out.begin(), out.end(), reshuffled) != out.end() );

	// Update the event record outgoing
	PList::iterator posOut = find(eventOutgoing.begin(), eventOutgoing.end(), *p);

	if ( posOut != eventOutgoing.end() ) {
	eventOutgoing.erase(posOut);
	eventOutgoing.push_back(reshuffled);
	}

	else {
	PList::iterator posHard = find(eventHard.begin(), eventHard.end(), *p);
	assert( posHard != eventHard.end() );
	eventHard.erase(posHard);
	eventHard.push_back(reshuffled);
	}

	// Replace the particle in the the decay process outgoing
	if ( decayProc ) {
	vector<pair<PPtr,PerturbativeProcessPtr> >::iterator decayOutIt = decayProc->outgoing().end();
	for ( decayOutIt = decayProc->outgoing().begin();
	decayOutIt!= decayProc->outgoing().end(); ++decayOutIt ) {

	if ( decayOutIt->first == *p ){
	break;
	}
	}
	assert( decayOutIt != decayProc->outgoing().end() );
	decayOutIt->first = reshuffled;
	}
	}
	}

	void ConstituentReshuffler::updateSpinInfo( PPtr& oldPart,
	PPtr& newPart ) {

	const Lorentz5Momentum& oldMom = oldPart->momentum();
	const Lorentz5Momentum& newMom = newPart->momentum();

	// Rotation from old momentum to +ve z-axis
	LorentzRotation oldToZAxis;
	Axis axisOld(oldMom.vect().unit());
	if( axisOld.perp2() > 1e-12 ) {
	double sinth(sqrt(1.-sqr(axisOld.z())));
	oldToZAxis.rotate( -acos(axisOld.z()),Axis(-axisOld.y()/sinth,axisOld.x()/sinth,0.));
	}

	// Rotation from new momentum to +ve z-axis
	LorentzRotation newToZAxis;
	Axis axisNew(newMom.vect().unit());
	if( axisNew.perp2() > 1e-12 ) {
	double sinth(sqrt(1.-sqr(axisNew.z())));
	newToZAxis.rotate( -acos(axisNew.z()),Axis(-axisNew.y()/sinth,axisNew.x()/sinth,0.));
	}

	// Boost from old momentum to new momentum along z-axis
	Lorentz5Momentum momOldRotated = oldToZAxis*Lorentz5Momentum(oldMom);
	Lorentz5Momentum momNewRotated = newToZAxis*Lorentz5Momentum(newMom);

	Energy2 a = sqr(momOldRotated.z()) + sqr(momNewRotated.t());
	Energy2 b = 2.momOldRotated.t()momOldRotated.z();
	Energy2 c = sqr(momOldRotated.t()) - sqr(momNewRotated.t());
	double beta;

	// The rotated momentum should always lie along the +ve z-axis
	if ( momOldRotated.z() > ZERO )
	beta = (-b + sqrt(sqr(b)-4.ac)) / 2. / a;
	else
	beta = (-b - sqrt(sqr(b)-4.ac)) / 2. / a;

	LorentzRotation boostOldToNew(0., 0., beta);

	// Total transform
	LorentzRotation transform = (newToZAxis.inverse())boostOldToNewoldToZAxis;

	// Assign the same spin info to the old and new particles
	newPart->spinInfo(oldPart->spinInfo());
	newPart->spinInfo()->transform(oldMom, transform);
	}

	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void ConstituentReshuffler::persistentOutput(PersistentOStream &) const {
	}

	void ConstituentReshuffler::persistentInput(PersistentIStream &, int) {
	}

	ClassDescription<ConstituentReshuffler> ConstituentReshuffler::initConstituentReshuffler;
	// Definition of the static class description member.

	void ConstituentReshuffler::Init() {

	static ClassDocumentation<ConstituentReshuffler> documentation
	("The ConstituentReshuffler class implements reshuffling "
	"of partons on their nominal mass shell to their constituent "
	"mass shells.");

	}

	diff --git a/Shower/Dipole/Utility/ConstituentReshuffler.h b/Shower/Dipole/Utility/ConstituentReshuffler.h
	--- a/Shower/Dipole/Utility/ConstituentReshuffler.h
	+++ b/Shower/Dipole/Utility/ConstituentReshuffler.h
	@@ -1,299 +1,274 @@
	// -- C++ --
	//
	// ConstituentReshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_ConstituentReshuffler_H
	#define HERWIG_ConstituentReshuffler_H
	//
	// This is the declaration of the ConstituentReshuffler class.
	//

	#include "ThePEG/Handlers/HandlerBase.h"
	#include "ThePEG/Utilities/Exception.h"
	#include "Herwig/Shower/PerturbativeProcess.h"
	+#include "Herwig/Utilities/Reshuffler.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer, Stephen Webster
	*
	* \brief The ConstituentReshuffler class implements reshuffling
	* of partons on their nominal mass shell to their constituent
	* mass shells.
	*
	*/
	-class ConstituentReshuffler: public HandlerBase {
	+class ConstituentReshuffler: public HandlerBase, public Reshuffler {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	ConstituentReshuffler();

	/**
	* The destructor.
	*/
	virtual ~ConstituentReshuffler();
	//@}

	public:

	/**
	* Reshuffle the outgoing partons to constituent
	* masses. Optionally, incoming partons are given
	* to absorb recoils. Add the non-reshuffled partons
	* to the intermediates list. Throw ConstituentReshufflerProblem
	* if a numerical problem prevents the solution of
	* the reshuffling equation.
	*/
	void reshuffle(PList& out,
	PPair& in,
	PList& intermediates,
	const bool decay,
	PList& decayPartons,
	PList& decayRecoilers);

	/**
	* Reshuffle the outgoing partons to constituent
	* masses. Optionally, incoming partons are given
	* to absorb recoils. Add the non-reshuffled partons
	* to the intermediates list. Throw ConstituentReshufflerProblem
	* if a numerical problem prevents the solution of
	* the reshuffling equation.
	*/
	void reshuffle(PList& out,
	PPair& in,
	PList& intermediates,
	const bool decay=false) {

	PList decayPartons;
	PList decayRecoilers;

	reshuffle(out,
	in,
	intermediates,
	decay,
	decayPartons,
	decayRecoilers);
	}


	/**
	* Reshuffle the outgoing partons following the showering
	* of the initial hard interaction to constituent masses,
	* for the case of outgoing decaying particles.
	* Throw ConstituentReshufflerProblem
	* if a numerical problem prevents the solution of
	* the reshuffling equation.
	*/
	void hardProcDecayReshuffle(PList& decaying,
	PList& eventOutgoing,
	PList& eventHard,
	PPair& eventIncoming,
	PList& eventIntermediates) ;

	/**
	* Reshuffle the outgoing partons following the showering
	* of a particle decay to constituent masses.
	* Throw ConstituentReshufflerProblem
	* if a numerical problem prevents the solution of
	* the reshuffling equation.
	*/
	void decayReshuffle(PerturbativeProcessPtr& decayProc,
	PList& eventOutgoing,
	PList& eventHard,
	PList& eventIntermediates) ;

	/**
	* Update the dipole event record and, if appropriate,
	* the relevant decay process.
	**/
	void updateEvent( PList& intermediates,
	PList& eventIntermediates,
	PList& out,
	PList& eventOutgoing,
	PList& eventHard,
	PerturbativeProcessPtr decayProc = PerturbativeProcessPtr() ) ;


	/**
	* Update the spinInfo of a particle following reshuffling
	* to take account of the change in momentum.
	* Used only for unstable particles that need to be dealt with.
	**/
	void updateSpinInfo( PPtr& oldPart,
	PPtr& newPart ) ;

	protected:

	/**
	* The function object defining the equation
	- * to be solved.
	- */
	- struct ReshuffleEquation {
	-
	- ReshuffleEquation (Energy q,
	- PList::iterator m_begin,
	- PList::iterator m_end)
	- : w(q), p_begin(m_begin), p_end(m_end) {}
	-
	- typedef double ArgType;
	- typedef double ValType;
	-
	- static double aUnit();
	- static double vUnit();
	-
	- double operator() (double xi) const;
	-
	- Energy w;
	-
	- PList::iterator p_begin;
	- PList::iterator p_end;
	-
	- };
	-
	- /**
	- * The function object defining the equation
	* to be solved in the case of separate recoilers
	* TODO - refine the whole implementation of separate partons and recoilers
	*/
	struct DecayReshuffleEquation {

	DecayReshuffleEquation (Energy q,
	PList::iterator m_begin,
	PList::iterator m_end,
	PList::iterator n_begin,
	PList::iterator n_end)
	: w(q), p_begin(m_begin), p_end(m_end), r_begin(n_begin), r_end(n_end) {}

	typedef double ArgType;
	typedef double ValType;

	static double aUnit();
	static double vUnit();

	double operator() (double xi) const;

	Energy w;

	PList::iterator p_begin;
	PList::iterator p_end;

	PList::iterator r_begin;
	PList::iterator r_end;
	};



	public:

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

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

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

	protected:

	/** @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;
	//@}


	// If needed, insert declarations of virtual function defined in the
	// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).


	private:

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

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the ConstituentReshuffler class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::ConstituentReshuffler>
	: public ClassTraitsBase<Herwig::ConstituentReshuffler> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::ConstituentReshuffler"; }
	/**
	* The name of a file containing the dynamic library where the class
	* ConstituentReshuffler is implemented. It may also include several, space-separated,
	* libraries if the class ConstituentReshuffler depends on other classes (base classes
	* excepted). In this case the listed libraries will be dynamically
	* linked in the order they are specified.
	*/
	static string library() { return "HwDipoleShower.so"; }
	};

	/** @endcond */

	}

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

	#include "SudakovFormFactor.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.h"
	#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "Herwig/Shower/QTilde/QTildeShowerHandler.h"
	#include "Herwig/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h"
	#include "Herwig/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h"
	#include "Herwig/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h"
	#include "Herwig/Shower/QTilde/Kinematics/KinematicHelpers.h"
	#include "SudakovCutOff.h"

	#include <array>
	using std::array;

	using namespace Herwig;

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

	void SudakovFormFactor::persistentOutput(PersistentOStream & os) const {
	os << splittingFn_ << alpha_ << pdfmax_ << particles_ << pdffactor_ << cutoff_;
	}

	void SudakovFormFactor::persistentInput(PersistentIStream & is, int) {
	is >> splittingFn_ >> alpha_ >> pdfmax_ >> particles_ >> pdffactor_ >> cutoff_;
	}

	void SudakovFormFactor::Init() {

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

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

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

	static Reference<SudakovFormFactor,SudakovCutOff>
	interfaceCutoff("Cutoff",
	"A reference to the SudakovCutOff object",
	&Herwig::SudakovFormFactor::cutoff_,
	false, false, true, false);

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

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


	}

	bool SudakovFormFactor::alphaSVeto(Energy2 pt2) const {
	double ratio=alphaSVetoRatio(pt2,1.);
	return UseRandom::rnd() > ratio;
	}

	double SudakovFormFactor::alphaSVetoRatio(Energy2 pt2, double factor) const {
	factor *= ShowerHandler::currentHandler()->renormalizationScaleFactor();
	- return alpha_->ratio(pt2, factor);
	+ return alpha_->showerRatio(pt2, factor);
	}


	bool SudakovFormFactor::PDFVeto(const Energy2 t, const double x,
	const tcPDPtr parton0, const tcPDPtr parton1,
	Ptr<BeamParticleData>::transient_const_pointer beam) const {
	double ratio=PDFVetoRatio(t,x,parton0,parton1,beam,1.);
	return UseRandom::rnd() > ratio;
	}

	double SudakovFormFactor::PDFVetoRatio(const Energy2 t, const double x,
	const tcPDPtr parton0, const tcPDPtr parton1,
	Ptr<BeamParticleData>::transient_const_pointer beam,double factor) const {
	assert(pdf_);
	Energy2 theScale = t * sqr(ShowerHandler::currentHandler()->factorizationScaleFactor()*factor);
	if (theScale < sqr(freeze_)) theScale = sqr(freeze_);

	const double newpdf=pdf_->xfx(beam,parton0,theScale,x/z());
	if(newpdf<=0.) return 0.;

	const double oldpdf=pdf_->xfx(beam,parton1,theScale,x);
	if(oldpdf<=0.) return 1.;

	const double ratio = newpdf/oldpdf;
	double maxpdf = pdfmax_;

	switch (pdffactor_) {
	case 0: break;
	case 1: maxpdf /= z(); break;
	case 2: maxpdf /= 1.-z(); break;
	case 3: maxpdf /= (z()*(1.-z())); break;
	case 4: maxpdf /= sqrt(z()); break;
	case 5: maxpdf *= sqrt(z()); break;
	default :
	throw Exception() << "SudakovFormFactor::PDFVetoRatio invalid PDFfactor = "
	<< pdffactor_ << Exception::runerror;

	}

	if (ratio > maxpdf) {
	generator()->log() << "PDFVeto warning: Ratio > " << name()
	<< ":PDFmax (by a factor of "
	<< ratio/maxpdf <<") for "
	<< parton0->PDGName() << " to "
	<< parton1->PDGName() << "\n";
	}
	return ratio/maxpdf ;
	}

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

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

	void SudakovFormFactor::guesstz(Energy2 t1,unsigned int iopt,
	const IdList &ids,
	double enhance,bool ident,
	double detune,
	Energy2 &t_main, double &z_main){
	unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_;
	double lower = splittingFn_->integOverP(zlimits_.first ,ids,pdfopt);
	double upper = splittingFn_->integOverP(zlimits_.second,ids,pdfopt);
	double c = 1./((upper - lower)
	- * alpha_->overestimateValue()/Constants::twopienhancedetune);
	+ * alpha_->showerOverestimateValue()/Constants::twopienhancedetune);
	double r = UseRandom::rnd();
	assert(iopt<=2);
	if(iopt==1) {
	c/=pdfmax_;
	//symmetry of FS gluon splitting
	if(ident) c*=0.5;
	}
	else if(iopt==2) c*=-1.;
	// guessing t
	if(iopt!=2 \|\| c*log(r)<log(Constants::MaxEnergy2/t1)) {
	t_main = t1*pow(r,c);
	}
	else
	t_main = Constants::MaxEnergy2;
	// guessing z
	z_main = splittingFn_->invIntegOverP(lower + UseRandom::rnd()
	*(upper - lower),ids,pdfopt);
	}

	bool SudakovFormFactor::guessTimeLike(Energy2 &t,Energy2 tmin,double enhance,
	double detune) {
	Energy2 told = t;
	// calculate limits on z and if lower>upper return
	if(!computeTimeLikeLimits(t)) return false;
	// guess values of t and z
	guesstz(told,0,ids_,enhance,ids_[1]==ids_[2],detune,t,z_);
	// actual values for z-limits
	if(!computeTimeLikeLimits(t)) return false;
	if(t<tmin) {
	t=-1.0*GeV2;
	return false;
	}
	else
	return true;
	}

	bool SudakovFormFactor::guessSpaceLike(Energy2 &t, Energy2 tmin, const double x,
	double enhance,
	double detune) {
	Energy2 told = t;
	// calculate limits on z if lower>upper return
	if(!computeSpaceLikeLimits(t,x)) return false;
	// guess values of t and z
	guesstz(told,1,ids_,enhance,ids_[1]==ids_[2],detune,t,z_);
	// actual values for z-limits
	if(!computeSpaceLikeLimits(t,x)) return false;
	if(t<tmin) {
	t=-1.0*GeV2;
	return false;
	}
	else
	return true;
	}

	bool SudakovFormFactor::PSVeto(const Energy2 t) {
	// still inside PS, return true if outside
	// check vs overestimated limits
	if (z() < zlimits_.first \|\| z() > zlimits_.second) return true;

	Energy2 m02 = (ids_[0]->id()!=ParticleID::g && ids_[0]->id()!=ParticleID::gamma) ?
	masssquared_[0] : Energy2();

	Energy2 pt2 = QTildeKinematics::pT2_FSR(t,z(),m02,masssquared_[1],masssquared_[2],
	masssquared_[1],masssquared_[2]);

	// if pt2<0 veto
	if (pt2<cutoff_->pT2min()) return true;
	// otherwise calculate pt and return
	pT_ = sqrt(pt2);
	return false;
	}



	ShoKinPtr SudakovFormFactor::generateNextTimeBranching(const Energy startingScale,
	const IdList &ids,
	const RhoDMatrix & rho,
	double enhance,
	double detuning) {
	// First reset the internal kinematics variables that can
	// have been eventually set in the previous call to the method.
	q_ = ZERO;
	z_ = 0.;
	phi_ = 0.;
	// perform initialization
	Energy2 tmax(sqr(startingScale)),tmin;
	initialize(ids,tmin);
	// check max > min
	if(tmax<=tmin) return ShoKinPtr();
	// calculate next value of t using veto algorithm
	Energy2 t(tmax);
	// no shower variations to calculate
	if(ShowerHandler::currentHandler()->showerVariations().empty()){
	// Without variations do the usual Veto algorithm
	// No need for more if-statements in this loop.
	do {
	if(!guessTimeLike(t,tmin,enhance,detuning)) break;
	}
	while(PSVeto(t) \|\|
	SplittingFnVeto(z()(1.-z())t,ids,true,rho,detuning) \|\|
	alphaSVeto(splittingFn()->pTScale() ? sqr(z()(1.-z()))t : z()(1.-z())t));
	}
	else {
	bool alphaRew(true),PSRew(true),SplitRew(true);
	do {
	if(!guessTimeLike(t,tmin,enhance,detuning)) break;
	PSRew=PSVeto(t);
	if (PSRew) continue;
	SplitRew=SplittingFnVeto(z()(1.-z())t,ids,true,rho,detuning);
	alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(z()(1.-z()))t : z()(1.-z())t);
	double factor=alphaSVetoRatio(splittingFn()->pTScale() ? sqr(z()(1.-z()))t : z()(1.-z())t,1.)*
	SplittingFnVetoRatio(z()(1.-z())t,ids,true,rho,detuning);

	tShowerHandlerPtr ch = ShowerHandler::currentHandler();

	if( !(SplitRew \|\| alphaRew) ) {
	//Emission
	q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV;
	if (q_ <= ZERO) break;
	}

	for ( map<string,ShowerVariation>::const_iterator var =
	ch->showerVariations().begin();
	var != ch->showerVariations().end(); ++var ) {
	if ( ( ch->firstInteraction() && var->second.firstInteraction ) \|\|
	( !ch->firstInteraction() && var->second.secondaryInteractions ) ) {

	double newfactor = alphaSVetoRatio(splittingFn()->pTScale() ?
	sqr(z()(1.-z()))t :
	z()(1.-z())t,var->second.renormalizationScaleFactor)
	* SplittingFnVetoRatio(z()(1.-z())t,ids,true,rho,detuning);

	double varied;
	if ( SplitRew \|\| alphaRew ) {
	// No Emission
	varied = (1. - newfactor) / (1. - factor);
	} else {
	// Emission
	varied = newfactor / factor;
	}

	map<string,double>::iterator wi = ch->currentWeights().find(var->first);
	if ( wi != ch->currentWeights().end() )
	wi->second *= varied;
	else {
	assert(false);
	//ch->currentWeights()[var->first] = varied;
	}
	}
	}

	}
	while(PSRew \|\| SplitRew \|\| alphaRew);
	}
	q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV;
	if(q_ < ZERO) return ShoKinPtr();

	// return the ShowerKinematics object
	return new_ptr(FS_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this));
	}

	ShoKinPtr SudakovFormFactor::
	generateNextSpaceBranching(const Energy startingQ,
	const IdList &ids,
	double x,
	const RhoDMatrix & rho,
	double enhance,
	Ptr<BeamParticleData>::transient_const_pointer beam,
	double detuning) {
	// First reset the internal kinematics variables that can
	// have been eventually set in the previous call to the method.
	q_ = ZERO;
	z_ = 0.;
	phi_ = 0.;
	// perform the initialization
	Energy2 tmax(sqr(startingQ)),tmin;
	initialize(ids,tmin);
	// check max > min
	if(tmax<=tmin) return ShoKinPtr();
	// calculate next value of t using veto algorithm
	Energy2 t(tmax),pt2(ZERO);
	// no shower variations
	if(ShowerHandler::currentHandler()->showerVariations().empty()){
	// Without variations do the usual Veto algorithm
	// No need for more if-statements in this loop.
	do {
	if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break;
	pt2 = QTildeKinematics::pT2_ISR(t,z(),masssquared_[2]);
	}
	while(pt2 < cutoff_->pT2min()\|\|
	z() > zlimits_.second\|\|
	SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning)\|\|
	alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())t : (1.-z())t)\|\|
	PDFVeto(t,x,ids[0],ids[1],beam));
	}
	// shower variations
	else
	{
	bool alphaRew(true),PDFRew(true),ptRew(true),zRew(true),SplitRew(true);
	do {
	if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break;
	pt2 = QTildeKinematics::pT2_ISR(t,z(),masssquared_[2]);
	ptRew=pt2 < cutoff_->pT2min();
	zRew=z() > zlimits_.second;
	if (ptRew\|\|zRew) continue;
	SplitRew=SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning);
	alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())t : (1.-z())t);
	PDFRew=PDFVeto(t,x,ids[0],ids[1],beam);
	double factor=PDFVetoRatio(t,x,ids[0],ids[1],beam,1.)*
	alphaSVetoRatio(splittingFn()->pTScale() ? sqr(1.-z())t : (1.-z())t,1.)*
	SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning);

	tShowerHandlerPtr ch = ShowerHandler::currentHandler();

	if( !(PDFRew \|\| SplitRew \|\| alphaRew) ) {
	//Emission
	q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV;
	if (q_ <= ZERO) break;
	}

	for ( map<string,ShowerVariation>::const_iterator var =
	ch->showerVariations().begin();
	var != ch->showerVariations().end(); ++var ) {
	if ( ( ch->firstInteraction() && var->second.firstInteraction ) \|\|
	( !ch->firstInteraction() && var->second.secondaryInteractions ) ) {



	double newfactor = PDFVetoRatio(t,x,ids[0],ids[1],beam,var->second.factorizationScaleFactor)*
	alphaSVetoRatio(splittingFn()->pTScale() ?
	sqr(1.-z())t : (1.-z())t,var->second.renormalizationScaleFactor)
	SplittingFnVetoRatio((1.-z())t/z(),ids,false,rho,detuning);

	double varied;
	if( PDFRew \|\| SplitRew \|\| alphaRew) {
	// No Emission
	varied = (1. - newfactor) / (1. - factor);
	} else {
	// Emission
	varied = newfactor / factor;
	}


	map<string,double>::iterator wi = ch->currentWeights().find(var->first);
	if ( wi != ch->currentWeights().end() )
	wi->second *= varied;
	else {
	assert(false);
	//ch->currentWeights()[var->first] = varied;
	}
	}
	}

	}
	while( PDFRew \|\| SplitRew \|\| alphaRew);
	}
	if(t > ZERO && zlimits_.first < zlimits_.second) q_ = sqrt(t);
	else return ShoKinPtr();

	pT_ = sqrt(pt2);
	// create the ShowerKinematics and return it
	return new_ptr(IS_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this));
	}

	void SudakovFormFactor::initialize(const IdList & ids, Energy2 & tmin) {
	ids_=ids;
	tmin = 4.*cutoff_->pT2min();
	masses_ = cutoff_->virtualMasses(ids);
	masssquared_.clear();
	for(unsigned int ix=0;ix<masses_.size();++ix) {
	masssquared_.push_back(sqr(masses_[ix]));
	if(ix>0) tmin=max(masssquared_[ix],tmin);
	}
	}

	ShoKinPtr SudakovFormFactor::generateNextDecayBranching(const Energy startingScale,
	const Energy stoppingScale,
	const Energy minmass,
	const IdList &ids,
	const RhoDMatrix & rho,
	double enhance,
	double detuning) {
	// First reset the internal kinematics variables that can
	// have been eventually set in the previous call to this method.
	q_ = Constants::MaxEnergy;
	z_ = 0.;
	phi_ = 0.;
	// perform initialisation
	Energy2 tmax(sqr(stoppingScale)),tmin;
	initialize(ids,tmin);
	tmin=sqr(startingScale);
	// check some branching possible
	if(tmax<=tmin) return ShoKinPtr();
	// perform the evolution
	Energy2 t(tmin),pt2(-MeV2);
	do {
	if(!guessDecay(t,tmax,minmass,enhance,detuning)) break;
	pt2 = QTildeKinematics::pT2_Decay(t,z(),masssquared_[0],masssquared_[2]);
	}
	while(SplittingFnVeto((1.-z())*t/z(),ids,true,rho,detuning)\|\|
	alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())t : (1.-z())t ) \|\|
	pt2<cutoff_->pT2min() \|\|
	t*(1.-z())>masssquared_[0]-sqr(minmass));
	if(t > ZERO) {
	q_ = sqrt(t);
	pT_ = sqrt(pt2);
	}
	else return ShoKinPtr();
	phi_ = 0.;
	// create the ShowerKinematics object
	return new_ptr(Decay_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this));
	}

	bool SudakovFormFactor::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass,
	double enhance, double detune) {
	minmass = max(minmass,GeV);
	// previous scale
	Energy2 told = t;
	// overestimated limits on z
	if(tmax<masssquared_[0]) {
	t=-1.0*GeV2;
	return false;
	}
	Energy2 tm2 = tmax-masssquared_[0];
	Energy tm = sqrt(tm2);
	zlimits_ = make_pair(sqr(minmass/masses_[0]),
	1.-sqrt(masssquared_[2]+cutoff_->pT2min()+
	0.25*sqr(masssquared_[2])/tm2)/tm
	+0.5*masssquared_[2]/tm2);
	if(zlimits_.second<zlimits_.first) {
	t=-1.0*GeV2;
	return false;
	}
	// guess values of t and z
	guesstz(told,2,ids_,enhance,ids_[1]==ids_[2],detune,t,z_);
	// actual values for z-limits
	if(t<masssquared_[0]) {
	t=-1.0*GeV2;
	return false;
	}
	tm2 = t-masssquared_[0];
	tm = sqrt(tm2);
	zlimits_ = make_pair(sqr(minmass/masses_[0]),
	1.-sqrt(masssquared_[2]+cutoff_->pT2min()+
	0.25*sqr(masssquared_[2])/tm2)/tm
	+0.5*masssquared_[2]/tm2);
	if(t>tmax\|\|zlimits_.second<zlimits_.first) {
	t=-1.0*GeV2;
	return false;
	}
	else
	return true;
	}

	bool SudakovFormFactor::computeTimeLikeLimits(Energy2 & t) {
	if (t < 1e-20 * GeV2) {
	t=-1.*GeV2;
	return false;
	}
	const Energy2 pT2min = cutoff_->pT2min();
	// special case for gluon radiating
	if(ids_[0]->id()==ParticleID::g\|\|ids_[0]->id()==ParticleID::gamma) {
	// no emission possible
	if(t<16.*(masssquared_[1]+pT2min)) {
	t=-1.*GeV2;
	return false;
	}
	// overestimate of the limits
	zlimits_.first = 0.5(1.-sqrt(1.-4.sqrt((masssquared_[1]+pT2min)/t)));
	zlimits_.second = 1.-zlimits_.first;
	}
	// special case for radiated particle is gluon
	else if(ids_[2]->id()==ParticleID::g\|\|ids_[2]->id()==ParticleID::gamma) {
	zlimits_.first = sqrt((masssquared_[1]+pT2min)/t);
	zlimits_.second = 1.-sqrt((masssquared_[2]+pT2min)/t);
	}
	else if(ids_[1]->id()==ParticleID::g\|\|ids_[1]->id()==ParticleID::gamma) {
	zlimits_.second = sqrt((masssquared_[2]+pT2min)/t);
	zlimits_.first = 1.-sqrt((masssquared_[1]+pT2min)/t);
	}
	else {
	zlimits_.first = (masssquared_[1]+pT2min)/t;
	zlimits_.second = 1.-(masssquared_[2]+pT2min)/t;
	}
	if(zlimits_.first>=zlimits_.second) {
	t=-1.*GeV2;
	return false;
	}
	return true;
	}

	bool SudakovFormFactor::computeSpaceLikeLimits(Energy2 & t, double x) {
	if (t < 1e-20 * GeV2) {
	t=-1.*GeV2;
	return false;
	}
	// compute the limits
	zlimits_.first = x;
	double yy = 1.+0.5*masssquared_[2]/t;
	zlimits_.second = yy - sqrt(sqr(yy)-1.+cutoff_->pT2min()/t);
	// return false if lower>upper
	if(zlimits_.second<zlimits_.first) {
	t=-1.*GeV2;
	return false;
	}
	else
	return true;
	}

	namespace {

	tShowerParticlePtr findCorrelationPartner(ShowerParticle & particle,
	bool forward,
	ShowerInteraction inter) {
	tPPtr child = &particle;
	tShowerParticlePtr mother;
	if(forward) {
	mother = !particle.parents().empty() ?
	dynamic_ptr_cast<tShowerParticlePtr>(particle.parents()[0]) : tShowerParticlePtr();
	}
	else {
	mother = particle.children().size()==2 ?
	dynamic_ptr_cast<tShowerParticlePtr>(&particle) : tShowerParticlePtr();
	}
	tShowerParticlePtr partner;
	while(mother) {
	tPPtr otherChild;
	if(forward) {
	for (unsigned int ix=0;ix<mother->children().size();++ix) {
	if(mother->children()[ix]!=child) {
	otherChild = mother->children()[ix];
	break;
	}
	}
	}
	else {
	otherChild = mother->children()[1];
	}
	tShowerParticlePtr other = dynamic_ptr_cast<tShowerParticlePtr>(otherChild);
	if((inter==ShowerInteraction::QCD && otherChild->dataPtr()->coloured()) \|\|
	(inter==ShowerInteraction::QED && otherChild->dataPtr()->charged())) {
	partner = other;
	break;
	}
	if(forward && !other->isFinalState()) {
	partner = dynamic_ptr_cast<tShowerParticlePtr>(mother);
	break;
	}
	child = mother;
	if(forward) {
	mother = ! mother->parents().empty() ?
	dynamic_ptr_cast<tShowerParticlePtr>(mother->parents()[0]) : tShowerParticlePtr();
	}
	else {
	if(mother->children()[0]->children().size()!=2)
	break;
	tShowerParticlePtr mtemp =
	dynamic_ptr_cast<tShowerParticlePtr>(mother->children()[0]);
	if(!mtemp)
	break;
	else
	mother=mtemp;
	}
	}
	if(!partner) {
	if(forward) {
	partner = dynamic_ptr_cast<tShowerParticlePtr>( child)->partner();
	}
	else {
	if(mother) {
	tShowerParticlePtr parent;
	if(!mother->children().empty()) {
	parent = dynamic_ptr_cast<tShowerParticlePtr>(mother->children()[0]);
	}
	if(!parent) {
	parent = dynamic_ptr_cast<tShowerParticlePtr>(mother);
	}
	partner = parent->partner();
	}
	else {
	partner = dynamic_ptr_cast<tShowerParticlePtr>(&particle)->partner();
	}
	}
	}
	return partner;
	}

	pair<double,double> softPhiMin(double phi0, double phi1, double A, double B, double C, double D) {
	double c01 = cos(phi0 - phi1);
	double s01 = sin(phi0 - phi1);
	double s012(sqr(s01)), c012(sqr(c01));
	double A2(AA), B2(BB), C2(CC), D2(DD);
	if(abs(B/A)<1e-10 && abs(D/C)<1e-10) return make_pair(phi0,phi0+Constants::pi);
	double root = sqr(B2)C2D2sqr(s012) + 2.AB2BC2CDc01s012 + 2.AB2BCD2Dc01*s012
	+ 4.A2B2C2D2c012 - A2B2C2D2s012 - A2B2sqr(D2)s012 - sqr(B2)sqr(C2)s012
	- sqr(B2)C2D2s012 - 4.A2ABCD2Dc01 - 4.AB2BC2CDc01 + sqr(A2)sqr(D2)
	+ 2.A2B2C2D2 + sqr(B2)*sqr(C2);
	if(root<0.) return make_pair(phi0,phi0+Constants::pi);
	root = sqrt(root);
	double denom = (-2.ABCDc01 + A2D2 + B2*C2);
	double denom2 = (-BCc01 + A*D);

	double num = B2CD*s012;
	double y1 = Bs01(-C(num + root) + Ddenom) / denom2;
	double y2 = Bs01(-C(num - root) + Ddenom) / denom2;
	double x1 = -(num + root );
	double x2 = -(num - root );
	if(denom<0.) {
	y1*=-1.;
	y2*=-1.;
	x1*=-1.;
	x2*=-1.;
	}
	return make_pair(atan2(y1,x1) + phi0,atan2(y2,x2) + phi0);
	}

	}

	double SudakovFormFactor::generatePhiForward(ShowerParticle & particle,
	const IdList & ids,
	ShoKinPtr kinematics,
	const RhoDMatrix & rho) {
	// no correlations, return flat phi
	if(! dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->correlations())
	return Constants::twopi*UseRandom::rnd();
	// get the kinematic variables
	double z = kinematics->z();
	Energy2 t = z(1.-z)sqr(kinematics->scale());
	Energy pT = kinematics->pT();
	// if soft correlations
	Energy2 pipj,pik;
	bool canBeSoft[2] = {ids[1]->id()==ParticleID::g \|\| ids[1]->id()==ParticleID::gamma,
	ids[2]->id()==ParticleID::g \|\| ids[2]->id()==ParticleID::gamma };
	array<Energy2,3> pjk;
	array<Energy,3> Ek;
	Energy Ei,Ej;
	Energy2 m12(ZERO),m22(ZERO);
	InvEnergy2 aziMax(ZERO);
	bool softAllowed = dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()&&
	(canBeSoft[0] \|\| canBeSoft[1]);
	if(softAllowed) {
	// find the partner for the soft correlations
	tShowerParticlePtr partner=findCorrelationPartner(particle,true,splittingFn()->interactionType());
	// remember we want the softer gluon
	bool swapOrder = !canBeSoft[1] \|\| (canBeSoft[0] && canBeSoft[1] && z < 0.5);
	double zFact = !swapOrder ? (1.-z) : z;
	// compute the transforms to the shower reference frame
	// first the boost
	Lorentz5Momentum pVect = particle.showerBasis()->pVector();
	Lorentz5Momentum nVect = particle.showerBasis()->nVector();
	Boost beta_bb;
	if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) {
	beta_bb = -(pVect + nVect).boostVector();
	}
	else if(particle.showerBasis()->frame()==ShowerBasis::Rest) {
	beta_bb = -pVect.boostVector();
	}
	else
	assert(false);
	pVect.boost(beta_bb);
	nVect.boost(beta_bb);
	Axis axis;
	if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) {
	axis = pVect.vect().unit();
	}
	else if(particle.showerBasis()->frame()==ShowerBasis::Rest) {
	axis = nVect.vect().unit();
	}
	else
	assert(false);
	// and then the rotation
	LorentzRotation rot;
	if(axis.perp2()>0.) {
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	}
	else if(axis.z()<0.) {
	rot.rotate(Constants::pi,Axis(1.,0.,0.));
	}
	rot.invert();
	pVect *= rot;
	nVect *= rot;
	// shower parameters
	Energy2 pn = pVect*nVect, m2 = pVect.m2();
	double alpha0 = particle.showerParameters().alpha;
	double beta0 = 0.5/alpha0/pn*
	(sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt));
	Lorentz5Momentum qperp0(particle.showerParameters().ptx,
	particle.showerParameters().pty,ZERO,ZERO);
	assert(partner);
	Lorentz5Momentum pj = partner->momentum();
	pj.boost(beta_bb);
	pj *= rot;
	// compute the two phi independent dot products
	pik = 0.5zFact(sqr(alpha0)m2 - sqr(particle.showerParameters().pt) + 2.alpha0beta0pn )
	+0.5*sqr(pT)/zFact;
	Energy2 dot1 = pj*pVect;
	Energy2 dot2 = pj*nVect;
	Energy2 dot3 = pj*qperp0;
	pipj = alpha0dot1+beta0dot2+dot3;
	// compute the constants for the phi dependent dot product
	pjk[0] = zFact(alpha0dot1+dot3-0.5dot2/pn(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
	+0.5sqr(pT)dot2/pn/zFact/alpha0;
	pjk[1] = (pj.x() - dot2/alpha0/pnqperp0.x())pT;
	pjk[2] = (pj.y() - dot2/alpha0/pnqperp0.y())pT;
	m12 = sqr(particle.dataPtr()->mass());
	m22 = sqr(partner->dataPtr()->mass());
	if(swapOrder) {
	pjk[1] *= -1.;
	pjk[2] *= -1.;
	}
	Ek[0] = zFact(alpha0pVect.t()-0.5nVect.t()/pn(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
	+0.5sqr(pT)nVect.t()/pn/zFact/alpha0;
	Ek[1] = -nVect.t()/alpha0/pnqperp0.x()pT;
	Ek[2] = -nVect.t()/alpha0/pnqperp0.y()pT;
	if(swapOrder) {
	Ek[1] *= -1.;
	Ek[2] *= -1.;
	}
	Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2]));
	Ei = alpha0pVect.t()+beta0nVect.t();
	Ej = pj.t();
	double phi0 = atan2(-pjk[2],-pjk[1]);
	if(phi0<0.) phi0 += Constants::twopi;
	double phi1 = atan2(-Ek[2],-Ek[1]);
	if(phi1<0.) phi1 += Constants::twopi;
	double xi_min = pik/Ei/(Ek[0]+mag2), xi_max = pik/Ei/(Ek[0]-mag2), xi_ij = pipj/Ei/Ej;
	if(xi_min>xi_max) swap(xi_min,xi_max);
	if(xi_min>xi_ij) softAllowed = false;
	Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
	if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==1) {
	aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
	}
	else if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==2) {
	double A = (pipjEk[0]- Ejpik)/Ej/sqr(Ej);
	double B = -sqrt(sqr(pipj)*(sqr(Ek[1])+sqr(Ek[2])))/Ej/sqr(Ej);
	double C = pjk[0]/sqr(Ej);
	double D = -sqrt(sqr(pjk[1])+sqr(pjk[2]))/sqr(Ej);
	pair<double,double> minima = softPhiMin(phi0,phi1,A,B,C,D);
	aziMax = 0.5/pik/(Ek[0]-mag2)(Ei-m12(Ek[0]-mag2)/pik + max(Ej(A+Bcos(minima.first -phi1))/(C+D*cos(minima.first -phi0)),
	Ej(A+Bcos(minima.second-phi1))/(C+D*cos(minima.second-phi0))));
	}
	else
	assert(false);
	}
	// if spin correlations
	vector<pair<int,Complex> > wgts;
	if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->spinCorrelations()) {
	// calculate the weights
	wgts = splittingFn()->generatePhiForward(z,t,ids,rho);
	}
	else {
	wgts = {{ {0, 1.} }};
	}
	// generate the azimuthal angle
	double phi,wgt;
	static const Complex ii(0.,1.);
	unsigned int ntry(0);
	double phiMax(0.),wgtMax(0.);
	do {
	phi = Constants::twopi*UseRandom::rnd();
	// first the spin correlations bit (gives 1 if correlations off)
	Complex spinWgt = 0.;
	for(unsigned int ix=0;ix<wgts.size();++ix) {
	if(wgts[ix].first==0)
	spinWgt += wgts[ix].second;
	else
	spinWgt += exp(double(wgts[ix].first)iiphi)*wgts[ix].second;
	}
	wgt = spinWgt.real();
	if(wgt-1.>1e-10) {
	generator()->log() << "Forward spin weight problem " << wgt << " " << wgt-1.
	<< " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n";
	generator()->log() << "Weights \n";
	for(unsigned int ix=0;ix<wgts.size();++ix)
	generator()->log() << wgts[ix].first << " " << wgts[ix].second << "\n";
	}
	// soft correlations bit
	double aziWgt = 1.;
	if(softAllowed) {
	Energy2 dot = pjk[0]+pjk[1]cos(phi)+pjk[2]sin(phi);
	Energy Eg = Ek[0]+Ek[1]cos(phi)+Ek[2]sin(phi);
	if(pipjEg>pikEj) {
	if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==1) {
	aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
	}
	else if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==2) {
	aziWgt = max(ZERO,0.5/pik/Eg(Ei-m12Eg/pik + (pipjEg - Ejpik)/dot)/aziMax);
	}
	if(aziWgt-1.>1e-10\|\|aziWgt<-1e-10) {
	generator()->log() << "Forward soft weight problem " << aziWgt << " " << aziWgt-1.
	<< " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n";
	}
	}
	else {
	aziWgt = 0.;
	}
	}
	wgt *= aziWgt;
	if(wgt>wgtMax) {
	phiMax = phi;
	wgtMax = wgt;
	}
	++ntry;
	}
	while(wgt<UseRandom::rnd()&&ntry<10000);
	if(ntry==10000) {
	generator()->log() << "Too many tries to generate phi in forward evolution\n";
	phi = phiMax;
	}
	// return the azimuthal angle
	return phi;
	}

	double SudakovFormFactor::generatePhiBackward(ShowerParticle & particle,
	const IdList & ids,
	ShoKinPtr kinematics,
	const RhoDMatrix & rho) {
	// no correlations, return flat phi
	if(! dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->correlations())
	return Constants::twopi*UseRandom::rnd();
	// get the kinematic variables
	double z = kinematics->z();
	Energy2 t = (1.-z)*sqr(kinematics->scale())/z;
	Energy pT = kinematics->pT();
	// if soft correlations
	bool softAllowed = dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations() &&
	(ids[2]->id()==ParticleID::g \|\| ids[2]->id()==ParticleID::gamma);
	Energy2 pipj,pik,m12(ZERO),m22(ZERO);
	array<Energy2,3> pjk;
	Energy Ei,Ej,Ek;
	InvEnergy2 aziMax(ZERO);
	if(softAllowed) {
	// find the partner for the soft correlations
	tShowerParticlePtr partner=findCorrelationPartner(particle,false,splittingFn()->interactionType());
	double zFact = (1.-z);
	// compute the transforms to the shower reference frame
	// first the boost
	Lorentz5Momentum pVect = particle.showerBasis()->pVector();
	Lorentz5Momentum nVect = particle.showerBasis()->nVector();
	assert(particle.showerBasis()->frame()==ShowerBasis::BackToBack);
	Boost beta_bb = -(pVect + nVect).boostVector();
	pVect.boost(beta_bb);
	nVect.boost(beta_bb);
	Axis axis = pVect.vect().unit();
	// and then the rotation
	LorentzRotation rot;
	if(axis.perp2()>0.) {
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	}
	else if(axis.z()<0.) {
	rot.rotate(Constants::pi,Axis(1.,0.,0.));
	}
	rot.invert();
	pVect *= rot;
	nVect *= rot;
	// shower parameters
	Energy2 pn = pVect*nVect;
	Energy2 m2 = pVect.m2();
	double alpha0 = particle.x();
	double beta0 = -0.5/alpha0/pnsqr(alpha0)m2;
	Lorentz5Momentum pj = partner->momentum();
	pj.boost(beta_bb);
	pj *= rot;
	double beta2 = 0.5(1.-zFact)(sqr(alpha0zFact/(1.-zFact))m2+sqr(pT))/alpha0/zFact/pn;
	// compute the two phi independent dot products
	Energy2 dot1 = pj*pVect;
	Energy2 dot2 = pj*nVect;
	pipj = alpha0dot1+beta0dot2;
	pik = alpha0(alpha0zFact/(1.-zFact)m2+pn(beta2+zFact/(1.-zFact)*beta0));
	// compute the constants for the phi dependent dot product
	pjk[0] = alpha0zFact/(1.-zFact)dot1+beta2*dot2;
	pjk[1] = pj.x()*pT;
	pjk[2] = pj.y()*pT;
	m12 = ZERO;
	m22 = sqr(partner->dataPtr()->mass());
	Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
	if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==1) {
	aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
	}
	else if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==2) {
	Ek = alpha0zFact/(1.-zFact)pVect.t()+beta2*nVect.t();
	Ei = alpha0pVect.t()+beta0nVect.t();
	Ej = pj.t();
	if(pipjEk> Ejpik) {
	aziMax = 0.5/pik/Ek(Ei-m12Ek/pik + (pipjEk- Ejpik)/(pjk[0]-mag));
	}
	else {
	aziMax = 0.5/pik/Ek(Ei-m12Ek/pik);
	}
	}
	else {
	assert(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==0);
	}
	}
	// if spin correlations
	vector<pair<int,Complex> > wgts;
	if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->spinCorrelations()) {
	// get the weights
	wgts = splittingFn()->generatePhiBackward(z,t,ids,rho);
	}
	else {
	wgts = {{ {0, 1.} }};
	}
	// generate the azimuthal angle
	double phi,wgt;
	static const Complex ii(0.,1.);
	unsigned int ntry(0);
	double phiMax(0.),wgtMax(0.);
	do {
	phi = Constants::twopi*UseRandom::rnd();
	Complex spinWgt = 0.;
	for(unsigned int ix=0;ix<wgts.size();++ix) {
	if(wgts[ix].first==0)
	spinWgt += wgts[ix].second;
	else
	spinWgt += exp(double(wgts[ix].first)iiphi)*wgts[ix].second;
	}
	wgt = spinWgt.real();
	if(wgt-1.>1e-10) {
	generator()->log() << "Backward weight problem " << wgt << " " << wgt-1.
	<< " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << z << " " << phi << "\n";
	generator()->log() << "Weights \n";
	for(unsigned int ix=0;ix<wgts.size();++ix)
	generator()->log() << wgts[ix].first << " " << wgts[ix].second << "\n";
	}
	// soft correlations bit
	double aziWgt = 1.;
	if(softAllowed) {
	Energy2 dot = pjk[0]+pjk[1]cos(phi)+pjk[2]sin(phi);
	if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==1) {
	aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
	}
	else if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==2) {
	aziWgt = max(ZERO,0.5/pik/Ek(Ei-m12Ek/pik + pipjEk/dot - Ejpik/dot)/aziMax);
	}
	if(aziWgt-1.>1e-10\|\|aziWgt<-1e-10) {
	generator()->log() << "Backward soft weight problem " << aziWgt << " " << aziWgt-1.
	<< " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n";
	}
	}
	wgt *= aziWgt;
	if(wgt>wgtMax) {
	phiMax = phi;
	wgtMax = wgt;
	}
	++ntry;
	}
	while(wgt<UseRandom::rnd()&&ntry<10000);
	if(ntry==10000) {
	generator()->log() << "Too many tries to generate phi in backward evolution\n";
	phi = phiMax;
	}
	// return the azimuthal angle
	return phi;
	}

	double SudakovFormFactor::generatePhiDecay(ShowerParticle & particle,
	const IdList & ids,
	ShoKinPtr kinematics,
	const RhoDMatrix &) {
	// only soft correlations in this case
	// no correlations, return flat phi
	if( !(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations() &&
	(ids[2]->id()==ParticleID::g \|\| ids[2]->id()==ParticleID::gamma )))
	return Constants::twopi*UseRandom::rnd();
	// get the kinematic variables
	double z = kinematics->z();
	Energy pT = kinematics->pT();
	// if soft correlations
	// find the partner for the soft correlations
	tShowerParticlePtr partner = findCorrelationPartner(particle,true,splittingFn()->interactionType());
	double zFact(1.-z);
	// compute the transforms to the shower reference frame
	// first the boost
	Lorentz5Momentum pVect = particle.showerBasis()->pVector();
	Lorentz5Momentum nVect = particle.showerBasis()->nVector();
	assert(particle.showerBasis()->frame()==ShowerBasis::Rest);
	Boost beta_bb = -pVect.boostVector();
	pVect.boost(beta_bb);
	nVect.boost(beta_bb);
	Axis axis = nVect.vect().unit();
	// and then the rotation
	LorentzRotation rot;
	if(axis.perp2()>0.) {
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	}
	else if(axis.z()<0.) {
	rot.rotate(Constants::pi,Axis(1.,0.,0.));
	}
	rot.invert();
	pVect *= rot;
	nVect *= rot;
	// shower parameters
	Energy2 pn = pVect*nVect;
	Energy2 m2 = pVect.m2();
	double alpha0 = particle.showerParameters().alpha;
	double beta0 = 0.5/alpha0/pn*
	(sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt));
	Lorentz5Momentum qperp0(particle.showerParameters().ptx,
	particle.showerParameters().pty,ZERO,ZERO);
	Lorentz5Momentum pj = partner->momentum();
	pj.boost(beta_bb);
	pj *= rot;
	// compute the two phi independent dot products
	Energy2 pik = 0.5zFact(sqr(alpha0)m2 - sqr(particle.showerParameters().pt) + 2.alpha0beta0pn )
	+0.5*sqr(pT)/zFact;
	Energy2 dot1 = pj*pVect;
	Energy2 dot2 = pj*nVect;
	Energy2 dot3 = pj*qperp0;
	Energy2 pipj = alpha0dot1+beta0dot2+dot3;
	// compute the constants for the phi dependent dot product
	array<Energy2,3> pjk;
	pjk[0] = zFact(alpha0dot1+dot3-0.5dot2/pn(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
	+0.5sqr(pT)dot2/pn/zFact/alpha0;
	pjk[1] = (pj.x() - dot2/alpha0/pnqperp0.x())pT;
	pjk[2] = (pj.y() - dot2/alpha0/pnqperp0.y())pT;
	Energy2 m12 = sqr(particle.dataPtr()->mass());
	Energy2 m22 = sqr(partner->dataPtr()->mass());
	Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
	InvEnergy2 aziMax;
	array<Energy,3> Ek;
	Energy Ei,Ej;
	if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==1) {
	aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
	}
	else if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==2) {
	Ek[0] = zFact(alpha0pVect.t()+-0.5nVect.t()/pn(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
	+0.5sqr(pT)nVect.t()/pn/zFact/alpha0;
	Ek[1] = -nVect.t()/alpha0/pnqperp0.x()pT;
	Ek[2] = -nVect.t()/alpha0/pnqperp0.y()pT;
	Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2]));
	Ei = alpha0pVect.t()+beta0nVect.t();
	Ej = pj.t();
	aziMax = 0.5/pik/(Ek[0]-mag2)(Ei-m12(Ek[0]-mag2)/pik + pipj(Ek[0]+mag2)/(pjk[0]-mag) - Ejpik/(pjk[0]-mag) );
	}
	else
	assert(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==0);
	// generate the azimuthal angle
	double phi,wgt(0.);
	unsigned int ntry(0);
	double phiMax(0.),wgtMax(0.);
	do {
	phi = Constants::twopi*UseRandom::rnd();
	Energy2 dot = pjk[0]+pjk[1]cos(phi)+pjk[2]sin(phi);
	if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==1) {
	wgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
	}
	else if(dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->softCorrelations()==2) {
	if(qperp0.m2()==ZERO) {
	wgt = 1.;
	}
	else {
	Energy Eg = Ek[0]+Ek[1]cos(phi)+Ek[2]sin(phi);
	wgt = max(ZERO,0.5/pik/Eg(Ei-m12Eg/pik + (pipjEg - Ejpik)/dot)/aziMax);
	}
	}
	if(wgt-1.>1e-10\|\|wgt<-1e-10) {
	generator()->log() << "Decay soft weight problem " << wgt << " " << wgt-1.
	<< " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n";
	}
	if(wgt>wgtMax) {
	phiMax = phi;
	wgtMax = wgt;
	}
	++ntry;
	}
	while(wgt<UseRandom::rnd()&&ntry<10000);
	if(ntry==10000) {
	phi = phiMax;
	generator()->log() << "Too many tries to generate phi\n";
	}
	// return the azimuthal angle
	return phi;
	}


	Energy SudakovFormFactor::calculateScale(double zin, Energy pt, IdList ids,
	unsigned int iopt) {
	Energy2 tmin;
	initialize(ids,tmin);
	// final-state branching
	if(iopt==0) {
	Energy2 scale=(sqr(pt)+masssquared_[1](1.-zin)+masssquared_[2]zin);
	if(ids[0]->id()!=ParticleID::g) scale -= zin(1.-zin)masssquared_[0];
	scale /= sqr(zin*(1-zin));
	return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
	}
	else if(iopt==1) {
	Energy2 scale=(sqr(pt)+zin*masssquared_[2])/sqr(1.-zin);
	return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
	}
	else if(iopt==2) {
	Energy2 scale = (sqr(pt)+zin*masssquared_[2])/sqr(1.-zin)+masssquared_[0];
	return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
	}
	else {
	throw Exception() << "Unknown option in SudakovFormFactor::calculateScale() "
	<< "iopt = " << iopt << Exception::runerror;
	}
	}
	diff --git a/Shower/ShowerAlpha.h b/Shower/ShowerAlpha.h
	--- a/Shower/ShowerAlpha.h
	+++ b/Shower/ShowerAlpha.h
	@@ -1,150 +1,179 @@
	// -- C++ --
	//
	// ShowerAlpha.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_ShowerAlpha_H
	#define HERWIG_ShowerAlpha_H
	//
	// This is the declaration of the ShowerAlpha class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "ShowerAlpha.fh"

	namespace Herwig {

	using namespace ThePEG;

	/** \ingroup Shower
	*
	* This class is the abstract class from which all types of running couplings
	* used in the Showering derive from.
	* The main purpose of this class, and the ones that derive from it, is
	* to allow systematic uncertainties for the initial-state radiation and,
	* independently, the final-state radiation effects, to be evaluated.
	*
	* This is achieved by allowing a multiplicative factor,
	* which is 1.0 for the "central value",
	* for the scale argument, \f$\mu^2\f$, of the running coupling. This
	* factor, \f$f\f$ is given by the scaleFactor() member and the coupling
	* returned by the value() member is such that
	* \f[\alpha(\mu^2)\to \alpha(f\times\mu^2).\f]
	* This scale factor is a parameter which is settable by the user, via the
	* interface.
	* Although, of course, it is not clear my how much we should scale
	* in order to get a \f$1\sigma\f$ systematic error (but factors:
	* 1/2 and 2 are quite common), this method provides a double side error,
	* and it appears more sensible than the rough and one-sided evaluation
	* obtained
	* via turning off the I.S.R. and/or F.S.R. (possibilities which are,
	* anyway, provided by Herwig).
	*
	* @see \ref ShowerAlphaInterfaces "The interfaces"
	* defined for ShowerAlpha.
	*/
	class ShowerAlpha: public Interfaced {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	ShowerAlpha() : _scaleFactor( 1.0 ) {}
	//@}

	public:

	/**
	* Methods to return the coupling and the scaleFactor
	*/
	//@{
	/**
	* Pure virtual method that is supposed to return the
	* running alpha value evaluated at the input scale.
	* @param scale The scale
	* @return The coupling
	*/
	virtual double value(const Energy2 scale) const = 0;

	/**
	* Virtual method, which
	* should be overridden in a derived class to provide an
	* overestimate approximation of the alpha value.
	*/
	virtual double overestimateValue() const = 0;

	/**
	* Virtual method which returns the ratio of the running alpha
	* value at the input scale to the overestimated value.
	* @param scale The scale
	* @return The ratio
	*/
	virtual double ratio(const Energy2 scale,double factor=1.) const = 0;

	/**
	* It returns the factor that multiplies the
	* scale argument, \f$\mu^2\f$, of the running coupling.
	* This is supposed to be <I>1.0</I> in normal conditions (central values)
	* whereas different values can be useful for systematics evaluation
	* for Initial State radiation or Final State radiation effects.
	*/
	double scaleFactor() const {return _scaleFactor;}

	/**
	+ * Virtual method that is supposed to return the
	+ * running alpha value evaluated at the input scale.
	+ * @param scale The scale
	+ * @return The coupling
	+ */
	+ virtual inline double showerValue(const Energy2 scale) const {
	+ return value(scale);
	+ }
	+
	+ /**
	+ * Virtual method, which
	+ * should be overridden in a derived class to provide an
	+ * overestimate approximation of the alpha value.
	+ */
	+ virtual inline double showerOverestimateValue() const {
	+ return overestimateValue();
	+ }
	+
	+ /**
	+ * Virtual method which returns the ratio of the running alpha
	+ * value at the input scale to the overestimated value.
	+ * @param scale The scale
	+ * @return The ratio
	+ */
	+ virtual inline double showerRatio(const Energy2 scale,double factor=1.) const {
	+ return ratio(scale,factor);
	+ }
	+
	+ /**
	* Initialize this coupling.
	*/
	virtual void initialize () {}
	//@}

	public:

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

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

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

	private:

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

	private:

	/**
	* The scale factor
	*/
	double _scaleFactor;

	};

	}

	#endif /* HERWIG_ShowerAlpha_H */
	diff --git a/Shower/ShowerHandler.cc b/Shower/ShowerHandler.cc
	--- a/Shower/ShowerHandler.cc
	+++ b/Shower/ShowerHandler.cc
	@@ -1,1132 +1,1147 @@
	// -- C++ --
	//
	// ShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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 ShowerHandler class.
	//

	#include "ShowerHandler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Command.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/PDF/PartonBinInstance.h"
	#include "Herwig/PDT/StandardMatchers.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Utilities/StringUtils.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig/Utilities/EnumParticles.h"
	#include "Herwig/PDF/MPIPDF.h"
	#include "Herwig/PDF/MinBiasPDF.h"
	#include "ThePEG/Handlers/EventHandler.h"
	#include "Herwig/PDF/HwRemDecayer.h"
	#include <cassert>
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "Herwig/Decay/DecayIntegrator.h"
	#include "Herwig/Decay/PhaseSpaceMode.h"

	using namespace Herwig;

	DescribeClass<ShowerHandler,CascadeHandler>
	describeShowerHandler ("Herwig::ShowerHandler","HwShower.so");

	ShowerHandler::~ShowerHandler() {}

	tShowerHandlerPtr ShowerHandler::currentHandler_ = tShowerHandlerPtr();

	void ShowerHandler::doinit() {
	CascadeHandler::doinit();
	// copy particles to decay before showering from input vector to the
	// set used in the simulation
	if ( particlesDecayInShower_.empty() ) {
	for(unsigned int ix=0;ix<inputparticlesDecayInShower_.size();++ix)
	particlesDecayInShower_.insert(abs(inputparticlesDecayInShower_[ix]));
	}
	if ( profileScales() ) {
	if ( profileScales()->unrestrictedPhasespace() &&
	restrictPhasespace() ) {
	generator()->log()
	<< "ShowerApproximation warning: The scale profile chosen requires an unrestricted phase space,\n"
	<< "however, the phase space was set to be restricted. Will switch to unrestricted phase space.\n"
	<< flush;
	restrictPhasespace_ = false;
	}
	}
	}

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

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

	ShowerHandler::ShowerHandler() :
	maxtry_(10),maxtryMPI_(10),maxtryDP_(10),maxtryDecay_(100),
	factorizationScaleFactor_(1.0),
	renormalizationScaleFactor_(1.0),
	hardScaleFactor_(1.0),
	restrictPhasespace_(true), maxPtIsMuF_(false),
	spinOpt_(1), pdfFreezingScale_(2.5*GeV),
	doFSR_(true), doISR_(true),
	splitHardProcess_(true),
	includeSpaceTime_(false), vMin_(0.1*GeV2),
	reweight_(1.0) {
	inputparticlesDecayInShower_.push_back( 6 ); // top
	inputparticlesDecayInShower_.push_back( 23 ); // Z0
	inputparticlesDecayInShower_.push_back( 24 ); // W+/-
	inputparticlesDecayInShower_.push_back( 25 ); // h0
	}

	void ShowerHandler::doinitrun(){
	CascadeHandler::doinitrun();
	//can't use isMPIOn here, because the EventHandler is not set at that stage
	if(MPIHandler_) {
	MPIHandler_->initialize();
	if(MPIHandler_->softInt())
	remDec_->initSoftInteractions(MPIHandler_->Ptmin(), MPIHandler_->beta());
	}
	}

	void ShowerHandler::dofinish() {
	CascadeHandler::dofinish();
	if(MPIHandler_) MPIHandler_->finalize();
	}

	void ShowerHandler::persistentOutput(PersistentOStream & os) const {
	os << remDec_ << ounit(pdfFreezingScale_,GeV) << maxtry_
	<< maxtryMPI_ << maxtryDP_ << maxtryDecay_
	<< inputparticlesDecayInShower_
	<< particlesDecayInShower_ << MPIHandler_ << PDFA_ << PDFB_
	<< PDFARemnant_ << PDFBRemnant_
	<< includeSpaceTime_ << ounit(vMin_,GeV2)
	<< factorizationScaleFactor_ << renormalizationScaleFactor_
	<< hardScaleFactor_
	<< restrictPhasespace_ << maxPtIsMuF_ << hardScaleProfile_
	<< showerVariations_ << doFSR_ << doISR_ << splitHardProcess_
	- << spinOpt_ << useConstituentMasses_;
	+ << spinOpt_ << useConstituentMasses_ << tagIntermediates_;
	}

	void ShowerHandler::persistentInput(PersistentIStream & is, int) {
	is >> remDec_ >> iunit(pdfFreezingScale_,GeV) >> maxtry_
	>> maxtryMPI_ >> maxtryDP_ >> maxtryDecay_
	>> inputparticlesDecayInShower_
	>> particlesDecayInShower_ >> MPIHandler_ >> PDFA_ >> PDFB_
	>> PDFARemnant_ >> PDFBRemnant_
	>> includeSpaceTime_ >> iunit(vMin_,GeV2)
	>> factorizationScaleFactor_ >> renormalizationScaleFactor_
	>> hardScaleFactor_
	>> restrictPhasespace_ >> maxPtIsMuF_ >> hardScaleProfile_
	>> showerVariations_ >> doFSR_ >> doISR_ >> splitHardProcess_
	- >> spinOpt_ >> useConstituentMasses_;
	+ >> spinOpt_ >> useConstituentMasses_ >> tagIntermediates_;
	}

	void ShowerHandler::Init() {

	static ClassDocumentation<ShowerHandler> documentation
	("Main driver class for the showering.");

	static Reference<ShowerHandler,HwRemDecayer>
	interfaceRemDecayer("RemDecayer",
	"A reference to the Remnant Decayer object",
	&Herwig::ShowerHandler::remDec_,
	false, false, true, false);

	static Parameter<ShowerHandler,Energy> interfacePDFFreezingScale
	("PDFFreezingScale",
	"The PDF freezing scale",
	&ShowerHandler::pdfFreezingScale_, GeV, 2.5GeV, 2.0GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<ShowerHandler,unsigned int> interfaceMaxTry
	("MaxTry",
	"The maximum number of attempts for the main showering loop",
	&ShowerHandler::maxtry_, 10, 1, 100,
	false, false, Interface::limited);

	static Parameter<ShowerHandler,unsigned int> interfaceMaxTryMPI
	("MaxTryMPI",
	"The maximum number of regeneration attempts for an additional scattering",
	&ShowerHandler::maxtryMPI_, 10, 0, 100,
	false, false, Interface::limited);

	static Parameter<ShowerHandler,unsigned int> interfaceMaxTryDP
	("MaxTryDP",
	"The maximum number of regeneration attempts for an additional hard scattering",
	&ShowerHandler::maxtryDP_, 10, 0, 100,
	false, false, Interface::limited);

	static ParVector<ShowerHandler,long> interfaceDecayInShower
	("DecayInShower",
	"PDG codes of the particles to be decayed in the shower",
	&ShowerHandler::inputparticlesDecayInShower_, -1, 0l, -10000000l, 10000000l,
	false, false, Interface::limited);

	static Reference<ShowerHandler,UEBase> interfaceMPIHandler
	("MPIHandler",
	"The object that administers all additional scatterings.",
	&ShowerHandler::MPIHandler_, false, false, true, true);

	static Reference<ShowerHandler,PDFBase> interfacePDFA
	("PDFA",
	"The PDF for beam particle A. Overrides the particle's own PDF setting."
	"By default used for both the shower and forced splitting in the remnant",
	&ShowerHandler::PDFA_, false, false, true, true, false);

	static Reference<ShowerHandler,PDFBase> interfacePDFB
	("PDFB",
	"The PDF for beam particle B. Overrides the particle's own PDF setting."
	"By default used for both the shower and forced splitting in the remnant",
	&ShowerHandler::PDFB_, false, false, true, true, false);

	static Reference<ShowerHandler,PDFBase> interfacePDFARemnant
	("PDFARemnant",
	"The PDF for beam particle A used to generate forced splittings of the remnant."
	" This overrides both the particle's own PDF setting and the value set by PDFA if used.",
	&ShowerHandler::PDFARemnant_, false, false, true, true, false);

	static Reference<ShowerHandler,PDFBase> interfacePDFBRemnant
	("PDFBRemnant",
	"The PDF for beam particle B used to generate forced splittings of the remnant."
	" This overrides both the particle's own PDF setting and the value set by PDFB if used.",
	&ShowerHandler::PDFBRemnant_, false, false, true, true, false);

	static Switch<ShowerHandler,bool> interfaceIncludeSpaceTime
	("IncludeSpaceTime",
	"Whether to include the model for the calculation of space-time distances",
	&ShowerHandler::includeSpaceTime_, false, false, false);
	static SwitchOption interfaceIncludeSpaceTimeYes
	(interfaceIncludeSpaceTime,
	"Yes",
	"Include the model",
	true);
	static SwitchOption interfaceIncludeSpaceTimeNo
	(interfaceIncludeSpaceTime,
	"No",
	"Only include the displacement from the particle-s lifetime for decaying particles",
	false);

	static Parameter<ShowerHandler,Energy2> interfaceMinimumVirtuality
	("MinimumVirtuality",
	"The minimum virtuality for the space-time model",
	&ShowerHandler::vMin_, GeV2, 0.1GeV2, 0.0GeV2, 1000.0*GeV2,
	false, false, Interface::limited);

	static Parameter<ShowerHandler,double> interfaceFactorizationScaleFactor
	("FactorizationScaleFactor",
	"The factorization scale factor.",
	&ShowerHandler::factorizationScaleFactor_, 1.0, 0.0, 0,
	false, false, Interface::lowerlim);

	static Parameter<ShowerHandler,double> interfaceRenormalizationScaleFactor
	("RenormalizationScaleFactor",
	"The renormalization scale factor.",
	&ShowerHandler::renormalizationScaleFactor_, 1.0, 0.0, 0,
	false, false, Interface::lowerlim);

	static Parameter<ShowerHandler,double> interfaceHardScaleFactor
	("HardScaleFactor",
	"The hard scale factor.",
	&ShowerHandler::hardScaleFactor_, 1.0, 0.0, 0,
	false, false, Interface::lowerlim);

	static Parameter<ShowerHandler,unsigned int> interfaceMaxTryDecay
	("MaxTryDecay",
	"The maximum number of attempts to generate a decay",
	&ShowerHandler::maxtryDecay_, 200, 10, 0,
	false, false, Interface::lowerlim);

	static Reference<ShowerHandler,HardScaleProfile> interfaceHardScaleProfile
	("HardScaleProfile",
	"The hard scale profile to use.",
	&ShowerHandler::hardScaleProfile_, false, false, true, true, false);

	static Switch<ShowerHandler,bool> interfaceMaxPtIsMuF
	("MaxPtIsMuF",
	"",
	&ShowerHandler::maxPtIsMuF_, false, false, false);
	static SwitchOption interfaceMaxPtIsMuFYes
	(interfaceMaxPtIsMuF,
	"Yes",
	"",
	true);
	static SwitchOption interfaceMaxPtIsMuFNo
	(interfaceMaxPtIsMuF,
	"No",
	"",
	false);

	static Switch<ShowerHandler,bool> interfaceRestrictPhasespace
	("RestrictPhasespace",
	"Switch on or off phasespace restrictions",
	&ShowerHandler::restrictPhasespace_, true, false, false);
	static SwitchOption interfaceRestrictPhasespaceYes
	(interfaceRestrictPhasespace,
	"Yes",
	"Perform phasespace restrictions",
	true);
	static SwitchOption interfaceRestrictPhasespaceNo
	(interfaceRestrictPhasespace,
	"No",
	"Do not perform phasespace restrictions",
	false);

	static Command<ShowerHandler> interfaceAddVariation
	("AddVariation",
	"Add a shower variation.",
	&ShowerHandler::doAddVariation, false);

	static Switch<ShowerHandler,bool> interfaceDoFSR
	("DoFSR",
	"Switch on or off final state radiation.",
	&ShowerHandler::doFSR_, true, false, false);
	static SwitchOption interfaceDoFSRYes
	(interfaceDoFSR,
	"Yes",
	"Switch on final state radiation.",
	true);
	static SwitchOption interfaceDoFSRNo
	(interfaceDoFSR,
	"No",
	"Switch off final state radiation.",
	false);

	static Switch<ShowerHandler,bool> interfaceDoISR
	("DoISR",
	"Switch on or off initial state radiation.",
	&ShowerHandler::doISR_, true, false, false);
	static SwitchOption interfaceDoISRYes
	(interfaceDoISR,
	"Yes",
	"Switch on initial state radiation.",
	true);
	static SwitchOption interfaceDoISRNo
	(interfaceDoISR,
	"No",
	"Switch off initial state radiation.",
	false);

	static Switch<ShowerHandler,bool> interfaceSplitHardProcess
	("SplitHardProcess",
	"Whether or not to try and split the hard process into production and decay processes",
	&ShowerHandler::splitHardProcess_, true, false, false);
	static SwitchOption interfaceSplitHardProcessYes
	(interfaceSplitHardProcess,
	"Yes",
	"Split the hard process",
	true);
	static SwitchOption interfaceSplitHardProcessNo
	(interfaceSplitHardProcess,
	"No",
	"Don't split the hard process",
	false);

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

	static Switch<ShowerHandler,bool> interfaceUseConstituentMasses
	("UseConstituentMasses",
	"Whether or not to use constituent masses for the reconstruction of the particle after showering.",
	&ShowerHandler::useConstituentMasses_, true, false, false);
	static SwitchOption interfaceUseConstituentMassesYes
	(interfaceUseConstituentMasses,
	"Yes",
	"Use constituent masses.",
	true);
	static SwitchOption interfaceUseConstituentMassesNo
	(interfaceUseConstituentMasses,
	"No",
	"Don't use constituent masses.",
	false);

	+ static Parameter<ShowerHandler,int> interfaceTagIntermediates
	+ ("TagIntermediates",
	+ "Tag particles after shower with the given status code; if zero, no tagging will be performed.",
	+ &ShowerHandler::tagIntermediates_, 0, 0, 0,
	+ false, false, Interface::lowerlim);
	+
	}

	Energy ShowerHandler::hardScale() const {
	assert(false);
	return ZERO;
	}

	void ShowerHandler::cascade() {
	useMe();
	// Initialise the weights in the event object
	// so that any variations are output regardless of
	// whether showering occurs for the given event
	initializeWeights();
	// get the PDF's from ThePEG (if locally overridden use the local versions)
	tcPDFPtr first = PDFA_ ? tcPDFPtr(PDFA_) : firstPDF().pdf();
	tcPDFPtr second = PDFB_ ? tcPDFPtr(PDFB_) : secondPDF().pdf();
	resetPDFs(make_pair(first,second));
	// set the PDFs for the remnant
	if( ! rempdfs_.first)
	rempdfs_.first = PDFARemnant_ ? PDFPtr(PDFARemnant_) : const_ptr_cast<PDFPtr>(first);
	if( ! rempdfs_.second)
	rempdfs_.second = PDFBRemnant_ ? PDFPtr(PDFBRemnant_) : const_ptr_cast<PDFPtr>(second);
	// get the incoming partons
	tPPair incomingPartons =
	eventHandler()->currentCollision()->primarySubProcess()->incoming();
	// and the parton bins
	PBIPair incomingBins =
	make_pair(lastExtractor()->partonBinInstance(incomingPartons.first),
	lastExtractor()->partonBinInstance(incomingPartons.second));
	// and the incoming hadrons
	tPPair incomingHadrons =
	eventHandler()->currentCollision()->incoming();
	remnantDecayer()->setHadronContent(incomingHadrons);
	// check if incoming hadron == incoming parton
	// and get the incoming hadron if exists or parton otherwise
	incoming_ = make_pair(incomingBins.first ?
	incomingBins.first ->particle() : incomingPartons.first,
	incomingBins.second ?
	incomingBins.second->particle() : incomingPartons.second);
	// check the collision is of the beam particles
	// and if not boost collision to the right frame
	// i.e. the hadron-hadron CMF of the collision
	bool btotal(false);
	LorentzRotation rtotal;
	if(incoming_.first != incomingHadrons.first \|\|
	incoming_.second != incomingHadrons.second ) {
	btotal = true;
	boostCollision(false);
	}
	// set the current ShowerHandler
	setCurrentHandler();
	// first shower the hard process
	try {
	SubProPtr sub = eventHandler()->currentCollision()->primarySubProcess();
	incomingPartons = cascade(sub,lastXCombPtr());
	}
	catch(ShowerTriesVeto &veto){
	throw Exception() << "Failed to generate the shower after "
	<< veto.tries
	<< " attempts in ShowerHandler::cascade()"
	<< Exception::eventerror;
	}
	if(showerHardProcessVeto()) throw Veto();
	+ // tag particles as after shower
	+ if ( tagIntermediates_ > 0 ) {
	+ ParticleSet original = newStep()->particles();
	+ for ( auto p : original ) {
	+ //if ( !p->data().coloured() ) continue;
	+ newStep()->copyParticle(p);
	+ p->status(tagIntermediates_);
	+ }
	+ }
	// if a non-hadron collision return (both incoming non-hadronic)
	if( ( !incomingBins.first\|\|
	!isResolvedHadron(incomingBins.first ->particle()))&&
	( !incomingBins.second\|\|
	!isResolvedHadron(incomingBins.second->particle()))) {
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// perform the reweighting for the hard process shower
	combineWeights();
	// unset the current ShowerHandler
	unSetCurrentHandler();
	return;
	}
	// get the remnants for hadronic collision
	pair<tRemPPtr,tRemPPtr> remnants(getRemnants(incomingBins));
	// set the starting scale of the forced splitting to the PDF freezing scale
	remnantDecayer()->initialize(remnants, incoming_, *currentStep(), pdfFreezingScale());
	// do the first forcedSplitting
	try {
	remnantDecayer()->doSplit(incomingPartons, make_pair(rempdfs_.first,rempdfs_.second), true);
	}
	catch (ExtraScatterVeto) {
	throw Exception() << "Remnant extraction failed in "
	<< "ShowerHandler::cascade() from primary interaction. "
	<< "Please check the PDFs you are using and set/unset them if necessary."
	<< Exception::eventerror;
	}
	// perform the reweighting for the hard process shower
	combineWeights();
	// if no MPI return
	if( !isMPIOn() ) {
	remnantDecayer()->finalize();
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// unset the current ShowerHandler
	unSetCurrentHandler();
	return;
	}
	// generate the multiple scatters use modified pdf's now:
	setMPIPDFs();
	// additional "hard" processes
	unsigned int tries(0);
	// This is the loop over additional hard scatters (most of the time
	// only one, but who knows...)
	for(unsigned int i=1; i <= getMPIHandler()->additionalHardProcs(); i++){
	//counter for regeneration
	unsigned int multSecond = 0;
	// generate the additional scatters
	while( multSecond < getMPIHandler()->multiplicity(i) ) {
	// generate the hard scatter
	tStdXCombPtr lastXC = getMPIHandler()->generate(i);
	SubProPtr sub = lastXC->construct();
	// add to the Step
	newStep()->addSubProcess(sub);
	// increment the counters
	tries++;
	multSecond++;
	if(tries == maxtryDP_)
	throw Exception() << "Failed to establish the requested number "
	<< "of additional hard processes. If this error "
	<< "occurs often, your selection of additional "
	<< "scatter is probably unphysical"
	<< Exception::eventerror;
	// Generate the shower. If not possible veto the event
	try {
	incomingPartons = cascade(sub,lastXC);
	}
	catch(ShowerTriesVeto &veto){
	throw Exception() << "Failed to generate the shower of "
	<< "a secondary hard process after "
	<< veto.tries
	<< " attempts in Evolver::showerHardProcess()"
	<< Exception::eventerror;
	}
	try {
	// do the forcedSplitting
	remnantDecayer()->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false);
	}
	catch(ExtraScatterVeto){
	//remove all particles associated with the subprocess
	newStep()->removeParticle(incomingPartons.first);
	newStep()->removeParticle(incomingPartons.second);
	//remove the subprocess from the list
	newStep()->removeSubProcess(sub);
	//regenerate the scattering
	multSecond--;
	continue;
	}
	// connect with the remnants but don't set Remnant colour,
	// because that causes problems due to the multiple colour lines.
	if ( !remnants.first ->extract(incomingPartons.first , false) \|\|
	!remnants.second->extract(incomingPartons.second, false) )
	throw Exception() << "Remnant extraction failed in "
	<< "ShowerHandler::cascade() for additional scatter"
	<< Exception::runerror;
	}
	// perform the reweighting for the additional hard scatter shower
	combineWeights();
	}
	// the underlying event processes
	unsigned int ptveto(1), veto(0);
	unsigned int max(getMPIHandler()->multiplicity());
	for(unsigned int i=0; i<max; i++) {
	// check how often this scattering has been regenerated
	if(veto > maxtryMPI_) break;
	//generate PSpoint
	tStdXCombPtr lastXC = getMPIHandler()->generate();
	SubProPtr sub = lastXC->construct();
	//If Algorithm=1 additional scatters of the signal type
	// with pt > ptmin have to be vetoed
	//with probability 1/(m+1), where m is the number of occurances in this event
	if( getMPIHandler()->Algorithm() == 1 ){
	//get the pT
	Energy pt = sub->outgoing().front()->momentum().perp();
	if(pt > getMPIHandler()->PtForVeto() && UseRandom::rnd() < 1./(ptveto+1) ){
	ptveto++;
	i--;
	continue;
	}
	}
	// add to the SubProcess to the step
	newStep()->addSubProcess(sub);
	// Run the Shower. If not possible veto the scattering
	try {
	incomingPartons = cascade(sub,lastXC);
	}
	// discard this extra scattering, but try the next one
	catch(ShowerTriesVeto) {
	newStep()->removeSubProcess(sub);
	//regenerate the scattering
	veto++;
	i--;
	continue;
	}
	try{
	//do the forcedSplitting
	remnantDecayer()->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false);
	}
	catch (ExtraScatterVeto) {
	//remove all particles associated with the subprocess
	newStep()->removeParticle(incomingPartons.first);
	newStep()->removeParticle(incomingPartons.second);
	//remove the subprocess from the list
	newStep()->removeSubProcess(sub);
	//regenerate the scattering
	veto++;
	i--;
	continue;
	}
	//connect with the remnants but don't set Remnant colour,
	//because that causes problems due to the multiple colour lines.
	if ( !remnants.first ->extract(incomingPartons.first , false) \|\|
	!remnants.second->extract(incomingPartons.second, false) )
	throw Exception() << "Remnant extraction failed in "
	<< "ShowerHandler::cascade() for MPI hard scattering"
	<< Exception::runerror;
	//reset veto counter
	veto = 0;
	// perform the reweighting for the MPI process shower
	combineWeights();
	}
	// finalize the remnants
	remnantDecayer()->finalize(getMPIHandler()->colourDisrupt(),
	getMPIHandler()->softMultiplicity());
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// unset the current ShowerHandler
	unSetCurrentHandler();
	getMPIHandler()->clean();
	resetPDFs(make_pair(first,second));
	}

	void ShowerHandler::initializeWeights() {
	if ( !showerVariations().empty() ) {

	tEventPtr event = eventHandler()->currentEvent();

	for ( map<string,ShowerVariation>::const_iterator var =
	showerVariations().begin();
	var != showerVariations().end(); ++var ) {

	// Check that this is behaving as intended
	//map<string,double>::iterator wi = event->optionalWeights().find(var->first);
	//assert(wi == event->optionalWeights().end() );

	event->optionalWeights()[var->first] = 1.0;
	currentWeights_[var->first] = 1.0;
	}
	}
	reweight_ = 1.0;
	}

	void ShowerHandler::resetWeights() {
	for ( map<string,double>::iterator w = currentWeights_.begin();
	w != currentWeights_.end(); ++w ) {
	w->second = 1.0;
	}
	reweight_ = 1.0;
	}

	void ShowerHandler::combineWeights() {
	tEventPtr event = eventHandler()->currentEvent();
	for ( map<string,double>::const_iterator w =
	currentWeights_.begin(); w != currentWeights_.end(); ++w ) {
	map<string,double>::iterator ew = event->optionalWeights().find(w->first);
	if ( ew != event->optionalWeights().end() )
	ew->second *= w->second;
	else {
	assert(false && "Weight name unknown.");
	//event->optionalWeights()[w->first] = w->second;
	}
	}
	if ( reweight_ != 1.0 ) {
	Ptr<StandardEventHandler>::tptr eh =
	dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(eventHandler());
	if ( !eh ) {
	throw Exception() << "ShowerHandler::combineWeights() : Cross section reweighting "
	<< "through the shower is currently only available with standard "
	<< "event generators" << Exception::runerror;
	}
	eh->reweight(reweight_);
	}
	}

	string ShowerHandler::doAddVariation(string in) {
	if ( in.empty() )
	return "expecting a name and a variation specification";
	string name = StringUtils::car(in);
	ShowerVariation var;
	string res = var.fromInFile(StringUtils::cdr(in));
	if ( res.empty() ) {
	if ( !var.firstInteraction && !var.secondaryInteractions ) {
	// TODO what about decay showers?
	return "variation does not apply to any shower";
	}
	if ( var.renormalizationScaleFactor == 1.0 &&
	var.factorizationScaleFactor == 1.0 ) {
	return "variation does not vary anything";
	}
	/*
	Repository::clog() << "adding a variation with tag '" << name << "' using\nxir = "
	<< var.renormalizationScaleFactor
	<< " xif = "
	<< var.factorizationScaleFactor
	<< "\napplying to:\n"
	<< "first interaction = " << var.firstInteraction << " "
	<< "secondary interactions = " << var.secondaryInteractions << "\n"
	<< flush;
	*/
	showerVariations()[name] = var;
	}
	return res;
	}

	tPPair ShowerHandler::cascade(tSubProPtr, XCPtr) {
	assert(false);
	return tPPair();
	}

	ShowerHandler::RemPair
	ShowerHandler::getRemnants(PBIPair incomingBins) {
	RemPair remnants;
	// first beam particle
	if(incomingBins.first&&!incomingBins.first->remnants().empty()) {
	remnants.first =
	dynamic_ptr_cast<tRemPPtr>(incomingBins.first->remnants()[0] );
	if(remnants.first) {
	ParticleVector children=remnants.first->children();
	for(unsigned int ix=0;ix<children.size();++ix) {
	if(children[ix]->dataPtr()==remnants.first->dataPtr())
	remnants.first = dynamic_ptr_cast<RemPPtr>(children[ix]);
	}
	//remove existing colour lines from the remnants
	if(remnants.first->colourLine())
	remnants.first->colourLine()->removeColoured(remnants.first);
	if(remnants.first->antiColourLine())
	remnants.first->antiColourLine()->removeAntiColoured(remnants.first);
	}
	}
	// seconnd beam particle
	if(incomingBins.second&&!incomingBins. second->remnants().empty()) {
	remnants.second =
	dynamic_ptr_cast<tRemPPtr>(incomingBins.second->remnants()[0] );
	if(remnants.second) {
	ParticleVector children=remnants.second->children();
	for(unsigned int ix=0;ix<children.size();++ix) {
	if(children[ix]->dataPtr()==remnants.second->dataPtr())
	remnants.second = dynamic_ptr_cast<RemPPtr>(children[ix]);
	}
	//remove existing colour lines from the remnants
	if(remnants.second->colourLine())
	remnants.second->colourLine()->removeColoured(remnants.second);
	if(remnants.second->antiColourLine())
	remnants.second->antiColourLine()->removeAntiColoured(remnants.second);
	}
	}
	assert(remnants.first \|\| remnants.second);
	return remnants;
	}

	namespace {

	void addChildren(tPPtr in,set<tPPtr> & particles) {
	particles.insert(in);
	for(unsigned int ix=0;ix<in->children().size();++ix)
	addChildren(in->children()[ix],particles);
	}
	}

	void ShowerHandler::boostCollision(bool boost) {
	// calculate boost from lab to rest
	if(!boost) {
	Lorentz5Momentum ptotal=incoming_.first ->momentum()+incoming_.second->momentum();
	boost_ = LorentzRotation(-ptotal.boostVector());
	Axis axis((boost_*incoming_.first ->momentum()).vect().unit());
	if(axis.perp2()>0.) {
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	boost_.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	}
	}
	// first call performs the boost and second inverse
	// get the particles to be boosted
	set<tPPtr> particles;
	addChildren(incoming_.first,particles);
	addChildren(incoming_.second,particles);
	// apply the boost
	for(set<tPPtr>::const_iterator cit=particles.begin();
	cit!=particles.end();++cit) {
	(*cit)->transform(boost_);
	}
	if(!boost) boost_.invert();
	}

	void ShowerHandler::setMPIPDFs() {

	if ( !mpipdfs_.first ) {
	// first have to check for MinBiasPDF
	tcMinBiasPDFPtr first = dynamic_ptr_cast<tcMinBiasPDFPtr>(firstPDF().pdf());
	if(first)
	mpipdfs_.first = new_ptr(MPIPDF(first->originalPDF()));
	else
	mpipdfs_.first = new_ptr(MPIPDF(firstPDF().pdf()));
	}

	if ( !mpipdfs_.second ) {
	tcMinBiasPDFPtr second = dynamic_ptr_cast<tcMinBiasPDFPtr>(secondPDF().pdf());
	if(second)
	mpipdfs_.second = new_ptr(MPIPDF(second->originalPDF()));
	else
	mpipdfs_.second = new_ptr(MPIPDF(secondPDF().pdf()));
	}

	if( !remmpipdfs_.first ) {
	tcMinBiasPDFPtr first = dynamic_ptr_cast<tcMinBiasPDFPtr>(rempdfs_.first);
	if(first)
	remmpipdfs_.first = new_ptr(MPIPDF(first->originalPDF()));
	else
	remmpipdfs_.first = new_ptr(MPIPDF(rempdfs_.first));
	}

	if( !remmpipdfs_.second ) {
	tcMinBiasPDFPtr second = dynamic_ptr_cast<tcMinBiasPDFPtr>(rempdfs_.second);
	if(second)
	remmpipdfs_.second = new_ptr(MPIPDF(second->originalPDF()));
	else
	remmpipdfs_.second = new_ptr(MPIPDF(rempdfs_.second));
	}
	// reset the PDFs stored in the base class
	resetPDFs(mpipdfs_);

	}

	bool ShowerHandler::isResolvedHadron(tPPtr particle) {
	if(!HadronMatcher::Check(particle->data())) return false;
	for(unsigned int ix=0;ix<particle->children().size();++ix) {
	if(particle->children()[ix]->id()==ParticleID::Remnant) return true;
	}
	return false;
	}

	namespace {

	bool decayProduct(tSubProPtr subProcess,
	tPPtr particle) {
	// must be time-like and not incoming
	if(particle->momentum().m2()<=ZERO\|\|
	particle == subProcess->incoming().first\|\|
	particle == subProcess->incoming().second) return false;
	// if only 1 outgoing and this is it
	if(subProcess->outgoing().size()==1 &&
	subProcess->outgoing()[0]==particle) return true;
	// must not be the s-channel intermediate otherwise
	if(find(subProcess->incoming().first->children().begin(),
	subProcess->incoming().first->children().end(),particle)!=
	subProcess->incoming().first->children().end()&&
	find(subProcess->incoming().second->children().begin(),
	subProcess->incoming().second->children().end(),particle)!=
	subProcess->incoming().second->children().end()&&
	subProcess->incoming().first ->children().size()==1&&
	subProcess->incoming().second->children().size()==1)
	return false;
	// if non-coloured this is enough
	if(!particle->dataPtr()->coloured()) return true;
	// if coloured must be unstable
	if(particle->dataPtr()->stable()) return false;
	// must not have same particle type as a child
	int id = particle->id();
	for(unsigned int ix=0;ix<particle->children().size();++ix)
	if(particle->children()[ix]->id()==id) return false;
	// otherwise its a decaying particle
	return true;
	}

	PPtr findParent(PPtr original, bool & isHard,
	set<PPtr> outgoingset,
	tSubProPtr subProcess) {
	PPtr parent=original;
	isHard \|=(outgoingset.find(original) != outgoingset.end());
	if(!original->parents().empty()) {
	PPtr orig=original->parents()[0];
	if(decayProduct(subProcess,orig))
	parent=findParent(orig,isHard,outgoingset,subProcess);
	}
	return parent;
	}
	}

	void ShowerHandler::findDecayProducts(PPtr in,PerturbativeProcessPtr hard,
	DecayProcessMap & decay) const {
	ParticleVector children=in->children();
	for(ParticleVector::const_iterator it=children.begin(); it!=children.end();++it) {
	// if decayed or should be decayed in shower make the PerturbaitveProcess
	bool radiates = false;
	if(!(**it).children().empty()) {
	// remove d,u,s,c,b quarks and leptons other than on-shell taus
	if( StandardQCDPartonMatcher::Check((**it).id()) \|\|
	( LeptonMatcher::Check((it).id()) && !(abs((it).id())==ParticleID::tauminus &&
	abs((it).mass()-(it).dataPtr()->mass())<MeV))) {
	radiates = true;
	}
	else {
	bool foundParticle(false),foundGauge(false);
	for(unsigned int iy=0;iy<(**it).children().size();++iy) {
	if((it).children()[iy]->id()==(it).id()) {
	foundParticle = true;
	}
	else if((**it).children()[iy]->id()==ParticleID::g \|\|
	(**it).children()[iy]->id()==ParticleID::gamma) {
	foundGauge = true;
	}
	}
	radiates = foundParticle && foundGauge;
	}
	}
	if(radiates) {
	findDecayProducts(*it,hard,decay);
	}
	else if(!(**it).children().empty()\|\|
	(decaysInShower((it).id())&&!(it).dataPtr()->stable())) {
	createDecayProcess(*it,hard,decay);
	}
	else {
	hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
	}
	}
	}

	void ShowerHandler::splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard,
	DecayProcessMap & decay) const {
	// temporary storage of the particles
	set<PPtr> hardParticles;
	// tagged particles in a set
	set<PPtr> outgoingset(tagged.begin(),tagged.end());
	bool isHard=false;
	// loop over the tagged particles
	for (tParticleVector::const_iterator taggedP = tagged.begin();
	taggedP != tagged.end(); ++taggedP) {
	// skip remnants
	if (eventHandler()->currentCollision()&&
	eventHandler()->currentCollision()->isRemnant(*taggedP)) continue;
	// find the parent and whether its a decaying particle
	bool isDecayProd=false;
	// check if hard
	isHard \|=(outgoingset.find(*taggedP) != outgoingset.end());
	if(splitHardProcess_) {
	tPPtr parent = *taggedP;
	// check if from s channel decaying colourless particle
	while(parent&&!parent->parents().empty()&&!isDecayProd) {
	parent = parent->parents()[0];
	if(parent == subProcess_->incoming().first \|\|
	parent == subProcess_->incoming().second ) break;
	isDecayProd = decayProduct(subProcess_,parent);
	}
	if (isDecayProd)
	hardParticles.insert(findParent(parent,isHard,outgoingset,subProcess_));
	}
	if (!isDecayProd)
	hardParticles.insert(*taggedP);
	}
	// there must be something to shower
	if(hardParticles.empty())
	throw Exception() << "No particles to shower in "
	<< "ShowerHandler::splitHardProcess()"
	<< Exception::eventerror;
	// must be a hard process
	if(!isHard)
	throw Exception() << "Starting on decay not yet implemented in "
	<< "ShowerHandler::splitHardProcess()"
	<< Exception::runerror;
	// create the hard process
	hard = new_ptr(PerturbativeProcess());
	// incoming particles
	hard->incoming().push_back(make_pair(subProcess_->incoming().first ,PerturbativeProcessPtr()));
	hard->incoming().push_back(make_pair(subProcess_->incoming().second,PerturbativeProcessPtr()));
	// outgoing particles
	for(set<PPtr>::const_iterator it=hardParticles.begin();it!=hardParticles.end();++it) {
	// if decayed or should be decayed in shower make the tree
	PPtr orig = *it;
	bool radiates = false;
	if(!orig->children().empty()) {
	// remove d,u,s,c,b quarks and leptons other than on-shell taus
	if( StandardQCDPartonMatcher::Check(orig->id()) \|\|
	( LeptonMatcher::Check(orig->id()) &&
	!(abs(orig->id())==ParticleID::tauminus && abs(orig->mass()-orig->dataPtr()->mass())<MeV))) {
	radiates = true;
	}
	else {
	bool foundParticle(false),foundGauge(false);
	for(unsigned int iy=0;iy<orig->children().size();++iy) {
	if(orig->children()[iy]->id()==orig->id()) {
	foundParticle = true;
	}
	else if(orig->children()[iy]->id()==ParticleID::g \|\|
	orig->children()[iy]->id()==ParticleID::gamma) {
	foundGauge = true;
	}
	}
	radiates = foundParticle && foundGauge;
	}
	}
	if(radiates) {
	findDecayProducts(orig,hard,decay);
	}
	else if(!(**it).children().empty()\|\|
	(decaysInShower((it).id())&&!(it).dataPtr()->stable())) {
	createDecayProcess(*it,hard,decay);
	}
	else {
	hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
	}
	}
	}

	void ShowerHandler::createDecayProcess(PPtr in,PerturbativeProcessPtr hard, DecayProcessMap & decay) const {
	// there must be an incoming particle
	assert(in);
	// create the new process and connect with the parent
	PerturbativeProcessPtr newDecay=new_ptr(PerturbativeProcess());
	newDecay->incoming().push_back(make_pair(in,hard));
	Energy width=in->dataPtr()->generateWidth(in->mass());
	decay.insert(make_pair(width,newDecay));
	hard->outgoing().push_back(make_pair(in,newDecay));
	// we need to deal with the decay products if decayed
	ParticleVector children = in->children();
	if(!children.empty()) {
	for(ParticleVector::const_iterator it = children.begin();
	it!= children.end(); ++it) {
	// if decayed or should be decayed in shower make the tree
	in->abandonChild(*it);
	bool radiates = false;
	if(!(**it).children().empty()) {
	if(StandardQCDPartonMatcher::Check((**it).id())\|\|
	(LeptonMatcher::Check((it).id())&& !(abs((it).id())==ParticleID::tauminus &&
	abs((it).mass()-(it).dataPtr()->mass())<MeV))) {
	radiates = true;
	}
	else {
	bool foundParticle(false),foundGauge(false);
	for(unsigned int iy=0;iy<(**it).children().size();++iy) {
	if((it).children()[iy]->id()==(it).id()) {
	foundParticle = true;
	}
	else if((**it).children()[iy]->id()==ParticleID::g \|\|
	(**it).children()[iy]->id()==ParticleID::gamma) {
	foundGauge = true;
	}
	}
	radiates = foundParticle && foundGauge;
	}
	// finally assume all non-decaying particles are in this class
	// pr 27/11/15 not sure about this bit
	// if(!radiates) {
	// radiates = !decaysInShower((**it).id());
	// }
	}
	if(radiates) {
	findDecayProducts(*it,newDecay,decay);
	}
	else if(!(**it).children().empty()\|\|
	(decaysInShower((it).id())&&!(it).dataPtr()->stable())) {
	createDecayProcess(*it,newDecay,decay);
	}
	else {
	newDecay->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
	}
	}
	}
	}

	tDMPtr ShowerHandler::decay(PerturbativeProcessPtr process,
	DecayProcessMap & decayMap,
	bool radPhotons ) const {
	PPtr parent = process->incoming()[0].first;
	assert(parent);
	if(parent->spinInfo()) parent->spinInfo()->decay(true);
	unsigned int ntry = 0;
	ParticleVector children;
	tDMPtr dm = DMPtr();
	while (true) {
	// exit if fails
	if (++ntry>=maxtryDecay_)
	throw Exception() << "Failed to perform decay in ShowerHandler::decay()"
	<< " after " << maxtryDecay_
	<< " attempts for " << parent->PDGName()
	<< Exception::eventerror;
	// select decay mode
	dm = parent->data().selectMode(*parent);

	if(!dm)
	throw Exception() << "Failed to select decay mode in ShowerHandler::decay()"
	<< "for " << parent->PDGName()
	<< Exception::eventerror;
	if(!dm->decayer())
	throw Exception() << "No Decayer for selected decay mode "
	<< " in ShowerHandler::decay()"
	<< Exception::runerror;
	// start of try block
	try {
	children = dm->decayer()->decay(dm, parent);
	// if no children have another go
	if(children.empty()) continue;

	if(radPhotons){
	// generate radiation in the decay
	tDecayIntegratorPtr hwdec=dynamic_ptr_cast<tDecayIntegratorPtr>(dm->decayer());
	if (hwdec && hwdec->canGeneratePhotons())
	children = hwdec->generatePhotons(*parent,children);
	}

	// 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]);
	}
	// if succeeded break out of loop
	break;
	}
	catch(Veto) {
	}
	}
	assert(!children.empty());
	for(ParticleVector::const_iterator it = children.begin();
	it!= children.end(); ++it) {
	if(!(**it).children().empty()\|\|
	(decaysInShower((it).id())&&!(it).dataPtr()->stable())) {
	createDecayProcess(*it,process,decayMap);
	}
	else {
	process->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
	}
	}
	return dm;
	}


	// Note: The tag must be constructed from an ordered particle container.
	tDMPtr ShowerHandler::findDecayMode(const string & tag) const {
	static map<string,DMPtr> cache;
	map<string,DMPtr>::const_iterator pos = cache.find(tag);

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

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


	/**
	* Operator for the particle ordering
	* @param p1 The first ParticleData object
	* @param p2 The second ParticleData object
	*/

	bool ShowerHandler::ParticleOrdering::operator() (tcPDPtr p1, tcPDPtr p2) const {
	return abs(p1->id()) > abs(p2->id()) \|\|
	( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) \|\|
	( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
	}
	diff --git a/Shower/ShowerHandler.h b/Shower/ShowerHandler.h
	--- a/Shower/ShowerHandler.h
	+++ b/Shower/ShowerHandler.h
	@@ -1,898 +1,908 @@
	// -- C++ --
	//
	// ShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_ShowerHandler_H
	#define HERWIG_ShowerHandler_H
	//
	// This is the declaration of the ShowerHandler class.
	//

	#include "ThePEG/Handlers/EventHandler.h"
	#include "ThePEG/Handlers/CascadeHandler.h"
	#include "ShowerVariation.h"
	#include "Herwig/PDF/HwRemDecayer.fh"
	#include "ThePEG/EventRecord/RemnantParticle.fh"
	#include "UEBase.h"
	#include "PerturbativeProcess.h"
	#include "Herwig/MatrixElement/Matchbox/Matching/HardScaleProfile.h"
	#include "ShowerHandler.fh"

	namespace Herwig {
	using namespace ThePEG;

	/** \ingroup Shower
	*
	* This class is the main driver of the shower: it is responsible for
	* the proper handling of all other specific collaborating classes
	* and for the storing of the produced particles in the event record.
	*
	* @see \ref ShowerHandlerInterfaces "The interfaces"
	*
	* @see ThePEG::CascadeHandler
	* @see MPIHandler
	* @see HwRemDecayer
	*/
	class ShowerHandler: public CascadeHandler {

	public:

	/**
	* Typedef for a pair of ThePEG::RemnantParticle pointers.
	*/
	typedef pair<tRemPPtr, tRemPPtr> RemPair;

	public:

	/**
	* Default constructor
	*/
	ShowerHandler();

	/**
	* Destructor
	*/
	virtual ~ShowerHandler();

	public:

	/**
	* The main method which manages the multiple interactions and starts
	* the shower by calling cascade(sub, lastXC).
	*/
	virtual void cascade();

	/**
	* pointer to "this", the current ShowerHandler.
	*/
	static const tShowerHandlerPtr currentHandler() {
	assert(currentHandler_);
	return currentHandler_;
	}

	/**
	* pointer to "this", the current ShowerHandler.
	*/
	static bool currentHandlerIsSet() {
	return currentHandler_;
	}

	public:

	/**
	* Hook to allow vetoing of event after showering hard sub-process
	* as in e.g. MLM merging.
	*/
	virtual bool showerHardProcessVeto() const { return false; }

	/**
	* Return true, if this cascade handler will perform reshuffling from hard
	* process masses.
	*/
	virtual bool isReshuffling() const { return true; }

	/**
	* Return true, if this cascade handler will put the final state
	* particles to their constituent mass. If false the nominal mass is used.
	*/
	virtual bool retConstituentMasses() const { return useConstituentMasses_; }



	/**
	* Return true, if the shower handler can generate a truncated
	* shower for POWHEG style events generated using Matchbox
	*/
	virtual bool canHandleMatchboxTrunc() const { return false; }

	/**
	* Get the PDF freezing scale
	*/
	Energy pdfFreezingScale() const { return pdfFreezingScale_; }

	/**
	* Get the local PDFs.
	*/
	PDFPtr getPDFA() const {return PDFA_;}

	/**
	* Get the local PDFs.
	*/
	PDFPtr getPDFB() const {return PDFB_;}

	/**
	* Return true if currently the primary subprocess is showered.
	*/
	bool firstInteraction() const {
	if (!eventHandler()->currentCollision())return true;
	return ( subProcess_ ==
	eventHandler()->currentCollision()->primarySubProcess() );
	}

	/**
	* Return the remnant decayer.
	*/
	tHwRemDecPtr remnantDecayer() const { return remDec_; }

	/**
	* Split the hard process into production and decays
	* @param tagged The tagged particles from the StepHandler
	* @param hard The hard perturbative process
	* @param decay The decay particles
	*/
	void splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard,
	DecayProcessMap & decay) const;

	/**
	* Information if the Showerhandler splits the hard process.
	*/
	bool doesSplitHardProcess()const {return splitHardProcess_;}

	/**
	* Decay a particle.
	* radPhotons switches the generation of photon
	* radiation on/off.
	* Required for Dipole Shower but not QTilde Shower.
	*/
	tDMPtr decay(PerturbativeProcessPtr,
	DecayProcessMap & decay,
	bool radPhotons = false) const;


	/**
	* Cached lookup of decay modes.
	* Generator::findDecayMode() is not efficient.
	*/
	tDMPtr findDecayMode(const string & tag) const;


	/**
	* A struct to order the particles in the same way as in the DecayMode's
	*/

	struct ParticleOrdering {

	bool operator() (tcPDPtr p1, tcPDPtr p2) const;

	};


	/**
	* A container for ordered particles required
	* for constructing tags for decay mode lookup.
	*/
	typedef multiset<tcPDPtr,ParticleOrdering> OrderedParticles;


	public:

	/**
	* @name Switches for initial- and final-state radiation
	*/
	//@{
	/**
	* Switch for any radiation
	*/
	bool doRadiation() const {return doFSR_ \|\| doISR_;}

	/**
	* Switch on or off final state radiation.
	*/
	bool doFSR() const { return doFSR_;}

	/**
	* Switch on or off initial state radiation.
	*/
	bool doISR() const { return doISR_;}
	//@}

	public:

	/**
	* @name Switches for scales
	*/
	//@{
	/**
	* Return true if maximum pt should be deduced from the factorization scale
	*/
	bool hardScaleIsMuF() const { return maxPtIsMuF_; }

	/**
	* The factorization scale factor.
	*/
	double factorizationScaleFactor() const {
	return factorizationScaleFactor_;
	}

	/**
	* The renormalization scale factor.
	*/
	double renFac() const {
	return renormalizationScaleFactor_;
	}

	/**
	* The factorization scale factor.
	*/
	double facFac() const {
	return factorizationScaleFactor_;
	}

	/**
	* The renormalization scale factor.
	*/
	double renormalizationScaleFactor() const {
	return renormalizationScaleFactor_;
	}

	/**
	* The scale factor for the hard scale
	*/
	double hardScaleFactor() const {
	return hardScaleFactor_;
	}
	/**
	* Return true, if the phase space restrictions of the dipole shower should
	* be applied.
	*/
	bool restrictPhasespace() const { return restrictPhasespace_; }

	/**
	* Return profile scales
	*/
	Ptr<HardScaleProfile>::tptr profileScales() const { return hardScaleProfile_; }

	/**
	* Return the relevant hard scale to be used in the profile scales
	*/
	virtual Energy hardScale() const;

	/**
	* Return information about shower phase space choices
	*/
	virtual int showerPhaseSpaceOption() const {
	assert(false && "not implemented in general");
	return -1;
	}
	//@}

	public:

	/**
	* Access the shower variations
	*/
	map<string,ShowerVariation>& showerVariations() {
	return showerVariations_;
	}

	/**
	* Return the shower variations
	*/
	const map<string,ShowerVariation>& showerVariations() const {
	return showerVariations_;
	}

	/**
	* Access the current Weights
	*/
	map<string,double>& currentWeights() {
	return currentWeights_;
	}

	/**
	* Return the current Weights
	*/
	const map<string,double>& currentWeights() const {
	return currentWeights_;
	}

	/**
	* Change the current reweighting factor
	*/
	void reweight(double w) {
	reweight_ = w;
	}

	/**
	* Return the current reweighting factor
	*/
	double reweight() const {
	return reweight_;
	}

	public :

	/**
	* Access to switches for spin correlations
	*/
	//@{
	/**
	* Spin Correlations
	*/
	unsigned int spinCorrelations() const {
	return spinOpt_;
	}

	/**
	* Any correlations
	*/
	virtual bool correlations() const {
	return spinOpt_!=0;
	}
	//@}

	public:

	/**
	* struct that is used to catch exceptions which are thrown
	* due to energy conservation issues of additional scatters
	*/
	struct ExtraScatterVeto {};

	/**
	* struct that is used to catch exceptions which are thrown
	* due to fact that the Shower has been invoked more than
	* a defined threshold on a certain configuration
	*/
	struct ShowerTriesVeto {
	/** variable to store the number of attempts */
	const int tries;

	/** constructor */
	ShowerTriesVeto(int t) : tries(t) {}
	};

	public:

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

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

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

	protected:

	/** @name Functions to perform the cascade
	*/
	//@{
	/**
	* The main method which manages the showering of a subprocess.
	*/
	virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb);

	/**
	* Set up for the cascade
	*/
	void prepareCascade(tSubProPtr sub) {
	current_ = currentStep();
	subProcess_ = sub;
	}

	/**
	* Boost all the particles in the collision so that the collision always occurs
	* in the rest frame with the incoming particles along the z axis
	*/
	void boostCollision(bool boost);
	//@}

	protected:

	/**
	* Set/unset the current shower handler
	*/
	//@{
	/**
	* Set the current handler
	*/
	void setCurrentHandler() {
	currentHandler_ = tShowerHandlerPtr(this);
	}

	/**
	* Unset the current handler
	*/
	void unSetCurrentHandler() {
	currentHandler_ = tShowerHandlerPtr();
	}
	//@}

	protected:

	/**
	* @name Members relating to the underlying event and MPI
	*/
	//@{
	/**
	* Return true if multiple parton interactions are switched on
	* and can be used for this beam setup.
	*/
	bool isMPIOn() const {
	return MPIHandler_ && MPIHandler_->beamOK();
	}

	/**
	* Access function for the MPIHandler, it should only be called after
	* checking with isMPIOn.
	*/
	tUEBasePtr getMPIHandler() const {
	assert(MPIHandler_);
	return MPIHandler_;
	}

	/**
	* Is a beam particle where hadronic structure is resolved
	*/
	bool isResolvedHadron(tPPtr);

	/**
	* Get the remnants from the ThePEG::PartonBinInstance es and
	* do some checks.
	*/
	RemPair getRemnants(PBIPair incbins);

	/**
	* Reset the PDF's after the hard collision has been showered
	*/
	void setMPIPDFs();
	//@}

	public:

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

	protected:

	/**
	* Members to handle splitting up of hard process and decays
	*/
	//@{
	/**
	* Find decay products from the hard process and create decay processes
	* @param parent The parent particle
	* @param hard The hard process
	* @param decay The decay processes
	*/
	void findDecayProducts(PPtr parent, PerturbativeProcessPtr hard, DecayProcessMap & decay) const;

	/**
	* Find decay products from the hard process and create decay processes
	* @param parent The parent particle
	* @param hard The parent hard process
	* @param decay The decay processes
	*/
	void createDecayProcess(PPtr parent,PerturbativeProcessPtr hard, DecayProcessMap & decay) const;
	//@}

	/**
	* @name Functions to return information relevant to the process being showered
	*/
	//@{
	/**
	* Return the currently used SubProcess.
	*/
	tSubProPtr currentSubProcess() const {
	assert(subProcess_);
	return subProcess_;
	}

	/**
	* Access to the incoming beam particles
	*/
	tPPair incomingBeams() const {
	return incoming_;
	}
	//@}

	protected:

	/**
	* Weight handling for shower variations
	*/
	//@
	/**
	* Combine the variation weights which have been encountered
	*/
	void combineWeights();

	/**
	* Initialise the weights in currentEvent()
	*/
	void initializeWeights();

	/**
	* Reset the current weights
	*/
	void resetWeights();
	//@}

	protected:

	/**
	* Return the maximum number of attempts for showering
	* a given subprocess.
	*/
	unsigned int maxtry() const { return maxtry_; }

	protected:

	/**
	* Parameters for the space-time model
	*/
	//@{
	/**
	* Whether or not to include spa-cetime distances in the shower
	*/
	bool includeSpaceTime() const {return includeSpaceTime_;}

	/**
	* The minimum virtuality for the space-time model
	*/
	Energy2 vMin() const {return vMin_;}
	//@}

	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;
	//@}

	protected:

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

	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();

	/**
	* Finalize this object. Called in the run phase just after a
	* run has ended. Used eg. to write out statistics.
	*/
	virtual void dofinish();
	//@}

	private:

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

	private:

	/**
	* pointer to "this", the current ShowerHandler.
	*/
	static tShowerHandlerPtr currentHandler_;

	/**
	* a MPIHandler to administer the creation of several (semihard)
	* partonic interactions.
	*/
	UEBasePtr MPIHandler_;

	/**
	* Pointer to the HwRemDecayer
	*/
	HwRemDecPtr remDec_;

	private:

	/**
	* Maximum tries for various stages of the showering process
	*/
	//@{
	/**
	* Maximum number of attempts for the
	* main showering loop
	*/
	unsigned int maxtry_;

	/**
	* Maximum number of attempts for the regeneration of an additional
	* scattering, before the number of scatters is reduced.
	*/
	unsigned int maxtryMPI_;

	/**
	* Maximum number of attempts for the regeneration of an additional
	* hard scattering, before this event is vetoed.
	*/
	unsigned int maxtryDP_;

	/**
	* Maximum number of attempts to generate a decay
	*/
	unsigned int maxtryDecay_;
	//@}

	private:

	/**
	* Factors for the various scales
	*/
	//@{
	/**
	* The factorization scale factor.
	*/
	double factorizationScaleFactor_;

	/**
	* The renormalization scale factor.
	*/
	double renormalizationScaleFactor_;

	/**
	* The scale factor for the hard scale
	*/
	double hardScaleFactor_;

	/**
	* True, if the phase space restrictions of the dipole shower should
	* be applied.
	*/
	bool restrictPhasespace_;

	/**
	* True if maximum pt should be deduced from the factorization scale
	*/
	bool maxPtIsMuF_;

	/**
	* The profile scales
	*/
	Ptr<HardScaleProfile>::ptr hardScaleProfile_;
	//@}

	/**
	* Option to include spin correlations
	*/
	unsigned int spinOpt_;

	private:

	/**
	* Storage of information about the current event
	*/
	//@{
	/**
	* The incoming beam particles for the current collision
	*/
	tPPair incoming_;

	/**
	* Boost to get back to the lab
	*/
	LorentzRotation boost_;

	/**
	* Const pointer to the currently handeled ThePEG::SubProcess
	*/
	tSubProPtr subProcess_;

	/**
	* Const pointer to the current step
	*/
	tcStepPtr current_;
	//@}

	private:

	/**
	* PDFs to be used for the various stages and related parameters
	*/
	//@{
	/**
	* The PDF freezing scale
	*/
	Energy pdfFreezingScale_;

	/**
	* PDFs to be used for the various stages and related parameters
	*/
	//@{
	/**
	* The PDF for beam particle A. Overrides the particle's own PDF setting.
	*/
	PDFPtr PDFA_;

	/**
	* The PDF for beam particle B. Overrides the particle's own PDF setting.
	*/
	PDFPtr PDFB_;

	/**
	* The PDF for beam particle A for remnant splitting. Overrides the particle's own PDF setting.
	*/
	PDFPtr PDFARemnant_;

	/**
	* The PDF for beam particle B for remnant splitting. Overrides the particle's own PDF setting.
	*/
	PDFPtr PDFBRemnant_;

	/**
	* The MPI PDF's to be used for secondary scatters.
	*/
	pair <PDFPtr, PDFPtr> mpipdfs_;

	/**
	* The MPI PDF's to be used for secondary scatters.
	*/
	pair <PDFPtr, PDFPtr> rempdfs_;

	/**
	* The MPI PDF's to be used for secondary scatters.
	*/
	pair <PDFPtr, PDFPtr> remmpipdfs_;
	//@}

	private:

	/**
	* @name Parameters for initial- and final-state radiation
	*/
	//@{
	/**
	* Switch on or off final state radiation.
	*/
	bool doFSR_;

	/**
	* Switch on or off initial state radiation.
	*/
	bool doISR_;
	//@}

	private:

	/**
	* @name Parameters for particle decays
	*/
	//@{
	/**
	* Whether or not to split into hard and decay trees
	*/
	bool splitHardProcess_;

	/**
	* PDG codes of the particles which decay during showering
	* this is fast storage for use during running
	*/
	set<long> particlesDecayInShower_;

	/**
	* PDG codes of the particles which decay during showering
	* this is a vector that is interfaced so they can be changed
	*/
	vector<long> inputparticlesDecayInShower_;
	//@}

	private:

	/**
	* Parameters for the space-time model
	*/
	//@{
	/**
	- * Whether or not to include spa-cetime distances in the shower
	+ * Whether or not to include space-time distances in the shower
	*/
	bool includeSpaceTime_;

	/**
	* The minimum virtuality for the space-time model
	*/
	Energy2 vMin_;
	//@}


	private:

	/**
	* Parameters for the constituent mass treatment.
	*/
	- //@{
	- // True if shower should return constituent masses.
	- bool useConstituentMasses_=true;
	+ //@{
	+
	+ /**
	+ * True if shower should return constituent masses.
	+ */
	+ bool useConstituentMasses_ = true;
	+
	+ /**
	+ * Status code to tag intermediate state as after the showering; if zero no tagging will be performed
	+ */
	+ int tagIntermediates_ = 0;
	+
	//@}
	+
	private:

	/**
	* Parameters relevant for reweight and variations
	*/
	//@{
	/**
	* The shower variations
	*/
	map<string,ShowerVariation> showerVariations_;

	/**
	* Command to add a shower variation
	*/
	string doAddVariation(string);

	/**
	* A reweighting factor applied by the showering
	*/
	double reweight_;

	/**
	* The shower variation weights
	*/
	map<string,double> currentWeights_;
	//@}
	};

	}

	#endif /* HERWIG_ShowerHandler_H */
	diff --git a/Utilities/Makefile.am b/Utilities/Makefile.am
	--- a/Utilities/Makefile.am
	+++ b/Utilities/Makefile.am
	@@ -1,46 +1,47 @@
	SUBDIRS = XML Statistics
	noinst_LTLIBRARIES = libHwUtils.la

	libHwUtils_la_SOURCES = \
	EnumParticles.h \
	Interpolator.tcc Interpolator.h \
	Kinematics.cc Kinematics.h \
	Progress.h Progress.cc \
	Maths.h Maths.cc \
	StandardSelectors.cc StandardSelectors.h\
	Histogram.cc Histogram.fh Histogram.h \
	GaussianIntegrator.cc GaussianIntegrator.h \
	GaussianIntegrator.tcc \
	Statistic.h HerwigStrategy.cc HerwigStrategy.h \
	GSLIntegrator.h GSLIntegrator.tcc \
	GSLBisection.h GSLBisection.tcc GSLHelper.h \
	+Reshuffler.h Reshuffler.cc \
	expm-1.h \
	HiggsLoopFunctions.h AlphaS.h

	nodist_libHwUtils_la_SOURCES = hgstamp.inc
	BUILT_SOURCES = hgstamp.inc
	CLEANFILES = hgstamp.inc

	HGVERSION := $(shell hg -R $(top_srcdir) parents --template '"Herwig {node\|short} ({branch})"' 2> /dev/null \|\| echo \"$(PACKAGE_STRING)\" \|\| true )

	.PHONY: update_hgstamp
	hgstamp.inc: update_hgstamp
	@[ -f $@ ] \|\| touch $@
	@echo '$(HGVERSION)' \| cmp -s $@ - \|\| echo '$(HGVERSION)' > $@

	libHwUtils_la_LIBADD = \
	XML/libHwXML.la \
	Statistics/libHwStatistics.la


	check_PROGRAMS = utilities_test
	utilities_test_SOURCES = \
	tests/utilitiesTestsMain.cc \
	tests/utilitiesTestsGlobalFixture.h \
	tests/utilitiesTestsKinematics.h \
	tests/utilitiesTestMaths.h \
	tests/utilitiesTestsStatistic.h
	utilities_test_LDADD = $(BOOST_UNIT_TEST_FRAMEWORK_LIBS) $(THEPEGLIB) -ldl libHwUtils.la
	utilities_test_LDFLAGS = $(AM_LDFLAGS) -export-dynamic $(BOOST_UNIT_TEST_FRAMEWORK_LDFLAGS) $(THEPEGLDFLAGS)
	utilities_test_CPPFLAGS = $(AM_CPPFLAGS) $(BOOST_CPPFLAGS)
	TESTS = utilities_test
	diff --git a/Utilities/Reshuffler.cc b/Utilities/Reshuffler.cc
	new file mode 100644
	--- /dev/null
	+++ b/Utilities/Reshuffler.cc
	@@ -0,0 +1,85 @@
	+// -- C++ --
	+//
	+// Reshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	+// Copyright (C) 2002-2019 The Herwig Collaboration
	+//
	+// Herwig is licenced under version 3 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 Reshuffler class.
	+//
	+
	+#include <config.h>
	+#include "Reshuffler.h"
	+#include "Herwig/Utilities/GSLBisection.h"
	+#include "Herwig/Shower/ShowerHandler.h"
	+
	+using namespace Herwig;
	+
	+void Reshuffler::reshuffle(const PVector& particles,
	+ const vector<Energy>& masses) const {
	+
	+ Energy zero (0.*GeV);
	+ Lorentz5Momentum Q (zero,zero,zero,zero);
	+
	+ for (PVector::const_iterator p = particles.begin();
	+ p != particles.end(); ++p) {
	+ Q += (**p).momentum();
	+ }
	+
	+ Boost beta = Q.findBoostToCM();
	+
	+ vector<Lorentz5Momentum> mbackup;
	+
	+ bool need_boost = (beta.mag2() > Constants::epsilon);
	+
	+ if (need_boost) {
	+
	+ for (PVector::const_iterator p = particles.begin();
	+ p != particles.end(); ++p) {
	+ Lorentz5Momentum mom = (**p).momentum();
	+ mbackup.push_back(mom);
	+ (**p).set5Momentum(mom.boost(beta));
	+ }
	+
	+ }
	+
	+ double xi;
	+
	+ ReshuffleEquation<PVector::const_iterator,vector<Energy>::const_iterator>
	+ solve (Q.m(),particles.begin(),particles.end(),
	+ masses.begin(),masses.end());
	+
	+ GSLBisection solver(1e-10,1e-8,10000);
	+
	+ try {
	+ xi = solver.value(solve,0.0,1.1);
	+ } catch (GSLBisection::GSLerror) {
	+ throw Exception("Failed to reshuffle.",Exception::eventerror);
	+ } catch (GSLBisection::IntervalError) {
	+ throw Exception("Failed to reshuffle.",Exception::eventerror);
	+ }
	+
	+ PVector::const_iterator p = particles.begin();
	+ vector<Energy>::const_iterator m = masses.begin();
	+ for (; p != particles.end(); ++p, ++m) {
	+
	+ Lorentz5Momentum rm = Lorentz5Momentum (xi(*p).momentum().x(),
	+ xi(*p).momentum().y(),
	+ xi(*p).momentum().z(),
	+ sqrt(sqr(*m) +
	+ xixi(sqr((p).momentum().t())-sqr((p).dataPtr()->mass()))));
	+
	+ rm.rescaleMass();
	+
	+ if (need_boost) {
	+ rm.boost(-beta);
	+ }
	+
	+ (**p).set5Momentum(rm);
	+
	+ }
	+
	+}
	diff --git a/Utilities/Reshuffler.h b/Utilities/Reshuffler.h
	new file mode 100644
	--- /dev/null
	+++ b/Utilities/Reshuffler.h
	@@ -0,0 +1,86 @@
	+// -- C++ --
	+//
	+// Reshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	+// Copyright (C) 2002-2019 The Herwig Collaboration
	+//
	+// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	+//
	+#ifndef HERWIG_Reshuffler_H
	+#define HERWIG_Reshuffler_H
	+//
	+// This is the declaration of the Reshuffler class.
	+//
	+
	+#include "ThePEG/Config/ThePEG.h"
	+
	+namespace Herwig {
	+
	+using namespace ThePEG;
	+
	+/**
	+ * \author Simon Platzer, Stephen Webster
	+ *
	+ * \brief The Reshuffler class implements reshuffling
	+ * of partons on their nominal mass shell to their constituent
	+ * mass shells.
	+ *
	+ */
	+class Reshuffler {
	+
	+protected:
	+
	+ /**
	+ * The function object defining the equation
	+ * to be solved.
	+ */
	+ template<class PIterator, class MIterator>
	+ struct ReshuffleEquation {
	+
	+ ReshuffleEquation (Energy q,
	+ PIterator pp_begin,
	+ PIterator pp_end,
	+ MIterator mm_begin,
	+ MIterator mm_end)
	+ : w(q), p_begin(pp_begin), p_end(pp_end),
	+ m_begin(mm_begin), m_end(mm_end) {}
	+
	+ typedef double ArgType;
	+ typedef double ValType;
	+
	+ static double aUnit() { return 1.; }
	+ static double vUnit() { return 1.; }
	+
	+ double operator() (double xi) const {
	+ double r = - w/GeV;
	+ PIterator p = p_begin;
	+ MIterator m = m_begin;
	+ for (; p != p_end; ++p, ++m) {
	+ r += sqrt(sqr(*m) +
	+ xixi(sqr((p).momentum().t())-sqr((p).dataPtr()->mass()))) / GeV;
	+ }
	+ return r;
	+ }
	+
	+ Energy w;
	+
	+ PIterator p_begin;
	+ PIterator p_end;
	+
	+ MIterator m_begin;
	+ MIterator m_end;
	+
	+ };
	+
	+ /**
	+ * Reshuffle to consitutent masses
	+ */
	+ void reshuffle(const PVector& particles,
	+ const vector<Energy>& masses) const;
	+
	+};
	+
	+}
	+
	+#endif /* HERWIG_Reshuffler_H */
	+

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 9:00 PM (1 d, 18 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806191
Default Alt Text: (239 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions