No OneTemporary
Actions

Size

68 KB

Subscribers

None

View Options

	diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc
	--- a/Hadronization/ClusterHadronizationHandler.cc
	+++ b/Hadronization/ClusterHadronizationHandler.cc
	@@ -1,303 +1,304 @@
	// -- C++ --
	//
	// ClusterHadronizationHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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","");

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

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

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


	void ClusterHadronizationHandler::persistentInput(PersistentIStream & is, int) {
	is >> _partonSplitter >> _clusterFinder >> _colourReconnector
	>> _clusterFissioner >> _lightClusterDecayer >> _clusterDecayer
	>> 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 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());
	// set the scale for coloured particles to just above the gluon mass squared
	// if less than this so they are classed as perturbative
	Energy2 Q02 = 1.01*sqr(getParticleData(ParticleID::g)->constituentMass());
	for(unsigned int ix=0;ix<currentlist.size();++ix) {
	if(currentlist[ix]->scale()<Q02) currentlist[ix]->scale(Q02);
	}

	// split the gluons
	_partonSplitter->split(currentlist);
	// form the clusters
	ClusterVector clusters =
	_clusterFinder->formClusters(currentlist);
	// 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(),
	mem_fun(&Particle::undecay));
	}
	}
	if (!lightOK) {
	- // currentHandler_ = 0;
	throw Exception("CluHad::handle(): tried LightClusterDecayer 10 times!",
	Exception::eventerror);
	}


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

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

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

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

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

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

	-void ClusterHadronizationHandler::_setChildren(ClusterVector clusters) const {
	-
	+void ClusterHadronizationHandler::_setChildren(const ClusterVector & clusters) const {
	// erase existing information about the partons' children
	tPVector partons;
	- for (ClusterVector::const_iterator cl = clusters.begin();
	- cl != clusters.end(); cl++) {
	- partons.push_back( (*cl)->colParticle() );
	- partons.push_back( (*cl)->antiColParticle() );
	+ for ( 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(), mem_fun(&Particle::undecay));

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

	diff --git a/Hadronization/ClusterHadronizationHandler.h b/Hadronization/ClusterHadronizationHandler.h
	--- a/Hadronization/ClusterHadronizationHandler.h
	+++ b/Hadronization/ClusterHadronizationHandler.h
	@@ -1,223 +1,209 @@
	// -- C++ --
	//
	// ClusterHadronizationHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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"

	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 {

	public:

	- /** @name Standard constructors and destructors. */
	- //@{
	- /**
	- * The default constructor.
	- */
	- ClusterHadronizationHandler()
	- : _minVirtuality2( 0.1GeV2 ), _maxDisplacement( 1.0e-10mm ),
	- _reduceToTwoComponents(true)
	- {}
	-
	- //@}
	-
	-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() {
	assert(currentHandler_);
	return currentHandler_;
	}

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

	/**
	* The minimum virtuality^2 of partons to use in calculating
	* distances.
	*/
	- Energy2 _minVirtuality2;
	+ 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;
	+ Length _maxDisplacement = 1.0e-10_mm;

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

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

	/**
	* Tag the constituents of the clusters as their parents
	*/
	- void _setChildren(ClusterVector clusters) const;
	+ void _setChildren(const ClusterVector & clusters) const;

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


	};


	}

	#endif /* HERWIG_ClusterHadronizationHandler_H */
	diff --git a/Hadronization/ColourReconnector.cc b/Hadronization/ColourReconnector.cc
	--- a/Hadronization/ColourReconnector.cc
	+++ b/Hadronization/ColourReconnector.cc
	@@ -1,488 +1,778 @@
	// -- C++ --
	//
	// ColourReconnector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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 ColourReconnector class.
	//

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

	+#include <ThePEG/Interface/Switch.h>
	+#include <ThePEG/Interface/Parameter.h>
	+
	+#include "Herwig/Utilities/Maths.h"

	using namespace Herwig;

	-typedef ClusterVector::iterator CluVecIt;
	+using CluVecIt = ColourReconnector::CluVecIt;
	+using Constants::pi;
	+using Constants::twopi;
	+

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

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

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

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

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

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


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


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


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

	const size_t nclusters = cv.size();

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

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

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

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

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

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

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

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

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

	// reduce temperature
	t *= _annealingFactor;
	}

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


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

	ClusterVector newcv = cv;

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

	// iterate over all clusters
	for (CluVecIt cit = newcv.begin(); cit != newcv.end(); cit++) {
	-
	// find the cluster which, if reconnected with *cit, would result in the
	// smallest sum of cluster masses
	// NB this method returns *cit if no reconnection partner can be found
	CluVecIt candidate = _findRecoPartner(cit, newcv);

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

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

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

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

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

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

	swap(cv,newcv);
	return;
	}

	+namespace {
	+ inline bool hasDiquark(CluVecIt cit) {
	+ for(int i = 0; i<(*cit)->numComponents(); i++) {
	+ if (DiquarkMatcher::Check(((cit)->particle(i)->dataPtr())))
	+ return true;
	+ }
	+ return false;
	+ }
	+}
	+
	+
	+// Implementation of the baryonic reconnection algorithm
	+void ColourReconnector::_doRecoBaryonic(ClusterVector & cv) const {
	+
	+ ClusterVector newcv = cv;
	+
	+ ClusterVector deleted; deleted.reserve(cv.size());
	+
	+ // try to avoid systematic errors by randomising the reconnection order
	+ long (*p_irnd)(long) = UseRandom::irnd;
	+ random_shuffle( newcv.begin(), newcv.end(), p_irnd );
	+
	+ // iterate over all clusters
	+ for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) {
	+ //avoid clusters already containing diuarks
	+ if (hasDiquark(cit)) continue;
	+
	+ //skip the cluster to be deleted later 3->2 cluster
	+ if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
	+ continue;
	+
	+ // Skip all found baryonic clusters, this biases the algorithm but implementing
	+ // something like re-reconnection is ongoing work
	+ if ((*cit)->numComponents()==3) continue;
	+
	+ // Find a candidate suitable for reconnection
	+ CluVecIt baryonic1, baryonic2;
	+ bool isBaryonicCandidate = false;
	+ CluVecIt candidate = _findPartnerBaryonic(cit, newcv,
	+ isBaryonicCandidate,
	+ deleted,
	+ baryonic1, baryonic2);
	+
	+ // skip this cluster if no possible reconnection partner can be found
	+ if ( !isBaryonicCandidate && candidate==cit )
	+ continue;
	+
	+ if ( isBaryonicCandidate
	+ && UseRandom::rnd() < _precoBaryonic ) {
	+ deleted.push_back(*baryonic2);
	+
	+ // Function that does the reconnection from 3 -> 2 clusters
	+ ClusterPtr b1, b2;
	+ _makeBaryonicClusters(cit,baryonic1,*baryonic2, b1, b2);
	+
	+ *cit = b1;
	+ *baryonic1 = b2;
	+
	+ // Baryonic2 is easily skipped in the next loop
	+ }
	+
	+ // Normal 2->2 Colour reconection
	+ if ( !isBaryonicCandidate
	+ && UseRandom::rnd() < _preco ) {
	+ auto reconnected = _reconnect(cit, candidate);
	+ *cit = reconnected.first;
	+ *candidate = reconnected.second;
	+ }
	+ }
	+
	+ // create a new vector of clusters except for the ones which are "deleted" during
	+ // baryonic reconnection
	+ ClusterVector clustervector;
	+ for ( const auto & cluster : newcv )
	+ if ( find(deleted.begin(),
	+ deleted.end(), cluster) == deleted.end() )
	+ clustervector.push_back(cluster);
	+
	+ swap(cv,clustervector);
	+}
	+
	+
	+
	+namespace {
	+
	+double calculateRapidityRF(const Lorentz5Momentum & q1,
	+ const Lorentz5Momentum & p2) {
	+ //calculate rapidity wrt the direction of q1
	+ //angle between the particles in the RF of cluster of q1
	+
	+ // calculate the z component of p2 w.r.t the direction of q1
	+ const Energy pz = p2.vect() * q1.vect().unit();
	+ if ( pz == ZERO ) return 0.;
	+
	+ // Transverse momentum of p2 w.r.t the direction of q1
	+ const Energy pt = sqrt(p2.vect().mag2() - sqr(pz));
	+
	+ // Transverse mass pf p2 w.r.t to the direction of q1
	+ const Energy mtrans = sqrt(p2.mass()p2.mass() + (ptpt));
	+
	+ // Correct formula
	+ const double y2 = log((p2.t() + abs(pz))/mtrans);
	+
	+ return ( pz < ZERO ) ? -y2 : y2;
	+}
	+
	+}
	+
	+
	+CluVecIt ColourReconnector::_findPartnerBaryonic(
	+ CluVecIt cl, ClusterVector & cv,
	+ bool & baryonicCand,
	+ const ClusterVector& deleted,
	+ CluVecIt &baryonic1,
	+ CluVecIt &baryonic2 ) const {
	+
	+ using Constants::pi;
	+ using Constants::twopi;
	+
	+ // Returns a candidate for possible reconnection
	+ CluVecIt candidate = cl;
	+
	+ bool bcand = false;
	+
	+ double maxrap = 0.0;
	+ double minrap = 0.0;
	+ double maxrapNormal = 0.0;
	+ double minrapNormal = 0.0;
	+ double maxsumnormal = 0.0;
	+
	+ double maxsum = 0.0;
	+ double secondsum = 0.0;
	+
	+ for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
	+ //avoid looping over clusters containing diquarks
	+ if ( hasDiquark(cit) ) continue;
	+ if ( (*cit)->numComponents()==3 ) continue;
	+ if ( cit==cl ) continue;
	+
	+ //skip the cluster to be deleted later 3->2 cluster
	+ if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() )
	+ continue;
	+
	+ if ( (cl)->isBeamCluster() && (cit)->isBeamCluster() )
	+ continue;
	+
	+ // stop it putting far apart clusters together
	+ if ( ( (cl).vertex()-(cit).vertex() ).m() >_maxDistance )
	+ continue;
	+
	+ const bool Colour8 =
	+ _isColour8( (cl)->colParticle(), (cit)->antiColParticle() )
	+ \|\|
	+ _isColour8( (cit)->colParticle(), (cl)->antiColParticle() ) ;
	+ if ( Colour8 ) continue;
	+
	+ // boost into restframe of cluster 1
	+ Lorentz5Momentum cl1 = (*cl)->momentum();
	+ Lorentz5Momentum cl2 = (*cit)->momentum();
	+ // boost cl1 constituents int restframe of cl1:
	+ Boost boostv(-cl1.boostVector());
	+ cl1.boost(boostv);
	+ cl2.boost(boostv);
	+
	+ Lorentz5Momentum p1col = (*cl)->colParticle()->momentum();
	+ Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum();
	+ Lorentz5Momentum p2col = (*cit)->colParticle()->momentum();
	+ Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum();
	+ p1col.boost(boostv);
	+ p1anticol.boost(boostv);
	+ // boost cl2 constituents int rf of cl1:
	+ p2col.boost(boostv);
	+ p2anticol.boost(boostv);
	+
	+ // calculate the rapidity of the other constituents of the clusters
	+ // w.r.t axis of p1anticol.vect.unit
	+ const double rapq = calculateRapidityRF(p1anticol,p2col);
	+ const double rapqbar = calculateRapidityRF(p1anticol,p2anticol);
	+
	+ // configuration for normal CR
	+ if ( rapq > 0.0 && rapqbar < 0.0
	+ && rapq > maxrap
	+ && rapqbar < minrap ) {
	+ maxrap = rapq;
	+ minrap = rapqbar;
	+ //sum of rapidities of quarks
	+ const double normalsum = abs(rapq) + abs(rapqbar);
	+ if ( normalsum > maxsumnormal ) {
	+ maxsumnormal = normalsum;
	+ maxrapNormal = rapq;
	+ minrapNormal = rapqbar;
	+ bcand = false;
	+ candidate = cit;
	+ }
	+ }
	+
	+ if ( rapq < 0.0 && rapqbar >0.0
	+ && rapqbar > maxrapNormal
	+ && rapq < minrapNormal ) {
	+ maxrap = rapqbar;
	+ minrap = rapq;
	+ const double sumrap = abs(rapqbar) + abs(rapq);
	+ // first candidate gets here. If second baryonic candidate has higher Ysum than the first
	+ // one, the second candidate becomes the first one and the first the second.
	+ if (sumrap > maxsum) {
	+ if(maxsum != 0){
	+ baryonic2 = baryonic1;
	+ baryonic1 = cit;
	+ bcand = true;
	+ } else {
	+ baryonic1 = cit;
	+ }
	+ maxsum = sumrap;
	+ } else {
	+ if (sumrap > secondsum && sumrap != maxsum) {
	+ secondsum = sumrap;
	+ bcand = true;
	+ baryonic2 = cit;
	+ }
	+ }
	+ }
	+
	+ }
	+
	+ if(bcand == true){
	+ baryonicCand = true;
	+ }
	+
	+ return candidate;
	+}
	+

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

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

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

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

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

	// stop it putting far apart clusters together
	if(((cl).vertex()-(cit).vertex()).m()>_maxDistance) continue;

	// momenta of the old clusters
	Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() +
	(*cl)->antiColParticle()->momentum();
	Lorentz5Momentum p2 = (*cit)->colParticle()->momentum() +
	(*cit)->antiColParticle()->momentum();

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

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

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

	+// forms two baryonic clusters from three clusters
	+void ColourReconnector::_makeBaryonicClusters(
	+ ClusterPtr &c1, ClusterPtr &c2,
	+ ClusterPtr &c3,
	+ ClusterPtr &newcluster1,
	+ ClusterPtr &newcluster2) const{
	+
	+ //make sure they all have 2 components
	+ assert(c1->numComponents()==2);
	+ assert(c2->numComponents()==2);
	+ assert(c3->numComponents()==2);
	+ //abandon childs
	+ c1->colParticle()->abandonChild(c1);
	+ c1->antiColParticle()->abandonChild(c1);
	+ c2->colParticle()->abandonChild(c2);
	+ c2->antiColParticle()->abandonChild(c2);
	+ c3->colParticle()->abandonChild(c3);
	+ c3->antiColParticle()->abandonChild(c3);
	+
	+ newcluster1 = new_ptr(Cluster(c1->colParticle(),c2->colParticle(), c3->colParticle()));
	+ c1->colParticle()->addChild(newcluster1);
	+ c2->colParticle()->addChild(newcluster1);
	+ c3->colParticle()->addChild(newcluster1);
	+ newcluster1->setVertex(LorentzPoint());
	+ newcluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->antiColParticle(),
	+ c3->antiColParticle()));
	+ c1->antiColParticle()->addChild(newcluster2);
	+ c2->antiColParticle()->addChild(newcluster2);
	+ c3->antiColParticle()->addChild(newcluster2);
	+ newcluster2->setVertex(LorentzPoint());
	+}
	+

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

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

	+c1->colParticle()->abandonChild(c1);
	+c2->antiColParticle()->abandonChild(c2);
	+
	ClusterPtr newCluster1
	= new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) );

	+ c1->colParticle()->addChild(newCluster1);
	+ c2->antiColParticle()->addChild(newCluster1);
	+
	newCluster1->setVertex(0.5*( c1->colParticle()->vertex() +
	- c2->antiColParticle()->vertex() ));
	+ c2->antiColParticle()->vertex() ));

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

	+ c1->antiColParticle()->addChild(newCluster2);
	+ c2->colParticle()->addChild(newCluster2);
	+
	newCluster2->setVertex(0.5*( c2->colParticle()->vertex() +
	- c1->antiColParticle()->vertex() ));
	+ c1->antiColParticle()->vertex() ));

	if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true);
	if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true);

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


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

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

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

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

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

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


	-bool ColourReconnector::isColour8(cPPtr p, cPPtr q) const {
	+
	+bool ColourReconnector::_isColour8(tcPPtr p, tcPPtr q) const {
	bool octet = false;

	// make sure we have a triplet and an anti-triplet
	if ( ( p->hasColour() && q->hasAntiColour() ) \|\|
	( p->hasAntiColour() && q->hasColour() ) ) {

	// true if p and q are originated from a colour octet
	if ( !p->parents().empty() && !q->parents().empty() ) {
	octet = ( p->parents()[0] == q->parents()[0] ) &&
	( p->parents()[0]->data().iColour() == PDT::Colour8 );
	}

	// (Final) option: check if same colour8 parent
	// or already found an octet.
	if(_octetOption==0\|\|octet) return octet;

	// (All) option handling more octets
	// by browsing particle history/colour lines.
	tColinePtr cline,aline;

	// Get colourlines form final states.
	if(p->hasColour() && q->hasAntiColour()) {
	cline = p-> colourLine();
	aline = q->antiColourLine();
	}
	else {
	cline = q-> colourLine();
	aline = p->antiColourLine();
	}

	// Follow the colourline of p.
	if ( !p->parents().empty() ) {
	tPPtr parent = p->parents()[0];
	while (parent) {
	if(parent->data().iColour() == PDT::Colour8) {
	// Coulour8 particles should have a colour
	// and an anticolour line. Currently the
	// remnant has none of those. Since the children
	// of the remnant are not allowed to emit currently,
	// the colour octet remnant is handled by the return
	// statement above. The assert also catches other
	// colour octets without clines. If the children of
	// a remnant should be allowed to emit, the remnant
	// should get appropriate colour lines and
	// colour states.
	// See Ticket: #407
	- // assert(parent->colourLine()&&parent->antiColourLine());
	+ // assert(parent->colourLine()&&parent->antiColourLine());
	octet = (parent-> colourLine()==cline &&
	parent->antiColourLine()==aline);
	}
	if(octet\|\|parent->parents().empty()) break;
	parent = parent->parents()[0];
	}
	}
	}

	return octet;
	}


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

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


	void ColourReconnector::Init() {

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


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


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

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

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


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


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

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

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

	static Parameter<ColourReconnector,Length> interfaceMaxDistance
	("MaxDistance",
	"Maximum distance between the clusters at which to consider rearrangement"
	" to avoid colour reconneections of displaced vertices",
	&ColourReconnector::_maxDistance, femtometer, 1000.femtometer, 0.0femtometer, 1e100*femtometer,
	false, false, Interface::limited);

	+
	static Switch<ColourReconnector,unsigned int> interfaceOctetTreatment
	("OctetTreatment",
	"Which octets are not allowed to be reconnected",
	&ColourReconnector::_octetOption, 0, false, false);
	static SwitchOption interfaceOctetTreatmentFinal
	(interfaceOctetTreatment,
	"Final",
	"Only prevent for the final (usuaslly non-perturbative) g -> q qbar splitting",
	0);
	static SwitchOption interfaceOctetTreatmentAll
	(interfaceOctetTreatment,
	"All",
	"Prevent for all octets",
	1);

	+}

	-}
	diff --git a/Hadronization/ColourReconnector.h b/Hadronization/ColourReconnector.h
	--- a/Hadronization/ColourReconnector.h
	+++ b/Hadronization/ColourReconnector.h
	@@ -1,247 +1,259 @@
	// -- C++ --
	//
	// ColourReconnector.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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_ColourReconnector_H
	#define HERWIG_ColourReconnector_H

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

	namespace Herwig {


	using namespace ThePEG;

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

	public:

	- /** @name Standard constructors and destructors. */
	- //@{
	- /**
	- * Default constructor.
	- */
	- ColourReconnector() :
	- _algorithm(0),
	- _annealingFactor(0.9),
	- _annealingSteps(50),
	- _clreco(0),
	- _initTemp(0.1),
	- _preco(0.5),
	- _triesPerStepFactor(5.0),
	- _maxDistance(1000.*femtometer),
	- _octetOption(0)
	- {}
	- //@}
	-
	/**
	* Does the colour rearrangement, starting out from the list of particles in
	* the event record and the collection of "usual" clusters passed as
	* arguments. If the actual rearrangement is accepted, the initial collection of
	* clusters is overridden by the old ones.
	*/
	void rearrange(ClusterVector & clusters);

	+ using CluVecIt = ClusterVector::iterator;

	private:

	/** PRIVATE MEMBER FUNCTIONS */

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


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

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

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

	+
	+ /**
	+ * Baryonic Colour Reconnection model
	+ */
	+ void _doRecoBaryonic(ClusterVector & cv) const;
	+
	+
	+ void _makeBaryonicClusters(ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3,
	+ ClusterPtr &newcluster1, ClusterPtr &newcluster2) const;
	+
	+
	/**
	* @brief Finds the cluster in cv which, if reconnected with the given
	* cluster cl, would result in the smallest sum of cluster masses.
	* If no reconnection partner can be found, a pointer to the
	* original Cluster cl is returned.
	* @arguments cv cluster vector
	* cl cluster iterator (must be from cv) which wants to have a reconnection partner
	* @return iterator to the found cluster, or the original cluster pointer if
	* no mass-reducing combination can be found
	*/
	- ClusterVector::iterator _findRecoPartner(ClusterVector::iterator cl,
	- ClusterVector & cv) const;
	+
	+
	+ CluVecIt _findRecoPartner(CluVecIt cl, ClusterVector & cv) const;
	+
	+ CluVecIt _findPartnerRapidity(CluVecIt cl, ClusterVector & cv) const;
	+
	+ CluVecIt _findPartnerBaryonic(CluVecIt cl, ClusterVector & cv,
	+ bool & tetraCand,
	+ const ClusterVector& a,
	+ CluVecIt &baryonic1,
	+ CluVecIt &baryonic2 ) const;
	+
	+
	+

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

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

	+
	+ /** DATA MEMBERS */
	+
	+ /**
	+ * Specifies the colour reconnection algorithm to be used.
	+ */
	+ int _algorithm = 0;
	+
	+ /**
	+ * The annealing factor is the ratio of two successive temperature steps:
	+ * T_n = _annealingFactor * T_(n-1)
	+ */
	+ double _annealingFactor = 0.9;
	+
	+ /**
	+ * Number of temperature steps in the statistical annealing algorithm
	+ */
	+ unsigned int _annealingSteps = 50;
	+
	+ /**
	+ * Do we do colour reconnections?
	+ */
	+ int _clreco = 0;
	+
	+ /**
	+ * Factor used to determine the initial temperature according to
	+ * InitialTemperature = _initTemp * median {energy changes in a few random
	+ * rearrangements}
	+ */
	+ double _initTemp = 0.1;
	+
	+ /**
	+ * Probability that a found reconnection possibility is actually accepted.
	+ */
	+ double _preco = 0.5;
	+
	+
	+ double _precoBaryonic = 0.5;
	+
	+ /**
	+ * The number of tries per temperature steps is the number of clusters times
	+ * this factor.
	+ */
	+ double _triesPerStepFactor = 5.0;
	+ /**
	+ * maximum allowed distance in the eta phi space for reconnection to occur
	+ */
	+
	+ /**
	+ * Maximium distance for reconnections
	+ */
	+ Length _maxDistance = picometer;
	+
	/**
	* @return true, if the two partons are splitting products of the same
	* gluon
	*/
	- bool isColour8(cPPtr p, cPPtr q) const;
	-
	+ bool _isColour8(tcPPtr p, tcPPtr q) const;
	+
	+ /**
	+ * Option for handling octets
	+ */
	+ unsigned int _octetOption = 0;

	public:

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

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

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

	protected:

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

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


	private:

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

	-private:
	-
	- /** DATA MEMBERS */
	-
	- /**
	- * Specifies the colour reconnection algorithm to be used.
	- */
	- int _algorithm;
	-
	- /**
	- * The annealing factor is the ratio of two successive temperature steps:
	- * T_n = _annealingFactor * T_(n-1)
	- */
	- double _annealingFactor;
	-
	- /**
	- * Number of temperature steps in the statistical annealing algorithm
	- */
	- unsigned _annealingSteps;
	-
	- /**
	- * Do we do colour reconnections?
	- */
	- int _clreco;
	-
	- /**
	- * Factor used to determine the initial temperature according to
	- * InitialTemperature = _initTemp * median {energy changes in a few random
	- * rearrangements}
	- */
	- double _initTemp;
	-
	- /**
	- * Probability that a found reconnection possibility is actually accepted.
	- */
	- double _preco;
	-
	- /**
	- * The number of tries per temperature steps is the number of clusters times
	- * this factor.
	- */
	- double _triesPerStepFactor;
	-
	- /**
	- * Maximium distance for reconnections
	- */
	- Length _maxDistance;
	-
	- /**
	- * Option for handling octets
	- */
	- unsigned int _octetOption;

	};


	}

	#endif /* HERWIG_ColourReconnector_H */
	+
	diff --git a/Hadronization/PartonSplitter.cc b/Hadronization/PartonSplitter.cc
	--- a/Hadronization/PartonSplitter.cc
	+++ b/Hadronization/PartonSplitter.cc
	@@ -1,149 +1,218 @@
	// -- C++ --
	//
	// PartonSplitter.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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 PartonSplitter class.
	//

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

	using namespace Herwig;

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

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

	void PartonSplitter::persistentOutput(PersistentOStream & os) const {
	- os << _quarkSelector << ounit(_gluonDistance,femtometer);
	+ os << _quarkSelector << ounit(_gluonDistance,femtometer)
	+ << _splitGluon << _splitPwtUquark << _splitPwtDquark << _splitPwtSquark;
	}

	void PartonSplitter::persistentInput(PersistentIStream & is, int) {
	- is >> _quarkSelector >> iunit(_gluonDistance,femtometer);
	+ is >> _quarkSelector >> iunit(_gluonDistance,femtometer)
	+ >> _splitGluon >> _splitPwtUquark >> _splitPwtDquark >> _splitPwtSquark;
	}

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

	+
	void PartonSplitter::Init() {

	static ClassDocumentation<PartonSplitter> documentation
	("This class is reponsible of the nonperturbative splitting of partons");
	-
	+
	+ static Switch<PartonSplitter,int> interfaceSplit
	+ ("Split",
	+ "Option for different splitting options",
	+ &PartonSplitter::_splitGluon, 1, false, false);
	+ static SwitchOption interfaceSplitDefault
	+ (interfaceSplit,
	+ "ud",
	+ "Normal cluster splitting where only u and d quarks are drawn is used.",
	+ 0);
	+ static SwitchOption interfaceSplitAll
	+ (interfaceSplit,
	+ "uds",
	+ "Alternative cluster splitting where all light quark pairs (u, d, s) can be drawn.",
	+ 1);
	+
	+ static Parameter<PartonSplitter,double> interfaceSplitPwtUquark
	+ ("SplitPwtUquark",
	+ "Weight for splitting in U quarks",
	+ &PartonSplitter::_splitPwtUquark, 1, 0.0, 1.0,
	+ false, false, Interface::limited);
	+ static Parameter<PartonSplitter,double> interfaceSplitPwtDquark
	+ ("SplitPwtDquark",
	+ "Weight for splitting in D quarks",
	+ &PartonSplitter::_splitPwtDquark, 1, 0.0, 1.0,
	+ false, false, Interface::limited);
	+ static Parameter<PartonSplitter,double> interfaceSplitPwtSquark
	+ ("SplitPwtSquark",
	+ "Weight for splitting in S quarks",
	+ &PartonSplitter::_splitPwtSquark, 0.5, 0.0, 1.0,
	+ false, false, Interface::limited);
	+
	}

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

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

	(*pit)->antiColourLine()->addAntiColoured(ptrQbar);
	(*pit)->addChild(ptrQbar);
	newtag.push_back(ptrQbar);

	// set the life length of gluon
	Length distance = UseRandom::rndExp(_gluonDistance);
	(pit).setLifeLength((distance/(pit).mass())(*pit).momentum());
	// assume quarks same position as gluon
	ptrQ ->setVertex((**pit).decayVertex());
	ptrQ ->setLifeLength(Lorentz5Distance());
	ptrQbar->setVertex((**pit).decayVertex());
	ptrQbar->setLifeLength(Lorentz5Distance());
	}
	swap(tagged,newtag);
	}

	void PartonSplitter::splitTimeLikeGluon(tcPPtr ptrGluon,
	PPtr & ptrQ,
	PPtr & ptrQbar){
	// select the quark flavour
	- tPDPtr quark = _quarkSelector.select(UseRandom::rnd());
	+ tPDPtr quark;
	+ long idNew=0;
	+
	+ switch(_splitGluon){
	+ case 0:
	+ quark = _quarkSelector.select(UseRandom::rnd());
	+ break;
	+ case 1:
	+ // Only allow light quarks u,d,s with the probabilities
	+ double prob_d = _splitPwtDquark;
	+ double prob_u = _splitPwtUquark;
	+ double prob_s = _splitPwtSquark;
	+
	+ int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
	+ switch(choice) {
	+ case 0: idNew = ThePEG::ParticleID::u; break;
	+ case 1: idNew = ThePEG::ParticleID::d; break;
	+ case 2: idNew = ThePEG::ParticleID::s; break;
	+ }
	+ ptrQ = getParticle(idNew);
	+ ptrQbar = getParticle(-idNew);
	+ break;
	+ }
	// Solve the kinematics of the two body decay G --> Q + Qbar
	Lorentz5Momentum momentumQ;
	Lorentz5Momentum momentumQbar;
	double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 );
	using Constants::pi;
	double phiStar = UseRandom::rnd( -pi , pi );
	- Energy constituentQmass = quark->constituentMass();

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

	void PartonSplitter::doinit() {
	Interfaced::doinit();
	// calculate the probabilties for the gluon to branch into each quark type
	// based on the available phase-space, as in fortran.
	Energy mg=getParticleData(ParticleID::g)->constituentMass();
	for( int ix=1; ix<6; ++ix ) {
	PDPtr quark = getParticleData(ix);
	Energy pcm = Kinematics::pstarTwoBodyDecay(mg,quark->constituentMass(),
	quark->constituentMass());
	if(pcm>ZERO) _quarkSelector.insert(pcm/GeV,quark);
	}
	if(_quarkSelector.empty())
	throw InitException() << "At least one quark must have constituent mass less "
	<< "then the constituent mass of the gluon in "
	<< "PartonSplitter::doinit()" << Exception::runerror;
	}
	diff --git a/Hadronization/PartonSplitter.h b/Hadronization/PartonSplitter.h
	--- a/Hadronization/PartonSplitter.h
	+++ b/Hadronization/PartonSplitter.h
	@@ -1,144 +1,172 @@
	// -- C++ --
	//
	// PartonSplitter.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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_PartonSplitter_H
	#define HERWIG_PartonSplitter_H

	#include "CluHadConfig.h"
	#include <ThePEG/Interface/Interfaced.h>
	#include <ThePEG/Utilities/Selector.h>
	#include "PartonSplitter.fh"

	namespace Herwig {


	using namespace ThePEG;


	/** \ingroup Hadronization
	* \class PartonSplitter
	* \brief This class splits the gluons from the end of the shower.
	* \author Philip Stephens
	* \author Alberto Ribon
	*
	* This class does all of the nonperturbative parton splittings needed
	* immediately after the end of the showering (both initial and final),
	* as very first step of the cluster hadronization.
	*
	- * \todo change so quark weights can be varied and quarks other
	- * than u and d can be produced
	+ * the quarks are attributed with different weights for the splitting
	+ * by default only the splitting in u and d quarks is allowed
	+ * the option "set /Herwig/Hadronization/PartonSplitter:Split 1"
	+ * allows for additional splitting into s quarks based on some weight
	+ * in order for that to work the mass of the strange quark has to be changed
	+ * from the default value s.t. m_g > 2m_s
	+ *
	*
	* * @see \ref PartonSplitterInterfaces "The interfaces"
	* defined for PartonSplitter.
	*/
	class PartonSplitter: public Interfaced {

	public:

	/**
	* Default constructor
	*/
	- PartonSplitter() : _gluonDistance(ZERO)
	+ PartonSplitter() :
	+ _splitPwtUquark(1),
	+ _splitPwtDquark(1),
	+ _splitPwtSquark(0.5),
	+ _gluonDistance(ZERO),
	+ _splitGluon(0)
	{}

	/**
	* This method does the nonperturbative splitting of:
	* time-like gluons. At the end of the shower the gluons should be
	* on a "physical" mass shell and should therefore be time-like.
	* @param tagged The tagged particles to be split
	* @return The particles which were not split and the products of splitting.
	*/
	void split(PVector & tagged);

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

	protected:

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

	private:

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

	/**
	* Non-perturbatively split a time-like gluon,
	* if something goes wrong null pointers are returned.
	* @param gluon The gluon to be split
	* @param quark The quark produced in the splitting
	* @param anti The antiquark produced in the splitting
	*/
	void splitTimeLikeGluon(tcPPtr gluon, PPtr & quark, PPtr & anti);

	+ // probabilities for the different quark types
	+ double _splitPwtUquark;
	+ double _splitPwtDquark;
	+ double _splitPwtSquark;
	+
	+
	private:

	/**
	* The selector to pick the type of quark
	*/
	Selector<PDPtr,double> _quarkSelector;

	/**
	+ * A pointer to a Herwig::HadronSelector object for generating hadrons.
	+ */
	+
	+ /**
	* c tau for gluon decays
	*/
	Length _gluonDistance;
	+
	+ /**
	+ * Flag used to determine between normal gluon splitting and alternative gluon splitting
	+ */
	+ int _splitGluon;
	+
	+
	};

	}

	#endif /* HERWIG_PartonSplitter_H */
	diff --git a/src/snippets/BaryonicReconnection.in b/src/snippets/BaryonicReconnection.in
	new file mode 100644
	--- /dev/null
	+++ b/src/snippets/BaryonicReconnection.in
	@@ -0,0 +1,17 @@
	+# Set strange quark mass to 0.45 in order to allow alternative gluon splitting
	+set /Herwig/Particles/s:ConstituentMass 0.45*GeV
	+set /Herwig/Particles/sbar:ConstituentMass 0.45*GeV
	+
	+# Use Baryonic Colour Reconnection Model
	+set /Herwig/Hadronization/ColourReconnector:Algorithm BaryonicReco
	+
	+# Allow alternative gluon splitting
	+set /Herwig/Hadronization/PartonSplitter:Split uds
	+
	+# Parameters for the Baryonic Reconnection Model
	+set /Herwig/Hadronization/ColourReconnector:ReconnectionProbability 0.772606
	+set /Herwig/Hadronization/ColourReconnector:ReconnectionProbabilityBaryonic 0.477612
	+set /Herwig/UnderlyingEvent/MPIHandler:pTmin0 3.053252
	+set /Herwig/UnderlyingEvent/MPIHandler:InvRadius 1.282032
	+set /Herwig/Hadronization/HadronSelector:PwtSquark 0.291717
	+set /Herwig/Hadronization/PartonSplitter:SplitPwtSquark 0.824135

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions