diff --git a/EventRecord/SpinInfo.cc b/EventRecord/SpinInfo.cc
--- a/EventRecord/SpinInfo.cc
+++ b/EventRecord/SpinInfo.cc
@@ -1,182 +1,187 @@
 // -*- C++ -*-
 //
 // SpinInfo.cc is a part of ThePEG - Toolkit for HEP Event Generation
 // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
 //
 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the SpinInfo class.
 //
 // Author: Peter Richardson
 //
 #include "SpinInfo.h"
 #include "ThePEG/Repository/CurrentGenerator.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "Event.h"
 
 using namespace ThePEG;
 
 const double SpinInfo::_eps=1.0e-8;
 
 SpinInfo::SpinInfo(const SpinInfo & x)
   : EventInfoBase(x), _production(x._production), _decay(x._decay),
     _timelike(x._timelike),
     _prodloc(x._prodloc), _decayloc(x._decayloc),
     _decayed(x._decayed), _developed(x._developed),_rhomatrix(x._rhomatrix),
     _Dmatrix(x._Dmatrix),_spin(x._spin),
     _productionmomentum(x._productionmomentum),
     _decaymomentum(x._decaymomentum),
     _currentmomentum(x._currentmomentum) 
 {
   x._production=VertexPtr();
   x._decay=VertexPtr();
   // set the vertex so it now points to the copy
   if(_production) {
     // for timelike
     if(_timelike) _production->resetOutgoing(this,_prodloc); 
     // for spacelike
     else          _production->resetIncoming(this,_prodloc);
   }
 }
 
 EIPtr SpinInfo::clone() const {
   tcSpinPtr temp=this;
   return const_ptr_cast<SpinPtr>(temp);
 }
 
 void SpinInfo::rebind(const EventTranslationMap & trans) {
   if(_production) _production = trans.translate(_production);
   if(_decay)      _decay      = trans.translate(_decay);
   EventInfoBase::rebind(trans);
 }
 
 
 NoPIOClassDescription<SpinInfo> SpinInfo::initSpinInfo;
 // Definition of the static class description member.
 
 void SpinInfo::Init() {}
 
 void SpinInfo::update() const {
   // number of instances fo this object
   int nref=referenceCount();
   if(nref<2||nref>3) return;
   // work out the number of references there should be
   int nmin=0;
   // check the production pointers
   if(_production) {
     if(_timelike) { if(_production->outgoing()[_prodloc]==this) ++nmin;}
     else          { if(_production->incoming()[_prodloc]==this) ++nmin;}
   }
   // check the decay pointers
   if(_decay) {
     if(_decay->incoming()[_decayloc]==this) ++nmin;
   }
   // delete the pointers
   SpinPtr temp;
   if(nmin+1==nref) {
     // delete the production pointers
     if(_production) {
       if(_timelike) {
 	if(_production->outgoing()[_prodloc]==this)
 	  _production->resetOutgoing(SpinPtr(),_prodloc);
       }
       else {
 	if(_production->incoming()[_prodloc]==this)
 	  _production->resetIncoming(SpinPtr(),_prodloc);
       }
     }
     // delete the decay pointers
     if(_decay) {
       if(_decay->incoming()[_decayloc]==this)
 	_decay->resetIncoming(SpinPtr(),_decayloc);
     }
   }
 }
 
 void SpinInfo::decay(bool recursive) const {
   // if the particle has already been decayed do nothing
   if(_decayed) return;
   // otherwise we need to obtain the correct rho (timelike) or D (spacelike) matrix
   assert(_developed!=NeedsUpdate);
   if(_timelike) {
     if(_developed==Developed&&iSpin()!=PDT::Spin0) _developed=NeedsUpdate;
     if(productionVertex()) {
       if(recursive) redecay();
       else _rhomatrix = productionVertex()->getRhoMatrix(_prodloc,true);
     }
   }
   else {
     if(_developed==Developed&&iSpin()!=PDT::Spin0) _developed=NeedsUpdate;
     if(_production) _Dmatrix = _production->getDMatrix(_prodloc);
   }
   _decaymomentum = _currentmomentum;
   _decayed=true;
 }
 
 void SpinInfo::redevelop() const {
   assert(developed()==NeedsUpdate);
   // calculate rho/D matrix
   if(_timelike) {
     _Dmatrix   = decayVertex() ? 
       decayVertex()->getDMatrix(decayLocation()) : RhoDMatrix(iSpin());
   }
   else {
     _rhomatrix = decayVertex() ? 
       decayVertex()->getRhoMatrix(decayLocation(),false) :  RhoDMatrix(iSpin());
   }
   // update the D matrix of this spininfo
   _developed = Developed;
   // update the parent if needed
   if(productionVertex() &&
      productionVertex()->incoming().size()==1) {
     tcSpinPtr parent = _timelike ? 
       productionVertex()->incoming()[0] : productionVertex()->outgoing()[0];
-    parent->needsUpdate();
-    parent->redevelop();
+    if ( parent->developed() != StopUpdate ) {
+      parent->needsUpdate();
+      parent->redevelop();
+    }
   }
 }
 
 void SpinInfo::develop() const {
   // if the particle has already been developed do nothing
   switch(_developed) {
   case Developed:
     return;
+  case StopUpdate:
+    return;
   case NeedsUpdate:
     redevelop();
     return;
   case Undeveloped:
     if(_timelike) {
       if(_decay) _Dmatrix = _decay->getDMatrix(_decayloc);
       else       _Dmatrix = RhoDMatrix(iSpin());
     }
     else {
       if(_decay) _rhomatrix = _decay->getRhoMatrix(_decayloc,false);
       else       _rhomatrix = RhoDMatrix(iSpin());
     }
     _developed=Developed;
     return;
   }
 }
 
 void SpinInfo::redecay() const {
   if(!productionVertex()) return;
   if(productionVertex()->incoming().size()==1) {
     tcSpinPtr parent;
     if(productionVertex()->incoming()[0]->timelike())
       parent = productionVertex()->incoming()[0];
     else {
       if(productionVertex()->outgoing()[0]!=this)
 	parent = productionVertex()->outgoing()[0];
       else
 	parent = productionVertex()->outgoing()[1];
     }
-    parent->redecay();
+    if ( parent->developed() != StopUpdate )
+      parent->redecay();
   }
   if(timelike())
     _rhomatrix = productionVertex()->getRhoMatrix(_prodloc,true);
   else
     _Dmatrix   = productionVertex()->getDMatrix(_prodloc);
 }
diff --git a/EventRecord/SpinInfo.h b/EventRecord/SpinInfo.h
--- a/EventRecord/SpinInfo.h
+++ b/EventRecord/SpinInfo.h
@@ -1,450 +1,457 @@
 // -*- C++ -*-
 //
 // SpinInfo.h is a part of ThePEG - Toolkit for HEP Event Generation
 // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
 //
 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
 //
 #ifndef ThePEG_SpinInfo_H
 #define ThePEG_SpinInfo_H
 // This is the declaration of the SpinInfo class.
 
 #include "ThePEG/EventRecord/EventInfoBase.h"
 #include "ThePEG/PDT/PDT.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "HelicityVertex.h"
 
 namespace ThePEG {
 
 /**
  *   The SpinInfo is the base class for the spin information for the
  *   spin correlation algorithm. The implementations for different spin
  *   states inherit from this.
  *
  *   The class contains pointers to the vertex where the particle is
  *   produced and where it decays, together with methods to set/get
  *   these.
  *
  *   There are two flags decayed which store information on the state
  *   of the particle.
  *
  *   The decayed() members provides access to the _decay data member
  *   which is true if the spin density matrix required to perform the
  *   decay of a timelike particle has been calculated (this would be a
  *   decay matrix for a spacelike particle.) This is set by the
  *   decay() method which calls a method from the production vertex to
  *   calculate this matrix. The decay() method should be called by a
  *   decayer which uses spin correlation method before it uses the
  *   spin density matrix to calculate the matrix element for the
  *   decay.
  *
  *   The developed() member provides access to the _developed data
  *   member which is true if the decay matrix required to perform the
  *   decays of the siblings of a particle has been calculated (this
  *   would a spin density matrix for a spacelike particle.) This is
  *   set by the developed() method which calls a method from the decay
  *   vertex to calculate the matrix. The developed() method is called
  *   by a DecayHandler which is capable of performing spin
  *   correlations after all the unstable particles produced by a
  *   decaying particle are decayed.
  *
  *   Methods are also provided to access the spin density and decay
  *   matrices for a particle.
  *
  * @author Peter Richardson
  *
  */
 class SpinInfo: public EventInfoBase {
 
 public:
 
   /**
    *  Status for the implementation of spin correlations
    */
   enum DevelopedStatus {
     Undeveloped=0, /**< Not developed. */
     Developed=1,   /**< Developed. */
-    NeedsUpdate=2  /**< Developed but needs recalculating due to some change. */
+    NeedsUpdate=2,  /**< Developed but needs recalculating due to some change. */
+    StopUpdate=3     /**< Stop recalculating at this spin info. */
   };
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * Default constructor.
    */
   SpinInfo() 
     : _timelike(false), _prodloc(-1), _decayloc(-1), 
       _decayed(false), _developed(Undeveloped) {}
 
   /**
    * Standard Constructor.
    * @param s the spin.
    * @param p the production momentum.
    * @param time true if the particle is time-like.
    */
   SpinInfo(PDT::Spin s, 
 	   const Lorentz5Momentum & p = Lorentz5Momentum(), 
 	   bool time = false)
     : _timelike(time), _prodloc(-1), _decayloc(-1),
       _decayed(false), _developed(Undeveloped),
       _rhomatrix(s), _Dmatrix(s), _spin(s),
       _productionmomentum(p), _currentmomentum(p) {}
 
   /**
    * Copy-constructor.
    */
   SpinInfo(const SpinInfo &);
   //@}
 
 public:
 
   /**
    * Returns true if the polarization() has been implemented in a
    * subclass. This default version returns false.
    */
   virtual bool hasPolarization() const { return false; }
 
   /**
    * Return the angles of the polarization vector as a pair of
    * doubles. first is the polar angle and second is the azimuth
    * wrt. the particles direction. This default version of the
    * function returns 0,0, and if a subclass implements a proper
    * function it should also implement 'hasPolarization()' to return
    * true.
    */
   virtual DPair polarization() const { return DPair(); }
 
 public:
 
   /**
    * Standard Init function.
    */
   static void Init();
 
   /**
    * Rebind to cloned objects. If a FermionSpinInfo is cloned together
    * with a whole Event and this has pointers to other event record
    * objects, these should be rebound to their clones in this
    * function.
    */
   virtual void rebind(const EventTranslationMap & trans);
 
   /**
    * Standard clone method.
    */
   virtual EIPtr clone() const;
 
   /**
    * Method to handle the delelation
    */
   void update() const;
 
   /**
    * Perform a lorentz rotation of the spin information
    */
   virtual void transform(const LorentzMomentum & m, const LorentzRotation & r) {
     _currentmomentum = m;
     _currentmomentum.transform(r);
   }
 
 public:
 
 
   /** @name Access the vertices. */
   //@{
   /**
    * Set the vertex at which the particle was produced.
    */
   void productionVertex(VertexPtr in) const {
     _production=in;
     // add to the list of outgoing if timelike
     int temp(-1);
     if(_timelike) in->addOutgoing(this,temp); 
     // or incoming if spacelike
     else          in->addIncoming(this,temp);
     _prodloc=temp;
   }
 
   /**
    * Get the vertex at which the particle was produced.
    */
   tcVertexPtr productionVertex() const { return _production; }
 
   /**
    * Set the vertex at which the particle decayed or branched.
    */
   void decayVertex(VertexPtr in) const {
     if(in) {
       _decay=in;
       if(_timelike) {
 	int temp(-1);
 	in->addIncoming(this,temp);
 	_decayloc=temp;
 	assert(temp==0);
       }
       else {
 	int temp(-1);
 	in->addOutgoing(this,temp);
 	_decayloc=temp;
       }
     }
     else {
       _decay=VertexPtr();
       _decayloc=-1;
     }
   }
 
   /**
    * Get the vertex at which the particle decayed or branched.
    */
   tcVertexPtr decayVertex() const { return _decay; }
   //@}
 
   /** @name Access information about the associated particle. */
   //@{
   /**
    * Has the particle decayed?
    */
   bool decayed() const { return _decayed; }
 
   /**
    * Set if the particle has decayed.
    */
   void decayed(bool b) const { _decayed = b; }
 
   /**
    * Return true if the decay matrix required to perform the decays of
    * the siblings of a particle has been calculated.
    */
   DevelopedStatus developed() const { return _developed; }
 
   /**
    * Calculate the rho matrix for the decay if not already done.
    */
   void decay(bool recursive=false) const ;
 
   /**
    * Set the developed flag and calculate the D matrix for the decay.
    */
   void develop() const ;
 
   /**
    *  Needs update
    */
   void needsUpdate() const {_developed=NeedsUpdate;}
 
   /**
+   *  Used for an unstable particle to *temporarily* stop
+   *  redevelop and redecay at that particle
+   */
+  void stopUpdate() const {_developed=StopUpdate;}
+
+  /**
    * Return 2s+1 for the particle
    */
   PDT::Spin iSpin() const { return _spin; }
 
   /**
    * Return the momentum of the particle when it was produced.
    */
   const Lorentz5Momentum & productionMomentum() const {
     return _productionmomentum;
   }
 
   /**
    *  The current momentum of the particle
    */
   const Lorentz5Momentum & currentMomentum() const {
     return _currentmomentum;
   }
 
   /**
    * Return true if particle is timelike (rather than spacelike).
    */
   bool timelike() const { return _timelike; }
   //@}
 
   /**
    *  Access to the locations
    */
   //@{
   /**
    *  Production Location
    */
   int productionLocation() const {return _prodloc;}
 
   /**
    *  Decay Location
    */
   int decayLocation() const {return _decayloc;}
   //@}
 
 public:
 
   /** @name Access the rho and D matrices. */
   //@{
   /**
    * Access the rho matrix.
    */
   RhoDMatrix rhoMatrix() const { return _rhomatrix; }
 
   /**
    * Access the rho matrix.
    */
   RhoDMatrix & rhoMatrix() { return _rhomatrix; }
 
   /**
    * Access the D matrix.
    */
   RhoDMatrix DMatrix() const { return _Dmatrix; }
 
   /**
    * Access the D matrix.
    */
   RhoDMatrix & DMatrix() { return _Dmatrix; }
   //@}
 
 protected:
 
   /**
    *  Check if momentum is near to the current momentum
    */
   bool isNear(const Lorentz5Momentum & p) {
     return currentMomentum().isNear(p,_eps);
   }
 
 private:
 
   /**
    * Describe a concrete class without persistent data.
    */
   static NoPIOClassDescription<SpinInfo> initSpinInfo;
 
   /**
    * Private and non-existent assignment operator.
    */
   SpinInfo & operator=(const SpinInfo &);
 
 private:
 
   /**
    * Set the developed flag and calculate the D matrix for the decay,
    * and all decays further up the chain.
    */
   void redevelop() const ;
 
   /**
    * Recursively recalulate all the rho matrices from the top of the chain
    */
   void redecay() const ;
 
 private:
 
   /**
    * Pointer to the production vertex for the particle
    */
   mutable VertexPtr _production;
 
   /**
    * Pointers to the decay vertex for the particle
    */
   mutable VertexPtr _decay;
 
   /**
    * Is this is timelike (true) or spacelike (false ) particle?  This
    * is used to decide if the particle is incoming or outgoing at the
    * production vertex
    */
   bool _timelike;
 
   /**
    * Location in the hard vertex array at production.
    */
   mutable int _prodloc;
 
   /**
    * Location in the hard vertex array at decay.
    */
   mutable int _decayloc;
 
   /**
    * Has the particle been decayed?  (I.e. has the rho matrix for the
    * decay been calculated.)
    */
   mutable bool _decayed;
 
   /**
    * Has the particle been developed?  (I.e. has the D matrix encoding
    * the info about the decay been calculated)
    */
   mutable DevelopedStatus _developed;
 
   /**
    * Storage of the rho matrix.
    */
   mutable RhoDMatrix _rhomatrix;
 
   /**
    * Storage of the decay matrix
    */
   mutable RhoDMatrix _Dmatrix;
 
   /**
    * The spin of the particle
    */
   PDT::Spin _spin;
 
   /**
    * Momentum of the particle when it was produced
    */
   Lorentz5Momentum _productionmomentum;
 
   /**
    * Momentum of the particle when it decayed
    */
   mutable Lorentz5Momentum _decaymomentum;
   
   /**
    * Current momentum of the particle
    */
   Lorentz5Momentum _currentmomentum;
 
   /**
    *  A small energy for comparing momenta to check if Lorentz Transformations
    *  should be performed
    */
   static const double _eps;
 };
 
 }
 
 
 namespace ThePEG {
 
 /** @cond TRAITSPECIALIZATIONS */
 
 /**
  * This template specialization informs ThePEG about the base class of
  * SpinInfo.
  */
 template <>
 struct BaseClassTrait<ThePEG::SpinInfo,1>: public ClassTraitsType {
   /** Typedef of the base class of SpinInfo. */
   typedef EventInfoBase NthBase;
 };
 
 /**
  * This template specialization informs ThePEG about the name of the
  * SpinInfo class and the shared object where it is defined.
  */
 template <>
 struct ClassTraits<ThePEG::SpinInfo>
   : public ClassTraitsBase<ThePEG::SpinInfo> {
   /**
    * Return the class name.
    */
   static string className() { return "ThePEG::SpinInfo"; }
 };
 
 /** @endcond */
 
 }
 
 #endif /* ThePEG_SpinInfo_H */
diff --git a/Utilities/SimplePhaseSpace.h b/Utilities/SimplePhaseSpace.h
--- a/Utilities/SimplePhaseSpace.h
+++ b/Utilities/SimplePhaseSpace.h
@@ -1,232 +1,232 @@
 // -*- C++ -*-
 //
 // SimplePhaseSpace.h is a part of ThePEG - Toolkit for HEP Event Generation
 // Copyright (C) 1999-2017 Leif Lonnblad
 //
 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
 //
 #ifndef ThePEG_SimplePhaseSpace_H
 #define ThePEG_SimplePhaseSpace_H
 
 #include "ThePEG/Config/ThePEG.h"
 
 #include "ThePEG/Vectors/LorentzRotation.h"
 #include "ThePEG/Vectors/LorentzRotation.h"
 #include "ThePEG/EventRecord/Particle.h"
 #include "ThePEG/EventRecord/ParticleTraits.h"
 #include "ThePEG/Repository/UseRandom.h"
 #include "SimplePhaseSpace.xh"
 #include <numeric>
 
 namespace ThePEG {
 
 /**
  * SimplePhaseSpace defines a set of static functions to be used for
  * distributing momenta evenly in phase space. In most cases pointers
  * and references to both particle and momentum objects can be used as
  * arguments as long as the ParticleTraits class is specialized
  * properly. When needed, random numbers are generated with the
  * generator given by the static UseRandom class.
  */
 namespace SimplePhaseSpace {
 
   /**
    * Set two momenta in their center of mass system. Their total
    * invariant mass squared is given by s, and their direction is
    * distributed isotropically.
    * @param s the total invariant mass squared.
    * @param p1 pointer or reference to the first momentum. Its
    * invariant mass will be preserved.
    * @param p2 pointer or reference to the second momentum. Its
    * invariant mass will be preserved.
    * @throw ImpossibleKinematics if the sum of the invariant masses was
    * larger than the given invariant mass (\f$\sqrt{s}\f$).
    */
   template <typename PType>
   void CMS(Energy2 s, PType & p1, PType & p2);
 
   /**
    * Set two momenta in their center of mass system. Their total
    * invariant mass squared is given by s, and their direction is
    * given in terms of the polar and azimuth angle of the first
    * momenta.
    * @param s the total invariant mass squared.
    * @param p1 pointer or reference to the first momentum. Its
    * invariant mass will be preserved.
    * @param p2 pointer or reference to the second momentum. Its
    * invariant mass will be preserved.
    * @param cosTheta cosine of the azimuth angle of the first momentum.
    * @param phi azimuth angle of the first momentum.
    * @throw ImpossibleKinematics if the sum of the invariant masses was
    * larger than the given invariant mass (\f$\sqrt{s}\f$).
    */
   template <typename PType>
   void CMS(PType & p1, PType & p2, Energy2 s,
 		  double cosTheta, double phi);
 
   /**
    * Set two momenta in their center of mass system. Their total
    * invariant mass squared is given by s. The helper momentum p0 is
    * used so that afterwards \f$t=(p0-p1)^2\f$ and p1 has the azimuth
    * angle phi around p0.
    * @param p1 pointer or reference to the first momentum. Its
    * invariant mass will be preserved.
    * @param p2 pointer or reference to the second momentum. Its
    * invariant mass will be preserved.
    * @param s the total invariant mass squared.
    * @param t \f$=(p0-p1)^2\f$.
    * @param phi azimuth angle of the first momentum around p0.
    * @param p0 pointer or reference to an auxiliary momentum.
    * @throw ImpossibleKinematics if the sum of the invariant masses was
    * larger than the given invariant mass (\f$\sqrt{s}\f$).
    */
   template <typename PType>
   void CMS(PType & p1, PType & p2, Energy2 s, Energy2 t, double phi,
 		  const PType & p0);
 
   /**
    * Set two momenta in their center of mass system. Their total
    * invariant mass squared is given by s. p1 will be along the z-axis.
    * @param p1 pointer or reference to the first momentum. Its
    * invariant mass will be preserved.
    * @param p2 pointer or reference to the second momentum. Its
    * invariant mass will be preserved.
    * @param s the total invariant mass squared.
    * @throw ImpossibleKinematics if the sum of the invariant masses was
    * larger than the given invariant mass (\f$\sqrt{s}\f$).
    */
   template <typename PType>
   void CMS(PType & p1, PType & p2, Energy2 s);
 
   /**
    * Set two momenta in their center of mass system. Their total
    * invariant mass squared is given by s. The first will be along the
    * z-axis.
    * @param p a pair of pointers or references to the two momenta. Their
    * invariant masses will be preserved.
    * @param s the total invariant mass squared.
    * @throw ImpossibleKinematics if the sum of the invariant masses was
    * larger than the given invariant mass (\f$\sqrt{s}\f$).
    */
   template <typename PPairType>
   void CMS(const PPairType & p, Energy2 s)
   {
     CMS(*p.first, *p.second, s);
   }
 
   /**
    * Set three momenta in their center of mass system. Their total
    * invariant mass squared is given by s. The energy fraction of
    * particle p1(3) is x1(3) of the total energy and the angles of the
    * system is distributed isotropically.
    * @param p1 pointer or reference to the first momentum. Its
    * invariant mass will be preserved.
    * @param p2 pointer or reference to the second momentum. Its
    * invariant mass will be preserved.
    * @param p3 pointer or reference to the second momentum. Its
    * invariant mass will be preserved.
    * @param s the total invariant mass squared.
    * @param x1 the energy fraction \f$2e_1/\sqrt{s}\f$.
    * @param x3 the energy fraction \f$2e_3/\sqrt{s}\f$.
    * @throw ImpossibleKinematics if the sum of the invariant masses was
    * larger than the given invariant mass (\f$\sqrt{s}\f$).
    */
   template <typename PType>
   void CMS(PType & p1, PType & p2, PType & p3, Energy2 s,
 		  double x1, double x3);
 
   /**
    * Set three momenta in their center of mass system. Their total
    * invariant mass squared is given by s. The energy fraction of
    * particle p1(3) is x1(3) of the total energy. Particle p1 is
    * initially placed along the z-axis and particle p2 is given
    * azimuth angle phii. Then the system is then rotated with
    * theta and phi respectively.
    * @param p1 pointer or reference to the first momentum. Its
    * invariant mass will be preserved.
    * @param p2 pointer or reference to the second momentum. Its
    * invariant mass will be preserved.
    * @param p3 pointer or reference to the second momentum. Its
    * invariant mass will be preserved.
    * @param s the total invariant mass squared.
    * @param x1 the energy fraction \f$2e_1/\sqrt{s}\f$.
    * @param x3 the energy fraction \f$2e_3/\sqrt{s}\f$.
    * @param phii the azimuth angle of p2 around p1.
    * @param theta the polar angle of p1.
    * @param phi the azimuth angle of p1.
    * @throw ImpossibleKinematics if the sum of the invariant masses was
    * larger than the given invariant mass (\f$\sqrt{s}\f$).
    */
   template <typename PType>
   void CMS(PType & p1, PType & p2, PType & p3, Energy2 s,
 		  double x1, double x3, double phii = 0.0,
 		  double theta = 0.0, double phi = 0.0);
 
   /**
    * Calculate the absolute magnitude of the momenta of two particles
    * with masses m1 and m2 when put in their CMS of total invariant
    * mass squared s.
    * @param s the total invariant mass squared.
    * @param m1 the mass of particle 1.
    * @param m2 the mass of particle 2.
    * @throw ImpossibleKinematics if the sum of the masses was
    * larger than the given invariant mass (\f$\sqrt{s}\f$).
    */
   Energy getMagnitude(Energy2 s, Energy m1, Energy m2);
 
   /**
    * Return a three-vector given the absolute momentum, cos(theta) and
    * phi.
    * @param p the magnitude of the momentum.
    * @param costheta the cosine of the polar angle.
    * @param phi the azimuth angle.
    */
   inline Momentum3 polar3Vector(Energy p, double costheta, double phi)
   {
     return Momentum3(p*sqrt(1.0 - sqr(costheta))*sin(phi),
 		     p*sqrt(1.0 - sqr(costheta))*cos(phi),
 		     p*costheta);
   }
 
   /**
    * Get a number of randomly distributed momenta.
    * Given a number specified invariant masses and a
    * total invariant mass m0, return corresponding four-momenta
    * randomly distributed according to phase space. 
    * @param m0 the
    * total invariant mass of the resulting momenta.
    * @param m a vector
    * of invariant masses of the resulting momenta.
    * @return a vector
    * of momenta with the given masses randomly distributed.
    * @throw ImpossibleKinematics if the sum of the masses was
    * larger than the given invariant mass (\f$\sqrt{s}\f$).
    */
   vector<LorentzMomentum>
   CMSn(Energy m0, const vector<Energy> & m);
 
   /**
    * Set the momentum of a number of particles. Given a number of
    * particles and a total invariant mass m0, distribute their
    * four-momenta randomly according to phase space.
    * @param particles a container of particles or pointers to
    * particles. The invariant mass of these particles will not be
    * chaned.
    * @param m0 the
    * total invariant mass of the resulting momenta.
    * @throw ImpossibleKinematics if the sum of the masses was
    * larger than the given invariant mass (\f$\sqrt{s}\f$).
    */
   template <typename Container>
   void CMSn(Container & particles, Energy m0);
 
-};
+}
 
 }
 
 #ifndef ThePEG_TEMPLATES_IN_CC_FILE
 #include "SimplePhaseSpace.tcc"
 #endif
 
 #endif /* ThePEG_SimplePhaseSpace_H */