diff --git a/EventRecord/SpinInfo.h b/EventRecord/SpinInfo.h
--- a/EventRecord/SpinInfo.h
+++ b/EventRecord/SpinInfo.h
@@ -1,473 +1,473 @@
 // -*- 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. */
     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),
       _oldDeveloped(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), _oldDeveloped(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 ;
 
   /**
    * Calculate the rho matrix for the decay if not already done.
    */
-  void undecay() const ;
+  virtual void undecay() const ;
 
   /**
    * Set the developed flag and calculate the D matrix for the decay.
    */
   void develop() const ;
 
   /**
    *  Needs update
    */
   void needsUpdate() const {
     if(_developed!=NeedsUpdate) _oldDeveloped = _developed;
     _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;
 
   /**
    * Has the particle been developed?  (I.e. has the D matrix encoding
    * the info about the decay been calculated)
    */
   mutable DevelopedStatus _oldDeveloped;
 
   /**
    * 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/Helicity/FermionSpinInfo.h b/Helicity/FermionSpinInfo.h
--- a/Helicity/FermionSpinInfo.h
+++ b/Helicity/FermionSpinInfo.h
@@ -1,175 +1,183 @@
 // -*- C++ -*-
 //
 // FermionSpinInfo.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_FermionSpinInfo_H
 #define ThePEG_FermionSpinInfo_H
 // This is the declaration of the FermionSpinInfo class.
 
 #include "ThePEG/EventRecord/SpinInfo.h"
 #include "ThePEG/Helicity/LorentzSpinor.h"
 #include "FermionSpinInfo.fh"
 #include <array>
 
 namespace ThePEG {
 namespace Helicity {
 
 /**
  *  The FermionSpinInfo class inherits from the SpinInfo class and
  *  implements the storage of the basis vectors for a spin-1/2
  *  particle.  The basis states are the u-spinors for a particle and
  *  the v-spinors for an antiparticle. The barred spinors can be
  *  obtained from these.
  *
  *  These basis states should be set by either matrixelements or
  *  decayers which are capable of generating spin correlation
  *  information.
  *
  *  The basis states in the rest frame of the particles can then be
  *  accessed by decayers to produce the correct correlations.
  *
  *  N.B. in our convention 0 is the \f$-\frac12\f$ helicity state and
  *  1 is the \f$+\frac12\f$ helicity state.
  *
  * @author Peter Richardson
  */
 class FermionSpinInfo: public SpinInfo {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * Default constructor.
    */
   FermionSpinInfo()
     : SpinInfo(PDT::Spin1Half), _decaycalc(false) {}
 
   /**
    * Standard Constructor.
    * @param p the production momentum.
    * @param time true if the particle is time-like.
    */
   FermionSpinInfo(const Lorentz5Momentum & p, bool time)
     : SpinInfo(PDT::Spin1Half, p, time), _decaycalc(false) {}
   //@}
 
 public:
 
   /** @name Set and get methods for the basis state. */
   //@{
   /**
    * Set the basis state, this is production state.
    * @param hel the helicity (0 or 1 as described above.)
    * @param in the LorentzSpinor for the given helicity.
    */
   void setBasisState(unsigned int hel,
 		     const LorentzSpinor<SqrtEnergy> & in) const {
     assert(hel<2);
     _productionstates[hel] = in;
     _currentstates   [hel] = in;
   }
 
   /**
    * Set the basis state for the decay.
    * @param hel the helicity (0 or 1 as described above.)
    * @param in the LorentzSpinor for the given helicity.
    */
   void setDecayState(unsigned int hel,
 		     const LorentzSpinor<SqrtEnergy> & in) const {
     assert(hel<2);
     _decaycalc = true;
     _decaystates[hel] = in;
   }
 
   /**
    * Get the basis state for the production for the given helicity, \a
    * hel (which is 0 or 1 as described above.)
    */
   const LorentzSpinor<SqrtEnergy> & getProductionBasisState(unsigned int hel) const {
     assert(hel<2);
     return _productionstates[hel];
   }
 
   /**
    * Get the current basis state for the given helicity, \a
    * hel (which is 0 or 1 as described above.)
    */
   const LorentzSpinor<SqrtEnergy> & getCurrentBasisState(unsigned int hel) const {
     assert(hel<2);
     return _currentstates[hel];
   }
 
   /**
    * Get the basis state for the decay for the given helicity, \a hel
    * (which is 0 or 1 as described above.)
    */
   const LorentzSpinor<SqrtEnergy> & getDecayBasisState(unsigned int hel) const {
     assert(hel<2);
     if(!_decaycalc) {
       for(unsigned int ix=0;ix<2;++ix) _decaystates[ix]=_currentstates[ix];
       _decaycalc=true;
     }
     return _decaystates[hel];
   }
   //@}
 
   /**
    * Perform a lorentz rotation of the spin information
    */
   virtual void transform(const LorentzMomentum &,const LorentzRotation &);
 
+  /**
+   *  Undecay
+   */
+  virtual void undecay() const {
+    _decaycalc=false;
+    SpinInfo::undecay();
+  }
+  
 public:
 
   /**
    * Standard Init function.
    */
   static void Init();
 
   /**
    * Standard clone method.
    */
   virtual EIPtr clone() const;
 
 private:
 
   /**
    * Private and non-existent assignment operator.
    */
   FermionSpinInfo & operator=(const FermionSpinInfo &);
 
 private:
 
   /**
    * basis states in the frame in which the particle was produced
    */
   mutable std::array<LorentzSpinor<SqrtEnergy>,2> _productionstates;
 
   /**
    * basis states in the current frame of the particle
    */
   mutable std::array<LorentzSpinor<SqrtEnergy>,2> _currentstates;
 
   /**
    * basis states in the frame in which the particle decays
    */
   mutable std::array<LorentzSpinor<SqrtEnergy>,2> _decaystates;
 
   /**
    * True if the decay state has been set.
    */
   mutable bool _decaycalc;
 
 };
 
 }
 }
 
 namespace ThePEG {
 
 }
 #endif /* ThePEG_FermionSpinInfo_H */
