diff --git a/Decay/Partonic/DarkQuarkoniumDecayer.cc b/Decay/Partonic/DarkQuarkoniumDecayer.cc
new file mode 100644
--- /dev/null
+++ b/Decay/Partonic/DarkQuarkoniumDecayer.cc
@@ -0,0 +1,152 @@
+// -*- C++ -*-
+//
+// DarkQuarkoniumDecayer.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 HeavyDecayer class.
+//
+
+#include "DarkQuarkoniumDecayer.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include <ThePEG/PDT/EnumParticles.h>
+#include <ThePEG/PDT/DecayMode.h>
+#include <ThePEG/Interface/ClassDocumentation.h>
+#include <ThePEG/Interface/Switch.h>
+#include "Herwig/Utilities/Kinematics.h"
+#include <ThePEG/Persistency/PersistentOStream.h>
+#include <ThePEG/Persistency/PersistentIStream.h>
+#include <ThePEG/Repository/UseRandom.h>
+#include <cassert>
+
+using namespace Herwig;
+
+DarkQuarkoniumDecayer::DarkQuarkoniumDecayer() : MECode(0) {}
+
+IBPtr DarkQuarkoniumDecayer::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr DarkQuarkoniumDecayer::fullclone() const {
+ return new_ptr(*this);
+}
+
+void DarkQuarkoniumDecayer::Init() {
+
+ static ClassDocumentation<DarkQuarkoniumDecayer> documentation
+ ("The DarkQuarkoniumDecayer performs partonic decays of quarkonium"
+ " resonances");
+
+ static Switch<DarkQuarkoniumDecayer,int> interfaceMECode
+ ("MECode",
+ "The code for the ME type to use in the decay",
+ &DarkQuarkoniumDecayer::MECode, 0, false, false);
+ static SwitchOption interfaceMECodePhaseSpace
+ (interfaceMECode,
+ "PhaseSpace",
+ "Use a phase-space distribution",
+ 0);
+ static SwitchOption interfaceMECodeOrePowell
+ (interfaceMECode,
+ "OrePowell",
+ "Use the Ore-Powell matrix element",
+ 130);
+
+}
+
+// The following static variable is needed for the type
+// description system in ThePEG.
+DescribeClass<DarkQuarkoniumDecayer,PartonicDecayerBase>
+describeHerwigDarkQuarkoniumDecayer("Herwig::DarkQuarkoniumDecayer", "HwShower.so HwPartonicDecay.so");
+
+bool DarkQuarkoniumDecayer::accept(tcPDPtr, const tPDVector & children) const {
+ return (children.size() == 3 || children.size() == 2);
+}
+
+ParticleVector DarkQuarkoniumDecayer::decay(const Particle & p,
+ const tPDVector & children) const {
+ ParticleVector partons;
+ for(unsigned int ix=0;ix<children.size();++ix) {
+ partons.push_back(children[ix]->produceParticle());
+ }
+ assert(partons.size()==2 || partons.size()==3);
+ Lorentz5Momentum products[3];
+ Energy gluMass = getParticleData(ParticleID::g)->constituentMass();
+ for(unsigned int i = 0; i<partons.size(); i++) {
+ if(partons[i]->id() == ParticleID::g) products[i].setMass(gluMass);
+ else products[i].setMass(partons[i]->mass());
+ }
+ if(partons.size() == 3) {
+ // 2 gluon or quark + dark pion decay
+ // Ore & Powell orthopositronium matrix element
+ if(MECode == 130) {
+ double x1, x2, x3, test;
+ do {
+ // if decay fails return empty vector.
+ if (! Kinematics::threeBodyDecay(p.momentum(), products[0],
+ products[1], products[2]))
+ return ParticleVector();
+ x1 = 2.*(p.momentum()*products[0])/sqr(p.mass());
+ x2 = 2.*(products[1]*products[2])/sqr(p.mass());
+ x3 = 2. - x1 - x2;
+ test = sqr(x1*(1.-x1)) + sqr(x2*(1.-x2)) + sqr(x3*(1.-x3));
+ test /= sqr(x1*x2*x3);
+ }
+ while(test < 2.*UseRandom::rnd());
+ }
+ else {
+ if (! Kinematics::threeBodyDecay(p.momentum(), products[0],
+ products[1],products[2]))
+ return ParticleVector();
+ }
+ // test the momenta
+ for(unsigned int i = 0; i<partons.size(); i++)
+ partons[i]->set5Momentum(products[i]);
+ // Now set colour connections
+ partons[0]->colourNeighbour(partons[1]);
+ if(abs(partons[0]->id()) == ParticleID::g)
+ partons[0]->colourNeighbour(partons[1]);
+ }
+ // two decay children
+ // 2 gluon or q-qbar decay
+ else {
+ double Theta, Phi;
+ Kinematics::generateAngles(Theta, Phi);
+ Energy p1 = partons[0]->mass();
+ Energy p2 = partons[1]->mass();
+ if(p1 == ZERO) p1 = gluMass;
+ if(p2 == ZERO) p2 = gluMass;
+ if (! Kinematics::twoBodyDecay(p.momentum(), p1, p2, Theta, Phi,
+ products[0], products[1]))
+ return ParticleVector();
+ for(unsigned int i = 0; i<partons.size(); i++)
+ partons[i]->set5Momentum(products[i]);
+ int first(0), second(1);
+ if(partons[0]->id() < 0) swap(first,second);
+ partons[first]->antiColourNeighbour(partons[second]);
+ if(abs(partons[first]->id()) == ParticleID::g)
+ partons[first]->colourNeighbour(partons[second]);
+ }
+ return partons;
+}
+
+void DarkQuarkoniumDecayer::persistentOutput(PersistentOStream &os) const {
+ os << MECode;
+}
+
+void DarkQuarkoniumDecayer::persistentInput(PersistentIStream &is, int) {
+ is >> MECode;
+}
+
+void DarkQuarkoniumDecayer::dataBaseOutput(ofstream & output, bool header) const {
+ if(header) output << "update decayers set parameters=\"";
+ // parameters for the PartonicDecayerBase base class
+ PartonicDecayerBase::dataBaseOutput(output,false);
+ output << "newdef " << name() << ":MECode " << MECode << " \n";
+ if(header) output << "\n\" where BINARY ThePEGName=\""
+ << fullName() << "\";" << endl;
+}
diff --git a/Decay/Partonic/DarkQuarkoniumDecayer.h b/Decay/Partonic/DarkQuarkoniumDecayer.h
new file mode 100644
--- /dev/null
+++ b/Decay/Partonic/DarkQuarkoniumDecayer.h
@@ -0,0 +1,133 @@
+// -*- C++ -*-
+//
+// DarkQuarkoniumDecayer.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_DarkQuarkoniumDecayer_H
+#define HERWIG_DarkQuarkoniumDecayer_H
+//
+// This is the declaration of the DarkQuarkoniumDecayer class.
+//
+
+#include "PartonicDecayerBase.h"
+
+namespace Herwig {
+using namespace ThePEG;
+
+/** \ingroup Decay
+ *
+ * The DarkQuarkoniumDecayer class is designed for the partonic decay of bottom and charmonium
+ * resonances. In general it is used for decays of the type:
+ * - \f$q,\bar{q}\f$ decay to a quark-antiquark pair generally using phase space,
+ * \e i.e. MECode=0.
+ *
+ * - \f$g,g\f$ decay to two gluons normally using phase space,
+ * \e i.e. MECode=0.
+ *
+ * - \f$g,g,g\f$ decay to three gluons, this will normally use the Ore-Powell
+ * matrix element, \e i.e. MECode=130.
+ *
+ * - \f$g,g,\gamma\f$ decay to two gluons and a photon, this will normally use
+ * the Ore-Powell matrix element, \e i.e. MECode=130.
+ *
+ *
+ * This class supports two values of the MECode variable which can be set using
+ * the interface
+ *
+ * - MECode=0 flat-phase space
+ * - MECode=130 The Ore-Powell onium matrix element.
+ *
+ * This is designed to be the same as the FORTRAN HERWIG routine.
+ *
+ * @see HeavyDecayer
+ * @see Hw64Decayer
+ * @see Decayer
+ *
+ */
+class DarkQuarkoniumDecayer: public PartonicDecayerBase {
+
+public:
+
+ /**
+ * Standard ctors and dtor
+ */
+ DarkQuarkoniumDecayer();
+
+ /**
+ * Check if this decayer can perfom the decay for a particular mode
+ * @param parent The decaying particle
+ * @param children The decay products
+ * @return true If this decayer can handle the given mode, otherwise false.
+ */
+ virtual bool accept(tcPDPtr parent, const tPDVector & children) const;
+
+ /**
+ * Perform the decay of the particle to the specified decay products
+ * @param parent The decaying particle
+ * @param children The decay products
+ * @return a ParticleVector containing the decay products.
+ */
+ virtual ParticleVector decay(const Particle & parent,
+ const tPDVector & children) const;
+
+
+ /**
+ * Output the setup information for the particle database
+ * @param os The stream to output the information to
+ * @param header Whether or not to output the information for MySQL
+ */
+ virtual void dataBaseOutput(ofstream & os,bool header) const;
+
+public:
+
+ /**
+ * Standard Init function used to initialize the interface.
+ */
+ static void Init();
+
+ /**
+ * Standard Persistent stream methods
+ */
+ void persistentOutput(PersistentOStream &) const;
+
+ /**
+ * Standard Persistent stream methods
+ */
+ void persistentInput(PersistentIStream &, int);
+
+ /**
+ * Standard clone methods
+ */
+protected:
+
+ /**
+ * Standard clone methods
+ */
+ virtual IBPtr clone() const;
+
+ /**
+ * Standard clone methods
+ */
+ virtual IBPtr fullclone() const;
+
+private:
+
+ /**
+ * Private and non-existent assignment operator.
+ */
+ const DarkQuarkoniumDecayer & operator=(const DarkQuarkoniumDecayer &) = delete;
+
+private:
+
+ /**
+ * The code for the type of matrix element being used.
+ */
+ int MECode;
+};
+
+}
+
+#endif /* HERWIG_DarkQuarkoniumDecayer_H */
diff --git a/Decay/Partonic/Makefile.am b/Decay/Partonic/Makefile.am
--- a/Decay/Partonic/Makefile.am
+++ b/Decay/Partonic/Makefile.am
@@ -1,24 +1,26 @@
BUILT_SOURCES = Partonic__all.cc
CLEANFILES = Partonic__all.cc
Partonic__all.cc : $(DIR_H_FILES) $(DIR_CC_FILES) Makefile
@echo "Concatenating .cc files into $@"
@$(top_srcdir)/cat_with_cpplines $(DIR_CC_FILES) > $@
EXTRA_DIST = $(ALL_H_FILES) $(ALL_CC_FILES)
DIR_H_FILES = $(addprefix $(srcdir)/,$(ALL_H_FILES))
ALL_H_FILES = \
QuarkoniumDecayer.h \
+DarkQuarkoniumDecayer.h \
HeavyDecayer.h \
WeakPartonicDecayer.h \
BtoSGammaDecayer.h \
PartonicDecayerBase.h
DIR_CC_FILES = $(addprefix $(srcdir)/,$(ALL_CC_FILES))
ALL_CC_FILES = \
+DarkQuarkoniumDecayer.cc \
QuarkoniumDecayer.cc \
HeavyDecayer.cc \
WeakPartonicDecayer.cc\
BtoSGammaDecayer.cc\
PartonicDecayerBase.cc
diff --git a/Decay/ScalarMeson/EtaPiGammaGammaDecayer.cc b/Decay/ScalarMeson/EtaPiGammaGammaDecayer.cc
--- a/Decay/ScalarMeson/EtaPiGammaGammaDecayer.cc
+++ b/Decay/ScalarMeson/EtaPiGammaGammaDecayer.cc
@@ -1,374 +1,374 @@
// -*- C++ -*-
//
// EtaPiGammaGammaDecayer.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 EtaPiGammaGammaDecayer class.
//
#include "EtaPiGammaGammaDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include "ThePEG/Helicity/HelicityFunctions.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
void EtaPiGammaGammaDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
_etamax = mode(0)->maxWeight();
_etapmax = mode(1)->maxWeight();
}
}
EtaPiGammaGammaDecayer::EtaPiGammaGammaDecayer()
: _grhoomega(12.924/GeV), _fpi(130.7*MeV),_rhomass(771.1*MeV),
_rhowidth(149.2*MeV),_grho(_rhomass/_fpi),_mpi(ZERO),_rhoconst(0.),
_localparameters(true),_ratiofpif8(1./1.3),_ratiofpif0(1./1.04),
_theta(-Constants::pi/9.),_etamax(2.36858),_etapmax(0.006),
_dconst(2), _econst(2) {
// intermediates
generateIntermediates(false);
}
void EtaPiGammaGammaDecayer::doinit() {
DecayIntegrator::doinit();
// set rho parameters if needed
tPDPtr rho(getParticleData(ParticleID::rho0));
if(!_localparameters) {
_rhomass = rho->mass();
_rhowidth = rho->width();
}
// constant for the running rho width
_mpi=getParticleData(ParticleID::pi0)->mass();
Energy pcm =Kinematics::pstarTwoBodyDecay(_rhomass,_mpi,_mpi);
_rhoconst=_rhomass*_rhomass*_rhowidth/(pcm*pcm*pcm);
// set the prefactors
double conv(sqrt(4.*Constants::pi*SM().alphaEM()));
conv *=_fpi*_fpi*_grho/_rhomass/_rhomass;
InvEnergy2 pre(2.*sqrt(3.)/9.*sqr(_grhoomega*conv));
double fact[2];
// constants for eta
fact[0] = _ratiofpif8*cos(_theta)-sqrt(2.)*_ratiofpif0*sin(_theta);
// constants for eta'
fact[1] = _ratiofpif8*sin(_theta)+sqrt(2.)*_ratiofpif0*cos(_theta);
for(unsigned int ix=0;ix<2;++ix) {
_dconst[ix]=fact[ix]*pre;
_econst[ix]=fact[ix]*pre;
}
- // set up the phsae space for the decays
+ // set up the phase space for the decays
tPDPtr eta[2]={getParticleData(ParticleID::eta),
getParticleData(ParticleID::etaprime)};
tPDVector out = {getParticleData(ParticleID::pi0),
getParticleData(ParticleID::gamma),
getParticleData(ParticleID::gamma)};
for(unsigned int ix=0;ix<2;++ix) {
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(eta[ix],out,
(ix==0 ? _etamax : _etapmax)));
PhaseSpaceChannel c1((PhaseSpaceChannel(mode),0,rho,0,2,1,1,1,3));
c1.weight(0.5);
mode->addChannel(c1);
PhaseSpaceChannel c2((PhaseSpaceChannel(mode),0,rho,0,3,1,1,1,2));
c2.weight(0.5);
mode->addChannel(c2);
addMode(mode);
}
}
int EtaPiGammaGammaDecayer::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
cc=false;
int id;
if(children.size()!=3) return -1;
tPDVector::const_iterator pit = children.begin();
unsigned int npi0(0),ngamma(0);
for( ;pit!=children.end();++pit) {
id=(**pit).id();
if(id==ParticleID::pi0) ++npi0;
else if(id==ParticleID::gamma) ++ngamma;
}
if(!(npi0==1&&ngamma==2)) return -1;
// number of the mode
switch (parent->id()) {
case ParticleID::eta : return 0;
case ParticleID::etaprime: return 1;
default: return -1;
}
}
void EtaPiGammaGammaDecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(_grhoomega,1/GeV)<< ounit(_fpi,GeV)<< _grho
<< ounit(_rhomass,GeV)<< ounit(_rhowidth,GeV)<< _localparameters
<< _ratiofpif8 << _ratiofpif0 << _theta << _etamax << _etapmax
<< _rhoconst << ounit(_mpi,GeV) << ounit(_dconst,1/GeV2)
<< ounit(_econst,1/GeV2);
}
void EtaPiGammaGammaDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(_grhoomega,1/GeV) >> iunit(_fpi,GeV)>> _grho
>> iunit(_rhomass,GeV)>> iunit(_rhowidth,GeV)>> _localparameters
>> _ratiofpif8 >> _ratiofpif0 >> _theta >> _etamax >> _etapmax
>> _rhoconst >> iunit(_mpi,GeV) >> iunit(_dconst,1/GeV2)
>> iunit(_econst,1/GeV2);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<EtaPiGammaGammaDecayer,DecayIntegrator>
describeHerwigEtaPiGammaGammaDecayer("Herwig::EtaPiGammaGammaDecayer", "HwSMDecay.so");
void EtaPiGammaGammaDecayer::Init() {
static ClassDocumentation<EtaPiGammaGammaDecayer> documentation
("The EtaPiGammaGammaDecayer class implements a VMD model for the"
" decay of the eta or etaprime to a pion and two photons.",
"The decays of $\\eta,\\eta'\\to\\pi^0\\gamma\\gamma$ were simulated using"
" the matrix elements of \\cite{Holstein:2001bt}",
"\\bibitem{Holstein:2001bt} B.~R.~Holstein,\n"
" Phys.\\ Scripta {\\bf T99} (2002) 55 [arXiv:hep-ph/0112150].\n"
"%%CITATION = PHSTB,T99,55;%%\n");
static Parameter<EtaPiGammaGammaDecayer,InvEnergy> interfacegrhoomega
("grhoomega",
"The couping of the rho, omega and a pion",
&EtaPiGammaGammaDecayer::_grhoomega, 1./GeV, 12.924/GeV, ZERO, 100./GeV,
false, false, true);
static Parameter<EtaPiGammaGammaDecayer,Energy> interfaceFpi
("Fpi",
"The pion decay constant",
&EtaPiGammaGammaDecayer::_fpi, MeV, 130.7*MeV, ZERO, 200.0*MeV,
false, false, true);
static Parameter<EtaPiGammaGammaDecayer,double> interfacegrho
("grho",
"Rho decay constant",
&EtaPiGammaGammaDecayer::_grho, 5.9, 0.0, 10.0,
false, false, true);
static Parameter<EtaPiGammaGammaDecayer,Energy> interfaceRhoMass
("RhoMass",
"The mass of the rho meson",
&EtaPiGammaGammaDecayer::_rhomass, MeV, 771.1*MeV, 500.0*MeV, 1000.0*MeV,
false, false, true);
static Parameter<EtaPiGammaGammaDecayer,Energy> interfaceRhoWidth
("RhoWidth",
"The width of the rho meson",
&EtaPiGammaGammaDecayer::_rhowidth, MeV, 149.2*MeV, 100.0*MeV, 200.0*MeV,
false, false, true);
static Parameter<EtaPiGammaGammaDecayer,double> interfaceRatioFpiF8
("RatioFpiF8",
"The ratio of the decay constant Fpi to F8",
&EtaPiGammaGammaDecayer::_ratiofpif8, 1./1.3, 0.0, 10.0,
false, false, true);
static Parameter<EtaPiGammaGammaDecayer,double> interfaceRatioFpiF0
("RatioFpiF0",
"The ratio of the decay constant Fpi to F0",
&EtaPiGammaGammaDecayer::_ratiofpif0, 1./1.04, 0.0, 10.0,
false, false, true);
static Parameter<EtaPiGammaGammaDecayer,double> interfaceTheta
("Theta",
"The eta etaprime mixing angle",
&EtaPiGammaGammaDecayer::_theta, -Constants::pi/9., -Constants::pi, Constants::pi,
false, false, true);
static Parameter<EtaPiGammaGammaDecayer,double> interfaceEtaMax
("EtaMax",
"THe maximum weight for the eta decay",
&EtaPiGammaGammaDecayer::_etamax, 1.35, -1.0e12, 1.0e12,
false, false, false);
static Parameter<EtaPiGammaGammaDecayer,double> interfaceEtaPrimeMax
("EtaPrimeMax",
"THe maximum weight for the eta prime decay",
&EtaPiGammaGammaDecayer::_etapmax, 0.006, -1.0e12, 1.0e12,
false, false, false);
static Switch<EtaPiGammaGammaDecayer,bool> interfaceLocalParameters
("LocalParameters",
"Use local values of the parameters",
&EtaPiGammaGammaDecayer::_localparameters, true, false, false);
static SwitchOption interfaceLocalParametersLocal
(interfaceLocalParameters,
"Local",
"Use local values",
true);
static SwitchOption interfaceLocalParametersParticleData
(interfaceLocalParameters,
"ParticleData",
"Use values from the particle data objects",
false);
}
void EtaPiGammaGammaDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
// set up the spin information for the decay products
ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
incoming,true);
ScalarWaveFunction::constructSpinInfo(decay[0],outgoing,true);
for(unsigned int ix=0;ix<2;++ix)
VectorWaveFunction::constructSpinInfo(_vectors[ix],decay[ix+1],
outgoing,true,true);
}
double EtaPiGammaGammaDecayer::me2(const int,const Particle & part,
const tPDVector &,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin1,PDT::Spin1)));
useMe();
if(meopt==Initialize) {
ScalarWaveFunction::
calculateWaveFunctions(_rho,const_ptr_cast<tPPtr>(&part),incoming);
}
for(unsigned int ix=0;ix<2;++ix) {
_vectors[ix].resize(3);
for(unsigned int ihel=0;ihel<3;ihel+=2) {
_vectors[ix][ihel] = HelicityFunctions::polarizationVector(-momenta[ix+1],ihel,Helicity::outgoing);
}
}
// dot products we need
Energy2 q1dotq2(momenta[1]*momenta[2]),
pdotq1(part.momentum()*momenta[1]),
pdotq2(part.momentum()*momenta[2]);
complex<Energy> e1dotq2[3],e1dotp[3],e2dotq1[3],e2dotp[3];
for(unsigned int ix=0;ix<3;++ix) {
if(ix==1) {
e1dotq2[ix]=ZERO;
e1dotp[ix] =ZERO;
e2dotq1[ix]=ZERO;
e2dotp[ix] =ZERO;
}
else {
e1dotq2[ix] =_vectors[0][ix]*momenta[2];
e1dotp[ix] =_vectors[0][ix]*part.momentum();
e2dotq1[ix] =_vectors[1][ix]*momenta[1];
e2dotp[ix] =_vectors[1][ix]*part.momentum();
}
}
// the momentum dependent pieces of the matrix element
Complex ii(0.,1.);
Energy2 mpi2(sqr(momenta[0].mass())),meta2(sqr(part.mass())),
mrho2(sqr(_rhomass)),
t(mpi2+2.*((momenta[0])*(momenta[1]))),
u(mpi2+2.*((momenta[0])*(momenta[2])));
Energy q(sqrt(t)),pcm(Kinematics::pstarTwoBodyDecay(q,_mpi,_mpi));
complex<Energy2> tgamma(ii*pcm*pcm*pcm*_rhoconst/q);
q=sqrt(u);pcm = Kinematics::pstarTwoBodyDecay(q,_mpi,_mpi);
complex<Energy2> ugamma(ii*pcm*pcm*pcm*_rhoconst/q);
complex<InvEnergy2> prop1(1./(mrho2-t-tgamma)),prop2(1./(mrho2-u-ugamma));
complex<InvEnergy2> Dfact(_dconst[imode()]*(prop1*(pdotq2-meta2)
+prop2*(pdotq1-meta2)));
complex<InvEnergy4> Efact(_econst[imode()]*(prop1+prop2));
Complex e1dote2;
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
if(ix==1||iy==1) (*ME())(0,0,ix,iy)=0.;
else {
e1dote2=_vectors[0][ix].dot(_vectors[1][iy]);
(*ME())(0,0,ix,iy) =
Complex(Dfact*complex<Energy2>(e1dote2*q1dotq2-
e1dotq2[ix]*e2dotq1[iy])
-Efact*complex<Energy4>(-e1dote2*pdotq1*pdotq2
-e1dotp[ix]*e2dotp[iy]*q1dotq2
+e1dotq2[ix]*e2dotp[iy]*pdotq1
+e1dotp[ix]*e2dotq1[iy]*pdotq2));
}
}
}
double me(ME()->contract(_rho).real());
// test of the me
// Energy M(part.mass());
// Energy2 M2(M*M);
// Energy2 s1(2.*(momenta[1]*momenta[2]));
// Energy2 s2(M2-2.*(part.momentum()*momenta[1]));
// Energy2 s3(M2-2.*(part.momentum()*momenta[2]));
// cout << "testing the matrix element " << (
// 2*(2*(Dfact*conj(Dfact)).real() + 2*(Dfact*conj(Efact)).real()*M2
// + (Efact*conj(Efact)).real()*M2*M2)*
// s1*s1 - 2*(Efact*conj(Efact)).real()*M2*s1*(M2 - s2)*
// (M2 - s3) +(Efact*conj(Efact)).real()*(M2 - s2)*(M2 - s2)*
// (M2-s3)*(M2-s3))/8. - me << endl;
return me;
}
double EtaPiGammaGammaDecayer::
threeBodyMatrixElement(const int imodeb, const Energy2 q2,const Energy2 s3,
const Energy2 s2,const Energy2 s1,const Energy ,
const Energy ,const Energy ) const {
// compute the prefactors
Energy2 mrho2 = sqr(_rhomass);
Energy q = sqrt(s3);
Energy pcm = Kinematics::pstarTwoBodyDecay(q,_mpi,_mpi);
Complex ii(0.,1.);
complex<Energy2> tgamma(ii*pcm*pcm*pcm*_rhoconst/q);
q = sqrt(s2);
pcm = Kinematics::pstarTwoBodyDecay(q,_mpi,_mpi);
complex<Energy2> ugamma(ii*pcm*pcm*pcm*_rhoconst/q);
complex<InvEnergy2> prop1(1./(mrho2-s3-tgamma)), prop2(1./(mrho2-s2-ugamma));
Energy2 pdotq2(0.5*(q2-s3)), pdotq1(0.5*(q2-s2));
complex<InvEnergy2> Dfact(_dconst[imodeb]*(prop1*(pdotq2-q2)+prop2*(pdotq1-q2)));
complex<InvEnergy4> Efact(_econst[imodeb]*(prop1+prop2));
InvEnergy4 D2 = (Dfact*conj(Dfact)).real();
InvEnergy8 E2((Efact*conj(Efact)).real());
InvEnergy6 ED((Efact*conj(Dfact)).real());
return (2 * (2*D2 + 2*ED*q2 + E2*sqr(q2)) * sqr(s1)
- double(2*E2*q2*s1*(q2-s2)*(q2-s3))
+ double(E2*sqr(q2-s2)*sqr(q2-s3))
)/8.;
}
WidthCalculatorBasePtr
EtaPiGammaGammaDecayer::threeBodyMEIntegrator(const DecayMode & dm) const {
// workout which mode we are doing
int id(dm.parent()->id()),imode(1);
if(id==ParticleID::eta){imode=0;}
// construct the integrator
vector<double> inweights; inweights.push_back(0.5); inweights.push_back(0.5);
Energy mrho(getParticleData(ParticleID::rhoplus)->mass());
Energy wrho(getParticleData(ParticleID::rhoplus)->width());
vector<Energy> inmass; inmass.push_back(mrho); inmass.push_back(mrho);
vector<Energy> inwidth; inwidth.push_back(wrho); inwidth.push_back(wrho);
vector<int> intype; intype.push_back(1); intype.push_back(2);
vector<double> inpow(2,0.0);
return new_ptr(ThreeBodyAllOnCalculator<EtaPiGammaGammaDecayer>
(inweights,intype,inmass,inwidth,inpow,*this,
imode,_mpi,ZERO,ZERO));
}
void EtaPiGammaGammaDecayer::dataBaseOutput(ofstream & output,
bool header) const {
if(header) output << "update decayers set parameters=\"";
DecayIntegrator::dataBaseOutput(output,false);
output << "newdef " << name() << ":grhoomega " << _grhoomega*GeV << "\n";
output << "newdef " << name() << ":Fpi " << _fpi/MeV << "\n";
output << "newdef " << name() << ":grho " << _grho << "\n";
output << "newdef " << name() << ":RhoMass " << _rhomass/MeV << "\n";
output << "newdef " << name() << ":RhoWidth " << _rhowidth/MeV << "\n";
output << "newdef " << name() << ":RatioFpiF8 " << _ratiofpif8 << "\n";
output << "newdef " << name() << ":RatioFpiF0 " << _ratiofpif0 << "\n";
output << "newdef " << name() << ":Theta " << _theta << "\n";
output << "newdef " << name() << ":EtaMax " << _etamax << "\n";
output << "newdef " << name() << ":EtaPrimeMax " << _etapmax << "\n";
output << "newdef " << name() << ":LocalParameters " << _localparameters << "\n";
if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
}
diff --git a/Hadronization/DarkHadronSpectrum.cc b/Hadronization/DarkHadronSpectrum.cc
--- a/Hadronization/DarkHadronSpectrum.cc
+++ b/Hadronization/DarkHadronSpectrum.cc
@@ -1,373 +1,376 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DarkHadronSpectrum class.
//
#include "DarkHadronSpectrum.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ParMap.h"
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Repository/EventGenerator.h>
#include <ThePEG/Repository/CurrentGenerator.h>
#include <ThePEG/Repository/Repository.h>
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
namespace {
bool weightIsLess (pair<long,double> a, pair<long,double> b) {
return a.second < b.second;
}
}
DarkHadronSpectrum::DarkHadronSpectrum(unsigned int opt)
: HadronSpectrum(),
_topt(opt),_trial(0),
_nlightquarks(2),
_nheavyquarks(0) {}
DarkHadronSpectrum::~DarkHadronSpectrum() {}
void DarkHadronSpectrum::persistentOutput(PersistentOStream & os) const {
os << _sngWt << _decWt << _repwt << _pwtDIquark
<< _nlightquarks << _nheavyquarks
<< belowThreshold_;
}
void DarkHadronSpectrum::persistentInput(PersistentIStream & is, int) {
is >> _sngWt >> _decWt >> _repwt >> _pwtDIquark
>> _nlightquarks >> _nheavyquarks
>> belowThreshold_;
}
// *** 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).
DescribeAbstractClass<DarkHadronSpectrum,HadronSpectrum>
describeHerwigDarkHadronSpectrum("Herwig::DarkHadronSpectrum", "Herwig.so");
void DarkHadronSpectrum::Init() {
static ClassDocumentation<DarkHadronSpectrum> documentation
("There is no documentation for the DarkHadronSpectrum class");
static ParMap<DarkHadronSpectrum,double>
interfacePwt("Pwt","Weights for choosing the quarks",
&DarkHadronSpectrum::_pwt, 0, -1, 1.0, 0.0, 10.0,
false, false, false);
static Parameter<DarkHadronSpectrum,double>
interfacePwtDIquark("PwtDIquark","Weight for choosing a DIquark",
&DarkHadronSpectrum::_pwtDIquark, 0, 1.0, 0.0, 100.0,
false,false,false);
static Parameter<DarkHadronSpectrum,int> interfaceNLightQuarks
("NumLightQuarks",
"The number of light quarks (roughly mass degenerate)",
&DarkHadronSpectrum::_nlightquarks, 0, 2, 0, 9, false,false,false);
static Parameter<DarkHadronSpectrum,int> interfaceNHeavyQuarks
("NumHeavyQuarks",
"The number of heavy quarks (non-mass degenerate)",
&DarkHadronSpectrum::_nheavyquarks, 0, 0, 0, 9, false,false,false);
//
// mixing angles
//
// the ideal mixing angle
/*
// the meson weights
//
static ParVector<DarkHadronSpectrum,double> interface1S0Weights
("1S0Weights",
"The weights for the 1S0 multiplets start with n=1.",
&DarkHadronSpectrum::_weight1S0, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<DarkHadronSpectrum,double> interface3S1Weights
("3S1Weights",
"The weights for the 3S1 multiplets start with n=1.",
&DarkHadronSpectrum::_weight3S1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
*/
static Switch<DarkHadronSpectrum,unsigned int> interfaceTrial
("Trial",
"A Debugging option to only produce certain types of hadrons",
&DarkHadronSpectrum::_trial, 0, false, false);
static SwitchOption interfaceTrialAll
(interfaceTrial,
"All",
"Produce all the hadrons",
0);
static SwitchOption interfaceTrialPions
(interfaceTrial,
"Pions",
"Only produce pions",
1);
static SwitchOption interfaceTrialSpin2
(interfaceTrial,
"Spin2",
"Only mesons with spin less than or equal to two are produced",
2);
static SwitchOption interfaceTrialSpin3
(interfaceTrial,
"Spin3",
"Only hadrons with spin less than or equal to three are produced",
3);
static Switch<DarkHadronSpectrum,unsigned int> interfaceBelowThreshold
("BelowThreshold",
"Option for the selection of the hadrons if the cluster is below the pair threshold",
&DarkHadronSpectrum::belowThreshold_, 0, false, false);
static SwitchOption interfaceBelowThresholdLightest
(interfaceBelowThreshold,
"Lightest",
"Force cluster to decay to the lightest hadron with the appropriate flavours",
0);
static SwitchOption interfaceBelowThresholdAll
(interfaceBelowThreshold,
"All",
"Select from all the hadrons below the two hadron threshold according to their spin weights",
1);
}
Energy DarkHadronSpectrum::hadronPairThreshold(tcPDPtr par1, tcPDPtr par2) const {
// Determine the sum of the nominal masses of the two lightest hadrons
// with the right flavour numbers as the cluster under consideration.
// Notice that we don't need real masses (drawn by a Breit-Wigner
// distribution) because the lightest pair of hadrons does not involve
// any broad resonance.
return massLightestHadronPair(par1,par2);
}
double DarkHadronSpectrum::mixingStateWeight(long id) const {
switch(id) {
//case ParticleID::eta: return 0.5*probabilityMixing(_etamix ,1);
//case ParticleID::etaprime: return 0.5*probabilityMixing(_etamix ,2);
default: return 1./3.;
}
}
void DarkHadronSpectrum::doinit() {
if (_nlightquarks + _nheavyquarks > 9) {
throw InitException() << "Can have a maximum of 9 dark quarks!";
}
// the default partons allowed
// the quarks
for (long ix : hadronizingQuarks()) {
_partons.push_back(getParticleData(ix));
}
// the diquarks
for(long ix : hadronizingQuarks()) {
for(long iy : hadronizingQuarks()) {
// Only add diquarks if they've been defined
if(ix==iy)
if (getParticleData(makeDiquarkID(ix,iy,long(3))))
_partons.push_back(getParticleData(makeDiquarkID(ix,iy,long(3))));
else
if (getParticleData(makeDiquarkID(ix,iy,long(1))))
_partons.push_back(getParticleData(makeDiquarkID(ix,iy,long(1))));
}
}
// set the weights for the various excited mesons
// set all to one to start with
for (int l = 0; l < Lmax; ++l ) {
for (int j = 0; j < Jmax; ++j) {
for (int n = 0; n < Nmax; ++n) {
_repwt[l][j][n] = 1.0;
}
}
}
// weights for the different quarks etc
// NB: Assuming max 9 quarks, first digit used in tagging hadrons and diquarks!
vector<long> light = lightHadronizingQuarks();
for(unsigned int ix=0; ix<light.size(); ++ix){
for(unsigned int iy=0; iy<ix; ++iy){
_pwt[_DarkHadOffset + (light[ix]%10)*1000 + (light[iy]%10)*100 + 1] =
0.5 * _pwtDIquark * _pwt[light[ix]] * _pwt[light[iy]];
}
_pwt[_DarkHadOffset + (light[ix]%10)*1100 + 3] =
_pwtDIquark * _pwt[light[ix]] * _pwt[light[ix]];
}
// Set any other weights to zero
// NB: Only producing light diquarks!
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
if (_pwt.find(_partons[ix]->id()) == _pwt.end()){
_pwt[_partons[ix]->id()]=0.;
}
}
// find the maximum
map<long,double>::iterator pit =
max_element(_pwt.begin(),_pwt.end(),weightIsLess);
const double pmax = pit->second;
for(pit=_pwt.begin(); pit!=_pwt.end(); ++pit) {
pit->second/=pmax;
}
HadronSpectrum::doinit();
}
void DarkHadronSpectrum::constructHadronTable() {
// initialise the table
_table.clear();
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
for(unsigned int iy=0; iy<_partons.size(); ++iy) {
if (!(DiquarkMatcher::Check(_partons[ix]->id())
&& DiquarkMatcher::Check(_partons[iy]->id())))
_table[make_pair(_partons[ix]->id(),_partons[iy]->id())] = KupcoData();
}
}
// get the particles from the event generator
ParticleMap particles = generator()->particles();
// loop over the particles
for(ParticleMap::iterator it=particles.begin();
it!=particles.end(); ++it) {
long pid = it->first;
tPDPtr particle = it->second;
int pspin = particle->iSpin();
// Don't include hadrons which are explicitly forbidden
if(find(_forbidden.begin(),_forbidden.end(),particle)!=_forbidden.end())
continue;
// Don't include non-hadrons or antiparticles
if(pid < 100) continue;
// remove diffractive particles
if(pspin == 0) continue;
// K_0S and K_0L not made make K0 and Kbar0
if(pid==ParticleID::K_S0||pid==ParticleID::K_L0) continue;
// Debugging options
// Only include those with 2J+1 less than...5
if(_trial==2 && pspin >= 5) continue;
// Only include those with 2J+1 less than...7
if(_trial==3 && pspin >= 7) continue;
// Only include pions
if(_trial==1 && pid!=111 && pid!=211) continue;
// shouldn't be coloured
if(particle->coloured()) continue;
// Require dark hadrons
if ((pid/100000)!=49) continue;
// Get the flavours
const int x4 = (pid/1000)%10;
const int x3 = (pid/100 )%10;
const int x2 = (pid/10 )%10;
int flav1;
int flav2;
// Skip non-hadrons (susy particles, etc...)
if(x3 == 0 || x2 == 0) continue;
else if(x4 == 0) { // meson
flav1 = x2;
flav2 = x3;
}
else { // baryon
// insert the spin 1 diquark, sort out the rest later
flav1 = makeDiquarkID(x2,x3,3);
flav2 = x4;
}
insertToHadronTable(particle,flav1,flav2);
}
// normalise weights to one for first option
if(_topt == 0) {
HadronTable::const_iterator tit;
KupcoData::iterator it;
for(tit=_table.begin();tit!=_table.end();++tit) {
double weight=0;
for(it = tit->second.begin(); it!=tit->second.end(); ++it)
weight=max(weight,it->overallWeight);
weight = 1./weight;
}
}
}
void DarkHadronSpectrum::insertMeson(HadronInfo a, int flav1, int flav2) {
// identical light flavours
if(flav1 == flav2 && flav1<=_nlightquarks) {
vector<long> light = lightHadronizingQuarks();
// light quark admixture states
- a.overallWeight *= 1 / light.size();
+ a.overallWeight *= 1. / light.size();
for(unsigned int ix=0; ix<light.size(); ++ix){
_table[make_pair(light[ix],light[ix])].insert(a);
}
}
else {
- _table[make_pair(90+flav1, 90+flav2)].insert(a);
- if(flav1 != flav2) _table[make_pair(90+flav2, 90+flav1)].insert(a);
+ _table[make_pair(_DarkHadOffset+flav1, _DarkHadOffset+flav2)].insert(a);
+ if(flav1 != flav2) _table[make_pair(_DarkHadOffset+flav2, _DarkHadOffset+flav1)].insert(a);
}
}
+
+double DarkHadronSpectrum::mesonWeight(long id) const {
+ // Don't currently have radial excitations (practically clashes with dark had offset
+ // in pdgId codes; theoretically doesn't make much sense to consider such complex
+ // states). For now just return 1
+ return 1.0;
+}
+
+
long DarkHadronSpectrum::makeDiquarkID(long id1, long id2, long pspin) const {
assert(id1 * id2 > 0);
// Currently not checking if these are dark quarks to allow giving either ID
// or 1 digit flav, is this a good idea long term?
// && DarkQuarkMatcher::Check(id1)
// && DarkQuarkMatcher::Check(id2)) ;
long ida = abs(id1)%10;
long idb = abs(id2)%10;
if (ida < idb) swap(ida,idb);
if (pspin != 1 && pspin != 3) assert(false);
long idnew = ida*1000 + idb*100 + pspin + _DarkHadOffset;
// Diquarks made of quarks of the same type: uu, dd, ss, cc, bb,
// have spin 1, and therefore the less significant digit (which
// corresponds to 2*J+1) is 3 rather than 1 as all other Diquarks.
if (id1 == id2 && pspin == 1) {
//cerr<<"WARNING: spin-0 diquiark of the same type cannot exist."
// <<" Switching to spin-1 diquark.\n";
- idnew = ida*1000 + idb*100 + 3;
+ idnew = ida*1000 + idb*100 + 3 + _DarkHadOffset;
}
return id1 > 0 ? idnew : -idnew;
}
bool DarkHadronSpectrum::isExotic(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) const {
- /// \todo make this more general
- long id1 = par1 ? par1->id(): 0;
- long id2 = par2 ? par2->id(): 0;
- long id3 = par3 ? par3->id(): 0;
-return
- ( (id1/1000000)% 10 != 0 && (id1/1000000)% 10 != 9 ) ||
- ( (id2/1000000)% 10 != 0 && (id2/1000000)% 10 != 9 ) ||
- ( (id3/1000000)% 10 != 0 && (id3/1000000)% 10 != 9 ) ||
- abs(id1)==6||abs(id2)==6;
+ // Don't list dark particles as exotic to allow them to be treated as either light
+ // or heavy in the hadronisation
+ return false;
}
bool DarkHadronSpectrum::canBeBaryon(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3) const {
assert(par1 && par2);
long id1 = par1->id(), id2 = par2->id();
if (!par3) {
if( id1*id2 < 0) return false;
if(DiquarkMatcher::Check(id1))
return abs(int(par2->iColour())) == 3 && !DiquarkMatcher::Check(id2);
if(DiquarkMatcher::Check(id2))
return abs(int(par1->iColour())) == 3;
return false;
}
else {
// In this case, to be a baryon, all three components must be (anti-)quarks
// and with the same sign.
return (par1->iColour() == 3 && par2->iColour() == 3 && par3->iColour() == 3) ||
(par1->iColour() == -3 && par2->iColour() == -3 && par3->iColour() == -3);
}
}
diff --git a/Hadronization/DarkHadronSpectrum.h b/Hadronization/DarkHadronSpectrum.h
--- a/Hadronization/DarkHadronSpectrum.h
+++ b/Hadronization/DarkHadronSpectrum.h
@@ -1,360 +1,365 @@
// -*- C++ -*-
#ifndef Herwig_DarkHadronSpectrum_H
#define Herwig_DarkHadronSpectrum_H
//
// This is the declaration of the DarkHadronSpectrum class.
//
#include "Herwig/Hadronization/HadronSpectrum.h"
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/Repository/EventGenerator.h>
#include <ThePEG/PDT/EnumParticles.h>
#include "ThePEG/Repository/CurrentGenerator.h"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the DarkHadronSpectrum class.
*
* @see \ref DarkHadronSpectrumInterfaces "The interfaces"
* defined for DarkHadronSpectrum.
*/
class DarkHadronSpectrum: public HadronSpectrum {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DarkHadronSpectrum(unsigned int opt);
/**
* The destructor.
*/
virtual ~DarkHadronSpectrum();
//@}
public:
/** @name Partonic content */
//@{
/**
* Return the id of the gluon
*/
virtual long gluonId() const { return ParticleID::darkGluon; }
/**
* Return the ids of all hadronizing quarks
*/
virtual const vector<long>& hadronizingQuarks() const {
static vector<long> hadronizing = lightHadronizingQuarks();
static vector<long> heavy = heavyHadronizingQuarks();
hadronizing.insert(hadronizing.end(), heavy.begin(), heavy.end());
return hadronizing;
}
/**
* The light hadronizing quarks
*/
virtual const vector<long>& lightHadronizingQuarks() const {
if (_lightquarks.size() != _nlightquarks) {
for (long il=0; il<_nlightquarks; il++) {
- _lightquarks.push_back(il+91);
+ _lightquarks.push_back(il+_DarkHadOffset+1);
}
}
return _lightquarks;
}
/**
* The heavy hadronizing quarks
*/
virtual const vector<long>& heavyHadronizingQuarks() const {
if (_heavyquarks.size() != _nheavyquarks) {
for (long il=0; il<_nheavyquarks; il++) {
- _heavyquarks.push_back(il+91+_nlightquarks);
+ _heavyquarks.push_back(il+_DarkHadOffset+1+_nlightquarks);
}
}
return _heavyquarks;
}
/**
* The lightest quarks, used for finding the lightest Hadron Pair
*/
virtual const vector<long>& lightestQuarks() const {
// May need to be updated in future for strange-like quarks
return lightHadronizingQuarks();
}
/**
* Return true if any of the possible three input particles contains
* the indicated heavy quark. false otherwise. In the case that
* only the first particle is specified, it can be: an (anti-)quark,
* an (anti-)diquark an (anti-)meson, an (anti-)baryon; in the other
* cases, each pointer is assumed to be either (anti-)quark or
* (anti-)diquark.
*/
virtual bool hasHeavy(long id, tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const {
//ToDo: this should work for the heavyHadronizingQuarks
return false;
}
//@}
/**
* Return the threshold for a cluster to split into a pair of hadrons.
* This is normally the mass of the lightest hadron Pair, but can be
* higher for heavy and exotic clusters
*/
virtual Energy hadronPairThreshold(tcPDPtr par1, tcPDPtr par2) const;
/**
* Return the weight for the given flavour
*/
virtual double pwtQuark(const long& id) const {
return pwt(id);
}
/**
* The diquark weight.
*/
double pwtDIquark() const {
return _pwtDIquark;
}
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 Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
*
* The array _repwt is initialized using the interfaces to set different
* weights for different meson multiplets and the constructHadronTable()
* method called to complete the construction of the hadron tables.
*
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
/**
* Return the id of the diquark (anti-diquark) made by the two
* quarks (antiquarks) of id specified in input (id1, id2).
* Caller must ensure that id1 and id2 are quarks.
*/
long makeDiquarkID(long id1, long id2, long pspin) const;
/**
+ * Weights for mesons
+ */
+ virtual double mesonWeight(long id) const;
+
+ /**
* Return true, if any of the possible input particle pointer is an exotic quark, e.g. Susy quark;
* false otherwise.
*/
bool isExotic(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const;
protected:
/**
* Construct the table of hadron data
* This is the main method to initialize the hadron data (mainly the
* weights associated to each hadron, taking into account its spin,
* eventual isoscalar-octect mixing, singlet-decuplet factor). This is
* the method that one should update when new or updated hadron data is
* available.
*
* This class implements the construction of the basic table but can be
* overridden if needed in inheriting classes.
*
* The rationale for factors used for diquarks involving different quarks can
* be can be explained by taking a prototype example that in the exact SU(2) limit,
* in which:
* \f[m_u=m_d\f]
* \f[M_p=M_n=M_\Delta\f]
* and we will have equal numbers of u and d quarks produced.
* Suppose that we weight 1 the diquarks made of the same
* quark and 1/2 those made of different quarks, the fractions
* of u and d baryons (p, n, Delta) we get are the following:
* - \f$\Delta^{++}\f$: 1 possibility only u uu with weight 1
* - \f$\Delta^- \f$: 1 possibility only d dd with weight 1
* - \f$p,\Delta^+ \f$: 2 possibilities u ud with weight 1/2
* d uu with weight 1
* - \f$n,\Delta^0 \f$: 2 possibilities d ud with weight 1/2
* u dd with weight 1
* In the latter two cases, we have to take into account the
* fact that p and n have spin 1/2 whereas Delta+ and Delta0
* have spin 3/2 therefore from phase space we get a double weight
* for Delta+ and Delta0 relative to p and n respectively.
* Therefore the relative amount of these baryons that is
* produced is the following:
* # p = # n = ( 1/2 + 1 ) * 1/3 = 1/2
* # Delta++ = # Delta- = 1 = ( 1/2 + 1) * 2/3 # Delta+ = # Delta0
* which is correct, and therefore the weight 1/2 for the
* diquarks of different types of quarks is justified (at least
* in this limit of exact SU(2) ).
*/
virtual void constructHadronTable();
/**
* Access the parton weights
*/
double pwt(long pid) const {
map<long,double>::const_iterator it = _pwt.find(abs(pid));
assert( it != _pwt.end() );
return it->second;
}
/**
* Insert a meson in the table
*/
virtual void insertMeson(HadronInfo a, int flav1, int flav2);
/**
* Methods for the mixing of \f$I=0\f$ mesons
*/
//@{
/**
* Return the probability of mixing for Octet-Singlet isoscalar mixing,
* the probability of the
* \f$\frac1{\sqrt{2}}(|u\bar{u}\rangle + |d\bar{d}\rangle)\f$ component
* is returned.
* @param angleMix The mixing angle in degrees (not radians)
* @param order is 0 for no mixing, 1 for the first resonance of a pair,
* 2 for the second one.
* The mixing is defined so that for example with \f$\eta-\eta'\f$ mixing where
* the mixing angle is \f$\theta=-23^0$ with $\eta\f$ as the first particle
* and \f$\eta'\f$ the second one.
* The convention used is
* \f[\eta = \cos\theta|\eta_{\rm octet }\rangle
* -\sin\theta|\eta_{\rm singlet}\rangle\f]
* \f[\eta' = \sin\theta|\eta_{\rm octet }\rangle
* -\cos\theta|\eta_{\rm singlet}\rangle\f]
* with
* \f[|\eta_{\rm singlet}\rangle = \frac1{\sqrt{3}}
* \left[|u\bar{u}\rangle + |d\bar{d}\rangle + |s\bar{s}\rangle\right]\f]
* \f[|\eta_{\rm octet }\rangle = \frac1{\sqrt{6}}
* \left[|u\bar{u}\rangle + |d\bar{d}\rangle - 2|s\bar{s}\rangle\right]\f]
*/
double probabilityMixing(const double angleMix,
const int order) const {
static double convert=Constants::pi/180.0;
if (order == 1)
return sqr( cos( angleMix*convert + atan( sqrt(2.0) ) ) );
else if (order == 2)
return sqr( sin( angleMix*convert + atan( sqrt(2.0) ) ) );
else
return 1.;
}
/**
* Returns the weight of given mixing state.
* @param id The PDG code of the meson
*/
virtual double mixingStateWeight(long id) const;
//@}
virtual double specialQuarkWeight(double quarkWeight, long id,
const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const {
return quarkWeight;
}
/**
* The probability of producting a diquark.
*/
double _pwtDIquark;
/**
* Singlet and Decuplet weights
*/
//@{
/**
* The singlet weight
*/
double _sngWt;
/**
* The decuplet weight
*/
double _decWt;
//@}
/**
* Return true if the two or three particles in input can be the components
* of a baryon; false otherwise.
*/
virtual bool canBeBaryon(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3 = PDPtr()) const;
private:
/**
* Option for the construction of the tables
*/
unsigned int _topt;
/**
* Which particles to produce for debugging purposes
*/
unsigned int _trial;
/**
* Prefix for Dark Hadron pdgID
*/
int _DarkHadOffset = 4900000;
/**
* The number of light quarks
*/
int _nlightquarks;
/**
* The number of heavy quarks
*/
int _nheavyquarks;
/**
* The pdgIds of the light quarks
*/
mutable vector<long> _lightquarks = {};
/**
* The pdgIds of the heavy quarks
*/
mutable vector<long> _heavyquarks = {};
};
}
#endif /* Herwig_DarkHadronSpectrum_H */
diff --git a/Shower/QTilde/Base/PartnerFinder.cc b/Shower/QTilde/Base/PartnerFinder.cc
--- a/Shower/QTilde/Base/PartnerFinder.cc
+++ b/Shower/QTilde/Base/PartnerFinder.cc
@@ -1,872 +1,870 @@
// -*- C++ -*-
//
// PartnerFinder.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 PartnerFinder class.
//
#include "PartnerFinder.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
DescribeClass<PartnerFinder,Interfaced>
describePartnerFinder ("Herwig::PartnerFinder","HwShower.so");
// some useful functions to avoid using #define
namespace {
// return bool if final-state particle
inline bool FS(const tShowerParticlePtr a) {
return a->isFinalState();
}
// return colour line pointer
inline Ptr<ThePEG::ColourLine>::transient_pointer
CL(const tShowerParticlePtr a, unsigned int index=0) {
return a->colourInfo()->colourLines().empty() ? ThePEG::tColinePtr() :
const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->colourLines()[index]);
}
// return colour line size
inline size_t
CLSIZE(const tShowerParticlePtr a) {
return a->colourInfo()->colourLines().size();
}
inline Ptr<ThePEG::ColourLine>::transient_pointer
ACL(const tShowerParticlePtr a, unsigned int index=0) {
return a->colourInfo()->antiColourLines().empty() ? ThePEG::tColinePtr() :
const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->antiColourLines()[index]);
}
inline size_t
ACLSIZE(const tShowerParticlePtr a) {
return a->colourInfo()->antiColourLines().size();
}
}
void PartnerFinder::persistentOutput(PersistentOStream & os) const {
os << partnerMethod_ << QEDPartner_ << scaleChoice_;
}
void PartnerFinder::persistentInput(PersistentIStream & is, int) {
is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_;
}
void PartnerFinder::Init() {
static ClassDocumentation<PartnerFinder> documentation
("This class is responsible for finding the partners for each interaction types ",
"and within the evolution scale range specified by the ShowerVariables ",
"then to determine the initial evolution scales for each pair of partners.");
static Switch<PartnerFinder,int> interfacePartnerMethod
("PartnerMethod",
"Choice of partner finding method for gluon evolution.",
&PartnerFinder::partnerMethod_, 0, false, false);
static SwitchOption interfacePartnerMethodRandom
(interfacePartnerMethod,
"Random",
"Choose partners of a gluon randomly.",
0);
static SwitchOption interfacePartnerMethodMaximum
(interfacePartnerMethod,
"Maximum",
"Choose partner of gluon with largest angle.",
1);
static Switch<PartnerFinder,int> interfaceQEDPartner
("QEDPartner",
"Control of which particles to use as the partner for QED radiation",
&PartnerFinder::QEDPartner_, 0, false, false);
static SwitchOption interfaceQEDPartnerAll
(interfaceQEDPartner,
"All",
"Consider all possible choices which give a positive contribution"
" in the soft limit.",
0);
static SwitchOption interfaceQEDPartnerIIandFF
(interfaceQEDPartner,
"IIandFF",
"Only allow initial-initial or final-final combinations",
1);
static SwitchOption interfaceQEDPartnerIF
(interfaceQEDPartner,
"IF",
"Only allow initial-final combinations",
2);
static Switch<PartnerFinder,int> interfaceScaleChoice
("ScaleChoice",
"The choice of the evolution scales",
&PartnerFinder::scaleChoice_, 0, false, false);
static SwitchOption interfaceScaleChoicePartner
(interfaceScaleChoice,
"Partner",
"Scale of all interactions is that of the evolution partner",
0);
static SwitchOption interfaceScaleChoiceDifferent
(interfaceScaleChoice,
"Different",
"Allow each interaction to have different scales",
1);
}
void PartnerFinder::setInitialEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
ShowerInteraction type,
bool darkInteraction,
const bool setPartners) {
// clear the existing partners
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) (*cit)->clearPartners();
// set them
if(type==ShowerInteraction::QCD) {
setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
}
else if(type==ShowerInteraction::QED) {
setInitialQEDEvolutionScales(particles,isDecayCase,setPartners);
}
else if(type==ShowerInteraction::EW) {
setInitialEWEvolutionScales(particles,isDecayCase,false);
}
else if(type==ShowerInteraction::DARK) {
setInitialDARKEvolutionScales(particles,isDecayCase,setPartners);
}
else if(type==ShowerInteraction::QEDQCD) {
setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
setInitialQEDEvolutionScales(particles,isDecayCase,false);
}
else if(type==ShowerInteraction::ALL) {
setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
setInitialQEDEvolutionScales(particles,isDecayCase,false);
setInitialEWEvolutionScales(particles,isDecayCase,false);
}
else
assert(false);
if (darkInteraction) {
setInitialDARKEvolutionScales(particles,isDecayCase,setPartners);
}
// \todo EW scales here
// print out for debugging
if(Debug::level>=10) {
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) {
generator()->log() << "Particle: " << **cit << "\n";
if(!(**cit).partner()) continue;
generator()->log() << "Primary partner: " << *(**cit).partner() << "\n";
for(vector<ShowerParticle::EvolutionPartner>::const_iterator it= (**cit).partners().begin();
it!=(**cit).partners().end();++it) {
generator()->log() << static_cast<long>(it->type) << " "
<< it->weight << " "
<< it->scale/GeV << " "
<< *(it->partner)
<< "\n";
}
}
generator()->log() << flush;
}
}
void PartnerFinder::setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
// Loop over particles and consider only coloured particles which don't
// have already their colour partner fixed and that don't have children
// (the latter requirement is relaxed in the case isDecayCase is true).
// Build a map which has as key one of these particles (i.e. a pointer
// to a ShowerParticle object) and as a corresponding value the vector
// of all its possible *normal* candidate colour partners, defined as follows:
// --- have colour, and no children (this is not required in the case
// isDecayCase is true);
// --- if both are initial (incoming) state particles, then the (non-null) colourLine()
// of one of them must match the (non-null) antiColourLine() of the other.
// --- if one is an initial (incoming) state particle and the other is
// a final (outgoing) state particle, then both must have the
// same (non-null) colourLine() or the same (non-null) antiColourLine();
// Notice that this definition exclude the special case of baryon-violating
// processes (as in R-parity Susy), which will show up as particles
// without candidate colour partners, and that we will be treated a part later
// (this means that no modifications of the following loop is needed!
for ( const auto & sp : particles ) {
// Skip colourless particles
- if(!sp->data().coloured()) continue;
- if(abs(sp->id())==90 && abs(sp->id())==91 && abs(sp->id())==92) continue;
+ if(!sp->data().coloured() || sp->data().colouredInteraction() != 0) continue;
// find the partners
auto partners = findQCDPartners(sp,particles);
// must have a partner
if(partners.empty()) {
throw Exception() << "`Failed to make colour connections in "
<< "PartnerFinder::setQCDInitialEvolutionScales"
<< *sp
<< Exception::eventerror;
}
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
int position = -1;
for(size_t ix=0; ix< partners.size(); ++ix) {
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp, partners[ix].second),isDecayCase));
if (!setPartners && partners[ix].second) position = ix;
}
assert(setPartners || position >= 0);
// set partners if required
if (setPartners) {
// In the case of more than one candidate colour partners,
// there are now two approaches to choosing the partner. The
// first method is based on two assumptions:
// 1) the choice of which is the colour partner is done
// *randomly* between the available candidates.
// 2) the choice of which is the colour partner is done
// *independently* from each particle: in other words,
// if for a particle "i" its selected colour partner is
// the particle "j", then the colour partner of "j"
// does not have to be necessarily "i".
// The second method always chooses the furthest partner
// for hard gluons and gluinos.
// random choice
if( partnerMethod_ == 0 ) {
// random choice of partner
position = UseRandom::irnd(partners.size());
}
// take the one with largest angle
else if (partnerMethod_ == 1 ) {
if (sp->perturbative() == 1 &&
sp->dataPtr()->iColour()==PDT::Colour8 ) {
assert(partners.size()==2);
// Determine largest angle
double maxAngle(0.);
for(unsigned int ix=0;ix<partners.size();++ix) {
double angle = sp->momentum().vect().
angle(partners[ix].second->momentum().vect());
if(angle>maxAngle) {
maxAngle = angle;
position = ix;
}
}
}
else position = UseRandom::irnd(partners.size());
}
else assert(false);
// set the evolution partner
sp->partner(partners[position].second);
}
// primary partner set, set the others and do the scale
for(size_t ix=0; ix<partners.size(); ++ix) {
sp->addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,1.,partners[ix].first,
scales[ix].first));
}
// set scales for all interactions to that of the partner, default
Energy scale = scales[position].first;
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first==ShowerPartnerType::QCDColourLine) {
sp->scales().QCD_c =
sp->scales().QCD_c_noAO =
(scaleChoice_==0 ? scale : scales[ix].first);
}
else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
sp->scales().QCD_ac =
sp->scales().QCD_ac_noAO =
(scaleChoice_==0 ? scale : scales[ix].first);
}
else assert(false);
}
}
}
void PartnerFinder::setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
// loop over all the particles
for(const auto & sp : particles) {
// not charged or photon continue
if(!sp->dataPtr()->charged() && sp->dataPtr()->id()!=ParticleID::gamma) continue;
- if(abs(sp->id())==90 || abs(sp->id())==91 || abs(sp->id())==92) continue;
// find the potential partners
vector<pair<double,tShowerParticlePtr> > partners = findQEDPartners(sp,particles,isDecayCase);
if(partners.empty()) {
throw Exception() << "Failed to find partner in "
<< "PartnerFinder::setQEDInitialEvolutionScales"
<< *sp << Exception::eventerror;
}
// calculate the probabilities
double prob(0.);
for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
// normalise
for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
// set the partner if required
int position(-1);
// use QCD partner if set
if(!setPartners&&sp->partner()) {
for(unsigned int ix=0;ix<partners.size();++ix) {
if(sp->partner()==partners[ix].second) {
position = ix;
break;
}
}
}
// set the partner
if(setPartners||!sp->partner()||position<0) {
prob = UseRandom::rnd();
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first>prob) {
position = ix;
break;
}
prob -= partners[ix].first;
}
if(position>=0&&(setPartners||!sp->partner())) {
sp->partner(partners[position].second);
}
}
// must have a partner
if(position<0) throw Exception() << "Failed to find partner in "
<< "PartnerFinder::setQEDInitialEvolutionScales"
<< *sp << Exception::eventerror;
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
for(unsigned int ix=0;ix< partners.size();++ix) {
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp,partners[ix].second),
isDecayCase));
}
// store all the possible partners
for(unsigned int ix=0;ix<partners.size();++ix) {
sp->addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
partners[ix].first,
ShowerPartnerType::QED,
scales[ix].first));
}
// set scales
sp->scales().QED = scales[position].first;
sp->scales().QED_noAO = scales[position].first;
}
}
pair<Energy,Energy> PartnerFinder::
calculateInitialEvolutionScales(const ShowerPPair &particlePair,
const bool isDecayCase, int key) {
bool FS1=FS(particlePair.first),FS2= FS(particlePair.second);
if(FS1 && FS2){
return calculateFinalFinalScales(particlePair.first->momentum(),particlePair.second->momentum(), key);
}
else if(FS1 && !FS2) {
pair<Energy,Energy> rval = calculateInitialFinalScales(particlePair.second->momentum(),
particlePair.first->momentum(),
isDecayCase);
return { rval.second, rval.first };
}
else if(!FS1 &&FS2)
return calculateInitialFinalScales(particlePair.first->momentum(),particlePair.second->momentum(),isDecayCase);
else
return calculateInitialInitialScales(particlePair.first->momentum(),particlePair.second->momentum());
}
vector< pair<ShowerPartnerType, tShowerParticlePtr> >
PartnerFinder::findQCDPartners(tShowerParticlePtr particle,
const ShowerParticleVector &particles) {
vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners;
for(const auto & sp : particles) {
if(!sp->data().coloured() || particle==sp) continue;
// one initial-state and one final-state particle
if(FS(particle) != FS(sp)) {
// loop over all the colours of both particles
for(size_t ix=0; ix<CLSIZE(particle); ++ix) {
for(size_t jx=0; jx<CLSIZE(sp); ++jx) {
if((CL(particle,ix) && CL(particle,ix)==CL(sp,jx))) {
partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
}
}
}
//loop over all the anti-colours of both particles
for(size_t ix=0; ix<ACLSIZE(particle); ++ix) {
for(size_t jx=0; jx<ACLSIZE(sp); ++jx) {
if((ACL(particle,ix) && ACL(particle,ix)==ACL(sp,jx))) {
partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
}
}
}
}
// two initial-state or two final-state particles
else {
//loop over the colours of the first particle and the anti-colours of the other
for(size_t ix=0; ix<CLSIZE(particle); ++ix){
for(size_t jx=0; jx<ACLSIZE(sp); ++jx){
if(CL(particle,ix) && CL(particle,ix)==ACL(sp,jx)) {
partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
}
}
}
//loop over the anti-colours of the first particle and the colours of the other
for(size_t ix=0; ix<ACLSIZE(particle); ++ix){
for(size_t jx=0; jx<CLSIZE(sp); jx++){
if(ACL(particle,ix) && ACL(particle,ix)==CL(sp,jx)) {
partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
}
}
}
}
}
// if we haven't found any partners look for RPV
if (partners.empty()) {
// special for RPV
tColinePtr col = CL(particle);
if(FS(particle)&&col&&col->sourceNeighbours().first) {
tColinePair cpair = col->sourceNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))||
(!FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second ))) {
partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
}
}
}
else if(col&&col->sinkNeighbours().first) {
tColinePair cpair = col->sinkNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second))||
(!FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))) {
partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
}
}
}
col = ACL(particle);
if(FS(particle)&&col&&col->sinkNeighbours().first) {
tColinePair cpair = col->sinkNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second))||
(!FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second ))) {
partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
}
}
}
else if(col&&col->sourceNeighbours().first) {
tColinePair cpair = col->sourceNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))||
(!FS(sp) && (ACL(sp) == cpair.first ||ACL(sp) == cpair.second))) {
partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
}
}
}
}
// return the partners
return partners;
}
vector< pair<ShowerPartnerType, tShowerParticlePtr> >
PartnerFinder::findDARKPartners(tShowerParticlePtr particle,
const ShowerParticleVector &particles) {
vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners;
for (const auto & sp : particles) {
if(!sp->data().coloured() || particle==sp) continue;
- if(abs(sp->id())!=90 && abs(sp->id())!=91 && abs(sp->id())!=92) continue;
+ if(sp->data().colouredInteraction() != 1) continue;
// one initial-state and one final-state particle
if(FS(particle) != FS(sp)) {
// loop over all the colours of both particles
for(size_t ix=0; ix<CLSIZE(particle); ++ix) {
for(size_t jx=0; jx<CLSIZE(sp); ++jx) {
if((CL(particle,ix) && CL(particle,ix)==CL(sp,jx))) {
partners.push_back({ ShowerPartnerType::DARKColourLine, sp });
}
}
}//loop over all the anti-colours of both particles
for(size_t ix=0; ix<ACLSIZE(particle); ++ix) {
for(size_t jx=0; jx<ACLSIZE(sp); ++jx) {
if((ACL(particle,ix) && ACL(particle,ix)==ACL(sp,jx))) {
partners.push_back({ ShowerPartnerType::DARKAntiColourLine, sp });
}
}
}
}
// two initial-state or two final-state particles
else {
//loop over the colours of the first particle and the anti-colours of the other
for(size_t ix=0; ix<CLSIZE(particle); ++ix){
for(size_t jx=0; jx<ACLSIZE(sp); ++jx){
if(CL(particle,ix) && CL(particle,ix)==ACL(sp,jx)) {
partners.push_back({ ShowerPartnerType::DARKColourLine, sp });
}
}
}
//loop over the anti-colours of the first particle and the colours of the other
for(size_t ix=0; ix<ACLSIZE(particle); ++ix){
for(size_t jx=0; jx<CLSIZE(sp); jx++){
if(ACL(particle,ix) && ACL(particle,ix)==CL(sp,jx)) {
partners.push_back({ ShowerPartnerType::DARKAntiColourLine, sp });
}
}
}
}
}
// if we haven't found any partners look for RPV
if (partners.empty()) {
// special for RPV
tColinePtr col = CL(particle);
if(FS(particle)&&col&&col->sourceNeighbours().first) {
tColinePair cpair = col->sourceNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))||
(!FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second ))) {
partners.push_back({ ShowerPartnerType::DARKColourLine, sp });
}
}
}
else if(col&&col->sinkNeighbours().first) {
tColinePair cpair = col->sinkNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second))||
(!FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))) {
partners.push_back({ ShowerPartnerType::DARKColourLine, sp });
}
}
}
col = ACL(particle);
if(FS(particle)&&col&&col->sinkNeighbours().first) {
tColinePair cpair = col->sinkNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second))||
(!FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second ))) {
partners.push_back({ ShowerPartnerType::DARKAntiColourLine, sp });
}
}
}
else if(col&&col->sourceNeighbours().first) {
tColinePair cpair = col->sourceNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))||
(!FS(sp) && (ACL(sp) == cpair.first ||ACL(sp) == cpair.second))) {
partners.push_back({ ShowerPartnerType::DARKAntiColourLine, sp });
}
}
}
}
// return the partners
return partners;
}
void PartnerFinder::setInitialEWEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
// loop over all the particles
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) {
// if not weakly interacting continue
if( !weaklyInteracting( (**cit).dataPtr()))
continue;
// find the potential partners
vector<pair<double,tShowerParticlePtr> > partners = findEWPartners(*cit,particles,isDecayCase);
if(partners.empty()) {
throw Exception() << "Failed to find partner in "
<< "PartnerFinder::setEWInitialEvolutionScales"
<< (**cit) << Exception::eventerror;
}
// calculate the probabilities
double prob(0.);
for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
// normalise
for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
// set the partner if required
int position(-1);
// use QCD partner if set
if(!setPartners&&(*cit)->partner()) {
for(unsigned int ix=0;ix<partners.size();++ix) {
if((*cit)->partner()==partners[ix].second) {
position = ix;
break;
}
}
}
// set the partner
if(setPartners||!(*cit)->partner()||position<0) {
prob = UseRandom::rnd();
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first>prob) {
position = ix;
break;
}
prob -= partners[ix].first;
}
if(position>=0&&(setPartners||!(*cit)->partner())) {
(*cit)->partner(partners[position].second);
}
}
// must have a partner
if(position<0) throw Exception() << "Failed to find partner in "
<< "PartnerFinder::setEWInitialEvolutionScales"
<< (**cit) << Exception::eventerror;
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
for(unsigned int ix=0;ix< partners.size();++ix) {
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
isDecayCase, 0));
}
// store all the possible partners
for(unsigned int ix=0;ix<partners.size();++ix) {
(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
partners[ix].first,
ShowerPartnerType::EW,
scales[ix].first));
}
// set scales
(**cit).scales().EW = scales[position].first;
}
}
void PartnerFinder::setInitialDARKEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
for ( const auto & sp : particles ) {
// Skip colourless particles
if(!sp->data().coloured()) continue;
- if(abs(sp->id())!=90 && abs(sp->id())!=91 && abs(sp->id())!=92) continue;
+ if(sp->data().colouredInteraction() != 1) continue;
// find the partners
auto partners = findDARKPartners(sp,particles);
// must have a partner
if(partners.empty()) {
throw Exception() << "`Failed to make colour connections in "
<< "PartnerFinder::setDARKInitialEvolutionScales"
<< *sp
<< Exception::eventerror;
}
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
int position = -1;
for(size_t ix=0; ix< partners.size(); ++ix) {
scales.push_back(
calculateInitialEvolutionScales(
ShowerPPair(sp, partners[ix].second),
isDecayCase
)
);
if (!setPartners && partners[ix].second) position = ix;
}
assert(setPartners || position >= 0);
// set partners if required
if (setPartners) {
// random choice of partner
position = UseRandom::irnd(partners.size());
// set the evolution partner
sp->partner(partners[position].second);
}
// primary partner set, set the others and do the scale
for(size_t ix=0; ix<partners.size(); ++ix) {
sp->addPartner(
ShowerParticle::EvolutionPartner(
partners[ix].second,
1.,
partners[ix].first,
scales[ix].first
)
);
}
// set scales for all interactions to that of the partner, default
Energy scale = scales[position].first;
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first==ShowerPartnerType::DARKColourLine) {
sp->scales().DARK_c =
sp->scales().DARK_c_noAO =
(scaleChoice_==0 ? scale : scales[ix].first);
}
else if(partners[ix].first==ShowerPartnerType::DARKAntiColourLine) {
sp->scales().DARK_ac =
sp->scales().DARK_ac_noAO =
(scaleChoice_==0 ? scale : scales[ix].first);
}
else assert(false);
}
}
}
vector< pair<double, tShowerParticlePtr> >
PartnerFinder::findEWPartners(tShowerParticlePtr particle,
const ShowerParticleVector &particles,
const bool ) {
vector< pair<double, tShowerParticlePtr> > partners;
for(ShowerParticlePtr partner : particles) {
if(particle==partner || !weaklyInteracting(partner->dataPtr()))
continue;
partners.push_back(make_pair(1.,partner));
}
return partners;
}
vector< pair<double, tShowerParticlePtr> >
PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
const ShowerParticleVector & particles,
const bool isDecayCase) {
vector< pair<double, tShowerParticlePtr> > partners;
const double pcharge =
particle->id()==ParticleID::gamma ? 1 : double(particle->data().iCharge());
vector< pair<double, tShowerParticlePtr> > photons;
for (const auto & sp : particles) {
if (particle == sp) continue;
if (sp->id()==ParticleID::gamma) photons.push_back(make_pair(1.,sp));
if (!sp->data().charged() ) continue;
double charge = pcharge*double((sp)->data().iCharge());
if ( FS(particle) != FS(sp) ) charge *=-1.;
if ( QEDPartner_ != 0 && !isDecayCase ) {
// only include II and FF as requested
if ( QEDPartner_ == 1 && FS(particle) != FS(sp) )
continue;
// only include IF is requested
else if (QEDPartner_ == 2 && FS(particle) == FS(sp) )
continue;
}
if (particle->id()==ParticleID::gamma) charge = -abs(charge);
// only keep positive dipoles
if (charge<0.) partners.push_back({ -charge, sp });
}
if (particle->id()==ParticleID::gamma && partners.empty())
return photons;
return partners;
}
pair<Energy,Energy>
PartnerFinder::calculateFinalFinalScales(
const Lorentz5Momentum & p1,
const Lorentz5Momentum & p2, int key)
{
static const double eps=1e-7;
// Using JHEP 12(2003)045 we find that we need ktilde = 1/2(1+b-c+lambda)
// ktilde = qtilde^2/Q^2 therefore qtilde = sqrt(ktilde*Q^2)
// find momenta in rest frame of system
// calculate quantities for the scales
Energy2 Q2 = (p1+p2).m2();
double b = p1.mass2()/Q2;
double c = p2.mass2()/Q2;
if(b<0.) {
if(b<-eps) {
throw Exception() << "Negative Mass squared b = " << b
<< "in PartnerFinder::calculateFinalFinalScales()"
<< Exception::eventerror;
}
b = 0.;
}
if(c<0.) {
if(c<-eps) {
throw Exception() << "Negative Mass squared c = " << c
<< "in PartnerFinder::calculateFinalFinalScales()"
<< Exception::eventerror;
}
c = 0.;
}
// KMH & PR - 16 May 2008 - swapped lambda calculation from
// double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)),
// which should be identical for p1 & p2 onshell in their COM
// but in the inverse construction for the Nason method, this
// was not the case, leading to misuse.
const double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c))
*(sqrt(b)-1.-sqrt(c))*(sqrt(c)-1.-sqrt(b)));
Energy firstQ, secondQ;
double kappab(0.), kappac(0.);
//key = 0; // symmetric case pre-selection
switch(key) {
case 0: // symmetric case
firstQ = sqrt(0.5*Q2*(1.+b-c+lam));
secondQ = sqrt(0.5*Q2*(1.-b+c+lam));
break;
case 1: // maximum emission from both legs
kappab=4.*(1.-2.*sqrt(c)-b+c);
kappac=4.*(1.-2.*sqrt(b)-c+b);
firstQ = sqrt(Q2*kappab);
secondQ = sqrt(Q2*kappac);
break;
default:
assert(false);
}
// calculate the scales
return pair<Energy,Energy>(firstQ, secondQ);
}
pair<Energy,Energy>
PartnerFinder::calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
const bool isDecayCase) {
if(!isDecayCase) {
// In this case from JHEP 12(2003)045 we find the conditions
// ktilde_b = (1+c) and ktilde_c = (1+2c)
// We also find that c = m_c^2/Q^2. The process is a+b->c where
// particle a is not colour connected (considered as a colour singlet).
// Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and
// q_c = sqrt(Q^2+2 m_c^2)
// We also assume that the first particle in the pair is the initial
// state particle and the second is the final state one c
const Energy2 mc2 = sqr(pc.mass());
const Energy2 Q2 = -(pb-pc).m2();
return { sqrt(Q2+mc2), sqrt(Q2+2*mc2) };
}
else {
// In this case from JHEP 12(2003)045 we find, for the decay
// process b->c+a(neutral), the condition
// (ktilde_b-1)*(ktilde_c-c)=(1/4)*sqr(1-a+c+lambda).
// We also assume that the first particle in the pair is the initial
// state particle (b) and the second is the final state one (c).
// - We find maximal phase space coverage through emissions from
// c if we set ktilde_c = 4.*(sqr(1.-sqrt(a))-c)
// - We find the most 'symmetric' way to populate the phase space
// occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda)
// - We find the most 'smooth' way to populate the phase space
// occurs for...
Energy2 mb2(sqr(pb.mass()));
double a=(pb-pc).m2()/mb2;
double c=sqr(pc.mass())/mb2;
double lambda = 1. + a*a + c*c - 2.*a - 2.*c - 2.*a*c;
lambda = sqrt(max(lambda,0.));
const double PROD = 0.25*sqr(1. - a + c + lambda);
int key = 0;
double ktilde_b, ktilde_c, cosi = 0.;
switch (key) {
case 0: // the 'symmetric' choice
ktilde_c = 0.5*(1-a+c+lambda) + c ;
ktilde_b = 1.+PROD/(ktilde_c-c) ;
break;
case 1: // the 'maximal' choice
ktilde_c = 4.0*(sqr(1.-sqrt(a))-c);
ktilde_b = 1.+PROD/(ktilde_c-c) ;
break;
case 2: // the 'smooth' choice
c = max(c,1.*GeV2/mb2);
cosi = (sqr(1-sqrt(c))-a)/lambda;
ktilde_b = 2.0/(1.0-cosi);
ktilde_c = (1.0-a+c+lambda)*(1.0+c-a-lambda*cosi)/(2.0*(1.0+cosi));
break;
}
return { sqrt(mb2*ktilde_b), sqrt(mb2*ktilde_c) };
}
}
pair<Energy,Energy>
PartnerFinder::calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2) {
// This case is quite simple. From JHEP 12(2003)045 we find the condition
// that ktilde_b = ktilde_c = 1. In this case we have the process
// b+c->a so we need merely boost to the CM frame of the two incoming
// particles and then qtilde is equal to the energy in that frame
const Energy Q = sqrt((p1+p2).m2());
return {Q,Q};
}
diff --git a/Shower/QTilde/Dark/HiddenValleyAlpha.cc b/Shower/QTilde/Dark/HiddenValleyAlpha.cc
--- a/Shower/QTilde/Dark/HiddenValleyAlpha.cc
+++ b/Shower/QTilde/Dark/HiddenValleyAlpha.cc
@@ -1,328 +1,370 @@
// -*- C++ -*-
//
// HiddenValleyAlpha.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HiddenValleyAlpha class.
//
#include "HiddenValleyAlpha.h"
#include "HiddenValleyModel.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/Throw.h"
using namespace Herwig;
IBPtr HiddenValleyAlpha::clone() const {
return new_ptr(*this);
}
IBPtr HiddenValleyAlpha::fullclone() const {
return new_ptr(*this);
}
void HiddenValleyAlpha::persistentOutput(PersistentOStream & os) const {
os << _asType << _asMaxNP << ounit(_qmin,GeV) << _nloop << _thresopt
<< ounit(_lambdain,GeV) << _alphain
<< _tolerance << _maxtry << _alphamin
<< ounit(_thresholds,GeV) << ounit(_lambda,GeV) << _ca << _cf << _tr;
}
void HiddenValleyAlpha::persistentInput(PersistentIStream & is, int) {
is >> _asType >> _asMaxNP >> iunit(_qmin,GeV) >> _nloop >> _thresopt
>> iunit(_lambdain,GeV) >> _alphain
>> _tolerance >> _maxtry >> _alphamin
>> iunit(_thresholds,GeV) >> iunit(_lambda,GeV) >> _ca >> _cf >> _tr;
}
ClassDescription<HiddenValleyAlpha> HiddenValleyAlpha::initHiddenValleyAlpha;
// Definition of the static class description member.
void HiddenValleyAlpha::Init() {
static ClassDocumentation<HiddenValleyAlpha> documentation
("This (concrete) class describes the QCD alpha running.");
static Switch<HiddenValleyAlpha, int> intAsType
("NPAlphaS",
"Behaviour of AlphaS in the NP region",
&HiddenValleyAlpha::_asType, 1, false, false);
static SwitchOption intAsTypeZero
(intAsType, "Zero","zero below Q_min", 1);
static SwitchOption intAsTypeConst
(intAsType, "Const","const as(qmin) below Q_min", 2);
static SwitchOption intAsTypeLin
(intAsType, "Linear","growing linearly below Q_min", 3);
static SwitchOption intAsTypeQuad
(intAsType, "Quadratic","growing quadratically below Q_min", 4);
static SwitchOption intAsTypeExx1
(intAsType, "Exx1", "quadratic from AlphaMaxNP down to as(Q_min)", 5);
static SwitchOption intAsTypeExx2
(intAsType, "Exx2", "const = AlphaMaxNP below Q_min", 6);
// default such that as(qmin) = 1 in the current parametrization.
// min = Lambda3
static Parameter<HiddenValleyAlpha,Energy> intQmin
("Qmin", "Q < Qmin is treated with NP parametrization as of (unit [GeV])",
&HiddenValleyAlpha::_qmin, GeV, 0.630882*GeV, 0.330445*GeV,
100.0*GeV,false,false,false);
static Parameter<HiddenValleyAlpha,double> interfaceAlphaMaxNP
("AlphaMaxNP",
"Max value of alpha in NP region, only relevant if NPAlphaS = 5,6",
&HiddenValleyAlpha::_asMaxNP, 1.0, 0., 100.0,
false, false, Interface::limited);
static Parameter<HiddenValleyAlpha,unsigned int> interfaceNumberOfLoops
("NumberOfLoops",
"The number of loops to use in the alpha_S calculation",
&HiddenValleyAlpha::_nloop, 2, 1, 2,
false, false, Interface::limited);
+ static Switch<HiddenValleyAlpha,bool> interfaceLambdaOption
+ ("LambdaOption",
+ "Option for the calculation of the Lambda used in the simulation from the input"
+ " Lambda_MSbar",
+ &HiddenValleyAlpha::_lambdaopt, false, false, false);
+ static SwitchOption interfaceLambdaOptionfalse
+ (interfaceLambdaOption,
+ "Same",
+ "Use the same value",
+ false);
+ static SwitchOption interfaceLambdaOptionConvert
+ (interfaceLambdaOption,
+ "Convert",
+ "Use the conversion to the Herwig scheme from NPB349, 635",
+ true);
+
static Parameter<HiddenValleyAlpha,Energy> interfaceLambdaQCD
("LambdaQCD",
"Input value of Lambda_MSBar",
- &HiddenValleyAlpha::_lambdain, MeV, 0.208364*GeV, 100.0*MeV, 500.0*MeV,
+ &HiddenValleyAlpha::_lambdain, GeV, 0.208364*GeV, 100.0*MeV, 100.0*GeV,
false, false, Interface::limited);
static Parameter<HiddenValleyAlpha,double> interfaceAlphaMZ
("AlphaMZ",
"The input value of the strong coupling at the Z mass ",
&HiddenValleyAlpha::_alphain, 0.118, 0.1, 0.2,
false, false, Interface::limited);
+ static Switch<HiddenValleyAlpha,bool> interfaceInputOption
+ ("InputOption",
+ "Option for inputing the initial value of the coupling",
+ &HiddenValleyAlpha::_inopt, true, false, false);
+ static SwitchOption interfaceInputOptionAlphaMZ
+ (interfaceInputOption,
+ "AlphaMZ",
+ "Use the value of alpha at MZ to calculate the coupling",
+ true);
+ static SwitchOption interfaceInputOptionLambdaQCD
+ (interfaceInputOption,
+ "LambdaQCD",
+ "Use the input value of Lambda to calculate the coupling",
+ false);
+
static Parameter<HiddenValleyAlpha,double> interfaceTolerance
("Tolerance",
"The tolerance for discontinuities in alphaS at thresholds.",
&HiddenValleyAlpha::_tolerance, 1e-10, 1e-20, 1e-4,
false, false, Interface::limited);
static Parameter<HiddenValleyAlpha,unsigned int> interfaceMaximumIterations
("MaximumIterations",
"The maximum number of iterations for the Newton-Raphson method to converge.",
&HiddenValleyAlpha::_maxtry, 100, 10, 1000,
false, false, Interface::limited);
static Switch<HiddenValleyAlpha,bool> interfaceThresholdOption
("ThresholdOption",
"Whether to use the consistuent or normal masses for the thresholds",
&HiddenValleyAlpha::_thresopt, true, false, false);
static SwitchOption interfaceThresholdOptionCurrent
(interfaceThresholdOption,
"Current",
"Use the current masses",
true);
static SwitchOption interfaceThresholdOptionConstituent
(interfaceThresholdOption,
"Constituent",
"Use the constitent masses.",
false);
-
}
double HiddenValleyAlpha::value(const Energy2 scale) const {
pair<short,Energy> nflam;
Energy q = sqrt(scale);
double val(0.);
// special handling if the scale is less than Qmin
if (q < _qmin) {
nflam = getLamNfTwoLoop(_qmin);
double val0 = alphaS(_qmin, nflam.second, nflam.first);
switch (_asType) {
case 1:
// flat, zero; the default type with no NP effects.
val = 0.;
break;
case 2:
// flat, non-zero alpha_s = alpha_s(q2min).
val = val0;
break;
case 3:
// linear in q
val = val0*q/_qmin;
break;
case 4:
// quadratic in q
val = val0*sqr(q/_qmin);
break;
case 5:
// quadratic in q, starting off at asMaxNP, ending on as(qmin)
val = (val0 - _asMaxNP)*sqr(q/_qmin) + _asMaxNP;
break;
case 6:
// just asMaxNP and constant
val = _asMaxNP;
break;
}
}
else {
// the 'ordinary' case
nflam = getLamNfTwoLoop(q);
val = alphaS(q, nflam.second, nflam.first);
}
return scaleFactor() * val;
}
double HiddenValleyAlpha::overestimateValue() const {
return scaleFactor() * _alphamin;
}
double HiddenValleyAlpha::ratio(const Energy2 scale,double) const {
pair<short,Energy> nflam;
Energy q = sqrt(scale);
double val(0.);
// special handling if the scale is less than Qmin
if (q < _qmin) {
nflam = getLamNfTwoLoop(_qmin);
double val0 = alphaS(_qmin, nflam.second, nflam.first);
switch (_asType) {
case 1:
// flat, zero; the default type with no NP effects.
val = 0.;
break;
case 2:
// flat, non-zero alpha_s = alpha_s(q2min).
val = val0;
break;
case 3:
// linear in q
val = val0*q/_qmin;
break;
case 4:
// quadratic in q
val = val0*sqr(q/_qmin);
break;
case 5:
// quadratic in q, starting off at asMaxNP, ending on as(qmin)
val = (val0 - _asMaxNP)*sqr(q/_qmin) + _asMaxNP;
break;
case 6:
// just asMaxNP and constant
val = _asMaxNP;
break;
}
} else {
// the 'ordinary' case
nflam = getLamNfTwoLoop(q);
val = alphaS(q, nflam.second, nflam.first);
}
// denominator
return val/_alphamin;
}
Energy HiddenValleyAlpha::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))/derivativealphaS(match,lamtest,nflav);
lamtest=match/exp(0.5*xtest);
}
while(abs(alpha-alphaS(match,lamtest,nflav)) > _tolerance && ntry < _maxtry);
return lamtest;
}
pair<short, Energy> HiddenValleyAlpha::getLamNfTwoLoop(Energy q) const {
unsigned int ix=1;
- for(;ix<_thresholds.size();++ix) {
+ for(;ix<_thresholds.size();ix++) {
if(q<_thresholds[ix]) break;
}
- if(ix==_thresholds.size()) --ix;
--ix;
- return pair<short,Energy>(ix, _lambda[ix]);
+ return pair<short,Energy>(ix+_nf_light, _lambda[ix]);
}
void HiddenValleyAlpha::doinit() {
ShowerAlpha::doinit();
// get the model for parameters
tcHiddenValleyPtr model = dynamic_ptr_cast<tcHiddenValleyPtr>
(generator()->standardModel());
if(!model) throw InitException() << "Must be using the HiddenValleyModel"
<< " in HiddenValleyAlpha::doinit()"
<< Exception::runerror;
// get the colour factors
_ca = model->CA();
_cf = model->CF();
_tr = model->TR();
// get the thresholds
_thresholds.push_back(_qmin);
+ _nf_light = 0;
for(unsigned int ix=1;ix<=model->NF();++ix) {
- _thresholds.push_back(getParticleData(ParticleID::darkGluon+long(ix))->mass());
+ Energy qmass = getParticleData(ParticleID::darkGluon+long(ix))->mass();
+ if (qmass > _qmin) _thresholds.push_back(qmass);
+ else _nf_light++;
}
_lambda.resize(_thresholds.size());
Energy mz = getParticleData(ThePEG::ParticleID::Z0)->mass();
- unsigned int nf;
- for(nf=0;nf<_thresholds.size();++nf) {
- if(mz<_thresholds[nf]) break;
+ unsigned int nf_heavy;
+ for(nf_heavy=0;nf_heavy<_thresholds.size();++nf_heavy) {
+ if(mz<_thresholds[nf_heavy]) break;
}
- nf-=1;
+ nf_heavy-=1;
+ unsigned int nf = _nf_light+nf_heavy;
+
// value of Lambda from alphas if needed using Newton-Raphson
- _lambda[nf] = computeLambda(mz,_alphain,nf-1);
+ if(_inopt) {
+ _lambda[nf_heavy] = computeLambda(mz,_alphain,nf-1);
+ }
+ // otherwise it was an input parameter
+ else{_lambda[nf_heavy]=_lambdain;}
+ // convert lambda to the Monte Carlo scheme if needed
+ using Constants::pi;
+ if(_lambdaopt) _lambda[nf_heavy] *=exp((0.5*_ca*(67./3.-sqr(pi))-_tr*nf*10./3.)/(11*_ca-2*nf))/sqrt(2.);
+
// compute the threshold matching
// above the Z mass
- for(int ix=nf;ix<int(_thresholds.size())-1;++ix) {
+ for(int ix=nf_heavy;ix<int(_thresholds.size())-1;++ix) {
_lambda[ix+1] = computeLambda(_thresholds[ix+1],alphaS(_thresholds[ix+1],
_lambda[ix],ix),ix+1);
}
// below Z mass
- for(int ix=nf-1;ix>=0;--ix) {
+ for(int ix=nf_heavy-1;ix>=0;--ix) {
_lambda[ix] = computeLambda(_thresholds[ix+1],alphaS(_thresholds[ix+1],
_lambda[ix+1],ix+1),ix);
}
// compute the maximum value of as
if ( _asType < 5 ) _alphamin = value(sqr(_qmin)+1.0e-8*sqr(MeV)); // approx as = 1
else _alphamin = max(_asMaxNP, value(sqr(_qmin)+1.0e-8*sqr(MeV)));
// check consistency lambda_3 < qmin
if(_lambda[0]>_qmin)
Throw<InitException>() << "The value of Qmin is less than Lambda in "
<< _qmin/GeV << " < " << _lambda[0]/GeV
<< " HiddenValleyAlpha::doinit " << Exception::abortnow;
}
double HiddenValleyAlpha::derivativealphaS(Energy q, Energy lam, int nf) const {
using Constants::pi;
double lx = log(sqr(q/lam));
// N.B. b_1 is divided by 2 due Hw++ convention
double b0 = 11./3.*_ca - 4./3.*_tr*nf;
double b1 = 17./3.*sqr(_ca) - nf*_tr*(10./3.*_ca+2.*_cf);
if(_nloop==1)
return -4.*pi/(b0*sqr(lx));
else if(_nloop==2)
return -4.*pi/(b0*sqr(lx))*(1.-2.*b1/sqr(b0)/lx*(1.-2.*log(lx)));
else
assert(false);
return 0.;
}
double HiddenValleyAlpha::alphaS(Energy q, Energy lam, int nf) const {
using Constants::pi;
double lx(log(sqr(q/lam)));
// N.B. b_1 is divided by 2 due Hw++ convention
double b0 = 11./3.*_ca - 4./3.*_tr*nf;
double b1 = 17./3.*sqr(_ca) - nf*_tr*(10./3.*_ca+2.*_cf);
// one loop
if(_nloop==1)
return 4.*pi/(b0*lx);
// two loop
else if(_nloop==2) {
return 4.*pi/(b0*lx)*(1.-2.*b1/sqr(b0)*log(lx)/lx);
}
else
assert(false);
return 0.;
}
diff --git a/Shower/QTilde/Dark/HiddenValleyAlpha.h b/Shower/QTilde/Dark/HiddenValleyAlpha.h
--- a/Shower/QTilde/Dark/HiddenValleyAlpha.h
+++ b/Shower/QTilde/Dark/HiddenValleyAlpha.h
@@ -1,323 +1,338 @@
// -*- C++ -*-
//
// HiddenValleyAlpha.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_HiddenValleyAlpha_H
#define HERWIG_HiddenValleyAlpha_H
//
// This is the declaration of the HiddenValleyAlpha class.
//
#include "Herwig/Shower/ShowerAlpha.h"
#include "ThePEG/Config/Constants.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
*
* This concrete class provides the definition of the
* pure virtual function value() and overestimateValue() for the
* strong coupling.
*
* A number of different options for the running of the coupling
* and its initial definition are supported.
*
* @see \ref HiddenValleyAlphaInterfaces "The interfaces"
* defined for HiddenValleyAlpha.
*/
class HiddenValleyAlpha: public ShowerAlpha {
public:
/**
* The default constructor.
*/
HiddenValleyAlpha() : ShowerAlpha(),
_qmin(0.630882*GeV), _asType(1), _asMaxNP(1.0),
_nloop(2),_thresopt(false),
_lambdain(0.208364*GeV),_alphain(0.118),
_tolerance(1e-10),
_maxtry(100),_alphamin(0.) {}
public:
/**
* Methods to return the coupling
*/
//@{
/**
* It returns the running coupling value evaluated at the input scale
* multiplied by the scale factor scaleFactor().
* @param scale The scale
* @return The coupling
*/
virtual double value(const Energy2 scale) const;
/**
* It returns the running coupling value evaluated at the input scale
* multiplied by the scale factor scaleFactor().
*/
virtual double overestimateValue() const;
/**
* Return the ratio of the coupling at the scale to the overestimated value
*/
virtual double ratio(const Energy2 scale,double factor=1.) const;
/**
* Initialize this coupling.
*/
virtual void initialize() { doinit(); }
//@}
/**
* Get the value of \f$\Lambda_{\rm QCd}\f$
* @param nf number of flavours
*/
Energy lambdaQCD(unsigned int nf) {
if (nf <= 3) return _lambda[0];
else if (nf==4 || nf==5) return _lambda[nf-3];
else return _lambda[3];
}
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:
/**
* Member functions which calculate the coupling
*/
//@{
/**
* The 1,2,3-loop parametrization of \f$\alpha_S\f$.
* @param q The scale
* @param lam \f$\Lambda_{\rm QCD}\f$
* @param nf The number of flavours
*/
double alphaS(Energy q, Energy lam, int nf) const;
/**
* The derivative of \f$\alpha_S\f$ with respect to \f$\ln(Q^2/\Lambda^2)\f$
* @param q The scale
* @param lam \f$\Lambda_{\rm QCD}\f$
* @param nf The number of flavours
*/
double derivativealphaS(Energy q, Energy lam, int nf) const;
/**
* Compute the value of \f$Lambda\f$ needed to get the input value of
* the strong coupling at the scale given for the given number of flavours
* using the Newton-Raphson method
* @param match The scale for the coupling
* @param alpha The input coupling
* @param nflav The number of flavours
*/
Energy computeLambda(Energy match, double alpha, unsigned int nflav) const;
/**
* Return the value of \f$\Lambda\f$ and the number of flavours at the scale.
* @param q The scale
* @return The number of flavours at the scale and \f$\Lambda\f$.
*/
pair<short, Energy> getLamNfTwoLoop(Energy q) const;
//@}
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<HiddenValleyAlpha> initHiddenValleyAlpha;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
HiddenValleyAlpha & operator=(const HiddenValleyAlpha &);
private:
/**
* Minimum value of the scale
*/
Energy _qmin;
/**
* Parameter controlling the behaviour of \f$\alpha_S\f$ in the
* non-perturbative region.
*/
int _asType;
/**
* Another parameter, a possible (maximum) value of alpha in the
* non-perturbative region.
*/
double _asMaxNP;
/**
* Thresholds for the different number of flavours
*/
vector<Energy> _thresholds;
/**
* \f$\Lambda\f$ for the different number of flavours
*/
vector<Energy> _lambda;
/**
* Option for the number of loops
*/
unsigned int _nloop;
/**
* Option for the threshold masses
*/
bool _thresopt;
/**
* Input value of Lambda
*/
Energy _lambdain;
/**
* Input value of \f$alpha_S(M_Z)\f$
*/
double _alphain;
/**
+ * Whether to convert lambda from MSbar scheme
+ */
+ bool _lambdaopt;
+
+ /**
+ * Whether to use input value of alphas(MZ) or lambda
+ */
+ bool _inopt;
+
+ /**
* Tolerance for discontinuities at the thresholds
*/
double _tolerance;
/**
* Maximum number of iterations for the Newton-Raphson method to converge
*/
unsigned int _maxtry;
/**
+ * Number of light flavours (below qmin)
+ */
+ unsigned int _nf_light;
+
+ /**
* The minimum value of the coupling
*/
double _alphamin;
/**
* Colour factors
*/
//@{
/**
* \f$C_A\f$
*/
double _ca;
/**
* \f$C_A\f$
*/
double _cf;
/**
* \f$T_R\f$
*/
double _tr;
//@}
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of HiddenValleyAlpha. */
template <>
struct BaseClassTrait<Herwig::HiddenValleyAlpha,1> {
/** Typedef of the first base class of HiddenValleyAlpha. */
typedef Herwig::ShowerAlpha NthBase;
};
/** This template specialization informs ThePEG about the name of
* the HiddenValleyAlpha class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::HiddenValleyAlpha>
: public ClassTraitsBase<Herwig::HiddenValleyAlpha> {
/** Return a platform-independent class name */
static string className() { return "Herwig::HiddenValleyAlpha"; }
/**
* The name of a file containing the dynamic library where the class
* HiddenValleyAlpha is implemented. It may also include several,
* space-separated, libraries if the class HiddenValleyAlpha depends on
* other classes (base classes excepted). In this case the listed
* libraries will be dynamically linked in the order they are
* specified.
*/
static string library() { return "HwShower.so HwHiddenValley.so"; }
};
/** @endcond */
}
#endif /* HERWIG_HiddenValleyAlpha_H */
diff --git a/Shower/QTilde/Dark/HiddenValleyModel.cc b/Shower/QTilde/Dark/HiddenValleyModel.cc
--- a/Shower/QTilde/Dark/HiddenValleyModel.cc
+++ b/Shower/QTilde/Dark/HiddenValleyModel.cc
@@ -1,130 +1,130 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HiddenValleyModel class.
//
#include "HiddenValleyModel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
HiddenValleyModel::HiddenValleyModel() : groupType_(SU), Nc_(3), Nf_(1),
qL_(-0.2), uR_(-0.2), dR_(0.6),
lL_(0.6), lR_(-0.2), gPrime_(1./7.),
qCharge_(1,-0.2)
{}
IBPtr HiddenValleyModel::clone() const {
return new_ptr(*this);
}
IBPtr HiddenValleyModel::fullclone() const {
return new_ptr(*this);
}
void HiddenValleyModel::persistentOutput(PersistentOStream & os) const {
os << oenum(groupType_) << Nc_ << Nf_ << FFZPVertex_
<< qL_ << uR_ << dR_ << lL_ << lR_ << qCharge_ << gPrime_;
}
void HiddenValleyModel::persistentInput(PersistentIStream & is, int) {
is >> ienum(groupType_) >> Nc_ >> Nf_ >> FFZPVertex_
>> qL_ >> uR_ >> dR_ >> lL_ >> lR_ >> qCharge_ >> gPrime_;
}
ClassDescription<HiddenValleyModel> HiddenValleyModel::initHiddenValleyModel;
// Definition of the static class description member.
void HiddenValleyModel::Init() {
static ClassDocumentation<HiddenValleyModel> documentation
("The HiddenValleyModel class is the main class for the Hidden Valley Model.");
static Switch<HiddenValleyModel,GroupType> interfaceGroupType
("GroupType",
"Type of the new unbroken group",
&HiddenValleyModel::groupType_, SU, false, false);
static SwitchOption interfaceGroupTypeSU
(interfaceGroupType,
"SU",
"Use an SU(N) group",
SU);
static Parameter<HiddenValleyModel,unsigned int> interfaceGroupOrder
("GroupOrder",
"The value of the number of 'colours' for the new symmetry group",
&HiddenValleyModel::Nc_, 3, 1, 1000,
false, false, Interface::limited);
static Parameter<HiddenValleyModel,unsigned int> interfaceNumberOfFermions
("NumberOfFermions",
"The number of fermions charged under the new group",
&HiddenValleyModel::Nf_, 1, 1, 100,
false, false, Interface::limited);
static Reference<HiddenValleyModel,AbstractFFVVertex> interfaceVertexFFZPrime
("Vertex/FFZPrime",
"The vertex coupling the Zprime to the fermions",
&HiddenValleyModel::FFZPVertex_, false, false, true, false, false);
static ParVector<HiddenValleyModel,double> interfaceQuirkCharges
("QuirkCharges",
"The charges under the new U(1) of the new quirks",
- &HiddenValleyModel::qCharge_, 1, -0.2, -10., 10.0,
+ &HiddenValleyModel::qCharge_, -1, -0.2, -10., 10.0,
false, false, Interface::limited);
static Parameter<HiddenValleyModel,double> interfaceZPrimeQL
("ZPrime/QL",
"The charge of the left-handed quarks under the new U(1)",
&HiddenValleyModel::qL_, -0.2, 10.0, 10.0,
false, false, Interface::limited);
static Parameter<HiddenValleyModel,double> interfaceZPrimeUR
("ZPrime/UR",
"The charge of the right-handed up quarks under the new U(1)",
&HiddenValleyModel::uR_, -0.2, 10.0, 10.0,
false, false, Interface::limited);
static Parameter<HiddenValleyModel,double> interfaceZPrimeDR
("ZPrime/DR",
"The charge of the right-handed down quarks under the new U(1)",
&HiddenValleyModel::uR_, 0.6, 10.0, 10.0,
false, false, Interface::limited);
static Parameter<HiddenValleyModel,double> interfaceZPrimeLL
("ZPrime/LL",
"The charge of the left-handed leptons under the new U(1)",
&HiddenValleyModel::lL_, 0.6, 10.0, 10.0,
false, false, Interface::limited);
static Parameter<HiddenValleyModel,double> interfaceZPrimeLR
("ZPrime/LR",
"The charge of the right-handed charged leptons under the new U(1)",
&HiddenValleyModel::lR_, -0.2, 10.0, 10.0,
false, false, Interface::limited);
static Parameter<HiddenValleyModel,double> interfaceGPrime
("GPrime",
"The new gauge coupling",
&HiddenValleyModel::gPrime_, 1./7., 0.0, 10.0,
false, false, Interface::limited);
}
void HiddenValleyModel::doinit() {
addVertex(FFZPVertex_);
BSMModel::doinit();
FFZPVertex_->init();
if(qCharge_.size()!=Nf_)
- throw InitException() << "Number of new fermions and number of new charges"
+ throw InitException() << "Number of new fermions and number of new charges "
<< "different in HiddenValleyModel::doinit()"
<< Exception::runerror;
}
diff --git a/Shower/QTilde/SplittingFunctions/PTCutOff.cc b/Shower/QTilde/SplittingFunctions/PTCutOff.cc
--- a/Shower/QTilde/SplittingFunctions/PTCutOff.cc
+++ b/Shower/QTilde/SplittingFunctions/PTCutOff.cc
@@ -1,66 +1,66 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the PTCutOff class.
//
#include "PTCutOff.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/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IBPtr PTCutOff::clone() const {
return new_ptr(*this);
}
IBPtr PTCutOff::fullclone() const {
return new_ptr(*this);
}
void PTCutOff::persistentOutput(PersistentOStream & os) const {
os << ounit(pTmin_,GeV) << ounit(pT2min_,GeV2);
}
void PTCutOff::persistentInput(PersistentIStream & is, int) {
is >> iunit(pTmin_,GeV) >> iunit(pT2min_,GeV2);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<PTCutOff,SudakovCutOff>
describeHerwigPTCutOff("Herwig::PTCutOff", "HwShower.so");
void PTCutOff::Init() {
static ClassDocumentation<PTCutOff> documentation
("There is no documentation for the PTCutOff class");
static Parameter<PTCutOff,Energy> interfacepTmin
("pTmin",
"The minimum pT if using a cut-off on the pT",
- &PTCutOff::pTmin_, GeV, 1.0*GeV, ZERO, 10.0*GeV,
+ &PTCutOff::pTmin_, GeV, 1.0*GeV, ZERO, 100.0*GeV,
false, false, Interface::limited);
}
void PTCutOff::doinit() {
pT2min_ = sqr(pTmin_);
SudakovCutOff::doinit();
}
const vector<Energy> & PTCutOff::virtualMasses(const IdList & ids) {
static vector<Energy> output;
output.clear();
for(auto id : ids)
output.push_back(id->mass());
return output;
}
diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in
--- a/src/defaults/Shower.in
+++ b/src/defaults/Shower.in
@@ -1,501 +1,501 @@
# -*- ThePEG-repository -*-
############################################################
# Setup of default parton shower
#
# Useful switches for users are marked near the top of
# this file.
#
# Don't edit this file directly, but reset the switches
# in your own input files!
############################################################
library HwMPI.so
library HwShower.so
library HwMatching.so
mkdir /Herwig/Shower
cd /Herwig/Shower
create Herwig::QTildeShowerHandler ShowerHandler
newdef ShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
newdef ShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
# use LO PDFs for Shower, can be changed later
newdef ShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
newdef ShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
newdef ShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
newdef ShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
#####################################
# initial setup, don't change these!
#####################################
create Herwig::SplittingGenerator SplittingGenerator
create Herwig::ShowerAlphaQCD AlphaQCD
create Herwig::ShowerAlphaQCD AlphaQCDFSR
create Herwig::ShowerAlphaQED AlphaQED
set AlphaQED:CouplingSource Thompson
create Herwig::ShowerAlphaQED AlphaEW
set AlphaEW:CouplingSource MZ
create Herwig::HiddenValleyAlpha AlphaDARK
create Herwig::PartnerFinder PartnerFinder
newdef PartnerFinder:PartnerMethod 0
newdef PartnerFinder:ScaleChoice 0
create Herwig::KinematicsReconstructor KinematicsReconstructor
newdef KinematicsReconstructor:ReconstructionOption Colour3
newdef KinematicsReconstructor:InitialStateReconOption SofterFraction
newdef KinematicsReconstructor:InitialInitialBoostOption LongTransBoost
newdef KinematicsReconstructor:FinalFinalWeight Yes
newdef /Herwig/Partons/RemnantDecayer:AlphaS AlphaQCD
newdef /Herwig/Partons/RemnantDecayer:AlphaEM AlphaQED
newdef ShowerHandler:PartnerFinder PartnerFinder
newdef ShowerHandler:KinematicsReconstructor KinematicsReconstructor
newdef ShowerHandler:SplittingGenerator SplittingGenerator
-newdef ShowerHandler:Interactions QEDQCD # TODO: is this a good choice?
+newdef ShowerHandler:Interactions QEDQCD
newdef ShowerHandler:SpinCorrelations Yes
newdef ShowerHandler:SoftCorrelations Singular
##################################################################
# Intrinsic pT
#
# Recommended:
# 1.9 GeV for Tevatron W/Z production.
# 2.1 GeV for LHC W/Z production at 10 TeV
# 2.2 GeV for LHC W/Z production at 14 TeV
#
# Set all parameters to 0 to disable
##################################################################
newdef ShowerHandler:IntrinsicPtGaussian 2.2*GeV
newdef ShowerHandler:IntrinsicPtBeta 0
newdef ShowerHandler:IntrinsicPtGamma 0*GeV
newdef ShowerHandler:IntrinsicPtIptmax 0*GeV
#############################################################
# Set up truncated shower handler.
#############################################################
create Herwig::PowhegShowerHandler PowhegShowerHandler
set PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
set PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
newdef PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:PartnerFinder PartnerFinder
newdef PowhegShowerHandler:KinematicsReconstructor KinematicsReconstructor
newdef PowhegShowerHandler:SplittingGenerator SplittingGenerator
newdef PowhegShowerHandler:Interactions QEDQCD
newdef PowhegShowerHandler:SpinCorrelations Yes
newdef PowhegShowerHandler:SoftCorrelations Singular
newdef PowhegShowerHandler:IntrinsicPtGaussian 2.2*GeV
newdef PowhegShowerHandler:IntrinsicPtBeta 0
newdef PowhegShowerHandler:IntrinsicPtGamma 0*GeV
newdef PowhegShowerHandler:IntrinsicPtIptmax 0*GeV
newdef PowhegShowerHandler:EvolutionScheme DotProduct
#############################################################
# End of interesting user servicable section.
#
# Anything that follows below should only be touched if you
# know what you're doing.
#
# Really.
#############################################################
#
# a few default values
newdef ShowerHandler:MECorrMode 1
newdef ShowerHandler:EvolutionScheme DotProduct
newdef AlphaQCD:ScaleFactor 1.0
newdef AlphaQCD:NPAlphaS 2
newdef AlphaQCD:Qmin 0.935
newdef AlphaQCD:NumberOfLoops 2
newdef AlphaQCD:AlphaIn 0.126
newdef AlphaQCDFSR:ScaleFactor 1.0
newdef AlphaQCDFSR:NPAlphaS 2
newdef AlphaQCDFSR:Qmin 0.935
newdef AlphaQCDFSR:NumberOfLoops 2
newdef AlphaQCDFSR:AlphaIn 0.1186
newdef AlphaDARK:ScaleFactor 1.0
newdef AlphaDARK:NPAlphaS 2
newdef AlphaDARK:Qmin 0.935
newdef AlphaDARK:NumberOfLoops 2
#newdef AlphaDARK:AlphaIn 0.1186
#
#
#
# Lets set up all the splittings
create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn
set QtoQGammaSplitFn:InteractionType QED
set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral
set QtoQGammaSplitFn:AngularOrdered Yes
set QtoQGammaSplitFn:StrictAO Yes
create Herwig::HalfHalfOneSplitFn QtoQGSplitFn
newdef QtoQGSplitFn:InteractionType QCD
newdef QtoQGSplitFn:ColourStructure TripletTripletOctet
set QtoQGSplitFn:AngularOrdered Yes
set QtoQGSplitFn:StrictAO Yes
create Herwig::HalfHalfOneDarkSplitFn QtoQGDARKSplitFn
newdef QtoQGDARKSplitFn:InteractionType DARK
newdef QtoQGDARKSplitFn:ColourStructure TripletTripletOctet
set QtoQGDARKSplitFn:AngularOrdered Yes
set QtoQGDARKSplitFn:StrictAO Yes
create Herwig::OneOneOneSplitFn GtoGGSplitFn
newdef GtoGGSplitFn:InteractionType QCD
newdef GtoGGSplitFn:ColourStructure OctetOctetOctet
set GtoGGSplitFn:AngularOrdered Yes
set GtoGGSplitFn:StrictAO Yes
create Herwig::OneOneOneDarkSplitFn GtoGGDARKSplitFn
newdef GtoGGDARKSplitFn:InteractionType DARK
newdef GtoGGDARKSplitFn:ColourStructure OctetOctetOctet
set GtoGGDARKSplitFn:AngularOrdered Yes
set GtoGGDARKSplitFn:StrictAO Yes
create Herwig::OneOneOneMassiveSplitFn WtoWGammaSplitFn
newdef WtoWGammaSplitFn:InteractionType QED
newdef WtoWGammaSplitFn:ColourStructure ChargedChargedNeutral
set WtoWGammaSplitFn:AngularOrdered Yes
set WtoWGammaSplitFn:StrictAO Yes
create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn
newdef GtoQQbarSplitFn:InteractionType QCD
newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet
set GtoQQbarSplitFn:AngularOrdered Yes
set GtoQQbarSplitFn:StrictAO Yes
create Herwig::OneHalfHalfDarkSplitFn GtoQQbarDARKSplitFn
newdef GtoQQbarDARKSplitFn:InteractionType DARK
newdef GtoQQbarDARKSplitFn:ColourStructure OctetTripletTriplet
set GtoQQbarDARKSplitFn:AngularOrdered Yes
set GtoQQbarDARKSplitFn:StrictAO Yes
create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn
newdef GammatoQQbarSplitFn:InteractionType QED
newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged
set GammatoQQbarSplitFn:AngularOrdered Yes
set GammatoQQbarSplitFn:StrictAO Yes
create Herwig::HalfOneHalfSplitFn QtoGQSplitFn
newdef QtoGQSplitFn:InteractionType QCD
newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet
set QtoGQSplitFn:AngularOrdered Yes
set QtoGQSplitFn:StrictAO Yes
create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn
newdef QtoGammaQSplitFn:InteractionType QED
newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged
set QtoGammaQSplitFn:AngularOrdered Yes
set QtoGammaQSplitFn:StrictAO Yes
create Herwig::HalfHalfOneEWSplitFn QtoQWZSplitFn
newdef QtoQWZSplitFn:InteractionType EW
newdef QtoQWZSplitFn:ColourStructure EW
create Herwig::HalfHalfOneEWSplitFn LtoLWZSplitFn
newdef LtoLWZSplitFn:InteractionType EW
newdef LtoLWZSplitFn:ColourStructure EW
create Herwig::HalfHalfZeroEWSplitFn QtoQHSplitFn
newdef QtoQHSplitFn:InteractionType EW
newdef QtoQHSplitFn:ColourStructure EW
create Herwig::OneOneOneEWSplitFn VtoVVSplitFn
newdef VtoVVSplitFn:InteractionType EW
newdef VtoVVSplitFn:ColourStructure EW
create Herwig::OneOneOneQEDSplitFn GammatoWWSplitFn
newdef GammatoWWSplitFn:InteractionType QED
newdef GammatoWWSplitFn:ColourStructure ChargedNeutralCharged
create Herwig::OneOneZeroEWSplitFn VtoVHSplitFn
newdef VtoVHSplitFn:InteractionType EW
newdef VtoVHSplitFn:ColourStructure EW
#
# Now the Sudakovs
create Herwig::PTCutOff PTCutOff
newdef PTCutOff:pTmin 0.958*GeV
create Herwig::SudakovFormFactor SudakovCommon
newdef SudakovCommon:Alpha AlphaQCD
newdef SudakovCommon:Cutoff PTCutOff
newdef SudakovCommon:PDFmax 1.0
cp SudakovCommon QtoQGSudakov
newdef QtoQGSudakov:SplittingFunction QtoQGSplitFn
newdef QtoQGSudakov:PDFmax 1.9
cp SudakovCommon QtoQGDARKSudakov
newdef QtoQGDARKSudakov:SplittingFunction QtoQGDARKSplitFn
newdef QtoQGDARKSudakov:PDFmax 1.9
cp SudakovCommon QtoQGammaSudakov
set QtoQGammaSudakov:SplittingFunction QtoQGammaSplitFn
set QtoQGammaSudakov:Alpha AlphaQED
set QtoQGammaSudakov:PDFmax 1.9
cp QtoQGammaSudakov LtoLGammaSudakov
cp PTCutOff LtoLGammaPTCutOff
# Technical parameter to stop evolution.
set LtoLGammaPTCutOff:pTmin 0.000001
set LtoLGammaSudakov:Cutoff LtoLGammaPTCutOff
cp SudakovCommon QtoQWZSudakov
set QtoQWZSudakov:SplittingFunction QtoQWZSplitFn
set QtoQWZSudakov:Alpha AlphaEW
set QtoQWZSudakov:PDFmax 1.9
cp SudakovCommon LtoLWZSudakov
set LtoLWZSudakov:SplittingFunction LtoLWZSplitFn
set LtoLWZSudakov:Alpha AlphaEW
set LtoLWZSudakov:PDFmax 1.9
cp SudakovCommon QtoQHSudakov
set QtoQHSudakov:SplittingFunction QtoQHSplitFn
set QtoQHSudakov:Alpha AlphaEW
set QtoQHSudakov:PDFmax 1.9
cp QtoQHSudakov LtoLHSudakov
cp SudakovCommon VtoVVSudakov
set VtoVVSudakov:SplittingFunction VtoVVSplitFn
set VtoVVSudakov:Alpha AlphaQED
cp SudakovCommon GammatoWWSudakov
set GammatoWWSudakov:SplittingFunction GammatoWWSplitFn
set GammatoWWSudakov:Alpha AlphaEW
set GammatoWWSudakov:PDFmax 1.9
cp SudakovCommon VtoVHSudakov
set VtoVHSudakov:SplittingFunction VtoVHSplitFn
set VtoVHSudakov:Alpha AlphaEW
set VtoVHSudakov:PDFmax 1.9
cp SudakovCommon GtoGGSudakov
newdef GtoGGSudakov:SplittingFunction GtoGGSplitFn
newdef GtoGGSudakov:PDFmax 2.0
cp SudakovCommon GtoGGDARKSudakov
newdef GtoGGDARKSudakov:SplittingFunction GtoGGDARKSplitFn
newdef GtoGGDARKSudakov:PDFmax 2.0
cp SudakovCommon WtoWGammaSudakov
newdef WtoWGammaSudakov:SplittingFunction WtoWGammaSplitFn
set WtoWGammaSudakov:Alpha AlphaQED
cp SudakovCommon GtoQQbarSudakov
newdef GtoQQbarSudakov:SplittingFunction GtoQQbarSplitFn
newdef GtoQQbarSudakov:PDFmax 120.0
cp SudakovCommon GtoQQbarDARKSudakov
newdef GtoQQbarDARKSudakov:SplittingFunction GtoQQbarDARKSplitFn
newdef GtoQQbarDARKSudakov:PDFmax 120.0
cp SudakovCommon GammatoQQbarSudakov
newdef GammatoQQbarSudakov:SplittingFunction GammatoQQbarSplitFn
set GammatoQQbarSudakov:Alpha AlphaQED
newdef GammatoQQbarSudakov:PDFmax 10000.0
cp SudakovCommon GtobbbarSudakov
newdef GtobbbarSudakov:SplittingFunction GtoQQbarSplitFn
newdef GtobbbarSudakov:PDFmax 40000.0
cp SudakovCommon GtoccbarSudakov
newdef GtoccbarSudakov:SplittingFunction GtoQQbarSplitFn
newdef GtoccbarSudakov:PDFmax 2000.0
cp SudakovCommon QtoGQSudakov
newdef QtoGQSudakov:SplittingFunction QtoGQSplitFn
cp SudakovCommon QtoGammaQSudakov
newdef QtoGammaQSudakov:SplittingFunction QtoGammaQSplitFn
set QtoGammaQSudakov:Alpha AlphaQED
newdef QtoGammaQSudakov:PDFmax 2000.0
cp SudakovCommon utoGuSudakov
newdef utoGuSudakov:SplittingFunction QtoGQSplitFn
newdef utoGuSudakov:PDFFactor OverOneMinusZ
newdef utoGuSudakov:PDFmax 5.0
cp SudakovCommon dtoGdSudakov
newdef dtoGdSudakov:SplittingFunction QtoGQSplitFn
newdef dtoGdSudakov:PDFFactor OverOneMinusZ
#Use a different Sudakov for FS QCD splittings in order to use a different value of alphaS
cp QtoQGSudakov QtoQGSudakovFSR
cp GtoGGSudakov GtoGGSudakovFSR
cp GtoQQbarSudakov GtoQQbarSudakovFSR
cp GtobbbarSudakov GtobbbarSudakovFSR
cp GtobbbarSudakov GtoccbarSudakovFSR
set QtoQGSudakovFSR:Alpha AlphaQCDFSR
set GtoGGSudakovFSR:Alpha AlphaQCDFSR
set GtoQQbarSudakovFSR:Alpha AlphaQCDFSR
set GtobbbarSudakovFSR:Alpha AlphaQCDFSR
set GtoccbarSudakovFSR:Alpha AlphaQCDFSR
#
# Now add the final splittings
#
do SplittingGenerator:AddFinalSplitting u->u,g; QtoQGSudakovFSR
do SplittingGenerator:AddFinalSplitting d->d,g; QtoQGSudakovFSR
do SplittingGenerator:AddFinalSplitting s->s,g; QtoQGSudakovFSR
do SplittingGenerator:AddFinalSplitting c->c,g; QtoQGSudakovFSR
do SplittingGenerator:AddFinalSplitting b->b,g; QtoQGSudakovFSR
do SplittingGenerator:AddFinalSplitting t->t,g; QtoQGSudakovFSR
#
do SplittingGenerator:AddFinalSplitting g->g,g; GtoGGSudakovFSR
#
do SplittingGenerator:AddFinalSplitting g->u,ubar; GtoQQbarSudakovFSR
do SplittingGenerator:AddFinalSplitting g->d,dbar; GtoQQbarSudakovFSR
do SplittingGenerator:AddFinalSplitting g->s,sbar; GtoQQbarSudakovFSR
do SplittingGenerator:AddFinalSplitting g->c,cbar; GtoccbarSudakovFSR
do SplittingGenerator:AddFinalSplitting g->b,bbar; GtobbbarSudakovFSR
do SplittingGenerator:AddFinalSplitting g->t,tbar; GtoQQbarSudakovFSR
#
do SplittingGenerator:AddFinalSplitting gamma->u,ubar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->d,dbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->s,sbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->c,cbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->b,bbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->t,tbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->e-,e+; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->mu-,mu+; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->tau-,tau+; GammatoQQbarSudakov
#
do SplittingGenerator:AddFinalSplitting u->u,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting d->d,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting s->s,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting c->c,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting b->b,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting t->t,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting e-->e-,gamma; LtoLGammaSudakov
do SplittingGenerator:AddFinalSplitting mu-->mu-,gamma; LtoLGammaSudakov
do SplittingGenerator:AddFinalSplitting tau-->tau-,gamma; LtoLGammaSudakov
do SplittingGenerator:AddFinalSplitting W+->W+,gamma; WtoWGammaSudakov
#
# Now lets add the initial splittings. Remember the form a->b,c; means
# that the current particle b is given and we backward branch to new
# particle a which is initial state and new particle c which is final state
#
do SplittingGenerator:AddInitialSplitting u->u,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting d->d,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting s->s,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting c->c,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting b->b,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting u->u,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting d->d,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting s->s,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting c->c,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting b->b,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting t->t,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting g->g,g; GtoGGSudakov
#
do SplittingGenerator:AddInitialSplitting g->d,dbar; GtoQQbarSudakov
do SplittingGenerator:AddInitialSplitting g->u,ubar; GtoQQbarSudakov
do SplittingGenerator:AddInitialSplitting g->s,sbar; GtoQQbarSudakov
do SplittingGenerator:AddInitialSplitting g->c,cbar; GtoccbarSudakov
do SplittingGenerator:AddInitialSplitting g->b,bbar; GtobbbarSudakov
#
#do SplittingGenerator:AddInitialSplitting gamma->d,dbar; GammatoQQbarSudakov
#do SplittingGenerator:AddInitialSplitting gamma->u,ubar; GammatoQQbarSudakov
#do SplittingGenerator:AddInitialSplitting gamma->s,sbar; GammatoQQbarSudakov
#do SplittingGenerator:AddInitialSplitting gamma->c,cbar; GammatoQQbarSudakov
#do SplittingGenerator:AddInitialSplitting gamma->b,bbar; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->e-,e+; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->mu-,mu+; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->tau-,tau+; GammatoQQbarSudakov
#
do SplittingGenerator:AddInitialSplitting d->g,d; dtoGdSudakov
do SplittingGenerator:AddInitialSplitting u->g,u; utoGuSudakov
do SplittingGenerator:AddInitialSplitting s->g,s; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting c->g,c; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting b->g,b; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting dbar->g,dbar; dtoGdSudakov
do SplittingGenerator:AddInitialSplitting ubar->g,ubar; utoGuSudakov
do SplittingGenerator:AddInitialSplitting sbar->g,sbar; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting cbar->g,cbar; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting bbar->g,bbar; QtoGQSudakov
#
do SplittingGenerator:AddInitialSplitting d->gamma,d; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting u->gamma,u; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting s->gamma,s; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting c->gamma,c; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting b->gamma,b; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting dbar->gamma,dbar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting ubar->gamma,ubar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting sbar->gamma,sbar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting cbar->gamma,cbar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting bbar->gamma,bbar; QtoGammaQSudakov
#
# Electroweak
#
do SplittingGenerator:AddFinalSplitting u->u,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting d->d,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting s->s,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting c->c,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting b->b,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting t->t,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting u->u,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting d->d,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting s->s,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting c->c,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting b->b,Z0; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting t->t,Z0; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting u->d,W+; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting c->s,W+; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting d->u,W-; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting s->c,W-; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting u->d,W+; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting c->s,W+; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting d->u,W-; QtoQWZSudakov
do SplittingGenerator:AddInitialSplitting s->c,W-; QtoQWZSudakov
do SplittingGenerator:AddFinalSplitting c->c,h0; QtoQHSudakov
do SplittingGenerator:AddFinalSplitting b->b,h0; QtoQHSudakov
do SplittingGenerator:AddFinalSplitting t->t,h0; QtoQHSudakov
do SplittingGenerator:AddInitialSplitting c->c,h0; QtoQHSudakov
do SplittingGenerator:AddInitialSplitting b->b,h0; QtoQHSudakov
do SplittingGenerator:AddInitialSplitting t->t,h0; QtoQHSudakov
do SplittingGenerator:AddFinalSplitting gamma->W+,W-; GammatoWWSudakov
do SplittingGenerator:AddFinalSplitting Z0->W+,W-; VtoVVSudakov
do SplittingGenerator:AddFinalSplitting W+->W+,gamma; VtoVVSudakov
do SplittingGenerator:AddFinalSplitting W+->W+,Z0; VtoVVSudakov
do SplittingGenerator:AddFinalSplitting W+->W+,h0; VtoVHSudakov
do SplittingGenerator:AddFinalSplitting Z0->Z0,h0; VtoVHSudakov
#
# Electroweak l -> l V
#
#do SplittingGenerator:AddFinalSplitting e-->e-,Z0; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting mu-->mu-,Z0; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting tau-->tau-,Z0; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting nu_e->nu_e,Z0; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting nu_mu->nu_mu,Z0; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting nu_tau->nu_tau,Z0; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting e-->nu_e,W-; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting mu-->nu_mu,W-; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting tau-->nu_tau,W-; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting nu_e->e-,W+; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting nu_mu->mu-,W+; LtoLWZSudakov
#do SplittingGenerator:AddFinalSplitting nu_tau->tau-,W+; LtoLWZSudakov