diff --git a/Decay/ScalarMeson/Makefile.am b/Decay/ScalarMeson/Makefile.am
--- a/Decay/ScalarMeson/Makefile.am
+++ b/Decay/ScalarMeson/Makefile.am
@@ -1,42 +1,44 @@
 BUILT_SOURCES  = SMDecayer__all.cc
 CLEANFILES = SMDecayer__all.cc
 
 SMDecayer__all.cc : $(DIR_H_FILES) $(DIR_CC_FILES) Makefile
 	@echo "Concatenating .cc files into $@"
 	@$(top_srcdir)/cat_with_cpplines $(DIR_CC_FILES) > $@
 
 EXTRA_DIST = $(ALL_H_FILES) $(ALL_CC_FILES)
 
 DIR_H_FILES = $(addprefix $(srcdir)/,$(ALL_H_FILES))
 ALL_H_FILES = \
 EtaPiGammaGammaDecayer.h\
 EtaPiPiGammaDecayer.h \
 EtaPiPiFermionsDecayer.h \
 EtaPiPiPiDecayer.h \
 PScalar4FermionsDecayer.h\
 PScalarLeptonNeutrinoDecayer.h\
 PScalarPScalarVectorDecayer.h  \
 PScalarVectorFermionsDecayer.h\
 PScalarVectorVectorDecayer.h\
 ScalarMesonTensorScalarDecayer.h\
 ScalarScalarScalarDecayer.h  \
 SemiLeptonicScalarDecayer.h  \
 ScalarMesonFactorizedDecayer.h \
-ScalarVectorVectorDecayer.h
+ScalarVectorVectorDecayer.h \
+PseudoScalar2FermionsDecayer.h 
 
 DIR_CC_FILES = $(addprefix $(srcdir)/,$(ALL_CC_FILES))
 ALL_CC_FILES = \
 EtaPiGammaGammaDecayer.cc \
 EtaPiPiGammaDecayer.cc \
 EtaPiPiFermionsDecayer.cc \
 EtaPiPiPiDecayer.cc \
 PScalar4FermionsDecayer.cc \
 PScalarLeptonNeutrinoDecayer.cc \
 PScalarPScalarVectorDecayer.cc \
 PScalarVectorFermionsDecayer.cc  \
 PScalarVectorVectorDecayer.cc \
 ScalarMesonTensorScalarDecayer.cc \
 ScalarScalarScalarDecayer.cc \
 SemiLeptonicScalarDecayer.cc \
 ScalarMesonFactorizedDecayer.cc \
