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.1*GeV2, ZERO, 10.0*GeV2,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.0*mm, 1.0e-9*mm,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.0*m0*n ){
+ 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.0*m0*(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