diff --git a/Helicity/LorentzSpinorBar.h b/Helicity/LorentzSpinorBar.h
--- a/Helicity/LorentzSpinorBar.h
+++ b/Helicity/LorentzSpinorBar.h
@@ -1,355 +1,355 @@
 // -*- C++ -*-
 //
 // LorentzSpinorBar.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_LorentzSpinorBar_H
 #define ThePEG_LorentzSpinorBar_H
 // This is the declaration of the LorentzSpinorBar class.
 
 #include "ThePEG/Config/ThePEG.h"
 #include "ThePEG/Vectors/LorentzRotation.h"
 #include "ThePEG/Vectors/ThreeVector.h"
 #include "HelicityDefinitions.h"
 #include "LorentzSpinor.fh"
 #include "LorentzSpinorBar.fh"
 
 namespace ThePEG {
 namespace Helicity {
 
 /**
  *  The LorentzSpinorBar class implements the storage of a barred
  *  LorentzSpinor. The design is based on that of the LorentzSpinor
  *  class where the details of the implemented are discussed in more
  *  detail.
  *
  * @see LorentzSpinor
  *
  * @author Peter Richardson
  */
 
 template<typename Value>
 class LorentzSpinorBar {
 public:
 
   /** @name Standard constructors. */
   //@{
   /**
    * Default zero constructor, optionally specifying \a t, the type
    */
   LorentzSpinorBar(SpinorType t = SpinorType::unknown) : _type(t), _spin() {}
 
   /**
    * Constructor with complex numbers specifying the components,
    * optionally specifying \a t, the type
    */
   LorentzSpinorBar(complex<Value> a, complex<Value> b,
 		   complex<Value> c, complex<Value> d,
 		   SpinorType t = SpinorType::unknown)
     : _type(t), _spin{{a,b,c,d}} {}
 
   template <typename U>
   LorentzSpinorBar(const LorentzSpinorBar<U> & other)
     : _type(other._type), _spin(other._spin) {}
 
   //@}
 
   /** @name Access the components. */
   //@{
   /**
    * Subscript operator to return spinor components
    */
   complex<Value> operator[](int i) const {
     assert( i>= 0 && i <= 3 );
     return _spin[i];
   }
 
   /**
    * Subscript operator to return spinor components
    */
   complex<Value> operator()(int i) const {
     assert( i>= 0 && i <= 3 );
     return _spin[i];
   }
 
   /**
    * Set components by index.
    */
   complex<Value> & operator()(int i) {
     assert( i>= 0 && i <= 3 );
     return _spin[i];
   }
 
   /**
    * Set components by index.
    */
   complex<Value> & operator[](int i) {
     assert( i>= 0 && i <= 3 );
     return _spin[i];
   }
 
   /**
    * Get first component.
    */
   complex<Value> s1() const {return _spin[0];}
 
   /**
    * Get second component.
    */
   complex<Value> s2() const {return _spin[1];}
 
   /**
    * Get third component.
    */
   complex<Value> s3() const {return _spin[2];}
 
   /**
    * Get fourth component.
    */
   complex<Value> s4() const {return _spin[3];}
 
   /**
    * Set first component.
    */
   void setS1(complex<Value> in) {_spin[0]=in;}
 
   /**
    * Set second component.
    */
   void setS2(complex<Value> in) {_spin[1]=in;}
 
   /**
    * Set third component.
    */
   void setS3(complex<Value> in) {_spin[2]=in;}
 
   /**
    * Set fourth component.
    */
   void setS4(complex<Value> in) {_spin[3]=in;}
   //@}
 
   /// @name Mathematical assignment operators.
   //@{
   template <typename ValueB>
   LorentzSpinorBar<Value> & operator+=(const LorentzSpinorBar<ValueB> & a) {
     for(unsigned int ix=0;ix<4;++ix) _spin[ix] += a._spin[ix];
     return *this;
   }
   
   template <typename ValueB>
   LorentzSpinorBar<Value> & operator-=(const LorentzSpinorBar<ValueB> & a) {
     for(unsigned int ix=0;ix<4;++ix) _spin[ix] -= a._spin[ix];
     return *this;
   }
   
   LorentzSpinorBar<Value> & operator*=(double a) {
     for(unsigned int ix=0;ix<4;++ix) _spin[ix] *=a;
     return *this;
   }
 
   LorentzSpinorBar<Value> & operator/=(double a) {
     for(unsigned int ix=0;ix<4;++ix) _spin[ix] /=a;
     return *this;
   }
   //@}
 
   /** @name Transformations. */
   //@{
   /**
    * Return the barred spinor
    */
   LorentzSpinor<Value> bar() const;
 
   /**
    * Return the conjugated spinor \f$u_c=C\bar{u}^T\f$. This operation
    * transforms u-spinors to v-spinors and vice-versa and is required when
    * dealing with majorana particles.
    */
   LorentzSpinorBar conjugate() const;
 
   /**
    * Standard Lorentz boost specifying the components of the beta vector.
    */
   LorentzSpinorBar & boost(double,double,double);
 
   /**
    * Standard Lorentz boost specifying the beta vector.
    */
   LorentzSpinorBar & boost(const Boost &);
 
   /**
    * General Lorentz transformation
    */
   LorentzSpinorBar & transform(const SpinHalfLorentzRotation &) ;
   
   /**
    * General Lorentz transformation
    */
   LorentzSpinorBar & transform(const LorentzRotation & r) {
     transform(r.half());
     return *this;
   }
   //@}
 
   /** @name Functions related to type. */
   //@{
   /**
    * Return the type of the spinor.
    */
   SpinorType Type() const {return _type;}
   //@}
 
   /**
    *  @name Functions to apply the projection operator
    */
   //@{
   /**
    *   Apply \f$p\!\!\!\!\!\not\,\,\,+m\f$
    */
   template<typename ValueB> 
   auto projectionOperator(const LorentzVector<ValueB> & p, 
                           const ValueB & m) const 
-  -> LorentzSpinorBar<decltype(m*s1())>
+  -> LorentzSpinorBar<decltype(m*Value())>
   {
-    typedef decltype(m*s1()) ResultT;
+    typedef decltype(m*Value()) ResultT;
     LorentzSpinorBar<ResultT> spin;
     static const Complex ii(0.,1.);
     complex<ValueB> p0pp3=p.t()+p.z();
     complex<ValueB> p0mp3=p.t()-p.z();
     complex<ValueB> p1pp2=p.x()+ii*p.y();
     complex<ValueB> p1mp2=p.x()-ii*p.y();
     spin.setS1(m*s1()+p0pp3*s3()+p1pp2*s4());
     spin.setS2(m*s2()+p0mp3*s4()+p1mp2*s3());
     spin.setS3(m*s3()+p0mp3*s1()-p1pp2*s2());
     spin.setS4(m*s4()+p0pp3*s2()-p1mp2*s1());
     return spin;
   }
 
   /**
    *  Apply \f$g^LP_L+g^RP_R\f$
    */
   LorentzSpinorBar
   helicityProjectionOperator(const Complex & gL, const Complex & gR) const {
     LorentzSpinorBar spin;
     spin.setS1(gL*s1());
     spin.setS2(gL*s2());
     spin.setS3(gR*s3());
     spin.setS4(gR*s4());
     return spin;
   }
   //@}
 
 private:
   /**
    * Type of spinor
    */
   SpinorType _type;
 
   /**
    * Storage of the components.
    */
   std::array<complex<Value>,4> _spin;
 };
 
 /// @name Basic mathematical operations
 //@{
 template <typename Value>
 inline LorentzSpinorBar<double>
 operator/(const LorentzSpinorBar<Value> & v, Value a) {
   return LorentzSpinorBar<double>(v.s1()/a, v.s2()/a, v.s3()/a, v.s4()/a,v.Type());
 }
 
 inline LorentzSpinorBar<double>
 operator/(const LorentzSpinorBar<double> & v, Complex a) {
   return LorentzSpinorBar<double>(v.s1()/a, v.s2()/a, v.s3()/a, v.s4()/a,v.Type());
 }
 
 template <typename Value>
 inline LorentzSpinorBar<Value> operator-(const LorentzSpinorBar<Value> & v) {
   return LorentzSpinorBar<Value>(-v.s1(),-v.s2(),-v.s3(),-v.s4(),v.Type());
 }
 
 template <typename ValueA, typename ValueB>
 inline LorentzSpinorBar<ValueA>
 operator+(LorentzSpinorBar<ValueA> a, const LorentzSpinorBar<ValueB> & b) {
   return a += b;
 }
 
 template <typename ValueA, typename ValueB>
 inline LorentzSpinorBar<ValueA>
 operator-(LorentzSpinorBar<ValueA> a, const LorentzSpinorBar<ValueB> & b) {
   return a -= b;
 }
 
 template <typename Value>
 inline LorentzSpinorBar<Value>
 operator*(const LorentzSpinorBar<Value> & a, double b) {
   return LorentzSpinorBar<Value>(a.s1()*b, a.s2()*b, a.s3()*b, a.s4()*b,a.Type());
 }
 
 template <typename Value>
 inline LorentzSpinorBar<Value>
 operator*(double b, LorentzSpinorBar<Value> a) {
   return a *= b;
 }
   
 template <typename Value>
 inline LorentzSpinorBar<Value>
 operator*(const LorentzSpinorBar<Value> & a, Complex b) {
   return LorentzSpinorBar<Value>(a.s1()*b, a.s2()*b, a.s3()*b, a.s4()*b,a.Type());
 }
 
 template <typename ValueA, typename ValueB>
 inline auto operator*(complex<ValueB> a, const LorentzSpinorBar<ValueA> & v) 
   -> LorentzSpinorBar<decltype(a.real()*v.s1().real())>
 {
   return {a*v.s1(), a*v.s2(), a*v.s3(), a*v.s4(),v.Type()};
 }
 
 template <typename ValueA, typename ValueB>
 inline auto operator*(const LorentzSpinorBar<ValueA> & v, complex<ValueB> b) 
   -> LorentzSpinorBar<decltype(b.real()*v.s1().real())>
 {
   return b*v;
 }
 
 template <typename ValueA, typename ValueB>
 inline auto operator/(const LorentzSpinorBar<ValueA> & v, complex<ValueB> b) 
   -> LorentzSpinorBar<decltype(v.s1().real()/b.real())>
 {
   return {v.s1()/b, v.s2()/b, v.s3()/b, v.s4()/b,v.Type()};
 }
   
 template <typename ValueA, typename ValueB>
 inline auto operator*(ValueB a, const LorentzSpinorBar<ValueA> & v) 
   -> LorentzSpinorBar<decltype(a*v.s1().real())>
 {
   return {a*v.s1(), a*v.s2(), a*v.s3(), a*v.s4(),v.Type()};
 }
 
 template <typename ValueA, typename ValueB>
 inline auto operator*(const LorentzSpinorBar<ValueA> & v, ValueB b) 
   -> LorentzSpinorBar<decltype(b*v.s1().real())>
 {
   return b*v;
 }
 
 template <typename ValueA, typename ValueB>
 inline auto operator/(const LorentzSpinorBar<ValueA> & v, ValueB b) 
   -> LorentzSpinorBar<decltype(v.s1().real()/b)>
 {
   return {v.s1()/b, v.s2()/b, v.s3()/b, v.s4()/b,v.Type()};
 }
 
 }
 }
 
 #ifndef ThePEG_TEMPLATES_IN_CC_FILE
 #include "LorentzSpinorBar.tcc"
 #endif 
 
 #endif
diff --git a/Helicity/RSFermionSpinInfo.h b/Helicity/RSFermionSpinInfo.h
--- a/Helicity/RSFermionSpinInfo.h
+++ b/Helicity/RSFermionSpinInfo.h
@@ -1,170 +1,178 @@
 // -*- C++ -*-
 //
 // RSFermionSpinInfo.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_RSFermionSpinInfo_H
 #define THEPEG_RSFermionSpinInfo_H
 // This is the declaration of the RSFermionSpinInfo class.
 
 #include "ThePEG/EventRecord/SpinInfo.h"
 #include "ThePEG/Helicity/LorentzRSSpinor.h"
 #include "RSFermionSpinInfo.fh"
 #include <array>
 
 namespace ThePEG {
 namespace Helicity {
 
 /**
  *  The RSFermionSpinInfo class inherits from the SpinInfo class and
  *  implements the storage of the basis vector for a spin-3/2 particle.
  *  The basis states are the vector u spinors for a particle and the vector
  *  v-spinors for an antiparticle. The barred spinors can be obtained from these.
  *
  *  These basis states should be set by either the matrixelements or decayers
  *  which are capable of generating spin correlation information.
  *
  *  The basis states in the rest frame of the particles can then be accessed by
  *  the decayers to produce the correct correlations.
  *
  *  N.B. in our convention 0 is the \f$-\frac32\f$ helicity state,
  *  1 is the \f$-\frac12\f$ helicity state,
  *  2 is the \f$+\frac12\f$ helicity state,
  *  3 is the \f$+\frac32\f$ helicity state.
  *
  * @see SpinInfo
  *
  * \author Peter Richardson
  *
  */
 class RSFermionSpinInfo: public SpinInfo {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * Default constructor.
    */
   RSFermionSpinInfo()  : SpinInfo(PDT::Spin3Half), _decaycalc(false) {}
 
   /**
    * Standard Constructor.
    * @param p the production momentum.
    * @param time true if the particle is time-like.
    */
   RSFermionSpinInfo(const Lorentz5Momentum & p,bool time)
     : SpinInfo(PDT::Spin3Half, p, time), _decaycalc(false) {}
   //@}
 
 public:
 
   /** @name Set and get methods for the basis state. */
   //@{
   /**
    * Set the basis state, this is production state.
    * @param hel the helicity (0,1,2,3 as described above.)
    * @param in the LorentzRSSpinor for the given helicity.
    */
   void setBasisState(unsigned int hel,
 		     const LorentzRSSpinor<SqrtEnergy> & in) const {
     assert(hel<4);
     _productionstates[hel] = in;
     _currentstates   [hel] = in;
   }
 
   /**
    * Set the basis state for the decay.
    * @param hel the helicity (0,1,2,3 as described above.)
    * @param in the LorentzRSSpinor for the given helicity.
    */
   void setDecayState(unsigned int hel,
 		     const LorentzRSSpinor<SqrtEnergy> & in) const {
     assert(hel<4);
     _decaycalc = true;
     _decaystates[hel] = in;
   }
 
   /**
    * Get the basis state for the production for the given helicity, \a
    * hel (0,1,2,3 as described above.)
    */
   const LorentzRSSpinor<SqrtEnergy> & getProductionBasisState(unsigned int hel) const {
     assert(hel<4);
     return _productionstates[hel];
   }
 
   /**
    * Get the basis state for the decay for the given helicity, \a hel
    * (0,1,2,3 as described above.)
    */
   const LorentzRSSpinor<SqrtEnergy> & getDecayBasisState(unsigned int hel) const {
     assert(hel<4);
     if(!_decaycalc) {
       for(unsigned int ix=0;ix<4;++ix) _decaystates[ix]=_currentstates[ix];
       _decaycalc=true;
     }
     return _decaystates[hel];
   }
 
   /**
    * Perform a lorentz rotation of the spin information
    */
   virtual void transform(const LorentzMomentum &,const LorentzRotation &);
   //@}
 
+  /**
+   *  Undecay
+   */
+  virtual void undecay() const {
+    _decaycalc=false;
+    SpinInfo::undecay();
+  }
+
 public:
 
   /**
    * Standard Init function used to initialize the interfaces.
    */
   static void Init();
 
   /**
    * Standard clone method.
    */
   virtual EIPtr clone() const;
 
 private:
 
   /**
    * Private and non-existent assignment operator.
    */
   RSFermionSpinInfo & operator=(const RSFermionSpinInfo &);
 
 private:
 
   /**
    * Basis states in the frame in which the particle was produced.
    */
   mutable std::array<LorentzRSSpinor<SqrtEnergy>,4> _productionstates;
 
   /**
    * Basis states in the frame in which the particle decays.
    */
   mutable std::array<LorentzRSSpinor<SqrtEnergy>,4> _decaystates;
 
   /**
    * Basis states in the current frame of the particle
    */
   mutable std::array<LorentzRSSpinor<SqrtEnergy>,4> _currentstates;
 
   /**
    * True if the decay state has been set.
    */
   mutable bool _decaycalc;
 
 };
 
 }
 }
 
 
 
 namespace ThePEG {
 
 }
 #endif /* THEPEG_RSFermionSpinInfo_H */
diff --git a/Helicity/TensorSpinInfo.h b/Helicity/TensorSpinInfo.h
--- a/Helicity/TensorSpinInfo.h
+++ b/Helicity/TensorSpinInfo.h
@@ -1,168 +1,176 @@
 // -*- C++ -*-
 //
 // TensorSpinInfo.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_TensorSpinInfo_H
 #define THEPEG_TensorSpinInfo_H
 // This is the declaration of the TensorSpinInfo class.
 
 #include "ThePEG/EventRecord/SpinInfo.h"
 #include "ThePEG/Helicity/LorentzTensor.h"
 #include "TensorSpinInfo.fh"
 // #include "TensorSpinInfo.xh"
 #include <array>
 
 namespace ThePEG {
 namespace Helicity {
 
 /**
  *  The TensorSpinInfo class is the implementation of the spin
  *  information for tensor particles.  It inherits from the SpinInfo
  *  class and implements the storage of the basis tensors.
  *
  *  These basis states should be set by either matrix elements or
  *  decayers which are capable of generating spin correlation
  *  information.
  *
  *  The basis states in the rest frame of the particles can then be
  *  accessed by decayers to produce the correct correlation.
  *
  *  N.B. in our convention 0 is the \f$-2\f$ helicity state,
  *  1 is the \f$-1\f$ helicity state,
  *  2 is the \f$0\f$ helicity state,
  *  3 is the \f$+1\f$ helicity state and
  *  4 is the \f$+2\f$ helicity state.
  *
  * @author Peter Richardson
  *
  */
 class TensorSpinInfo: public SpinInfo {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * Default constructor.
    */
   TensorSpinInfo() : SpinInfo(PDT::Spin2), _decaycalc(false) {}
   
   /**
    * Standard Constructor.
    * @param p the production momentum.
    * @param time true if the particle is time-like.
    */
   TensorSpinInfo(const Lorentz5Momentum & p, bool time)
     : SpinInfo(PDT::Spin2, p, time), _decaycalc(false) {}
   //@}
 
 public:
 
   /** @name Access the basis states. */
   //@{
   /**
    * Set the basis state, this is production state.
    * @param hel the helicity (0,1,2,3,4 as described above.)
    * @param in the LorentzTensor for the given helicity.
    */
   void setBasisState(unsigned int hel, LorentzTensor<double> in) const {
     assert(hel<5);
     _productionstates[hel]=in;
     _currentstates   [hel]=in;
   }
 
   /**
    * Set the basis state for the decay.
    * @param hel the helicity (0,1,2,3,4 as described above.)
    * @param in the LorentzTensor for the given helicity.
    */
   void setDecayState(unsigned int hel, LorentzTensor<double> in) const {
     assert(hel<5);
     _decaycalc = true;
     _decaystates[hel] = in;
   }
 
   /**
    * Get the basis state for the production for the given helicity, \a
    * hel  (0,1,2,3,4 as described above.)
    */
   const LorentzTensor<double> & getProductionBasisState(unsigned int hel) const {
     assert(hel<5);
     return _productionstates[hel];
   }
 
   /**
    * Get the basis state for the decay for the given helicity, \a hel
    * (0,1,2,3,4 as described above.)
    */
   const LorentzTensor<double> & getDecayBasisState(unsigned int hel) const {
     assert(hel<5);
     if(!_decaycalc) {
       for(unsigned int ix=0;ix<5;++ix)
 	_decaystates[ix]=_currentstates[ix].conjugate();
       _decaycalc=true;
     }
     return _decaystates[hel];
   }
   //@}
 
   /**
    * Perform a lorentz rotation of the spin information
    */
   virtual void transform(const LorentzMomentum &,const LorentzRotation &);
 
+  /**
+   *  Undecay
+   */
+  virtual void undecay() const {
+    _decaycalc=false;
+    SpinInfo::undecay();
+  }
+
 public:
 
   /**
    * Standard Init function.
    */
   static void Init();
 
   /**
    * Standard clone method.
    */
   virtual EIPtr clone() const;
 
 private:
 
   /**
    * Private and non-existent assignment operator.
    */
   TensorSpinInfo & operator=(const TensorSpinInfo &);
 
 private:
 
   /**
    * Basis states in the frame in which the particle was produced.
    */
   mutable std::array<LorentzTensor<double>,5> _productionstates;
 
   /**
    * Basis states in the frame in which the particle decays.
    */
   mutable std::array<LorentzTensor<double>,5> _decaystates;
 
   /**
    * Basis states in the current frame of the particle
    */
   mutable std::array<LorentzTensor<double>,5> _currentstates;
 
   /**
    * True if the decay state has been set.
    */
   mutable bool _decaycalc;
 
 };
 
 }
 }
 
 
 namespace ThePEG {
 
 }
 #endif /* THEPEG_TensorSpinInfo_H */
