No OneTemporary
Actions

Size

58 KB

Subscribers

None

View Options

	diff --git a/Hadronization/DynamicClusterFissioner.cc b/Hadronization/DynamicClusterFissioner.cc
	new file mode 100644
	--- /dev/null
	+++ b/Hadronization/DynamicClusterFissioner.cc
	@@ -0,0 +1,404 @@
	+// -- C++ --
	+//
	+// ClusterFissioner.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	+// Copyright (C) 2002-2019 The Herwig Collaboration
	+//
	+// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	+//
	+//
	+// Thisk is the implementation of the non-inlined, non-templated member
	+// functions of the ClusterFissioner class.
	+//
	+
	+#include "DynamicClusterFissioner.h"
	+#include <ThePEG/Interface/ClassDocumentation.h>
	+#include <ThePEG/Interface/Reference.h>
	+#include <ThePEG/Interface/Parameter.h>
	+#include <ThePEG/Interface/Switch.h>
	+#include <ThePEG/Persistency/PersistentOStream.h>
	+#include <ThePEG/Persistency/PersistentIStream.h>
	+#include <ThePEG/PDT/EnumParticles.h>
	+#include "Herwig/Utilities/Kinematics.h"
	+#include "ThePEG/Repository/UseRandom.h"
	+#include "ThePEG/Repository/EventGenerator.h"
	+#include <ThePEG/Utilities/DescribeClass.h>
	+
	+using namespace Herwig;
	+
	+DescribeClass<DynamicClusterFissioner,ClusterFissioner>
	+describeDynamicClusterFissioner("Herwig::DynamicClusterFissioner","Herwig.so");
	+
	+DynamicClusterFissioner::DynamicClusterFissioner() :
	+ _Qqtilde(4*GeV),
	+ _Qg2tilde(4*GeV),
	+ _AngOrdFission(1),
	+ _restrictGluon(1)
	+{}
	+
	+IBPtr DynamicClusterFissioner::clone() const {
	+ return new_ptr(*this);
	+}
	+
	+IBPtr DynamicClusterFissioner::fullclone() const {
	+ return new_ptr(*this);
	+}
	+
	+void DynamicClusterFissioner::persistentOutput(PersistentOStream & os) const {
	+ os << _AngOrdFission << _restrictGluon << ounit(_Qqtilde,GeV) << ounit(_Qg2tilde,GeV) << _dynamicGluonMassGenerator << _dynamicPartonSplitter;
	+}
	+
	+void DynamicClusterFissioner::persistentInput(PersistentIStream & is, int) {
	+ is >> _AngOrdFission >> _restrictGluon >> iunit(_Qqtilde,GeV) >> iunit(_Qg2tilde,GeV) >> _dynamicGluonMassGenerator >> _dynamicPartonSplitter;
	+}
	+
	+void DynamicClusterFissioner::Init() {
	+
	+ static Parameter<DynamicClusterFissioner,Energy> interfaceQqtilde
	+ ("Qqtilde",
	+ "Scale of the non-pert. quark splitting in the cluster fission",
	+ &DynamicClusterFissioner::_Qqtilde, GeV, 4.0GeV, 1.5GeV, 10000.0*GeV,false,false,false); //Set upper/lower bound and default later
	+
	+ static Parameter<DynamicClusterFissioner,Energy> interfaceQg2tilde
	+ ("Qg2tilde",
	+ "Upper scale for the Sudakov for the non-pert. gluon splitting in the cluster fission",
	+ &DynamicClusterFissioner::_Qg2tilde, GeV, 4.0GeV, 0.5GeV, 10000.0*GeV,false,false,false); //Set upper/lower bound and default later
	+
	+ static Switch<DynamicClusterFissioner,int> interfaceAngOrdFission
	+ ("AngularOrderedFission",
	+ "Option to switch on angular ordering in the dynamic fission model",
	+ &DynamicClusterFissioner::_AngOrdFission,0,false,false);
	+ static SwitchOption interfaceAngOrdFissionNo
	+ (interfaceAngOrdFission,
	+ "No",
	+ "No (Qqtilde and Qg2tilde independent)",
	+ 0);
	+ static SwitchOption interfaceAngOrdFissionYes
	+ (interfaceAngOrdFission,
	+ "Yes",
	+ "Qg2tilde=(1-z)*Qqtilde",
	+ 1);
	+
	+
	+ static Switch<DynamicClusterFissioner,int> interfacerestrictGluon
	+ ("restrictGluon",
	+ "Option to restrict the gluon virtuality in the dynamic fission in order to make sure that the splitting works out kinematically",
	+ &DynamicClusterFissioner::_restrictGluon,0,false,false);
	+ static SwitchOption interfacerestrictGluonNo
	+ (interfacerestrictGluon,
	+ "No",
	+ "Gluon can have virtuality up to Qg2tilde. Might be necessary to start fission again with new z and qtilde for the quark",
	+ 0);
	+ static SwitchOption interfacerestrictGluonYes
	+ (interfacerestrictGluon,
	+ "Yes",
	+ "Gluon virtuality is restricted to the range that the kinematics allow once that z and qitlde are fixed for the quark",
	+ 1);
	+
	+ static Reference<DynamicClusterFissioner,DynamicGluonMassGenerator> interfaceDynamicGluonMassGenerator
	+ ("GluonMassGenerator",
	+ "Set a reference to a gluon mass generator.",
	+ &DynamicClusterFissioner::_dynamicGluonMassGenerator, false, false, true, true, false);
	+
	+ static Reference<DynamicClusterFissioner,DynamicPartonSplitter> interfaceDynamicPartonSplitter
	+ ("PartonSplitter",
	+ "Set a reference to a parton splitter",
	+ &DynamicClusterFissioner::_dynamicPartonSplitter, false, false, true, true, false);
	+
	+}
	+
	+
	+
	+
	+InvEnergy DynamicClusterFissioner::Pqzproposal(double z,Energy qtilde, Energy Qqtilde, Energy mq) const{
	+ return _dynamicGluonMassGenerator->ClusterAlphaS(0.GeV2)2.(1.-sqr(mq/Qqtilde))/(qtilde(1.-z));
	+}
	+
	+InvEnergy DynamicClusterFissioner::Pqz(double z, Energy qtilde, Energy mq) const {
	+ if ( z*qtilde > mq){
	+ return _dynamicGluonMassGenerator->ClusterAlphaS(sqr(z(1.0-z)qtilde))( (1.+zz)/(1.-z) - ((2.sqr(mq))/( z(1.-z)*sqr(qtilde) )))/qtilde;}
	+ else {
	+ return 0./(1.*GeV);
	+ }
	+}
	+
	+
	+void DynamicClusterFissioner::dynamicFission(tPPtr ptrP, tPPtr ptrPbar, PPtr& ptrQ, PPtr& ptrQbar) const {
	+
	+ LorentzMomentum P1, P2;
	+ bool fromcolor;
	+
	+ Energy mq2, mg, MP1, Qqtilde, Qgtilde, qtilde, mq;
	+
	+ double z;
	+
	+ //minimal quark mass
	+ Energy m0, mu, ms, md;
	+ 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;}
	+
	+
	+ bool repeat2=true;
	+
	+
	+ LorentzMomentum Pcl(ptrP->momentum()+ptrPbar->momentum());
	+ Energy Mcl = Pcl.m();
	+ //boost to rest frame of the cluster
	+ Boost bv = Pcl.boostVector();
	+
	+
	+ //upper scale of the quark splitting
	+ Qqtilde = _Qqtilde;
	+
	+
	+while(repeat2){
	+
	+ //choose from which constituent we radiate the gluon
	+ if ( UseRandom::rnd(0.0, 1.0)<0.5 ){
	+ P1 = ptrP->momentum();
	+ P2 = ptrPbar->momentum();
	+ mq = ptrP->data().constituentMass();
	+ mq2 = ptrPbar->data().constituentMass();
	+ fromcolor=true;
	+ }
	+ else{
	+ P2 = ptrP->momentum();
	+ P1 = ptrPbar->momentum();
	+ mq = ptrPbar->data().constituentMass();
	+ mq2 = ptrP->data().constituentMass();
	+ fromcolor=false;
	+ }
	+
	+
	+
	+ //and check if it is kinematically possible for that constituent
	+ if (!canSplitMinimally(Mcl,mq,mq2,m0)){
	+ swap(P1,P2);
	+ swap(mq,mq2);
	+ fromcolor=(!fromcolor);
	+ }
	+
	+
	+ P1.boost(-bv);
	+ P2.boost(-bv);
	+
	+
	+
	+ //Qqtilde = Qqtilde + mq;
	+ //get z and qtilde
	+ double zmin=mq/Qqtilde;
	+ double zmax=1-(4*sqr(m0)/( sqr(Mcl-mq2)-sqr(mq) ));
	+ //if ((_AngOrdFission==1) && (1-(4*m0/Qqtilde)<zmax)){
	+ // zmax=1-(4*m0/Qqtilde);
	+// }
	+
	+
	+
	+
	+ double rzmin=-log(1-zmin);
	+ double rzmax=-log(1-zmax);
	+
	+ bool repeat = true;
	+ while (repeat){
	+ //generate z and qtilde from a proposal distribution and check that at least with minimal quark mass the kinematics can work out
	+ bool repeat3 = true;
	+ double ymin, ymax;
	+ while (repeat3){
	+ z = 1.0 - exp(-UseRandom::rnd(rzmin,rzmax));
	+
	+ ymin =mq/(z*Qqtilde);
	+ if ((_AngOrdFission==1)&&(4m0/((1.-z)Qqtilde)>ymin)){
	+ ymin=4m0/((1.-z)Qqtilde);
	+ }
	+ ymax =sqrt( ( sqr(Mcl-mq2) -sqr(mq) -(4sqr(m0)/(1-z)) )/( z(1-z) ) )/Qqtilde;
	+ if (ymax >1.){
	+ ymax=1.;
	+ }
	+ if (ymin<ymax){
	+ repeat3=false;
	+ }
	+ }
	+
	+ double rymin=-log(ymax);
	+ double rymax=-log(ymin);
	+
	+ qtilde=Qqtilde*exp(-UseRandom::rnd(rymin,rymax));
	+
	+ //compare with actual distrubition
	+ if (Pqzproposal(z,qtilde,Qqtilde,mq)*UseRandom::rnd(0.0,1.0)<Pqz(z,qtilde,mq)){
	+ repeat = false;
	+ }
	+ }
	+
	+
	+ //scale for the Sudakov in the gluon splitting
	+ if (_AngOrdFission==1){
	+ Qgtilde=(1.-z)*qtilde;
	+ }
	+ else{
	+ Qgtilde=_Qg2tilde;
	+ }
	+
	+
	+ //gluon virtuality
	+ if (_restrictGluon==1){
	+ //restirct the gluon virtuality to the range where kinematic reconstruction is possible, i.e. Mcl>MP1+mq2 is always true
	+ Energy mgmax=sqrt((1.-z)(sqr(Mcl-mq2) - sqr(mq) -z(1.-z)*sqr(qtilde)));
	+ mg=_dynamicGluonMassGenerator->generate(Qgtilde,mgmax);
	+ }
	+ else{
	+ //no additional restriction on gluon virtuality. Mcl>MP1+mq2 might fail and new z and qtilde for the quark are generated
	+ mg=_dynamicGluonMassGenerator->generate(Qgtilde);
	+ }
	+
	+
	+
	+
	+
	+
	+
	+ //kinematic reconstruction / get momenta P1prime and P2prime before radiating
	+ //invariant mass of quark before radiating
	+ MP1 = sqrt(sqr(mq)+z(1.0-z)sqr(qtilde)+(sqr(mg)/(1.-z)));
	+
	+
	+
	+ if(Mcl>MP1+mq2){
	+ repeat2 = false;
	+ }
	+
	+ }//end while(repeat2)
	+
	+
	+
	+
	+ Energy absP = sqrt(sqr(sqr(Mcl))+sqr(sqr(MP1))+sqr(sqr(mq2))-2.0(sqr(MclMP1)+sqr(Mclmq2)+sqr(MP1mq2)))/(2.0*Mcl);
	+
	+
	+
	+ Lorentz5Momentum P1prime(absP*P1.vect().unit(),sqrt(sqr(absP)+sqr(MP1)),MP1);
	+ Lorentz5Momentum P2prime(absP*P2.vect().unit(),sqrt(sqr(absP)+sqr(mq2)),mq2);
	+
	+
	+
	+
	+
	+ //splitting
	+ //construct the backwards light-like vector
	+ LorentzVector<double> nbar(-P1prime.vect().unit(),1.0);
	+
	+ //construct the transverse momentum with random azimuthal angle
	+ //construct an orthonormal basis of transverse 4-vectors
	+ ThreeVector<double> Three_qperp1 = nbar.vect().orthogonal().unit();
	+ ThreeVector<double> Three_qperp2 = nbar.vect().unit().cross(Three_qperp1);
	+ LorentzVector<double> qperp1(Three_qperp1,0.0);
	+ LorentzVector<double> qperp2(Three_qperp2,0.0);
	+
	+
	+
	+ using Constants::pi;
	+ double phi = UseRandom::rnd( -pi , pi );
	+ Energy gperpabs=sqrt(sqr(1.-z)(sqr(zqtilde)-sqr(mq)) );
	+
	+ LorentzMomentum gperp = gperpabscos(phi)qperp1+gperpabssin(phi)qperp2;
	+
	+ //quark and gluon momenta
	+ Lorentz5Momentum k(zP1prime + ((sqr(mq)-sqr(zMP1)+sqr(gperpabs))/(2.0zP1prime.dot(nbar)))*nbar - gperp,mq);
	+ Lorentz5Momentum g((1.-z)P1prime + ((sqr(mg)-sqr((1.-z)MP1)+sqr(gperpabs))/(2.0(1.-z)P1prime.dot(nbar)))*nbar +gperp,mg);
	+
	+
	+
	+
	+ //create particles and give them to PartonSplitter
	+
	+
	+ PPtr ptrPprime = PPtr();
	+ PPtr ptrG = PPtr();
	+ ptrG = getParticle(ThePEG::ParticleID::g);
	+ ptrG->set5Momentum(g);
	+ if (fromcolor){
	+ ptrP->set5Momentum(k);
	+ ptrPbar->set5Momentum(P2prime);
	+ ptrPprime = getParticle(ptrP->id());
	+ ptrPprime->set5Momentum(P1prime);
	+ //_dynamicPartonSplitter->splitTimeLikeGluon(ptrG,ptrPprime,ptrPbar,ptrQ,ptrQbar,Qgtilde,nbar,fromcolor);
	+ _dynamicPartonSplitter->splitTimeLikeGluon(ptrG,ptrQ, ptrQbar, ptrPprime, ptrPbar,Qgtilde,nbar,fromcolor);
	+ }
	+ else{
	+ ptrP->set5Momentum(P2prime);
	+ ptrPbar->set5Momentum(k);
	+ ptrPprime = getParticle(ptrPbar->id());
	+ ptrPprime->set5Momentum(P1prime);
	+ //_dynamicPartonSplitter->splitTimeLikeGluon(ptrG,ptrP,ptrPprime,ptrQ,ptrQbar,Qgtilde,nbar,fromcolor);
	+ _dynamicPartonSplitter->splitTimeLikeGluon(ptrG,ptrQ, ptrQbar, ptrP, ptrPprime,Qgtilde,nbar,fromcolor);
	+ }
	+
	+
	+
	+
	+ //boost back to lab frame
	+ //(this must be possible in a nicer way)
	+ Lorentz5Momentum temp;
	+ temp = ptrP->momentum();
	+ temp.boost(bv);
	+ ptrP->set5Momentum(temp);
	+
	+ temp = ptrPbar->momentum();
	+ temp.boost(bv);
	+ ptrPbar->set5Momentum(temp);
	+
	+ temp = ptrQ->momentum();
	+ temp.boost(bv);
	+ ptrQ->set5Momentum(temp);
	+
	+ temp = ptrQbar->momentum();
	+ temp.boost(bv);
	+ ptrQbar->set5Momentum(temp);
	+
	+}
	+
	+
	+bool DynamicClusterFissioner::canSplitMinimally(Energy Mcl, Energy m1, Energy m2, Energy m0) const {
	+
	+ if (Mcl<m1+m2+2*m0){ return false; }
	+
	+ Energy Qtilde=_Qqtilde;
	+
	+
	+ if (Qtilde>m1){
	+
	+ if (_AngOrdFission==0){
	+ if (Qtilde>2*m0+m1){
	+ return Mcl>m2+m1+2*m0;
	+ }
	+ else{
	+ return Mcl>m2+sqrt(Qtilde(m1- (4sqr(m0)/(Qtilde-m1)) ));
	+ }
	+ }
	+ else{
	+ if (Qtilde>4*m0+m1){
	+ return Mcl>m2+sqrt((m0+m1)(4m0+m1));
	+ }
	+ else{
	+ return Mcl<m2+sqrt( sqr(m1) + (4sqr(m0)(4*m1+Qtilde))/(Qtilde-m1) );
	+ }
	+ }
	+
	+ }
	+ else{
	+ return false;
	+ }
	+
	+
	+}
	+
	+bool DynamicClusterFissioner::canSplitMinimally(tcClusterPtr clu, Energy minmass) const {
	+ Energy Mcl=clu->mass();
	+ Energy m1 = clu->particle(0)->data().constituentMass();
	+ Energy m2 = clu->particle(1)->data().constituentMass();
	+ return ((canSplitMinimally(Mcl,m1,m2,minmass))\|\|(canSplitMinimally(Mcl,m2,m1,minmass)));
	+}
	diff --git a/Hadronization/DynamicClusterFissioner.h b/Hadronization/DynamicClusterFissioner.h
	new file mode 100644
	--- /dev/null
	+++ b/Hadronization/DynamicClusterFissioner.h
	@@ -0,0 +1,210 @@
	+// -- C++ --
	+//
	+// ClusterFissioner.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	+// Copyright (C) 2002-2019 The Herwig Collaboration
	+//
	+// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	+//
	+#ifndef HERWIG_DynamicClusterFissioner_H
	+#define HERWIG_DynamicClusterFissioner_H
	+
	+#include <ThePEG/Interface/Interfaced.h>
	+#include "CluHadConfig.h"
	+#include "ClusterFissioner.h"
	+#include "DynamicPartonSplitter.h"
	+#include "DynamicGluonMassGenerator.h"
	+
	+namespace Herwig {
	+using namespace ThePEG;
	+
	+
	+class DynamicClusterFissioner: public ClusterFissioner {
	+
	+public:
	+
	+ /** @name Standard constructors and destructors. */
	+ //@{
	+ /**
	+ * Default constructor.
	+ */
	+ DynamicClusterFissioner();
	+
	+ //@}
	+
	+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.
	+ */
	+ DynamicClusterFissioner & operator=(const DynamicClusterFissioner &) = delete;
	+
	+
	+
	+public:
	+
	+ //@{
	+
	+ /**
	+ * Split three-component cluster
	+ */
	+ virtual cutType cutThree(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn){
	+ assert(false);//implementation missing
	+ }
	+ //@}
	+
	+protected:
	+
	+
	+
	+
	+ InvEnergy Pqzproposal(double z, Energy qtilde, Energy Qqtilde, Energy mq) const;
	+ InvEnergy Pqz(double z, Energy Qqtilde, Energy mq) const;
	+
	+
	+ /**
	+ * The actual implementation of the dynamic fission
	+ */
	+ virtual void dynamicFission(tPPtr ptrP, tPPtr ptrPbar, PPtr & ptrQ, PPtr & ptrQbar) const;
	+
	+
	+ /**
	+ * call the dynamic fission function here
	+ */
	+ virtual void drawNewFlavour(PPtr& newPtr1, PPtr& newPtr2,const ClusterPtr& cluster) const {
	+ if(!cluster->particle(0)->hasAntiColour()){
	+ dynamicFission(cluster->particle(0), cluster->particle(1) ,newPtr2,newPtr1);
	+ }
	+ else{
	+ dynamicFission(cluster->particle(1), cluster->particle(0) ,newPtr1,newPtr2);
	+ }
	+ }
	+
	+
	+
	+ bool canSplitMinimally(tcClusterPtr clu, Energy minmass) const;
	+ bool canSplitMinimally(Energy Mcl, Energy m1, Energy m2, Energy m0) const;
	+
	+
	+ /**
	+ * set the momenta here
	+ */
	+ virtual pair<Energy,Energy> drawNewMasses(Energy Mc, bool soft1, bool soft2,
	+ Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
	+ tPPtr ptrQ1, Lorentz5Momentum& pQ1,
	+ tPPtr newPtr1, Lorentz5Momentum& pQone,
	+ tPPtr newPtr2, Lorentz5Momentum& pQtwo,
	+ tPPtr ptrQ2, Lorentz5Momentum& pQ2) const
	+ {
	+ pQ1 = ptrQ1->momentum();
	+ pQ2 = ptrQ2->momentum();
	+ pQone = newPtr1->momentum();
	+ pQtwo = newPtr2->momentum();
	+ pClu1 = pQ1+pQone;
	+ pClu2 = pQ2 + pQtwo;
	+
	+ pClu1.setMass(pClu1.m());
	+ pClu2.setMass(pClu2.m());
	+
	+ pQ1.setMass(ptrQ1->data().constituentMass());
	+ pQ2.setMass(ptrQ2->data().constituentMass());
	+ pQone.setMass(newPtr1->data().constituentMass());
	+ pQtwo.setMass(newPtr2->data().constituentMass());
	+
	+ pair<Energy,Energy> result;
	+ result.first = pClu1.m();
	+ result.second = pClu2.m();
	+ return result;
	+ }
	+
	+ /**
	+ * does nothing
	+ * kinematics already done in drawNewMasses
	+ */
	+ virtual void calculateKinematics(const Lorentz5Momentum &pClu,
	+ const Lorentz5Momentum &p0Q1,
	+ const bool toHadron1, const bool toHadron2,
	+ Lorentz5Momentum &pClu1, Lorentz5Momentum &pClu2,
	+ Lorentz5Momentum &pQ1, Lorentz5Momentum &pQb,
	+ Lorentz5Momentum &pQ2, Lorentz5Momentum &pQ2b) const {}
	+
	+
	+
	+ /**
	+ * pointer to the dynamic parton splitter
	+ */
	+ /*Ptr<DynamicPartonSplitter>::tptr partonSplitter() const {
	+ Ptr<DynamicPartonSplitter>::tptr partonSplitter_ = dynamic_ptr_cast<Ptr<DynamicPartonSplitter>::tptr>(ClusterHadronizationHandler::currentHandler()->partonSplitter());
	+ return partonSplitter_;
	+ }*/
	+
	+ /**
	+ * pointer to the dynamic gluon mass generator
	+ */
	+ /*Ptr<DynamicGluonMassGenerator>::tptr gluonMassGenerator() const {
	+ Ptr<DynamicGluonMassGenerator>::tptr gluonMassGenerator_ = dynamic_ptr_cast<Ptr<DynamicGluonMassGenerator>::tptr>(ClusterHadronizationHandler::currentHandler()->gluonMassGenerator());
	+ return gluonMassGenerator_;
	+ }*/
	+
	+
	+private:
	+
	+ //scales for quark and gluon splitting in the fission
	+ Energy _Qqtilde;
	+ Energy _Qg2tilde;
	+
	+ //switch for angular ordering in dynamic fission
	+ int _AngOrdFission;
	+
	+ //switch for restricting the gluon virtuality in dynamic fission
	+ int _restrictGluon;
	+
	+ Ptr<DynamicGluonMassGenerator>::ptr _dynamicGluonMassGenerator;
	+ Ptr<DynamicPartonSplitter>::ptr _dynamicPartonSplitter;
	+
	+
	+
	+};
	+}
	+
	+#endif /* HERWIG_DynamicClusterFissioner_H */
	diff --git a/Hadronization/DynamicGluonMassGenerator.cc b/Hadronization/DynamicGluonMassGenerator.cc
	new file mode 100644
	--- /dev/null
	+++ b/Hadronization/DynamicGluonMassGenerator.cc
	@@ -0,0 +1,153 @@
	+// -- C++ --
	+//
	+// This is the implementation of the non-inlined, non-templated member
	+// functions of the GluonMassGenerator class.
	+//
	+
	+#include "DynamicGluonMassGenerator.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"
	+#include <ThePEG/Interface/Parameter.h>
	+
	+using namespace Herwig;
	+
	+DynamicGluonMassGenerator::DynamicGluonMassGenerator() {}
	+
	+DynamicGluonMassGenerator::~DynamicGluonMassGenerator() {}
	+
	+IBPtr DynamicGluonMassGenerator::clone() const {
	+ return new_ptr(*this);
	+}
	+
	+IBPtr DynamicGluonMassGenerator::fullclone() const {
	+ return new_ptr(*this);
	+}
	+
	+
	+
	+InvEnergy DynamicGluonMassGenerator::PmgProposal(Energy, Energy mq) const {
	+ return (ClusterAlphaS(4mqmq)*0.32/mq);
	+}
	+
	+
	+InvEnergy DynamicGluonMassGenerator::PmgProposal(Energy mg) const {
	+ Energy mu, ms, md;
	+ mu=getParticleData(ThePEG::ParticleID::u)->constituentMass();
	+ md=getParticleData(ThePEG::ParticleID::d)->constituentMass();
	+ ms=getParticleData(ThePEG::ParticleID::s)->constituentMass();
	+ return PmgProposal(mg, mu)+PmgProposal(mg, md)+PmgProposal(mg, ms);
	+}
	+
	+
	+
	+InvEnergy DynamicGluonMassGenerator::Pmg(Energy mg, Energy mq, Energy Qtilde) const {
	+ if((2mq<mg) && (mg <= sqrt(mqQtilde))){
	+ return (ClusterAlphaS(mgmg)/mg)sqrt(1-pow(2mq/mg,2))(1+(2*pow(mq/mg,2)));
	+ }
	+ else if ((sqrt(mq*Qtilde)<mg) &&(mg<=Qtilde/2.0)){
	+ return (ClusterAlphaS(mgmg)/mg)sqrt(1-pow(2mg/Qtilde,2))(1+(3*pow(mq/mg,2))-pow(mg/Qtilde,2) );
	+ }
	+ else{
	+ return 0.0*InvGeV;
	+ }
	+}
	+
	+InvEnergy DynamicGluonMassGenerator::Pmg(Energy mg,Energy Qtilde) const {
	+ Energy mu, ms, md;
	+ mu=getParticleData(ThePEG::ParticleID::u)->constituentMass();
	+ md=getParticleData(ThePEG::ParticleID::d)->constituentMass();
	+ ms=getParticleData(ThePEG::ParticleID::s)->constituentMass();
	+ return Pmg(mg,mu,Qtilde)+Pmg(mg,md,Qtilde)+Pmg(mg,ms,Qtilde);
	+}
	+
	+
	+
	+Energy DynamicGluonMassGenerator::generateProposal(Energy mgmax) const {
	+ Energy m0, mu, ms, md;
	+ 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;}
	+ return UseRandom::rnd(2.*m0,mgmax);
	+}
	+
	+Energy DynamicGluonMassGenerator::generate(Energy Qtilde, Energy mgmax) const {
	+ Energy mg;
	+ double r;
	+ bool repeat;
	+ Energy max=mgmax;
	+ if(0.5Qtilde<max){max=0.5Qtilde;};
	+ repeat=true;
	+ while(repeat){
	+ mg=generateProposal(max);
	+ r=UseRandom::rnd(0.0, 1.0);
	+ if(r*PmgProposal(mg)<Pmg(mg,Qtilde)){
	+ repeat=false;
	+ }
	+ }
	+ return mg;
	+}
	+
	+
	+Energy DynamicGluonMassGenerator::generate(Energy Qtilde) const {
	+ return generate(Qtilde,Qtilde);
	+}
	+
	+
	+Energy DynamicGluonMassGenerator::generate() const {
	+ return generate(Qgtilde());
	+}
	+
	+
	+
	+// If needed, insert default implementations of virtual function defined
	+// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
	+
	+
	+
	+void DynamicGluonMassGenerator::persistentOutput(PersistentOStream & os) const {
	+ // * ATTENTION * os << ; // Add all member variable which should be written persistently here.
	+ os << ounit(_Qgtilde,GeV) << ounit(_clusteralphasfreeze,GeV);
	+}
	+
	+void DynamicGluonMassGenerator::persistentInput(PersistentIStream & is, int) {
	+ // * ATTENTION * is >> ; // Add all member variable which should be read persistently here.
	+ is >> iunit(_Qgtilde,GeV) >> iunit(_clusteralphasfreeze,GeV);
	+}
	+
	+
	+// * 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<DynamicGluonMassGenerator,GluonMassGenerator>
	+ describeHerwigDynamicGluonMassGenerator("Herwig::DynamicGluonMassGenerator", "");
	+
	+void DynamicGluonMassGenerator::Init() {
	+
	+ static ClassDocumentation<DynamicGluonMassGenerator> documentation
	+ ("There is no documentation for the DynamicGluonMassGenerator class");
	+
	+
	+ static Parameter<DynamicGluonMassGenerator,Energy> interfaceQgtilde
	+ ("Qgtilde",
	+ "Upper scale for the Sudakov for the non-pert. gluon splitting",
	+ &DynamicGluonMassGenerator::_Qgtilde, GeV, 4.0GeV, 0.5GeV, 10000.0*GeV,false,false,false);
	+
	+ static Parameter<DynamicGluonMassGenerator,Energy> interfaceClusterAlphaSFreeze
	+ ("ClusterAlphaSFreeze",
	+ "Freeze-out scale for the non-pert. coupling in gluon splitting and cluster fission",
	+ &DynamicGluonMassGenerator::_clusteralphasfreeze, GeV, 1.0GeV, 0.1GeV, 10.0*GeV,false,false,false);
	+
	+}
	+
	diff --git a/Hadronization/DynamicGluonMassGenerator.h b/Hadronization/DynamicGluonMassGenerator.h
	new file mode 100644
	--- /dev/null
	+++ b/Hadronization/DynamicGluonMassGenerator.h
	@@ -0,0 +1,174 @@
	+// -- C++ --
	+#ifndef Herwig_DynamicGluonMassGenerator_H
	+#define Herwig_DynamicGluonMassGenerator_H
	+//
	+// This is the declaration of the DynamicGluonMassGenerator class.
	+//
	+
	+
	+
	+#include "Herwig/Hadronization/GluonMassGenerator.h"
	+#include "ClusterHadronizationHandler.h"
	+#include "ThePEG/StandardModel/StandardModelBase.h"
	+#include "Cluster.h"
	+
	+
	+namespace Herwig {
	+
	+using namespace ThePEG;
	+
	+/**
	+ * Here is the documentation of the GluonMassGenerator class.
	+ *
	+ * @see \ref GluonMassGeneratorInterfaces "The interfaces"
	+ * defined for GluonMassGenerator.
	+ */
	+class DynamicGluonMassGenerator: public GluonMassGenerator {
	+
	+public:
	+
	+ /** @name Standard constructors and destructors. */
	+ //@{
	+ /**
	+ * The default constructor.
	+ */
	+ DynamicGluonMassGenerator();
	+
	+ /**
	+ * The destructor.
	+ */
	+ virtual ~DynamicGluonMassGenerator();
	+ //@}
	+
	+public:
	+
	+
	+ /**
	+ * Returns the scale Qgtilde of the Sudakov in the gluon splitting
	+ */
	+ virtual Energy Qgtilde() const {
	+ return _Qgtilde;
	+ }
	+
	+ /**
	+ * Returns the minimum mass of a gluon that can by generated
	+ */
	+ virtual Energy minGluonMass() const { //ToDo: get it from hadron spectrum
	+ Energy m0, mu, md, ms;
	+ 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;}
	+ return m0;
	+ }
	+
	+
	+ /**
	+ * Generate a gluon mass drawn from the proposal distribution (overestimate)
	+ */
	+ virtual Energy generateProposal(Energy Qtilde) const;
	+
	+ /**
	+ * The proposal distributions for the gluon mass (overestimate)
	+ */
	+ virtual InvEnergy PmgProposal(Energy mg, Energy mq) const;
	+ virtual InvEnergy PmgProposal(Energy mg) const;
	+
	+
	+ /**
	+ * The gluon mass distribution
	+ */
	+ virtual InvEnergy Pmg(Energy mg, Energy mq, Energy Qtilde) const;
	+ virtual InvEnergy Pmg(Energy mg, Energy Qtilde) const;
	+
	+
	+ /**
	+ * Generate a single gluon mass
	+ */
	+ virtual Energy generate(Energy Qtilde, Energy mgmax) const;
	+ virtual Energy generate(Energy Qtilde) const;
	+ virtual Energy generate() const;
	+
	+
	+ /**
	+ * Return the the strong coupling used in ClusterFission
	+ */
	+ double ClusterAlphaS(Energy2 q2) const {
	+ if(q2>sqr(_clusteralphasfreeze)){return ClusterHadronizationHandler::currentHandler()->SM().alphaS(q2);}
	+ else{return ClusterHadronizationHandler::currentHandler()->SM().alphaS(sqr(_clusteralphasfreeze));}
	+ }
	+
	+
	+
	+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:
	+
	+ /**
	+ * Private and non-existent assignment operator.
	+ */
	+ DynamicGluonMassGenerator & operator=(const DynamicGluonMassGenerator &) = delete;
	+
	+ /**
	+ * Scale of the Sudakov in the gluon splitting
	+ */
	+ Energy _Qgtilde = 4.0*GeV;
	+
	+ /**
	+ * freezing scale for the non-pert. alphas
	+ */
	+ Energy _clusteralphasfreeze = 1.0*GeV;
	+
	+};
	+
	+}
	+
	+#endif /* Herwig_DynamicGluonMassGenerator_H */
	diff --git a/Hadronization/DynamicPartonSplitter.cc b/Hadronization/DynamicPartonSplitter.cc
	new file mode 100644
	--- /dev/null
	+++ b/Hadronization/DynamicPartonSplitter.cc
	@@ -0,0 +1,312 @@
	+// -- C++ --
	+//
	+// PartonSplitter.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 PartonSplitter class.
	+//
	+
	+#include "DynamicPartonSplitter.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 <ThePEG/EventRecord/Particle.h>
	+#include <ThePEG/PDT/PDT.h>
	+
	+
	+using namespace Herwig;
	+
	+IBPtr DynamicPartonSplitter::clone() const {
	+ return new_ptr(*this);
	+}
	+
	+IBPtr DynamicPartonSplitter::fullclone() const {
	+ return new_ptr(*this);
	+}
	+
	+void DynamicPartonSplitter::persistentOutput(PersistentOStream & os) const {
	+ os << _findProgenitor << _restrictZ << _dynamicGluonMassGenerator;
	+}
	+
	+void DynamicPartonSplitter::persistentInput(PersistentIStream & is, int) {
	+ is >> _findProgenitor >> _restrictZ >> _dynamicGluonMassGenerator;
	+}
	+
	+DescribeClass<DynamicPartonSplitter,PartonSplitter>
	+describeDynamicPartonSplitter("Herwig::DynamicPartonSplitter","Herwig.so");
	+
	+
	+void DynamicPartonSplitter::Init() {
	+
	+ static ClassDocumentation<DynamicPartonSplitter> documentation
	+ ("This class is reponsible of the nonperturbative splitting of partons");
	+
	+
	+ static Switch<DynamicPartonSplitter,int> interfaceRestrictZ
	+ ("RestrictZ",
	+ "Option to restrict z in gluon splitting to z > 0.5",
	+ &DynamicPartonSplitter::_restrictZ,0,false,false);
	+ static SwitchOption interfaceRestrictZNo
	+ (interfaceRestrictZ,
	+ "No",
	+ "No restriction",
	+ 0);
	+ static SwitchOption interfaceRestrictZYes
	+ (interfaceRestrictZ,
	+ "Yes",
	+ "Restric z > 0.5",
	+ 1);
	+
	+
	+ static Switch<DynamicPartonSplitter,int> interfaceFindProgenitor
	+ ("FindProgenitor",
	+ "Option for how to choose which color partner is the progenitor of the gluon",
	+ &DynamicPartonSplitter::_findProgenitor,0,false,false);
	+ static SwitchOption interfaceFindProgenitorAngle
	+ (interfaceFindProgenitor,
	+ "Angle",
	+ "choose color partner with maximum cos(theta)",
	+ 0);
	+ static SwitchOption interfaceFindProgenitorMomentum
	+ (interfaceFindProgenitor,
	+ "Momentum",
	+ "choose color partner with maximum \|p\|*cos(theta)",
	+ 1);
	+ static SwitchOption interfaceFindProgenitorpT
	+ (interfaceFindProgenitor,
	+ "pT",
	+ "choose color partner with maximum cos^2(theta) as backwards direction",
	+ 2);
	+ static SwitchOption interfaceFindProgenitorpTforward
	+ (interfaceFindProgenitor,
	+ "pTforward",
	+ "choose color partner with maximum cos^2(theta), gluon always in forward direction",
	+ 3);
	+ static SwitchOption interfaceFindProgenitorGluon
	+ (interfaceFindProgenitor,
	+ "Gluon",
	+ "Gluon as progenitor",
	+ 4);
	+ static SwitchOption interfaceFindProgenitorIsotropic
	+ (interfaceFindProgenitor,
	+ "Isotropic",
	+ "Isotropic gluon decay",
	+ 5);
	+
	+ static Reference<DynamicPartonSplitter,DynamicGluonMassGenerator> interfaceDynamicGluonMassGenerator
	+ ("GluonMassGenerator",
	+ "Set a reference to a gluon mass generator.",
	+ &DynamicPartonSplitter::_dynamicGluonMassGenerator, false, false, true, true, false);
	+}
	+
	+
	+void DynamicPartonSplitter::drawNewFlavour(PPtr & ptrQ, PPtr & ptrQbar, Energy mg, Energy Qtilde){
	+ static const Energy md = getParticleData(ThePEG::ParticleID::d)->constituentMass();
	+ static const Energy mu = getParticleData(ThePEG::ParticleID::u)->constituentMass();
	+ static const Energy ms = getParticleData(ThePEG::ParticleID::s)->constituentMass();
	+
	+ double prob_d=_dynamicGluonMassGenerator->Pmg(mg,md,Qtilde)*UnitRemoval::E;
	+ double prob_u=_dynamicGluonMassGenerator->Pmg(mg,mu,Qtilde)*UnitRemoval::E;
	+ double prob_s=_dynamicGluonMassGenerator->Pmg(mg,ms,Qtilde)*UnitRemoval::E;
	+
	+
	+ long idNew=0;
	+ 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);
	+ assert (ptrQ);
	+ assert(ptrQbar);
	+ assert (ptrQ->dataPtr());
	+ assert(ptrQbar->dataPtr());
	+
	+}
	+
	+
	+void DynamicPartonSplitter::splitTimeLikeGluon(tcPPtr ptrGluon, PPtr & ptrQ, PPtr & ptrQbar, tcPPtr ptrP, tcPPtr ptrPbar, Energy Qtilde, LorentzVector<double> nbar, bool ProgenitorHasColor){
	+
	+ LorentzMomentum momentumG = ptrGluon->momentum();
	+ LorentzMomentum momentumP = ptrP->momentum();
	+ LorentzMomentum momentumPbar = ptrPbar->momentum();
	+ Energy mg = momentumG.m();
	+
	+ drawNewFlavour(ptrQ,ptrQbar,mg,Qtilde);
	+
	+ Energy mQ = ptrQ->data().constituentMass();
	+ LorentzMomentum momentumQ, momentumQbar;
	+
	+ if(_findProgenitor == 5){
	+
	+ Lorentz5Momentum momentum5Qtemp, momentum5Qbartemp;
	+ double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 );
	+ using Constants::pi;
	+ double phiStar = UseRandom::rnd( -pi , pi );
	+ Kinematics::twoBodyDecay(momentumG, mQ,
	+ mQ, cosThetaStar, phiStar, momentum5Qtemp,
	+ momentum5Qbartemp );
	+ momentumQ = momentum5Qtemp;
	+ momentumQbar = momentum5Qbartemp;
	+ }
	+ else{
	+ //construct an orthonormal basis of transverse 4-vectors
	+ ThreeVector<double> Three_qperp1 = nbar.vect().orthogonal().unit();
	+ ThreeVector<double> Three_qperp2 = nbar.vect().unit().cross(Three_qperp1);
	+ LorentzVector<double> qperp1(Three_qperp1,0.0);
	+ LorentzVector<double> qperp2(Three_qperp2,0.0);
	+
	+ //get the z-value from the gluon splitting function
	+ bool repeat;
	+ double z,r,x;
	+ if (mg < sqrt(mQ*Qtilde)){x = pow(mQ/mg,2);}
	+ else {x = pow(mg/Qtilde,2);}
	+ repeat = true;
	+ while(repeat){
	+ r = UseRandom::rnd(0.0, 1.0);
	+ z = (r-0.5)sqrt(1-4x)+0.5;
	+ r = UseRandom::rnd(0.0, 1.0-2.0x+2.0pow(mQ/mg,2));
	+ if( r < (1-(2z(1-z))+2*pow(mQ/mg,2)) ){
	+ repeat = false;
	+ }
	+ }
	+
	+ if (_restrictZ==1){//without this ProgenitorHasColor would not be necessary
	+ if (ProgenitorHasColor){
	+ if (z>0.5){z = 1.-z;} //anti-quark in direction of progenitor
	+ }
	+ else{
	+ if (z<0.5){z = 1.-z;} //quark in direction of progenitor
	+ }
	+ }
	+
	+ //magniuted of relative transverse momentum
	+ Energy qT = sqrt(sqr(mg)z(1.-z)-sqr(mQ));
	+ //azimuthal angle of relative tranverse momentum
	+ using Constants::pi;
	+ double phi = UseRandom::rnd( -pi , pi );
	+
	+ //set the relative transverse momentum
	+ LorentzMomentum qperp = qTcos(phi)qperp1+qTsin(phi)qperp2;
	+
	+ //calculate the quarks' momenta
	+ momentumQ = zmomentumG + ( ( (mgmg)(1-2.z) -2momentumG.dot(qperp) )/( 2.momentumG.dot(nbar) ) )*nbar + qperp;
	+ momentumQbar = momentumG - momentumQ;
	+
	+ } //else _findProgenitor Isotropic
	+
	+ //without this I would not need to give also P and Pbar
	+ if (((_findProgenitor == 4)\|\|(_findProgenitor==5))&&(_restrictZ==1)){
	+ LorentzMomentum c1 = momentumQ+momentumP;
	+ LorentzMomentum c2 = momentumQbar+momentumPbar;
	+ LorentzMomentum c3 = momentumQ+momentumPbar;
	+ LorentzMomentum c4 = momentumQbar+momentumP;
	+
	+ //linear
	+ Energy M12=c1.m()+c2.m();
	+ Energy M34=c3.m()+c4.m();
	+
	+ //quadratic
	+ //Energy2 M12 = c1.m()c1.m() + c2.m()c2.m();
	+ //Energy2 M34 = c3.m()c3.m() + c4.m()c4.m();
	+
	+ //product
	+ //Energy2 M12=c1.m()*c2.m();
	+ //Energy2 M34=c3.m()*c4.m();
	+
	+ //maximum
	+ //Energy M12 = std::max(c1.m(),c2.m());
	+ //Energy M34 = std::max(c3.m(),c4.m());
	+
	+ if (M12 < M34){
	+ swap(momentumQ,momentumQbar);
	+ }
	+ }
	+
	+ Lorentz5Momentum momentum5Q(momentumQ,mQ);
	+ Lorentz5Momentum momentum5Qbar(momentumQbar,mQ);
	+
	+ ptrQ ->set5Momentum( momentum5Q );
	+ ptrQbar ->set5Momentum( momentum5Qbar );
	+
	+}
	+
	+
	+
	+void DynamicPartonSplitter::splitTimeLikeGluon(tcPPtr ptrGluon, PPtr & ptrQ, PPtr & ptrQbar){
	+ LorentzVector<double> nbar;
	+ bool ProgenitorHasColor;
	+
	+ tPPtr ColorPartner = ptrGluon->antiColourLine()->endParticle();
	+ tPPtr AntiColorPartner = ptrGluon->colourLine()->startParticle();
	+ Findnbar(ptrGluon,ColorPartner,AntiColorPartner,nbar,ProgenitorHasColor);
	+
	+ splitTimeLikeGluon(ptrGluon,ptrQ,ptrQbar,ColorPartner,AntiColorPartner,_dynamicGluonMassGenerator->Qgtilde(),nbar,ProgenitorHasColor);
	+}
	+
	+
	+
	+void DynamicPartonSplitter::Findnbar(tcPPtr ptrGluon, tcPPtr ptrP, tcPPtr ptrPbar, LorentzVector<double> & nbarfinal, bool & ProgenitorHasColor){
	+
	+ LorentzMomentum momentumG = ptrGluon->momentum();
	+ LorentzMomentum momentumP = ptrP->momentum();
	+ LorentzMomentum momentumPbar = ptrPbar->momentum();
	+
	+ LorentzMomentum nbar = momentumP;
	+ ProgenitorHasColor = true;
	+
	+ if (_findProgenitor==0){
	+ if( momentumG.vect().unit().dot( momentumPbar.vect().unit() ) > momentumG.vect().unit().dot( momentumP.vect().unit() )){
	+ nbar = momentumPbar; //minimum angle theta
	+ ProgenitorHasColor = false; //has anticolor
	+ }
	+ nbar.setVect( (-1.)*nbar.vect() );
	+ }
	+ else if (_findProgenitor==1){
	+ if ( momentumG.vect().dot( momentumPbar.vect() ) > momentumG.vect().dot( momentumP.vect() ) ){
	+ nbar = momentumPbar; //maximum \|p\|*cos(theta)
	+ ProgenitorHasColor = false; //has anticolor
	+ }
	+ nbar.setVect( (-1.)*nbar.vect() );
	+ }
	+ else if ((_findProgenitor==2)\|\|(_findProgenitor==3)){
	+ ProgenitorHasColor = false;
	+ if( pow(momentumG.vect().unit().dot( momentumPbar.vect().unit()),2) > pow(momentumG.vect().unit().dot( momentumP.vect().unit() ),2)){
	+ nbar = momentumPbar; //maximum cos^2(theta) --> minimum gluon transverse momentum
	+ ProgenitorHasColor = true;
	+ }
	+ if (_findProgenitor==3){
	+ if ( momentumG.vect().unit().dot( nbar.vect().unit() ) > 0. ){
	+ nbar.setVect( (-1.)*nbar.vect() ); //gluon always in forward direction
	+ ProgenitorHasColor = !ProgenitorHasColor;
	+ }
	+ }
	+ }
	+ else if (_findProgenitor==4){
	+ nbar = momentumG;
	+ nbar.setVect( (-1.)*nbar.vect() );
	+ }
	+
	+ nbarfinal.setT(1.0);
	+ nbarfinal.setVect(nbar.vect().unit());
	+}
	+
	+
	diff --git a/Hadronization/DynamicPartonSplitter.h b/Hadronization/DynamicPartonSplitter.h
	new file mode 100644
	--- /dev/null
	+++ b/Hadronization/DynamicPartonSplitter.h
	@@ -0,0 +1,180 @@
	+// -- C++ --
	+//
	+// PartonSplitter.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_DynamicPartonSplitter_H
	+#define HERWIG_DynamicPartonSplitter_H
	+
	+#include "CluHadConfig.h"
	+#include <ThePEG/Interface/Interfaced.h>
	+#include <ThePEG/Utilities/Selector.h>
	+#include "PartonSplitter.h"
	+#include "DynamicGluonMassGenerator.h"
	+
	+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.
	+ *
	+ * 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 DynamicPartonSplitter: public PartonSplitter {
	+
	+public:
	+
	+ /**
	+ * Default constructor
	+ */
	+ DynamicPartonSplitter() :
	+ PartonSplitter(),
	+ _findProgenitor(0),
	+ _restrictZ(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);
	+ //@}
	+
	+ /**
	+ * 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;
	+ //@}
	+
	+private:
	+
	+ /**
	+ * Private and non-existent assignment operator.
	+ */
	+ DynamicPartonSplitter & operator=(const DynamicPartonSplitter &) = delete;
	+
	+
	+public:
	+ /**
	+ * 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
	+ */
	+ virtual void splitTimeLikeGluon(tcPPtr gluon, PPtr & quark, PPtr & anti);
	+
	+
	+ /**
	+ * 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
	+ * @colorpartner1 The color partner of the gluon
	+ * @colorpartner2 The anticolor partner of the gluon
	+ * @Qtilde The scale of the Sudakov in the splitting
	+ * @nbar Backwards lightlike direction of the splitting
	+ * @ProgenitorHasColor If the (anti)quark that is identified as the "progenitor" of the gluon has color or anticolor
	+ * @
	+ */
	+ void splitTimeLikeGluon(tcPPtr gluon, PPtr & quark, PPtr & anti, tcPPtr colorpartner1, tcPPtr colorpartner2, Energy Qtilde, LorentzVector<double> nbar, bool ProgenitorHasColor);
	+
	+private:
	+ /**
	+ * Finding the backwards light like direction for the splitting and defining the gluon "progenitor"
	+ * @gluon The gluon to be split
	+ * @colorpartner1 The color partner of the gluon
	+ * @colorpartner2 The anticolor partner of the gluon
	+ * @nbar The lightlike vector that will be defined
	+ * @ProgenitorHasColor Returns whether the "progenitor" of the gluon has color or anticolor
	+ * */
	+ void Findnbar(tcPPtr gluon, tcPPtr colorpartner1, tcPPtr colorpartner2, LorentzVector<double> & nbar, bool & ProgenitorHasColor);
	+
	+
	+ /**
	+ * draw the new light flavour according to probabilities given by gluon mass distribution
	+ */
	+ void drawNewFlavour(PPtr & ptrQ, PPtr & ptrQbar, Energy mg, Energy Qtilde);
	+
	+
	+
	+
	+private:
	+
	+ /**
	+ * switch for different options for choosing the gluon's progenitor
	+ */
	+ int _findProgenitor;
	+
	+
	+ /**
	+ * switch for restricting z > 0.5
	+ */
	+ int _restrictZ;
	+
	+
	+ /**
	+ * pointer to the dynamic gluon mass generator
	+ */
	+ Ptr<DynamicGluonMassGenerator>::ptr _dynamicGluonMassGenerator;
	+
	+
	+};
	+
	+}
	+
	+#endif /* HERWIG_DynamicPartonSplitter_H */
	diff --git a/Hadronization/Makefile.am b/Hadronization/Makefile.am
	--- a/Hadronization/Makefile.am
	+++ b/Hadronization/Makefile.am
	@@ -1,23 +1,26 @@
	noinst_LTLIBRARIES = libHwHadronization.la
	libHwHadronization_la_SOURCES = \
	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 \
	MatrixElementClusterFissioner.cc MatrixElementClusterFissioner.h \
	ClusterHadronizationHandler.cc ClusterHadronizationHandler.h \
	ClusterHadronizationHandler.fh \
	ColourReconnector.cc ColourReconnector.h ColourReconnector.fh \
	+DynamicGluonMassGenerator.h DynamicGluonMassGenerator.cc \
	+DynamicPartonSplitter.h DynamicPartonSplitter.cc \
	+DynamicClusterFissioner.h DynamicClusterFissioner.cc \
	GluonMassGenerator.h GluonMassGenerator.cc \
	Hw64Selector.cc Hw64Selector.h \
	HwppSelector.cc HwppSelector.h \
	Hw7Selector.cc Hw7Selector.h \
	LightClusterDecayer.cc LightClusterDecayer.h LightClusterDecayer.fh \
	PartonSplitter.cc PartonSplitter.h PartonSplitter.fh \
	SpinHadronizer.h SpinHadronizer.cc \
	HadronSpectrum.h HadronSpectrum.fh HadronSpectrum.cc \
	StandardModelHadronSpectrum.h StandardModelHadronSpectrum.cc \
	DarkHadronSpectrum.h DarkHadronSpectrum.cc \
	DarkHwppSelector.cc DarkHwppSelector.h

	diff --git a/Shower/QTilde/Couplings/ShowerAlphaQCDNP.cc b/Shower/QTilde/Couplings/ShowerAlphaQCDNP.cc
	--- a/Shower/QTilde/Couplings/ShowerAlphaQCDNP.cc
	+++ b/Shower/QTilde/Couplings/ShowerAlphaQCDNP.cc
	@@ -1,419 +1,422 @@
	// -- C++ --
	//
	// ShowerAlphaQCDNP.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 ShowerAlphaQCDNP class.
	//
	#include "ShowerAlphaQCDNP.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Command.h"
	#include "ThePEG/Interface/Deleted.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Config/Constants.h"
	#include "Herwig/Utilities/AlphaS.h"
	#include "gsl/gsl_sf_lambert.h"

	using namespace Herwig;
	using Herwig::Math::alphaS;
	using Herwig::Math::derivativeAlphaS;

	DescribeClass<ShowerAlphaQCDNP,ShowerAlpha>
	describeShowerAlphaQCDNP("Herwig::ShowerAlphaQCDNP","HwShower.so");

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

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

	void ShowerAlphaQCDNP::persistentOutput(PersistentOStream & os) const {
	os << ounit(_qmin,GeV) << _asMaxNP << _nfMaxNP
	<< ounit(_thresholds,GeV) << ounit(_lambda,GeV)
	<< _nloop << _thresopt
	<< _alphain << _tolerance << _maxtry
	<< _npPower << ounit(_optInputScale,GeV) << _quarkBranching
	<< ounit(_quarkMasses,GeV);
	}

	void ShowerAlphaQCDNP::persistentInput(PersistentIStream & is, int) {
	is >> iunit(_qmin,GeV) >> _asMaxNP >> _nfMaxNP
	>> iunit(_thresholds,GeV) >> iunit(_lambda,GeV)
	>> _nloop >> _thresopt
	>> _alphain >> _tolerance >> _maxtry
	>> _npPower >> iunit(_optInputScale,GeV) >> _quarkBranching
	>> iunit(_quarkMasses,GeV);
	}

	void ShowerAlphaQCDNP::Init() {

	static ClassDocumentation<ShowerAlphaQCDNP> documentation
	("This (concrete) class describes the QCD alpha running.");

	// default such that as(qmin) = 1 in the current parametrization.
	- static Parameter<ShowerAlphaQCDNP,Energy> intQmin
	+ static Parameter<ShowerAlphaQCDNP,Energy> interfaceQmin
	("Qmin", "Q < Qmin is treated with NP parametrization as of (unit [GeV])",
	- &ShowerAlphaQCDNP::_qmin, GeV, 1GeV, 1GeV,
	+ &ShowerAlphaQCDNP::_qmin, GeV, 1GeV, 0GeV,
	100.0*GeV,false,false,false);

	static Parameter<ShowerAlphaQCDNP,unsigned int> interfaceNumberOfLoops
	("NumberOfLoops",
	"The number of loops to use in the alpha_S calculation",
	&ShowerAlphaQCDNP::_nloop, 3, 1, 3,
	false, false, Interface::limited);

	static Parameter<ShowerAlphaQCDNP,double> interfaceAlphaMZ
	("AlphaIn",
	"The input value of the strong coupling at the chosen InputScale (default: MZ)",
	&ShowerAlphaQCDNP::_alphain, 0.118, 0.1, 0.2,
	false, false, Interface::limited);

	static Parameter<ShowerAlphaQCDNP,double> interfaceTolerance
	("Tolerance",
	"The tolerance for discontinuities in alphaS at thresholds.",
	&ShowerAlphaQCDNP::_tolerance, 1e-10, 1e-20, 1e-4,
	false, false, Interface::limited);

	static Parameter<ShowerAlphaQCDNP,unsigned int> interfaceMaximumIterations
	("MaximumIterations",
	"The maximum number of iterations for the Newton-Raphson method to converge.",
	&ShowerAlphaQCDNP::_maxtry, 100, 10, 1000,
	false, false, Interface::limited);

	static Switch<ShowerAlphaQCDNP,bool> interfaceThresholdOption
	("ThresholdOption",
	"Whether to use the consistuent or normal masses for the thresholds",
	&ShowerAlphaQCDNP::_thresopt, true, false, false);
	static SwitchOption interfaceThresholdOptionCurrent
	(interfaceThresholdOption,
	"Current",
	"Use the current masses",
	true);
	static SwitchOption interfaceThresholdOptionConstituent
	(interfaceThresholdOption,
	"Constituent",
	"Use the constitent masses.",
	false);

	static Command<ShowerAlphaQCDNP> interfaceValue
	("Value",
	"",
	&ShowerAlphaQCDNP::value, false);

	static Command<ShowerAlphaQCDNP> interfacecheck
	("check",
	"check",
	&ShowerAlphaQCDNP::check, false);

	static Command<ShowerAlphaQCDNP> interfacecheckNf
	("checkNf",
	"checkNf",
	&ShowerAlphaQCDNP::checkNf, false);

	static Command<ShowerAlphaQCDNP> interfacecheckscaleNP
	("checkscaleNP",
	"checkscaleNP",
	&ShowerAlphaQCDNP::checkscaleNP, false);

	static Parameter<ShowerAlphaQCDNP,Energy> interfaceInputScale
	("InputScale",
	"An optional input scale. MZ will be used if not set.",
	&ShowerAlphaQCDNP::_optInputScale, GeV, 91.1876_GeV, ZERO, 0*GeV,
	false, false, Interface::lowerlim);

	static ParVector<ShowerAlphaQCDNP,Energy> interfaceQuarkMasses
	("QuarkMasses",
	"The quark masses to be used instead of the masses set in the particle data.",
	&ShowerAlphaQCDNP::_quarkMasses, GeV, -1, 0.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);


	static Switch<ShowerAlphaQCDNP,bool> interfaceQuarkBranching
	("QuarkBranching",
	"True, if this coupling is used in a gluon to qqbar branching.",
	&ShowerAlphaQCDNP::_quarkBranching, false, false, false);
	static SwitchOption interfaceQuarkBranchingYes
	(interfaceQuarkBranching,
	"Yes",
	"Use in gluon to qqbar branching.",
	true);
	static SwitchOption interfaceQuarkBranchingNo
	(interfaceQuarkBranching,
	"No",
	"Use in gluon emission.",
	false);

	static Parameter<ShowerAlphaQCDNP,double> interfaceNPPower
	("NPPower",
	"The non-perturbative power law",
	- &ShowerAlphaQCDNP::_npPower, 2.0, 1.0, 0,
	- false, false, Interface::lowerlim);
	+ &ShowerAlphaQCDNP::_npPower, 2.0, 1.0, 10.,
	+ false, false, false);
	+

	}

	void ShowerAlphaQCDNP::doinit() {
	ShowerAlpha::doinit();
	// calculate the value of 5-flavour lambda
	// evaluate the initial
	// value of Lambda from alphas if needed using Newton-Raphson
	_lambda[2]=computeLambda(_optInputScale,_alphain,5);

	// compute the threshold matching
	// top threshold
	for(int ix=1;ix<4;++ix) {
	if ( _quarkMasses.empty() ) {
	if(_thresopt)
	_thresholds[ix]=getParticleData(ix+3)->mass();
	else
	_thresholds[ix]=getParticleData(ix+3)->constituentMass();
	} else {
	// starting at zero rather than one, cf the other alphas's
	_thresholds[ix] = _quarkMasses[ix+2];
	}
	}
	// compute 6 flavour lambda by matching at top mass using Newton Raphson
	_lambda[3]=computeLambda(_thresholds[3],alphaS(_thresholds[3],_lambda[2],5,_nloop),6);
	// bottom threshold
	// compute 4 flavour lambda by matching at bottom mass using Newton Raphson
	_lambda[1]=computeLambda(_thresholds[2],alphaS(_thresholds[2],_lambda[2],5,_nloop),4);
	// charm threshold
	// compute 3 flavour lambda by matching at charm mass using Newton Raphson
	_lambda[0]=computeLambda(_thresholds[1],alphaS(_thresholds[1],_lambda[1],4,_nloop),3);

	//Energy q = scaleNPMin();
	//auto nflam = getLamNf(q);
	//_asMaxNP = alphaS(q, nflam.second, nflam.first, _nloop);

	Energy q = scaleThatmMinimizesScaleNP() ;
	_asMaxNP = value(qq)1.05;

	//_absoluteCutoff=_qminpow(100.log(10.),-1./_npPower);
	- _absoluteCutoff=2.0*GeV;
	+ _absoluteCutoff=_qminpow(log(10e100GeV/_qmin),-1./_npPower);

	_nfMaxNP = nfNP(sqr(_absoluteCutoff));

	}

	string ShowerAlphaQCDNP::check(string args) {

	doinit();

	istringstream argin(args);

	double Q_low, Q_high;
	long n_steps;

	argin >> Q_low >> Q_high >> n_steps;

	string fname;
	argin >> fname;

	Repository::clog() << "checking alpha_s in range [" << Q_low << "," << Q_high << "] GeV in "
	<< n_steps << " steps.\nResults are written to " << fname << "\n";

	double step_width = (Q_high-Q_low)/n_steps;

	ofstream out (fname.c_str());

	for (long k = 0; k <= n_steps; ++k) {

	Energy Q = Q_lowGeV + kstep_width*GeV;

	//out << (Q/GeV) << " " << value(Q*Q) << " " << _asMaxNP << "\n";
	out << (Q/GeV) << " " << value(Q*Q) << "\n";


	}

	return "alpha_s check finished";

	}

	string ShowerAlphaQCDNP::checkNf(string args) {

	doinit();

	istringstream argin(args);

	double Q_low, Q_high;
	long n_steps;

	argin >> Q_low >> Q_high >> n_steps;

	string fname;
	argin >> fname;

	Repository::clog() << "checking nfNP in range [" << Q_low << "," << Q_high << "] GeV in "
	<< n_steps << " steps.\nResults are written to " << fname << "\n";

	double step_width = (Q_high-Q_low)/n_steps;

	ofstream out (fname.c_str());

	for (long k = 0; k <= n_steps; ++k) {

	Energy Q = Q_lowGeV + kstep_width*GeV;

	//out << (Q/GeV) << " " << nfNP(Q*Q) << " " << _nfMaxNP << "\n";
	out << (Q/GeV) << " " << nfNP(Q*Q) << "\n";

	}

	return "nfNP check finished";

	}



	string ShowerAlphaQCDNP::checkscaleNP(string args) {

	doinit();

	istringstream argin(args);

	double Q_low, Q_high;
	long n_steps;

	argin >> Q_low >> Q_high >> n_steps;

	string fname;
	argin >> fname;

	Repository::clog() << "checking scaleNP (in GeV) in range [" << Q_low << "," << Q_high << "] GeV in "
	<< n_steps << " steps.\nResults are written to " << fname << "\n";

	double step_width = (Q_high-Q_low)/n_steps;

	ofstream out (fname.c_str());

	for (long k = 0; k <= n_steps; ++k) {

	Energy Q = Q_lowGeV + kstep_width*GeV;

	out << (Q/GeV) << " " << scaleNP(scaleFactor()Q)/(1.GeV) << "\n";

	}

	return "scaleNP check finished";

	}


	double ShowerAlphaQCDNP::value(const Energy2 scale) const {

	Energy q = scaleFactor()*sqrt(scale);
	+
	+ if (q<_absoluteCutoff){
	+ auto nflam = getLamNf(_absoluteCutoff);
	+ double as = alphaS(scaleNP(_absoluteCutoff), nflam.second, nflam.first, _nloop);
	+ return as*pow(q/_absoluteCutoff,_npPower);
	+ } else {
	+ auto nflam = getLamNf(q);
	+ return alphaS(scaleNP(q), nflam.second, nflam.first, _nloop);
	+ }

	- if ( q > _absoluteCutoff ) {
	- auto nflam = getLamNf(q);
	- return alphaS(scaleNP(q), nflam.second, nflam.first, _nloop);
	- }
	-
	-
	- return 0.;

	}


	string ShowerAlphaQCDNP::value(string scale) {
	istringstream readscale(scale);
	double inScale; readscale >> inScale;
	Energy theScale = inScale * GeV;
	initialize();
	ostringstream showvalue ("");
	showvalue << "alpha_s (" << theScale/GeV << " GeV) = "
	<< value (sqr(theScale));
	return showvalue.str();
	}

	pair<short, Energy> ShowerAlphaQCDNP::getLamNf(Energy q) const {
	short nf = 6;
	// get lambda and nf according to the thresholds
	if (q < _thresholds[1]) nf = 3;
	else if (q < _thresholds[2]) nf = 4;
	else if (q < _thresholds[3]) nf = 5;
	return pair<short,Energy>(nf, _lambda[nf-3]);
	}

	Energy ShowerAlphaQCDNP::computeLambda(Energy match,
	double alpha,
	unsigned int nflav) const {

	Energy lamtest=200.0*MeV;
	double xtest;
	unsigned int ntry=0;
	do {
	++ntry;
	xtest = log(sqr(match/lamtest));
	xtest += (alpha-alphaS(match,lamtest,nflav,_nloop))/derivativeAlphaS(match,lamtest,nflav,_nloop);
	Energy newLambda = match/exp(0.5 *xtest);
	lamtest = newLambda<match ? newLambda : 0.5*(lamtest+match);
	}
	while(abs(alpha-alphaS(match,lamtest,nflav,_nloop)) > _tolerance && ntry < _maxtry);
	return lamtest;

	}

	// --------------------------------------------------------------------------------
	// main modification for nonperturbative running in terms of changed argument and
	// nonperturbative number of flavours
	// --------------------------------------------------------------------------------

	double ShowerAlphaQCDNP::nfNP(const Energy2 q2) const {
	Energy q=sqrt(q2);

	if (q<_absoluteCutoff){
	- return 0.;
	+ return nfNP(_absoluteCutoff_absoluteCutoff)pow(_absoluteCutoff/q,_npPower);
	}
	else{

	//get p-nf
	double nf=getLamNf(q).first;

	//beta function coefficients; b=sum_i (alphaS/(4Pi))^i * sum_j b_ij *nf^j
	using Constants::pi;
	double b00 = 11.;
	double b01 = - 2./3.;
	double b10 = 51.;
	double b11 = - 19./3.;
	double c = (b00+b10value(q2)/(4.pi))/(b01+b11value(q2)/(4.pi));

	//absolute number of np-flavors
	double nfNPabs=(c+nf)(qscaleNPDerivative(q)/scaleNP(q))-c;

	//return ratio of np-flavors over p-flavors
	return nfNPabs/nf;
	}
	}

	Energy ShowerAlphaQCDNP::scaleNP(const Energy q) const {
	return q(1.+(_qmin/q)(exp(pow(_qmin/q,_npPower))-1.));
	}

	double ShowerAlphaQCDNP::scaleNPDerivative(const Energy q) const {
	return 1.-exp(pow(_qmin/q,_npPower))_npPowerpow(_qmin/q,1.+_npPower);
	}

	Energy ShowerAlphaQCDNP::scaleThatmMinimizesScaleNP() const {
	double x=(1./_npPower)*gsl_sf_lambert_W0(pow(_npPower,1./(1.+_npPower))/(1.+_npPower));
	return _qminexp(x)pow(_npPower,1./(1.+_npPower));

	}

	Energy ShowerAlphaQCDNP::scaleNPMin() const {
	double x=((1.+_npPower)/_npPower)*gsl_sf_lambert_W0(pow(_npPower,1./(1.+_npPower))/(1.+_npPower));
	return _qmin*(-1.+exp(x)+pow(x,-1./_npPower));
	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 5:53 AM (17 h, 31 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4972828
Default Alt Text: (58 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions