Page MenuHomeHEPForge
Files F10881477

No OneTemporary
Actions

View Options

diff --git a/Models/Susy/Makefile.am b/Models/Susy/Makefile.am
--- a/Models/Susy/Makefile.am
+++ b/Models/Susy/Makefile.am
@@ -1,31 +1,32 @@
+SUBDIRS = NMSSM
if WANT_MSSM
pkglib_LTLIBRARIES = HwSusy.la
endif
HwSusy_la_SOURCES = SusyBase.cc SusyBase.h SusyBase.fh \
MSSM.cc MSSM.h MSSM.fh\
MixingMatrix.h MixingMatrix.fh MixingMatrix.cc\
SSCFSVertex.cc SSCFSVertex.h \
SSGFSVertex.cc SSGFSVertex.h \
SSHSFSFVertex.cc SSHSFSFVertex.h \
SSNFSVertex.cc SSNFSVertex.h \
SSWSSVertex.cc SSWSSVertex.h \
SSGSSVertex.cc SSGSSVertex.h \
SSGGSQSQVertex.cc SSGGSQSQVertex.h \
SSGSGSGVertex.cc SSGSGSGVertex.h \
SSNNZVertex.cc SSNNZVertex.h \
SSCCZVertex.cc SSCCZVertex.h \
SSCNWVertex.cc SSCNWVertex.h \
SSFFHVertex.cc SSFFHVertex.h \
SSGOGOHVertex.cc SSGOGOHVertex.h \
SSWWHVertex.cc SSWWHVertex.h \
SSWHHVertex.cc SSWHHVertex.h \
SSHHHVertex.cc SSHHHVertex.h \
SSHGGVertex.cc SSHGGVertex.h \
SSHPPVertex.cc SSHPPVertex.h \
SSNNPVertex.h SSNNPVertex.cc \
SSGNGVertex.h SSGNGVertex.cc \
SSNCTVertex.h SSNCTVertex.cc \
SSGVNHVertex.h SSGVNHVertex.cc\
SSGVNVVertex.h SSGVNVVertex.cc\
SSGVFSVertex.h SSGVFSVertex.cc
HwSusy_la_LDFLAGS = -module -version-info 6:0:0
diff --git a/Models/Susy/NMSSM/Makefile.am b/Models/Susy/NMSSM/Makefile.am
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/Makefile.am
@@ -0,0 +1,14 @@
+if WANT_NMSSM
+pkglib_LTLIBRARIES = HwNMSSM.la
+endif
+HwNMSSM_la_SOURCES = \
+NMSSM.cc NMSSM.h NMSSM.fh \
+NMSSMFFHVertex.h NMSSMFFHVertex.cc \
+NMSSMWWHVertex.h NMSSMWWHVertex.cc \
+NMSSMWHHVertex.h NMSSMWHHVertex.cc \
+NMSSMHSFSFVertex.h NMSSMHSFSFVertex.cc\
+NMSSMGOGOHVertex.h NMSSMGOGOHVertex.cc \
+NMSSMHHHVertex.h NMSSMHHHVertex.cc \
+NMSSMGGHVertex.h NMSSMGGHVertex.cc \
+NMSSMPPHVertex.h NMSSMPPHVertex.cc
+HwNMSSM_la_LDFLAGS = -module
diff --git a/Models/Susy/NMSSM/NMSSM.cc b/Models/Susy/NMSSM/NMSSM.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSM.cc
@@ -0,0 +1,141 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the NMSSM class.
+//
+
+#include "NMSSM.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/Interface/Reference.h"
+
+using namespace Herwig;
+
+void NMSSM::persistentOutput(PersistentOStream & os) const {
+ os << theHiggsAMix << _lambda << _kappa << ounit(_theAlambda,GeV)
+ << ounit(_theAkappa, GeV) << ounit(_lambdaVEV, GeV)
+ << ounit(_MQ3, GeV) << ounit(_MU2, GeV);
+}
+
+void NMSSM::persistentInput(PersistentIStream & is, int) {
+ is >> theHiggsAMix >> _lambda >> _kappa >> iunit(_theAlambda,GeV)
+ >> iunit(_theAkappa, GeV) >> iunit(_lambdaVEV, GeV)
+ >> iunit(_MQ3, GeV) >> iunit(_MU2, GeV);
+}
+
+ClassDescription<NMSSM> NMSSM::initNMSSM;
+// Definition of the static class description member.
+
+void NMSSM::Init() {
+
+ static ClassDocumentation<NMSSM> documentation
+ ("The NMSSM class is the base class for the NMSSM model");
+
+}
+
+void NMSSM::extractParameters(bool checkmodel) {
+ MSSM::extractParameters(false);
+ if(checkmodel) {
+ map<string,ParamMap>::const_iterator pit;
+ pit = parameters().find("modsel");
+ if(pit == parameters().end()) return;
+ ParamMap::const_iterator it;
+ // nmssm or mssm
+ it = pit->second.find(3);
+ int inmssm = (it != pit->second.end()) ? int(it->second) : 0;
+ if(inmssm == 0)
+ throw Exception() << "R-parity, CP and flavour conserving NMSSM model"
+ << " used but MSSM read in." << Exception::runerror;
+ // RPV
+ it = pit->second.find(4);
+ int irpv = (it != pit->second.end()) ? int(it->second) : 0;
+ if(irpv != 0) throw Exception() << "NMSSM model does not support RPV"
+ << Exception::runerror;
+ // CPV
+ it = pit->second.find(5);
+ int icpv = (it != pit->second.end()) ? int(it->second) : 0;
+ if(icpv != 0) throw Exception() << "NMSSM model does not support CPV"
+ << Exception::runerror;
+ // flavour violation
+ it = pit->second.find(6);
+ int ifv = (it != pit->second.end()) ? int(it->second) : 0;
+ if(ifv != 0) throw Exception() << "NMSSM model does not support "
+ << "flavour violation"
+ << Exception::runerror;
+ }
+ // get the NMSSM parameters
+ map<string,ParamMap>::const_iterator pit;
+ pit=parameters().find("msoft");
+ if( pit != parameters().end() ) {
+ ParamMap::const_iterator it;
+ it = pit->second.find(43);
+ if(it != pit->second.end()) _MQ3 = it->second*GeV;
+ it = pit->second.find(46);
+ if(it != pit->second.end()) _MU2 = it->second*GeV;
+ }
+ pit=parameters().find("nmssmrun");
+ if( pit != parameters().end() ) {
+ ParamMap::const_iterator it = pit->second.find(1);
+ if(it != pit->second.end()) _lambda = it->second;
+ it = pit->second.find(2);
+ if(it != pit->second.end()) _kappa = it->second;
+ it = pit->second.find(3);
+ if(it != pit->second.end()) _theAlambda = it->second*GeV;
+ it = pit->second.find(4);
+ if(it != pit->second.end()) _theAkappa = it->second*GeV;
+ it = pit->second.find(5);
+ if(it != pit->second.end()) _lambdaVEV = it->second*GeV;
+ }
+ pit=parameters().find("extpar");
+ if( pit != parameters().end() ) {
+ ParamMap::const_iterator it = pit->second.find(61);
+ if(_lambda==ZERO && it != pit->second.end()) _lambda = it->second;
+ it = pit->second.find(62);
+ if(_kappa==ZERO && it != pit->second.end()) _kappa = it->second;
+ it = pit->second.find(63);
+ if(_theAlambda==ZERO && it != pit->second.end()) _theAlambda = it->second*GeV;
+ it = pit->second.find(64);
+ if(_theAkappa==ZERO && it != pit->second.end()) _theAkappa = it->second*GeV;
+ it = pit->second.find(65);
+ if(_lambdaVEV==ZERO && it != pit->second.end()) _lambdaVEV = it->second*GeV;
+ it = pit->second.find(43);
+ if(_MQ3==ZERO && it != pit->second.end()) _MQ3 = it->second*GeV;
+ it = pit->second.find(46);
+ if(_MU2==ZERO && it != pit->second.end()) _MU2 = it->second*GeV;
+ }
+ else {
+ throw Exception() << "NMSSM::extractParameters - There was no EXTPAR block "
+ << "in the extracted parameters list. The model cannot "
+ << "be used without these." << Exception::runerror;
+ }
+
+
+
+
+ pit=parameters().find("msoft");
+ if( pit != parameters().end() ) {
+ ParamMap::const_iterator it;
+ if(_MQ3==ZERO) {
+ it = pit->second.find(43);
+ if(it != pit->second.end()) _MQ3 = it->second*GeV;
+ }
+ if(_MU2==ZERO) {
+ it = pit->second.find(46);
+ if(it != pit->second.end()) _MU2 = it->second*GeV;
+ }
+ }
+}
+
+void NMSSM::createMixingMatrices() {
+ map<string,pair<MatrixSize, MixingVector> >::const_iterator it;
+ for(it=mixings().begin();it!=mixings().end();++it) {
+ string name=it->first;
+ // pseudo-scalar higgs mixing
+ if (name == "nmamix") {
+ createMixingMatrix(theHiggsAMix,name,it->second.second,it->second.first);
+ }
+ }
+ // base class for neutralinos and charginos
+ MSSM::createMixingMatrices();
+}
diff --git a/Models/Susy/NMSSM/NMSSM.fh b/Models/Susy/NMSSM/NMSSM.fh
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSM.fh
@@ -0,0 +1,22 @@
+// -*- C++ -*-
+//
+// This is the forward declaration of the NMSSM class.
+//
+#ifndef HERWIG_NMSSM_FH
+#define HERWIG_NMSSM_FH
+
+#include "ThePEG/Config/Pointers.h"
+
+namespace Herwig {
+
+class NMSSM;
+
+}
+
+namespace ThePEG {
+
+ThePEG_DECLARE_POINTERS(Herwig::NMSSM,NMSSMPtr);
+
+}
+
+#endif
diff --git a/Models/Susy/NMSSM/NMSSM.h b/Models/Susy/NMSSM/NMSSM.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSM.h
@@ -0,0 +1,249 @@
+// -*- C++ -*-
+#ifndef HERWIG_NMSSM_H
+#define HERWIG_NMSSM_H
+//
+// This is the declaration of the NMSSM class.
+//
+
+#include "Herwig++/Models/Susy/MSSM.h"
+#include "NMSSM.fh"
+
+namespace Herwig {
+using namespace ThePEG;
+
+/**
+ * Here is the documentation of the NMSSM class.
+ *
+ * @see \ref NMSSMInterfaces "The interfaces"
+ * defined for NMSSM.
+ */
+class NMSSM: public MSSM {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ NMSSM() : _lambda(0.), _kappa(0.), _theAlambda(0.*MeV),
+ _theAkappa(0.*MeV), _lambdaVEV(0.*MeV),
+ _MQ3(0.*MeV), _MU2(0.*MeV)
+ {}
+
+public:
+
+ /**
+ * Mixing matrix for the neutral CP-odd Higgs bosons
+ */
+ const MixingMatrixPtr & CPoddHiggsMix() const {
+ return theHiggsAMix;
+ }
+
+ /**
+ * The NMSSM couplings
+ */
+ //@{
+ /**
+ * Superpotential \f$\lambda\f$ term
+ */
+ double lambda() const {
+ return _lambda;
+ }
+
+ /**
+ * Superpotential \f$\kappa\f$ coupling
+ */
+ double kappa() const {
+ return _kappa;
+ }
+
+ /**
+ * The V.E.V of the extra singlet field scaled
+ * by \f$ lambda\f$,
+ */
+ Energy lambdaVEV() const {
+ return _lambdaVEV;
+ }
+
+ /**
+ * Soft trilinear \f$SH_2 H_1\f$ coupling
+ */
+ Energy trilinearLambda() const {
+ return _theAlambda;
+ }
+
+ /**
+ * Soft cubic \f$S\f$ coupling
+ */
+ Energy trilinearKappa() const {
+ return _theAkappa;
+ }
+ /**
+ * left 3rd generation scalar quark mass
+ */
+ Energy MQ3() const {
+ return _MQ3;
+ }
+
+ /**
+ * right scalar top mass
+ */
+ Energy MU2() const {
+ return _MU2;
+ }
+ //@}
+
+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:
+
+ /**
+ * Extract the parameters from the input blocks
+ */
+ virtual void extractParameters(bool checkModel=true);
+
+ /**
+ * Create the mixing matrices for the model
+ */
+ virtual void createMixingMatrices();
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const {return new_ptr(*this);}
+
+ /** 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 {return new_ptr(*this);}
+ //@}
+
+private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is a concrete class with persistent data.
+ */
+ static ClassDescription<NMSSM> initNMSSM;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ NMSSM & operator=(const NMSSM &);
+
+private:
+
+ /**
+ * Higgs mixing matrix
+ */
+ MixingMatrixPtr theHiggsAMix;
+
+ /**
+ * The NMSSM couplings
+ */
+ //@{
+ /**
+ * Superpotential \f$\lambda\f$ term
+ */
+ double _lambda;
+
+ /**
+ * Superpotential \f$\kappa\f$ coupling
+ */
+ double _kappa;
+
+ /**
+ * Soft trilinear \f$SH_2 H_1\f$ coupling
+ */
+ Energy _theAlambda;
+
+ /**
+ * Soft cubic \f$S\f$ coupling
+ */
+ Energy _theAkappa;
+
+ /**
+ * The V.E.V of the extra singlet field scaled
+ * by \f$ lambda\f$
+ */
+ Energy _lambdaVEV;
+ /**
+ * left 3rd generation scalar quark mass
+ */
+ Energy _MQ3;
+
+ /**
+ * right scalar top mass
+ */
+ Energy _MU2;
+ //@}
+};
+
+}
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of NMSSM. */
+template <>
+struct BaseClassTrait<Herwig::NMSSM,1> {
+ /** Typedef of the first base class of NMSSM. */
+ typedef Herwig::MSSM NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the NMSSM class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::NMSSM>
+ : public ClassTraitsBase<Herwig::NMSSM> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::NMSSM"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * NMSSM is implemented. It may also include several, space-separated,
+ * libraries if the class NMSSM 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 "HwSusy.so HwNMSSM.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_NMSSM_H */
diff --git a/Models/Susy/NMSSM/NMSSMFFHVertex.cc b/Models/Susy/NMSSM/NMSSMFFHVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMFFHVertex.cc
@@ -0,0 +1,159 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the NMSSMFFHVertex class.
+//
+
+#include "NMSSMFFHVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "NMSSM.h"
+
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+NMSSMFFHVertex::NMSSMFFHVertex() : _mw(0.*MeV), _sinb(0.), _cosb(0.),
+ _tanb(0.), _idlast(make_pair(0,0)),
+ _q2last(0.*MeV2),
+ _masslast(make_pair(0.*MeV,0*MeV)),
+ _couplast(0.) {
+ // the quarks and neutral higgs
+ int in[5]={25,35,45,36,46};
+ for(unsigned int iy=0;iy<5;++iy)
+ for(int ix=1;ix<7;++ix)
+ addToList( -ix, ix, in[iy] );
+
+ // leptons and neutral higgs
+ for(unsigned int iy=0;iy<5;++iy)
+ for(int ix=11;ix<17;ix+=2)
+ addToList( -ix, ix, in[iy] );
+
+ // the quarks and the charged higgs
+ //H-
+ for(int ix=0;ix<3;++ix)
+ addToList(2*ix+2, -2*ix-1, -37);
+
+ //H+
+ for(int ix=0;ix<3;++ix)
+ addToList(-(2*ix+2), 2*ix+1, 37);
+
+ // the leptons and the charged higgs
+ //H-
+ for(int ix=0;ix<3;++ix)
+ addToList( 2*ix+12, -2*ix-11, -37 );
+
+ //H+
+ for(int ix=0;ix<3;++ix)
+ addToList( -(2*ix+12), 2*ix+11, 37 );
+}
+
+void NMSSMFFHVertex::persistentOutput(PersistentOStream & os) const {
+ os << _mixS << _mixP << ounit(_mw,GeV)
+ << _sinb << _cosb << _tanb << _sw << _theSM;
+}
+
+void NMSSMFFHVertex::persistentInput(PersistentIStream & is, int) {
+ is >> _mixS >> _mixP >> iunit(_mw,GeV)
+ >> _sinb >> _cosb >> _tanb >> _sw >> _theSM;
+}
+
+void NMSSMFFHVertex::doinit() {
+ // cast to NMSSM model
+ tcNMSSMPtr model=dynamic_ptr_cast<tcNMSSMPtr>(generator()->standardModel());
+ if(!model)
+ throw InitException() << "Must have the NMSSM Model in NMSSMFFHVertex::doinit()"
+ << Exception::runerror;
+ _theSM = model;
+ // sin theta_W
+ double sw2=_theSM->sin2ThetaW();
+ _sw = sqrt(sw2);
+ // get the mixing matrices
+ _mixS=model->CPevenHiggsMix();
+ if(!_mixS) throw InitException() << "Mixing matrix for CP-even neutral Higgs"
+ << " bosons is not set in NMSSMFFHVertex::doinit()"
+ << Exception::runerror;
+ _mixP=model->CPoddHiggsMix();
+ if(!_mixP) throw InitException() << "Mixing matrix for CP-odd neutral Higgs"
+ << " bosons is not set in NMSSMFFHVertex::doinit()"
+ << Exception::runerror;
+ // Mass of the W boson
+ _mw=getParticleData(ParticleID::Wplus)->mass();
+ // sin and cos beta
+ _tanb = model->tanBeta();
+ double beta = atan(_tanb);
+ _sinb=sin(beta);
+ _cosb=cos(beta);
+ // order in couplings
+ orderInGem(1);
+ orderInGs(0);
+ // base class
+ FFSVertex::doinit();
+}
+
+ClassDescription<NMSSMFFHVertex> NMSSMFFHVertex::initNMSSMFFHVertex;
+// Definition of the static class description member.
+
+void NMSSMFFHVertex::Init() {
+
+ static ClassDocumentation<NMSSMFFHVertex> documentation
+ ("The NMSSMFFHVertex class implements the vertex for the couplings"
+ " of the Higgs bosons of the NMSSM to Standard Model fermions");
+
+}
+//calulate the couplings
+void NMSSMFFHVertex::setCoupling(Energy2 q2,tcPDPtr a,tcPDPtr b, tcPDPtr c) {
+ int ihiggs=c->id();
+ int id(abs(a->id()));
+ Complex output(1.);
+ // neutral Higgs
+ if(ihiggs==25||ihiggs==35||ihiggs==45||ihiggs==36||ihiggs==46) {
+ if(_idlast.first!=id||q2!=_q2last) {
+ _idlast.first=id;
+ _masslast.first = _theSM->mass(q2,a);
+ }
+ output = _masslast.first/_mw;
+ // CP-even
+ if(ihiggs==25||ihiggs==35||ihiggs==45) {
+ int iloc = (ihiggs-25)/10;
+ output *= (id%2==0) ? (*_mixS)(iloc,1)/_sinb : (*_mixS)(iloc,0)/_cosb;
+ left(1.); right(1.);
+ }
+ // CP-odd
+ else {
+ int iloc = (ihiggs-36)/10;
+ output *= (id%2==0) ? (*_mixP)(iloc,1)/_sinb : (*_mixP)(iloc,0)/_cosb;
+ left(1.); right(-1.);
+ output *= Complex(0., 1.);
+ }
+ }
+ // Charged higgs
+ else if(abs(ihiggs)==37) {
+ output *= -sqrt(2.);
+ int id2=abs(b->id());
+ if(id2<id) swap(id,id2);
+ if(_idlast.first!=id||_idlast.second!=id2||q2!=_q2last) {
+ _idlast.first =id ;
+ _idlast.second=id2;
+ _masslast.first = _theSM->mass(q2,a);
+ _masslast.second = _theSM->mass(q2,b);
+ }
+ double rgt = _masslast.first *_tanb/_mw;
+ double lft = _masslast.second/_tanb/_mw;
+ if(ihiggs<0) swap(lft,rgt);
+ right(rgt);
+ left (lft);
+ }
+ else {
+ throw Exception() << "Unknown Higgs boson, PDG code = " << ihiggs
+ << "in NMSSMFFHVertex::setCoupling()"
+ << Exception::runerror;
+ }
+ // prefactor
+ if(q2!=_q2last) {
+ _couplast = 0.5*weakCoupling(q2);
+ _q2last=q2;
+ }
+ norm(-_couplast*output);
+}
diff --git a/Models/Susy/NMSSM/NMSSMFFHVertex.h b/Models/Susy/NMSSM/NMSSMFFHVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMFFHVertex.h
@@ -0,0 +1,209 @@
+// -*- C++ -*-
+#ifndef HERWIG_NMSSMFFHVertex_H
+#define HERWIG_NMSSMFFHVertex_H
+//
+// This is the declaration of the NMSSMFFHVertex class.
+//
+
+#include "ThePEG/Helicity/Vertex/Scalar/FFSVertex.h"
+#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/Susy/MixingMatrix.h"
+
+namespace Herwig {
+using namespace ThePEG;
+
+/** \ingroup Helicity
+ *
+ * The NMSSMFFHVertex class implements the interactions of the NMSSM Higgs bosons
+ * with the Standard Model fermions.
+ *
+ * @see \ref NMSSMFFHVertexInterfaces "The interfaces"
+ * defined for NMSSMFFHVertex.
+ */
+class NMSSMFFHVertex: public FFSVertex {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ NMSSMFFHVertex();
+
+ /** @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();
+
+ /**
+ * Calculate the couplings.
+ * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
+ * @param part1 The ParticleData pointer for the first particle.
+ * @param part2 The ParticleData pointer for the second particle.
+ * @param part3 The ParticleData pointer for the third particle.
+ * @param ioff Which particle is off-shell
+ */
+ virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,tcPDPtr part3);
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const {return new_ptr(*this);}
+
+ /** 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 {return new_ptr(*this);}
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is a concrete class with persistent data.
+ */
+ static ClassDescription<NMSSMFFHVertex> initNMSSMFFHVertex;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ NMSSMFFHVertex & operator=(const NMSSMFFHVertex &);
+
+private:
+
+ /**
+ * Mixing matrix for the CP-even Higgs bosons
+ */
+ MixingMatrixPtr _mixS;
+
+ /**
+ * Mixing matrix for the CP-odd Higgs bosons
+ */
+ MixingMatrixPtr _mixP;
+
+ /**
+ * Pointer to the SM object.
+ */
+ tcHwSMPtr _theSM;
+
+ /**
+ * Mass of the \f$W\f$ boson
+ */
+ Energy _mw;
+
+ /**
+ * \f$\sin\beta\f$
+ */
+ double _sinb;
+
+ /**
+ * \f$\cos\beta\f$
+ */
+ double _cosb;
+
+ /**
+ * \f$\tan\beta\f$
+ */
+ double _tanb;
+
+ /**
+ * \f$\sin\theta_W\f$
+ */
+ double _sw;
+
+ /**
+ * The PDG code of the last fermion the coupling was evaluated for.
+ */
+ pair<int,int> _idlast;
+
+ /**
+ * The last \f$q^2\f$ the coupling was evaluated at.
+ */
+ Energy2 _q2last;
+
+ /**
+ * The mass of the last fermion for which the coupling was evaluated.
+ */
+ pair<Energy,Energy> _masslast;
+
+ /**
+ * The last value of the coupling
+ */
+ double _couplast;
+
+};
+}
+
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of NMSSMFFHVertex. */
+template <>
+struct BaseClassTrait<Herwig::NMSSMFFHVertex,1> {
+ /** Typedef of the first base class of NMSSMFFHVertex. */
+ typedef ThePEG::Helicity::FFSVertex NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the NMSSMFFHVertex class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::NMSSMFFHVertex>
+ : public ClassTraitsBase<Herwig::NMSSMFFHVertex> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::NMSSMFFHVertex"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * NMSSMFFHVertex is implemented. It may also include several, space-separated,
+ * libraries if the class NMSSMFFHVertex 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 "HwSusy.so HwNMSSM.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_NMSSMFFHVertex_H */
diff --git a/Models/Susy/NMSSM/NMSSMGGHVertex.cc b/Models/Susy/NMSSM/NMSSMGGHVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMGGHVertex.cc
@@ -0,0 +1,203 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the NMSSMGGHVertex class.
+//
+
+#include "NMSSMGGHVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "Herwig++/Models/Susy/NMSSM/NMSSM.h"
+
+using namespace Herwig;
+
+NMSSMGGHVertex::NMSSMGGHVertex() : _sw(0.), _cw(0.), _mw(0.*MeV),
+ _mz(0.*MeV),_lambdaVEV(0.*MeV), _lambda(0.), _v1(0.*MeV),
+ _v2(0.*MeV), _triTp(0.*MeV), _triBt(0.*MeV),
+ _sb(0.), _cb(0.), _masslast(make_pair(0.*MeV,0.*MeV)),
+ _q2last(0.*MeV2), _couplast(0.), _coup(0.),
+ _hlast(0), _recalc(true) {
+ addToList(21,21,25);
+ addToList(21,21,35);
+ addToList(21,21,36);
+ addToList(21,21,45);
+ addToList(21,21,46);
+}
+
+void NMSSMGGHVertex::doinit() {
+ _theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
+ if( !_theSM ) {
+ throw InitException() << "NMSSMGGHVertex::doinit - The SM pointer is null!"
+ << Exception::abortnow;
+ }
+ // SM parameters
+ _sw = sqrt(sin2ThetaW());
+ _cw = sqrt(1. - sin2ThetaW());
+ _mw = getParticleData(24)->mass();
+ _mz = getParticleData(23)->mass();
+ _top = getParticleData(6);
+ _bt = getParticleData(5);
+
+ //NMSSM parameters
+ tcNMSSMPtr nmssm = dynamic_ptr_cast<tcNMSSMPtr>(_theSM);
+ _mixS = nmssm->CPevenHiggsMix();
+ _mixP = nmssm->CPoddHiggsMix();
+ _mixQt = nmssm->stopMix();
+ _mixQb = nmssm->sbottomMix();
+
+ double beta = atan(nmssm->tanBeta());
+ _sb = sin(beta);
+ _cb = cos(beta);
+
+ _v1 = sqrt(2.)*_mw*_cb;
+ _v2 = sqrt(2.)*_mw*_sb;
+
+ _lambda = nmssm->lambda();
+ _lambdaVEV = nmssm->lambdaVEV();
+
+ _triTp = nmssm->topTrilinear();
+ _triBt = nmssm->bottomTrilinear();
+
+ // resize vectors here and use setNParticles method
+ // to the set the actual number in the loop.
+ // Also only the top mass hass to be calculated at runtime
+ masses.resize(6, Energy());
+ masses[0] = getParticleData(6)->mass();
+ masses[1] = getParticleData(5)->mass();
+
+ masses[2] = getParticleData(1000005)->mass();
+ masses[3] = getParticleData(2000005)->mass();
+
+ masses[4] = getParticleData(1000006)->mass();
+ masses[5] = getParticleData(2000006)->mass();
+
+ type.resize(6, PDT::Spin0);
+ type[0] = PDT::Spin1Half;
+ type[1] = PDT::Spin1Half;
+ couplings.resize(6);
+
+ orderInGem(1);
+ orderInGs(2);
+
+ VVSLoopVertex::doinit();
+}
+
+void NMSSMGGHVertex::persistentOutput(PersistentOStream & os) const {
+ os << _theSM << _sw << _cw << ounit(_mw, GeV) << ounit(_mz, GeV)
+ << ounit(_lambdaVEV,GeV) << _lambda << ounit(_v1,GeV) << ounit(_v2,GeV)
+ << ounit(_triTp,GeV) << ounit(_triBt,GeV)
+ << _top << _bt << _mixS << _mixP << _mixQt << _mixQb << _sb << _cb;
+}
+
+
+void NMSSMGGHVertex::persistentInput(PersistentIStream & is, int) {
+ is >> _theSM >> _sw >> _cw >> iunit(_mw, GeV) >> iunit(_mz, GeV)
+ >> iunit(_lambdaVEV,GeV) >> _lambda >> iunit(_v1,GeV) >> iunit(_v2,GeV)
+ >> iunit(_triTp,GeV) >> iunit(_triBt,GeV)
+ >> _top >> _bt >> _mixS >> _mixP >> _mixQt >> _mixQb >> _sb >> _cb;
+}
+
+ClassDescription<NMSSMGGHVertex> NMSSMGGHVertex::initNMSSMGGHVertex;
+// Definition of the static class description member.
+
+void NMSSMGGHVertex::Init() {
+
+ static ClassDocumentation<NMSSMGGHVertex> documentation
+ ("The effective coupling of a higgs to a pair of gluons in the "
+ "NMSSM.");
+
+}
+
+void NMSSMGGHVertex::setCoupling(Energy2 q2, tcPDPtr p1, tcPDPtr p2,
+ tcPDPtr p3) {
+ long hid(p3->id());
+ if( q2 != _q2last ) {
+ _couplast = sqr(strongCoupling(q2));
+ _coup = weakCoupling(q2);
+ _q2last = q2;
+ _recalc = true;
+ }
+ norm(_couplast*_coup);
+ // scalar higgs bosons
+ if( hid != _hlast ) {
+ _hlast = hid;
+ _recalc = true;
+ if( hid % 5 == 0 ) {
+ // location of the higgs
+ int iloc = (hid - 25)/10;
+ // 6 particles in the loop
+ setNParticles(6);
+ // top and bottom quark masses
+ Energy mt = _theSM->mass(q2, _top);
+ Energy mb = _theSM->mass(q2, _bt);
+ Complex c(0.);
+ // couplings for the top quark loop
+ c = -0.25*mt*(*_mixS)(iloc, 1)/_sb/_mw;
+
+ couplings[0] = make_pair(c,c);
+ masses[0] = mt;
+ // couplings for the bottom quark loop
+ c = -0.25*mb*(*_mixS)(iloc, 0)/_cb/_mw;
+
+ couplings[1] = make_pair(c,c);
+ masses[1] = mb;
+ // sbottoms
+ double f1 = mb/_mw/_cb;
+ complex<Energy> f2 = 0.5*_mz/_cw*
+ ( - _cb*(*_mixS)(iloc,0) + _sb*(*_mixS)(iloc,1));
+ complex<Energy> cpl;
+ for(unsigned int ix=0;ix<2;++ix) {
+ cpl = -f2*( (1. - 2.*sqr(_sw)/3.)*(*_mixQb)(ix, 0)*(*_mixQb)(ix, 0)
+ + 2.*sqr(_sw)*(*_mixQb)(ix, 1)*(*_mixQb)(ix, 1)/3.)
+ - f1*mb*(*_mixS)(iloc,0)
+ *((*_mixQb)(ix, 0)*(*_mixQb)(ix, 0) + (*_mixQb)(ix, 1)*(*_mixQb)(ix, 1))
+ - 0.5*f1*(-_lambdaVEV*(*_mixS)(iloc,1) - _lambda*_v2*(*_mixS)(iloc,2)/_coup
+ + _triBt*(*_mixS)(iloc,0))*((*_mixQb)(ix, 1)*(*_mixQb)(ix, 0)
+ + (*_mixQb)(ix, 0)*(*_mixQb)(ix, 1));
+
+ couplings[2+ix] = make_pair(0.5*cpl*UnitRemoval::InvE,0.5*cpl*UnitRemoval::InvE);
+ }
+ // stop
+ f1 = mt/_mw/_sb;
+ for(unsigned int ix=0;ix<2;++ix) {
+ cpl =+f2*( (1. - 4.*sqr(_sw)/3.)*(*_mixQt)(ix, 0)*(*_mixQt)(ix, 0)
+ + 4.*sqr(_sw)*(*_mixQt)(ix, 1)*(*_mixQt)(ix, 1)/3.)
+ - f1*mt*(*_mixS)(iloc,1)
+ *((*_mixQt)(ix, 0)*(*_mixQt)(ix, 0)
+ + (*_mixQt)(ix, 1)*(*_mixQt)(ix, 1))
+ - 0.5*f1*(-_lambdaVEV*(*_mixS)(iloc,0) - _lambda*_v1*(*_mixS)(iloc,2)/_coup
+ + _triTp*(*_mixS)(iloc,1))*((*_mixQt)(ix, 1)*(*_mixQt)(ix, 0)
+ + (*_mixQt)(ix, 0)*(*_mixQt)(ix, 1));
+
+ couplings[4+ix] = make_pair(0.5*cpl*UnitRemoval::InvE,0.5*cpl*UnitRemoval::InvE);
+ }
+ }
+ // pseudoscalar higgs bosons
+ else {
+ // location of the higgs
+ int iloc = (hid - 36)/10;
+ // 2 particles in the loop
+ setNParticles(2);
+ // top and bottom quark masses
+ Energy mt = _theSM->mass(q2, _top);
+ Energy mb = _theSM->mass(q2, _bt);
+ Complex c(0.);
+ // top quark couplings
+ c = Complex(0.,-1.)*0.25*mt*(*_mixP)(iloc, 1)/_sb/_mw;
+
+ couplings[0] = make_pair(-c,c);
+ masses[0] = mt;
+ // bottom quark couplings
+ c = Complex(0., -1.)*0.25*mb*(*_mixP)(iloc, 0)/_cb/_mw;
+
+ couplings[1] = make_pair(-c,c);
+ masses[1] = mb;
+ }
+ }
+
+ if( _recalc ) {
+ VVSLoopVertex::setCoupling(q2, p1, p2, p3);
+ _recalc = false;
+ }
+}
diff --git a/Models/Susy/NMSSM/NMSSMGGHVertex.h b/Models/Susy/NMSSM/NMSSMGGHVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMGGHVertex.h
@@ -0,0 +1,289 @@
+// -*- C++ -*-
+#ifndef HERWIG_NMSSMGGHVertex_H
+#define HERWIG_NMSSMGGHVertex_H
+//
+// This is the declaration of the NMSSMGGHVertex class.
+//
+
+#include "Herwig++/Models/General/VVSLoopVertex.h"
+#include "Herwig++/Models/StandardModel/StandardModel.fh"
+#include "Herwig++/Models/Susy/MixingMatrix.h"
+
+namespace Herwig {
+using namespace ThePEG;
+
+/**
+ * This class implements the effective vertex for a higgs coupling
+ * to a pair of gluons in the NMSSM.
+ *
+ * @see \ref NMSSMGGHVertexInterfaces "The interfaces"
+ * defined for NMSSMGGHVertex.
+ */
+class NMSSMGGHVertex: public VVSLoopVertex {
+
+public:
+
+ /** @name Standard constructors and destructors. */
+ //@{
+ /**
+ * The default constructor.
+ */
+ NMSSMGGHVertex();
+ //@}
+
+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();
+
+ /**
+ * Calculate couplings
+ *@param q2 Scale at which to evaluate coupling
+ *@param p1 ParticleData pointer to first particle
+ *@param p2 ParticleData pointer to second particle
+ *@param p3 ParticleData pointer to third particle
+ */
+ virtual void setCoupling(Energy2 q2, tcPDPtr p1, tcPDPtr p2,
+ tcPDPtr p3);
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const {return new_ptr(*this);}
+
+ /** 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 {return new_ptr(*this);}
+ //@}
+
+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:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is a concrete class with persistent data.
+ */
+ static ClassDescription<NMSSMGGHVertex> initNMSSMGGHVertex;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ NMSSMGGHVertex & operator=(const NMSSMGGHVertex &);
+
+private:
+
+ /** @name Stored parameters. */
+ //@{
+ /**
+ * The SM pointer
+ */
+ tcHwSMPtr _theSM;
+
+ /**
+ * \f$ \sin\theta_W\f$
+ */
+ double _sw;
+
+ /**
+ * \f$ \cos\theta_W\f$
+ */
+ double _cw;
+
+ /**
+ * \f$ M_W\f$
+ */
+ Energy _mw;
+
+ /**
+ * \f$ M_Z \f$
+ */
+ Energy _mz;
+
+ /**
+ * The product \f$\lambda \langle S\rangle \f$.
+ */
+
+ Energy _lambdaVEV;
+
+ /**
+ * The coefficient of the trilinear \f$SH_2 H_1\f$ term in the superpotential
+ */
+
+ double _lambda;
+
+ /**
+ * The value of the VEV of the higgs that couples to the up-type sector
+ * \f$ g*sqrt(2)M_W\cos\beta \f$
+ */
+
+ Energy _v1;
+
+ /**
+ * The value of the VEV of the higgs that couples to the down-type sector
+ * \f$ g*sqrt(2)M_W\cos\beta \f$
+ */
+ Energy _v2;
+
+ /**
+ * The top quark trilinear coupling
+ */
+ complex<Energy> _triTp;
+
+ /**
+ * The bottom quark trilinear coupling
+ */
+ complex<Energy> _triBt;
+
+ /**
+ * A pointer to the top quark
+ */
+ tcPDPtr _top;
+
+ /**
+ * A pointer to the bottom quark
+ */
+ tcPDPtr _bt;
+
+ /**
+ * CP-even Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixS;
+
+ /**
+ * CP-even Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixP;
+
+ /**
+ * \f$\tilde{t}\f$ mixing matrix
+ */
+ MixingMatrixPtr _mixQt;
+
+ /**
+ * \f$\tilde{b}\f$ mixing matrix
+ */
+ MixingMatrixPtr _mixQb;
+
+ /**
+ * \f$ \sin\beta\f$
+ */
+ double _sb;
+
+ /**
+ * \f$ \cos\beta\f$
+ */
+ double _cb;
+
+ /**
+ * The top and bottom quark masses calculated at the last value
+ * of \f$q^2\f$
+ */
+ pair<Energy, Energy> _masslast;
+
+ /**
+ * The scale at which the coupling was last evaluated
+ */
+ Energy2 _q2last;
+
+ /**
+ * The value of the overall normalisation when the coupling was last
+ * evaluated.
+ */
+ double _couplast;
+
+ /**
+ * The value of the weak coupling was last
+ * evaluated.
+ */
+ double _coup;
+
+ /**
+ * The PDG code of the Higgs particle when the vertex was last evaluated
+ */
+ long _hlast;
+
+ /**
+ * Whether the tensor coefficient need recalculating or not
+ */
+ bool _recalc;
+ //@}
+};
+
+}
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of NMSSMGGHVertex. */
+template <>
+struct BaseClassTrait<Herwig::NMSSMGGHVertex,1> {
+ /** Typedef of the first base class of NMSSMGGHVertex. */
+ typedef Herwig::VVSLoopVertex NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the NMSSMGGHVertex class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::NMSSMGGHVertex>
+ : public ClassTraitsBase<Herwig::NMSSMGGHVertex> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::NMSSMGGHVertex"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * NMSSMGGHVertex is implemented. It may also include several, space-separated,
+ * libraries if the class NMSSMGGHVertex 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 "HwSusy.so HwNMSSM.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_NMSSMGGHVertex_H */
diff --git a/Models/Susy/NMSSM/NMSSMGOGOHVertex.cc b/Models/Susy/NMSSM/NMSSMGOGOHVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMGOGOHVertex.cc
@@ -0,0 +1,290 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the NMSSMGOGOHVertex class.
+//
+
+#include "NMSSMGOGOHVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "NMSSM.h"
+
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+NMSSMGOGOHVertex::NMSSMGOGOHVertex() : _lambda(0.), _kappa(0.), _sinb(0.),
+ _cosb(0.), _sw(0.), _cw(0.),
+ _q2last(0.*MeV2), _couplast(0.) {
+ int ieven[3]={25,35,45};
+ int iodd [2]={36,46};
+ long ichar[2]={1000024,1000037};
+ long ineut[5]={1000022,1000023,1000025,1000035,1000045};
+ // CP-even charginos
+ for(unsigned int ix=0;ix<2;++ix) {
+ for(unsigned int iy=0;iy<2;++iy) {
+ for(unsigned int iz=0;iz<3;++iz) {
+ addToList(-ichar[ix], ichar[iy], ieven[iz]);
+ }
+ }
+ }
+ // CP-odd charginos
+ for(unsigned int ix=0;ix<2;++ix) {
+ for(unsigned int iy=0;iy<2;++iy) {
+ for(unsigned int iz=0;iz<2;++iz) {
+ addToList(-ichar[ix], ichar[iy], iodd [iz]);
+
+ }
+ }
+ }
+ // CP-even neutralinos
+ for(unsigned int ix=0;ix<5;++ix) {
+ for(unsigned int iy=0;iy<5;++iy) {
+ for(unsigned int iz=0;iz<3;++iz) {
+ addToList( ineut[ix], ineut[iy], ieven[iz]);
+ }
+ }
+ }
+ // CP-odd neutralinos
+ for(unsigned int ix=0;ix<5;++ix) {
+ for(unsigned int iy=0;iy<5;++iy) {
+ for(unsigned int iz=0;iz<2;++iz) {
+ addToList( ineut[ix], ineut[iy], iodd[iz]);
+ }
+ }
+ }
+
+ // charged higgs
+ for(unsigned int ix=0;ix<5;++ix) {
+ for(unsigned int iy=0;iy<2;++iy) {
+ addToList(ineut[ix], -ichar[iy], 37);
+
+ addToList(ineut[ix], ichar[iy], -37);
+
+ }
+ }
+}
+
+void NMSSMGOGOHVertex::persistentOutput(PersistentOStream & os) const {
+ os << _mixV << _mixU << _mixN << _mixS << _mixP << _lambda << _kappa << _sinb
+ << _cosb << _sw << _cw;
+}
+
+void NMSSMGOGOHVertex::persistentInput(PersistentIStream & is, int) {
+ is >> _mixV >> _mixU >> _mixN >> _mixS >> _mixP >> _lambda >> _kappa >> _sinb
+ >> _cosb >> _sw >> _cw;
+}
+
+void NMSSMGOGOHVertex::doinit() {
+
+ tcNMSSMPtr model=dynamic_ptr_cast<tcNMSSMPtr>(generator()->standardModel());
+ // SM parameters
+ // sin theta_W
+
+ double sw2=sin2ThetaW();
+ _cw=sqrt(1.0 - sw2);
+ _sw=sqrt(sw2);
+ if(!model)
+ throw InitException() << "Must have the NMSSM Model in "
+ << "NMSSMGOGOHVertex::doinit()"
+ << Exception::runerror;
+ // get the mixing matrices
+ // higgs
+ _mixS=model->CPevenHiggsMix();
+ if(!_mixS)
+ throw InitException() << "Mixing matrix for CP-even neutral Higgs"
+ << " bosons is not set in NMSSMGOGOHVertex::doinit()"
+ << Exception::runerror;
+ _mixP=model->CPoddHiggsMix();
+ if(!_mixP)
+ throw InitException() << "Mixing matrix for CP-odd neutral Higgs"
+ << " bosons is not set in NMSSMGOGOHVertex::doinit()"
+ << Exception::runerror;
+ // charginos
+ _mixU = model->charginoUMix();
+ _mixV = model->charginoVMix();
+ if(!_mixU || !_mixV)
+ throw InitException() << "NMSSMGOGOHVertex::doinit - "
+ << "A mixing matrix pointer is null. U: "
+ << _mixU << " V: " << _mixV
+ << Exception::abortnow;
+ // neutralinos
+ _mixN = model->neutralinoMix();
+ if(!_mixN)
+ throw InitException() << "NMSSMGOGOHVertex::doinit - The neutralino "
+ << "mixing matrix pointer is null."
+ << Exception::abortnow;
+ // kappa and lambda couplings
+ _lambda = model->lambda();
+ _kappa = model->kappa();
+ // sin and cos beta
+ double beta = atan(model->tanBeta());
+ _sinb=sin(beta);
+ _cosb=cos(beta);
+ // order in the couplings
+ orderInGem(1);
+ orderInGs(0);
+ FFSVertex::doinit();
+}
+
+ClassDescription<NMSSMGOGOHVertex> NMSSMGOGOHVertex::initNMSSMGOGOHVertex;
+// Definition of the static class description member.
+
+void NMSSMGOGOHVertex::Init() {
+
+ static ClassDocumentation<NMSSMGOGOHVertex> documentation
+ ("The NMSSMGOGOHVertex class implements the couplings of the Higgs bosons"
+ " of the NMSSM and the electroweak gauginos");
+
+}
+
+void NMSSMGOGOHVertex::setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,
+ tcPDPtr part3) {
+ long id1(part1->id()), id2(part2->id()),
+ id3(part3->id()), ihigg(0), ig1(0), ig2(0);
+ if( abs(id1) == 25 || abs(id1) == 35 || abs(id1) == 45 ||
+ abs(id1) == 36 || abs(id1) == 46 || abs(id1) == 37 ) {
+ ihigg = id1;
+ ig1 = id2;
+ ig2 = id3;
+ }
+ else if( abs(id2) == 25 || abs(id2) == 35 ||
+ abs(id2) == 45 ||abs(id2) == 36 ||abs(id2) == 46 || abs(id2) == 37 ) {
+ ihigg = id2;
+ ig1 = id1;
+ ig2 = id3;
+ }
+ else if( abs(id3) ==25 || abs(id3) == 35 ||
+ abs(id3) == 45 ||abs(id3) == 36 ||abs(id3) == 46 || abs(id3) == 37 ) {
+ ihigg = id3;
+ ig1 = id1;
+ ig2 = id2;
+ }
+ else {
+ throw HelicityConsistencyError()
+ << "NMSSMGOGOHVertex::setCoupling - There is no higgs particle in "
+ << "this vertex. Particles: " << id1 << " " << id2 << " " << id3
+ << Exception::runerror;
+ return;
+ }
+
+ // weak coupling
+ if(q2!=_q2last) {
+ _couplast = weakCoupling(q2);
+ _q2last = q2;
+ }
+ double rt = sqrt(0.5);
+ // CP-even neutral higgs
+ if(ihigg == 25 || ihigg == 35 || ihigg == 45) {
+ int iloc = (ihigg - 25)/10;
+ // chargino
+ if(abs(ig1) == 1000024 || abs(ig1) == 1000037) {
+ if( ig1 < 0 ) swap(ig1, ig2);
+ int ic1 = abs(ig1)==1000024 ? 0 : 1;
+ int ic2 = abs(ig2)==1000024 ? 0 : 1;
+ Complex coupL = -_lambda*rt*conj((*_mixS)(iloc,2)*(*_mixU)(ic1,1)*(*_mixV)(ic2,1))
+ -_couplast*rt*(conj((*_mixS)(iloc,0)*(*_mixU)(ic1,1)*(*_mixV)(ic2,0) +
+ (*_mixS)(iloc,1)*(*_mixU)(ic1,0)*(*_mixV)(ic2,1)));
+ Complex coupR = -_lambda*rt*(*_mixS)(iloc,2)*(*_mixU)(ic2,1)*(*_mixV)(ic1,1)
+ -_couplast*rt*((*_mixS)(iloc,0)*(*_mixU)(ic2,1)*(*_mixV)(ic1,0)+
+ (*_mixS)(iloc,1)*(*_mixU)(ic2,0)*(*_mixV)(ic1,1));
+ left(coupL);
+ right(coupR);
+ norm(1.0);
+ }
+ // neutralino
+ else {
+ int in1 = (ig1 < 1000024) ? (ig1 - 1000022) : (ig1 - 1000005)/10;
+ int in2 = (ig2 < 1000024) ? (ig2 - 1000022) : (ig2 - 1000005)/10;
+ Complex us1 = (*_mixS)(iloc, 0), us2 = (*_mixS)(iloc, 1);
+ Complex us3 = (*_mixS)(iloc, 2);
+ Complex ni1 = (*_mixN)(in1,0), nj1 = (*_mixN)(in2,0);
+ Complex ni2 = (*_mixN)(in1,1), nj2 = (*_mixN)(in2,1);
+ Complex ni3 = (*_mixN)(in1,3), nj3 = (*_mixN)(in2,3);
+ Complex ni4 = (*_mixN)(in1,2), nj4 = (*_mixN)(in2,2);
+ Complex ni5 = (*_mixN)(in1,4), nj5 = (*_mixN)(in2,4);
+ Complex YL =
+ - _lambda*rt*(us2*(ni4*nj5 + ni5*nj4) +
+ us1*(ni3*nj5 + ni5*nj3) +
+ us3*(ni3*nj4 + ni4*nj3))
+ + sqrt(2.)*_kappa*us3*ni5*nj5
+ - _couplast*0.5*(us2*(ni2*nj3 + ni3*nj2) -
+ us1*(ni2*nj4 + ni4*nj2))
+ + _couplast*0.5*_sw*(us2*(ni1*nj3 + ni3*nj1) -
+ us1*(ni1*nj4 + ni4*nj1) )/_cw;
+ left(-conj(YL));
+ right(-YL);
+ norm(1.0);
+ }
+ }
+ // CP-odd neutral higgs
+ else if(ihigg==36||ihigg==46) {
+ int iloc = (ihigg-36)/10;
+ // chargino
+ if(abs(ig1)==1000024||abs(ig1)==1000037) {
+ if( ig1 < 0 ) swap(ig1, ig2);
+ int ic1 = abs(ig1)==1000024 ? 0 : 1;
+ int ic2 = abs(ig2)==1000024 ? 0 : 1;
+ Complex QL = Complex(0,-1.0)*
+ (_lambda*rt*conj((*_mixP)(iloc,2)*(*_mixU)(ic1,1)*(*_mixV)(ic2,1))
+ -_couplast*rt*conj(((*_mixP)(iloc,0)*(*_mixU)(ic1,1)*(*_mixV)(ic2,0) +
+ (*_mixP)(iloc,1)*(*_mixU)(ic1,0)*(*_mixV)(ic2,1))));
+ Complex QR = Complex(0,-1.0)*
+ (_lambda*rt*(*_mixP)(iloc,2)*(*_mixU)(ic2,1)*(*_mixV)(ic1,1)
+ -_couplast*rt*((*_mixP)(iloc,0)*(*_mixU)(ic2,1)*(*_mixV)(ic1,0) +
+ (*_mixP)(iloc,1)*(*_mixU)(ic2,0)*(*_mixV)(ic1,1)));
+ left(QL);
+ right(-QR);
+ norm(1.);
+ }
+ // neutralino
+ else {
+ int in1 = (ig1 < 1000024) ? (ig1 - 1000022) : (ig1 - 1000005)/10;
+ int in2 = (ig2 < 1000024) ? (ig2 - 1000022) : (ig2 - 1000005)/10;
+ Complex up1 = (*_mixP)(iloc, 0), up2 = (*_mixP)(iloc, 1);
+ Complex up3 = (*_mixP)(iloc, 2);
+ Complex ni1 = (*_mixN)(in1,0), nj1 = (*_mixN)(in2,0);
+ Complex ni2 = (*_mixN)(in1,1), nj2 = (*_mixN)(in2,1);
+ Complex ni3 = (*_mixN)(in1,2), nj3 = (*_mixN)(in2,2);
+ Complex ni4 = (*_mixN)(in1,3), nj4 = (*_mixN)(in2,3);
+ Complex ni5 = (*_mixN)(in1,4), nj5 = (*_mixN)(in2,4);
+ Complex AL =
+ _lambda*rt*(up2*(ni3*nj5 + ni5*nj3) +
+ up1*(ni4*nj5 + ni5*nj4) +
+ up3*(ni3*nj4 + ni4*nj3))
+ - sqrt(2.)*_kappa*up3*ni5*nj5
+ - _couplast*0.5*(up2*(ni2*nj4 + ni4*nj2) -
+ up1*(ni2*nj3 + ni3*nj2))
+ + _couplast*0.5*_sw*(up2*(ni1*nj4 + ni4*nj1) -
+ up1*(ni1*nj3 + ni3*nj1))/_cw;
+ AL *= Complex(0.0, -1.0);
+ left(conj(AL));
+ right(AL);
+ norm(1.);
+ }
+ }
+ // charged higgs
+ else {
+ if (abs(ig1) == 1000024 || abs(ig1) == 1000037) swap (ig1,ig2);
+ int in = (abs(ig1) < 1000024) ? (ig1-1000022) : (ig1-1000005)/10;
+ int ic = (abs(ig2) == 1000024) ? 0 : 1;
+ Complex QpR = _lambda*_cosb*(*_mixU)(ic,1)*(*_mixN)(in,4)
+ -_sinb*_couplast*(rt*(*_mixU)(ic,1)*(_sw*(*_mixN)(in,0)/_cw + (*_mixN)(in,1))
+ - (*_mixU)(ic,0)*(*_mixN)(in,2));
+ Complex QpL = _lambda*_sinb*(*_mixV)(ic,1)*(*_mixN)(in,4)
+ + _couplast*_cosb*(rt*(*_mixV)(ic,1)
+ *(_sw*(*_mixN)(in,0)/_cw + (*_mixN)(in,1))
+ + (*_mixV)(ic,0)*(*_mixN)(in,3));
+ QpL = conj(QpL);
+ if(ihigg > 0) {
+ left (QpL);
+ right(QpR);
+ norm(-1.);
+ }
+ else {
+ left (conj(QpR));
+ right(conj(QpL));
+ norm(-1.);
+ }
+ }
+}
diff --git a/Models/Susy/NMSSM/NMSSMGOGOHVertex.h b/Models/Susy/NMSSM/NMSSMGOGOHVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMGOGOHVertex.h
@@ -0,0 +1,221 @@
+// -*- C++ -*-
+#ifndef HERWIG_NMSSMGOGOHVertex_H
+#define HERWIG_NMSSMGOGOHVertex_H
+//
+// This is the declaration of the NMSSMGOGOHVertex class.
+//
+
+#include "ThePEG/Helicity/Vertex/Scalar/FFSVertex.h"
+#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/Susy/MixingMatrix.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * Here is the documentation of the NMSSMGOGOHVertex class.
+ *
+ * @see \ref NMSSMGOGOHVertexInterfaces "The interfaces"
+ * defined for NMSSMGOGOHVertex.
+ */
+class NMSSMGOGOHVertex: public FFSVertex {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ inline NMSSMGOGOHVertex();
+
+ /** @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();
+
+ /**
+ * Calculate the couplings. This method is virtual and must be implemented in
+ * classes inheriting from this.
+ * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
+ * @param part1 The ParticleData pointer for the first particle.
+ * @param part2 The ParticleData pointer for the second particle.
+ * @param part3 The ParticleData pointer for the third particle.
+ * @param ioff An integer referring to which particle in the list is
+ * offshell if applicable.
+ */
+ virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,
+ tcPDPtr part3);
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const {return new_ptr(*this);}
+
+ /** 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 {return new_ptr(*this);}
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is a concrete class with persistent data.
+ */
+ static ClassDescription<NMSSMGOGOHVertex> initNMSSMGOGOHVertex;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ NMSSMGOGOHVertex & operator=(const NMSSMGOGOHVertex &);
+
+private:
+
+ /**
+ * The various mixing matrices and couplings
+ */
+ //@{
+ /**
+ * The V chargino mixing matrix
+ */
+ MixingMatrixPtr _mixV;
+
+ /**
+ * The U chargino mixing matrix
+ */
+ MixingMatrixPtr _mixU;
+
+ /**
+ * The neutralino mixing matrix
+ */
+ MixingMatrixPtr _mixN;
+
+ /**
+ * The CP-even neutral Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixS;
+
+ /**
+ * The CP-odd neutral Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixP;
+
+ /**
+ * The tri-linear \f$\lambda\f$ coupling
+ */
+ double _lambda;
+
+ /**
+ * The tri-linear \f$\kappa\f$ coupling
+ */
+ double _kappa;
+
+ /**
+ * \f$\sin\beta\f$
+ */
+ double _sinb;
+
+ /**
+ * \f$\cos\beta\f$
+ */
+ double _cosb;
+
+ /**
+ * \f$\sin\theta_W\f$
+ */
+ double _sw;
+
+ /**
+ * \f$\cos\theta_W\f$
+ */
+ double _cw;
+
+ /**
+ * The last \f$q^2\f$ the coupling was evaluated at.
+ */
+ Energy2 _q2last;
+
+ /**
+ * The last value of the coupling
+ */
+ double _couplast;
+ //@}
+
+};
+}
+
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of NMSSMGOGOHVertex. */
+template <>
+struct BaseClassTrait<Herwig::NMSSMGOGOHVertex,1> {
+ /** Typedef of the first base class of NMSSMGOGOHVertex. */
+ typedef ThePEG::Helicity::FFSVertex NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the NMSSMGOGOHVertex class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::NMSSMGOGOHVertex>
+ : public ClassTraitsBase<Herwig::NMSSMGOGOHVertex> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::NMSSMGOGOHVertex"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * NMSSMGOGOHVertex is implemented. It may also include several, space-separated,
+ * libraries if the class NMSSMGOGOHVertex 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 "HwSusy.so HwNMSSM.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_NMSSMGOGOHVertex_H */
diff --git a/Models/Susy/NMSSM/NMSSMHHHVertex.cc b/Models/Susy/NMSSM/NMSSMHHHVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMHHHVertex.cc
@@ -0,0 +1,275 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the NMSSMHHHVertex class.
+//
+
+#include "NMSSMHHHVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "NMSSM.h"
+
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+NMSSMHHHVertex::NMSSMHHHVertex() : _mw(0.*MeV), _mz(0.*MeV), _sw2(0.),
+ _cw(0.), _lambda(0.), _kappa(0.) ,
+ _lambdaVEV(0.*MeV), _theAl(0.*MeV),
+ _theAk(0.*MeV), _sb(0.), _cb(0.),
+ _s2b(0.), _c2b(0.), _vu(0.*MeV),
+ _vd(0.*MeV), _s(0.*MeV), _q2last(0.*MeV2),
+ _glast(0.), _MQ3(0.*MeV), _MU2(0.*MeV),
+ _masslast(make_pair(0.*MeV,0*MeV)),
+ _includeRadiative(true) {
+ // PDG codes for the particles in vertex _vd
+ //CP-even Higgs
+ addToList(25, 35, 45);
+ for( unsigned int i = 25; i <= 45; i += 10 ) {
+ addToList(i, i, 25);
+ addToList(i, i, 35);
+ addToList(i, i, 45);
+ //Charged Higgs
+ addToList(i, 37, -37);
+ //CP-odd Higgs
+ addToList(i, 36, 36);
+ addToList(i, 36, 46);
+ addToList(i, 46, 36);
+ addToList(i, 46, 46);
+ }
+}
+
+void NMSSMHHHVertex::doinit() {
+ _theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
+ tcNMSSMPtr nmssm = dynamic_ptr_cast<tcNMSSMPtr>(_theSM);
+ if( !nmssm )
+ throw InitException() << "NMSSMHHHVertex::doinit - The model object is"
+ << "not the NMSSM object."
+ << Exception::runerror;
+ //SM parameters
+ _mw = getParticleData(24)->mass();
+ _mz = getParticleData(23)->mass();
+ _sw2 = sin2ThetaW();
+ _cw = sqrt(1. - _sw2);
+ //NMSSM parameters
+ _mixS = nmssm->CPevenHiggsMix();
+ _mixP = nmssm->CPoddHiggsMix();
+ if( !_mixS || !_mixP )
+ throw InitException() << "NMSSMHHHVertex::doinit - One of the mixing matrix "
+ << "pointers is null, cannot continue. S: "
+ << _mixS << " P: " << _mixP << Exception::runerror;
+ _lambda = nmssm->lambda();
+ _kappa = nmssm->kappa();
+ _lambdaVEV = nmssm->lambdaVEV();
+ _theAl = nmssm->trilinearLambda();
+ _theAk = nmssm->trilinearKappa();
+ _MQ3 = nmssm->MQ3();
+ _MU2 = nmssm->MU2();
+ double beta = atan(nmssm->tanBeta());
+ _sb = sin(beta);
+ _cb = cos(beta);
+ _vd = sqrt(2)*_mw*_cb;
+ _vu = sqrt(2)*_mw*_sb;
+ _s = _lambdaVEV/_lambda;
+ // order in couplings
+ orderInGem(1);
+ orderInGs(0);
+ SSSVertex::doinit();
+}
+
+void NMSSMHHHVertex::persistentOutput(PersistentOStream & os) const {
+ os << ounit(_mw, GeV) << ounit(_mz,GeV) << ounit(_mb, GeV) << ounit(_mt,GeV)
+ << _sw2 << _cw << _lambda << _includeRadiative
+ << _kappa << ounit(_lambdaVEV,GeV) << ounit(_theAl, GeV)
+ << ounit(_theAk,GeV) << _sb << _cb << _s2b << _c2b
+ << ounit(_vu,GeV) << ounit(_vu,GeV) << ounit(_s,GeV) << _mixS << _mixP
+ << ounit(_MQ3,GeV) << ounit(_MU2,GeV) << _theSM;
+}
+
+void NMSSMHHHVertex::persistentInput(PersistentIStream & is, int) {
+ is >> iunit(_mw, GeV) >> iunit(_mz,GeV)>> iunit(_mb, GeV) >> iunit(_mt,GeV)
+ >> _sw2 >> _cw >> _lambda >> _includeRadiative
+ >> _kappa >> iunit(_lambdaVEV,GeV) >> iunit(_theAl, GeV)
+ >> iunit(_theAk,GeV) >> _sb >> _cb >> _s2b >> _c2b
+ >> iunit(_vu,GeV) >> iunit(_vd,GeV) >> iunit(_s,GeV)>> _mixS >> _mixP
+ >> iunit(_MQ3,GeV) >> iunit(_MU2,GeV) >> _theSM;
+}
+
+ClassDescription<NMSSMHHHVertex> NMSSMHHHVertex::initNMSSMHHHVertex;
+// Definition of the static class description member.
+
+void NMSSMHHHVertex::Init() {
+
+ static ClassDocumentation<NMSSMHHHVertex> documentation
+ ("This is the triple Higgs coupling in the NMSSM.");
+
+ static Switch<NMSSMHHHVertex,bool> interfaceIncludeRadiativeCorrections
+ ("IncludeRadiativeCorrections",
+ "Include radiative corrections in the vertex",
+ &NMSSMHHHVertex::_includeRadiative, true, false, false);
+ static SwitchOption interfaceIncludeRadiativeCorrectionsYes
+ (interfaceIncludeRadiativeCorrections,
+ "Yes",
+ "Include the radiative terms",
+ true);
+ static SwitchOption interfaceIncludeRadiativeCorrectionsNo
+ (interfaceIncludeRadiativeCorrections,
+ "No",
+ "Don't include them",
+ false);
+
+}
+
+//calulate couplings
+void NMSSMHHHVertex::setCoupling(Energy2 q2,tcPDPtr p1,tcPDPtr p2,
+ tcPDPtr p3) {
+ using Constants::pi;
+ long higgs[3] = {p1->id(), p2->id(), p3->id()};
+ unsigned int ns(0), np(0), nc(0);
+ for( int i = 0; i < 3; ++i ) {
+ if( higgs[i] == 25 || higgs[i] == 35 || higgs[i] == 45 )
+ ++ns;
+ else if( higgs[i] == 36 || higgs[i] == 46 )
+ ++np;
+ else if( abs(higgs[i]) == 37 )
+ ++nc;
+ }
+ //check three Higgs in vertex
+ assert( ns + np + nc == 3 );
+ if( q2 != _q2last ) {
+ _q2last = q2;
+ _glast = weakCoupling(q2);
+ _masslast.first = _theSM->mass(q2,getParticleData(5));
+ _masslast.second = _theSM->mass(q2,getParticleData(6));
+ }
+ //define VEV's
+ double rt = sqrt(0.5);
+ _mb = _masslast.first;
+ _mt = _masslast.second;
+ Energy _mtpole = getParticleData(6)->mass();
+ Energy2 Qstsb = _MQ3*_MU2;
+ double temp = Qstsb/sqr(_mtpole);
+ assert(temp!=0.);
+ double radlog = log(temp);
+ complex<Energy> coupling;
+ //CP even Higgs
+ if( ns == 3 ) {
+ unsigned int a = (higgs[0] - 25)/10;
+ unsigned int b = (higgs[1] - 25)/10;
+ unsigned int c = (higgs[2] - 25)/10;
+ coupling =
+ sqr(_lambda)*rt*(_vu*(usMix(a,b,c,1,0,0) + usMix(a,b,c,1,2,2))/_glast +
+ _vd*(usMix(a,b,c,0,1,1) + usMix(a,b,c,0,2,2))/_glast +
+ _s *(usMix(a,b,c,2,1,1) + usMix(a,b,c,2,0,0)))
+ - _lambda*_kappa*rt*(_vu*usMix(a,b,c,0,2,2)/_glast +
+ _vd*usMix(a,b,c,2,1,2)/_glast + 2.*_s*usMix(a,b,c,1,0,2))
+ + sqr(_kappa)/rt*_s*usMix(a,b,c,2,2,2)
+ - _lambda*_theAl*rt*usMix(a,b,c,1,0,2)
+ + _kappa*_theAk*rt/3.*usMix(a,b,c,2,2,2)
+ + sqr(_glast)*0.25*rt/sqr(_cw)*(_vu*(usMix(a,b,c,1,1,1) -
+ usMix(a,b,c,1,0,0))/_glast -
+ _vd*(usMix(a,b,c,0,1,1) -
+ usMix(a,b,c,0,0,0))/_glast);
+ // additional radiative terms
+ if(_includeRadiative) {
+ complex <Energy> radtop = usMix(a,b,c,1,1,1)*3.0*sqrt(2.0)*radlog
+ *sqr(_mt)*sqr(_mt)*sqr(_glast)*_glast/
+ (16.0*sqr(pi)*_vu*_vu*_vu);
+ complex <Energy> radbot= usMix(a,b,c,0,0,0)*3.0*sqrt(2.0)*radlog
+ *sqr(_mb)*sqr(_mb)*sqr(_glast)*_glast
+ /(16.0*sqr(pi)*_vd*_vd*_vd);
+ coupling += radbot + radtop;
+ }
+ }
+ //CP even, CP odd Vertex
+ else if(ns == 1 && np == 2) {
+ unsigned int a(0), b(0), c(0);
+ if( higgs[0] == 25 || higgs[0] == 35 || higgs[0] == 45 ) {
+ a = (higgs[0] - 25)/10;
+ b = (higgs[1] - 36)/10;
+ c = (higgs[2] - 36)/10;
+ }
+ else if(higgs[1] == 25 || higgs[1] == 35 || higgs[1] == 45 ) {
+ a = (higgs[1] - 25)/10;
+ b = (higgs[0] - 36)/10;
+ c = (higgs[2] - 36)/10;
+ }
+ else {
+ a = (higgs[2] - 25)/10;
+ b = (higgs[0] - 36)/10;
+ c = (higgs[1] - 36)/10;
+ }
+ coupling =
+ sqr(_lambda)*rt*(_vu*(upMix(a,b,c,1,0,0) + upMix(a,b,c,1,2,2))/_glast +
+ _vd*(upMix(a,b,c,0,1,1) + upMix(a,b,c,0,2,2))/_glast +
+ _s *(upMix(a,b,c,2,1,1) + upMix(a,b,c,2,0,0)))
+ + _lambda*_kappa*rt*(_vu*(upMix(a,b,c,0,2,2)
+ - 2.*upMix(a,b,c,2,0,2))/_glast +
+ _vd*(upMix(a,b,c,1,2,2)
+ - 2.*upMix(a,b,c,2,1,2))/_glast
+ + 2.*_s*(upMix(a,b,c,2,1,0)
+ - upMix(a,b,c,1,0,2) - upMix(a,b,c,0,1,2)))
+ + sqr(_kappa)/rt*_s*upMix(a,b,c,2,2,2)
+ +_lambda*_theAl*rt*(upMix(a,b,c,1,0,2)
+ + upMix(a,b,c,0,1,2) + upMix(a,b,c,2,1,0))
+ - _kappa*_theAk*rt*upMix(a,b,c,2,2,2)
+ + sqr(_glast)*0.25*rt/sqr(_cw)*(_vu*(upMix(a,b,c,1,1,1) -
+ upMix(a,b,c,1,0,0))/_glast -
+ _vd*(upMix(a,b,c,0,1,1) -
+ upMix(a,b,c,0,0,0))/_glast);
+ if(_includeRadiative) {
+ complex <Energy> radtop = upMix(a,b,c,1,1,1)*3.0*sqrt(2.0)*radlog*
+ sqr(_mt)*sqr(_mt)*sqr(_glast)*_glast/
+ (16.0*sqr(pi)*_vu*_vu*_vu);
+ complex <Energy> radbot= upMix(a,b,c,0,0,0)*3.0*sqrt(2.0)*radlog*
+ sqr(_mb)*sqr(_mb)*sqr(_glast)*_glast
+ /(16.0*sqr(pi)*_vd*_vd*_vd);
+ coupling += radbot + radtop;
+ }
+ }
+ //Charged Higgs
+ else {
+ unsigned int a(0);
+ if( higgs[0] == 25 || higgs[0] == 35 || higgs[0] == 45 )
+ a = (higgs[0] - 25)/10;
+ else if(higgs[1] == 25 || higgs[1] == 35 || higgs[1] == 45 )
+ a = (higgs[1] - 25)/10;
+ else
+ a = (higgs[2] - 25)/10;
+ coupling =
+ sqr(_lambda)*rt*2.*(_s*((*_mixS)(a,2)*sqr(_cb) + (*_mixS)(a,2)*sqr(_sb))
+ - (_vu*(*_mixS)(a,0)/_glast +
+ _vd*(*_mixS)(a,1)/_glast)*_sb*_cb)
+ +_lambda*_sb*_cb*2.*(*_mixS)(a,2)*(_kappa*_s/rt + rt*_theAl)
+ + sqr(_glast)*0.5*rt*_sw2/sqr(_cw)*((_vu*(*_mixS)(a,1)/_glast -
+ _vd*(*_mixS)(a,0)/_glast)*sqr(_cb) +
+ (_vd*(*_mixS)(a,0)/_glast -
+ _vu*(*_mixS)(a,1)/_glast)*sqr(_sb))
+ + sqr(_glast)*0.5*rt*(_vu*((*_mixS)(a,1)*sqr(_cb) +
+ (*_mixS)(a,1)*sqr(_sb) +
+ 2.*(*_mixS)(a,0)*_cb*_sb)/_glast
+ + _vd*((*_mixS)(a,0)*sqr(_cb) +
+ (*_mixS)(a,0)*sqr(_sb)
+ + 2.*(*_mixS)(a,1)*_sb*_cb)/_glast);
+ if(_includeRadiative) {
+ complex <Energy> radtop =(*_mixS)(a,1)*sqr(_sb)*6.0*sqrt(2.0)*radlog*
+ sqr(_mt)*sqr(_mt)*sqr(_glast)*_glast/
+ (16.0*sqr(pi)*_vu*_vu*_vu);
+ complex <Energy> radbot=(*_mixS)(a,0)*sqr(_cb)*6.0*sqrt(2.0)*radlog*
+ sqr(_mb)*sqr(_mb)*sqr(_glast)*_glast
+ /(16.0*sqr(pi)*_vd*_vd*_vd);
+ complex <Energy> temp2 = _vu*((*_mixS)(a,1)*sqr(_cb) +
+ (*_mixS)(a,0)*_sb*_cb)/_glast+
+ _vd*((*_mixS)(a,1)*_sb*_cb +
+ (*_mixS)(a,0)*sqr(_sb))/_glast;
+
+ complex <Energy> radtopbot= temp2*6.0*sqrt(2.0)*radlog*
+ sqr(_mt)*sqr(_mb)*sqr(_glast)*sqr(_glast)
+ /(16.0*sqr(pi)*sqr(_vu)*sqr(_vd));
+ coupling += radbot + radtop + radtopbot;
+ }
+ }
+ norm(-coupling * UnitRemoval::InvE);
+}
diff --git a/Models/Susy/NMSSM/NMSSMHHHVertex.h b/Models/Susy/NMSSM/NMSSMHHHVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMHHHVertex.h
@@ -0,0 +1,324 @@
+// -*- C++ -*-
+#ifndef HERWIG_NMSSMHHHVertex_H
+#define HERWIG_NMSSMHHHVertex_H
+//
+// This is the declaration of the NMSSMHHHVertex class.
+//
+#include "ThePEG/Helicity/Vertex/Scalar/SSSVertex.h"
+#include "ThePEG/StandardModel/StandardModelBase.h"
+#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/Susy/MixingMatrix.h"
+
+namespace Herwig {
+using namespace ThePEG;
+using namespace ThePEG::Helicity;
+
+/** \ingroup Helicity
+ * The <code>NMSSMHHHVertex<\code> defines the triple
+ * Higgs coupling in the NMSSM.
+ *
+ * @see \ref NMSSMHHHVertexInterfaces "The interfaces"
+ * defined for NMSSMHHHVertex.
+ */
+class NMSSMHHHVertex: public SSSVertex {
+
+public:
+
+ /** @name Standard constructors and destructors. */
+ //@{
+ /**
+ * The default constructor.
+ */
+ NMSSMHHHVertex();
+
+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();
+
+ /**
+ * Calculate the couplings. This method is virtual and must be implemented in
+ * classes inheriting from this.
+ * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
+ * @param part1 The ParticleData pointer for the first particle.
+ * @param part2 The ParticleData pointer for the second particle.
+ * @param part3 The ParticleData pointer for the third particle.
+ */
+ virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,
+ tcPDPtr part3);
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const {return new_ptr(*this);}
+
+ /** 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 {return new_ptr(*this);}
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is a concrete class with persistent data.
+ */
+ static ClassDescription<NMSSMHHHVertex> initNMSSMHHHVertex;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ NMSSMHHHVertex & operator=(const NMSSMHHHVertex &);
+
+private:
+
+ /**
+ * The mixing matrix combination \f$U^S_{ai}U^S_{bj}U^S_{ck}\f$
+ * @param a The row element of the first CP-even mixing matrix
+ * @param b The column element of the first CP-even mixing matrix
+ * @param c The row element of the second CP-even mixing matrix
+ * @param i The column element of the second CP-even mixing matrix
+ * @param j The row element of the third CP-even mixing matrix
+ * @param k The column element of the third CP-even mixing matrix
+ */
+ Complex usMix(unsigned int a, unsigned int b, unsigned int c,
+ unsigned int i, unsigned int j, unsigned int k) const {
+ return (*_mixS)(a,i)*(*_mixS)(b,j)*(*_mixS)(c,k) +
+ (*_mixS)(a,i)*(*_mixS)(c,j)*(*_mixS)(b,k) +
+ (*_mixS)(b,i)*(*_mixS)(a,j)*(*_mixS)(c,k) +
+ (*_mixS)(b,i)*(*_mixS)(c,j)*(*_mixS)(a,k) +
+ (*_mixS)(c,i)*(*_mixS)(a,j)*(*_mixS)(b,k) +
+ (*_mixS)(c,i)*(*_mixS)(b,j)*(*_mixS)(a,k);
+ }
+
+ /**
+ * The mixing matrix combination \f$U^S_{ai}U^P_{bj}U^P_{ck}\f$
+ * @param a The row element of the first CP-even mixing matrix
+ * @param b The column element of the first CP-even mixing matrix
+ * @param c The row element of the second CP-even mixing matrix
+ * @param i The column element of the second CP-even mixing matrix
+ * @param j The row element of the third CP-even mixing matrix
+ * @param k The column element of the third CP-even mixing matrix
+ */
+ Complex upMix(unsigned int a, unsigned int b, unsigned int c,
+ unsigned int i, unsigned int j, unsigned int k) const {
+ return (*_mixS)(a,i)*((*_mixP)(b,j)*(*_mixP)(c,k) +
+ (*_mixP)(c,j)*(*_mixP)(b,k));
+ }
+
+private:
+
+ /**
+ * A pointer to the object containing the SM parameters
+ */
+ tcHwSMPtr _theSM;
+
+ /**
+ * The \f$W\f$ mass
+ */
+ Energy _mw;
+
+ /**
+ * The \f$Z\f$ mass
+ */
+ Energy _mz;
+ /**
+ * The \f$b\f$ mass
+ */
+ Energy _mb;
+
+ /**
+ * The \f$t\f$ mass
+ */
+ Energy _mt;
+
+ /**
+ * \f$\sin^2\theta_W\f$
+ */
+ double _sw2;
+
+ /**
+ * \f$\cos\theta_W\f$
+ */
+ double _cw;
+
+ /**
+ * The CP-even Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixS;
+
+ /**
+ * The CP-odd Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixP;
+
+ /**
+ * The coefficient of the trilinear \f$SH_2 H_1\f$ term in the superpotential
+ */
+ double _lambda;
+
+ /**
+ * The coefficient of the cubic singlet term in the superpotential
+ */
+ double _kappa;
+
+ /**
+ * The product \f$\lambda \langle S\rangle \f$.
+ */
+ Energy _lambdaVEV;
+
+ /**
+ * The soft trilinear \f$SH_2 H_1\f$ coupling
+ */
+ Energy _theAl;
+
+ /**
+ * The soft cubic \f$S\f$ coupling
+ */
+ Energy _theAk;
+
+ /**
+ * \f$\sin\beta\f$
+ */
+ double _sb;
+
+ /**
+ * \f$\cos\beta\f$
+ */
+ double _cb;
+
+ /**
+ * \f$\sin2\beta\f$
+ */
+ double _s2b;
+
+ /**
+ * \f$\cos2\beta\f$
+ */
+ double _c2b;
+
+ /**
+ * The value of the VEV of the higgs that couples to the down-type sector
+ * \f$ g*sqrt(2)M_W\cos\beta \f$
+ */
+ Energy _vu;
+
+ /**
+ * The value of the VEV of the higgs that couples to the up-type sector
+ * i.e. \f$ g*sqrt(2)M_W\sin\beta \f$
+ */
+ Energy _vd;
+
+ /**
+ * The value of the VEV of the singlet higgs
+ */
+ Energy _s;
+
+ /**
+ * The scale at which this vertex was last evaluated
+ */
+ Energy2 _q2last;
+
+ /**
+ * The value of the EW coupling when it was last evaluated
+ */
+ double _glast;
+ /**
+ * left 3rd generation scalar quark mass
+ */
+ Energy _MQ3;
+
+ /**
+ * right scalar top mass
+ */
+ Energy _MU2;
+ /**
+ * The mass of the last fermion for which the coupling was evaluated.
+ */
+ pair<Energy,Energy> _masslast;
+
+ /**
+ * Whether or onto to include the radiative terms
+ */
+ bool _includeRadiative;
+
+};
+
+}
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of NMSSMHHHVertex. */
+template <>
+struct BaseClassTrait<Herwig::NMSSMHHHVertex,1> {
+ /** Typedef of the first base class of NMSSMHHHVertex. */
+ typedef Helicity::SSSVertex NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the NMSSMHHHVertex class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::NMSSMHHHVertex>
+ : public ClassTraitsBase<Herwig::NMSSMHHHVertex> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::NMSSMHHHVertex"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * NMSSMHHHVertex is implemented. It may also include several, space-separated,
+ * libraries if the class NMSSMHHHVertex 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 "HwSusy.so HwNMSSM.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_NMSSMHHHVertex_H */
diff --git a/Models/Susy/NMSSM/NMSSMHSFSFVertex.cc b/Models/Susy/NMSSM/NMSSMHSFSFVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMHSFSFVertex.cc
@@ -0,0 +1,386 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the NMSSMHSFSFVertex class.
+//
+
+#include "NMSSMHSFSFVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "NMSSM.h"
+
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+NMSSMHSFSFVertex::NMSSMHSFSFVertex() :
+ _triTp(0.*MeV), _triBt(0.*MeV), _triTa(0.*MeV), _lambda(0.),
+ _lambdaVEV(0.*MeV), _v1(0.*MeV), _v2(0.*MeV), _sw(0.), _cw(0.),
+ _mw(0.*MeV), _mz(0.*MeV), _sb(0.), _cb(0.), _tb(0.), _q2last(0.*MeV2),
+ _couplast(0.), _masslast(make_pair(0.*MeV,0.*MeV)), _idlast(make_pair(0,0)) {
+ //CP even
+ int even[3] = {25, 35, 45};
+ for(size_t h = 0; h < 3; ++h ) {
+ //squarks
+ for(long q = 1; q < 7; ++q) {
+ //11
+ addToList(even[h], -1000000 - q, 1000000 + q);
+ //22
+ addToList(even[h], -2000000 - q, 2000000 + q);
+ //12
+ addToList(even[h], -1000000 - q, 2000000 + q);
+ //21
+ addToList(even[h], -2000000 - q, 1000000 + q);
+ }
+ //sleptons
+ for(long l = 11; l < 17; ++l) {
+ //11
+ addToList(even[h], -1000000 - l, 1000000 + l);
+ //no right handed sneutrinos
+ if( l % 2 != 0 ) {
+ //22
+ addToList(even[h], -2000000 - l, 2000000 + l);
+ //12
+ addToList(even[h], -1000000 - l, 2000000 + l);
+ //21
+ addToList(even[h], -2000000 - l, 1000000 + l);
+ }
+ }
+ }
+ //CP odd
+ int odd[2] = {36, 46};
+ for(size_t h = 0; h < 2; ++h ) {
+ //squarks
+ for(long q = 1; q < 7; ++q) {
+ //12
+ addToList(odd[h], -1000000 - q, 2000000 + q);
+ //21
+ addToList(odd[h], -2000000 - q, 1000000 + q);
+ }
+ //sleptons
+ for(long l = 11; l < 16; l += 2) {
+ //12
+ addToList(odd[h], -1000000 - l, 2000000 + l);
+ //21
+ addToList(odd[h], -2000000 - l, 1000000 + l);
+ }
+ }
+ //charged higgs
+ //squarks
+ for(long q = 1; q < 3; ++q ) {
+ //H-
+ //LL
+ addToList(-37, -2*q - 999999, 2*q + 1000000);
+ //RR
+ addToList(-37, -2*q - 1999999, 2*q + 2000000);
+ //LR
+ addToList(-37, -2*q - 999999, 2*q + 2000000);
+ //RL
+ addToList(-37, -2*q - 1999999, 2*q + 1000000);
+ //H+
+ //LL
+ addToList(37, -2*q - 1000000, 2*q + 999999);
+ //RR
+ addToList(37, -2*q - 2000000, 2*q + 1999999);
+ //LR
+ addToList(37, -2*q - 1000000, 2*q + 1999999);
+ //RL
+ addToList(37, -2*q - 2000000, 2*q + 999999);
+ }
+ //sleptons
+ //easier as there are no right handed sneutrinos
+ for(long l = 11; l <= 15; l +=2 ) {
+ //H-
+ //LL
+ addToList(-37, -l - 1000000, l + 1000001);
+ //RL
+ addToList(-37, -l - 2000000, l + 1000001);
+ //H+
+ //LL
+ addToList(+37, -l - 1000001, l + 1000000);
+ //RL
+ addToList(+37, -l - 1000001, l + 2000000);
+ }
+}
+
+void NMSSMHSFSFVertex::persistentOutput(PersistentOStream & os) const {
+ os << _theSM << _mixS << _mixP << _mixTp << _mixBt << _mixTa
+ << ounit(_triTp,GeV) << ounit(_triBt,GeV) << ounit(_triTa,GeV)
+ << _lambda << ounit(_lambdaVEV,GeV) << ounit(_v1,GeV) << ounit(_v2,GeV)
+ << _sw << _cw << ounit(_mw,GeV) << ounit(_mz,GeV) << _sb << _cb
+ << _tb;
+}
+
+
+void NMSSMHSFSFVertex::persistentInput(PersistentIStream & is, int) {
+ is >> _theSM >> _mixS >> _mixP >> _mixTp >> _mixBt >> _mixTa
+ >> iunit(_triTp,GeV) >> iunit(_triBt,GeV) >> iunit(_triTa,GeV)
+ >> _lambda >> iunit(_lambdaVEV,GeV) >> iunit(_v1,GeV) >> iunit(_v2,GeV)
+ >> _sw >> _cw >> iunit(_mw,GeV) >> iunit(_mz,GeV) >> _sb >> _cb >> _tb;
+}
+
+ClassDescription<NMSSMHSFSFVertex> NMSSMHSFSFVertex::initNMSSMHSFSFVertex;
+// Definition of the static class description member.
+
+void NMSSMHSFSFVertex::Init() {
+
+ static ClassDocumentation<NMSSMHSFSFVertex> documentation
+ ("The coupling of Higgs bosons to sfermions in the MSSM.");
+
+}
+
+void NMSSMHSFSFVertex::doinit() {
+ _theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
+ tcNMSSMPtr nmssm = dynamic_ptr_cast<tcNMSSMPtr>(_theSM);
+ if( !nmssm )
+ throw InitException() << "NMSSMHSFSFVertex::doinit() - The model pointer "
+ << "in this vertex is not an NMSSM one as it "
+ << "should be." << Exception::runerror;
+ _mixS = nmssm->CPevenHiggsMix();
+ _mixP = nmssm->CPoddHiggsMix();
+ _mixTp = nmssm->stopMix();
+ _mixBt = nmssm->sbottomMix();
+ _mixTa = nmssm->stauMix();
+ if( !_mixS || !_mixP || !_mixTp || !_mixBt || !_mixTa )
+ throw InitException()
+ << "NMSSMHSFSFVertex::doinit() - One of the mixing matrix pointers is "
+ << "null, cannot continue. CP-even: " << _mixS << " CP-odd: " << _mixP
+ << " ~t: " << _mixTp << " ~b: " << _mixBt << " ~tau: " << _mixTa
+ << Exception::runerror;
+
+ _triTp = nmssm->topTrilinear();
+ _triBt = nmssm->bottomTrilinear();
+ _triTa = nmssm->tauTrilinear();
+
+ _lambda = nmssm->lambda();
+ _lambdaVEV = nmssm->lambdaVEV();
+
+ _sw = sin2ThetaW();
+ _cw = sqrt( 1. - _sw);
+ _sw = sqrt(_sw);
+ _mw = getParticleData(24)->mass();
+ _mz = getParticleData(23)->mass();
+
+ _tb = nmssm->tanBeta();
+ double beta = atan(_tb);
+ _sb = sin(beta);
+ _cb = cos(beta);
+
+ _v1 = sqrt(2.)*_mw*_cb;
+ _v2 = sqrt(2.)*_mw*_sb;
+
+ orderInGem(1);
+ orderInGs(0);
+ SSSVertex::doinit();
+
+}
+
+void NMSSMHSFSFVertex::setCoupling(Energy2 q2,tcPDPtr part1,
+ tcPDPtr part2, tcPDPtr part3) {
+ // extract particle ids
+ long higgs(part1->id()), isf1(part2->id()), isf2(part3->id());
+ // higgs first
+ if(abs(isf1)<100) swap(higgs,isf1);
+ if(abs(isf2)<100) swap(higgs,isf2);
+ // squark second
+ if(isf1<0) swap(isf1,isf2);
+ // check higgs
+ assert( higgs == 25 || higgs == 35 || higgs == 45 ||
+ higgs == 36 || higgs == 46 || abs(higgs) == 37 );
+ // abs of antisquark and check
+ isf2 *=-1;
+ assert(isf1>0&&isf2>0);
+ // running coupling
+ if( q2 != _q2last ) {
+ _q2last = q2;
+ _couplast = weakCoupling(q2);
+ }
+
+ //charged higgs
+ if( abs(higgs) == 37 ) {
+ norm(_couplast*chargedHiggs(q2, isf1, isf2));
+ return;
+ }
+ // neutral higgs
+ // L/R states of the sfermions
+ unsigned int alpha = ( isf1 > 2000000 ) ? 1 : 0;
+ unsigned int beta = ( isf2 > 2000000 ) ? 1 : 0;
+ // id nad mass of corresponding SM fermion
+ long smid = ( alpha == 0 ) ? isf1 - 1000000 : isf1 - 2000000;
+ if( q2 != _q2last || smid != _idlast.first) {
+ _idlast.first = smid;
+ _masslast.first = _theSM->mass(q2, getParticleData(smid));
+ }
+ double f1 = _masslast.first/_mw;
+ complex<Energy> af(ZERO);
+ Complex m1a(0.), m1b(0.), m2a(0.), m2b(0.);
+ // mixing for down type squarks and charged sleptons
+ if( smid % 2 != 0 ) {
+ f1 /= _cb;
+ // sbottom
+ if( smid == 5 ) {
+ m1a = (*_mixBt)(alpha, 0);
+ m1b = (*_mixBt)(alpha, 1);
+ m2a = (*_mixBt)(beta , 0) ;
+ m2b = (*_mixBt)(beta , 1);
+ af = _triBt;
+ }
+ // stau
+ else if( smid == 15 ) {
+ m1a = (*_mixTa)(alpha, 0);
+ m1b = (*_mixTa)(alpha, 1);
+ m2a = (*_mixTa)(beta , 0) ;
+ m2b = (*_mixTa)(beta , 1);
+ af = _triTa;
+ }
+ // 1st 2 generations
+ else {
+ m1a = (alpha == 0) ? 1. : 0.;
+ m1b = (alpha == 0) ? 0. : 1.;
+ m2a = (beta == 0) ? 1. : 0.;
+ m2b = (beta == 0) ? 0. : 1.;
+ af = ZERO;
+ }
+ }
+ // mixing for up type squarks and sneutrions
+ else {
+ f1 /= _sb;
+ // stop
+ if( smid == 6 ) {
+ m1a = (*_mixTp)(alpha, 0);
+ m1b = (*_mixTp)(alpha, 1);
+ m2a = (*_mixTp)(beta , 0);
+ m2b = (*_mixTp)(beta , 1);
+ af = _triTp;
+ }
+ // everything else
+ else {
+ m1a = (alpha == 0) ? 1. : 0.;
+ m1b = (alpha == 0) ? 0. : 1.;
+ m2a = (beta == 0) ? 1. : 0.;
+ m2b = (beta == 0) ? 0. : 1.;
+ af = 0.*MeV;
+ }
+ }
+ // scalar higgs bosons
+ complex<Energy> fact(ZERO);
+ if( higgs == 25 || higgs == 35 || higgs == 45 ) {
+ int iloc = (higgs - 25)/10;
+ complex<Energy> f2 = 0.5*_mz*( - _cb*(*_mixS)(iloc,0)
+ + _sb*(*_mixS)(iloc,1))/_cw;
+
+ // down type squarks and charged sleptons
+ if( smid % 2 != 0 ) {
+ double ef = (smid < 7) ? -1./3. : -1.;
+ fact = - f2*( (1. + 2.*ef*sqr(_sw))*m1a*m2a - 2.*ef*sqr(_sw)*m1b*m2b)
+ - f1*_masslast.first*(*_mixS)(iloc,0)*(m1a*m2a + m1b*m2b)
+ - 0.5*f1*(( - _lambdaVEV*(*_mixS)(iloc,1)
+ - _lambda*_v2*(*_mixS)(iloc,2)/_couplast
+ + af*(*_mixS)(iloc,0)) *
+ (m2a*m1b + m1a*m2b) );
+ }
+ // up type squarks and sneutrinos
+ else {
+ double ef = (smid < 7) ? 2./3. : 0.;
+ fact = +f2*( (1. - 2.*ef*sqr(_sw))*m1a*m2a + 2.*ef*sqr(_sw)*m1b*m2b )
+ - f1*_masslast.first*(*_mixS)(iloc,1)*(m1a*m2a + m1b*m2b)
+ - 0.5*f1*(( - _lambdaVEV*(*_mixS)(iloc,0)
+ - _lambda*_v1*(*_mixS)(iloc,2)/_couplast
+ + af*(*_mixS)(iloc,1) ) *
+ (m2a*m1b + m1a*m2b));
+ }
+ }
+ // pseudo scalar
+ else if( higgs == 36 || higgs == 46 ) {
+ int iloc = (higgs - 36)/10;
+ // down type squarks and charged sleptons
+ if( smid % 2 != 0 ) {
+ fact = -0.5*f1*Complex(0.0,1.0)*
+ ( _lambdaVEV*(*_mixP)(iloc,1) +
+ _lambda*_v2*(*_mixP)(iloc,2)/_couplast +
+ af*(*_mixP)(iloc,0) );
+ }
+ // up-type squarks and sneutrinos
+ else {
+ fact =-0.5*f1*Complex(0.0,1.0)*
+ ( _lambdaVEV *(*_mixP)(iloc,0) +
+ _lambda*_v1*(*_mixP)(iloc,2)/_couplast +
+ af*(*_mixP)(iloc,1));
+ }
+ if(alpha<beta) fact *= -1.;
+ }
+ norm(_couplast*fact*UnitRemoval::InvE);
+}
+
+Complex NMSSMHSFSFVertex::chargedHiggs(Energy2 q2, long id1, long id2) {
+ //have id1 as up-type
+ if( id1 % 2 != 0) swap(id1, id2);
+ // sfermion L/R states
+ unsigned int alpha = ( id1/1000000 == 2 ) ? 1 : 0;
+ unsigned int beta = ( id2/1000000 == 2 ) ? 1 : 0;
+ // type of quarks
+ long utype = (alpha == 0) ? id1 - 1000000 : id1 - 2000000;
+ long dtype = ( beta == 0) ? id2 - 1000000 : id2 - 2000000;
+ // compute the running masses
+ if( q2 != _q2last || id1 != _idlast.first || id2 != _idlast.second) {
+ _idlast.first = id1;
+ _idlast.second = id2;
+ _masslast.first = _theSM->mass(q2, getParticleData(utype) );
+ _masslast.second = _theSM->mass(q2, getParticleData(dtype) );
+ }
+ Energy2 facta = 2.*sqr(_mw)*_sb*_cb;
+ complex<Energy2> coupling(ZERO);
+ // sleptons
+ if( dtype == 11 || dtype == 13 || dtype == 15) {
+ Complex l1b = 0., l2b = 0.;
+ complex<Energy> tri(ZERO);
+ // 1st 2 generations
+ if (dtype == 11 || dtype == 13) {
+ l1b = (beta == 0) ? 1.0 : 0.0;
+ l2b = (beta == 0) ? 0.0 : 1.0;
+ }
+ // stau
+ else {
+ l1b = (*_mixTa)(beta, 0) ;
+ l2b = (*_mixTa)(beta, 1);
+ tri = _triTa;
+ }
+ coupling = ( l1b*(sqr(_masslast.second)*_tb - facta) +
+ l2b*_masslast.second*(tri*_tb + _lambdaVEV) );
+ }
+ // squarks
+ else {
+ Complex q1a(0.0), q1b(0.0), q2a(0.0), q2b(0.0);
+ complex<Energy> triD(ZERO), triU(ZERO);
+ // up-type bit
+ // stop
+ if(utype == 6){
+ q1a = (*_mixTp)(alpha, 0) ;
+ q2a = (*_mixTp)(alpha, 1);
+ triU = _triTp;
+ }
+ // light
+ else{
+ q1a = (alpha == 0) ? 1.0 : 0.0;
+ q2a = (alpha == 0) ? 0.0 : 1.0;
+ }
+ // down-type bit
+ // sbottom
+ if(utype == 5){
+ q1b = (*_mixBt)(beta, 0) ;
+ q2b = (*_mixBt)(beta, 1);
+ }
+ // light
+ else{
+ q1b = (beta == 0) ? 1.0 : 0.0;
+ q2b = (beta == 0) ? 0.0 : 1.0;
+ }
+ Energy mfu = _masslast.first;
+ Energy mfd = _masslast.second;
+ coupling = ( q1a*q1b*((sqr(mfd)*_tb + sqr(mfu)/_tb) - facta)
+ + q2a*q2b*mfu*mfd*(_tb + (1./_tb))
+ + q1a*q2b*mfd*(triD*_tb + _lambdaVEV)
+ + q2a*q1b*mfu*(triU/_tb + _lambdaVEV));
+ }
+ return coupling * UnitRemoval::InvE/_mw/sqrt(2.);
+}
diff --git a/Models/Susy/NMSSM/NMSSMHSFSFVertex.h b/Models/Susy/NMSSM/NMSSMHSFSFVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMHSFSFVertex.h
@@ -0,0 +1,287 @@
+// -*- C++ -*-
+#ifndef HERWIG_NMSSMHSFSFVertex_H
+#define HERWIG_NMSSMHSFSFVertex_H
+//
+// This is the declaration of the NMSSMHSFSFVertex class.
+//
+
+#include "ThePEG/Helicity/Vertex/Scalar/SSSVertex.h"
+#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/Susy/MixingMatrix.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/** \ingroup Helicity
+ * This class defines the coupling of Higgs bosons to the sfermions
+ * in the NMSSM.
+ *
+ * @see \ref NMSSMHSFSFVertexInterfaces "The interfaces"
+ * defined for NMSSMHSFSFVertex.
+ */
+class NMSSMHSFSFVertex: public Helicity::SSSVertex {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ NMSSMHSFSFVertex();
+
+ /** @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();
+
+ /**
+ * Calculate the couplings. This method is virtual and must be implemented in
+ * classes inheriting from this.
+ * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
+ * @param part1 The ParticleData pointer for the first particle.
+ * @param part2 The ParticleData pointer for the second particle.
+ * @param part3 The ParticleData pointer for the third particle.
+ */
+ virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,
+ tcPDPtr part3);
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const {return new_ptr(*this);}
+
+ /** 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 {return new_ptr(*this);}
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is a concrete class with persistent data.
+ */
+ static ClassDescription<NMSSMHSFSFVertex> initNMSSMHSFSFVertex;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ NMSSMHSFSFVertex & operator=(const NMSSMHSFSFVertex &);
+
+private:
+
+ /**
+ * Return the coupling of the charged higgs to the sfermions
+ */
+ Complex chargedHiggs(Energy2 q2, long id1, long id2);
+
+private:
+
+ /** @name Stored parameters for fast access.*/
+ //@{
+ /**
+ * A pointer to the Standard Model object
+ */
+ tcHwSMPtr _theSM;
+
+ /**
+ * The CP-even higgs mixing matrix
+ */
+ MixingMatrixPtr _mixS;
+
+ /**
+ * The CP-odd higgs mixing matrix
+ */
+ MixingMatrixPtr _mixP;
+
+ /**
+ * The \f$ \tilde{t}\f$ mixing matrix
+ */
+ MixingMatrixPtr _mixTp;
+
+ /**
+ * The \f$ \tilde{b}\f$ mixing matrix
+ */
+ MixingMatrixPtr _mixBt;
+
+ /**
+ * The \f$ \tilde{\tau}\f$ mixing matrix
+ */
+ MixingMatrixPtr _mixTa;
+
+ /**
+ * The top quark trilinear coupling
+ */
+ complex<Energy> _triTp;
+
+ /**
+ * The bottom quark trilinear coupling
+ */
+ complex<Energy> _triBt;
+
+ /**
+ * The tau lepton trilinear coupling
+ */
+ complex<Energy> _triTa;
+
+ /**
+ * The value of \f$\lambda\f$.
+ */
+ double _lambda;
+
+ /**
+ * The value of \f$ \lambda <S> \f$, the
+ * V.E.V of the extra gauge singlet scaled
+ * by \f$\lambda\f$
+ */
+ Energy _lambdaVEV;
+
+ /**
+ * The value of the V.E.V \f$ v_1 \f$
+ */
+ Energy _v1;
+
+ /**
+ * The value of the V.E.V \f$ v_2 \f$
+ */
+ Energy _v2;
+
+ /**
+ * The value of \f$ \sin\theta_W \f$
+ */
+ double _sw;
+
+ /**
+ * The value of \f$ \cos\theta_W \f$
+ */
+ double _cw;
+
+ /**
+ * The value of \f$ M_W \f$
+ */
+ Energy _mw;
+
+ /**
+ * The value of \f$ M_Z \f$
+ */
+ Energy _mz;
+
+ /**
+ * The value of \f$ \sin\beta \f$
+ */
+ double _sb;
+
+ /**
+ * The value of \f$ \cos\beta \f$
+ */
+ double _cb;
+
+ /**
+ * The value of \f$ \tan\beta \f$
+ */
+ double _tb;
+
+ //@}
+
+ /** @name Store previously calculated values for speed. */
+ //@{
+ /**
+ * The scale at which the last calculation took place.
+ */
+ Energy2 _q2last;
+
+ /**
+ * The value of the dimensionless coupling \f$g_W\f$ when
+ * last calculated.
+ */
+ double _couplast;
+
+ /**
+ * The value of mass of the counterpart SM fermion when
+ * last calculated.
+ */
+ pair<Energy, Energy> _masslast;
+
+ /**
+ * The PDG codes of the particles in the vertex when it was last evaluated
+ */
+ pair<long, long> _idlast;
+ //@}
+};
+}
+
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of NMSSMHSFSFVertex. */
+template <>
+struct BaseClassTrait<Herwig::NMSSMHSFSFVertex,1> {
+ /** Typedef of the first base class of NMSSMHSFSFVertex. */
+ typedef ThePEG::Helicity::SSSVertex NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the NMSSMHSFSFVertex class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::NMSSMHSFSFVertex>
+ : public ClassTraitsBase<Herwig::NMSSMHSFSFVertex> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::NMSSMHSFSFVertex"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * NMSSMHSFSFVertex is implemented. It may also include several, space-separated,
+ * libraries if the class NMSSMHSFSFVertex 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 "HwSusy.so HwNMSSM.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_NMSSMHSFSFVertex_H */
diff --git a/Models/Susy/NMSSM/NMSSMPPHVertex.cc b/Models/Susy/NMSSM/NMSSMPPHVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMPPHVertex.cc
@@ -0,0 +1,277 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the NMSSMPPHVertex class.
+//
+
+#include "NMSSMPPHVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "Herwig++/Models/Susy/NMSSM/NMSSM.h"
+
+using namespace Herwig;
+
+NMSSMPPHVertex::NMSSMPPHVertex()
+ : _sw(0.), _cw(0.), _mw(0.*MeV),
+ _mz(0.*MeV),_lambdaVEV(0.*MeV), _lambda(0.),
+ _triTp(0.*MeV), _triBt(0.*MeV),
+ _sb(0.), _cb(0.),
+ _kappa(0.),_vu(ZERO),_vd(ZERO),_s(ZERO),_theAl(ZERO),
+ _masslast(make_pair(0.*MeV,0.*MeV)),_q2last(0.*MeV2),
+ _couplast(0.), _coup(0.), _hlast(0), _recalc(true) {
+ addToList(22,22,25);
+ addToList(22,22,35);
+ addToList(22,22,36);
+ addToList(22,22,45);
+ addToList(22,22,46);
+}
+
+void NMSSMPPHVertex::doinit() {
+ _theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
+ if( !_theSM ) {
+ throw InitException() << "NMSSMPPHVertex::doinit - The SM pointer is null!"
+ << Exception::abortnow;
+ }
+ // SM parameters
+ _sw = sqrt(sin2ThetaW());
+ _cw = sqrt(1. - sin2ThetaW());
+ _mw = getParticleData(24)->mass();
+ _mz = getParticleData(23)->mass();
+ _top = getParticleData(6);
+ _bt = getParticleData(5);
+ _tau = getParticleData(15);
+
+ //NMSSM parameters
+ tcNMSSMPtr nmssm = dynamic_ptr_cast<tcNMSSMPtr>(_theSM);
+ _mixS = nmssm->CPevenHiggsMix();
+ _mixP = nmssm->CPoddHiggsMix();
+ _mixQt = nmssm->stopMix();
+ _mixQb = nmssm->sbottomMix();
+ _mixLt = nmssm->stauMix();
+
+ double beta = atan(nmssm->tanBeta());
+ _sb = sin(beta);
+ _cb = cos(beta);
+
+ _lambda = nmssm->lambda();
+ _lambdaVEV = nmssm->lambdaVEV();
+
+ _triTp = nmssm->topTrilinear();
+ _triBt = nmssm->bottomTrilinear();
+ _triTa = nmssm->tauTrilinear();
+
+ _vd = sqrt(2)*_mw*_cb;
+ _vu = sqrt(2)*_mw*_sb;
+ _s = _lambdaVEV/_lambda;
+ _theAl = nmssm->trilinearLambda();
+ _kappa = nmssm->kappa();
+
+ _mixU = nmssm->charginoUMix();
+ _mixV = nmssm->charginoVMix();
+
+ // resize vectors here and use setNParticles method
+ // to the set the actual number in the loop.
+ // Also only the top mass hass to be calculated at runtime
+ masses.resize(13, Energy());
+ masses[ 0] = getParticleData( 6)->mass();
+ masses[ 1] = getParticleData( 5)->mass();
+ masses[ 2] = getParticleData(15)->mass();
+ masses[ 3] = getParticleData(ParticleID::SUSY_chi_1plus)->mass();
+ masses[ 4] = getParticleData(ParticleID::SUSY_chi_2plus)->mass();
+ masses[ 5] = _mw;
+ masses[ 6] = getParticleData(ParticleID::Hplus)->mass();
+ masses[ 7] = getParticleData(1000005)->mass();
+ masses[ 8] = getParticleData(2000005)->mass();
+ masses[ 9] = getParticleData(1000006)->mass();
+ masses[10] = getParticleData(2000006)->mass();
+ masses[11] = getParticleData(1000015)->mass();
+ masses[12] = getParticleData(2000015)->mass();
+ type.resize(13, PDT::Spin0);
+ type[0] = PDT::Spin1Half;
+ type[1] = PDT::Spin1Half;
+ type[2] = PDT::Spin1Half;
+ type[3] = PDT::Spin1Half;
+ type[4] = PDT::Spin1Half;
+ type[5] = PDT::Spin1;
+ couplings.resize(13);
+ orderInGem(3);
+ orderInGs(0);
+ VVSLoopVertex::doinit();
+}
+
+void NMSSMPPHVertex::persistentOutput(PersistentOStream & os) const {
+ os << _theSM << _sw << _cw << ounit(_mw, GeV) << ounit(_mz, GeV)
+ << ounit(_lambdaVEV,GeV) << _lambda
+ << ounit(_triTp,GeV) << ounit(_triBt,GeV) << ounit(_triTa,GeV)
+ << _top << _bt << _tau << _mixS << _mixP << _mixU << _mixV
+ << _mixQt << _mixQb << _mixLt << _sb << _cb << _kappa
+ << ounit(_vu,GeV) << ounit(_vd,GeV) << ounit(_s,GeV) << ounit(_theAl,GeV);
+}
+
+void NMSSMPPHVertex::persistentInput(PersistentIStream & is, int) {
+ is >> _theSM >> _sw >> _cw >> iunit(_mw, GeV) >> iunit(_mz, GeV)
+ >> iunit(_lambdaVEV,GeV) >> _lambda
+ >> iunit(_triTp,GeV) >> iunit(_triBt,GeV) >> iunit(_triTa,GeV)
+ >> _top >> _bt >> _tau >> _mixS >> _mixP >> _mixU >> _mixV
+ >> _mixQt >> _mixQb >> _mixLt >> _sb >> _cb >> _kappa
+ >> iunit(_vu,GeV) >> iunit(_vd,GeV) >> iunit(_s,GeV) >> iunit(_theAl,GeV);
+}
+
+ClassDescription<NMSSMPPHVertex> NMSSMPPHVertex::initNMSSMPPHVertex;
+// Definition of the static class description member.
+
+void NMSSMPPHVertex::Init() {
+
+ static ClassDocumentation<NMSSMPPHVertex> documentation
+ ("The effective coupling of a higgs to a pair of gluons in the "
+ "NMSSM.");
+
+}
+
+void NMSSMPPHVertex::setCoupling(Energy2 q2, tcPDPtr p1, tcPDPtr p2,
+ tcPDPtr p3) {
+ long hid(p3->id());
+ double rt = sqrt(0.5);
+ if( q2 != _q2last ) {
+ _couplast = sqr(electroMagneticCoupling(q2));
+ _coup = weakCoupling(q2);
+ _q2last = q2;
+ _recalc = true;
+ }
+ norm(_couplast*_coup);
+ // scalar higgs bosons
+ if( hid != _hlast ) {
+ _hlast = hid;
+ _recalc = true;
+ // top and bottom quark masses
+ Energy mt = _theSM->mass(q2, _top);
+ Energy mb = _theSM->mass(q2, _bt);
+ Energy mtau = _theSM->mass(q2, _tau);
+ // scalar
+ if( hid % 5 == 0 ) {
+ // location of the higgs
+ int iloc = (hid - 25)/10;
+ // 6 particles in the loop
+ setNParticles(13);
+ Complex c(0.);
+ // couplings for the top quark loop
+ c = -1.5*sqr(_theSM->eu())* mt*(*_mixS)(iloc, 1)/_sb/_mw;
+ couplings[0] = make_pair(c,c);
+ masses[0] = mt;
+ // couplings for the bottom quark loop
+ c = -1.5*sqr(_theSM->ed())* mb*(*_mixS)(iloc, 0)/_cb/_mw;
+ couplings[1] = make_pair(c,c);
+ masses[1] = mb;
+ // couplings for the tau lepton loop
+ c = -0.5*sqr(_theSM->ee())*mtau*(*_mixS)(iloc, 0)/_cb/_mw;
+ couplings[2] = make_pair(c,c);
+ masses[2] = mtau;
+ // charginos
+ for(unsigned int ic=0;ic<2;++ic) {
+ c = -_lambda/_coup*rt*(*_mixS)(iloc,2)*(*_mixU)(ic,1)*(*_mixV)(ic,1)
+ -rt*((*_mixS)(iloc,0)*(*_mixU)(ic,1)*(*_mixV)(ic,0) +
+ (*_mixS)(iloc,1)*(*_mixU)(ic,0)*(*_mixV)(ic,1));
+ couplings[3+ic] = make_pair(c,c);
+ }
+ // W boson
+ c = UnitRemoval::InvE*_mw*
+ (_cb*(*_mixS)(iloc,0)+_sb*(*_mixS)(iloc,1));
+ couplings[5] = make_pair(c,c);
+ // charged Higgs
+ complex<Energy> cpl;
+ cpl = sqr(_lambda)*rt*2.*(_s*((*_mixS)(iloc,2)*sqr(_cb) + (*_mixS)(iloc,2)*sqr(_sb))
+ - (_vu*(*_mixS)(iloc,0)/_coup +
+ _vd*(*_mixS)(iloc,1)/_coup)*_sb*_cb)
+ +_lambda*_sb*_cb*2*(*_mixS)(iloc,2)*(_kappa*_s/rt + rt*_theAl)
+ + sqr(_coup)*0.5*rt*sqr(_sw)/sqr(_cw)*((_vu*(*_mixS)(iloc,1)/_coup -
+ _vd*(*_mixS)(iloc,0)/_coup)*sqr(_cb) +
+ (_vd*(*_mixS)(iloc,0)/_coup -
+ _vu*(*_mixS)(iloc,1)/_coup)*sqr(_sb))
+ + sqr(_coup)*0.5*rt*(_vu*((*_mixS)(iloc,1)*sqr(_cb) +
+ (*_mixS)(iloc,1)*sqr(_sb) +
+ 2.*(*_mixS)(iloc,0)*_cb*_sb)/_coup
+ + _vd*((*_mixS)(iloc,0)*sqr(_cb) +
+ (*_mixS)(iloc,0)*sqr(_sb)
+ + 2.*(*_mixS)(iloc,1)*_sb*_cb)/_coup);
+ cpl /= -_coup;
+ couplings[6] = make_pair(cpl*UnitRemoval::InvE,cpl*UnitRemoval::InvE);
+ // sbottoms
+ double f1 = mb/_mw/_cb;
+ complex<Energy> f2 = 0.5*_mz/_cw*
+ ( - _cb*(*_mixS)(iloc,0) + _sb*(*_mixS)(iloc,1));
+ for(unsigned int ix=0;ix<2;++ix) {
+ cpl = -f2*( (1. - 2.*sqr(_sw)/3.)*(*_mixQb)(ix, 0)*(*_mixQb)(ix, 0)
+ + 2.*sqr(_sw)*(*_mixQb)(ix, 1)*(*_mixQb)(ix, 1)/3.)
+ - f1*mb*(*_mixS)(iloc,0)
+ *((*_mixQb)(ix, 0)*(*_mixQb)(ix, 0) + (*_mixQb)(ix, 1)*(*_mixQb)(ix, 1))
+ - 0.5*f1*(-_lambdaVEV*(*_mixS)(iloc,1) - _lambda*_vu*(*_mixS)(iloc,2)/_coup
+ + _triBt*(*_mixS)(iloc,0))*((*_mixQb)(ix, 1)*(*_mixQb)(ix, 0)
+ + (*_mixQb)(ix, 0)*(*_mixQb)(ix, 1));
+ cpl *= 3.*sqr(_theSM->ed());
+ couplings[7+ix] = make_pair(cpl*UnitRemoval::InvE,cpl*UnitRemoval::InvE);
+ }
+ // stop
+ f1 = mt/_mw/_sb;
+ for(unsigned int ix=0;ix<2;++ix) {
+ cpl =+f2*( (1. - 4.*sqr(_sw)/3.)*(*_mixQt)(ix, 0)*(*_mixQt)(ix, 0)
+ + 4.*sqr(_sw)*(*_mixQt)(ix, 1)*(*_mixQt)(ix, 1)/3.)
+ - f1*mt*(*_mixS)(iloc,1)
+ *((*_mixQt)(ix, 0)*(*_mixQt)(ix, 0)
+ + (*_mixQt)(ix, 1)*(*_mixQt)(ix, 1))
+ - 0.5*f1*(-_lambdaVEV*(*_mixS)(iloc,0) - _lambda*_vd*(*_mixS)(iloc,2)/_coup
+ + _triTp*(*_mixS)(iloc,1))*((*_mixQt)(ix, 1)*(*_mixQt)(ix, 0)
+ + (*_mixQt)(ix, 0)*(*_mixQt)(ix, 1));
+ cpl *= 3.*sqr(_theSM->eu());
+ couplings[9+ix] = make_pair(cpl*UnitRemoval::InvE,cpl*UnitRemoval::InvE);
+ } // sbottoms
+ f1 = mtau/_mw/_cb;
+ for(unsigned int ix=0;ix<2;++ix) {
+ cpl = -f2*( (1. - 2.*sqr(_sw))*(*_mixLt)(ix, 0)*(*_mixLt)(ix, 0)
+ + 2.*sqr(_sw)*(*_mixLt)(ix, 1)*(*_mixLt)(ix, 1))
+ - f1*mtau*(*_mixS)(iloc,0)
+ *((*_mixLt)(ix, 0)*(*_mixLt)(ix, 0) + (*_mixLt)(ix, 1)*(*_mixLt)(ix, 1))
+ - 0.5*f1*(-_lambdaVEV*(*_mixS)(iloc,1) - _lambda*_vu*(*_mixS)(iloc,2)/_coup
+ + _triTa*(*_mixS)(iloc,0))*((*_mixLt)(ix, 1)*(*_mixLt)(ix, 0)
+ + (*_mixLt)(ix, 0)*(*_mixLt)(ix, 1));
+ cpl *= sqr(_theSM->ee());
+ couplings[11+ix] = make_pair(cpl*UnitRemoval::InvE,cpl*UnitRemoval::InvE);
+ }
+ }
+ // pseudoscalar higgs bosons
+ else {
+ // location of the higgs
+ int iloc = (hid - 36)/10;
+ // 2 particles in the loop
+ setNParticles(5);
+ Complex c(0.);
+ // top quark couplings
+ c = Complex(0., 1.)*1.5*sqr(_theSM->eu())* mt*(*_mixP)(iloc, 1)/_sb/_mw;
+ couplings[0] = make_pair(c,-c);
+ masses[0] = mt;
+ // bottom quark couplings
+ c = Complex(0., 1.)*1.5*sqr(_theSM->ed())* mb*(*_mixP)(iloc, 0)/_cb/_mw;
+ couplings[1] = make_pair(c,-c);
+ masses[1] = mb;
+ // tau lepton couplings
+ c = Complex(0., 1.)*0.5*sqr(_theSM->ee())*mtau*(*_mixP)(iloc, 0)/_cb/_mw;
+ couplings[2] = make_pair(c,-c);
+ masses[2] = mtau;
+ // charginos
+ for(unsigned int ic=0;ic<2;++ic) {
+ c = Complex(0,-1.0)*
+ (_lambda/_coup*rt*(*_mixP)(iloc,2)*(*_mixU)(ic,1)*(*_mixV)(ic,1)
+ -rt*((*_mixP)(iloc,0)*(*_mixU)(ic,1)*(*_mixV)(ic,0)
+ + (*_mixP)(iloc,1)*(*_mixU)(ic,0)*(*_mixV)(ic,1)));
+ couplings[3+ic] = make_pair(-c,c);
+
+ }
+ }
+ }
+
+ if( _recalc ) {
+ VVSLoopVertex::setCoupling(q2, p1, p2, p3);
+ _recalc = false;
+ }
+}
diff --git a/Models/Susy/NMSSM/NMSSMPPHVertex.h b/Models/Susy/NMSSM/NMSSMPPHVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMPPHVertex.h
@@ -0,0 +1,341 @@
+// -*- C++ -*-
+#ifndef HERWIG_NMSSMPPHVertex_H
+#define HERWIG_NMSSMPPHVertex_H
+//
+// This is the declaration of the NMSSMPPHVertex class.
+//
+
+#include "Herwig++/Models/General/VVSLoopVertex.h"
+#include "Herwig++/Models/StandardModel/StandardModel.fh"
+#include "Herwig++/Models/Susy/MixingMatrix.h"
+
+namespace Herwig {
+using namespace ThePEG;
+
+/**
+ * This class implements the effective vertex for a higgs coupling
+ * to a pair of photons in the NMSSM.
+ *
+ * @see \ref NMSSMPPHVertexInterfaces "The interfaces"
+ * defined for NMSSMPPHVertex.
+ */
+class NMSSMPPHVertex: public VVSLoopVertex {
+
+public:
+
+ /** @name Standard constructors and destructors. */
+ //@{
+ /**
+ * The default constructor.
+ */
+ NMSSMPPHVertex();
+ //@}
+
+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();
+
+ /**
+ * Calculate couplings
+ *@param q2 Scale at which to evaluate coupling
+ *@param p1 ParticleData pointer to first particle
+ *@param p2 ParticleData pointer to second particle
+ *@param p3 ParticleData pointer to third particle
+ */
+ virtual void setCoupling(Energy2 q2, tcPDPtr p1, tcPDPtr p2,
+ tcPDPtr p3);
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const {return new_ptr(*this);}
+
+ /** 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 {return new_ptr(*this);}
+ //@}
+
+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:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is a concrete class with persistent data.
+ */
+ static ClassDescription<NMSSMPPHVertex> initNMSSMPPHVertex;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ NMSSMPPHVertex & operator=(const NMSSMPPHVertex &);
+
+private:
+
+ /** @name Stored parameters. */
+ //@{
+ /**
+ * The SM pointer
+ */
+ tcHwSMPtr _theSM;
+
+ /**
+ * \f$ \sin\theta_W\f$
+ */
+ double _sw;
+
+ /**
+ * \f$ \cos\theta_W\f$
+ */
+ double _cw;
+
+ /**
+ * \f$ M_W\f$
+ */
+ Energy _mw;
+
+ /**
+ * \f$ M_Z \f$
+ */
+ Energy _mz;
+
+ /**
+ * The product \f$\lambda \langle S\rangle \f$.
+ */
+
+ Energy _lambdaVEV;
+
+ /**
+ * The coefficient of the trilinear \f$SH_2 H_1\f$ term in the superpotential
+ */
+
+ double _lambda;
+
+ /**
+ * The value of the VEV of the higgs that couples to the up-type sector
+ * \f$ g*sqrt(2)M_W\cos\beta \f$
+ */
+
+ Energy _v1;
+
+ /**
+ * The value of the VEV of the higgs that couples to the down-type sector
+ * \f$ g*sqrt(2)M_W\cos\beta \f$
+ */
+ Energy _v2;
+
+ /**
+ * The top quark trilinear coupling
+ */
+ complex<Energy> _triTp;
+
+ /**
+ * The bottom quark trilinear coupling
+ */
+ complex<Energy> _triBt;
+
+ /**
+ * The tau quark trilinear coupling
+ */
+ complex<Energy> _triTa;
+
+ /**
+ * A pointer to the top quark
+ */
+ tcPDPtr _top;
+
+ /**
+ * A pointer to the bottom quark
+ */
+ tcPDPtr _bt;
+
+ /**
+ * A pointer to the tau lepton
+ */
+ tcPDPtr _tau;
+
+ /**
+ * CP-even Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixS;
+
+ /**
+ * CP-even Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixP;
+
+ /**
+ * \f$\tilde{t}\f$ mixing matrix
+ */
+ MixingMatrixPtr _mixQt;
+
+ /**
+ * \f$\tilde{b}\f$ mixing matrix
+ */
+ MixingMatrixPtr _mixQb;
+
+ /**
+ * \f$\tilde{\tau}\f$ mixing matrix
+ */
+ MixingMatrixPtr _mixLt;
+
+ /**
+ * \f$ \sin\beta\f$
+ */
+ double _sb;
+
+ /**
+ * \f$ \cos\beta\f$
+ */
+ double _cb;
+
+ /**
+ * The coefficient of the cubic singlet term in the superpotential
+ */
+ double _kappa;
+
+ /**
+ * The value of the VEV of the higgs that couples to the down-type sector
+ * \f$ g*sqrt(2)M_W\cos\beta \f$
+ */
+ Energy _vu;
+
+ /**
+ * The value of the VEV of the higgs that couples to the up-type sector
+ * i.e. \f$ g*sqrt(2)M_W\sin\beta \f$
+ */
+ Energy _vd;
+
+ /**
+ * The value of the VEV of the singlet higgs
+ */
+ Energy _s;
+
+ /**
+ * The soft trilinear \f$SH_2 H_1\f$ coupling
+ */
+ Energy _theAl;
+
+ /**
+ * The U mixing matrix
+ */
+ tMixingMatrixPtr _mixU;
+
+ /**
+ * The V mixing matrix
+ */
+ tMixingMatrixPtr _mixV;
+
+ /**
+ * The top and bottom quark masses calculated at the last value
+ * of \f$q^2\f$
+ */
+ pair<Energy, Energy> _masslast;
+
+ /**
+ * The scale at which the coupling was last evaluated
+ */
+ Energy2 _q2last;
+
+ /**
+ * The value of the overall normalisation when the coupling was last
+ * evaluated.
+ */
+ double _couplast;
+
+ /**
+ * The value of the weak coupling was last
+ * evaluated.
+ */
+ double _coup;
+
+ /**
+ * The PDG code of the Higgs particle when the vertex was last evaluated
+ */
+ long _hlast;
+
+ /**
+ * Whether the tensor coefficient need recalculating or not
+ */
+ bool _recalc;
+ //@}
+};
+
+}
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of NMSSMPPHVertex. */
+template <>
+struct BaseClassTrait<Herwig::NMSSMPPHVertex,1> {
+ /** Typedef of the first base class of NMSSMPPHVertex. */
+ typedef Herwig::VVSLoopVertex NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the NMSSMPPHVertex class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::NMSSMPPHVertex>
+ : public ClassTraitsBase<Herwig::NMSSMPPHVertex> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::NMSSMPPHVertex"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * NMSSMPPHVertex is implemented. It may also include several, space-separated,
+ * libraries if the class NMSSMPPHVertex 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 "HwSusy.so HwNMSSM.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_NMSSMPPHVertex_H */
diff --git a/Models/Susy/NMSSM/NMSSMWHHVertex.cc b/Models/Susy/NMSSM/NMSSMWHHVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMWHHVertex.cc
@@ -0,0 +1,156 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the NMSSMWHHVertex class.
+//
+
+#include "NMSSMWHHVertex.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "NMSSM.h"
+
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+NMSSMWHHVertex::NMSSMWHHVertex() : _sinb(0.), _cosb(0.), _sw(0.), _cw(0.),
+ _q2last(0.*MeV2), _couplast(0.) {
+ // codes for the neutral higgs
+ //CP even
+ int ieven[3]={25,35,45};
+ //CP odd
+ int iodd [2]={36,46};
+ // Z CP even CP odd
+ for(unsigned int ix=0;ix<3;++ix)
+ for(unsigned int iy=0;iy<2;++iy)
+ addToList( 23, ieven[ix], iodd[iy] );
+
+ // W H+ CP even
+ for(unsigned int ix=0;ix<3;++ix)
+ addToList( -24, 37, ieven[ix] );
+
+ // W+ H- CP even
+ for(unsigned int ix=0;ix<3;++ix)
+ addToList( 24, -37, ieven[ix] );
+
+ // W H+ CP odd
+ for(unsigned int ix=0;ix<2;++ix)
+ addToList( -24, 37, iodd[ix] );
+
+ //W+ H- CP odd
+ for(unsigned int ix=0;ix<2;++ix)
+ addToList( 24, -37, iodd[ix] );
+
+ // Charged higgs Z/gamma
+ addToList( 22, 37, -37 );
+ addToList( 23, 37, -37 );
+}
+
+void NMSSMWHHVertex::doinit() {
+ // cast to NMSSM model
+ tcNMSSMPtr model=dynamic_ptr_cast<tcNMSSMPtr>(generator()->standardModel());
+ if(!model)
+ throw InitException() << "Must have the NMSSM Model in NMSSMFFHVertex::doinit()"
+ << Exception::runerror;
+ // sin theta_W
+ double sw2 = sin2ThetaW();
+ _sw = sqrt(sw2);
+ _cw = sqrt(1.-sw2);
+ // get the mixing matrices
+ _mixS=model->CPevenHiggsMix();
+ if(!_mixS) throw InitException() << "Mixing matrix for CP-even neutral Higgs"
+ << " bosons is not set in NMSSMWHHVertex::doinit()"
+ << Exception::runerror;
+ _mixP=model->CPoddHiggsMix();
+ if(!_mixP) throw InitException() << "Mixing matrix for CP-odd neutral Higgs"
+ << " bosons is not set in NMSSMWHHVertex::doinit()"
+ << Exception::runerror;
+ // sin and cos beta
+ double beta = atan(model->tanBeta());
+ _sinb = sin(beta);
+ _cosb = cos(beta);
+ // order in the couplings
+ orderInGem(1);
+ orderInGs(0);
+ // base class
+ VSSVertex::doinit();
+}
+
+void NMSSMWHHVertex::persistentOutput(PersistentOStream & os) const {
+ os << _sinb << _cosb << _sw << _cw << _mixS << _mixP;
+}
+
+void NMSSMWHHVertex::persistentInput(PersistentIStream & is, int) {
+ is >> _sinb >> _cosb >> _sw >> _cw >> _mixS >> _mixP;
+}
+
+ClassDescription<NMSSMWHHVertex> NMSSMWHHVertex::initNMSSMWHHVertex;
+// Definition of the static class description member.
+
+void NMSSMWHHVertex::Init() {
+
+ static ClassDocumentation<NMSSMWHHVertex> documentation
+ ("The NMSSMWHHVertex class implements the coupling of an electroweak"
+ " gauge boson with two Higgs bosons in the NMSSM.");
+
+}
+
+//calulate the couplings
+void NMSSMWHHVertex::setCoupling(Energy2 q2,tcPDPtr a,tcPDPtr b,tcPDPtr c) {
+ // weak coupling
+ if(q2!=_q2last) {
+ _couplast = weakCoupling(q2);
+ _q2last=q2;
+ }
+ // gauge bosons
+ int ibos= a->id();
+ int ih1 = b->id();
+ int ih2 = c->id();
+ Complex fact;
+ if(ibos==ParticleID::Z0) {
+ fact = 0.5/_cw;
+ // Z H+ H-
+ if(abs(ih1)==37) {
+ fact *= (sqr(_cw)-sqr(_sw));
+ if(ih1<0) fact *=-1.;
+ }
+ // Z CP even CP odd
+ else {
+ fact *= -1.;
+ if(ih1%10==6) {
+ fact *= -1.;
+ swap(ih1,ih2);
+ }
+ int is = (ih1-25)/10;
+ int ip = (ih2-36)/10;
+ fact *= Complex(0.,1.)*((*_mixS)(is,1)*(*_mixP)(ip,1)-
+ (*_mixS)(is,0)*(*_mixP)(ip,0));
+ }
+ }
+ // gamma CP even CP odd
+ else if(ibos==ParticleID::gamma) {
+ fact = ih1>0 ? _sw : -_sw;
+ }
+ // W boson
+ else {
+ fact = 0.5;
+ if(abs(ih2)==37) {
+ swap(ih1,ih2);
+ fact*=-1;
+ }
+ if(ibos<0&&ih2%5==0) fact*=-1;
+ // H+ CP even
+ if(ih2%5==0) {
+ int is = (ih2-25)/10;
+ fact *= (_cosb*(*_mixS)(is,1)-_sinb*(*_mixS)(is,0));
+ }
+ // H+ CP odd
+ else {
+ int ip = (ih2-36)/10;
+ fact *= Complex(0.,1.)*(_cosb*(*_mixP)(ip,1)+_sinb*(*_mixP)(ip,0));
+ }
+ }
+ //output the coupling
+ norm(_couplast*fact);
+}
diff --git a/Models/Susy/NMSSM/NMSSMWHHVertex.h b/Models/Susy/NMSSM/NMSSMWHHVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMWHHVertex.h
@@ -0,0 +1,189 @@
+// -*- C++ -*-
+#ifndef HERWIG_NMSSMWHHVertex_H
+#define HERWIG_NMSSMWHHVertex_H
+//
+// This is the declaration of the NMSSMWHHVertex class.
+//
+
+#include "ThePEG/Helicity/Vertex/Scalar/VSSVertex.h"
+#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/Susy/MixingMatrix.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * The NMSSMWHHVertex class implements the coupling of an electroweak"
+ * gauge boson with two Higgs bosons in the NMSSM.
+ *
+ * @see \ref NMSSMWHHVertexInterfaces "The interfaces"
+ * defined for NMSSMWHHVertex.
+ */
+class NMSSMWHHVertex: public VSSVertex {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ NMSSMWHHVertex();
+
+ /** @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();
+
+ /**
+ * Calculate the couplings. This method is virtual and must be implemented in
+ * classes inheriting from this.
+ * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
+ * @param part1 The ParticleData pointer for the first particle.
+ * @param part2 The ParticleData pointer for the second particle.
+ * @param part3 The ParticleData pointer for the third particle.
+ */
+ virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,tcPDPtr part3);
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const {return new_ptr(*this);}
+
+ /** 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 {return new_ptr(*this);}
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is a concrete class with persistent data.
+ */
+ static ClassDescription<NMSSMWHHVertex> initNMSSMWHHVertex;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ NMSSMWHHVertex & operator=(const NMSSMWHHVertex &);
+
+private:
+
+ /**
+ * \f$\sin\beta\f$
+ */
+ double _sinb;
+
+ /**
+ * \f$\cos\beta\f$
+ */
+ double _cosb;
+
+ /**
+ * \f$\sin\theta_W\f$
+ */
+ double _sw;
+
+ /**
+ * \f$\cos\theta_W\f$
+ */
+ double _cw;
+
+ /**
+ * The last \f$q^2\f$ the coupling was evaluated at.
+ */
+ Energy2 _q2last;
+
+ /**
+ * The last value of the coupling
+ */
+ double _couplast;
+
+ /**
+ * Mixing matrix for the CP-even Higgs bosons
+ */
+ MixingMatrixPtr _mixS;
+
+ /**
+ * Mixing matrix for the CP-odd Higgs bosons
+ */
+ MixingMatrixPtr _mixP;
+
+};
+}
+
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of NMSSMWHHVertex. */
+template <>
+struct BaseClassTrait<Herwig::NMSSMWHHVertex,1> {
+ /** Typedef of the first base class of NMSSMWHHVertex. */
+ typedef ThePEG::Helicity::VSSVertex NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the NMSSMWHHVertex class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::NMSSMWHHVertex>
+ : public ClassTraitsBase<Herwig::NMSSMWHHVertex> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::NMSSMWHHVertex"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * NMSSMWHHVertex is implemented. It may also include several, space-separated,
+ * libraries if the class NMSSMWHHVertex 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 "HwSusy.so HwNMSSM.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_NMSSMWHHVertex_H */
diff --git a/Models/Susy/NMSSM/NMSSMWWHVertex.cc b/Models/Susy/NMSSM/NMSSMWWHVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMWWHVertex.cc
@@ -0,0 +1,91 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the NMSSMWWHVertex class.
+//
+
+#include "NMSSMWWHVertex.h"
+#include "ThePEG/PDT/EnumParticles.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "NMSSM.h"
+
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+NMSSMWWHVertex::NMSSMWWHVertex()
+ : _couplast(0.), _q2last(), _mw(), _zfact(0.), _sinb(0.),_cosb(0.) {
+ int id[3]={25,35,45};
+ // PDG codes for the particles in the vertex
+ for(unsigned int ix=0;ix<3;++ix) {
+ // Higgs WW
+ addToList( 24, -24, id[ix] );
+ //Higgs ZZ
+ addToList( 23, 23, id[ix] );
+ }
+}
+
+void NMSSMWWHVertex::doinit() {
+ // SM parameters
+ _mw = getParticleData(ThePEG::ParticleID::Wplus)->mass();
+ double sw = sin2ThetaW();
+ _zfact = 1./(1.-sw);
+ sw = sqrt(sw);
+ // NMSSM parameters
+ tcNMSSMPtr model=dynamic_ptr_cast<tcNMSSMPtr>(generator()->standardModel());
+ if(!model)
+ throw InitException() << "Must have the NMSSM Model in NMSSMWWHVertex::doinit()"
+ << Exception::runerror;
+ // get the mixing matrices
+ _mixS=model->CPevenHiggsMix();
+ if(!_mixS) throw InitException() << "Mixing matrix for CP-even neutral Higgs"
+ << " bosons is not set in NMSSMWWHVertex::doinit()"
+ << Exception::runerror;
+ // sin and cos beta
+ double beta = atan(model->tanBeta());
+ _sinb=sin(beta);
+ _cosb=cos(beta);
+ // order in the couplings
+ orderInGem(1);
+ orderInGs(0);
+ // base class
+ VVSVertex::doinit();
+}
+
+void NMSSMWWHVertex::persistentOutput(PersistentOStream & os) const {
+ os << ounit(_mw,GeV) << _zfact << _sinb << _cosb << _mixS;
+}
+
+void NMSSMWWHVertex::persistentInput(PersistentIStream & is, int) {
+ is >> iunit(_mw,GeV) >> _zfact >> _sinb >> _cosb >> _mixS;
+}
+
+ClassDescription<NMSSMWWHVertex> NMSSMWWHVertex::initNMSSMWWHVertex;
+// Definition of the static class description member.
+
+void NMSSMWWHVertex::Init() {
+
+ static ClassDocumentation<NMSSMWWHVertex> documentation
+ ("The NMSSMWWHVertex class implements the coupling of two electroweak gauge"
+ " bosons with the Higgs bosons of the NMSSM");
+
+}
+//calulate couplings
+void NMSSMWWHVertex::setCoupling(Energy2 q2,tcPDPtr a,tcPDPtr, tcPDPtr c) {
+ // ID of gauge bosons
+ int ibos=abs(a->id());
+ // ID of Higgs
+ int ihiggs = (c->id()-25)/10;
+ // first the overall normalisation
+ if(q2!=_q2last) {
+ _couplast = weakCoupling(q2)*_mw*UnitRemoval::InvE;
+ _q2last=q2;
+ }
+ // higgs mixing factor
+ Complex hmix = _cosb*(*_mixS)(ihiggs,0)+_sinb*(*_mixS)(ihiggs,1);
+ // couplings
+ if(ibos==24) norm(_couplast*hmix);
+ else if(ibos==23) norm(_couplast*hmix*_zfact);
+ else assert(false);
+}
diff --git a/Models/Susy/NMSSM/NMSSMWWHVertex.h b/Models/Susy/NMSSM/NMSSMWWHVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/NMSSM/NMSSMWWHVertex.h
@@ -0,0 +1,193 @@
+// -*- C++ -*-
+#ifndef HERWIG_NMSSMWWHVertex_H
+#define HERWIG_NMSSMWWHVertex_H
+//
+// This is the declaration of the NMSSMWWHVertex class.
+//
+
+#include "ThePEG/Helicity/Vertex/Scalar/VVSVertex.h"
+#include "Herwig++/Models/Susy/MixingMatrix.h"
+#include "Herwig++/Models/StandardModel/StandardModel.h"
+
+namespace Herwig {
+using namespace ThePEG;
+
+/** \ingroup Helicity
+ *
+ * The NMSSMWWHVertex class is the implementation of the coupling of two electroweak
+ * gauge bosons to the Higgs bosons of the NMSSM. It inherits from VVSVertex and
+ * implements the setCoupling member.
+ *
+ * @see \ref NMSSMWWHVertexInterfaces "The interfaces"
+ * @see VVSVertex
+ * @see SMWWHVertex
+ * defined for NMSSMWWHVertex.
+ */
+class NMSSMWWHVertex: public VVSVertex {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ NMSSMWWHVertex();
+
+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();
+
+ /**
+ * Calculate the couplings.
+ * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
+ * @param part1 The ParticleData pointer for the first particle.
+ * @param part2 The ParticleData pointer for the second particle.
+ * @param part3 The ParticleData pointer for the third particle.
+ */
+ virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,tcPDPtr part3);
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const {return new_ptr(*this);}
+
+ /** 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 {return new_ptr(*this);}
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The static object used to initialize the description of this class.
+ * Indicates that this is a concrete class with persistent data.
+ */
+ static ClassDescription<NMSSMWWHVertex> initNMSSMWWHVertex;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ NMSSMWWHVertex & operator=(const NMSSMWWHVertex &);
+
+private:
+
+ /**
+ * Storage of the couplings.
+ */
+ //@{
+ /**
+ * The last value of the electroweak coupling calculated.
+ */
+ Complex _couplast;
+
+ /**
+ * The scale \f$q^2\f$ at which the coupling was last evaluated.
+ */
+ Energy2 _q2last;
+
+ /**
+ * The mass of the \f$W\f$ boson.
+ */
+ Energy _mw;
+
+ /**
+ * The factor for the \f$Z\f$ vertex.
+ */
+ double _zfact;
+
+ /**
+ * \f$\sin\beta\f$
+ */
+ double _sinb;
+
+ /**
+ * \f$\cos\beta\f$
+ */
+ double _cosb;
+
+ /**
+ * Mixing matrix for the CP-even Higgs bosons
+ */
+ MixingMatrixPtr _mixS;
+ //@}
+
+};
+}
+
+
+#include "ThePEG/Utilities/ClassTraits.h"
+
+namespace ThePEG {
+
+/** @cond TRAITSPECIALIZATIONS */
+
+/** This template specialization informs ThePEG about the
+ * base classes of NMSSMWWHVertex. */
+template <>
+struct BaseClassTrait<Herwig::NMSSMWWHVertex,1> {
+ /** Typedef of the first base class of NMSSMWWHVertex. */
+ typedef ThePEG::Helicity::VVSVertex NthBase;
+};
+
+/** This template specialization informs ThePEG about the name of
+ * the NMSSMWWHVertex class and the shared object where it is defined. */
+template <>
+struct ClassTraits<Herwig::NMSSMWWHVertex>
+ : public ClassTraitsBase<Herwig::NMSSMWWHVertex> {
+ /** Return a platform-independent class name */
+ static string className() { return "Herwig::NMSSMWWHVertex"; }
+ /**
+ * The name of a file containing the dynamic library where the class
+ * NMSSMWWHVertex is implemented. It may also include several, space-separated,
+ * libraries if the class NMSSMWWHVertex 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 "HwSusy.so HwNMSSM.so"; }
+};
+
+/** @endcond */
+
+}
+
+#endif /* HERWIG_NMSSMWWHVertex_H */
diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -1,162 +1,163 @@
dnl Process this file with autoconf to produce a configure script.
AC_PREREQ([2.59])
AC_INIT([Herwig++],[SVN],[herwig@projects.hepforge.org],[Herwig++])
AC_CONFIG_SRCDIR([Utilities/HerwigStrategy.cc])
AC_CONFIG_AUX_DIR([Config])
AC_CONFIG_MACRO_DIR([m4])
AC_CONFIG_HEADERS([Config/config.h])
dnl AC_PRESERVE_HELP_ORDER
AC_CANONICAL_HOST
case "${host}" in
*-darwin[[0156]].*)
AC_MSG_ERROR([Herwig++ requires OS X 10.3 or later])
;;
*-darwin7.*)
if test "x$MACOSX_DEPLOYMENT_TARGET" != "x10.3"; then
AC_MSG_ERROR(
[Please add 'MACOSX_DEPLOYMENT_TARGET=10.3' to the configure line.])
fi
;;
esac
dnl === disable debug symbols by default =====
if test "x$CXXFLAGS" = "x"; then
CXXFLAGS=-O3
fi
if test "x$CFLAGS" = "x"; then
CFLAGS=-O3
fi
dnl Looptools manual requires optimization off
if test "x$FCFLAGS" = "x"; then
FCFLAGS=-O0
fi
dnl ==========================================
AC_LANG([C++])
AM_INIT_AUTOMAKE([1.9 gnu dist-bzip2 -Wall])
dnl Checks for programs.
AC_PROG_CXX([g++])
AC_PROG_INSTALL
AC_PROG_MAKE_SET
AC_PROG_LN_S
dnl modified search order
AC_PROG_FC([gfortran g95 g77])
dnl xlf95 f95 fort ifort ifc efc pgf95 lf95 ftn xlf90 f90 pgf90 pghpf epcf90 xlf f77 frt pgf77 cf77 fort77 fl32 af77])
AC_LANG_PUSH([Fortran])
AC_MSG_CHECKING([if the Fortran compiler ($FC) works])
AC_COMPILE_IFELSE(
AC_LANG_PROGRAM([],[ print *[,]"Hello"]),
[AC_MSG_RESULT([yes])],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([A Fortran compiler is required to build Herwig++.])
]
)
AC_LANG_POP([Fortran])
LT_PREREQ([2.2])
LT_INIT([disable-static dlopen pic-only])
dnl ####################################
dnl ####################################
dnl for Doc/fixinterfaces.pl
AC_PATH_PROG(PERL, perl)
HERWIG_CHECK_GSL
HERWIG_CHECK_THEPEG
HERWIG_COMPILERFLAGS
HERWIG_LOOPTOOLS
HERWIG_PDF_PATH
HERWIG_CHECK_FASTJET
HERWIG_VERSIONSTRING
HERWIG_CHECK_ABS_BUG
HERWIG_ENABLE_MODELS
SHARED_FLAG=-shared
AM_CONDITIONAL(NEED_APPLE_FIXES,
[test "xx${host/darwin/foundit}xx" != "xx${host}xx"])
if test "xx${host/darwin/foundit}xx" != "xx${host}xx"; then
APPLE_DSO_FLAGS=-Wl,-undefined,dynamic_lookup
SHARED_FLAG=-bundle
fi
AC_SUBST([APPLE_DSO_FLAGS])
AC_SUBST([SHARED_FLAG])
AC_CONFIG_FILES([UnderlyingEvent/Makefile
Models/Makefile
Models/StandardModel/Makefile
Models/RSModel/Makefile
Models/General/Makefile
Models/Susy/Makefile
+ Models/Susy/NMSSM/Makefile
Models/UED/Makefile
Models/Transplanckian/Makefile
Models/ADD/Makefile
Decay/Makefile
Decay/FormFactors/Makefile
Decay/Tau/Makefile
Decay/Baryon/Makefile
Decay/VectorMeson/Makefile
Decay/Perturbative/Makefile
Decay/ScalarMeson/Makefile
Decay/TensorMeson/Makefile
Decay/WeakCurrents/Makefile
Decay/Partonic/Makefile
Decay/General/Makefile
Decay/Radiation/Makefile
Doc/refman.conf
Doc/refman.h
PDT/Makefile
PDF/Makefile
MatrixElement/Makefile
MatrixElement/General/Makefile
MatrixElement/Lepton/Makefile
MatrixElement/Hadron/Makefile
MatrixElement/DIS/Makefile
MatrixElement/Powheg/Makefile
MatrixElement/Gamma/Makefile
Shower/SplittingFunctions/Makefile
Shower/Couplings/Makefile
Shower/Default/MECorrections/Makefile
Shower/Default/Makefile
Shower/Base/Makefile
Shower/Powheg/Makefile
Shower/Powheg/HardGenerators/Makefile
Shower/Makefile
Utilities/Makefile
Hadronization/Makefile
lib/Makefile
include/Makefile
src/Makefile
src/defaults/Makefile
src/herwig-config
Doc/Makefile
Doc/HerwigDefaults.in
Looptools/Makefile
Analysis/Makefile
src/Makefile-UserModules
src/defaults/Analysis.in
Contrib/Makefile
Contrib/make_makefiles.sh
Tests/Makefile
Makefile])
AC_CONFIG_FILES([Doc/fixinterfaces.pl],[chmod +x Doc/fixinterfaces.pl])
HERWIG_OVERVIEW
AC_CONFIG_COMMANDS([summary],[cat config.herwig])
AC_OUTPUT
diff --git a/m4/herwig.m4 b/m4/herwig.m4
--- a/m4/herwig.m4
+++ b/m4/herwig.m4
@@ -1,435 +1,440 @@
# check for gcc bug http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34130
AC_DEFUN([HERWIG_CHECK_ABS_BUG],
[
AC_REQUIRE([HERWIG_COMPILERFLAGS])
if test "$GCC" = "yes"; then
AC_MSG_CHECKING([for gcc abs bug])
AC_RUN_IFELSE(
AC_LANG_PROGRAM(
[ int foo (int i) { return -2 * __builtin_abs(i - 2); } ],
[ if ( foo(1) != -2 || foo(3) != -2 ) return 1; ]
),
[ AC_MSG_RESULT([not found. Compiler is ok.]) ],
[
AC_MSG_RESULT([found. Builtin abs() is buggy.])
AC_MSG_CHECKING([if -fno-builtin-abs works])
oldcxxflags=$CXXFLAGS
CXXFLAGS="$CXXFLAGS -fno-builtin-abs"
AC_RUN_IFELSE(
AC_LANG_PROGRAM(
[
#include <cstdlib>
int foo (int i) { return -2 * std::abs(i - 2); }
],
[
if (foo(1) != -2 || foo(3) != -2) return 1;
]
),
[
AC_MSG_RESULT([yes. Setting -fno-builtin-abs.])
AM_CXXFLAGS="$AM_CXXFLAGS -fno-builtin-abs"
AM_CFLAGS="$AM_CFLAGS -fno-builtin-abs"
],
[
AC_MSG_RESULT([no. Setting -fno-builtin.])
AC_MSG_WARN([
*****************************************************************************
For this version of gcc, -fno-builtin-abs alone did not work to avoid the
gcc abs() bug. Instead, all gcc builtin functions are now disabled.
Update gcc if possible.
*****************************************************************************])
AM_CXXFLAGS="$AM_CXXFLAGS -fno-builtin"
AM_CFLAGS="$AM_CFLAGS -fno-builtin"
]
)
CXXFLAGS=$oldcxxflags
]
)
fi
])
dnl ##### THEPEG #####
AC_DEFUN([HERWIG_CHECK_THEPEG],
[
defaultlocation="${prefix}"
test "x$defaultlocation" = xNONE && defaultlocation="${ac_default_prefix}"
AC_MSG_CHECKING([for libThePEG in])
AC_ARG_WITH(thepeg,
AC_HELP_STRING([--with-thepeg=DIR],[location of ThePEG installation]),
[],
[with_thepeg="${defaultlocation}"])
AC_MSG_RESULT([$with_thepeg])
if test "x$with_thepeg" = "xno"; then
AC_MSG_ERROR([Cannot build Herwig++ without ThePEG. Please set --with-thepeg.])
fi
THEPEGLDFLAGS="-L${with_thepeg}/lib/ThePEG"
THEPEGPATH="${with_thepeg}"
oldldflags="$LDFLAGS"
oldlibs="$LIBS"
LDFLAGS="$LDFLAGS $THEPEGLDFLAGS"
AC_CHECK_LIB([ThePEG],[debugThePEG],[],
[AC_MSG_ERROR([No ThePEG libraries in $THEPEGLDFLAGS. Please set --with-thepeg.])])
AC_SUBST([THEPEGLIB],[-lThePEG])
AC_SUBST(THEPEGLDFLAGS)
AC_SUBST(THEPEGPATH)
LIBS="$oldlibs"
LDFLAGS="$oldldflags"
AC_MSG_CHECKING([for ThePEG headers in])
AC_ARG_WITH([thepeg-headers],
AC_HELP_STRING([--with-thepeg-headers=DIR],[location of ThePEG include directory]),
[],
[with_thepeg_headers="${with_thepeg}/include"])
AC_MSG_RESULT([$with_thepeg_headers])
if test "x$with_thepeg_headers" = "xno"; then
AC_MSG_ERROR([Cannot build Herwig++ without ThePEG headers. Please set --with-thepeg-headers.])
fi
THEPEGINCLUDE="-I$with_thepeg_headers"
oldcppflags="$CPPFLAGS"
CPPFLAGS="$CPPFLAGS $THEPEGINCLUDE"
AC_CHECK_HEADER([ThePEG/Config/ThePEG.h],[],
[AC_MSG_ERROR([No ThePEG headers in $with_thepeg_headers. Please set --with-thepeg-headers.])])
CPPFLAGS="$oldcppflags"
AC_SUBST(THEPEGINCLUDE)
AC_MSG_CHECKING([for HepMCAnalysis.so in ThePEG])
if test -x "$THEPEGPATH/lib/ThePEG/HepMCAnalysis.so" ; then
CREATE_HEPMC="create"
AC_MSG_RESULT([found])
else
CREATE_HEPMC="# create"
AC_MSG_RESULT([not found])
fi
AC_SUBST([CREATE_HEPMC])
])
dnl ##### FastJet #####
AC_DEFUN([HERWIG_CHECK_FASTJET],[
FASTJETPATH=""
FASTJETLIBS=""
FASTJETINCLUDE=""
LOAD_FASTJET=""
CREATE_FASTJET="#create"
AC_MSG_CHECKING([for FastJet location])
AC_ARG_WITH(fastjet,
AC_HELP_STRING([--with-fastjet=DIR],[Location of FastJet installation @<:@default=system libs@:>@]),
[],
[with_fastjet=system])
if test "x$with_fastjet" = "xno"; then
AC_MSG_RESULT([FastJet support disabled.])
elif test "x$with_fastjet" = "xsystem"; then
AC_MSG_RESULT([in system libraries])
oldlibs="$LIBS"
AC_CHECK_LIB(fastjet,main,
[],
[with_fastjet=no
AC_MSG_WARN([
FastJet not found in system libraries])
])
FASTJETLIBS="$LIBS"
LIBS=$oldlibs
else
AC_MSG_RESULT([$with_fastjet])
FASTJETPATH="$with_fastjet"
FASTJETINCLUDE="-I$with_fastjet/include"
FASTJETLIBS="-L$with_fastjet/lib -R$with_fastjet/lib -lfastjet"
fi
if test "x$with_fastjet" != "xno"; then
# Now lets see if the libraries work properly
oldLIBS="$LIBS"
oldLDFLAGS="$LDFLAGS"
oldCPPFLAGS="$CPPFLAGS"
LIBS="$LIBS $FASTJETLIBS"
LDFLAGS="$LDFLAGS"
CPPFLAGS="$CPPFLAGS $FASTJETINCLUDE"
AC_MSG_CHECKING([that FastJet works])
AC_LINK_IFELSE([AC_LANG_PROGRAM([[#include <fastjet/ClusterSequence.hh>
]],[[fastjet::JetDefinition jet_def(fastjet::ee_kt_algorithm, fastjet::E_scheme, fastjet::Best); ]])],
[AC_MSG_RESULT([yes])],[AC_MSG_RESULT([no])
AC_MSG_ERROR([Use '--with-fastjet=' to set a path or use '--without-fastjet'.])
])
AC_CHECK_HEADERS([fastjet/ClusterSequence.hh])
LIBS="$oldLIBS"
LDFLAGS="$oldLDFLAGS"
CPPFLAGS="$oldCPPFLAGS"
CREATE_FASTJET="create"
LOAD_FASTJET="library HwLEPJetAnalysis.so"
fi
AC_SUBST(FASTJETINCLUDE)
AC_SUBST(CREATE_FASTJET)
AC_SUBST(LOAD_FASTJET)
AC_SUBST(FASTJETLIBS)
AM_CONDITIONAL(WANT_LIBFASTJET,[test "x$CREATE_FASTJET" = "xcreate"])
])
dnl ##### LOOPTOOLS #####
AC_DEFUN([HERWIG_LOOPTOOLS],
[
AC_REQUIRE([AC_PROG_FC])
AC_REQUIRE([AC_FC_LIBRARY_LDFLAGS])
AC_REQUIRE([AC_PROG_CC])
AC_REQUIRE([HERWIG_COMPILERFLAGS])
AC_MSG_CHECKING([if Looptools build works])
enable_looptools=yes
if test "x$GCC" = "xyes"; then
case "${host}" in
x86_64-*|*-darwin10.*)
AM_FCFLAGS="$AM_FCFLAGS -fdefault-integer-8"
;;
esac
AC_LANG_PUSH([Fortran])
oldFCFLAGS="$FCFLAGS"
FCFLAGS="$AM_FCFLAGS"
AC_COMPILE_IFELSE(
AC_LANG_PROGRAM([],[ print *[,]"Hello"]),
[],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([needs gfortran on 64bit machines])]
)
FCFLAGS="$oldFCFLAGS"
AC_LANG_POP([Fortran])
fi
AC_MSG_RESULT([$enable_looptools])
AC_SUBST([F77],[$FC])
AC_SUBST([FFLAGS],[$FCFLAGS])
AC_SUBST([AM_FFLAGS],[$AM_FCFLAGS])
AC_SUBST([FLIBS],[$FCLIBS])
])
dnl ##### PDF PATH #####
AC_DEFUN([HERWIG_PDF_PATH],
[
AC_MSG_CHECKING([which Herwig++ PDF data to use])
AC_ARG_WITH(pdf,
AC_HELP_STRING([--with-pdf=DIR],[installation path of Herwig++PDF data tarball]),
[],
[with_pdf=${prefix}]
)
HERWIG_PDF_PREFIX=${with_pdf}/share/Herwig++PDF/mrst
if test -f "${HERWIG_PDF_PREFIX}/2008/mrstMCal.dat"; then
AC_MSG_RESULT([$with_pdf])
localPDFneeded=false
else
AC_MSG_RESULT([Using built-in PDF data set. For other data sets, set --with-pdf.])
HERWIG_PDF_PREFIX=PDF/mrst
localPDFneeded=true
fi
HERWIG_PDF_DEFAULT=${HERWIG_PDF_PREFIX}/2008/mrstMCal.dat
HERWIG_PDF_NLO=${HERWIG_PDF_PREFIX}/2002/mrst2002nlo.dat
AM_CONDITIONAL(WANT_LOCAL_PDF,[test "x$localPDFneeded" = "xtrue"])
AC_SUBST(HERWIG_PDF_DEFAULT)
AC_SUBST(HERWIG_PDF_NLO)
])
dnl ###### GSL ######
AC_DEFUN([HERWIG_CHECK_GSL],
[
AC_MSG_CHECKING([for gsl location])
GSLINCLUDE=""
GSLLIBS=""
AC_ARG_WITH(gsl,
AC_HELP_STRING([--with-gsl=DIR],[location of gsl installation @<:@default=system libs@:>@]),
[],
[with_gsl=system])
if test "x$with_gsl" = "xno"; then
AC_MSG_ERROR([libgsl is required. Please install the GNU scientific library and header files.])
fi
if test "x$with_gsl" = "xsystem"; then
AC_MSG_RESULT([in system libraries])
oldlibs="$LIBS"
AC_CHECK_LIB(m,main)
AC_CHECK_LIB(gslcblas,main)
AC_CHECK_LIB(gsl,main,[],
[
AC_MSG_ERROR([Cannot find libgsl. Please install the GNU scientific library and header files or use --with-gsl=.])
]
)
GSLLIBS="$LIBS"
LIBS=$oldlibs
else
if test "`uname -m`" = "x86_64" -a -e "$with_gsl/lib64/libgsl.a" -a -d "$with_gsl/include/gsl"; then
AC_MSG_RESULT([found in $with_gsl])
GSLLIBS="-L$with_gsl/lib64 -R$with_gsl/lib64 -lgslcblas -lgsl"
GSLINCLUDE="-I$with_gsl/include"
elif test -e "$with_gsl/lib/libgsl.a" -a -d "$with_gsl/include/gsl"; then
AC_MSG_RESULT([found in $with_gsl])
GSLLIBS="-L$with_gsl/lib -R$with_gsl/lib -lgslcblas -lgsl"
GSLINCLUDE="-I$with_gsl/include"
else
AC_MSG_RESULT([not found])
AC_MSG_ERROR([Can't find $with_gsl/lib/libgsl.a or the headers in $with_gsl/include])
fi
fi
AC_SUBST(GSLINCLUDE)
AC_SUBST(GSLLIBS)
])
AC_DEFUN([HERWIG_VERSIONSTRING],
[
if test -d $srcdir/.svn; then
AC_CHECK_PROG(have_svnversion,[svnversion],[yes],[no])
fi
AM_CONDITIONAL(USE_SVNVERSION,[test "x$have_svnversion" = "xyes"])
])
dnl ##### COMPILERFLAGS #####
AC_DEFUN([HERWIG_COMPILERFLAGS],
[
AC_REQUIRE([HERWIG_CHECK_THEPEG])
AM_CPPFLAGS="-I\$(top_builddir)/include $THEPEGINCLUDE \$(GSLINCLUDE)"
AC_MSG_CHECKING([for debugging mode])
AC_ARG_ENABLE(debug,
AC_HELP_STRING([--enable-debug],[debug mode, use --enable-debug=slow for additional options that slow down the run.]),
[],
[enable_debug=no]
)
AC_MSG_RESULT([$enable_debug])
if test "x$enable_debug" = "xno"; then
AM_CPPFLAGS="$AM_CPPFLAGS -DNDEBUG"
else
debugflags="-g"
fi
dnl -Wfloat-equal -fvisibility-inlines-hidden -Wctor-dtor-privacy -Weffc++
if test -n "$GCC"; then
warnflags="-ansi -pedantic -Wall -W"
if test "x$enable_debug" = "xslow"; then
debugflags="$debugflags -fno-inline"
AM_CPPFLAGS="$AM_CPPFLAGS -D_GLIBCXX_DEBUG"
fi
fi
AC_SUBST(AM_CPPFLAGS)
AC_SUBST(AM_CFLAGS, ["$warnflags $debugflags"])
AC_SUBST(AM_CXXFLAGS,["$warnflags $debugflags"])
AC_SUBST(AM_FCFLAGS, ["$debugflags"])
AC_SUBST(AM_LDFLAGS)
])
AC_DEFUN([HERWIG_ENABLE_MODELS],
[
AC_MSG_CHECKING([for BSM models to include])
LOAD_RS=""
LOAD_SUSY=""
LOAD_TRP=""
LOAD_UED=""
LOAD_ADD=""
AC_ARG_ENABLE(models,
- AC_HELP_STRING([--enable-models=LIST],[Comma-separated list of BSM models to enable. Options are (mssm ued rs trp add) or --disable-models to turn them all off.]),
+ AC_HELP_STRING([--enable-models=LIST],[Comma-separated list of BSM models to enable. Options are (mssm nmssm ued rs trp add) or --disable-models to turn them all off.]),
[],
[enable_models=all]
)
if test "x$enable_models" = "xyes" -o "x$enable_models" = "xall"; then
all=yes
fi
AC_MSG_RESULT([$enable_models])
if test ! "$all"; then
oldIFS="$IFS"
IFS=","
for i in $enable_models; do
declare $i=yes
done
IFS="$oldIFS"
fi
+if test "$nmssm"; then
+ mssm=yes
+fi
+
if test "$rs" -o "$all" ; then
LOAD_RS="library HwRSModel.so"
fi
AC_SUBST(LOAD_RS)
if test "$mssm" -o "$all"; then
LOAD_SUSY="library HwSusy.so"
fi
AC_SUBST(LOAD_SUSY)
if test "$trp" -o "$all"; then
LOAD_TRP="library HwTransplanck.so"
fi
AC_SUBST(LOAD_TRP)
if test "$ued" -o "$all"; then
LOAD_UED="library HwUED.so"
fi
AC_SUBST(LOAD_UED)
if test "$add" -o "$all"; then
LOAD_ADD="library HwADDModel.so"
fi
AC_SUBST(LOAD_ADD)
AM_CONDITIONAL(WANT_MSSM,[test "$mssm" -o "$all"])
+AM_CONDITIONAL(WANT_NMSSM,[test "$nmssm" -o "$all"])
AM_CONDITIONAL(WANT_UED,[test "$ued" -o "$all"])
AM_CONDITIONAL(WANT_RS,[test "$rs" -o "$all"])
AM_CONDITIONAL(WANT_TRP,[test "$trp" -o "$all"])
AM_CONDITIONAL(WANT_ADD,[test "$add" -o "$all"])
])
AC_DEFUN([HERWIG_OVERVIEW],
[
FCSTRING=`$FC --version | head -1`
CXXSTRING=`$CXX --version | head -1`
CCSTRING=`$CC --version | head -1`
cat << _HW_EOF_ > config.herwig
*****************************************************
*** $PACKAGE_STRING configuration summary
*** Please include this information in bug reports!
***--------------------------------------------------
*** Prefix: $prefix
***
*** BSM models: $enable_models
*** Herwig debug mode: $enable_debug
***
*** GSL: $with_gsl
***
*** ThePEG: $with_thepeg
*** ThePEG headers: $with_thepeg_headers
***
*** Fastjet: $with_fastjet
***
*** Host: $host
*** CXX: $CXXSTRING
*** FC: $FCSTRING
*** CC: $CCSTRING
*****************************************************
_HW_EOF_
])
diff --git a/src/LHC-NMSSM.in b/src/LHC-NMSSM.in
new file mode 100644
--- /dev/null
+++ b/src/LHC-NMSSM.in
@@ -0,0 +1,60 @@
+##################################################
+# Example generator for the NMSSM
+# The best way to use this is to make your own
+# copy of this file and edit that as you require.
+#
+# The first section loads the model file which
+# does not contain anything that users need to touch.
+#
+# The second section contains the user settings.
+###################################################
+
+read NMSSM.model
+cd /Herwig/NewPhysics
+
+##################################################
+#
+# This section contains the user defined settings
+#
+##################################################
+# --- Hard Process ----
+# The particle name can be found in the relevant model file
+# by searching for its PDG code and noting the text
+# '/Herwig/Particles/###' where the hashes denote the name
+
+# Switch to decide whether to include EW diagrams in the
+# hard process (On by default)
+set HPConstructor:IncludeEW No
+
+# Example hard process: Incoming proton, outgoing squarks
+insert HPConstructor:Incoming 0 /Herwig/Particles/g
+insert HPConstructor:Incoming 1 /Herwig/Particles/u
+insert HPConstructor:Incoming 2 /Herwig/Particles/ubar
+insert HPConstructor:Incoming 3 /Herwig/Particles/d
+insert HPConstructor:Incoming 4 /Herwig/Particles/dbar
+
+insert HPConstructor:Outgoing 0 /Herwig/Particles/~u_L
+insert HPConstructor:Outgoing 1 /Herwig/Particles/~u_Lbar
+insert HPConstructor:Outgoing 2 /Herwig/Particles/~d_L
+insert HPConstructor:Outgoing 3 /Herwig/Particles/~d_Lbar
+
+# --- Perturbative Decays ---
+# Read in the spectrum file and optional decay table.
+# If a decay table is in a separate file
+# then add another 'setup' line with that
+# file as the argument. The provided
+# spectrum file is an example using NMHDecay-1.2.1
+setup NMSSM/Model NMSSM.spc
+
+# Intrinsic pT tune extrapolated to LHC energy
+set /Herwig/Shower/Evolver:IntrinsicPtGaussian 2.2*GeV
+
+# Other parameters for run
+cd /Herwig/Generators
+set LHCGenerator:NumberOfEvents 10000000
+set LHCGenerator:RandomNumberGenerator:Seed 31122001
+set LHCGenerator:DebugLevel 1
+set LHCGenerator:PrintEvent 10
+set LHCGenerator:MaxErrors 10000
+
+saverun LHC-NMSSM LHCGenerator
diff --git a/src/Makefile.am b/src/Makefile.am
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,118 +1,122 @@
SUBDIRS = defaults
AUTOMAKE_OPTIONS = -Wno-portability
defaultsdir = ${pkgdatadir}/defaults
bin_PROGRAMS = Herwig++
Herwig___SOURCES = Herwig++.cc HerwigRun.cc HerwigRun.h
Herwig___LDFLAGS = -export-dynamic $(THEPEGLDFLAGS)
Herwig___LDADD = $(THEPEGLIB) -ldl
Herwig___CPPFLAGS = $(AM_CPPFLAGS) \
-DHERWIG_PKGDATADIR="\"$(pkgdatadir)\"" \
-DHERWIG_PKGLIBDIR="\"$(pkglibdir)\"" \
-DTHEPEG_PKGLIBDIR="\"$(THEPEGPATH)/lib/ThePEG\""
bin_SCRIPTS = herwig-config
-HELPERFILES = SPhenoSPS1a.spc MUED.model MSSM.model RS.model ADD.model
+HELPERFILES = SPhenoSPS1a.spc NMSSM.spc MUED.model MSSM.model NMSSM.model RS.model ADD.model
-INPUTFILES = ILC.in LEP.in LHC.in TVT.in LHC-MSSM.in LHC-MUED.in LHC-RS.in LHC-ADD.in \
+INPUTFILES = ILC.in LEP.in LHC.in TVT.in LHC-MSSM.in \
+ LHC-MUED.in LHC-NMSSM.in LHC-RS.in LHC-ADD.in\
ILC-MSSM.in ILC-MUED.in ILC-RS.in DIS.in LHC-Powheg.in \
TVT-Powheg.in GammaGamma.in LHC-TRP.in
## provide -no3body for faster 'make check'
dist_noinst_DATA = LHC-MUED-no3body.in
dist_pkgdata_DATA = $(INPUTFILES) $(HELPERFILES)
pkgdata_DATA = Makefile-UserModules
CLEANFILES = HerwigDefaults.rpo \
*.run *.log *.out *.tex \
multi.test *.output probs.test chisq.value
## checking targets ##
HerwigDefaults.rpo: Herwig++ $(srcdir)/defaults/*.in defaults/PDF.in defaults/Analysis.in $(top_builddir)/lib/*.so
./Herwig++ init -L$(top_builddir)/lib -i defaults/HerwigDefaults.in --exitonerror
check_BSM_Full=
check_BSM=
if WANT_MSSM
-check_BSM += check-LHC-MSSM
+check_BSM += check-LHC-MSSM
check_BSM_Full += check-LHC-MSSM
endif
+if WANT_NMSSM
+check_BSM_Full += check-LHC-NMSSM
+endif
if WANT_UED
check_BSM += check-LHC-MUED-no3body
check_BSM_Full += check-LHC-MUED
endif
if WANT_RS
check_BSM += check-LHC-RS
check_BSM_Full += check-LHC-RS
endif
if WANT_TRP
check_BSM_Full += check-LHC-TRP
endif
if WANT_ADD
check_BSM_Full += check-LHC-ADD
endif
check-local: check-LHC check-LEP check-DIS check-ILC check-GammaGamma $(check_BSM) check-LHC-Powheg
check-Powheg: check-LHC-Powheg check-TVT-Powheg
check-BSM: $(check_BSM_Full)
link-helper-files:
@for i in $(HELPERFILES); do \
if test -f $(srcdir)/$$i -a ! -e $$i; then \
$(LN_S) -f $(srcdir)/$$i; fi; done
check-%: $(srcdir)/%.in HerwigDefaults.rpo link-helper-files
./Herwig++ read $< --exitonerror
./Herwig++ run $(notdir $(subst .in,.run,$<)) -N500 -d1 --exitonerror
## valgrind targets ##
VALGRIND=valgrind --leak-check=full --num-callers=25 --freelist-vol=100000000 --leak-resolution=med --trace-children=yes
valgrind: valgrind-init valgrind-read valgrind-run
valgrind-init:
$(VALGRIND) ./Herwig++ init -d 1 -L$(top_builddir)/lib -i defaults/HerwigDefaults.in \
&> /tmp/valgrind-init.log
valgrind-read:
$(VALGRIND) ./Herwig++ read -d 1 LHC.in &> /tmp/valgrind-read.log
valgrind-run:
$(VALGRIND) ./Herwig++ run -d 1 -N5 LHC.run &> /tmp/valgrind-run.log
SETUPTHEPEG=$(THEPEGPATH)/bin/setupThePEG
THEPEGREPO=$(THEPEGPATH)/lib/ThePEG/ThePEGDefaults.rpo
install-data-hook:
@echo Creating repository
@./Herwig++ init -L$(DESTDIR)$(pkglibdir) -i $(DESTDIR)$(defaultsdir)/HerwigDefaults.in -r $(DESTDIR)$(pkgdatadir)/HerwigDefaults.rpo --exitonerror
uninstall-hook:
rm -f $(DESTDIR)$(pkgdatadir)/HerwigDefaults.rpo
register: register-with-thepeg-repo
register-with-thepeg-repo:
@if test -x "$(SETUPTHEPEG)" -a -w "$(THEPEGREPO)"; \
then echo Registering with ThePEG; \
"$(SETUPTHEPEG)" --init \
$(DESTDIR)$(defaultsdir)/HerwigDefaults.in \
-r "$(THEPEGREPO)" -o "$(THEPEGREPO)" \
-l$(DESTDIR)$(pkglibdir) --exitonerror ; \
fi
unregister : unregister-from-thepeg-repo
unregister-from-thepeg-repo:
@if test -x "$(SETUPTHEPEG)" -a -w "$(THEPEGREPO)"; \
then echo Unregistering with ThePEG; \
"$(SETUPTHEPEG)" --init defaults/HerwigCleanup.in \
-r "$(THEPEGREPO)" -o "$(THEPEGREPO)" \
-l$(DESTDIR)$(pkglibdir) --exitonerror ; \
fi
include $(top_srcdir)/Utilities/Makefile.am.versionstring
diff --git a/src/NMSSM.model b/src/NMSSM.model
new file mode 100644
--- /dev/null
+++ b/src/NMSSM.model
@@ -0,0 +1,130 @@
+##################################################
+# Common setup for the NMSSM
+#
+# See LHC-NMSSM.in or ILC-NMSSM.in for example usage
+#
+# This file does not contain anything that
+# users need to touch. User settings are in
+# ???-NMSSM.in
+#
+###################################################
+
+###################################################
+#
+# Read the MSSM.
+#
+##################################################
+read MSSM.model
+
+###################################################
+#
+# Create the additional particle content in the
+# NMSSM
+#
+##################################################
+cd /Herwig/Particles
+create ThePEG::ParticleData H30
+setup H30 45 H30 500.0 0.0 0.0 0.0 0 0 1 0
+create ThePEG::ParticleData A20
+setup A20 46 A20 500.0 0.0 0.0 0.0 0 0 1 0
+create ThePEG::ParticleData ~chi_50
+setup ~chi_50 1000045 ~chi_50 500.0 0.0 0.0 0.0 0 0 2 0
+###################################################
+#
+# The main directory and model object
+#
+##################################################
+mkdir /Herwig/NewPhysics/NMSSM
+cd /Herwig/NewPhysics/NMSSM
+create Herwig::NMSSM Model HwNMSSM.so
+# SM couplings
+set Model:QCD/RunningAlphaS /Herwig/AlphaS
+set Model:EW/RunningAlphaEM /Herwig/AlphaEM
+set Model:EW/CKM /Herwig/CKM
+set Model:RunningMass /Herwig/RunningMass
+###################################################
+#
+# Vertices
+#
+##################################################
+# Create additional NMSSM vertices
+mkdir /Herwig/Vertices/NMSSM
+cd /Herwig/Vertices/NMSSM
+create Herwig::NMSSMFFHVertex NMSSM_FFH
+create Herwig::NMSSMGOGOHVertex NMSSM_GOGOH
+create Herwig::NMSSMHSFSFVertex NMSSM_HSS
+create Herwig::NMSSMWHHVertex NMSSM_WHH
+create Herwig::NMSSMWWHVertex NMSSM_WWH
+create Herwig::NMSSMHHHVertex NMSSM_HHH
+create Herwig::NMSSMGGHVertex NMSSM_GGH
+create Herwig::NMSSMPPHVertex NMSSM_PPH
+cd /Herwig/NewPhysics/NMSSM
+# SM vertices
+set Model:Vertex/FFZ /Herwig/Vertices/FFZVertex
+set Model:Vertex/FFW /Herwig/Vertices/FFWVertex
+set Model:Vertex/FFG /Herwig/Vertices/FFGVertex
+set Model:Vertex/FFP /Herwig/Vertices/FFPVertex
+set Model:Vertex/GGG /Herwig/Vertices/GGGVertex
+set Model:Vertex/GGGG /Herwig/Vertices/GGGGVertex
+set Model:Vertex/WWW /Herwig/Vertices/WWWVertex
+set Model:Vertex/WWWW /Herwig/Vertices/WWWWVertex
+set Model:Vertex/WWHH NULL
+set Model:Vertex/HHHH NULL
+# MSSM feynman rules
+set Model:Vertex/WSFSF /Herwig/Vertices/MSSM/MSSM_WSS
+set Model:Vertex/NFSF /Herwig/Vertices/MSSM/MSSM_NFS
+set Model:Vertex/GFSF /Herwig/Vertices/MSSM/MSSM_GFS
+set Model:Vertex/CFSF /Herwig/Vertices/MSSM/MSSM_CFS
+set Model:Vertex/GSFSF /Herwig/Vertices/MSSM/MSSM_GSS
+set Model:Vertex/GGSQSQ /Herwig/Vertices/MSSM/MSSM_GGSS
+set Model:Vertex/GSGSG /Herwig/Vertices/MSSM/MSSM_GGOGO
+set Model:Vertex/GNG NULL
+set Model:Vertex/NNZ /Herwig/Vertices/MSSM/MSSM_NNZ
+set Model:Vertex/NNP NULL
+set Model:Vertex/CCZ /Herwig/Vertices/MSSM/MSSM_CCZ
+set Model:Vertex/CNW /Herwig/Vertices/MSSM/MSSM_CNW
+set Model:Vertex/SSWHH /Herwig/Vertices/MSSM/MSSM_WHH
+set Model:Vertex/NCT /Herwig/Vertices/MSSM/MSSM_NCT
+# NMSSM feynman rules
+set Model:Vertex/FFH /Herwig/Vertices/NMSSM/NMSSM_FFH
+set Model:Vertex/GOGOH /Herwig/Vertices/NMSSM/NMSSM_GOGOH
+set Model:Vertex/HSFSF /Herwig/Vertices/NMSSM/NMSSM_HSS
+set Model:Vertex/SSWHH /Herwig/Vertices/NMSSM/NMSSM_WHH
+set Model:Vertex/WWH /Herwig/Vertices/NMSSM/NMSSM_WWH
+set Model:Vertex/HHH /Herwig/Vertices/NMSSM/NMSSM_HHH
+set Model:Vertex/HGG /Herwig/Vertices/NMSSM/NMSSM_GGH
+set Model:Vertex/HPP /Herwig/Vertices/NMSSM/NMSSM_PPH
+
+###################################################
+#
+# Set up spin correlation Decayers
+#
+###################################################
+cd /Herwig/NewPhysics
+
+set TwoBodyDC:CreateDecayModes Yes
+set ThreeBodyDC:CreateDecayModes Yes
+
+insert NewModel:DecayParticles 32 /Herwig/Particles/~chi_50
+insert NewModel:DecayParticles 33 /Herwig/Particles/H30
+insert NewModel:DecayParticles 34 /Herwig/Particles/A20
+
+insert /Herwig/Shower/ShowerHandler:DecayInShower 0 1000045 # SUSY_chi_50
+insert /Herwig/Shower/ShowerHandler:DecayInShower 0 45 # h30
+insert /Herwig/Shower/ShowerHandler:DecayInShower 0 46 # a20
+
+###################################################
+#
+# Generate processes and decays in the NMSSM
+#
+###################################################
+cd /Herwig/NewPhysics
+set NMSSM/Model:ModelGenerator NewModel
+###################################################
+#
+# Choose NMSSM over SM
+#
+###################################################
+cd /Herwig/Generators
+set LEPGenerator:StandardModelParameters /Herwig/NewPhysics/NMSSM/Model
+set LHCGenerator:StandardModelParameters /Herwig/NewPhysics/NMSSM/Model

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:25 AM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983021
Default Alt Text
(150 KB)

Event Timeline