diff --git a/Helicity/VectorSpinInfo.h b/Helicity/VectorSpinInfo.h
--- a/Helicity/VectorSpinInfo.h
+++ b/Helicity/VectorSpinInfo.h
@@ -1,167 +1,175 @@
 // -*- C++ -*-
 //
 // VectorSpinInfo.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_VectorSpinInfo_H
 #define THEPEG_VectorSpinInfo_H
 // This is the declaration of the VectorSpinInfo class.
 
 #include "ThePEG/EventRecord/SpinInfo.h"
 #include "ThePEG/Helicity/LorentzPolarizationVector.h"
 #include "VectorSpinInfo.fh"
 
 namespace ThePEG {
 namespace Helicity {
 
 /**
  *  The VectorSpinInfo class is the implementation of the spin
  *  information for vector particles.  It inherits from the SpinInfo
  *  class and implements the storage of the basis vectors.
  *
  *  These basis states should be set by either matrixelements or
  *  decayers which are capable of generating spin correlation
  *  information.
  *
  *  The basis states in the rest frame of the particles can then be
  *  accessed by decayers to produce the correct correlation.
  *
  *  N.B. in our convention 0 is the \f$-1\f$ helicity state,
  *  1 is the \f$0\f$ helicity state and
  *  2 is the \f$+1\f$ helicity state.
  *
  * @author Peter Richardson
  *
  */
 class VectorSpinInfo: public SpinInfo {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * Default constructor.
    */
   VectorSpinInfo() : SpinInfo(PDT::Spin1), _decaycalc(false) {}
   
   /**
    * Standard Constructor.
    * @param p the production momentum.
    * @param time true if the particle is time-like.
    */
   VectorSpinInfo(const Lorentz5Momentum & p, bool time)
     : SpinInfo(PDT::Spin1, p, time), _decaycalc(false) {}
   //@}
   
 public:
 
   /** @name Set and get methods for the basis state. */
   //@{
   /**
    * Set the basis state, this is production state.
    * @param hel the helicity (0,1,2 as described above.)
    * @param in the LorentzPolarizationVector for the given helicity.
    */
   void setBasisState(unsigned int hel, 
 		     const LorentzPolarizationVector & in) const {
     assert(hel<3);
     _productionstates[hel] = in;
     _currentstates   [hel] = in;
   }
 
   /**
    * Set the basis state for the decay.
    * @param hel the helicity (0,1,2 as described above.)
    * @param in the LorentzPolarizationVector for the given helicity.
    */
   void setDecayState(unsigned int hel, 
 		     const LorentzPolarizationVector & in) const {
     assert(hel<3);
     _decaycalc = true;
     _decaystates[hel] = in;;
   }
 
   /**
    * Get the basis state for the production for the given helicity, \a
    * hel (0,1,2 as described above.)
    */
   const LorentzPolarizationVector & getProductionBasisState(unsigned int hel) const {
     assert(hel<3);
     return _productionstates[hel];
   }
 
   /**
    * Get the basis state for the decay for the given helicity, \a hel 
    * (0,1,2 as described above.)
    */
   const LorentzPolarizationVector & getDecayBasisState(unsigned int hel) const {
     assert(hel<3);
     if(!_decaycalc) {
       for(unsigned int ix=0;ix<3;++ix)
 	_decaystates[ix]=_currentstates[ix].conjugate();
       _decaycalc=true;
     }
     // return the basis function
     return _decaystates[hel];
   }
   //@}
 
   /**
    * Perform a Lorentz rotation of the spin information
    */
   virtual void transform(const LorentzMomentum &,const LorentzRotation & );
 
+  /**
+   *  Undecay
+   */
+  virtual void undecay() const {
+    _decaycalc=false;
+    SpinInfo::undecay();
+  }
+
 public:
 
   /**
    * Standard Init function used to initialize the interfaces.
    */
   static void Init();
 
   /**
    * Standard clone method.
    */
   virtual EIPtr clone() const;
 
 private:
 
   /**
    * Private and non-existent assignment operator.
    */
   VectorSpinInfo & operator=(const VectorSpinInfo &);
 
 private:
 
   /**
    * Basis states in the frame in which the particle was produced.
    */
   mutable std::array<LorentzPolarizationVector,3> _productionstates;
 
   /**
    * Basis states in the frame in which the particle decays.
    */
   mutable std::array<LorentzPolarizationVector,3> _decaystates;
 
   /**
    * Basis states in the current frame of the particle
    */
   mutable std::array<LorentzPolarizationVector,3> _currentstates;
 
   /**
    * True if the decay state has been set.
    */
   mutable bool _decaycalc;
 
 };
 
 }
 }
 
 
 namespace ThePEG {
 
 }
 #endif /* THEPEG_VectorSpinInfo_H */