-ScalarVectorVectorDecayer.cc 
+ScalarVectorVectorDecayer.cc \
+PseudoScalar2FermionsDecayer.cc 
diff --git a/Decay/ScalarMeson/PseudoScalar2FermionsDecayer.cc b/Decay/ScalarMeson/PseudoScalar2FermionsDecayer.cc
new file mode 100644
--- /dev/null
+++ b/Decay/ScalarMeson/PseudoScalar2FermionsDecayer.cc
@@ -0,0 +1,278 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the PseudoScalar2FermionsDecayer class.
+//
+
+#include "PseudoScalar2FermionsDecayer.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Command.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+using namespace Herwig;
+using namespace ThePEG::Helicity;
+
+void PseudoScalar2FermionsDecayer::doinitrun() {
+  DecayIntegrator::doinitrun();
+  if(initialize()) {
+    for(unsigned int ix=0;ix<incoming_.size();++ix) {
+      if(mode(ix)) maxweight_[ix] = mode(ix)->maxWeight();
+    }
+  }
+}
+
+
+void PseudoScalar2FermionsDecayer::doinit() {
+  DecayIntegrator::doinit();
+  // check the parameters arew consistent
+  unsigned int isize=coupling_.size();
+  if(isize!=incoming_.size()  || isize!=outgoing_.size() ||
+     isize!=maxweight_.size())
+    throw InitException() << "Inconsistent parameters in VectorMeson2"
+			   << "FermionDecayer::doiin() " << Exception::runerror;
+  // set up the integration channels
+  PhaseSpaceModePtr mode;
+  for(unsigned int ix=0;ix<incoming_.size();++ix) {
+    tPDPtr    in  =  getParticleData(incoming_[ix]);
+    tPDVector out = {getParticleData(outgoing_[ix].first),
+		     getParticleData(outgoing_[ix].second)};
+    if(in&&out[0]&&out[1]) 
+      mode = new_ptr(PhaseSpaceMode(in,out,maxweight_[ix]));
+    else
+      mode=PhaseSpaceModePtr();
+    addMode(mode);
+  }
+}
+
+PseudoScalar2FermionsDecayer::PseudoScalar2FermionsDecayer() {
+  // don't include intermediates
+  generateIntermediates(false);
+}
+
+int PseudoScalar2FermionsDecayer::modeNumber(bool & cc,tcPDPtr parent,
+					     const tPDVector & children) const {
+  if(children.size()!=2) return -1;
+  int id(parent->id());
+  int idbar = parent->CC() ? parent->CC()->id() : id;
+  int id1(children[0]->id());
+  int id1bar = children[0]->CC() ? children[0]->CC()->id() : id1;
+  int id2(children[1]->id());
+  int id2bar = children[1]->CC() ? children[1]->CC()->id() : id2;
+  int imode(-1);
+  unsigned int ix(0);
+  cc=false;
+  do {
+    if(incoming_[ix]==id   ) {
+      if((id1   ==outgoing_[ix].first&&id2   ==outgoing_[ix].second)||
+	 (id2   ==outgoing_[ix].first&&id1   ==outgoing_[ix].second)) imode=ix;
+    }
+    if(incoming_[ix]==idbar) {
+      if((id1bar==outgoing_[ix].first&&id2bar==outgoing_[ix].second)||
+	 (id2bar==outgoing_[ix].first&&id1bar==outgoing_[ix].second)) {
+	imode=ix;
+	cc=true;
+      }
+    }
+    ++ix;
+  }
+  while(imode<0&&ix<incoming_.size());
+  return imode;
+}
+
+IBPtr PseudoScalar2FermionsDecayer::clone() const {
+  return new_ptr(*this);
+}
+
+IBPtr PseudoScalar2FermionsDecayer::fullclone() const {
+  return new_ptr(*this);
+}
+
+void PseudoScalar2FermionsDecayer::persistentOutput(PersistentOStream & os) const {
+  os << coupling_ << incoming_ << outgoing_ << maxweight_;
+}
+
+void PseudoScalar2FermionsDecayer::persistentInput(PersistentIStream & is, int) {
+  is >> coupling_ >> incoming_ >> outgoing_ >> maxweight_;
+}
+
+
+// The following static variable is needed for the type
+// description system in ThePEG.
+DescribeClass<PseudoScalar2FermionsDecayer,DecayIntegrator>
+describeHerwigPseudoScalar2FermionsDecayer("Herwig::PseudoScalar2FermionsDecayer", "HwSMDecay.so");
+
+void PseudoScalar2FermionsDecayer::Init() {
+
+  static ClassDocumentation<PseudoScalar2FermionsDecayer> documentation
+    ("The PseudoScalar2FermionsDecayer class implements the decay of a pseudoscalar meson "
+     "to a fermion and antifermion.");
+  
+  static Command<PseudoScalar2FermionsDecayer> interfaceSetUpDecayMode
+    ("SetUpDecayMode",
+     "Set up the particles (incoming, fermion, antifermion), coupling and max weight for a decay",
+     &PseudoScalar2FermionsDecayer::setUpDecayMode, false);
+
+}
+
+void PseudoScalar2FermionsDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+  unsigned int iferm(0),ianti(1);
+  // set up the spin information for the decay products
+  if(outgoing_[imode()].first!=decay[iferm]->id()) swap(iferm,ianti);
+  ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
+					incoming,true);
+  SpinorBarWaveFunction::
+    constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
+  SpinorWaveFunction::
+    constructSpinInfo(wave_   ,decay[ianti],outgoing,true);
+}
+
+
+double PseudoScalar2FermionsDecayer::me2(const int,const Particle & part,
+				       const tPDVector & outgoing,
+				       const vector<Lorentz5Momentum> & momenta,
+				       MEOption meopt) const {
+  if(!ME())
+    ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin0,PDT::Spin1Half,PDT::Spin1Half)));
+  // fermion and antifermion
+  unsigned int iferm(0),ianti(1);
+  if(outgoing_[imode()].first!=outgoing[iferm]->id()) swap(iferm,ianti);
+  // initialization
+  if(meopt==Initialize) {
+    ScalarWaveFunction::calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&part),incoming);
+  }
+  wave_.resize(2);
+  wavebar_.resize(2);
+  for(unsigned int ix=0;ix<2;++ix) {
+    wavebar_[ix] = HelicityFunctions::dimensionedSpinorBar(-momenta[iferm],ix,Helicity::outgoing);
+    wave_   [ix] = HelicityFunctions::dimensionedSpinor   (-momenta[ianti],ix,Helicity::outgoing);
+  }
+  // prefactor
+  InvEnergy pre(coupling_[imode()]/part.mass());
+  // now compute the ME
+  for(unsigned ix=0;ix<2;++ix) {
+    for(unsigned iy=0;iy<2;++iy) {
+      Complex temp = pre*wave_[ix].pseudoScalar(wavebar_[iy]);
+      if(iferm>ianti) (*ME())(0,ix,iy)=temp;
+      else            (*ME())(0,iy,ix)=temp;
+    }
+  }
+  double me = ME()->contract(rho_).real();
+  // test of the matrix element
+  // double test = 2*sqr(coupling_[imode()])*(1.-sqr((momenta[0].mass()-momenta[1].mass())/part.mass()));
+  // cerr << "testing matrix element for " << part.PDGName() << " -> " 
+  //      << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " "
+  //      << me << " " << test << " " << (me-test)/(me+test) << "\n";
+  // return the answer
+  return me;
+}
+
+bool PseudoScalar2FermionsDecayer::twoBodyMEcode(const DecayMode & dm,int & mecode,
+						double & coupling) const {
+  int imode(-1);
+  int id(dm.parent()->id()),idbar(id);
+  if(dm.parent()->CC()){idbar=dm.parent()->CC()->id();}
+  ParticleMSet::const_iterator pit(dm.products().begin());
+  int id1((**pit).id()),id1bar(id1);
+  if((**pit).CC()){id1bar=(**pit).CC()->id();}
+  ++pit;
+  int id2((**pit).id()),id2bar(id2);
+  if((**pit).CC()){id2bar=(**pit).CC()->id();}
+  unsigned int ix(0); bool order(false);
+  do {
+    if(id   ==incoming_[ix]) {
+      if(id1==outgoing_[ix].first&&id2==outgoing_[ix].second) {
+	imode=ix;
+	order=true;
+      }
+      if(id2==outgoing_[ix].first&&id1==outgoing_[ix].second) {
+	imode=ix;
+	order=false;
+      }
+    }
+    if(idbar==incoming_[ix]&&imode<0) {
+      if(id1bar==outgoing_[ix].first&&id2bar==outgoing_[ix].second) {
+	imode=ix;
+	order=true;
+      }
+      if(id2bar==outgoing_[ix].first&&id1bar==outgoing_[ix].second) {
+	imode=ix;
+	order=false;
+      }
+    }
+    ++ix;
+  }
+  while(ix<incoming_.size()&&imode<0);
+  coupling=coupling_[imode];
+  mecode=2;
+  return order;
+}
+
+// output the setup information for the particle database
+void PseudoScalar2FermionsDecayer::dataBaseOutput(ofstream & output,
+						bool header) const {
+  if(header) output << "update decayers set parameters=\"";
+  // parameters for the DecayIntegrator base class
+  DecayIntegrator::dataBaseOutput(output,false);
+  // the rest of the parameters
+  for(unsigned int ix=0;ix<incoming_.size();++ix) {
+    output << "do " << name() << ":SetUpDecayMode " << incoming_[ix] << " "
+	   << outgoing_[ix].first << " " << outgoing_[ix].second << " "
+	   << coupling_[ix] << " " << maxweight_[ix] << "\n";
+  }
+  if(header) output << "\n\" where BINARY ThePEGName=\"" 
+		    << fullName() << "\";" << endl;
+}
+
+string PseudoScalar2FermionsDecayer::setUpDecayMode(string arg) {
+  // parse first bit of the string
+  string stype = StringUtils::car(arg);
+  arg          = StringUtils::cdr(arg);
+  // extract PDG code for the incoming particle
+  long in = stoi(stype);
+  tcPDPtr pData = getParticleData(in);
+  if(!pData)
+    return "Incoming particle with id " + std::to_string(in) + "does not exist";
+  if(pData->iSpin()!=PDT::Spin0)
+    return "Incoming particle with id " + std::to_string(in) + "does not have spin 1";
+  // and outgoing particles
+  stype = StringUtils::car(arg);
+  arg   = StringUtils::cdr(arg);
+  pair<long,long> out;
+  out.first = stoi(stype);
+  pData = getParticleData(out.first);
+  if(!pData)
+    return "First outgoing particle with id " + std::to_string(out.first) + "does not exist";
+  if(pData->iSpin()!=PDT::Spin1Half)
+    return "First outgoing particle with id " + std::to_string(out.first) + "does not have spin 1/2";
+  stype = StringUtils::car(arg);
+  arg   = StringUtils::cdr(arg);
+  out.second = stoi(stype);
+  pData = getParticleData(out.second);
+  if(!pData)
+    return "Second outgoing particle with id " + std::to_string(out.second) + "does not exist";
+  if(pData->iSpin()!=PDT::Spin1Half)
+    return "Second outgoing particle with id " + std::to_string(out.second) + "does not have spin 1/2";
+  // get the coupling
+  stype = StringUtils::car(arg);
+  arg   = StringUtils::cdr(arg);
+  double g = stof(stype);
+  // and the maximum weight
+  stype = StringUtils::car(arg);
+  arg   = StringUtils::cdr(arg);
+  double wgt = stof(stype);
+  // store the information
+  incoming_.push_back(in);
+  outgoing_.push_back(out);
+  coupling_.push_back(g);
+  maxweight_.push_back(wgt);
+  // success
+  return "";
+}
diff --git a/Decay/ScalarMeson/PseudoScalar2FermionsDecayer.h b/Decay/ScalarMeson/PseudoScalar2FermionsDecayer.h
new file mode 100644
--- /dev/null
+++ b/Decay/ScalarMeson/PseudoScalar2FermionsDecayer.h
@@ -0,0 +1,195 @@
+// -*- C++ -*-
+#ifndef Herwig_PseudoScalar2FermionsDecayer_H
+#define Herwig_PseudoScalar2FermionsDecayer_H
+//
+// This is the declaration of the PseudoScalar2FermionsDecayer class.
+//
+
+#include "Herwig/Decay/DecayIntegrator.h"
+#include "Herwig/Decay/PhaseSpaceMode.h"
+#include "ThePEG/Helicity/LorentzSpinorBar.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * The PseudoScalar2FermionsDecayer class is designed for the decay of a pseudoscalar meson to
+ * two fermions, in reality this is always the decay of the \f$\eta_c\f$ to a baryon
+ * antibaryon pair.
+ *
+ * @see \ref PseudoScalar2FermionsDecayerInterfaces "The interfaces"
+ * defined for PseudoScalar2FermionsDecayer.
+ */
+class PseudoScalar2FermionsDecayer: public DecayIntegrator {
+
+public:
+
+  /**
+   * The default constructor.
+   */
+  PseudoScalar2FermionsDecayer();
+
+  /**
+   * Which of the possible decays is required
+   * @param cc Is this mode the charge conjugate
+   * @param parent The decaying particle
+   * @param children The decay products
+   */
+  virtual int modeNumber(bool & cc, tcPDPtr parent, 
+			 const tPDVector & children) const;
+
+  /**
+   * Return the matrix element squared for a given mode and phase-space channel.
+   * @param ichan The channel we are calculating the matrix element for. 
+   * @param part The decaying Particle.
+   * @param outgoing The particles produced in the decay
+   * @param momenta  The momenta of the particles produced in the decay
+   * @param meopt Option for the calculation of the matrix element
+   * @return The matrix element squared for the phase-space configuration.
+   */
+  double me2(const int ichan,const Particle & part,
+	     const tPDVector & outgoing,
+	     const vector<Lorentz5Momentum> & momenta,
+	     MEOption meopt) const;
+
+  /**
+   *   Construct the SpinInfos for the particles produced in the decay
+   */
+  virtual void constructSpinInfo(const Particle & part,
+				 ParticleVector outgoing) const;
+
+  /**
+   * Specify the \f$1\to2\f$ matrix element to be used in the running width calculation.
+   * @param dm The DecayMode
+   * @param mecode The code for the matrix element as described
+   *               in the GenericWidthGenerator class, in this case 2.
+   * @param coupling The coupling for the matrix element.
+   * @return True or False if this mode can be handled.
+   */
+  bool twoBodyMEcode(const DecayMode & dm, int & mecode, double & coupling) const;
+
+  /**
+   * Output the setup information for the particle database
+   * @param os The stream to output the information to
+   * @param header Whether or not to output the information for MySQL
+   */
+  virtual void dataBaseOutput(ofstream & os,bool header) const;
+
+public:
+
+  /** @name Functions used by the persistent I/O system. */
+  //@{
+  /**
+   * Function used to write out object persistently.
+   * @param os the persistent output stream written to.
+   */
+  void persistentOutput(PersistentOStream & os) const;
+
+  /**
+   * Function used to read in object persistently.
+   * @param is the persistent input stream read from.
+   * @param version the version number of the object when written.
+   */
+  void persistentInput(PersistentIStream & is, int version);
+  //@}
+
+  /**
+   * The standard Init function used to initialize the interfaces.
+   * Called exactly once for each class by the class description system
+   * before the main function starts or
+   * when this class is dynamically loaded.
+   */
+  static void Init();
+
+protected:
+
+  /** @name Clone Methods. */
+  //@{
+  /**
+   * Make a simple clone of this object.
+   * @return a pointer to the new object.
+   */
+  virtual IBPtr clone() const;
+
+  /** Make a clone of this object, possibly modifying the cloned object
+   * to make it sane.
+   * @return a pointer to the new object.
+   */
+  virtual IBPtr fullclone() const;
+  //@}
+
+protected:
+
+  /** @name Standard Interfaced functions. */
+  //@{
+  /**
+   * Initialize this object after the setup phase before saving and
+   * EventGenerator to disk.
+   * @throws InitException if object could not be initialized properly.
+   */
+  virtual void doinit();
+
+  /**
+   * Initialize this object to the begining of the run phase.
+   */
+  virtual void doinitrun();
+  //@}
+
+private:
+
+  /**
+   * The assignment operator is private and must never be called.
+   * In fact, it should not even be implemented.
+   */
+  PseudoScalar2FermionsDecayer & operator=(const PseudoScalar2FermionsDecayer &) = delete;
+
+public:
+
+  /**
+   *   Set the parameters for a decay mode
+   */
+  string setUpDecayMode(string arg);
+
+private:
+
+  /**
+   * coupling for a decay
+   */
+  vector<double> coupling_;
+
+  /**
+   * the PDG codes for the incoming particles
+   */
+  vector<long> incoming_;
+
+  /**
+   * the PDG codes for the outgoing fermion-antifermion
+   */
+  vector<pair<long,long>> outgoing_;
+
+  /**
+   * maximum weight for a decay
+   */
+  vector<double> maxweight_;
+
+  /**
+   *  Spin density matrix
+   */
+  mutable RhoDMatrix rho_;
+
+  /**
+   *  Spinors for the decay products
+   */
+  mutable vector<Helicity::LorentzSpinor   <SqrtEnergy> > wave_;
+
+  /**
+   *  barred spinors for the decay products
+   */
+  mutable vector<Helicity::LorentzSpinorBar<SqrtEnergy> > wavebar_;
+
+};
+
+}
+
+#endif /* Herwig_PseudoScalar2FermionsDecayer_H */