diff --git a/Config/std.h b/Config/std.h
--- a/Config/std.h
+++ b/Config/std.h
@@ -1,185 +1,187 @@
 // -*- C++ -*-
 //
 // std.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_std_H
 #define ThePEG_std_H
 
 /** \file
 
  * This file introduces a number of <code>std::</code> classes into
  * the ThePEG namespace. Also introduces some useful functions for
  * standard library classes.
  *
  * Do not make changes in this file. If you want to use alternatives
  * to the <code>std::</code> classes in ThePEG, edit a copy of this
  * file and include it in an alternative config file which can be
  * included in the main ThePEG.h config file using the macro
  * <code>ThePEG_ALTERNATE_CONFIG</code>.
  */
 
 #include <map>
 #include <set>
 #include <string>
+#include <array>
 #include <vector>
 #include <list>
 #include <iostream>
 #include <fstream>
 #include <sstream>
 #include <iomanip>
 #include <stack>
 #include <utility>
 #include <typeinfo>
 #include <stdexcept>
 #include <cmath>
 
 namespace std {
 
 /** @cond TRAITSPECIALIZATIONS */
 /**
  * This specialization of the std::less class is needed in order to be
  * able use put pointers to type_info objects as keys in maps and
  * sets.
  */
 template <>
 struct less<const type_info *> :
     public binary_function<const type_info *, const type_info *, bool> 
 {
   /**
    * This is the function called when comparing two pointers to
    * type_info.
    */
   bool operator()(const type_info * x, const type_info * y) const {
     return x->before(*y); }
 };
 /** @endcond */
 
 }
 
 namespace ThePEG {
 
+using std::array;
 using std::deque;
 using std::stack;
 using std::vector;
 using std::multiset;
 using std::set;
 using std::map;
 using std::list;
 using std::multimap;
 using std::pair;
 using std::make_pair;
 using std::less;
 using std::string;
 using std::type_info;
 using std::exception;
 using std::range_error;
 using std::ios;
 using std::ostream;
 using std::istream;
 using std::ofstream;
 using std::ifstream;
 using std::ostringstream;
 using std::istringstream;
 using std::cin;
 using std::cout;
 using std::cerr;
 using std::endl;
 using std::flush;
 using std::setprecision;
 using std::setw;
 using std::swap;
 using std::min;
 using std::max;
 using std::mem_fun;
 using std::sqrt;
 //using std::pow;
 using std::abs;
 using std::atan2;
 using std::isfinite;
 
 /** Powers - standard or non-standard */
 template <class ExponentT>
 inline constexpr double pow(double x, ExponentT p) {
   return std::pow(x,double(p));
 }
 
 /** Square root of an integer. */
 inline double sqrt(int x) {
   return std::sqrt(double(x));
 }
 
 /** factorial */
 inline constexpr long double factorial(unsigned int n) {
   return (n < 2) ? 1.0 : n * factorial(n - 1);
 }
 
 /** Check if a given object is a part of a container. */
 template <typename Container, typename Key>
 inline bool member(const Container & c, const Key & k) {
   return c.find(k) != c.end();
 }
 
 /** Check if a given object is a part of a vector. */
 template <typename T, typename Key>
 inline bool member(const vector<T> & v, const Key & k) {
   for ( typename vector<T>::const_iterator i = v.begin(); i != v.end(); ++i )
     if ( *i == k ) return true;
   return false;
   // return find(v.begin(), v.end(), k) != v.end();
 }
 
 /** Return an insert iterator for a given container. */
 template <typename Cont>
 inline std::insert_iterator<Cont> inserter(Cont & c) {
   return std::insert_iterator<Cont>(c, c.end());
 }
 
 
 /** Return an insert iterator for a given vector. Overrides the
  *  general version. */
 template <typename T, typename A>
 inline std::back_insert_iterator< vector<T,A> > inserter(vector<T,A> & v) {
   return back_inserter(v);
 }
 
 /** Return an insert iterator for a given vector. Overrides the
  *  general version. */
 template <typename T, typename A>
 inline std::back_insert_iterator< deque<T,A> > inserter(deque<T,A> & v) {
   return back_inserter(v);
 }
 
 /** Stream manipulator setting an ostream to left-adjust its ouput. */
 inline ostream& left(ostream& os) {
   os.setf(ios::left, ios::adjustfield);
   return os;
 }
 
 /** Stream manipulator setting an ostream to right-adjust its ouput. */
 inline ostream& right(ostream& os) {
   os.setf(ios::right, ios::adjustfield);
   return os;
 }
 
 }
 
 /** Macro for declaring a set. */
 #define ThePEG_DECLARE_SET(VALTYPE,NAME)                               \
   /** A set of VALTYPE. */                                             \
   typedef set<VALTYPE, less<VALTYPE> > NAME
 
 /** Macro for declaring a multiset. */
 #define ThePEG_DECLARE_MULTISET(VALTYPE,NAME)                          \
   /** A multiset of VALTYPE. */                                        \
   typedef multiset<VALTYPE, less<VALTYPE> > NAME
 
 /** Macro for declaring a map. */
 #define ThePEG_DECLARE_MAP(KEYTYPE,VALTYPE,NAME)                       \
   /** A map of VALTYPE indexed by KEYTYPE. */                          \
   typedef map<KEYTYPE, VALTYPE, less<KEYTYPE> > NAME
 
 #endif /* ThePEG_std_H */
diff --git a/Helicity/Vertex/Vector/VVVVVertex.h b/Helicity/Vertex/Vector/VVVVVertex.h
--- a/Helicity/Vertex/Vector/VVVVVertex.h
+++ b/Helicity/Vertex/Vector/VVVVVertex.h
@@ -1,169 +1,169 @@
 // -*- C++ -*-
 //
 // VVVVVertex.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_VVVVVertex_H
 #define ThePEG_VVVVVertex_H
 //
 // This is the declaration of the VVVVVertex class.
 
 #include "ThePEG/Helicity/Vertex/AbstractVVVVVertex.h"
 #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
 #include "VVVVVertex.fh"
 
 namespace ThePEG {
 namespace Helicity{
 
 /** \ingroup Helicity
  *
  * This is the implementation of the four vector vertex. 
  * It is based on the AbstractVVVVVertex class for the storage of particles 
  * which are allowed to interact at the vertex.
  * Classes implementation a specific vertex should inherit from this 
  * one and implement the virtual setCoupling member.
  *
  * The form of the vertex is
  * \f[ic^2\left[
  *   2\epsilon_1\cdot\epsilon_2\epsilon_3\cdot\epsilon_4-
  *   \epsilon_1\cdot\epsilon_3\epsilon_2\cdot\epsilon_4-
  *   \epsilon_1\cdot\epsilon_4\epsilon_2\cdot\epsilon_3
  * \right]\f]
  *  optional the additional diagrams from the three point vertices can be included.
  *
  * @see AbstractVVVVVertex
  */
 class VVVVVertex: public AbstractVVVVVertex {
   
 public:
 
   /**
    * Standard Init function used to initialize the interfaces.
    */
   static void Init();
   
 public:
   
   /**
    * Members to calculate the helicity amplitude expressions for vertices
    * and off-shell particles.
    */
   //@{
   /**
    * Evaluate the vertex.
    * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
    * @param iopt Evaluation option, 0 just evaluate the four point vertex, 1
    * include all the three point diagrams as well.
    * @param vec1 The wavefunction for the first  vector.
    * @param vec2 The wavefunction for the second vector.
    * @param vec3 The wavefunction for the third  vector.
    * @param vec4 The wavefunction for the fourth vector.
    */
   Complex evaluate(Energy2 q2, int iopt,
 		   const VectorWaveFunction & vec1, const VectorWaveFunction & vec2,
 		   const VectorWaveFunction & vec3, const VectorWaveFunction & vec4);
   //@}
 
   /**
    *   Set coupling methods
    */
   //@{
   /**
    * Dummy for a three point interaction.
    */
   virtual void setCoupling(Energy2,tcPDPtr,tcPDPtr,tcPDPtr) {
     assert(false);
   }
 
   /**
    * Calculate the couplings for a four point interaction.
    * 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 part4 The ParticleData pointer for the fourth particle.
    */
   virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,tcPDPtr part3,
 			   tcPDPtr part4)=0;
   //@}
   
 protected:
   
   /**
    * Set the order of the particles.
    * @param id1 The PDG code of the first  particle.
    * @param id2 The PDG code of the second particle.
    * @param id3 The PDG code of the third  particle.
    * @param id4 The PDG code of the fourth particle.
    */
   void setOrder(int id1,int id2,int id3,int id4) {
     _iorder[0]=id1;
     _iorder[1]=id2;
     _iorder[2]=id3;
     _iorder[3]=id4;
   }
 
   /**
    * Set the type of the vertex.
    * @param itype The type of vertex (QCD=1 or electroweak=2).
    */
   void setType(int itype) {
     _itype=itype;
   }
 
   /**
    * Set the intermediate particles if including s/u/t channel terms.
    * @param part1 The ParticleData pointer for the first  particle.
    * @param part2 The ParticleData pointer for the second particle.
    * @param c1 The coupling for the first  particle.
    * @param c2 The coupling for the second particle.
    */
   void setIntermediate(tcPDPtr part1,tcPDPtr part2,Complex c1,Complex c2) {
     _inter[0]=part1;
     _inter[1]=part2;
     _coup[0]=c1;
     _coup[1]=c2;
   }
   
 private:
   
   /**
    * Private and non-existent assignment operator.
    */
   VVVVVertex & operator=(const VVVVVertex &);
   
 private:
 
   /**
    * Type of vertex 1=QCD 2=EW.
    */
   int _itype;
 
   /**  
    * Order of the particles.
    */
-  int _iorder[4];
+  array<int,4> _iorder;
 
   /**
    *  Intermediate particles
    */
-  tcPDPtr _inter[2];
+  array<tcPDPtr,2> _inter;
 
   /**
    * Couplings of the intermediate particles.
    */
-  Complex _coup[2];
+  array<Complex,2> _coup;
 
 };
 }
 }
 
 namespace ThePEG {
 
 }
 #endif /* ThePEG_VVVVVertex_H */
diff --git a/LesHouches/LesHouches.h b/LesHouches/LesHouches.h
--- a/LesHouches/LesHouches.h
+++ b/LesHouches/LesHouches.h
@@ -1,253 +1,253 @@
 // -*- C++ -*-
 //
 // LesHouches.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_LesHouches_H
 #define THEPEG_LesHouches_H
 //
 // This is the declaration of the LesHouches class.
 //
 
 #include "ThePEG/Config/ThePEG.h"
 
 namespace ThePEG {
 
 /**
  * The HEPRUP class is a simple container corresponding to the Les
  * Houches accord (hep-ph/0109068) common block with the same
  * name. The members are named in the same way as in the common
  * block. However, fortran arrays are represented by vectors, except
  * for the arrays of length two which are represented by pair objects.
  */
 class HEPRUP {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * Default constructor.
    */
   HEPRUP() : IDWTUP(0), NPRUP(0) {}
   //@}
 
 public:
 
   /**
    * Set the NPRUP variable, corresponding to the number of
    * sub-processes, to \a nrup, and resize all relevant vectors
    * accordingly.
    */
   void resize(int nrup) {
     NPRUP = nrup;
     resize();
   }
 
   /**
    * Assuming the NPRUP variable, corresponding to the number of
    * sub-processes, is correctly set, resize the relevant vectors
    * accordingly.
    */
   void resize() {
     XSECUP.resize(NPRUP);
     XERRUP.resize(NPRUP);
     XMAXUP.resize(NPRUP);
     LPRUP.resize(NPRUP);
   }
 
   /**
    * PDG id's of beam particles. (first/second is in +/-z direction).
    */
   pair<long,long> IDBMUP;
 
   /**
    * Energy of beam particles given in GeV.
    */
   pair<double,double> EBMUP;
 
   /**
    * The author group for the PDF used for the beams according to the
    * PDFLib specification.
    */
   pair<int,int> PDFGUP;
 
   /**
    * The id number the PDF used for the beams according to the
    * PDFLib specification.
    */
   pair<int,int> PDFSUP;
 
   /**
    * Master switch indicating how the ME generator envisages the
    * events weights should be interpreted according to the Les Houches
    * accord.
    */
   int IDWTUP;
 
   /**
    * The number of different subprocesses in this file (should
    * typically be just one)
    */
   int NPRUP;
 
   /**
    * The cross sections for the different subprocesses in pb.
    */
   vector<double> XSECUP;
 
   /**
    * The statistical error in the cross sections for the different
    * subprocesses in pb.
    */
   vector<double> XERRUP;
 
   /**
    * The maximum event weights (in XWGTUP) for different subprocesses.
    */
   vector<double> XMAXUP;
 
   /**
    * The subprocess code for the different subprocesses.
    */
   vector<int> LPRUP;
 
 };
 
 
 /**
  * The HEPEUP class is a simple container corresponding to the Les
  * Houches accord (hep-ph/0109068) common block with the same
  * name. The members are named in the same way as in the common
  * block. However, fortran arrays are represented by vectors, except
  * for the arrays of length two which are represented by pair objects.
  */
 class HEPEUP {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * Default constructor.
    */
   HEPEUP()
     : NUP(0), IDPRUP(0), XWGTUP(0.0), XPDWUP(0.0, 0.0),
       SCALUP(0.0), AQEDUP(0.0), AQCDUP(0.0) {}
   //@}
 
 public:
 
   /**
    * Set the NUP variable, corresponding to the number of particles in
    * the current event, to \a nup, and resize all relevant vectors
    * accordingly.
    */
   void resize(int nup) {
     NUP = nup;
     resize();
   }
 
   /**
    * Assuming the NUP variable, corresponding to the number of
    * particles in the current event, is correctly set, resize the
    * relevant vectors accordingly.
    */
   void resize() {
     IDUP.resize(NUP);
     ISTUP.resize(NUP);
     MOTHUP.resize(NUP);
     ICOLUP.resize(NUP);
-    PUP.resize(NUP, std::vector<double>(5));
+    PUP.resize(NUP);
     VTIMUP.resize(NUP);
     SPINUP.resize(NUP);
   }
 
   /**
    * The number of particle entries in the current event.
    */
   int NUP;
 
   /**
    * The subprocess code for this event (as given in LPRUP).
    */
   int IDPRUP;
 
   /**
    * The weight for this event.
    */
   double XWGTUP;
 
   /**
    * The PDF weights for the two incoming partons. Note that this
    * variable is not present in the current LesHouches accord
    * (hep-ph/0109068), hopefully it will be present in a future
    * accord.
    */
   pair<double,double> XPDWUP;
 
   /**
    * The scale in GeV used in the calculation of the PDF's in this
    * event.
    */
   double SCALUP;
 
   /**
    * The value of the QED coupling used in this event.
    */
   double AQEDUP;
 
   /**
    * The value of the QCD coupling used in this event.
    */
   double AQCDUP;
 
   /**
    * The PDG id's for the particle entries in this event.
    */
   vector<long> IDUP;
 
   /**
    * The status codes for the particle entries in this event.
    */
   vector<int> ISTUP;
 
   /**
    * Indices for the first and last mother for the particle entries in
    * this event.
    */
   vector< pair<int,int> > MOTHUP;
 
   /**
    * The colour-line indices (first(second) is (anti)colour) for the
    * particle entries in this event.
    */
   vector< pair<int,int> > ICOLUP;
 
   /**
    * Lab frame momentum (Px, Py, Pz, E and M in GeV) for the particle
    * entries in this event.
    */
-  vector< vector<double> > PUP;
+  vector< array<double,5> > PUP;
 
   /**
    * Invariant lifetime (c*tau, distance from production to decay im
    * mm) for the particle entries in this event.
    */
   vector<double> VTIMUP;
 
   /**
    * Spin info for the particle entries in this event given as the
    * cosine of the angle between the spin vector of a particle and the
    * 3-momentum of the decaying particle, specified in the lab frame.
    */
   vector<double> SPINUP;
 
 };
 
 }
 
 #endif /* THEPEG_LesHouches_H */
diff --git a/LesHouches/LesHouchesReader.cc b/LesHouches/LesHouchesReader.cc
--- a/LesHouches/LesHouchesReader.cc
+++ b/LesHouches/LesHouchesReader.cc
@@ -1,1592 +1,1592 @@
 // -*- C++ -*-
 //
 // LesHouchesReader.cc 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.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the LesHouchesReader class.
 //
 
 #include "LesHouchesReader.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/RefVector.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/Command.h"
 #include "config.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/PDF/PartonExtractor.h"
 #include "ThePEG/PDF/NoPDF.h"
 #include "ThePEG/Cuts/Cuts.h"
 #include "ThePEG/EventRecord/TmpTransform.h"
 #include "ThePEG/Utilities/UtilityBase.h"
 #include "ThePEG/Handlers/XComb.h"
 #include "ThePEG/Handlers/CascadeHandler.h"
 #include "LesHouchesEventHandler.h"
 #include "ThePEG/Utilities/Throw.h"
 #include "ThePEG/Utilities/HoldFlag.h"
 #include "ThePEG/Utilities/Debug.h"
 #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
 
 using namespace ThePEG;
 
 LesHouchesReader::LesHouchesReader(bool active)
   : theNEvents(0), position(0), reopened(0), theMaxScan(-1), scanning(false),
     isActive(active), theCacheFileName(""), doCutEarly(true),
     preweight(1.0), reweightPDF(false), doInitPDFs(false),
     theMaxMultCKKW(0), theMinMultCKKW(0), lastweight(1.0), maxFactor(1.0), optionalnpLO(0), optionalnpNLO(0),
     weightScale(1.0*picobarn), skipping(false), theMomentumTreatment(0),
     useWeightWarnings(true),theReOpenAllowed(true), theIncludeSpin(true) {}
 
 LesHouchesReader::LesHouchesReader(const LesHouchesReader & x)
   : HandlerBase(x), LastXCombInfo<>(x), heprup(x.heprup), hepeup(x.hepeup),
     inData(x.inData), inPDF(x.inPDF), outPDF(x.outPDF),
     thePartonExtractor(x.thePartonExtractor), thePartonBins(x.thePartonBins),
     theXCombs(x.theXCombs), theCuts(x.theCuts),
     theNEvents(x.theNEvents), position(x.position), reopened(x.reopened),
     theMaxScan(x.theMaxScan), scanning(false),
     isActive(x.isActive),
     theCacheFileName(x.theCacheFileName), doCutEarly(x.doCutEarly),
     stats(x.stats), statmap(x.statmap),
     thePartonBinInstances(x.thePartonBinInstances),
     reweights(x.reweights), preweights(x.preweights),
     preweight(x.preweight), reweightPDF(x.reweightPDF),
     doInitPDFs(x.doInitPDFs),
     theMaxMultCKKW(x.theMaxMultCKKW), theMinMultCKKW(x.theMinMultCKKW),
     lastweight(x.lastweight), maxFactor(x.maxFactor),
     weightScale(x.weightScale), xSecWeights(x.xSecWeights),
     maxWeights(x.maxWeights), skipping(x.skipping),
     theMomentumTreatment(x.theMomentumTreatment),
     useWeightWarnings(x.useWeightWarnings),
     theReOpenAllowed(x.theReOpenAllowed),
     theIncludeSpin(x.theIncludeSpin) {}
 
 LesHouchesReader::~LesHouchesReader() {}
 
 void LesHouchesReader::doinitrun() {
   HandlerBase::doinitrun();
   stats.reset();
   for ( StatMap::iterator i = statmap.begin(); i != statmap.end(); ++i )
     i->second.reset();
   open();
   if ( cacheFileName().length() ) openReadCacheFile();
   position = 0;
   reopened = 0;
 }
 
 bool LesHouchesReader::preInitialize() const {
   if ( HandlerBase::preInitialize() ) return true;
   if ( doInitPDFs && ! ( inPDF.first && inPDF.second ) ) return true;
   return false;
 }
 
 void LesHouchesReader::doinit() {
   HandlerBase::doinit();
   open();
   close();
   if ( !heprup.IDBMUP.first || !heprup.IDBMUP.second )
     Throw<LesHouchesInitError>()
       << "No information about incoming particles were found in "
       << "LesHouchesReader '" << name() << "'." << Exception::warning;
   inData = make_pair(getParticleData(heprup.IDBMUP.first),
                      getParticleData(heprup.IDBMUP.second));
   if ( heprup.EBMUP.first <= 0.0 || heprup.EBMUP.second <= 0.0 )
     Throw<LesHouchesInitError>()
     << "No information about the energy of incoming particles were found in "
     << "LesHouchesReader '" << name() << "'." << Exception::warning;
   
   if ( doInitPDFs && ! ( inPDF.first && inPDF.second ) ) {
     initPDFs();
     if ( ! ( inPDF.first && inPDF.second ) ) Throw<InitException>()
       << "LesHouchesReader '" << name()
       << "' could not create PDFBase objects in pre-initialization."
       << Exception::warning;
   }
   else if ( !inPDF.first || !inPDF.second ) Throw<LesHouchesInitError>()
     << "No information about the PDFs of incoming particles were found in "
     << "LesHouchesReader '" << name() << "'." << Exception::warning;
 }
 
 void LesHouchesReader::initPDFs() {
   if ( inPDF.first && inPDF.second ) return;
 
   string remhname;
   if ( heprup.PDFSUP.first && !inPDF.first) {
     inPDF.first = dynamic_ptr_cast<PDFPtr>
       (generator()->preinitCreate("ThePEG::LHAPDF", fullName() + "/PDFA",
                                   "ThePEGLHAPDF.so"));
     if ( !inPDF.first ) {
       Throw<InitException>()
         << "LesHouchesReader '" << name() << "' could not use information "
         << "about the PDFs used because the LHAPDF library was not properly "
         "defined." << Exception::warning;
       return;
     }
     remhname = fullName() + "/DummyRemH";
     generator()->preinitCreate("ThePEG::NoRemnants", remhname);
     generator()->preinitInterface(inPDF.first, "RemnantHandler",
                                   "set", remhname);
     if ( heprup.PDFGUP.first > 0 && heprup.PDFGUP.first < 10 ) {
       ostringstream os;
       os << heprup.PDFGUP.first << " " << heprup.PDFSUP.first;
       generator()->preinitInterface(inPDF.first, "PDFLIBNumbers",
                                     "set", os.str());
     } else {
       ostringstream os;
       os << heprup.PDFGUP.first*1000 + heprup.PDFSUP.first;
       generator()->preinitInterface(inPDF.first, "PDFNumber",
                                     "set", os.str());
     }
     generator()->preinitInterface(inPDF.first, "RangeException",
                                     "newdef", "Freeze");
   }
 
   if ( heprup.PDFSUP.second && !inPDF.second) {
     inPDF.second = dynamic_ptr_cast<PDFPtr>
       (generator()->preinitCreate("ThePEG::LHAPDF", fullName() + "/PDFB",
                                   "ThePEGLHAPDF.so"));
     if ( !inPDF.second ) {
       Throw<InitException>()
         << "LesHouchesReader '" << name() << "' could not use information "
         << "about the PDFs used because the LHAPDF library was not properly "
         "defined." << Exception::warning;
       return;
     }
     if ( remhname == "" ) {
       remhname = fullName() + "/DummyRemH";
       generator()->preinitCreate("ThePEG::NoRemnants", remhname);
     }
     generator()->preinitInterface(inPDF.second, "RemnantHandler",
                                   "set", remhname);      
     if ( heprup.PDFGUP.second > 0 && heprup.PDFGUP.second < 10 ) {
       ostringstream os;
       os << heprup.PDFGUP.second << " " << heprup.PDFSUP.second;
       generator()->preinitInterface(inPDF.second, "PDFLIBNumbers",
                                     "set", os.str());
     } else {
       ostringstream os;
       os << heprup.PDFGUP.second*1000 + heprup.PDFSUP.second;
       generator()->preinitInterface(inPDF.second, "PDFNumber",
                                     "set", os.str());
     }
     generator()->preinitInterface(inPDF.second, "RangeException",
                                     "newdef", "Freeze");
   }
   
   if ( ! ( inPDF.first && inPDF.second ) ) Throw<InitException>()
     << "LesHouchesReader '" << name()
     << "' could not find information about the PDFs used."
     << Exception::warning;
 }
 
 void LesHouchesReader::initialize(LesHouchesEventHandler & eh) {
   Energy2 Smax = ZERO;
   double Y = 0.0;
   if ( !theCuts ) {
     theCuts = eh.cuts();
     if ( !theCuts ) Throw<LesHouchesInitError>()
       << "No Cuts object was assigned to the LesHouchesReader '"
       << name() << "' nor was one\nassigned to the controlling "
       << "LesHouchesEventHandler '" << eh.name() << "'.\nAt least one of them "
       << "needs to have a Cuts object." << Exception::runerror;
 
     Smax = cuts().SMax();
     Y = cuts().Y();
   }
 
   theCKKW = eh.CKKWHandler();
 
   if ( !partonExtractor() ) {
     thePartonExtractor = eh.partonExtractor();
     if ( !partonExtractor() )  Throw<LesHouchesInitError>()
       << "No PartonExtractor object was assigned to the LesHouchesReader '"
       << name() << "' nor was one\nassigned to the controlling "
       << "LesHouchesEventHandler '" << eh.name() << "'.\nAt least one of them "
       << "needs to have a PartonExtractor object." << Exception::runerror;
   }
   open();
 
   Energy emax = 2.0*sqrt(heprup.EBMUP.first*heprup.EBMUP.second)*GeV;
   theCuts->initialize(sqr(emax),
                       0.5*log(heprup.EBMUP.first/heprup.EBMUP.second));
   if ( Smax > ZERO && ( Smax != cuts().SMax() || Y != cuts().Y() ) )
     Throw<LesHouchesInitError>()
       << "The LesHouchesReader '" << name() << "' uses the same Cuts object "
       << "as another LesHouchesReader which has not got the same energies of "
       << "the colliding particles. For the generation to work properly "
       << "different LesHouchesReader object with different colliding particles "
       << "must be assigned different (although possibly identical) Cuts "
       << "objects." << Exception::warning;
 
   thePartonBins = partonExtractor()->getPartons(emax, inData, cuts());
   for ( int i = 0, N = partonBins().size(); i < N; ++i ) {
     theXCombs[partonBins()[i]] =
       new_ptr(XComb(emax, inData, &eh, partonExtractor(), CKKWHandler(),
                     partonBins()[i], theCuts));
     partonExtractor()->nDims(partonBins()[i]);
   }
   outPDF = make_pair(partonExtractor()->getPDF(inData.first),
                      partonExtractor()->getPDF(inData.second));
 
   close();
 
   if ( !heprup.IDWTUP && useWeightWarnings )
     Throw<LesHouchesInitError>()
       << "No information about the weighting scheme was found. The events "
       << "produced by LesHouchesReader " << name()
       << " may not be sampled correctly." << Exception::warning;
 
   if ( abs(heprup.IDWTUP) > 1 && !eh.weighted() && useWeightWarnings  )
     Throw<LesHouchesInitError>()
       << "LesHouchesReader " << name() << " has the IDWTUP flag set to "
       << heprup.IDWTUP << " which is not supported by this reader, the "
       << "produced events may not be sampled correctly. It is up to "
       << "sub-classes of LesHouchesReader to correctly convert to match IDWTUP "
       << "+/- 1. Will try to make intelligent guesses to get "
       << "correct statistics.\nIn most cases this should be sufficient. "
       << "Unset <interface>WeightWarnings</interface> to avoid this message,"
       << "or set <interface>Weighted</interface> to on."
       << Exception::warning;
 
   if ( heprup.IDWTUP != eh.weightOption() && abs(heprup.IDWTUP) < 3 &&
        useWeightWarnings  )
     Throw<LesHouchesInitError>()
       << "LesHouchesReader " << name() << " has the IDWTUP flag set to "
       << heprup.IDWTUP 
       << ", which does not correspond\nto the weight option "
       << eh.weightOption() << " set in "
       << "the LesHouchesEventHandler " << eh.name() << ".\n\n"
       << "Use the following handler setting instead:\n"
       << "  set " << eh.name() << ":WeightOption " << heprup.IDWTUP
       << "\nWill try to make intelligent guesses to get "
       << "correct statistics. In most cases this should be sufficient. "
       << "Unset <interface>WeightWarnings</interface> to avoid this message"
       << Exception::warning;
 
   scan();
   initStat();
 }
 
 
 long LesHouchesReader::scan() {
   
   open();
 
   // Shall we write the events to a cache file for fast reading? If so
   // we write to a temporary file if the caches events should be
   // randomized.
   if ( cacheFileName().length() ) openWriteCacheFile();
 
   // Keep track of the number of events scanned.
   long neve = 0;
   long cuteve = 0;
   bool negw = false;
 
   // If the open() has not already gotten information about subprocesses
   // and cross sections we have to scan through the events.
   if ( !heprup.NPRUP || cacheFile() || abs(heprup.IDWTUP) != 1 ) { // why scan if IDWTUP != 1?
 
     HoldFlag<> isScanning(scanning);
 
     double oldsum = 0.0;
     vector<int> lprup;
     vector<double> newmax;
     vector<long> oldeve;
     vector<long> neweve;
     vector<double> sumlprup;
     vector<double> sumsqlprup;  
     vector<long> nscanned;
     for ( int i = 0; ( maxScan() < 0 || i < maxScan() ) && readEvent(); ++i ) {
       if ( !checkPartonBin() ) Throw<LesHouchesInitError>()
         << "Found event in LesHouchesReader '" << name()
         << "' which cannot be handeled by the assigned PartonExtractor '"
         << partonExtractor()->name() << "'." << Exception::runerror;
       vector<int>::iterator idit =
         find(lprup.begin(), lprup.end(), hepeup.IDPRUP);
       int id = lprup.size();
       if ( idit == lprup.end() ) {
         lprup.push_back(hepeup.IDPRUP);
         newmax.push_back(0.0);
         neweve.push_back(0);
         oldeve.push_back(0);
         sumlprup.push_back(0.);
         sumsqlprup.push_back(0.);
         nscanned.push_back(0);
       } else {
         id = idit - lprup.begin();
       }
       ++neve;
       ++oldeve[id];
       oldsum += hepeup.XWGTUP;
       sumlprup[id] += hepeup.XWGTUP;
       sumsqlprup[id] += sqr(hepeup.XWGTUP);
       ++nscanned[id];
       if ( cacheFile() ) {
         if ( eventWeight() == 0.0 ) {
           ++cuteve;
           continue;
         }
         cacheEvent();
       }
       ++neweve[id];
       newmax[id] = max(newmax[id], abs(eventWeight()));
       if ( eventWeight() < 0.0 ) negw = true;
     } //end of scanning events
     xSecWeights.resize(oldeve.size(), 1.0);
     for ( int i = 0, N = oldeve.size(); i < N; ++i )
       if ( oldeve[i] ) xSecWeights[i] = double(neweve[i])/double(oldeve[i]);
 
     if ( maxScan() < 0 || neve > NEvents() ) NEvents(neve - cuteve);
 
     if ( lprup.size() == heprup.LPRUP.size() ) {
       for ( int id = 0, N = lprup.size(); id < N; ++id ) {
         vector<int>::iterator idit =
           find(heprup.LPRUP.begin(), heprup.LPRUP.end(), hepeup.IDPRUP);
         if ( idit == heprup.LPRUP.end() ) {
           Throw<LesHouchesInitError>()
             << "When scanning events, the LesHouschesReader '" << name()
             << "' found undeclared processes."  << Exception::warning;
           heprup.NPRUP = 0;
           break;
         }
         int idh = idit - heprup.LPRUP.begin();
         heprup.XMAXUP[idh] = newmax[id];
       } 
     }
     if ( heprup.NPRUP == 0 ) {
       // No heprup block was supplied or something went wrong.
       heprup.NPRUP = lprup.size();
       heprup.LPRUP.resize(lprup.size());
       heprup.XMAXUP.resize(lprup.size());
       for ( int id = 0, N = lprup.size(); id < N; ++id ) {
         heprup.LPRUP[id] = lprup[id];
         heprup.XMAXUP[id] = newmax[id];
       }
     }
     if ( abs(heprup.IDWTUP) != 1 ) {
       // Try to fix things if abs(heprup.IDWTUP) != 1.
       double sumxsec = 0.0;
       if(abs(heprup.IDWTUP)==3) {
 	for ( int id = 0; id < heprup.NPRUP; ++id ) sumxsec += heprup.XSECUP[id];
       }
       else {
 	for ( int id = 0; id < heprup.NPRUP; ++id )  {
 	  //set the cross section directly from the event weights read
 	  heprup.XSECUP[id] = sumlprup[id]/nscanned[id];
 	  heprup.XERRUP[id] = (sumsqlprup[id]/nscanned[id] - sqr(sumlprup[id]/nscanned[id])) / nscanned[id];
 	  if(heprup.XERRUP[id] < 0.) {
 	    if( heprup.XERRUP[id]/(sumsqlprup[id]/nscanned[id])>-1e-10)
 	      heprup.XERRUP[id] = 0.;
 	    else {
 	      Throw<LesHouchesInitError>()
 		<< "Negative error when scanning events in LesHouschesReader '" << name()
 		<< Exception::warning;
 	      heprup.XERRUP[id] = 0.;
 	    }
 	  }
 	  heprup.XERRUP[id] = sqrt( heprup.XERRUP[id] );
 	  heprup.XMAXUP[id] = newmax[id];
 	  sumxsec += heprup.XSECUP[id];
 	}
       }
       weightScale = picobarn*neve*sumxsec/oldsum;
     }
   }
 
   if ( cacheFile() ) closeCacheFile();
 
   if ( negw ) heprup.IDWTUP = min(-abs(heprup.IDWTUP), -1);
  
   return neve;
 
 }
 
 void LesHouchesReader::setWeightScale(long) {}
 
 void LesHouchesReader::initStat() {
 
   stats.reset();
   statmap.clear();
   if ( heprup.NPRUP <= 0 ) return;
 
   double sumx = 0.0;
   xSecWeights.resize(heprup.NPRUP, 1.0);
   maxWeights.clear();
   for ( int ip = 0; ip < heprup.NPRUP; ++ip ) {
     sumx = max(heprup.XMAXUP[ip]*xSecWeights[ip], sumx);
     statmap[heprup.LPRUP[ip]] =
       XSecStat(heprup.XMAXUP[ip]*weightScale*xSecWeights[ip]);
     maxWeights[heprup.LPRUP[ip]] = heprup.XMAXUP[ip];
   }
   stats.maxXSec(sumx*weightScale);
   maxFactor = 1.0;
 }
 
 void LesHouchesReader::increaseMaxXSec(CrossSection maxxsec) {
   for ( int i = 0; i < heprup.NPRUP; ++i )
     statmap[heprup.LPRUP[i]].maxXSec(statmap[heprup.LPRUP[i]].maxXSec()*
       maxxsec/stats.maxXSec());
   maxFactor *= maxxsec/stats.maxXSec();
   stats.maxXSec(maxxsec);
 }
 
 tXCombPtr LesHouchesReader::getXComb() {
   if ( lastXCombPtr() ) return lastXCombPtr();
   fillEvent();
   connectMothers();
   tcPBPair sel = createPartonBinInstances();
   tXCombPtr lastXC = xCombs()[sel];
   // clean up the old XComb object before switching to a new one
   if ( theLastXComb && theLastXComb != lastXC ) 
     theLastXComb->clean();
   theLastXComb = lastXC;
   lastXCombPtr()->subProcess(SubProPtr());
   lastXCombPtr()->setPartonBinInstances(partonBinInstances(),
                                         sqr(hepeup.SCALUP)*GeV2);
   lastXCombPtr()->lastAlphaS(hepeup.AQCDUP);
   lastXCombPtr()->lastAlphaEM(hepeup.AQEDUP);
   return lastXCombPtr();
 }
 
 tSubProPtr LesHouchesReader::getSubProcess() {
   getXComb();
   if ( subProcess() ) return subProcess();
   lastXCombPtr()->subProcess(new_ptr(SubProcess(lastPartons(), tCollPtr(), this)));
   subProcess()->setOutgoing(outgoing().begin(), outgoing().end());
   subProcess()->setIntermediates(intermediates().begin(),
                                  intermediates().end());
   return subProcess();
 }
 
 void LesHouchesReader::fillEvent() {
   if ( !particleIndex.empty() ) return;
   particleIndex.clear();
   colourIndex.clear();
   colourIndex(0, tColinePtr());
   createParticles();
   createBeams();
     
 }
 
 void LesHouchesReader::reopen() {
   // If we didn't know how many events there were, we know now.
   if ( NEvents() <= 0 ) NEvents(position);
   ++reopened;
   // How large fraction of the events have we actually used? And how
   // large will we have used if we go through the file again?
   double frac = double(stats.attempts())/double(NEvents());
   if ( frac*double(reopened + 1)/double(reopened) > 1.0 &&
     NEvents() - stats.attempts() <
        generator()->N() - generator()->currentEventNumber() ) {
     if(theReOpenAllowed)
       generator()->logWarning(LesHouchesReopenWarning()
                               << "Reopening LesHouchesReader '" << name()
                               << "' after accessing " << stats.attempts() 
                               << " events out of "
                               << NEvents() << Exception::warning);
     else
       throw LesHouchesReopenWarning()
         << "More events requested than available in LesHouchesReader "
         << name() << Exception::runerror;
   }
   if ( cacheFile() ) {
     closeCacheFile();
     openReadCacheFile();
     if ( !uncacheEvent() ) Throw<LesHouchesReopenError>()
       << "Could not reopen LesHouchesReader '" << name()
       << "'." << Exception::runerror;
   } else {  
     close();
     open();
     if ( !readEvent() ) Throw<LesHouchesReopenError>()
       << "Could not reopen LesHouchesReader '" << name()
       << "'." << Exception::runerror;
   }
 }
 
 void LesHouchesReader::reset() {
   particleIndex.clear();
   colourIndex.clear();
   if ( theLastXComb ) theLastXComb->clean();
   theLastXComb = tXCombPtr();
 }
 
 bool LesHouchesReader::readEvent() {
 
   reset();
 
   if ( !doReadEvent() ) return false;
 
   // If we are just skipping event we do not need to reweight or do
   // anything fancy.
   if ( skipping ) { return true; }
 
   if ( cacheFile() && !scanning ) { return true;  }
 
   // Reweight according to the re- and pre-weights objects in the
   // LesHouchesReader base class.
   lastweight = reweight();
   
   if ( !reweightPDF && !cutEarly() ) return true;
   // We should try to reweight the PDFs or make early cuts here.
 
   fillEvent();
 
   double x1 = incoming().first->momentum().plus()/
     beams().first->momentum().plus();
 
   if ( reweightPDF &&
        inPDF.first && outPDF.first && inPDF.first != outPDF.first ) {
     if ( hepeup.XPDWUP.first <= 0.0 )
       hepeup.XPDWUP.first =
         inPDF.first->xfx(inData.first, incoming().first->dataPtr(),
                          sqr(hepeup.SCALUP*GeV), x1);
     double xf = outPDF.first->xfx(inData.first, incoming().first->dataPtr(),
                                   sqr(hepeup.SCALUP*GeV), x1);
     lastweight *= xf/hepeup.XPDWUP.first;
     hepeup.XPDWUP.first = xf;
   }
 
   double x2 = incoming().second->momentum().minus()/
     beams().second->momentum().minus();
 
   if ( reweightPDF &&
        inPDF.second && outPDF.second && inPDF.second != outPDF.second ) {
     if ( hepeup.XPDWUP.second <= 0.0 )
       hepeup.XPDWUP.second =
         inPDF.second->xfx(inData.second, incoming().second->dataPtr(),
                          sqr(hepeup.SCALUP*GeV), x2);
     double xf =
       outPDF.second->xfx(inData.second, incoming().second->dataPtr(),
                          sqr(hepeup.SCALUP*GeV), x2);
     lastweight *= xf/hepeup.XPDWUP.second;
     hepeup.XPDWUP.second = xf;
   }
 
   if ( cutEarly() ) {
     if ( !cuts().initSubProcess((incoming().first->momentum() +
                                  incoming().second->momentum()).m2(),
                                 0.5*log(x1/x2)) ) lastweight = 0.0;
     tSubProPtr sub = getSubProcess();
     TmpTransform<tSubProPtr> tmp(sub, Utilities::getBoostToCM(sub->incoming()));
     if ( !cuts().passCuts(*sub) ) lastweight = 0.0;
   }
 
   return true;
 }
 
 double LesHouchesReader::getEvent() {
   if ( cacheFile() ) {
     if ( !uncacheEvent() ) reopen();
   } else {
     if ( !readEvent() ) reopen();
   }
   ++position;
 
   double max = maxWeights[hepeup.IDPRUP]*maxFactor;
   // normalize all the weights to the max weight
   for(map<string,double>::iterator it=optionalWeights.begin();
       it!=optionalWeights.end();++it) {
     it->second = (max != 0.0) ? it->second/max : 0.0;
   }
   return (max != 0.0) ? eventWeight()/max : 0.0; 
 }
 
 void LesHouchesReader::skip(long n) {
   HoldFlag<> skipflag(skipping);
   while ( n-- ) getEvent();
 }
 
 double LesHouchesReader::reweight() {
   preweight = 1.0;
   if ( reweights.empty() && preweights.empty() &&
        !( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) )
     return 1.0;
   fillEvent();
   getSubProcess();
   for ( int i = 0, N = preweights.size(); i < N; ++i ) {
     preweights[i]->setXComb(lastXCombPtr());
     preweight *= preweights[i]->weight();
   }
   double weight = preweight;
   for ( int i = 0, N = reweights.size(); i < N; ++i ) {
     reweights[i]->setXComb(lastXCombPtr());
     weight *= reweights[i]->weight();
   }
 
   // If we are caching events we do not want to do CKKW reweighting.
   if ( cacheFile() ) return weight;
 
   if ( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) {
     CKKWHandler()->setXComb(lastXCombPtr());
     weight *= CKKWHandler()->reweightCKKW(minMultCKKW(), maxMultCKKW());
   }
   return weight;
 }
 
 bool LesHouchesReader::checkPartonBin() {
 
   // First find the positions of the incoming partons.
   pair< vector<int>, vector<int> > inc;
   for ( int i = 0; i < hepeup.NUP; ++i ) {
     if ( hepeup.ISTUP[i] == -9 ) {
       if ( inc.first.empty() ) inc.first.push_back(i);
       else if ( inc.second.empty() ) inc.second.push_back(i);
     }
     else if ( hepeup.ISTUP[i] == -1 ) {
       if ( inc.first.size() &&
            hepeup.MOTHUP[i].first == inc.first.back() + 1 )
         inc.first.push_back(i);
       else if ( inc.second.size() &&
                 hepeup.MOTHUP[i].first == inc.second.back() + 1 )
         inc.second.push_back(i);
       else if ( inc.first.empty() ) {
         inc.first.push_back(-1);
         inc.first.push_back(i);
       }
       else if ( inc.second.empty() ) {
         inc.second.push_back(-1);
         inc.second.push_back(i);
       }
       else if ( inc.first.size() <= inc.second.size() )
         inc.first.push_back(i);
       else
         inc.second.push_back(i);
     }
   }
 
   // Now store the corresponding id numbers
   pair< vector<long>, vector<long> > ids;
   ids.first.push_back(inc.first[0] < 0? heprup.IDBMUP.first:
                       hepeup.IDUP[inc.first[0]]);
   for ( int i = 1, N = inc.first.size(); i < N; ++i )
     ids.first.push_back(hepeup.IDUP[inc.first[i]]);
   ids.second.push_back(inc.second[0] < 0? heprup.IDBMUP.second:
                        hepeup.IDUP[inc.second[0]]);
   for ( int i = 1, N = inc.second.size(); i < N; ++i )
     ids.second.push_back(hepeup.IDUP[inc.second[i]]);
 
   // Find the correct pair of parton bins.
   PBPair pbp;
   for ( int i = 0, N = partonBins().size(); i < N; ++i ) {
     tcPBPtr curr = partonBins()[i].first;
     int icurr = inc.first.size() - 1;
     while ( curr && icurr >= 0 ) {
       if ( curr->parton()->id () != ids.first[icurr] ) break;
       curr = curr->incoming();
       --icurr;
     }
     if(!(!partonBins()[i].first->incoming() &&
          !partonBins()[i].first->particle() &&  
          partonBins()[i].first->parton()->id () == ids.first[0] &&
          ( inc.first.size()==1 ||
            (inc.first.size()==2 && ids.first[0]==ids.first[1]))) &&
        ( curr || icurr >= 0 ) ) continue;
 
     curr = partonBins()[i].second;
     icurr = inc.second.size() - 1;
     while ( curr && icurr >= 0 ) {
       if ( curr->parton()->id () != ids.second[icurr] ) break;
       curr = curr->incoming();
       --icurr;
     }
     if(!(!partonBins()[i].second->incoming() &&
          !partonBins()[i].second->particle() &&  
          partonBins()[i].second->parton()->id () == ids.second[0] &&
          ( inc.second.size()==1 ||
            (inc.second.size()==2 && ids.second[0]==ids.second[1]))) &&
        ( curr || icurr >= 0 ) ) continue;
 
     pbp = partonBins()[i];
   }
 
   // If we are only checking we return here.
   return ( pbp.first && pbp.second );
 
 }
 
 namespace {
   bool recursionNotNull(tcPBPtr bin, tcPPtr p) {
     while ( bin && p ) {
       if ( p->dataPtr() != bin->parton() ) break;
       bin = bin->incoming();
       p = p->parents().size()? p->parents()[0]: tPPtr();
     }
     return bin || p;
   }
 }
 
 
 tcPBPair LesHouchesReader::createPartonBinInstances() {
   tcPBPair sel;
   for ( int i = 0, N = partonBins().size(); i < N; ++i ) {
     tcPBPtr bin = partonBins()[i].first;
     tcPPtr p = incoming().first;
     if ( recursionNotNull(bin,p) ) continue;
     bin = partonBins()[i].second;
     p = incoming().second;
     if ( recursionNotNull(bin,p) ) continue;
     sel = partonBins()[i];
     break;
   }
   if ( !sel.first || !sel.second ) Throw<LesHouchesInconsistencyError>()
     << "Could not find appropriate PartonBin objects for event produced by "
     << "LesHouchesReader '" << name() << "'." << Exception::runerror;
 
   Direction<0> dir(true);
   thePartonBinInstances.first =
     new_ptr(PartonBinInstance(incoming().first, sel.first,
                               -sqr(hepeup.SCALUP*GeV)));
   if ( thePartonBinInstances.first->xi() > 1.00001 ) {
     Throw<LesHouchesInconsistencyError>()
       << "Found an event with momentum fraction larger than unity (x1="
       << thePartonBinInstances.first->xi()
       << "). The event will be skipped." << Exception::warning;
     throw Veto();
   }
   dir.reverse();
   thePartonBinInstances.second =
     new_ptr(PartonBinInstance(incoming().second, sel.second,
                               -sqr(hepeup.SCALUP*GeV)));
 
   if ( thePartonBinInstances.second->xi() > 1.00001 ) {
     Throw<LesHouchesInconsistencyError>()
       << "Found an event with momentum fraction larger than unity (x2="
       << thePartonBinInstances.second->xi()
       << "). The event will be skipped." << Exception::warning;
     throw Veto();
   }
   return sel;
 
 }
                      
 
 void LesHouchesReader::createParticles() {
   theBeams = PPair();
   theIncoming = PPair();
   theOutgoing = PVector();
   theIntermediates = PVector();
   for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
     if ( !hepeup.IDUP[i] ) continue;
     Lorentz5Momentum mom(hepeup.PUP[i][0]*GeV, hepeup.PUP[i][1]*GeV,
                          hepeup.PUP[i][2]*GeV, hepeup.PUP[i][3]*GeV,
                          hepeup.PUP[i][4]*GeV);
     //   cout << hepeup.PUP[i][0] << " " << hepeup.PUP[i][1] << " " << hepeup.PUP[i][2] << " " << hepeup.PUP[i][3] << " " << hepeup.PUP[i][4] << endl;
     if(theMomentumTreatment == 1)      mom.rescaleEnergy();
     else if(theMomentumTreatment == 2) mom.rescaleMass();
     // cout << hepeup.PUP[i][0] << " " << hepeup.PUP[i][1] << " " << hepeup.PUP[i][2] << " " << hepeup.PUP[i][3] << " " << hepeup.PUP[i][4] << endl;
     PDPtr pd = getParticleData(hepeup.IDUP[i]);
     if (!pd) {
       Throw<LesHouchesInitError>()
         << "LesHouchesReader '" << name() << "' found unknown particle ID "
         << hepeup.IDUP[i]
         << " in Les Houches common block structure.\n"
         << "You need to define the new particle in an input file.\n"
         << Exception::runerror;
     }
     if ( ! pd->coloured() 
          && ( hepeup.ICOLUP[i].first != 0 || hepeup.ICOLUP[i].second != 0 ) ) {
       Throw<LesHouchesInconsistencyError>()
         << "LesHouchesReader " << name() << ": " << pd->PDGName() 
         << " is not a coloured particle.\nIt should not have "
         << "(anti-)colour lines " << hepeup.ICOLUP[i].first
         << ' ' << hepeup.ICOLUP[i].second
         << " set; the event file needs to be fixed."
         << Exception::runerror;
     }
     PPtr p = pd->produceParticle(mom);
     if(hepeup.ICOLUP[i].first>=0 && hepeup.ICOLUP[i].second >=0) {
       tColinePtr c = colourIndex(hepeup.ICOLUP[i].first);
       if ( c ) {
         c->addColoured(p);
       }
       c = colourIndex(hepeup.ICOLUP[i].second);
       if ( c ) c->addAntiColoured(p);
     }
     else {
       tColinePtr c1 = colourIndex(abs(hepeup.ICOLUP[i].first ));
       tColinePtr c2 = colourIndex(abs(hepeup.ICOLUP[i].second));
       if(pd->hasColour()) {
         c1->addColouredIndexed(p,1);
         c2->addColouredIndexed(p,2);
       }
       else {
         c1->addAntiColouredIndexed(p,1);
         c2->addAntiColouredIndexed(p,2);
       }
     }
     particleIndex(i + 1, p);
     switch ( hepeup.ISTUP[i] ) {
     case -9:
       if ( !theBeams.first ) theBeams.first = p;
       else if ( !theBeams.second ) theBeams.second = p;
       else Throw<LesHouchesInconsistencyError>()
         << "To many incoming beam particles in the LesHouchesReader '"
         << name() << "'." << Exception::runerror;
       break;
     case -1:
       if ( !theIncoming.first ) theIncoming.first = p;
       else if ( !theIncoming.second ) theIncoming.second = p;
       else if ( particleIndex(theIncoming.first) == hepeup.MOTHUP[i].first )
         theIncoming.first = p;
       else if ( particleIndex(theIncoming.second) == hepeup.MOTHUP[i].first )
         theIncoming.second = p;
       else Throw<LesHouchesInconsistencyError>()
         << "To many incoming particles to hard subprocess in the "
         << "LesHouchesReader '" << name() << "'." << Exception::runerror;
       p->scale(sqr(hepeup.SCALUP*GeV));
       break;
     case 1:
       theOutgoing.push_back(p);
       p->scale(sqr(hepeup.SCALUP*GeV));
       break;
     case -2:
     case 2:
     case 3:
       theIntermediates.push_back(p);
       break;
     default:
       Throw<LesHouchesInconsistencyError>()
         << "Unknown status code (" << hepeup.ISTUP[i]
         << ") in the LesHouchesReader '" << name() << "'."
         << Exception::runerror;
     }
 
     // value 9 is defined as "Unknown or unpolarized particles"
     double spinup = hepeup.SPINUP[i];
     if ( abs(spinup - 9) < 1.0e-3 )
       spinup = 0.;
     if ( spinup < -1. || spinup > 1. ) {
       Throw<LesHouchesInconsistencyError>()
         << "Polarization must be between -1 and 1, not "
         << spinup << " as found in the "
         << "LesHouches event file.\nThe event file needs to be fixed." 
         << Exception::runerror;
     }
     if( theIncludeSpin 
         && abs(pd->id()) == ParticleID::tauminus 
         && spinup !=0) {
       if(pd->iSpin() == PDT::Spin1Half ) {
         vector<Helicity::SpinorWaveFunction> wave;
         Helicity::SpinorWaveFunction(wave,p,Helicity::outgoing,true);
         RhoDMatrix rho(pd->iSpin(),true);
         rho(0,0) = 0.5*(1.-spinup);
         rho(1,1) = 0.5*(1.+spinup);
         p->spinInfo()->rhoMatrix() = rho;
         p->spinInfo()->  DMatrix() = rho;
       }
     }
   }
 
   // check the colour flows, and if necessary create any sources/sinks
   // hard process
   // get the particles in the hard process
   PVector external;
   for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
     unsigned int moth; 
     switch ( hepeup.ISTUP[i] ) {
     case -1:
       external.push_back(particleIndex.find(i+1));
       break;
     case 1: case 2: case 3:
       moth = hepeup.MOTHUP[i].first;
       if(moth!=0 && (hepeup.ISTUP[moth]==-1||hepeup.ISTUP[moth]==-2||
                    hepeup.ISTUP[moth]==-9))
         external.push_back(particleIndex.find(i+1));
       moth = hepeup.MOTHUP[i].second;
       if(moth!=0 && (hepeup.ISTUP[moth]==-1||hepeup.ISTUP[moth]==-2||
                      hepeup.ISTUP[moth]==-9))
         external.push_back(particleIndex.find(i+1));
       break;
     case -2: case -9: default: 
       break;
     }
   }
   // check the incoming/outgoing lines match
   vector<tColinePtr> unMatchedColour,unMatchedAntiColour;
   for(unsigned int ix=0;ix<external.size();++ix) {
     vector<tcColinePtr> 
       col  = external[ix]->colourInfo()->    colourLines();
     vector<tcColinePtr> 
       anti = external[ix]->colourInfo()->antiColourLines();
     if(hepeup.ISTUP[particleIndex(external[ix])-1]<0)
       swap(col,anti);
     if(!col.empty()) {
       for(unsigned int ic1=0;ic1<col.size();++ic1) {
         bool matched=false;
         for(unsigned int iy=0;iy<external.size();++iy) {
           vector<tcColinePtr> col2;
           if(hepeup.ISTUP[particleIndex(external[iy])-1]<0) {
             if(external[iy]->colourInfo()->colourLines().empty()) continue;
             col2 = external[iy]->colourInfo()->colourLines();
           } 
           else if(hepeup.ISTUP[particleIndex(external[iy])-1]>0) {
             if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
             col2 = external[iy]->colourInfo()->antiColourLines();
           }
           for(unsigned int ic2=0;ic2<col2.size();++ic2) {
             if(col[ic1]==col2[ic2]) {
               matched=true;
               break;
             }
           }
           if(matched) break;
         }
         if(!matched) unMatchedColour.push_back(const_ptr_cast<tColinePtr>(col[ic1]));
       }
     }
     if(!anti.empty()) {
       for(unsigned int ic1=0;ic1<col.size();++ic1) {
         bool matched=false;
         for(unsigned int iy=0;iy<external.size();++iy) {
           vector<tcColinePtr> anti2;
           if(hepeup.ISTUP[particleIndex(external[iy])-1]<0) {
             if(external[iy]->colourInfo()->colourLines().empty()) continue;
             anti2 = external[iy]->colourInfo()->antiColourLines();
           } 
           else if(hepeup.ISTUP[particleIndex(external[iy])-1]>0) {
             if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
             anti2 = external[iy]->colourInfo()->colourLines();
           }
           for(unsigned int ic2=0;ic2<anti2.size();++ic2) {
             if(col[ic1]==anti2[ic2]) {
               matched=true;
               break;
             }
           }
           if(matched) break;
         }
         if(!matched) unMatchedAntiColour.push_back(const_ptr_cast<tColinePtr>(anti[ic1]));
       }
     }
   }
     
   // might have source/sink
   if( unMatchedColour.size() + unMatchedAntiColour.size() != 0) {
     if(unMatchedColour.size() == 3 ) {
       unMatchedColour[0]->setSourceNeighbours(unMatchedColour[1],
                                               unMatchedColour[2]);
     }
     else if(unMatchedColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
       Throw<LesHouchesInconsistencyError>()
         << "LesHouchesReader '" << name() << "' found inconsistent colour "
         << "flow in Les Houches common block structure for hard process.\n"
         << hepeup << Exception::runerror;
     }
     if(unMatchedAntiColour.size() == 3 ) {
       unMatchedAntiColour[0]->setSinkNeighbours(unMatchedAntiColour[1],
                                                 unMatchedAntiColour[2]);
     }
     else if(unMatchedAntiColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
       Throw<LesHouchesInconsistencyError>()
         << "LesHouchesReader '" << name() << "' found inconsistent colour "
         << "flow in Les Houches common block structure for hard process.\n"
         << hepeup << Exception::runerror;
     }
   }
 
   // any subsequent decays
   for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
     if(hepeup.ISTUP[i] !=2 && hepeup.ISTUP[i] !=3) continue;
     PVector external;
     external.push_back(particleIndex.find(i+1));
     for ( int j = 0; j < N; ++j ) {
       if(hepeup.MOTHUP[j].first==i+1||  hepeup.MOTHUP[j].second==i+1)
         external.push_back(particleIndex.find(j+1));
     }
     // check the incoming/outgoing lines match
     vector<tColinePtr> unMatchedColour,unMatchedAntiColour;
     for(unsigned int ix=0;ix<external.size();++ix) {
       vector<tcColinePtr> 
         col  = external[ix]->colourInfo()->    colourLines();
       vector<tcColinePtr> 
         anti = external[ix]->colourInfo()->antiColourLines();
       if(ix==0) swap(col,anti);
       if(!col.empty()) {
         for(unsigned int ic1=0;ic1<col.size();++ic1) {
           bool matched=false;
           for(unsigned int iy=0;iy<external.size();++iy) {
             if(iy==ix) continue;
             vector<tcColinePtr> col2;
             if(iy==0) {
               if(external[iy]->colourInfo()->colourLines().empty()) continue;
               col2 = external[iy]->colourInfo()->colourLines();
             } 
             else {
               if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
               col2 = external[iy]->colourInfo()->antiColourLines();
             }
             for(unsigned int ic2=0;ic2<col2.size();++ic2) {
               if(col[ic1]==col2[ic2]) {
                 matched=true;
                 break;
               }
             }
             if(matched) break;
           }
           if(!matched) unMatchedColour.push_back(const_ptr_cast<tColinePtr>(col[ic1]));
         }
       }
       if(!anti.empty()) {
         for(unsigned int ic1=0;ic1<anti.size();++ic1) {
           bool matched=false;
           for(unsigned int iy=0;iy<external.size();++iy) {
             if(iy==ix) continue;
             vector<tcColinePtr> anti2;
             if(iy==0) {
               if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
               anti2 = external[iy]->colourInfo()->antiColourLines();
             } 
             else {
               if(external[iy]->colourInfo()->colourLines().empty()) continue;
               anti2 = external[iy]->colourInfo()->colourLines();
             }
             for(unsigned int ic2=0;ic2<anti2.size();++ic2) {
               if(anti[ic1]==anti2[ic2]) {
                 matched=true;
                 break;
               }
             }
             if(matched) break;
           }
           if(!matched) unMatchedAntiColour.push_back(const_ptr_cast<tColinePtr>(anti[ic1]));
         }
       }
     }
     // might have source/sink
     if( unMatchedColour.size() + unMatchedAntiColour.size() != 0) {
       if(unMatchedColour.size() == 3 ) {
         unMatchedColour[0]->setSourceNeighbours(unMatchedColour[1],
                                                 unMatchedColour[2]);
       }
       else if(unMatchedColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
         Throw<LesHouchesInconsistencyError>()
           << "LesHouchesReader '" << name() << "' found inconsistent colour "
           << "flow in Les Houches common block structure for decay of \n"
           << *external[0] << "\n"
           << hepeup << Exception::runerror;
       }
       if(unMatchedAntiColour.size() == 3 ) {
         unMatchedAntiColour[0]->setSinkNeighbours(unMatchedAntiColour[1],
                                                   unMatchedAntiColour[2]);
       }
       else if(unMatchedAntiColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
         Throw<LesHouchesInconsistencyError>()
           << "LesHouchesReader '" << name() << "' found inconsistent colour "
           << "flow in Les Houches common block structure for decay of\n"
           << *external[0] << "\n"
           << hepeup << Exception::runerror;
       }
     }
   }
 
 }
 
 void LesHouchesReader::createBeams() {
 
   if ( !theBeams.first && dynamic_ptr_cast<Ptr<NoPDF>::tcp>(inPDF.first) ) {
     theBeams.first = theIncoming.first;
   }
   else if ( !theBeams.first ) {
     theBeams.first = getParticleData(heprup.IDBMUP.first)->produceParticle();
     double m = theBeams.first->mass()/GeV;
     theBeams.first->set5Momentum
       (Lorentz5Momentum(ZERO, ZERO,
                         sqrt(sqr(heprup.EBMUP.first) - sqr(m))*GeV,
                         heprup.EBMUP.first*GeV, m*GeV));
     hepeup.IDUP.push_back(heprup.IDBMUP.first);
     hepeup.ISTUP.push_back(-9);
     hepeup.MOTHUP.push_back(make_pair(0, 0));
     hepeup.ICOLUP.push_back(make_pair(0, 0));
     hepeup.VTIMUP.push_back(0.0);
     hepeup.SPINUP.push_back(0.0);
     particleIndex(hepeup.IDUP.size(), theBeams.first);
     hepeup.MOTHUP[particleIndex(theIncoming.first) - 1].first =
       hepeup.IDUP.size();
   }
   if ( !theBeams.second && dynamic_ptr_cast<Ptr<NoPDF>::tcp>(inPDF.second) ) {
     theBeams.second = theIncoming.second;
   }
   else if ( !theBeams.second ) {
     theBeams.second = getParticleData(heprup.IDBMUP.second)->produceParticle();
     double m = theBeams.second->mass()/GeV;
     theBeams.second->set5Momentum
       (Lorentz5Momentum(ZERO, ZERO,
                         -sqrt(sqr(heprup.EBMUP.second) - sqr(m))*GeV,
                         heprup.EBMUP.second*GeV, m*GeV));
     hepeup.IDUP.push_back(heprup.IDBMUP.second);
     hepeup.ISTUP.push_back(-9);
     hepeup.MOTHUP.push_back(make_pair(0, 0));
     hepeup.ICOLUP.push_back(make_pair(0, 0));
     hepeup.VTIMUP.push_back(0.0);
     hepeup.SPINUP.push_back(0.0);
     particleIndex(hepeup.IDUP.size(), theBeams.second);
     hepeup.MOTHUP[particleIndex(theIncoming.second) - 1].first =
       hepeup.IDUP.size();
   }
 }
 
 void LesHouchesReader::connectMothers() {
   const ObjectIndexer<long,Particle> & pi = particleIndex;
   for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
     if ( pi(hepeup.MOTHUP[i].first) ) 
       pi(hepeup.MOTHUP[i].first)->addChild(pi(i + 1));
     if ( pi(hepeup.MOTHUP[i].second) 
          && hepeup.MOTHUP[i].second != hepeup.MOTHUP[i].first ) 
       pi(hepeup.MOTHUP[i].second)->addChild(pi(i + 1));
   }
 }
 
 void LesHouchesReader::openReadCacheFile() {
   if ( cacheFile() ) closeCacheFile();
   cacheFile().open(cacheFileName(), "r");
   position = 0;
 }
 
 void LesHouchesReader::openWriteCacheFile() {
   if ( cacheFile() ) closeCacheFile();
   cacheFile().open(cacheFileName(), "w");
 }
 
 void LesHouchesReader::closeCacheFile() {
   cacheFile().close();
 }
 
 void LesHouchesReader::cacheEvent() const {
   static vector<char> buff;
   cacheFile().write(&hepeup.NUP, sizeof(hepeup.NUP));
   buff.resize(eventSize(hepeup.NUP));
   char * pos = &buff[0];
   pos = mwrite(pos, hepeup.IDPRUP);
   pos = mwrite(pos, hepeup.XWGTUP);
   pos = mwrite(pos, hepeup.XPDWUP);
   pos = mwrite(pos, hepeup.SCALUP);
   pos = mwrite(pos, hepeup.AQEDUP);
   pos = mwrite(pos, hepeup.AQCDUP);
   pos = mwrite(pos, hepeup.IDUP[0], hepeup.NUP);
   pos = mwrite(pos, hepeup.ISTUP[0], hepeup.NUP);
   pos = mwrite(pos, hepeup.MOTHUP[0], hepeup.NUP);
   pos = mwrite(pos, hepeup.ICOLUP[0], hepeup.NUP);
   for ( int i = 0; i < hepeup.NUP; ++i )
     pos = mwrite(pos, hepeup.PUP[i][0], 5);
   pos = mwrite(pos, hepeup.VTIMUP[0], hepeup.NUP);
   pos = mwrite(pos, hepeup.SPINUP[0], hepeup.NUP);
   pos = mwrite(pos, lastweight);
   pos = mwrite(pos, optionalWeights);
   for(size_t ff = 0; ff < optionalWeightsNames.size(); ff++) {
     pos = mwrite(pos, optionalWeightsNames[ff]);
   }
   pos = mwrite(pos, optionalnpLO);
   pos = mwrite(pos, optionalnpNLO);
   pos = mwrite(pos, preweight);
   cacheFile().write(&buff[0], buff.size(), 1);
 }
 
 bool LesHouchesReader::uncacheEvent() {
   reset();
   static vector<char> buff;
   if ( cacheFile().read(&hepeup.NUP, sizeof(hepeup.NUP)) != 1 )
     return false;
   buff.resize(eventSize(hepeup.NUP));
   if ( cacheFile().read(&buff[0], buff.size()) != 1 ) return false;
   const char * pos = &buff[0];
   pos = mread(pos, hepeup.IDPRUP);
   pos = mread(pos, hepeup.XWGTUP);
   pos = mread(pos, hepeup.XPDWUP);
   pos = mread(pos, hepeup.SCALUP);
   pos = mread(pos, hepeup.AQEDUP);
   pos = mread(pos, hepeup.AQCDUP);
   hepeup.IDUP.resize(hepeup.NUP);
   pos = mread(pos, hepeup.IDUP[0], hepeup.NUP);
   hepeup.ISTUP.resize(hepeup.NUP);
   pos = mread(pos, hepeup.ISTUP[0], hepeup.NUP);
   hepeup.MOTHUP.resize(hepeup.NUP);
   pos = mread(pos, hepeup.MOTHUP[0], hepeup.NUP);
   hepeup.ICOLUP.resize(hepeup.NUP);
   pos = mread(pos, hepeup.ICOLUP[0], hepeup.NUP);
-  hepeup.PUP.resize(hepeup.NUP, vector<double>(5));
+  hepeup.PUP.resize(hepeup.NUP);
   for ( int i = 0; i < hepeup.NUP; ++i ) 
     pos = mread(pos, hepeup.PUP[i][0], 5);
   hepeup.VTIMUP.resize(hepeup.NUP);
   pos = mread(pos, hepeup.VTIMUP[0], hepeup.NUP);
   hepeup.SPINUP.resize(hepeup.NUP);
   pos = mread(pos, hepeup.SPINUP[0], hepeup.NUP);
   pos = mread(pos, lastweight);
   pos = mread(pos, optionalWeights);
   for(size_t ff = 0; ff < optionalWeightsNames.size(); ff++) {
     pos = mread(pos, optionalWeightsNames[ff]);
   }
   pos = mread(pos, optionalnpLO);
   pos = mread(pos, optionalnpNLO);
   pos = mread(pos, preweight);
 
   // If we are skipping, we do not have to do anything else.
   if ( skipping ) return true;
 
   if ( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) {
     // The cached event has not been submitted to CKKW reweighting, so
     // we do that now.
     fillEvent();
     getSubProcess();
     CKKWHandler()->setXComb(lastXCombPtr());
     lastweight *= CKKWHandler()->reweightCKKW(minMultCKKW(), maxMultCKKW());
   }
   return true;
 }
 
 void LesHouchesReader::persistentOutput(PersistentOStream & os) const {
   os << heprup.IDBMUP << heprup.EBMUP << heprup.PDFGUP << heprup.PDFSUP
      << heprup.IDWTUP << heprup.NPRUP << heprup.XSECUP << heprup.XERRUP
      << heprup.XMAXUP << heprup.LPRUP << hepeup.NUP << hepeup.IDPRUP
      << hepeup.XWGTUP << hepeup.XPDWUP << hepeup.SCALUP << hepeup.AQEDUP
      << hepeup.AQCDUP << hepeup.IDUP << hepeup.ISTUP << hepeup.MOTHUP
      << hepeup.ICOLUP << hepeup.PUP << hepeup.VTIMUP << hepeup.SPINUP
      << inData << inPDF << outPDF << thePartonExtractor << theCKKW
      << thePartonBins << theXCombs << theCuts << theNEvents << position
      << reopened << theMaxScan << isActive
      << theCacheFileName << doCutEarly << stats << statmap
      << thePartonBinInstances
      << theBeams << theIncoming << theOutgoing << theIntermediates
      << reweights << preweights << preweight << reweightPDF << doInitPDFs
      << theLastXComb << theMaxMultCKKW << theMinMultCKKW << lastweight << optionalWeights << optionalnpLO << optionalnpNLO
      << maxFactor << ounit(weightScale, picobarn) << xSecWeights << maxWeights
      << theMomentumTreatment << useWeightWarnings << theReOpenAllowed
      << theIncludeSpin;
 }
 
 void LesHouchesReader::persistentInput(PersistentIStream & is, int) {
   if ( cacheFile() ) closeCacheFile();
   is >> heprup.IDBMUP >> heprup.EBMUP >> heprup.PDFGUP >> heprup.PDFSUP
      >> heprup.IDWTUP >> heprup.NPRUP >> heprup.XSECUP >> heprup.XERRUP
      >> heprup.XMAXUP >> heprup.LPRUP >> hepeup.NUP >> hepeup.IDPRUP
      >> hepeup.XWGTUP >> hepeup.XPDWUP >> hepeup.SCALUP >> hepeup.AQEDUP
      >> hepeup.AQCDUP >> hepeup.IDUP >> hepeup.ISTUP >> hepeup.MOTHUP
      >> hepeup.ICOLUP >> hepeup.PUP >> hepeup.VTIMUP >> hepeup.SPINUP
      >> inData >> inPDF >> outPDF >> thePartonExtractor >> theCKKW
      >> thePartonBins >> theXCombs >> theCuts >> theNEvents >> position
      >> reopened >> theMaxScan >> isActive
      >> theCacheFileName >> doCutEarly >> stats >> statmap
      >> thePartonBinInstances
      >> theBeams >> theIncoming >> theOutgoing >> theIntermediates
      >> reweights >> preweights >> preweight >> reweightPDF >> doInitPDFs
      >> theLastXComb >> theMaxMultCKKW >> theMinMultCKKW >> lastweight >> optionalWeights >> optionalnpLO >> optionalnpNLO
      >> maxFactor >> iunit(weightScale, picobarn) >> xSecWeights >> maxWeights
      >> theMomentumTreatment >> useWeightWarnings >> theReOpenAllowed
      >> theIncludeSpin;
 }
 
 AbstractClassDescription<LesHouchesReader>
 LesHouchesReader::initLesHouchesReader;
 // Definition of the static class description member.
 
 void LesHouchesReader::setBeamA(long id) { heprup.IDBMUP.first = id; }
 long LesHouchesReader::getBeamA() const { return heprup.IDBMUP.first; }
 void LesHouchesReader::setBeamB(long id) { heprup.IDBMUP.second = id; }
 long LesHouchesReader::getBeamB() const { return heprup.IDBMUP.second; }
 void LesHouchesReader::setEBeamA(Energy e) { heprup.EBMUP.first = e/GeV; }
 Energy LesHouchesReader::getEBeamA() const { return heprup.EBMUP.first*GeV; }
 void LesHouchesReader::setEBeamB(Energy e) { heprup.EBMUP.second = e/GeV; }
 Energy LesHouchesReader::getEBeamB() const { return heprup.EBMUP.second*GeV; }
 void LesHouchesReader::setPDFA(PDFPtr pdf) { inPDF.first = pdf; }
 PDFPtr LesHouchesReader::getPDFA() const { return inPDF.first; }
 void LesHouchesReader::setPDFB(PDFPtr pdf) { inPDF.second = pdf; }
 PDFPtr LesHouchesReader::getPDFB() const { return inPDF.second; }
 
 void LesHouchesReader::Init() {
 
   static ClassDocumentation<LesHouchesReader> documentation
     ("ThePEG::LesHouchesReader is an abstract base class to be used "
      "for objects which reads event files or streams from matrix element "
      "generators.");
 
   static Parameter<LesHouchesReader,long> interfaceBeamA
     ("BeamA",
      "The PDG id of the incoming particle along the positive z-axis. "
      "If zero the corresponding information is to be deduced from the "
      "event stream/file.",
      0, 0, 0, 0,
      true, false, false,
      &LesHouchesReader::setBeamA,
      &LesHouchesReader::getBeamA, 0, 0, 0);
 
   static Parameter<LesHouchesReader,long> interfaceBeamB
     ("BeamB",
      "The PDG id of the incoming particle along the negative z-axis. "
      "If zero the corresponding information is to be deduced from the "
      "event stream/file.",
      0, 0, 0, 0,
      true, false, false,
      &LesHouchesReader::setBeamB,
      &LesHouchesReader::getBeamB, 0, 0, 0);
 
   static Parameter<LesHouchesReader,Energy> interfaceEBeamA
     ("EBeamA",
      "The energy of the incoming particle along the positive z-axis. "
      "If zero the corresponding information is to be deduced from the "
      "event stream/file.",
      0, GeV, ZERO, ZERO, 1000000000.0*GeV,
      true, false, true,
      &LesHouchesReader::setEBeamA, &LesHouchesReader::getEBeamA, 0, 0, 0);
 
   static Parameter<LesHouchesReader,Energy> interfaceEBeamB
     ("EBeamB",
      "The energy of the incoming particle along the negative z-axis. "
      "If zero the corresponding information is to be deduced from the "
      "event stream/file.",
      0, GeV, ZERO, ZERO, 1000000000.0*GeV,
      true, false, true,
      &LesHouchesReader::setEBeamB, &LesHouchesReader::getEBeamB, 0, 0, 0);
 
   static Reference<LesHouchesReader,PDFBase> interfacePDFA
     ("PDFA",
      "The PDF used for incoming particle along the positive z-axis. "
      "If null the corresponding information is to be deduced from the "
      "event stream/file.",
      0, true, false, true, true, false,
      &LesHouchesReader::setPDFA, &LesHouchesReader::getPDFA, 0);
 
   static Reference<LesHouchesReader,PDFBase> interfacePDFB
     ("PDFB",
      "The PDF used for incoming particle along the negative z-axis. "
      "If null the corresponding information is to be deduced from the "
      "event stream/file.",
      0, true, false, true, true, false,
      &LesHouchesReader::setPDFB, &LesHouchesReader::getPDFB, 0);
 
   static Parameter<LesHouchesReader,long> interfaceMaxScan
     ("MaxScan",
      "The maximum number of events to scan to obtain information about "
      "processes and cross section in the intialization.",
      &LesHouchesReader::theMaxScan, -1, 0, 0,
      true, false, false);
 
   static Parameter<LesHouchesReader,string> interfaceCacheFileName
     ("CacheFileName",
      "Name of file used to cache the events from the reader in a fast-readable "
      "form. If empty, no cache file will be generated.",
      &LesHouchesReader::theCacheFileName, "",
      true, false);
   interfaceCacheFileName.fileType();
 
   static Switch<LesHouchesReader,bool> interfaceCutEarly
     ("CutEarly",
      "Determines whether to apply cuts to events before converting to "
      "ThePEG format.",
      &LesHouchesReader::doCutEarly, true, true, false);
   static SwitchOption interfaceCutEarlyYes
     (interfaceCutEarly,
      "Yes",
      "Event are cut before converted.",
      true);
   static SwitchOption interfaceCutEarlyNo
     (interfaceCutEarly,
      "No",
      "Events are not cut before converted.",
      false);
 
   static Reference<LesHouchesReader,PartonExtractor> interfacePartonExtractor
     ("PartonExtractor",
      "The PartonExtractor object used to construct remnants. If no object is "
      "provided the LesHouchesEventHandler object must provide one instead.",
      &LesHouchesReader::thePartonExtractor, true, false, true, true, false);
 
 
   static Reference<LesHouchesReader,Cuts> interfaceCuts
     ("Cuts",
      "The Cuts object to be used for this reader. Note that these "
      "must not be looser cuts than those used in the actual generation. "
      "If no object is provided the LesHouchesEventHandler object must "
      "provide one instead.",
      &LesHouchesReader::theCuts, true, false, true, true, false);
 
   static RefVector<LesHouchesReader,ReweightBase> interfaceReweights
     ("Reweights",
      "A list of ThePEG::ReweightBase objects to modify this the weight of "
      "this reader.",
      &LesHouchesReader::reweights, 0, false, false, true, false);
 
   static RefVector<LesHouchesReader,ReweightBase> interfacePreweights
     ("Preweights",
      "A list of ThePEG::ReweightBase objects to bias the phase space for this "
      "reader without influencing the actual cross section.",
      &LesHouchesReader::preweights, 0, false, false, true, false);
 
   static Switch<LesHouchesReader,bool> interfaceReweightPDF
     ("ReweightPDF",
      "If the PDFs used in the generation for this reader is different "
      "from the ones assumed by the associated PartonExtractor object, "
      "should the events be reweighted to fit the latter?",
      &LesHouchesReader::reweightPDF, false, true, false);
   static SwitchOption interfaceReweightPDFNo
     (interfaceReweightPDF, "No", "The event weights are kept as they are.",
      false);
   static SwitchOption interfaceReweightPDFYes
     (interfaceReweightPDF,
      "Yes", "The events are reweighted.", true);
 
   static Switch<LesHouchesReader,bool> interfaceInitPDFs
     ("InitPDFs",
      "If no PDFs were specified in <interface>PDFA</interface> or "
      "<interface>PDFB</interface>for this reader, try to extract the "
      "information from the event file and assign the relevant PDFBase"
      "objects when the reader is initialized.",
      &LesHouchesReader::doInitPDFs, false, true, false);
   static SwitchOption interfaceInitPDFsYes
     (interfaceInitPDFs,
      "Yes",
      "Extract PDFs during initialization.",
      true);
   static SwitchOption interfaceInitPDFsNo
     (interfaceInitPDFs,
      "No",
      "Do not extract PDFs during initialization.",
      false);
 
   static Parameter<LesHouchesReader,int> interfaceMaxMultCKKW
     ("MaxMultCKKW",
      "If this reader is to be used (possibly together with others) for CKKW-"
      "reweighting and veto, this should give the multiplicity of outgoing "
      "particles in the highest multiplicity matrix element in the group. "
      "If set to zero, no CKKW procedure should be applied.",
      &LesHouchesReader::theMaxMultCKKW, 0, 0, 0,
      true, false, Interface::lowerlim);
 
   static Parameter<LesHouchesReader,int> interfaceMinMultCKKW
     ("MinMultCKKW",
      "If this reader is to be used (possibly together with others) for CKKW-"
      "reweighting and veto, this should give the multiplicity of outgoing "
      "particles in the lowest multiplicity matrix element in the group. If "
      "larger or equal to <interface>MaxMultCKKW</interface>, no CKKW "
      "procedure should be applied.",
      &LesHouchesReader::theMinMultCKKW, 0, 0, 0,
      true, false, Interface::lowerlim);
 
 
   static Switch<LesHouchesReader,unsigned int> interfaceMomentumTreatment
     ("MomentumTreatment",
      "Treatment of the momenta supplied by the interface",
      &LesHouchesReader::theMomentumTreatment, 0, false, false);
   static SwitchOption interfaceMomentumTreatmentAccept
     (interfaceMomentumTreatment,
      "Accept",
      "Just accept the momenta given",
      0);
   static SwitchOption interfaceMomentumTreatmentRescaleEnergy
     (interfaceMomentumTreatment,
      "RescaleEnergy",
      "Rescale the energy supplied so it is consistent with the mass",
      1);
   static SwitchOption interfaceMomentumTreatmentRescaleMass
     (interfaceMomentumTreatment,
      "RescaleMass",
      "Rescale the mass supplied so it is consistent with the"
      " energy and momentum",
      2);
 
 
   static Switch<LesHouchesReader,bool> interfaceWeightWarnings
     ("WeightWarnings",
      "Determines if warnings about possible weight incompatibilities should "
      "be issued when this reader is initialized.",
      &LesHouchesReader::useWeightWarnings, true, true, false);
   static SwitchOption interfaceWeightWarningsWarnAboutWeights
     (interfaceWeightWarnings,
      "WarnAboutWeights",
      "Warn about possible incompatibilities with the weight option in the "
      "Les Houches common block and the requested weight treatment.",
      true);
   static SwitchOption interfaceWeightWarningsDontWarnAboutWeights
     (interfaceWeightWarnings,
      "DontWarnAboutWeights",
      "Do not warn about possible incompatibilities with the weight option "
      "in the Les Houches common block and the requested weight treatment.",
      false);
 
   static Switch<LesHouchesReader,bool> interfaceAllowedTopReOpen
     ("AllowedToReOpen",
      "Can the file be reopened if more events are requested than the file contains?",
      &LesHouchesReader::theReOpenAllowed, true, false, false);
   static SwitchOption interfaceAllowedTopReOpenYes
     (interfaceAllowedTopReOpen,
      "Yes",
      "Allowed to reopen the file",
      true);
   static SwitchOption interfaceAllowedTopReOpenNo
     (interfaceAllowedTopReOpen,
      "No",
      "Not allowed to reopen the file",
      false);
 
   static Switch<LesHouchesReader,bool> interfaceIncludeSpin
     ("IncludeSpin",
      "Use the spin information present in the event file, for tau leptons"
      " only as this is the only case which makes any sense",
      &LesHouchesReader::theIncludeSpin, true, false, false);
   static SwitchOption interfaceIncludeSpinYes
     (interfaceIncludeSpin,
      "Yes",
      "Use the spin information",
      true);
   static SwitchOption interfaceIncludeSpinNo
     (interfaceIncludeSpin,
      "No",
      "Don't use the spin information",
      false);
 
 
 
   interfaceCuts.rank(8);
   interfacePartonExtractor.rank(7);
   interfaceBeamA.rank(5);
   interfaceBeamB.rank(4);
   interfaceEBeamA.rank(3);
   interfaceEBeamB.rank(2);
   interfaceMaxMultCKKW.rank(1.5);
   interfaceMinMultCKKW.rank(1.0);
 
   interfaceBeamA.setHasDefault(false);
   interfaceBeamB.setHasDefault(false);
   interfaceEBeamA.setHasDefault(false);
   interfaceEBeamB.setHasDefault(false);
   interfaceMaxMultCKKW.setHasDefault(false);
   interfaceMinMultCKKW.setHasDefault(false);
 
 }
 
 namespace ThePEG {
 
 ostream & operator<<(ostream & os, const HEPEUP & h) {
   os << "<event>\n"
      << " " << setw(4) << h.NUP
      << " " << setw(6) << h.IDPRUP
      << " " << setw(14) << h.XWGTUP
      << " " << setw(14) << h.SCALUP
      << " " << setw(14) << h.AQEDUP
      << " " << setw(14) << h.AQCDUP << "\n";
 
   for ( int i = 0; i < h.NUP; ++i )
     os << " " << setw(8) << h.IDUP[i]
        << " " << setw(2) << h.ISTUP[i]
        << " " << setw(4) << h.MOTHUP[i].first
        << " " << setw(4) << h.MOTHUP[i].second
        << " " << setw(4) << h.ICOLUP[i].first
        << " " << setw(4) << h.ICOLUP[i].second
        << " " << setw(14) << h.PUP[i][0]
        << " " << setw(14) << h.PUP[i][1]
        << " " << setw(14) << h.PUP[i][2]
        << " " << setw(14) << h.PUP[i][3]
        << " " << setw(14) << h.PUP[i][4]
        << " " << setw(1) << h.VTIMUP[i]
        << " " << setw(1) << h.SPINUP[i] << std::endl;
   os << "</event>" << std::endl;
   return os;
 }
 
 }
 
diff --git a/LesHouches/MadGraphOneCut.cc b/LesHouches/MadGraphOneCut.cc
--- a/LesHouches/MadGraphOneCut.cc
+++ b/LesHouches/MadGraphOneCut.cc
@@ -1,170 +1,170 @@
 // -*- C++ -*-
 //
 // MadGraphOneCut.cc 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.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the MadGraphOneCut class.
 //
 
 #include "MadGraphOneCut.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/PDT/ParticleData.h"
 #include "ThePEG/PDT/EnumParticles.h"
 #include "ThePEG/Utilities/EnumIO.h"
 #include "ThePEG/Cuts/Cuts.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace ThePEG;
 
 IBPtr MadGraphOneCut::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr MadGraphOneCut::fullclone() const {
   return new_ptr(*this);
 }
 
 Energy MadGraphOneCut::minKT(tcPDPtr p) const {
-  if ( cutType != PT || !checkType(p) ) return ZERO;
+  if ( cutType != Cut::PT || !checkType(p) ) return ZERO;
   return theCut*GeV;
 }
 
 double MadGraphOneCut::minEta(tcPDPtr p) const {
-  if ( cutType != ETA || !checkType(p) ) return -Constants::MaxRapidity;
+  if ( cutType != Cut::ETA || !checkType(p) ) return -Constants::MaxRapidity;
   return -theCut;
 }
 
 double MadGraphOneCut::maxEta(tcPDPtr p) const {
-  if ( cutType != ETA || !checkType(p) ) return Constants::MaxRapidity;
+  if ( cutType != Cut::ETA || !checkType(p) ) return Constants::MaxRapidity;
   return theCut;
 }
 
 Energy MadGraphOneCut::minMaxKT(tcPDPtr p) const {
-  if ( cutType != XPT || !checkType(p) ) return ZERO;
+  if ( cutType != Cut::XPT || !checkType(p) ) return ZERO;
   return theCut*GeV;  
 }
 
 bool MadGraphOneCut::passCuts(tcCutsPtr parent,
 			      tcPDPtr ptype, LorentzMomentum p) const {
   if ( !checkType(ptype) ) return true;
-  if ( cutType == PT ) return p.perp() > theCut*GeV;
-  if ( cutType == ETA ) {
+  if ( cutType == Cut::PT ) return p.perp() > theCut*GeV;
+  if ( cutType == Cut::ETA ) {
     double y = p.rapidity() + parent->Y() + parent->currentYHat();
     return abs(p.mt()*sinh(y)) < p.perp()*sinh(theCut);
   }
   return true;
 }
 
 bool MadGraphOneCut::checkType(tcPDPtr p) const {
   switch ( abs(p->id()) ) {
   case ParticleID::d:
   case ParticleID::u:
   case ParticleID::s:
   case ParticleID::c:
   case ParticleID::g:
-    return particleType == JET;
+    return particleType == P::JET;
   case ParticleID::b:
-    return particleType == JET || particleType == BOT;
+    return particleType == P::JET || particleType == P::BOT;
   case ParticleID::gamma:
-    return particleType == PHO;
+    return particleType == P::PHO;
   case ParticleID::eminus:
   case ParticleID::nu_e:
   case ParticleID::muminus:
   case ParticleID::nu_mu:
   case ParticleID::tauminus:
   case ParticleID::nu_tau:
-    return particleType == LEP;
+    return particleType == P::LEP;
   default:
     return false;
   }
 }
 
 void MadGraphOneCut::persistentOutput(PersistentOStream & os) const {
   os << oenum(cutType) << oenum(particleType) << theCut;
 }
 
 void MadGraphOneCut::persistentInput(PersistentIStream & is, int) {
   is >> ienum(cutType) >> ienum(particleType) >> theCut;
 }
 
 ClassDescription<MadGraphOneCut> MadGraphOneCut::initMadGraphOneCut;
 // Definition of the static class description member.
 
 void MadGraphOneCut::Init() {
 
   static ClassDocumentation<MadGraphOneCut> documentation
     ("Objects of the MadGraphOneCut class can be created automatically by "
      "the MadGraphReader class when scanning event files for information "
      "about cuts. It is also possible to create objects by hand and use "
      "it as any other OneCutBase object.");
 
-  static Switch<MadGraphOneCut,CutType> interfaceCutType
+  static Switch<MadGraphOneCut,Cut> interfaceCutType
     ("CutType",
      "The type of cut this object will do.",
-     &MadGraphOneCut::cutType, PT, true, false);
+     &MadGraphOneCut::cutType, Cut::PT, true, false);
   static SwitchOption interfaceCutTypePT
     (interfaceCutType,
      "MinPT",
      "The minimum transverse momentum of a particle.",
-     PT);
+     Cut::PT);
   static SwitchOption interfaceCutTypeMaxEta
     (interfaceCutType,
      "MaxEta",
      "The maximum (absolute value of) pseudo-rapidity of a particle.",
-     ETA);
+     Cut::ETA);
   static SwitchOption interfaceCutTypeMinMaxPT
     (interfaceCutType,
      "MinMaxPT",
      "The minimum transverse momentum of the particle with largest "
      "transverse momentum.",
-     XPT);
+     Cut::XPT);
 
-  static Switch<MadGraphOneCut,PType> interfaceParticleType
+  static Switch<MadGraphOneCut,P> interfaceParticleType
     ("ParticleType",
      "The types of particles this cut is applied to.",
-     &MadGraphOneCut::particleType, JET, true, false);
+     &MadGraphOneCut::particleType, P::JET, true, false);
   static SwitchOption interfaceParticleTypeJets
     (interfaceParticleType,
      "Jets",
      "The cut applies only to coloured particles (jets).",
-     JET);
+     P::JET);
   static SwitchOption interfaceParticleTypeLeptons
     (interfaceParticleType,
      "Leptons",
      "The cut applies only to leptons.",
-     LEP);
+     P::LEP);
   static SwitchOption interfaceParticleTypePhotons
     (interfaceParticleType,
      "Photons",
      "The cut applies only to photons.",
-     PHO);
+     P::PHO);
   static SwitchOption interfaceParticleTypeBottom
     (interfaceParticleType,
      "Bottom",
      "The cut applies only to bottom quarks.",
-     BOT);
+     P::BOT);
 
   static Parameter<MadGraphOneCut,double> interfaceCut
     ("Cut",
      "The value of the cut to be applied (in units of GeV in case of a "
      "transverse momentum).",
      &MadGraphOneCut::theCut, 0.0, 0.0, 0,
      true, false, Interface::lowerlim);
 
   interfaceCut.rank(10);
   interfaceCutType.rank(9);
   interfaceParticleType.rank(8);
   interfaceCut.setHasDefault(false);
   interfaceCutType.setHasDefault(false);
   interfaceParticleType.setHasDefault(false);
 
 }
 
diff --git a/LesHouches/MadGraphOneCut.h b/LesHouches/MadGraphOneCut.h
--- a/LesHouches/MadGraphOneCut.h
+++ b/LesHouches/MadGraphOneCut.h
@@ -1,227 +1,227 @@
 // -*- C++ -*-
 //
 // MadGraphOneCut.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_MadGraphOneCut_H
 #define THEPEG_MadGraphOneCut_H
 //
 // This is the declaration of the MadGraphOneCut class.
 //
 
 #include "ThePEG/Cuts/OneCutBase.h"
 
 namespace ThePEG {
 
 /**
  * Objects of the MadGraphOneCut class can be created automatically by
  * the MadGraphReader class when scanning event files for information
  * about cuts. It is also possible to create objects by hand and use
  * it as any other OneCutBase object.
  *
  * @see \ref MadGraphOneCutInterfaces "The interfaces"
  * defined for MadGraphOneCut.
  */
 class MadGraphOneCut: public OneCutBase {
 
 public:
 
   /**
    * Enumerate the different kinds of cuts made by MadGraph.
    */
-  enum CutType {
+  enum class Cut {
     PT,  /**< The minimum transverse momentum of a particle. */
     ETA, /**< The maximum (absolute value of) pseudo-rapidity of a particle. */
     XPT  /**< The minimum transverse momentum of the particle with
 	      largest transverse momentum. */
   };
 
   /**
    * Enumerate the types of particles the cut is made on.
    */
-  enum PType {
+  enum class P {
     JET, /**< The cut applies only to coloured particles. */
     LEP, /**< The cut applies only to leptons. */
     PHO, /**< The cut applies only to photons. */
     BOT  /**< The cut applies only to bottom quarks. */
   };
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * The default constructor.
    */
-  MadGraphOneCut() : cutType(PT), particleType(JET), theCut(0.0) {}
+  MadGraphOneCut() : cutType(Cut::PT), particleType(P::JET), theCut(0.0) {}
 
   /**
    * The constructor used by the MadGraphReader.
    * @param t is the type of the cut.
    * @param p is the type of particles the cut is applied to.
    * @param c is the value of the cut (in units of GeV where applicable).
    */
-  MadGraphOneCut(CutType t, PType p, double c)
+  MadGraphOneCut(Cut t, P p, double c)
     : cutType(t), particleType(p), theCut(c) {}
   //@}
 
 public:
 
   /** @name Virtual functions mandated by the base class. */
   //@{
   /**
    * Return the minimum allowed value of the transverse momentum of an
    * outgoing parton.
    */
   virtual Energy minKT(tcPDPtr p) const;
 
   /**
    * Return the minimum allowed pseudo-rapidity of an outgoing parton
    * of the given type. The pseudo-rapidity is measured in the lab
    * system.
    */
   virtual double minEta(tcPDPtr p) const;
 
   /**
    * Return the maximum allowed pseudo-rapidity of an outgoing parton
    * of the given type. The pseudo-rapidity is measured in the lab
    * system.
    */
   virtual double maxEta(tcPDPtr p) const;
 
   /**
    * Return the minimum allowed value of the transverse momentum of
    * the outgoing parton with the lagrest transverse momentum. This
    * version simply returns minKt().
    */
   virtual Energy minMaxKT(tcPDPtr p) const;
 
   /**
    * Return true if a particle with type \a ptype and momentum \a p
    * passes the cuts. The \a parent contains information about the
    * kinematics of the hard sub-process.
    */
   virtual bool passCuts(tcCutsPtr parent,
 			tcPDPtr ptype, LorentzMomentum p) const;
   //@}
 
 protected:
 
   /**
    * Returns true if cut should be applied to a particle of type \a p.
    */
   bool checkType(tcPDPtr p) 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;
   //@}
 
 private:
 
   /**
    * The type of this cut.
    */
-  CutType cutType;
+  Cut cutType;
 
   /**
    * The type of particles this cut applies to.
    */
-  PType particleType;
+  P particleType;
 
   /**
    * The value of the cut to be applied.
    */
   double theCut;
 
 private:
 
   /**
    * The static object used to initialize the description of this class.
    * Indicates that this is a concrete class with persistent data.
    */
   static ClassDescription<MadGraphOneCut> initMadGraphOneCut;
 
   /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
   MadGraphOneCut & operator=(const MadGraphOneCut &);
 
 };
 
 }
 
 #include "ThePEG/Utilities/ClassTraits.h"
 
 namespace ThePEG {
 
 /** @cond TRAITSPECIALIZATIONS */
 
 /** This template specialization informs ThePEG about the
  *  base classes of MadGraphOneCut. */
 template <>
 struct BaseClassTrait<MadGraphOneCut,1> {
   /** Typedef of the first base class of MadGraphOneCut. */
   typedef OneCutBase NthBase;
 };
 
 /** This template specialization informs ThePEG about the name of
  *  the MadGraphOneCut class and the shared object where it is defined. */
 template <>
 struct ClassTraits<MadGraphOneCut>
   : public ClassTraitsBase<MadGraphOneCut> {
   /** Return a platform-independent class name */
   static string className() { return "ThePEG::MadGraphOneCut"; }
   /** Return the name(s) of the shared library (or libraries) be loaded to get
    *  access to the MadGraphOneCut class and any other class on which it depends
    *  (except the base class). */
   static string library() { return "MadGraphReader.so"; }
 };
 
 /** @endcond */
 
 }
 
 #endif /* THEPEG_MadGraphOneCut_H */
diff --git a/LesHouches/MadGraphReader.cc b/LesHouches/MadGraphReader.cc
--- a/LesHouches/MadGraphReader.cc
+++ b/LesHouches/MadGraphReader.cc
@@ -1,618 +1,618 @@
 // -*- C++ -*-
 //
 // MadGraphReader.cc 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.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the MadGraphReader class.
 //
 
 #include "MadGraphReader.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Interface/Command.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/PDT/EnumParticles.h"
 #include "ThePEG/PDT/ParticleData.h"
 #include "ThePEG/PDF/PDFBase.h"
 #include "ThePEG/LesHouches/MadGraphOneCut.h"
 #include "ThePEG/LesHouches/MadGraphTwoCut.h"
 #include "ThePEG/Cuts/Cuts.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Utilities/Throw.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace ThePEG;
 using std::fgetc;
 using std::fgets;
 
 IBPtr MadGraphReader::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr MadGraphReader::fullclone() const {
   return new_ptr(*this);
 }
 
 void MadGraphReader::open() {
   LesHouchesFileReader::open();
 
   heprup.IDWTUP = 1;
 
   static const int ntags = 33;
   static const char * cuttags[] = {"ptj", "ptb", "pta", "ptl",
 				   "etaj", "etab", "etaa", "etal",
 				   "drjj", "drbb", "draa", "drll",
 				   "drbj", "draj", "drjl", "drab",
 				   "drbl", "dral", "mmjj", "mmbb",
 				   "mmaa", "mmll", "mmbj", "mmaj",
 				   "mmjl", "mmab", "mmbl", "mmal", 
 				   "xptj", "xptb", "xpta", "xptl", "xetamin"};
 
   ieve = neve = 0;
 
   // If we are reading a LHF formatted file things are rather easy.
 
   if ( LHFVersion.size() ) {
     // extract the number of events
     neve = numberOfEvents(outsideBlock);
     string madheader = outsideBlock;
     if ( neve == 0 ) {
       neve = numberOfEvents(headerBlock);
       madheader = headerBlock;
     }
     if ( neve == 0 )
       Throw<LesHouchesFileError>()
       << "The MadGraphReader '" << name() << "' expected the LHE file '"
       << filename() << "' to include MadGraph-specific header information, "
       << "but did not find any. The events may not be properly sampled."
       << Exception::warning;
 
     if ( neve != 0 ) NEvents(neve);
 
     // MadEvent has gives wrong values for XMAXUP and XWGTUP, they
     // need to be multiplied by the number of events to be LHEF
     // compliant.
     weightScale = neve*picobarn;
 
     // Scan information about cuts.
     for ( int itag = 0; itag < ntags; ++itag ) {
      string::size_type pos = madheader.find(string("= ") + cuttags[itag]);
       if ( pos != string::npos ) {
 	string::size_type beg = max(madheader.rfind("#", pos) + 1,
 				    madheader.rfind("\n", pos) + 1);
 	string value = madheader.substr(beg, pos - beg);
 	for ( string::size_type i = 0; i < value.length(); ++i )
 	  if ( value[i] == 'd' || value[i] == 'D' ) value[i] = 'e';
 	cuts[cuttags[itag]] = std::strtod(value.c_str(), NULL);
       }
     }
 
     return;
 
   }
 
   double xsec = -1.0;
   double maxw = -1.0;
   double ebeam1 = -1.0;
   double ebeam2 = -1.0;
   int lpp1 = 0;
   int lpp2 = 0;
   string pdftag;
   // First scan banner to extract some information
   // (LesHoushesFileReader::open has already read in the first line).
   cfile.resetline();
   do {
     if ( !cfile ) break;
     if ( cfile.getc() != '#' ) {
       int test;
       cfile >> test;
       if ( cfile ) {
 	cfile.resetline();
 	break;
       }
     }
     if ( cfile.find("#  Number of Events") ) {
       cfile.skip(':');
       cfile >> neve;
     } else if ( cfile.find("Integrated weight") ) {
       cfile.skip(':');
       cfile >> xsec;
     } else if ( cfile.find("Max wgt") ) {
       cfile.skip(':');
       cfile >> maxw;
     } else if ( cfile.find("ebeam(1)") || cfile.find("ebeam1") ) {
       cfile >> ebeam1;
     } else if ( cfile.find("ebeam(2)") || cfile.find("ebeam2") ) {
       cfile >> ebeam2;
     } else if ( cfile.find("lpp(1)") || cfile.find("lpp1") ) {
       cfile >> lpp1;
     } else if ( cfile.find("lpp(2)") || cfile.find("lpp2") ) {
       cfile >> lpp2;
     } else if ( cfile.find("PDF set") ) {
       cfile.skip('\'');
       cfile >> pdftag;
       pdftag = pdftag.substr(0, 7);
     } else if ( cfile.find("Number of Events Written ") ) {
       cfile.skip(':');
       cfile >> neve;
       maxw = xsec/double(neve);
     } else {
       for ( int itag = 0; itag < ntags; ++itag ) {
 	if ( cfile.find(string("= ") + cuttags[itag]) ) {
 	  
 	  cfile >> cuts[cuttags[itag]];
 	  if ( cfile.getc() == 'd' ) {
 	    long x = 0;
 	    cfile >> x;
 	    cuts[cuttags[itag]] *= pow(10.0, double(x));
 	  }
 	  break;
 	}
       }
     }
   } while ( cfile.readline() );
 
   // Return here if no comment block was found.
   if ( neve <= 0 ) return;
 
   // Convert the extracted information to LesHouches format.
   heprup.NPRUP = 1;
   heprup.LPRUP.push_back(0);
   heprup.XSECUP.push_back(xsec);
   heprup.XERRUP.push_back(0.0);
   heprup.XMAXUP.push_back(maxw);
   NEvents(neve);
   // MadEvent has gives wrong values for XMAXUP and XWGTUP, they
   // need to be multiplied by the number of events to be LHEF
   // compliant.
   weightScale = neve*picobarn;
 
   if ( !heprup.IDBMUP.first ) {
     if ( lpp1 == 1 ) heprup.IDBMUP.first = ParticleID::pplus;
     else if ( lpp1 == -1 ) heprup.IDBMUP.first = ParticleID::pbarminus;
   }
   if ( !heprup.IDBMUP.second ) {
     if ( lpp2 == 1 ) heprup.IDBMUP.second = ParticleID::pplus;
     else if ( lpp2 == -1 ) heprup.IDBMUP.second = ParticleID::pbarminus;
   }
 
   if ( heprup.EBMUP.first <= 0.0 ) heprup.EBMUP.first = ebeam1;
   if ( heprup.EBMUP.second <= 0.0 ) heprup.EBMUP.second = ebeam2;
 
   if ( !cfile )
     throw LesHouchesFileError()
       << "An error occurred while '" << name() << "' was reading the file '"
       << filename() << "'." << Exception::runerror;
 
   if ( heprup.PDFSUP.first != 0 || heprup.PDFSUP.first != 0 ) return;
   // If we have an old MadGraph we have to try to figure out which PDF
   // codes to use.
   heprup.PDFGUP.first = heprup.PDFGUP.second = 0;
   if ( pdftag == "mrs02nl" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 20200;
   else if ( pdftag == "mrs02nn" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 20270;
   else if ( pdftag == "cteq6_m" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 10050;
   else if ( pdftag == "cteq6_l" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 10041;
   else if ( pdftag == "cteq6l1" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 10042;
   else if ( pdftag == "cteq5_m" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 19050;
   else if ( pdftag == "cteq5_d" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 19060;
   else if ( pdftag == "cteq5_l" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 19070;
   else if ( pdftag == "cteq4_m" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 19150;
   else if ( pdftag == "cteq4_d" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 19160;
   else if ( pdftag == "cteq4_l" )
     heprup.PDFSUP.first = heprup.PDFSUP.second = 19170;
 }
 
 
 long MadGraphReader::scan() {
   bool fixscale = !NEvents();
   long neve = LesHouchesFileReader::scan();
   if ( fixscale ) {
     // MadEvent has gives wrong values for XMAXUP and XWGTUP, they
     // need to be multiplied by the number of events to be LHEF
     // compliant.
     weightScale = neve*picobarn;
     if ( heprup.NPRUP > 1 ) weightScale /= heprup.NPRUP;
   }
   return neve;
 }
 
 bool MadGraphReader::doReadEvent() {
 
   if ( LesHouchesFileReader::doReadEvent() )  return true;
 
   if ( !cfile ) return false;
 
   hepeup.NUP = 0;
   ieve = 0;
   long evno = 0;
   hepeup.XWGTUP = 0.0;
   double scale = 0.0;
   double aEM = 0.0;
   double aS = 0.0;
   bool oldformat = false;
 
   cfile >> hepeup.NUP >> evno >> hepeup.XWGTUP >> scale >> aEM >> aS;
   if ( !cfile ) {
     hepeup.IDPRUP = evno;
     hepeup.SCALUP = fixedScale/GeV;
     hepeup.AQEDUP = fixedAEM;
     hepeup.AQCDUP = fixedAS;
     ++ieve;
     oldformat = true;
   } else {
     hepeup.IDPRUP = 0;
     ieve = evno;
     hepeup.SCALUP = scale;
     hepeup.AQEDUP = aEM;
     hepeup.AQCDUP = aS;
   }
 
   hepeup.IDUP.resize(hepeup.NUP);
   if ( !cfile.readline() ) return false;
   for ( int i = 0; i < hepeup.NUP; ++i ) cfile >> hepeup.IDUP[i];
 
   hepeup.MOTHUP.resize(hepeup.NUP);
   if ( !cfile.readline() ) return false;
   for ( int i = 0; i < hepeup.NUP; ++i ) cfile >> hepeup.MOTHUP[i].first;
   if ( !cfile.readline() ) return false;
   for ( int i = 0; i < hepeup.NUP; ++i ) cfile >> hepeup.MOTHUP[i].second;
 
   hepeup.ICOLUP.resize(hepeup.NUP);
   if ( !cfile.readline() ) return false;
   for ( int i = 0; i < hepeup.NUP; ++i ) cfile >> hepeup.ICOLUP[i].first;
   if ( !cfile.readline() ) return false;
   for ( int i = 0; i < hepeup.NUP; ++i ) cfile >> hepeup.ICOLUP[i].second;
 
   // Try to figure out if the colour lines are reversed
   bool colrev = false;
   for ( int i = 0; i < hepeup.NUP; ++i )
     if ( abs(hepeup.IDUP[i]) < 10 && hepeup.IDUP[i] < 0 &&
 	 !hepeup.ICOLUP[i].second ) colrev = true;
   if ( colrev ) for ( int i = 0; i < hepeup.NUP; ++i )
     swap(hepeup.ICOLUP[i].first, hepeup.ICOLUP[i].second);
 
   if ( oldformat ) {
     hepeup.ISTUP.assign(hepeup.NUP, 1);
     hepeup.ISTUP[0] = hepeup.ISTUP[1] = -1;
     hepeup.SPINUP.assign(hepeup.NUP, 9);
   } else {
     hepeup.ISTUP.resize(hepeup.NUP);
     if ( !cfile.readline() ) return false;
     for ( int i = 0; i < hepeup.NUP; ++i ) cfile >> hepeup.ISTUP[i];
     hepeup.SPINUP.resize(hepeup.NUP, 9);
     if ( !cfile.readline() ) return false;
     for ( int i = 0; i < hepeup.NUP; ++i ) cfile >> hepeup.SPINUP[i];
   }
 
-  hepeup.PUP.resize(hepeup.NUP, vector<double>(5));
+  hepeup.PUP.resize(hepeup.NUP);
   for ( int i = 0; i < hepeup.NUP; ++i ) {
     if ( !cfile.readline() ) return false;
     int dummy = 0;
     cfile  >> dummy >> hepeup.PUP[i][3]
 	   >> hepeup.PUP[i][0] >> hepeup.PUP[i][1] >> hepeup.PUP[i][2];
     hepeup.PUP[i][4] =
       sqrt(max(sqr(hepeup.PUP[i][3]) - sqr(hepeup.PUP[i][0]) -
 	       sqr(hepeup.PUP[i][1]) - sqr(hepeup.PUP[i][2]), 0.0));
   }
 
   if ( !cfile ) return false;
 
   // Set info not obtained from MadGraph.
   hepeup.VTIMUP = vector<double>(hepeup.NUP, -1.0);
 
   // Deduce positions of incoming beams and corresponding partons.
   pair<int,int> beampos(-1, -1);
   for ( int i = 0; i < hepeup.NUP; ++i ) {
     if ( hepeup.ISTUP[i] != -9 ) continue;
     if ( beampos.first < 0 ) beampos.first = i;
     else if ( beampos.second < 0 ) beampos.second = i;
   }
   pair<int,int> partpos(-1, -1);
   for ( int i = hepeup.NUP - 1; i >= 0; --i ) {
     if ( hepeup.ISTUP[i] != -1 ) continue;
     if ( hepeup.MOTHUP[i].first > 0 &&
 	 hepeup.MOTHUP[i].first == beampos.first )
       partpos.first = i;
     else if ( hepeup.MOTHUP[i].first > 0 &&
 	      hepeup.MOTHUP[i].first == beampos.second )
       partpos.second = i;
     else if ( partpos.second < 0 ) partpos.second = i;
     else if ( partpos.first < 0 ) partpos.first = i;
   }
 
   // We set these to -1 to let the base class do the work.
   hepeup.XPDWUP.first = -1.0;
   hepeup.XPDWUP.second = -1.0;
 
   cfile.readline();
 
   // Return true even if last read failed.
   return true;
 
 }
 
 void MadGraphReader::persistentOutput(PersistentOStream & os) const {
   os << ounit(fixedScale, GeV) << fixedAEM << fixedAS << cuts
      << doInitCuts;
 }
 
 void MadGraphReader::persistentInput(PersistentIStream & is, int) {
   is >> iunit(fixedScale, GeV) >> fixedAEM >> fixedAS >> cuts
      >> doInitCuts;
 }
 
 bool MadGraphReader::preInitialize() const {
   if ( LesHouchesFileReader::preInitialize() ) return true;
   if ( doInitCuts && !theCuts ) return true;
   return false;
 }
 
 void MadGraphReader::doinit() {
   LesHouchesFileReader::doinit();
   if ( doInitCuts && !theCuts ) {
     theCuts = initCuts();
     if ( !theCuts ) Throw<Exception>()
       << "MadGraphReader '" << name()
       << "' could not create cut objects in pre-initialization."
       << Exception::warning;
   }
 }
 
 void MadGraphReader::initPDFs() {
   LesHouchesFileReader::initPDFs();
 }
 
 CutsPtr MadGraphReader::initCuts() {
   CutsPtr newCuts;
   open();
   close();
   if ( cuts.empty() ) return CutsPtr();
   vector<OneCutPtr> ones;
   vector<TwoCutPtr> twos;
   vector<string> onames;
   vector<string> tnames;
   for ( map<string,double>::iterator i = cuts.begin(); i != cuts.end(); ++i ) {
     if ( i->second <= 0.0 ) continue;
-    MadGraphOneCut::CutType t = MadGraphOneCut::PT;
+    MadGraphOneCut::Cut t = MadGraphOneCut::Cut::PT;
     char p = 0;
     if ( i->first.substr(0, 2) == "pt" ) {
-      t = MadGraphOneCut::PT;
+      t = MadGraphOneCut::Cut::PT;
       p = i->first[2];
     }
     else if  ( i->first.substr(0, 3) == "eta" ) {
-      t = MadGraphOneCut::ETA;
+      t = MadGraphOneCut::Cut::ETA;
       p = i->first[3];
     }
     else if  ( i->first.substr(0, 3) == "xpt" ) {
-      t = MadGraphOneCut::XPT;
+      t = MadGraphOneCut::Cut::XPT;
       p = i->first[3];
     }
     if ( p ) {
-      MadGraphOneCut::PType pt = MadGraphOneCut::JET;
+      MadGraphOneCut::P pt = MadGraphOneCut::P::JET;
       switch ( p ) {
-      case 'j':	pt = MadGraphOneCut::JET; break;
-      case 'b':	pt = MadGraphOneCut::BOT; break;
-      case 'a':	pt = MadGraphOneCut::PHO; break;
-      case 'l':	pt = MadGraphOneCut::LEP; break;
+      case 'j':	pt = MadGraphOneCut::P::JET; break;
+      case 'b':	pt = MadGraphOneCut::P::BOT; break;
+      case 'a':	pt = MadGraphOneCut::P::PHO; break;
+      case 'l':	pt = MadGraphOneCut::P::LEP; break;
       }
       ones.push_back(new_ptr(MadGraphOneCut(t, pt, i->second)));
       onames.push_back(i->first);
       continue;
     }
     if ( i->first.substr(0, 2) == "dr" || i->first.substr(0, 2) == "mm" ) {
-      MadGraphTwoCut::CutType tt = MadGraphTwoCut::DELTAR;
-      if ( i->first.substr(0, 2) == "mm" ) tt = MadGraphTwoCut::INVMASS;
-      MadGraphTwoCut::PPType pp = MadGraphTwoCut::JETJET;
-      if ( i->first.substr(2, 2) == "jj" ) pp = MadGraphTwoCut::JETJET;
-      else if ( i->first.substr(2, 2) == "bb" )	pp = MadGraphTwoCut::BOTBOT;
-      else if ( i->first.substr(2, 2) == "aa" )	pp = MadGraphTwoCut::PHOPHO;
-      else if ( i->first.substr(2, 2) == "ll" )	pp = MadGraphTwoCut::LEPLEP;
-      else if ( i->first.substr(2, 2) == "bj" )	pp = MadGraphTwoCut::BOTJET;
-      else if ( i->first.substr(2, 2) == "aj" )	pp = MadGraphTwoCut::PHOJET;
-      else if ( i->first.substr(2, 2) == "jl" )	pp = MadGraphTwoCut::JETLEP;
-      else if ( i->first.substr(2, 2) == "ab" )	pp = MadGraphTwoCut::PHOBOT;
-      else if ( i->first.substr(2, 2) == "bl" )	pp = MadGraphTwoCut::BOTLEP;
-      else if ( i->first.substr(2, 2) == "al" )	pp = MadGraphTwoCut::PHOLEP;
+      MadGraphTwoCut::Cut tt = MadGraphTwoCut::Cut::DELTAR;
+      if ( i->first.substr(0, 2) == "mm" ) tt = MadGraphTwoCut::Cut::INVMASS;
+      MadGraphTwoCut::PP pp = MadGraphTwoCut::PP::JETJET;
+      if ( i->first.substr(2, 2) == "jj" ) pp = MadGraphTwoCut::PP::JETJET;
+      else if ( i->first.substr(2, 2) == "bb" )	pp = MadGraphTwoCut::PP::BOTBOT;
+      else if ( i->first.substr(2, 2) == "aa" )	pp = MadGraphTwoCut::PP::PHOPHO;
+      else if ( i->first.substr(2, 2) == "ll" )	pp = MadGraphTwoCut::PP::LEPLEP;
+      else if ( i->first.substr(2, 2) == "bj" )	pp = MadGraphTwoCut::PP::BOTJET;
+      else if ( i->first.substr(2, 2) == "aj" )	pp = MadGraphTwoCut::PP::PHOJET;
+      else if ( i->first.substr(2, 2) == "jl" )	pp = MadGraphTwoCut::PP::JETLEP;
+      else if ( i->first.substr(2, 2) == "ab" )	pp = MadGraphTwoCut::PP::PHOBOT;
+      else if ( i->first.substr(2, 2) == "bl" )	pp = MadGraphTwoCut::PP::BOTLEP;
+      else if ( i->first.substr(2, 2) == "al" )	pp = MadGraphTwoCut::PP::PHOLEP;
       twos.push_back(new_ptr(MadGraphTwoCut(tt, pp, i->second)));
       tnames.push_back(i->first);
     }
   }
   if ( ones.empty() && twos.empty() ) return CutsPtr();
 
   newCuts = new_ptr(Cuts());
   generator()->preinitRegister(newCuts, fullName() + "/ExtractedCuts");
 
   for ( int i = 0, N = ones.size(); i < N; ++i ) {
     generator()->preinitRegister(ones[i], fullName() + "/" + onames[i]);
     generator()->preinitInterface
       (newCuts, "OneCuts", 0, "insert",  ones[i]->fullName());
     //    newCuts->add(tOneCutPtr(ones[i]));
   }
 
   for ( int i = 0, N = twos.size(); i < N; ++i ) {
     reporeg(twos[i], tnames[i]);
     generator()->preinitInterface
       (newCuts, "TwoCuts", 0, "insert",  twos[i]->fullName());
     //    newCuts->add(tTwoCutPtr(twos[i]));
   }
 
   return newCuts;
 
 
 }
 
 string MadGraphReader::scanCuts(string) {
   if ( theCuts )
     return "A Cuts object has already been assigned to this reader.";
   open();
   close();
   if ( cuts.empty() ) return "No information about cuts were found. "
 			"Maybe the file was from an old version of MadGraph";
   vector<OneCutPtr> ones;
   vector<TwoCutPtr> twos;
   vector<string> onames;
   vector<string> tnames;
   for ( map<string,double>::iterator i = cuts.begin(); i != cuts.end(); ++i ) {
     if ( i->second <= 0.0 ) continue;
-    MadGraphOneCut::CutType t = MadGraphOneCut::PT;
+    MadGraphOneCut::Cut t = MadGraphOneCut::Cut::PT;
     char p = 0;
     if ( i->first.substr(0, 2) == "pt" ) {
-      t = MadGraphOneCut::PT;
+      t = MadGraphOneCut::Cut::PT;
       p = i->first[2];
     }
     else if  ( i->first.substr(0, 3) == "eta" ) {
-      t = MadGraphOneCut::ETA;
+      t = MadGraphOneCut::Cut::ETA;
       p = i->first[3];
     }
     else if  ( i->first.substr(0, 3) == "xpt" ) {
-      t = MadGraphOneCut::XPT;
+      t = MadGraphOneCut::Cut::XPT;
       p = i->first[3];
     }
     if ( p ) {
-      MadGraphOneCut::PType pt = MadGraphOneCut::JET;
+      MadGraphOneCut::P pt = MadGraphOneCut::P::JET;
       switch ( p ) {
-      case 'j':	pt = MadGraphOneCut::JET; break;
-      case 'b':	pt = MadGraphOneCut::BOT; break;
-      case 'a':	pt = MadGraphOneCut::PHO; break;
-      case 'l':	pt = MadGraphOneCut::LEP; break;
+      case 'j':	pt = MadGraphOneCut::P::JET; break;
+      case 'b':	pt = MadGraphOneCut::P::BOT; break;
+      case 'a':	pt = MadGraphOneCut::P::PHO; break;
+      case 'l':	pt = MadGraphOneCut::P::LEP; break;
       }
       ones.push_back(new_ptr(MadGraphOneCut(t, pt, i->second)));
       onames.push_back(i->first);
       continue;
     }
     if ( i->first.substr(0, 2) == "dr" || i->first.substr(0, 2) == "mm" ) {
-      MadGraphTwoCut::CutType tt = MadGraphTwoCut::DELTAR;
-      if ( i->first.substr(0, 2) == "mm" ) tt = MadGraphTwoCut::INVMASS;
-      MadGraphTwoCut::PPType pp = MadGraphTwoCut::JETJET;
-      if ( i->first.substr(2, 2) == "jj" ) pp = MadGraphTwoCut::JETJET;
-      else if ( i->first.substr(2, 2) == "bb" )	pp = MadGraphTwoCut::BOTBOT;
-      else if ( i->first.substr(2, 2) == "aa" )	pp = MadGraphTwoCut::PHOPHO;
-      else if ( i->first.substr(2, 2) == "ll" )	pp = MadGraphTwoCut::LEPLEP;
-      else if ( i->first.substr(2, 2) == "bj" )	pp = MadGraphTwoCut::BOTJET;
-      else if ( i->first.substr(2, 2) == "aj" )	pp = MadGraphTwoCut::PHOJET;
-      else if ( i->first.substr(2, 2) == "jl" )	pp = MadGraphTwoCut::JETLEP;
-      else if ( i->first.substr(2, 2) == "ab" )	pp = MadGraphTwoCut::PHOBOT;
-      else if ( i->first.substr(2, 2) == "bl" )	pp = MadGraphTwoCut::BOTLEP;
-      else if ( i->first.substr(2, 2) == "al" )	pp = MadGraphTwoCut::PHOLEP;
+      MadGraphTwoCut::Cut tt = MadGraphTwoCut::Cut::DELTAR;
+      if ( i->first.substr(0, 2) == "mm" ) tt = MadGraphTwoCut::Cut::INVMASS;
+      MadGraphTwoCut::PP pp = MadGraphTwoCut::PP::JETJET;
+      if ( i->first.substr(2, 2) == "jj" ) pp = MadGraphTwoCut::PP::JETJET;
+      else if ( i->first.substr(2, 2) == "bb" )	pp = MadGraphTwoCut::PP::BOTBOT;
+      else if ( i->first.substr(2, 2) == "aa" )	pp = MadGraphTwoCut::PP::PHOPHO;
+      else if ( i->first.substr(2, 2) == "ll" )	pp = MadGraphTwoCut::PP::LEPLEP;
+      else if ( i->first.substr(2, 2) == "bj" )	pp = MadGraphTwoCut::PP::BOTJET;
+      else if ( i->first.substr(2, 2) == "aj" )	pp = MadGraphTwoCut::PP::PHOJET;
+      else if ( i->first.substr(2, 2) == "jl" )	pp = MadGraphTwoCut::PP::JETLEP;
+      else if ( i->first.substr(2, 2) == "ab" )	pp = MadGraphTwoCut::PP::PHOBOT;
+      else if ( i->first.substr(2, 2) == "bl" )	pp = MadGraphTwoCut::PP::BOTLEP;
+      else if ( i->first.substr(2, 2) == "al" )	pp = MadGraphTwoCut::PP::PHOLEP;
       twos.push_back(new_ptr(MadGraphTwoCut(tt, pp, i->second)));
       tnames.push_back(i->first);
     }
   }
   if ( ones.empty() && twos.empty() ) return "No non-zero cuts found.";
 
   theCuts = new_ptr(Cuts());
   reporeg(theCuts, "ExtractedCuts");
 
   for ( int i = 0, N = ones.size(); i < N; ++i ) {
     reporeg(ones[i], onames[i]);
     theCuts->add(tOneCutPtr(ones[i]));
   }
 
   for ( int i = 0, N = twos.size(); i < N; ++i ) {
     reporeg(twos[i], tnames[i]);
     theCuts->add(tTwoCutPtr(twos[i]));
   }
 
   return "";
 }
 
 ClassDescription<MadGraphReader> MadGraphReader::initMadGraphReader;
 // Definition of the static class description member.
 
 void MadGraphReader::Init() {
 
   static ClassDocumentation<MadGraphReader> documentation
     ("ThePEG::MadGraphReader is used together with the LesHouchesEventHandler "
      "to read event files generated with the MadGraph/MadEvent program.",
      "Events were read from event files generated "
      "with the MadGraph/MadEvent\\cite{ThePEG::MadGraph} program.",
      "\\bibitem{ThePEG::MadGraph} F. Maltoni and T. Stelzer, "
      "hep-ph/0208156;\\\\"
      "T. Stelzer and W.F. Long, \\textit{Comput.~Phys.~Commun.} "
      "\\textbf{81} (1994) 357-371.");
 
   static Parameter<MadGraphReader,Energy> interfaceFixedScale
     ("FixedScale",
      "Old MadGraph files do not necessarily contain information about "
      "the factorization (or renormalization) scale. In this case this "
      "is used instead.",
      &MadGraphReader::fixedScale, GeV, 91.188*GeV, ZERO, 1000.0*GeV,
      true, false, true);
   interfaceFixedScale.setHasDefault(false);
 
   static Parameter<MadGraphReader,double> interfaceFixedAlphaEM
     ("FixedAlphaEM",
      "Old MadGraph files do not necessarily contain information about "
      "the value of \\f$\\alpha_{EM}\\f$. In this case this is used instead.",
      &MadGraphReader::fixedAEM, 0.007546772, 0.0, 1.0,
      true, false, true);
   interfaceFixedAlphaEM.setHasDefault(false);
 
   static Parameter<MadGraphReader,double> interfaceFixedAlphaS
     ("FixedAlphaS",
      "Old MadGraph files do not necessarily contain information about "
      "the value of \\f$\\alpha_S\\f$. In this case this is used instead.",
      &MadGraphReader::fixedAS, 0.12, 0.0, 1.0,
      true, false, true);
   interfaceFixedAlphaS.setHasDefault(false);
 
   static Command<MadGraphReader> interfaceScanCuts
     ("ScanCuts",
      "If no <interface>LesHouchesReader::Cuts</interface> has been assigned, "
      "the event file is scanned for information about generation cuts. If cuts "
      "are found, the corresponding objects will be created in a sub-directory "
      "with the same name as this object and assigned as the "
      "<interface>LesHouchesReader::Cuts</interface> of this reader.",
      &MadGraphReader::scanCuts, true);
 
   static Switch<MadGraphReader,bool> interfaceInitCuts
     ("InitCuts",
      "If no cuts were specified for this reader, try to extract cut "
      "information from the MadGraph file and assign the relevant cut "
      "objects when the reader is initialized.",
      &MadGraphReader::doInitCuts, false, true, false);
   static SwitchOption interfaceInitCutsYes
     (interfaceInitCuts,
      "Yes",
      "Extract cuts during initialization.",
      true);
   static SwitchOption interfaceInitCutsNo
     (interfaceInitCuts,
      "No",
      "Do not extract cuts during initialization.",
      false);
 
 
   interfaceScanCuts.rank(10.5);
   interfaceInitCuts.rank(10.6);
 
 }
 
 long MadGraphReader::numberOfEvents(string block) {
   long output(0);
   // Check for number of events in the file.
   string::size_type pos = block.find("##  Number of Events       :");
   if ( pos == string::npos )
     pos = block.find("#  Number of Events        :");
   if ( pos != string::npos ) {
     pos += 28;
     output = std::strtol(block.c_str() + pos, NULL, 0);
   }
   return output;
 }
diff --git a/LesHouches/MadGraphTwoCut.cc b/LesHouches/MadGraphTwoCut.cc
--- a/LesHouches/MadGraphTwoCut.cc
+++ b/LesHouches/MadGraphTwoCut.cc
@@ -1,242 +1,242 @@
 // -*- C++ -*-
 //
 // MadGraphTwoCut.cc 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.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the MadGraphTwoCut class.
 //
 
 #include "MadGraphTwoCut.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/PDT/ParticleData.h"
 #include "ThePEG/PDT/EnumParticles.h"
 #include "ThePEG/Utilities/EnumIO.h"
 #include "ThePEG/Cuts/Cuts.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace ThePEG;
 
 IBPtr MadGraphTwoCut::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr MadGraphTwoCut::fullclone() const {
   return new_ptr(*this);
 }
 
 Energy2 MadGraphTwoCut::minSij(tcPDPtr pi, tcPDPtr pj) const {
-  if ( !checkType(pi, pj) || cutType != INVMASS ) return ZERO;
+  if ( !checkType(pi, pj) || cutType != Cut::INVMASS ) return ZERO;
   return sqr(theCut*GeV);
 }
 
 Energy2 MadGraphTwoCut::minTij(tcPDPtr, tcPDPtr) const {
   return ZERO;
 }
 
 double MadGraphTwoCut::minDeltaR(tcPDPtr pi, tcPDPtr pj) const {
-  if ( !checkType(pi, pj) || cutType != DELTAR ) return 0.0;
+  if ( !checkType(pi, pj) || cutType != Cut::DELTAR ) return 0.0;
   return theCut;
 }
 
 Energy MadGraphTwoCut::minKTClus(tcPDPtr, tcPDPtr) const {
   return ZERO;
 }
 
 double MadGraphTwoCut::minDurham(tcPDPtr, tcPDPtr) const {
   return 0.0;
 }
 
 bool MadGraphTwoCut::passCuts(tcCutsPtr, tcPDPtr pitype, tcPDPtr pjtype,
 			      LorentzMomentum pi, LorentzMomentum pj,
 			      bool inci, bool incj) const {
   if ( inci || incj || !checkType(pitype, pjtype) ) return true;
-  if ( cutType == INVMASS ) return (pi + pj).m2() > sqr(theCut*GeV);
-  if ( cutType == DELTAR ) {
+  if ( cutType == Cut::INVMASS ) return (pi + pj).m2() > sqr(theCut*GeV);
+  if ( cutType == Cut::DELTAR ) {
     double deta2 = sqr(pi.eta() - pj.eta());
     double dphi = abs(pi.phi() - pj.phi());
     if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
     return sqrt(deta2 + sqr(dphi)) > theCut;
   }
   return true;
 }
 
-MadGraphTwoCut::PType MadGraphTwoCut::getType(tcPDPtr p) const {
+MadGraphTwoCut::P MadGraphTwoCut::getType(tcPDPtr p) const {
   switch ( abs(p->id()) ) {
   case ParticleID::d:
   case ParticleID::u:
   case ParticleID::s:
   case ParticleID::c:
   case ParticleID::g:
-    return JET;
+    return P::JET;
   case ParticleID::b:
-    return BOT;
+    return P::BOT;
   case ParticleID::gamma:
-    return PHO;
+    return P::PHO;
   case ParticleID::eminus:
   case ParticleID::nu_e:
   case ParticleID::muminus:
   case ParticleID::nu_mu:
   case ParticleID::tauminus:
   case ParticleID::nu_tau:
-    return LEP;
+    return P::LEP;
   default:
-    return NOT;
+    return P::NOT;
   }
 }
 
 bool MadGraphTwoCut::checkType(tcPDPtr pi, tcPDPtr pj) const {
   switch ( pairType ) {
-  case JETJET:
-    return getType(pi) == JET && getType(pj) == JET;
-  case LEPLEP:
-    if ( getType(pi) != LEP || getType(pj) != LEP ) return false;
-    if ( cutType == DELTAR ) return true;
-    // Special treatment for INVMASS.
+  case PP::JETJET:
+    return getType(pi) == P::JET && getType(pj) == P::JET;
+  case PP::LEPLEP:
+    if ( getType(pi) != P::LEP || getType(pj) != P::LEP ) return false;
+    if ( cutType == Cut::DELTAR ) return true;
+    // Special treatment for Cut::INVMASS.
     if ( pi->id()*pj->id() >= 0 ) return false;
     // OK we have a lepton-anti-lepton pair. I it the same lepton
     if ( pi->id() == -pj->id() ) return true;
     // NO, well is it the same family?
     if ( max(abs(pi->id()), abs(pj->id()))%2 ) return false;
     return abs(pi->id() + pj->id()) == 1 ;
-  case PHOPHO:
-    return getType(pi) == PHO && getType(pj) == PHO;
-  case BOTBOT:
-    return getType(pi) == BOT && getType(pj) == BOT;
-  case BOTJET:
-    return ( getType(pi) == BOT && getType(pj) == JET ) ||
-           ( getType(pi) == JET && getType(pj) == BOT );
-  case PHOJET:
-    return ( getType(pi) == PHO && getType(pj) == JET ) ||
-           ( getType(pi) == JET && getType(pj) == PHO );
-  case JETLEP:
-    return ( getType(pi) == LEP && getType(pj) == JET ) ||
-           ( getType(pi) == JET && getType(pj) == LEP );
-  case PHOBOT:
-    return ( getType(pi) == PHO && getType(pj) == BOT ) ||
-           ( getType(pi) == BOT && getType(pj) == PHO );
-  case BOTLEP:
-    return ( getType(pi) == BOT && getType(pj) == LEP ) ||
-           ( getType(pi) == LEP && getType(pj) == BOT );
-  case PHOLEP:
-    return ( getType(pi) == PHO && getType(pj) == LEP ) ||
-           ( getType(pi) == LEP && getType(pj) == PHO );
+  case PP::PHOPHO:
+    return getType(pi) == P::PHO && getType(pj) == P::PHO;
+  case PP::BOTBOT:
+    return getType(pi) == P::BOT && getType(pj) == P::BOT;
+  case PP::BOTJET:
+    return ( getType(pi) == P::BOT && getType(pj) == P::JET ) ||
+           ( getType(pi) == P::JET && getType(pj) == P::BOT );
+  case PP::PHOJET:
+    return ( getType(pi) == P::PHO && getType(pj) == P::JET ) ||
+           ( getType(pi) == P::JET && getType(pj) == P::PHO );
+  case PP::JETLEP:
+    return ( getType(pi) == P::LEP && getType(pj) == P::JET ) ||
+           ( getType(pi) == P::JET && getType(pj) == P::LEP );
+  case PP::PHOBOT:
+    return ( getType(pi) == P::PHO && getType(pj) == P::BOT ) ||
+           ( getType(pi) == P::BOT && getType(pj) == P::PHO );
+  case PP::BOTLEP:
+    return ( getType(pi) == P::BOT && getType(pj) == P::LEP ) ||
+           ( getType(pi) == P::LEP && getType(pj) == P::BOT );
+  case PP::PHOLEP:
+    return ( getType(pi) == P::PHO && getType(pj) == P::LEP ) ||
+           ( getType(pi) == P::LEP && getType(pj) == P::PHO );
   }
   return false;
 }
 
 void MadGraphTwoCut::persistentOutput(PersistentOStream & os) const {
   os << oenum(cutType) << oenum(pairType) << theCut;
 }
 
 void MadGraphTwoCut::persistentInput(PersistentIStream & is, int) {
   is >> ienum(cutType) >> ienum(pairType) >> theCut;
 }
 
 ClassDescription<MadGraphTwoCut> MadGraphTwoCut::initMadGraphTwoCut;
 // Definition of the static class description member.
 
 void MadGraphTwoCut::Init() {
 
   static ClassDocumentation<MadGraphTwoCut> documentation
     ("Objects of the MadGraphTwoCut class can be created automatically by "
      "the MadGraphReader class when scanning event files for information "
      "about cuts. It is also possible to create objects by hand and use "
      "it as any other MadGraphTwoCut object.");
 
-  static Switch<MadGraphTwoCut,CutType> interfaceCutType
+  static Switch<MadGraphTwoCut,Cut> interfaceCutType
     ("CutType",
      "The kind of cut this object will do.",
-     &MadGraphTwoCut::cutType, DELTAR, true, false);
+     &MadGraphTwoCut::cutType, Cut::DELTAR, true, false);
   static SwitchOption interfaceCutTypeInvariantMass
     (interfaceCutType,
      "InvariantMass",
      "The minimum invariant mass of two particles.",
-     INVMASS);
+     Cut::INVMASS);
   static SwitchOption interfaceCutTypeDeltaR
     (interfaceCutType,
      "DeltaR",
      "The minimum pseudo-rapidity--azimuth-angle distance between two "
      "particles.",
-     DELTAR);
+     Cut::DELTAR);
 
-  static Switch<MadGraphTwoCut,PPType> interfacePairType
+  static Switch<MadGraphTwoCut,PP> interfacePairType
     ("PairType",
      "The type of particle pairs this cut is applied to.",
-     &MadGraphTwoCut::pairType, JETJET, true, false);
+     &MadGraphTwoCut::pairType, PP::JETJET, true, false);
   static SwitchOption interfacePairTypeJetJet
     (interfacePairType,
      "JetJet",
      "The cut applies only to pairs of coloured particles (jets).",
-     JETJET);
+     PP::JETJET);
   static SwitchOption interfacePairTypeLeptonLepton
     (interfacePairType,
      "LeptonLepton",
      "The cut applies only to lepton pairs (in case of invariant mass, "
      "lepton--anti-lepton pairs of same flavour).",
-     LEPLEP);
+     PP::LEPLEP);
   static SwitchOption interfacePairTypePhotonPhoton
     (interfacePairType,
      "PhotonPhoton",
      "The cut applies only to pairs photons.",
-     PHOPHO);
+     PP::PHOPHO);
   static SwitchOption interfacePairTypeBottomPairs
     (interfacePairType,
      "BottomPairs",
      "The cut applies only to pairs of bottom quarks.",
-     BOTBOT);
+     PP::BOTBOT);
   static SwitchOption interfacePairTypeJetBottom
     (interfacePairType,
      "JetBottom",
      "The cut applies only to bottom quarks paired with another coloured "
      "particle (jet).",
-     BOTJET);
+     PP::BOTJET);
   static SwitchOption interfacePairTypePhotonJet
     (interfacePairType,
      "PhotonJet",
      "The cut applies only to a photon paired with a coloured particle (jet).",
-     PHOJET);
+     PP::PHOJET);
   static SwitchOption interfacePairTypeJetLepton
     (interfacePairType,
      "JetLepton",
      "The cut applies only to a coloured particle (jet) paired with a lepton.",
-     JETLEP);
+     PP::JETLEP);
   static SwitchOption interfacePairTypePhotonBottom
     (interfacePairType,
      "PhotonBottom",
      "The cut applies only to a photon paired with a bottom quark.",
-     PHOBOT);
+     PP::PHOBOT);
   static SwitchOption interfacePairTypeBottomLepton
     (interfacePairType,
      "BottomLepton",
      "The cut applies only to bottom quarks paired with a lepton.",
-     BOTLEP);
+     PP::BOTLEP);
   static SwitchOption interfacePairTypePhotonLepton
     (interfacePairType,
      "PhotonLepton",
      "The cut applies only to a photon paired with a lepton.",
-     PHOLEP);
+     PP::PHOLEP);
 
   static Parameter<MadGraphTwoCut,double> interfaceCut
     ("Cut",
      "The value of the cut to be applied (in units of GeV in case of "
      "minimum invariant mass).",
      &MadGraphTwoCut::theCut, 0.0, 0.0, 0,
      true, false, Interface::lowerlim);
 
   interfaceCut.rank(10);
   interfaceCutType.rank(9);
   interfacePairType.rank(8);
   interfaceCut.setHasDefault(false);
   interfaceCutType.setHasDefault(false);
   interfacePairType.setHasDefault(false);
 
 }
 
diff --git a/LesHouches/MadGraphTwoCut.h b/LesHouches/MadGraphTwoCut.h
--- a/LesHouches/MadGraphTwoCut.h
+++ b/LesHouches/MadGraphTwoCut.h
@@ -1,270 +1,270 @@
 // -*- C++ -*-
 //
 // MadGraphTwoCut.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_MadGraphTwoCut_H
 #define THEPEG_MadGraphTwoCut_H
 //
 // This is the declaration of the MadGraphTwoCut class.
 //
 
 #include "ThePEG/Cuts/TwoCutBase.h"
 
 namespace ThePEG {
 
 /**
  * Objects of the MadGraphTwoCut class can be created automatically by
  * the MadGraphReader class when scanning event files for information
  * about cuts. It is also possible to create objects by hand and use
  * it as any other OneCutBase object.
  *
  * @see \ref MadGraphTwoCutInterfaces "The interfaces"
  * defined for MadGraphTwoCut.
  */
 class MadGraphTwoCut: public TwoCutBase {
 
 public:
 
   /**
    * Enumerate the different kinds of cuts made by MadGraph.
    */
-  enum CutType {
+  enum class Cut {
     INVMASS, /**< The minimum invariant mass of two particles. */
     DELTAR   /**< The minimum pseudo-rapidity--azimuth-angle distance
                   between two particles. */
   };
 
   /**
    * Enumerate the types of particles the cut is made on.
    */
-  enum PType {
+  enum class P {
     JET, /**< Coloured particles (jets). */
     LEP, /**< Leptons. */
     PHO, /**< Photons. */
     BOT,  /**< Bottom quarks. */
     NOT  /**< Other types not cut on. */
   };
 
   /**
    * Enumerate the types of particles pairs the cut is made on.
    */
-  enum PPType {
+  enum class PP {
     JETJET, /**< The cut applies only to pairs of coloured particles (jets). */
     LEPLEP, /**< The cut applies only to lepton pairs (in case of INVMASS
                  lepton--anti-lepton pairs of same flavour). */
     PHOPHO, /**< The cut applies only to pairs photons. */
     BOTBOT, /**< The cut applies only to pairs of bottom quarks. */
     BOTJET, /**< The cut applies only to bottom quarks paired with another
                  coloured particle (jet). */
     PHOJET, /**< The cut applies only to a photon paired with a coloured
                  particle (jet). */
     JETLEP, /**< The cut applies only to a coloured particle (jet) paired
                  with a lepton. */
     PHOBOT, /**< The cut applies only to a photon paired with a bottom quark. */
     BOTLEP, /**< The cut applies only to bottom quarks paired with a lepton. */
     PHOLEP  /**< The cut applies only to a photon paired with a lepton. */
   };
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * The default constructor.
    */
   MadGraphTwoCut()
-    : cutType(DELTAR), pairType(JETJET), theCut(0.0) {}
+    : cutType(Cut::DELTAR), pairType(PP::JETJET), theCut(0.0) {}
 
   /**
    * The constructor used by the MadGraphReader.
    * @param t is the type of the cut.
    * @param p is the type of particles the cut is applied to.
    * @param c is the value of the cut (in units of GeV where applicable).
    */
-  MadGraphTwoCut(CutType t, PPType p, double c)
+  MadGraphTwoCut(Cut t, PP p, double c)
     : cutType(t), pairType(p), theCut(c) {}
 
   //@}
 
 public:
 
   /** @name Virtual functions mandated by the base class. */
   //@{
   /**
    * Return the minimum allowed squared invariant mass of two outgoing
    * partons of type \a pi and \a pj.
    */
   virtual Energy2 minSij(tcPDPtr pi, tcPDPtr pj) const;
 
   /**
    * Return the minimum allowed value of the negative of the squared
    * invariant mass of an incoming parton of type \a pi and an
    * outgoing parton of type \a po.
    */
   virtual Energy2 minTij(tcPDPtr pi, tcPDPtr po) const;
 
   /**
    * Return the minimum allowed value of \f$\Delta
    * R_{ij}=\sqrt{\Delta\eta_{ij}^2+\Delta\phi_{ij}^2}\f$ of two
    * outgoing partons of type \a pi and \a pj.
    */
   virtual double minDeltaR(tcPDPtr pi, tcPDPtr pj) const;
 
   /**
    * Return the minimum allowed value of the longitudinally invariant
    * \f$k_\perp\f$-algorithms distance measure. This is defined as
    * \f$\min(p_{\perp i}, p_{\perp
    * j})\sqrt{\Delta\eta_{ij}^2+\Delta\phi_{ij}^2}\f$ for two outgoing
    * partons, or simply \f$p_{\perp i}\f$ or \f$p_{\perp j}\f$ for a
    * single outgoing parton. Returns 0 if both partons are incoming. A
    * null pointer indicates an incoming parton, hence the type of the
    * incoming parton is irrelevant.
    */
   virtual Energy minKTClus(tcPDPtr pi, tcPDPtr pj) const;
 
   /**
    * Return the minimum allowed value of the Durham
    * \f$k_\perp\f$-algorithms distance measure. This is defined as
    * \f$2\min(E_j^2, E_j^2)(1-\cos\theta_{ij})/\hat{s}\f$ for two
    * outgoing partons.
    */
   virtual double minDurham(tcPDPtr pi, tcPDPtr pj) const;
 
   /**
    * Return true if a pair of particles with type \a pitype and \a
    * pjtype and momenta \a pi and \a pj respectively passes the
    * cuts. \a inci and \a inj indicates if the corresponding particles
    * are incoming.
    */
   virtual bool passCuts(tcCutsPtr parent, tcPDPtr pitype, tcPDPtr pjtype,
 			LorentzMomentum pi, LorentzMomentum pj,
 			bool inci = false, bool incj = false) const;
   //@}
 
 protected:
 
   /**
    * Returns true if cut should be applied to pair of particles of
    * type \a pi and \a pj.
    */
   bool checkType(tcPDPtr pi, tcPDPtr pj) const;
 
   /**
    * Get the type of particle \a p.
    */
-  PType getType(tcPDPtr p) const;
+  P getType(tcPDPtr p) 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;
   //@}
 
 private:
 
   /**
    * The type of this cut.
    */
-  CutType cutType;
+  Cut cutType;
 
   /**
    * The type of particle pairs this cut applies to.
    */
-  PPType pairType;
+  PP pairType;
 
   /**
    * The value of the cut to be applied.
    */
   double theCut;
 
 private:
 
   /**
    * The static object used to initialize the description of this class.
    * Indicates that this is a concrete class with persistent data.
    */
   static ClassDescription<MadGraphTwoCut> initMadGraphTwoCut;
 
   /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
   MadGraphTwoCut & operator=(const MadGraphTwoCut &);
 
 };
 
 }
 
 #include "ThePEG/Utilities/ClassTraits.h"
 
 namespace ThePEG {
 
 /** @cond TRAITSPECIALIZATIONS */
 
 /** This template specialization informs ThePEG about the
  *  base classes of MadGraphTwoCut. */
 template <>
 struct BaseClassTrait<MadGraphTwoCut,1> {
   /** Typedef of the first base class of MadGraphTwoCut. */
   typedef TwoCutBase NthBase;
 };
 
 /** This template specialization informs ThePEG about the name of
  *  the MadGraphTwoCut class and the shared object where it is defined. */
 template <>
 struct ClassTraits<MadGraphTwoCut>
   : public ClassTraitsBase<MadGraphTwoCut> {
   /** Return a platform-independent class name */
   static string className() { return "ThePEG::MadGraphTwoCut"; }
   /** Return the name(s) of the shared library (or libraries) be loaded to get
    *  access to the MadGraphTwoCut class and any other class on which it depends
    *  (except the base class). */
   static string library() { return "MadGraphReader.so"; }
 };
 
 /** @endcond */
 
 }
 
 #endif /* THEPEG_MadGraphTwoCut_H */
diff --git a/Persistency/PersistentIStream.h b/Persistency/PersistentIStream.h
--- a/Persistency/PersistentIStream.h
+++ b/Persistency/PersistentIStream.h
@@ -1,659 +1,667 @@
 // -*- C++ -*-
 //
 // PersistentIStream.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_PersistentIStream_H
 #define ThePEG_PersistentIStream_H
 // This is the declaration of the PersistentIStream class.
 
 #include "ThePEG/Config/ThePEG.h"
 #include "InputDescription.h"
 #include "PersistentIStream.fh"
 #include "ThePEG/Utilities/Exception.h"
 #include <climits>
 #include <valarray>
 
 namespace ThePEG {
 
 /** @ingroup Persistency
  * PersistentIStream is used to read persistent objects from a stream
  * where they were previously written using PersistentOStream. Basic
  * types and pointers to objects derived from
  * <code>PersistentBase</code> should be read in the same order they
  * were written out. If <code>pedantic()</code> is true the same
  * classes that were written out must be present in the current
  * program. If <code>pedantic()</code> is false and if an object is
  * read for which only a base class is present in the current program,
  * only the parts corresponding to the base class will be read, and
  * the rest will be gracefully skipped.
  *
  * Each base class of a given object will be asked to read its
  * members from the stream starting from the least derived class going to
  * the most derived one. Members may be pointers to other persistent
  * objects or basic types or containers of these. The output for each
  * object part should be implemented by specializing the
  * ClassTraits<T>::input method, which otherwise
  * will call the non-virtual <code>persistentInput</code> function of
  * the class. Note that for diamond-shaped multiple inheritance
  * structures, the virtual base classes will be written out several
  * times for the same object.
  *
  * @see PersistentOStream
  * @see ClassTraits
  */
 class PersistentIStream {
 
 public:
 
   ThePEG_DECLARE_POINTERS(PersistentBase,BPtr);
 
   /** A vector of pointers to persistent objects */
   typedef vector<BPtr> ObjectVector;
 
   /** A vector of bare pointers to InputDescription objects. */
   typedef InputDescription::DescriptionVector DescriptionVector;
 
 public:
 
   /**
    * Constuctor giving an input stream to be used as an underlying
    * istream.
    */
   PersistentIStream(istream & is) 
     : theIStream(&is), isPedantic(true), 
       allocStream(false), badState(false)
   {
     init();
   }
 
 
 
   /**
    * Constuctor giving a file name to read from. If the first
    * character in the string is a '|', the corresponding program is
    * run and its standard output is used instead. If the filename ends
    * in ".gz" the file is uncompressed with gzip.
    */
   PersistentIStream(string);
 
   /**
    * The destructor.
    */
   ~PersistentIStream();
 
   /**
    * Operator for extracting persistent objects from the stream.
    * @param ptr this pointer will refer to the extracted object.
    * @return a reference to the stream.
    */
   template <typename T>
   PersistentIStream & operator>>(RCPtr<T> & ptr) {
     BPtr b = getObject();
     ptr = dynamic_ptr_cast< RCPtr<T> >(b);
     if ( b && !ptr ) setBadState();
     return *this;
   }
 
   /**
    * Operator for extracting persistent objects from the stream.
    * @param ptr this pointer will refer to the extracted object.
    * @return a reference to the stream.
    */
   template <typename T>
   PersistentIStream & operator>>(ConstRCPtr<T> & ptr) {
     BPtr b = getObject();
     ptr = dynamic_ptr_cast< ConstRCPtr<T> >(b);
     if ( b && !ptr ) setBadState();
     return *this;
   }
 
   /**
    * Operator for extracting persistent objects from the stream.
    * @param ptr this pointer will refer to the extracted object.
    * @return a reference to the stream.
    */
   template <typename T>
   PersistentIStream & operator>>(TransientRCPtr<T> & ptr) {
     BPtr b = getObject();
     ptr = dynamic_ptr_cast< TransientRCPtr<T> >(b);
     if ( b && !ptr ) setBadState();
     return *this;
   }
 
   /**
    * Operator for extracting persistent objects from the stream.
    * @param ptr this pointer will refer to the extracted object.
    * @return a reference to the stream.
    */
   template <typename T>
   PersistentIStream & operator>>(TransientConstRCPtr<T> & ptr) {
     BPtr b = getObject();
     ptr = dynamic_ptr_cast< TransientConstRCPtr<T> >(b);
     if ( b && !ptr ) setBadState();
     return *this;
   }
 
 
   /** @name Operators for extracting built-in types from the stream. */
   //@{
   /**
    * Read a character string.
    */
   PersistentIStream & operator>>(string &);
 
   /**
    * Read a character.
    */
   PersistentIStream & operator>>(char &);
 
   /**
    * Read a signed character.
    */
   PersistentIStream & operator>>(signed char &);
 
   /**
    * Read an unsigned character.
    */
   PersistentIStream & operator>>(unsigned char &);
 
   /**
    * Read an integer.
    */
   PersistentIStream & operator>>(int & i) {
     is() >> i;
     getSep();
     return *this;
   }
 
   /**
    * Read an unsigned integer.
    */
   PersistentIStream & operator>>(unsigned int & i) {
     is() >> i;
     getSep();
     return *this;
   }
 
   /**
    * Read a long integer.
    */
   PersistentIStream & operator>>(long & i) {
     is() >> i;
     getSep();
     return *this;
   }
 
   /**
    * Read an unsigned long integer.
    */
   PersistentIStream & operator>>(unsigned long & i) {
     is() >> i;
     getSep();
     return *this;
   }
 
   /**
    * Read a short integer.
    */
   PersistentIStream & operator>>(short & i) {
     is() >> i;
     getSep();
     return *this;
   }
 
   /**
    * Read an unsigned short integer.
    */
   PersistentIStream & operator>>(unsigned short & i) {
     is() >> i;
     getSep();
     return *this;
   }
 
   /**
    * Read a double.
    */
   PersistentIStream & operator>>(double & d) {
     is() >> d;
     getSep();
     return *this;
   }
 
   /**
    * Read a float.
    */
   PersistentIStream & operator>>(float & f) {
     is() >> f;
     getSep();
     return *this;
   }
 
   /**
    * Read a bool.
    */
   PersistentIStream & operator>>(bool &);
 
   /**
    * Read a Complex.
    */
   PersistentIStream & operator>>(Complex &);
   //@}
 
   /**
    * Intput of containers streamable objects.
    * @param c the container into which objects are added.
    */
   template <typename Container> void getContainer(Container & c) {
     long size;
     typename Container::value_type val;
     c.clear();
     *this >> size;
     while ( size-- && good() ) {
       *this >> val;
       c.insert(c.end(), val);
     }
   }
 
   /**
    * Read in an object. Create an object and read its data from the
    * stream.
    * @return a pointer to the read object.
    */
   BPtr getObject();
 
   /**
    * For a given object, read the member variables corresponding to a
    * given InputDescription object.
    * @param obj the object to be read into.
    * @param pid a pointer to an InputDescription describing the
    * (sub)class to be read.
    */
   void getObjectPart(tBPtr obj, const InputDescription * pid);
 
   /**
    * Read a class description from the underlying stream and return a
    * corresponding InputDescription object
    */
   const InputDescription * getClass();
   
   /**
    * Set pedantic mode.  If the stream is set to be tolerant it is
    * allowed to read an object from the stream even if the
    * corresponding class is not known to the running executable, under
    * the condition that a public base class of the unknown class is
    * known. If the stream is set to be pedantic this is not allowed.
    * By default, the stream is pedantic.
    */
   PersistentIStream & setPedantic() {
     isPedantic = true;
     return *this;
   }
 
   /**
    * Set tolerant mode.  If the stream is set to be tolerant it is
    * allowed to read an object from the stream even if the
    * corresponding class is not known to the running executable, under
    * the condition that a public base class of the unknown class is
    * known. If the stream is set to be pedantic this is not allowed.
    * By default, the stream is pedantic.
    */
   PersistentIStream & setTolerant() {
     isPedantic = false;
     return *this;
   }
 
   /**
    * Check the state of the stream.
    */
   bool good() const { return !badState && is(); }
   
   /**
    * Check the state of the stream.
    */
   bool operator!() const { return !good(); }
 
   /**
    * Check the state of the stream.
    */
   operator bool() const { return good(); }
 
   /**
    * Check the tolerance. Returns true if setPedantic() has been
    * called or if not setTolerant() has been called.
    */
   bool pedantic() const { return isPedantic; }
 
   /**
    * The global libraries loaded on initialization.
    */
   const vector<string> & globalLibraries() const {
     return theGlobalLibraries;
   }
 
 private:
 
   /** @cond EXCEPTIONCLASSES */
   /** @ingroup Persistency
       Thrown if a class is missing */
   struct MissingClass: public Exception {};
 
   /** @ingroup Persistency Thrown if an object which should have been
       read in is missing. */
   struct MissingObject: public Exception {};
 
   /** @ingroup Persistency Thrown if reading from the stream failed
       for some reason. */
   struct ReadFailure: public Exception {};
   /** @endcond */
 
   /**
    * Internal initialization.
    */
   void init();
 
   /**
    * Get the next character from the associated istream.
    */
   char get() { return is().get(); }
 
   /**
    * Get the next character from the associated istream and decode it
    * if it is escaped.
    */
   char escaped() {
     char c = get();
     return c == tNoSep? tSep: c;
   }
 
   /**
    * Set the stream in a bad state
    */
   void setBadState() {
     breakThePEG();
     badState = true;
   }
 
   /**
    * Read a field separator from the stream.
    */
   void getSep() {
     if ( !pedantic() ) skipField();
     else if ( get() != tSep ) setBadState();
   }
 
   /**
    * Scan the stream for the next field separator.
    */
   void skipField() {
     is().ignore(INT_MAX, tSep);
     if ( !is() ) setBadState();
   }
 
 
   /**
    * Check if the next char to be read is a tBegin marker.
    */
   bool beginObject() { return is().peek() == tBegin; }
 
   /**
    * Scan the stream to the end of the current object. If any new object are
    * found these are read prom the stream to ensure that the pointer structure
    * is preserved.
    */
   void endObject();
 
   /**
    * Scan stream for "end base class" marker. The \a classname is the
    * name of the class currently being read and is only used for
    * documenting exceptions.
    */
   void endBase(string classname);
 
   /**
    * Return a reference to the associated stream.
    */
   istream & is() { return *theIStream; }
 
   /**
    * Return a const reference to the associated stream.
    */
   const istream & is() const { return *theIStream; }
 
   /**
    * Lists of objects that have been read.
    */
   ObjectVector readObjects;
 
   /**
    * Lists of classes and corresponding version strings that have been read.
    */
   DescriptionVector readClasses;
 
   /**
    * A pointer to the associated istream.
    */
   istream * theIStream;
 
   /**
    * Pedantic or tolerant. See description of the setPedantic() and
    * setTolerant() methods.
    */
   bool isPedantic;
 
   /**
    * True if the associated istream should be deleted when the PersistentIStream
    * is destroyed.
    */
   bool allocStream;
 
   /**
    * False if no errors has occurred.
    */
   bool badState;
 
   /** Version number of the PersistentOStream which has written the
    *  file being read. */
   int version;
 
   /** Subversion number of the PersistentOStream which has written the
    *  file being read. */
   int subVersion;
 
   /**
    * Global libraries loaded in the initialization.
    */
   vector<string> theGlobalLibraries;
 
   /** @name Special marker characters */
   //@{
   /**
    * The special marker character indicating the beginning of an object.
    */
   static const char tBegin = '{';
 
   /**
    * The special marker character indicating the end of an object.
    */
   static const char tEnd = '}';
 
   /**
    * The marker character indicating the beginning of the next base
    * class in case of multiple inheritance.
    */
   static const char tNext = '|';
 
   /**
    * The special marker character indicating an escaped marker character.
    */
   static const char tNull = '\\';
 
   /**
    * The special marker character indicating the end of a value.
    */
   static const char tSep = '\n';
 
   /**
    * The special marker character used to avoid confusion with escaped
    * tSep markers.
    */
   static const char tNoSep = 'n';
 
   /**
    * The special marker character indicating a true boolean value.
    */
   static const char tYes = 'y';
 
   /**
    * The special marker character indicating a false boolean value.
    */
   static const char tNo = 'n';
   //@}
 
 private:
 
   /**
    * Standard ctors and assignment are private and not implemented.
    */
   PersistentIStream();
 
   /**
    * Standard ctors and assignment are private and not implemented.
    */
   PersistentIStream(const PersistentIStream &);
 
   /**
    * Standard ctors and assignment are private and not implemented.
    */
   PersistentIStream & operator=(const PersistentIStream &);
 
 };
 
 
 /**
  * Operator for applying manipulators to the stream.
  */
 inline PersistentIStream & operator>>(PersistentIStream & is, 
 				      PersistentIManip func) {
   return (*func)(is);
 }
   
 /**
  * The manipulator for setting pedantic mode.
  */
 inline PersistentIStream & pedantic(PersistentIStream & is) {
   return is.setPedantic();
 }
 
 
 /**
  * The manipulator for setting tolerant mode.
  */
 inline PersistentIStream & tolerant(PersistentIStream & is) {
   return is.setTolerant();
 }
 
 
 /**
  * @name Partial specializations of operator>> for input of
  * std::containers.
  */
 //@{
 /** Input a pair of objects. */
 template <typename T1, typename T2>
 inline PersistentIStream & operator>>(PersistentIStream & is, pair<T1,T2> & p) {
   return is >> p.first >> p.second;
 }
 
 /** Input a map of key/objects pairs. */
 template <typename Key, typename T, typename Cmp, typename A>
 inline PersistentIStream & operator>>(PersistentIStream & is, map<Key,T,Cmp,A> & m) {
   m.clear();
   long size;
   Key k;
   is >> size;
   while ( size-- && is ) {
     is >> k;
     is >> m[k];
   }
   return is;
 }
 
 /** Input a multimap of key/objects pairs. */
 template <typename Key, typename T, typename Cmp, typename A>
 inline PersistentIStream & operator>>(PersistentIStream & is,
 				      multimap<Key,T,Cmp,A> & m) {
   m.clear();
   long size;
   Key k;
   T t;
   is >> size;
   while ( size-- && is ) {
     is >> k;
     is >> t;
     m.insert(make_pair(k, t));
   }
   return is;
 }
 
 
 /** Input a set of objects. */
 template <typename Key, typename Cmp, typename A>
 inline PersistentIStream & operator>>(PersistentIStream & is, 
 				      set<Key,Cmp,A> & s) {
   is.getContainer(s);
   return is;
 }
 
 /** Input a multoset of objects. */
 template <typename Key, typename Cmp, typename A>
 inline PersistentIStream & operator>>(PersistentIStream & is,
 				      multiset<Key,Cmp,A> & s) {
   is.getContainer(s);
   return is;
 }
 
 
 /** Input a list of objects. */
 template <typename T, typename A>
 inline PersistentIStream & operator>>(PersistentIStream & is, 
 				      list<T,A> & l) {
   is.getContainer(l);
   return is;
 }
 
 
 /** Input a vector of objects. */
 template <typename T, typename A>
 inline PersistentIStream & operator>>(PersistentIStream & is, 
-				      vector<T,A> & v) {
+              vector<T,A> & v) {
   is.getContainer(v);
   return is;
 }
 
+/** Input an array of objects. */
+template <typename T, size_t N>
+inline PersistentIStream & operator>>(PersistentIStream & is, 
+              array<T,N> & a) {
+  for ( size_t i = 0; i < N && is.good(); ++i ) 
+    is >> a[i];
+  return is;
+}
 
 /** Input a deque of objects. */
 template <typename T, typename A>
 inline PersistentIStream & operator>>(PersistentIStream & is, 
 				      deque<T,A> & d) {
   is.getContainer(d);
   return is;
 }
 
-/** Input a deque of objects. */
+/** Input a valarray. */
 template <typename T>
 inline PersistentIStream & operator>>(PersistentIStream & is, 
 				      std::valarray<T> & v) {
   long size;
   is >> size;
   v = std::valarray<T>(size);
   for ( int i = 0; i < size && is.good(); ++i ) is >> v[i];
   return is;
 }
 
 }
 
 #endif /* ThePEG_PersistentIStream_H */
diff --git a/Persistency/PersistentOStream.h b/Persistency/PersistentOStream.h
--- a/Persistency/PersistentOStream.h
+++ b/Persistency/PersistentOStream.h
@@ -1,660 +1,671 @@
 // -*- C++ -*-
 //
 // PersistentOStream.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_PersistentOStream_H
 #define ThePEG_PersistentOStream_H
 // This is the declaration of the PersistentOStream class.
 
 #include "ThePEG/Config/ThePEG.h"
 #include "ThePEG/Utilities/ClassDescription.h"
 #include "ThePEG/Utilities/Exception.h"
 #include "ThePEG/Utilities/Debug.h"
 #include "PersistentOStream.fh"
 #include "PersistentOStream.xh"
 #include <valarray>
 
 namespace ThePEG {
 
 /** @ingroup Persistency
  * PersistentOStream is used to write objects persistently
  * to a stream from which they can be read in again with a
  * PersistentIStream. Pointers to objects of classes
  * derived from <code>PersistentBase</code> may be written out if a
  * static ClassDescription object is present for the
  * class. Also basic types may be written to the stream, as well as
  * containers of pointers to persistent objects and basic types.
  *
  * The <code>PersistentOStream</code> keeps a list of all pointers to
  * written persistent objects, so that if several pointers to the
  * smame object is written, the object will only be written once.
  *
  * Each base class of a given object will be asked to write its
  * members to the stream starting from the least derived class going to
  * the most derived one. Members may be pointers to other persistent
  * objects or basic types or containers of these. The output for each
  * object part should be implemented by specializing the
  * ClassTraits<T>::output method, which otherwise
  * will call the non-virtual <code>persistentOutput</code> function of
  * the class. Note that for diamond-shaped multiple inheritance
  * structures, the virtual base classes will be written out several
  * times for the same object.
  *
  * @see PersistentIStream
  * @see ClassDescription
  * @see ClassTraits
  */
 class PersistentOStream {
 
 public:
 
   ThePEG_DECLARE_POINTERS(PersistentBase,BPtr);
 
   /** A map of objects indexed by integers */
   ThePEG_DECLARE_MAP(cBPtr, int, ObjectMap);
 
   /** A map relating class descriptions to integers. */
   ThePEG_DECLARE_MAP(const ClassDescriptionBase *, int, ClassMap);
 
   /** A vector of bare pointers to InputDescription objects. */
   typedef ClassDescriptionBase::DescriptionVector DescriptionVector;
 
 public:
 
   /**
    * Constuctor giving an output stream. Optionally a vector of
    * libraries to be loaded before the resulting file can be read in
    * again can be given in \a libs.
    */
   PersistentOStream(ostream &, const vector<string> & libs = vector<string>());
 
   /**
    * Constuctor giving a file name to read. If the first
    * character in the string is a '|', the corresponding program is
    * run and its standard input is used instead. If the filename ends
    * in ".gz" the file is compressed with gzip. Optionally a vector of
    * libraries to be loaded before the resulting file can be read in
    * again can be given in \a libs.
    */
   PersistentOStream(string, const vector<string> & libs = vector<string>());
 
   /**
    * The destructor
    */
   ~PersistentOStream();
 
   /**
    * Operator for writing persistent objects to the stream.
    * @param p a pointer to the object to be written.
    * @return a reference to the stream.
    */
   template <typename T>
   PersistentOStream & operator<<(const RCPtr<T> & p) { 
     return outputPointer(p); 
   }
 
   /**
    * Operator for writing persistent objects to the stream.
    * @param p a pointer to the object to be written.
    * @return a reference to the stream.
    */
   template <typename T>
   PersistentOStream & operator<<(const ConstRCPtr<T> & p) { 
     return outputPointer(p); 
   }
 
   /**
    * Operator for writing persistent objects to the stream.
    * @param p a pointer to the object to be written.
    * @return a reference to the stream.
    */
   template <typename T>
   PersistentOStream & operator<<(const TransientRCPtr<T> & p) { 
     return outputPointer(p); 
   }
 
   /**
    * Operator for writing persistent objects to the stream.
    * @param p a pointer to the object to be written.
    * @return a reference to the stream.
    */
   template <typename T>
   PersistentOStream & operator<<(const TransientConstRCPtr<T> & p) {
     return outputPointer(p);
   }
 
 
   /** @name Operators for extracting built-in types from the stream. */
   //@{
   /**
    * Write a character string.
    */
   PersistentOStream & operator<<(string s) {
     for ( string::const_iterator i = s.begin(); i < s.end(); ++i ) escape(*i);
     put(tSep);
     return *this;
   }
 
   /**
    * Write a character.
    */
   PersistentOStream & operator<<(char c) {
     escape(c);
     put(tSep);
     return *this;
   }
 
   /**
    * Write a signed character.
    */
   PersistentOStream & operator<<(signed char c) {
     return (*this) << static_cast<char>(c);
   }
   
   /**
    * Write an unsigned character.
    */
   PersistentOStream & operator<<(unsigned char c) {
     return (*this) << static_cast<char>(c);
   }
   
   /**
    * Write an integer.
    */
   PersistentOStream & operator<<(int i) {
     os() << i;
     put(tSep);
     return *this;
   }
 
   /**
    * Write an unsigned integer.
    */
   PersistentOStream & operator<<(unsigned int i) {
     os() << i;
     put(tSep);
     return *this;
   }
 
   /**
    * Write a long integer.
    */
   PersistentOStream & operator<<(long i) {
     os() << i;
     put(tSep);
     return *this;
   }
 
   /**
    * Write an unsigned long integer.
    */
   PersistentOStream & operator<<(unsigned long i) {
     os() << i;
     put(tSep);
     return *this;
   }
 
   /**
    * Write a short integer.
    */
   PersistentOStream & operator<<(short i) {
     os() << i;
     put(tSep);
     return *this;
   }
 
   /**
    * Write an unsigned short integer.
    */
   PersistentOStream & operator<<(unsigned short i) {
     os() << i;
     put(tSep);
     return *this;
   }
 
   /**
    * Write a double.
    */
   PersistentOStream & operator<<(double d) {
     if ( ! isfinite(d) )
       throw WriteError()
 	<< "Tried to write a NaN or Inf double to a persistent stream."
 	<< Exception::runerror;
     os() << setprecision(18) << d;
     put(tSep);
     return *this;
   }
 
   /**
    * Write a float.
    */
   PersistentOStream & operator<<(float f) {
     if ( ! isfinite(f) )
       throw WriteError()
 	<< "Tried to write a NaN or Inf float to a persistent stream."
 	<< Exception::runerror;
     os() << setprecision(9) << f;
     put(tSep);
     return *this;
   }
 
   /**
    * Write a boolean.
    */
   PersistentOStream & operator<<(bool t) {
     if (t) put(tYes);
     else put(tNo);
     // This is a workaround for a possible bug in gcc 4.0.0
     // which inserts tYes and tNo as global symbols although
     // they are private
     //  put(t? tYes: tNo);
     put(tSep);
     return *this;
   }
 
   /**
    * Write a c-style character string (to be read in as a std::string).
    */
   PersistentOStream & operator<<(const char * s) {
     *this << string(s);
     return *this;
   }
 
   /**
    * Write a Complex.
    */
   PersistentOStream & operator<<(Complex z) {
     *this << z.real() << z.imag();
     return *this;
   }
   //@}
 
   /**
    * Output of containers of streamable objects.
    */
   template <typename Container>
   void putContainer(const Container & c) {
     *this << c.size();
     for ( typename Container::const_iterator it = c.begin();
 	  it != c.end() && good() ; ++it )
       *this << *it;
   }
 
   /**
    * Write out a persistent object given a pointer to it.
    */
   PersistentOStream & outputPointer(tcBPtr);
 
   /**
    * For a given object, write the member variables corresponding to a
    * given ClassDescriptionBase object.
    * @param obj the object to be written.
    * @param cd a pointer to a ClassDescriptionBase describing the
    * (sub)class to written.
    */
   void putObjectPart(tcBPtr obj, const ClassDescriptionBase * cd);
 
   /**
    * Remove all objects that have been written, except those which are
    * to be saved, from the list of written objects.
    */
   PersistentOStream & flush();
   
   /**
    * Instuct the stream to save the following objects (protecting them from
    * being flushed).
    */
   PersistentOStream & push() {
     lastSavedObject.push(writtenObjects.size() - 1);
     return *this;
   }
 
   /**
    * Instuct the stream not to save the following objects.
    */
   PersistentOStream & pop() {
     lastSavedObject.pop();
     return *this;
   }
 
   /**
    * Check the state of the stream.
    */
   bool good() const { return !badState && os(); }
 
   /**
    * Check the state of the stream.
    */
   operator bool() const { return good(); }
 
   /**
    * Check the state of the stream.
    */
   bool operator!() const { return !good(); }
 
 private:
 
   /** @cond EXCEPTIONCLASSES */
   /** @ingroup Persistency
    * Internal exception class.
    */
   struct MissingClass: public Exception {};
   /** @ingroup Persistency
    * Internal exception class.
    */
   struct WriteError: public Exception {};
   /** @endcond */
 
   /**
    * The version of this PersistentOStream implementation.
    */
   static const int version = 0;
 
   /**
    * The subversion of this PersistentOStream implementation.
    */
   static const int subVersion = 3;
 
   /** @name Special marker characters */
   //@{
   /**
    * The special marker character indicating the beginning of an object.
    */
   static const char tBegin = '{';
 
   /**
    * The special marker character indicating the end of an object.
    */
   static const char tEnd = '}';
 
   /**
    * The marker character indicating the beginning of the next base
    * class in case of multiple inheritance.
    */
   /**
    * The special marker character indicating an escaped marker character.
    */
   static const char tNext = '|';
 
   /**
    * The special marker character indicating an escaped marker character.
    */
   static const char tNull = '\\';
 
   /**
    * The special marker character indicating the end of a value.
    */
   static const char tSep = '\n';
 
   /**
    * The special marker character used to avoid confusion with escaped
    * tSep markers.
    */
   static const char tNoSep = 'n';
 
   /**
    * The special marker character indicating a true boolean value.
    */
   static const char tYes = 'y';
 
   /**
    * The special marker character indicating a false boolean value.
    */
   static const char tNo = 'n';
   //@}
 
   /**
    * Return true if the given character is aspecial marker character.
    */
   bool isToken(char c) const {
     return c == tBegin || c == tEnd || c == tNext || c == tSep || c == tNull;
   }
 
   /**
    * Set the stream in a bad state.
    */
   void setBadState() {
     breakThePEG();
     badState = true;
   }
 
   /**
    * Check if the state is ok.
    */
   void checkState() { if ( ! os() ) badState = true; }
 
   /**
    * Write out class information to the associated ostream.
    */
   const ClassDescriptionBase * writeClassId(tcBPtr);
 
   /**
    * write out class information to the associated ostream.
    */
   void writeClassDescription(const ClassDescriptionBase *);
 
   /**
    * Put a "begin object" marker on the associated ostream
    */
   void beginObject() { put(tBegin); }
 
   /**
    * Put a "end of object" marker on the associated ostream
    */
   void endObject() { put(tEnd); }
 
   /**
    * Put an "next base class" marker on the associated ostream
    */
   void endBase() { put(tNext); }
 
   /**
    * Put a character on the associated ostream
    */
   void put(char c) { os().put(c); }
 
   /**
    * Put a character on the associated ostream but escape it if it is
    * a token.
    */
   void escape(char c) { 
     if ( isToken(c) ) {
       put(tNull);
       put( c == tSep? tNoSep: c );
     } else
       put(c);
   }
 
   /**
    * Return a reference to the associated ostream.
    */
   ostream & os() { return *theOStream; }
 
   /**
    * Return a const reference to the associated ostream.
    */
   const ostream & os() const { return *theOStream; }
 
   /**
    * Write out initial metainfo on the stream.
    */
   void init(const vector<string> & libs);
 
   /**
    * List of written objects.
    */
   ObjectMap writtenObjects;
 
   /**
    * List of written objects that are to be saved.
    */
   stack<int> lastSavedObject;
 
   /**
    * List of written classes.
    */
   ClassMap writtenClasses;
 
   /**
    * A pointer to the associated ostream.
    */
   ostream * theOStream;
 
   /**
    * True if no errors has occurred.
    */
   bool badState;
 
   /**
    * True if the associated ostream should be deleted in the destructor.
    */
   bool allocStream;
 
 private:
 
   /**
    * Standard ctors and assignment are private and not implemented.
    */
   PersistentOStream();
 
   /**
    * Standard ctors and assignment are private and not implemented.
    */
   PersistentOStream(const PersistentOStream &);
 
   /**
    * Standard ctors and assignment are private and not implemented.
    */
   PersistentOStream & operator=(const PersistentOStream &);
 
 };
 
 /**
  * Operator for applying manipulators to the stream.
  */
 inline PersistentOStream &
 operator<<(PersistentOStream & os, PersistentOManip func) {
   return (*func)(os);
 }
 
 
 /**
  * The manipulator for calling PersistentOStream::flush().
  */
 inline PersistentOStream & flush(PersistentOStream & os) { return os.flush(); }
 
 /**
  * The manipulator for calling PersistentOStream::push().
  */
 inline PersistentOStream & push(PersistentOStream & os) { return os.push(); }
 
 /**
  * The manipulator for calling PersistentOStream::pop().
  */
 inline PersistentOStream & pop(PersistentOStream & os) { return os.pop(); }
 
 
 /**
  * @name Partial specializations of operator<< for output of
  * std::containers.
  */
 //@{
 /** Output a pair of objects. */
 template <typename T1, typename T2>
 inline PersistentOStream & operator<<(PersistentOStream & os,
 				      const pair<T1,T2> & p) {
   return os << p.first << p.second;
 }
 
 /**
  * Output a multimap of key/object pairs.
  */
 template <typename Key, typename T, typename Cmp, typename A>
 inline PersistentOStream & operator<<(PersistentOStream & os,
 				      const multimap<Key,T,Cmp,A> & m) {
   os.putContainer(m);
   return os;
 }
 
 /**
  * Output a map of key/object pairs.
  */
 template <typename Key, typename T, typename Cmp, typename A>
 inline PersistentOStream & operator<<(PersistentOStream & os,
 				      const map<Key,T,Cmp,A> & m) {
   os.putContainer(m);
   return os;
 }
 
 /**
  * Output a set of objects.
  */
 template <typename Key, typename Cmp, typename A>
 inline PersistentOStream & operator<<(PersistentOStream & os,
 				      const set<Key,Cmp,A> & s) {
   os.putContainer(s);
   return os;
 }
 
 
 /**
  * Output a multiset of objects.
  */
 template <typename Key, typename Cmp, typename A>
 inline PersistentOStream & operator<<(PersistentOStream & os,
 				      const multiset<Key,Cmp,A> & s) {
   os.putContainer(s);
   return os;
 }
 
 
 /**
  * Output a list of objects.
  */
 template <typename T, typename A>
 inline PersistentOStream & operator<<(PersistentOStream & os,
 				      const list<T,A> & l) {
   os.putContainer(l);
   return os;
 }
 
 
 /**
  * Output a vector of objects.
  */
 template <typename T, typename A>
 inline PersistentOStream & operator<<(PersistentOStream & os,
-				      const vector<T,A> & v) {
+              const vector<T,A> & v) {
   os.putContainer(v);
   return os;
 }
 
 /**
+ * Output an array of objects.
+ */
+template <typename T, size_t N>
+inline PersistentOStream & operator<<(PersistentOStream & os,
+              const array<T,N> & a) {
+  for ( auto it = a.cbegin(); it != a.cend() && os.good() ; ++it )
+      os << *it;
+  return os;
+}
+
+/**
  * Output a deque of objects.
  */
 template <typename T, typename A>
 inline PersistentOStream & operator<<(PersistentOStream & os,
 				      const deque<T,A> & d) {
   os.putContainer(d);
   return os;
 }
 
 /**
- * Output a valarray of objects.
+ * Output a valarray.
  */
 template <typename T>
 inline PersistentOStream & operator<<(PersistentOStream & os,
 				      const std::valarray<T> & v) {
   os << v.size();
   for ( int i = 0, N = v.size(); i < N; ++i ) os << v[i];
   return os;
 }
 //@}
 
 }
 
 #endif /* ThePEG_PersistentOStream_H */
diff --git a/Repository/StandardRandom.h b/Repository/StandardRandom.h
--- a/Repository/StandardRandom.h
+++ b/Repository/StandardRandom.h
@@ -1,162 +1,162 @@
 // -*- C++ -*-
 //
 // StandardRandom.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_StandardRandom_H
 #define ThePEG_StandardRandom_H
 // This is the declaration of the StandardRandom class.
 
 #include "RandomGenerator.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 namespace ThePEG {
 
 /**
  * StandardRandom inherits from the RandomGenerator class and is an
  * interface to the CLHEP::JamesRandom engine class.
  *
  * @see \ref StandardRandomInterfaces "The interfaces"
  * defined for StandardRandom.
  */
 class StandardRandom: public RandomGenerator {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * Default constructor.
    */
-  StandardRandom() : u(97) { if ( theSeed != 0 ) setSeed(theSeed); }
+  StandardRandom() : u() { if ( theSeed != 0 ) setSeed(theSeed); }
   //@}
 
 public:
 
   /**
    * Reset the underlying random algorithm with the given seed. If the
    * \a seed is set to -1 a standard seed will be used.
    */
   virtual void setSeed(long seed);
 
 protected:
 
   /**
    * Fill the cache with random numbers.
    */
   virtual void fill();
 
 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);
   //@}
 
   /**
    * Standard Init function used to initialize the interface.
    */
   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;
   //@}
 
 private:
 
   /**
    * The internal state vector.
    */
-  vector<double> u;
+  array<double,97> u;
 
   /**
    * Parameter for the internal state.
    */
   double c;
 
   /**
    * Parameter for the internal state.
    */
   double cd;
 
   /**
    * Parameter for the internal state.
    */
   double cm;
 
   /**
    * Index for the internal state.
    */
   int i97;
 
   /**
    * Index for the internal state.
    */
   int j97;
 
 private:
 
   /**
    * Describe a concrete class with persistent data.
    */
   static ClassDescription<StandardRandom> initStandardRandom;
 
   /**
    *  Private and non-existent assignment operator.
    */
   StandardRandom & operator=(const StandardRandom &);
 
 };
 
 /** @cond TRAITSPECIALIZATIONS */
 
 /** This template specialization informs ThePEG about the base classes
  *  of StandardRandom. */
 template <>
 struct BaseClassTrait<StandardRandom,1>: public ClassTraitsType {
   /** Typedef of the first base class of StandardRandom. */
   typedef RandomGenerator NthBase;
 };
 
 /** This template specialization informs ThePEG about the name of the
  *  StandardRandom class. */
 template <>
 struct ClassTraits<StandardRandom>: public ClassTraitsBase<StandardRandom> {
   /** Return a platform-independent class name */
   static string className() { return "ThePEG::StandardRandom"; }
 };
 
 /** @endcond */
 
 }
 
 #endif /* ThePEG_StandardRandom_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.
  */
-struct SimplePhaseSpace {
+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>
-  static void CMS(Energy2 s, PType & p1, PType & p2);
+  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>
-  static void CMS(PType & p1, PType & p2, Energy2 s,
+  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>
-  static void CMS(PType & p1, PType & p2, Energy2 s, Energy2 t, double phi,
+  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>
-  static void CMS(PType & p1, PType & p2, Energy2 s);
+  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>
-  static void CMS(const PPairType & p, Energy2 s)
+  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>
-  static void CMS(PType & p1, PType & p2, PType & p3, Energy2 s,
+  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>
-  static void CMS(PType & p1, PType & p2, PType & p3, Energy2 s,
+  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$).
    */
-  static Energy getMagnitude(Energy2 s, Energy m1, Energy m2);
+  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.
    */
-  static Momentum3 polar3Vector(Energy p, double costheta, double phi)
+  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$).
    */
-  static vector<LorentzMomentum>
+  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>
-  static void CMSn(Container & particles, Energy m0);
+  void CMSn(Container & particles, Energy m0);
 
 };
 
 }
 
 #ifndef ThePEG_TEMPLATES_IN_CC_FILE
 #include "SimplePhaseSpace.tcc"
 #endif
 
 #endif /* ThePEG_SimplePhaseSpace_H */
diff --git a/Utilities/XSecStat.h b/Utilities/XSecStat.h
--- a/Utilities/XSecStat.h
+++ b/Utilities/XSecStat.h
@@ -1,330 +1,321 @@
 // -*- C++ -*-
 //
 // XSecStat.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_XSecStat_H
 #define THEPEG_XSecStat_H
 //
 // This is the declaration of the XSecStat class.
 //
 
 #include "ThePEG/Config/ThePEG.h"
 
 namespace ThePEG {
 
 /**
  * XSecStat is a concrete helper class used to collect statistics
  * about the cross section for a specific process or group of
  * processes. It contains an overestimated cross section and
  * information about the number of times the process has been used to
  * generate an event and how many times this event has been accepted.
  *
  * An object of this class must initially be given an overestimated
  * cross section in the constructor or with the maxXSec(CrossSection)
  * function. Each time the corresponding process is selected
  * (according to maxXSec()), the select(double) function should be
  * called giving the weight with which the event will be accepted as
  * argument. If the event is then accepted, the accept() function
  * should be called. If an event is later vetoed, the reject()
  * function should be called.
  * 
  */
 class XSecStat {
 
 public:
 
   /**
    * Enumerate the different weight classes
    */
   enum {
     plainWeights = 0,
     plainVetoedWeights,
     reweightedWeights,
     reweightedVetoedWeights
   };
 
   /** @name Standard constructors, destructor and assignment operator. */
   //@{
   /**
    * The default constructor.
    */
   XSecStat() 
     : theMaxXSec(ZERO), theAttempts(0), theAccepted(0), theVetoed(0),
-      theSumWeights (vector<double>(4,0.0)),
-      theSumWeights2(vector<double>(4,0.0)), theLastWeight(0.0) {}
+      theSumWeights (),
+      theSumWeights2(), theLastWeight(0.0) {}
 
   /**
    * Constructor taking the overestimated cross section, \a xsecmax,
    * as argument.
    */
   explicit XSecStat(CrossSection xsecmax) 
     : theMaxXSec(xsecmax), theAttempts(0), theAccepted(0), theVetoed(0),
-      theSumWeights (vector<double>(4,0.0)),
-      theSumWeights2(vector<double>(4,0.0)), theLastWeight(0.0) {}
+      theSumWeights (),
+      theSumWeights2(), theLastWeight(0.0) {}
 
   /**
    * The assignment operator.
    */
-  XSecStat & operator=(const XSecStat & x) {
-    theMaxXSec     = x.theMaxXSec;
-    theAttempts    = x.theAttempts;
-    theAccepted    = x.theAccepted;
-    theVetoed      = x.theVetoed;
-    theSumWeights  = x.theSumWeights;
-    theSumWeights2 = x.theSumWeights2;
-    theLastWeight  = x.theLastWeight;
-    return *this;
-  }
+  XSecStat & operator=(const XSecStat & x) = default;
 
   /**
    * Add the contents of another XSecStat.
    */
   XSecStat & operator+=(const XSecStat & x) {
     theAttempts    += x.theAttempts;
     theAccepted    += x.theAccepted;
     theVetoed      += x.theVetoed;
     for( unsigned int ix = 0; ix < 4; ++ix ) {
       theSumWeights [ix] += x.theSumWeights [ix];
       theSumWeights2[ix] += x.theSumWeights2[ix];
     }
     theLastWeight = 0.0;
     return *this;
   }
 
   /**
    * Reset the statistics.
    */
   void reset() {
     theAttempts = theAccepted = theVetoed = 0;
-    theSumWeights = theSumWeights2 = vector<double>(4,0.0);
+    theSumWeights = theSumWeights2 = {};
     theLastWeight = 0.0;
   }
 
   //@}
 
 public:
 
   /** @name Simple access functions */
   //@{
 
   /**
    * An event of the corresponding class has been accepted. The
    * select() method must have been called before.
    */
   void accept() { 
     theAccepted += 1;
   }
 
   /**
    * An event of the corresponding class has been attempted. It will
    * subsequently be accepted with the given \a weight.
    */
   void select(double weight) {
     theAttempts += 1;
     theSumWeights [reweightedWeights] +=     weight ;
     theSumWeights2[reweightedWeights] += sqr(weight);
     theSumWeights [plainWeights]      +=     weight ;
     theSumWeights2[plainWeights]      += sqr(weight);
     theLastWeight = weight;
   }
 
   /**
    * Reweight a selected and accepted event.
    */
   void reweight(double oldWeight, double newWeight) {
     theSumWeights [reweightedWeights] +=     newWeight  -     oldWeight ;
     theSumWeights2[reweightedWeights] += sqr(newWeight) - sqr(oldWeight);
   }
 
   /**
    * Reject the event which was last accepted with accept() or
    * selected with select(double). The \a weight should be set to the
    * value, \f$w\f$, used in the previous call to select(double),
    * except if the event has been accepted with the probability
    * \f$w\f$, in which case \a weight should be set to \f$sign(1,
    * w)\f$.
    */
   void reject(double weight = 1.0) {
     theSumWeights [reweightedVetoedWeights] +=            weight ;
     theSumWeights2[reweightedVetoedWeights] +=        sqr(weight);
     theSumWeights [plainVetoedWeights]      +=     theLastWeight ;
     theSumWeights2[plainVetoedWeights]      += sqr(theLastWeight);
     theVetoed += 1;
   }
 
   /**
    * The overestimated cross section.
    */
   CrossSection maxXSec() const { return theMaxXSec; }
 
   /**
    * The sum of the weights so far.
    */
   double sumWeights() const { 
     return theSumWeights[reweightedWeights] - theSumWeights[reweightedVetoedWeights];
   }
 
   /**
    * The sum of the squared weights so far.
    */
   double sumWeights2() const {
     return theSumWeights2[reweightedWeights] + theSumWeights2[reweightedVetoedWeights];
   }
 
   /**
    * The sum of the weights so far, excluding reweighting.
    */
   double sumWeightsNoReweight() const { 
     return theSumWeights[plainWeights] - theSumWeights[plainVetoedWeights];
   }
 
   /**
    * The sum of the squared weights so far, excluding reweighting.
    */
   double sumWeights2NoReweight() const { 
     return theSumWeights2[plainWeights] + theSumWeights2[plainVetoedWeights];
   }
 
   /**
    * The current estimate of the cross section for the corresponding
    * class of events. If no events have been generated, maxXSec() will
    * be returned.
    */
   CrossSection xSec(double att = 0) const {
     double n = (att == 0.0 ? attempts() : att);
     return n ? maxXSec()*sumWeights()/n : maxXSec();
   }
 
   /**
    * The current estimate of the error in the cross section for the
    * corresponding class of events. If no events have been generated,
    * maxXSec() will be returned.
    */
   CrossSection xSecErr(double att = 0) const {
     double n = (att == 0.0 ? attempts() : att);
     if ( n < 2 )
       return maxXSec();
     double sw = sumWeights(); double sw2 = sumWeights2();
     return
       maxXSec()*sqrt(abs(sw2/n-sqr(sw/n))/(n-1));
   }
 
   /**
    * The current estimate of the cross section for the corresponding
    * class of events, excluding reweighting. If no events have been
    * generated, maxXSec() will be returned.
    */
   CrossSection xSecNoReweight(double att = 0) const {
     double n = (att == 0.0 ? attempts() : att);
     return n ? maxXSec()*sumWeightsNoReweight()/n : maxXSec();
   }
 
   /**
    * The current estimate of the error in the cross section for the
    * corresponding class of events, excluding reweighting. If no
    * events have been generated, maxXSec() will be returned.
    */
   CrossSection xSecErrNoReweight(double att = 0) const {
     double n = (att == 0.0 ? attempts() : att);
     if ( n < 2 )
       return maxXSec();
     double sw = sumWeightsNoReweight(); 
     double sw2 = sumWeights2NoReweight();
     return
       maxXSec()*sqrt(abs(sw2/n-sqr(sw/n))/(n-1));
   }
 
   /**
    * Number of attempts so far.
    */
   double attempts() const { return theAttempts; }
 
   /**
    * Number of accepts so far.
    */
   double accepted() const { return theAccepted-theVetoed; }
 
   /**
    * Number of vetoes so far.
    */
   double vetoed() const { return theVetoed; }
 
   /**
    * Set the overestimated cross section.
    */
   void maxXSec(CrossSection x) { theMaxXSec = x; }
   //@}
 
 public:
 
   /** @name I/O functions */
   //@{
   /**
    * Output to a persistent stream.
    */
   void output(PersistentOStream & os) const;
 
   /**
    * Input from a persistent stream.
    */
   void input(PersistentIStream & is);
   //@}
 
 private:
 
   /**
    * The overestimated cross section.
    */
   CrossSection theMaxXSec;
 
   /**
    * Number of attempts so far.
    */
   double theAttempts;
 
   /**
    * Number of accepted events so far.
    */
   double theAccepted;
 
   /**
    * Number of events vetoed after being accepted
    */
   double theVetoed;
 
   /**
    * The sum of the weights so far.
    */
-  vector<double> theSumWeights;
+  array<double,4> theSumWeights;
 
   /**
    * The sum of the squared weights so far.
    */
-  vector<double> theSumWeights2;
+  array<double,4> theSumWeights2;
 
   /**
    * The last selected weight, ignoring reweighting.
    */
   double theLastWeight;
 
 };
 
 /** Ouptut an XSecStat to a persistent stream. */
 PersistentOStream & operator<<(PersistentOStream &, const XSecStat &);
 
 /** Input an XSecStat from a persistent stream. */
 PersistentIStream & operator>>(PersistentIStream &, XSecStat &);
 
 /** Add the contents of two XSecStat objects. */
 inline XSecStat operator+(const XSecStat & x1, const XSecStat & x2) {
   XSecStat x = x1;
   return x += x2;
 }
 
 }
 
 #endif /* THEPEG_XSecStat_H */
diff --git a/Vectors/SpinHalfLorentzRotation.cc b/Vectors/SpinHalfLorentzRotation.cc
--- a/Vectors/SpinHalfLorentzRotation.cc
+++ b/Vectors/SpinHalfLorentzRotation.cc
@@ -1,445 +1,412 @@
 // -*- C++ -*-
 //
 // SpinHalfLorentzRotation.cc 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.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the SpinHalfLorentzRotation class.
 //
 
 #include "SpinHalfLorentzRotation.h"
 
 using namespace ThePEG;
 
 // default constructor
 SpinHalfLorentzRotation::SpinHalfLorentzRotation() 
-  : _mx(4,vector<Complex>(4)) {
-  unsigned ix;
-  for(ix=0;ix<4;++ix) {
-    _mx[ix][ix]=1.0;
-  }
-}
+: _mx{{ {{1., 0., 0., 0.}},
+        {{0., 1., 0., 0.}},
+        {{0., 0., 1., 0.}},
+        {{0., 0., 0., 1.}} }}
+{}
 
 // constructor giving the components of a Lorentz boost
 SpinHalfLorentzRotation::
 SpinHalfLorentzRotation(double bx, double by, double bz, double gamma) 
-  : _mx(4,vector<Complex>(4)) 
+: _mx{}
 {
-  setBoost (bx,by,bz,gamma);
+  setBoost(bx,by,bz,gamma);
 }
 
 // constructor with boost vector
 SpinHalfLorentzRotation::
-SpinHalfLorentzRotation (const Boost & b, double gamma) 
-: _mx(4,vector<Complex>(4)) {
+SpinHalfLorentzRotation (const Boost & b, double gamma)
+: _mx{}
+{
   setBoost(b,gamma);
 }
 
 // protected set all elements constructor
 SpinHalfLorentzRotation::
 SpinHalfLorentzRotation(Complex c1c1,Complex c1c2,Complex c1c3,Complex c1c4,
 			Complex c2c1,Complex c2c2,Complex c2c3,Complex c2c4,
 			Complex c3c1,Complex c3c2,Complex c3c3,Complex c3c4,
 			Complex c4c1,Complex c4c2,Complex c4c3,Complex c4c4) 
-: _mx(4,vector<Complex>(4)) {
-  _mx[0][0]=c1c1;_mx[1][0]=c2c1;_mx[2][0]=c3c1;_mx[3][0]=c4c1;
-  _mx[0][1]=c1c2;_mx[1][1]=c2c2;_mx[2][1]=c3c2;_mx[3][1]=c4c2;
-  _mx[0][2]=c1c3;_mx[1][2]=c2c3;_mx[2][2]=c3c3;_mx[3][2]=c4c3;
-  _mx[0][3]=c1c4;_mx[1][3]=c2c4;_mx[2][3]=c3c4;_mx[3][3]=c4c4;
-}
+: _mx{{ {{c1c1, c1c2, c1c3, c1c4}},
+        {{c2c1, c2c2, c2c3, c2c4}},
+        {{c3c1, c3c2, c3c3, c3c4}},
+        {{c4c1, c4c2, c4c3, c4c4}} }}
+{}
 
 // check for identity matrix
 bool SpinHalfLorentzRotation::isIdentity() const {
-  return (_mx[0][0]==1.0&&_mx[0][1]==0.0&&_mx[0][2]==0.0&&_mx[0][3]==0.0&& 
-          _mx[1][0]==0.0&&_mx[1][1]==1.0&&_mx[1][2]==0.0&&_mx[1][3]==0.0&& 
-          _mx[2][0]==0.0&&_mx[2][1]==0.0&&_mx[2][2]==1.0&&_mx[2][3]==0.0&& 
-          _mx[3][0]==0.0&&_mx[3][1]==0.0&&_mx[3][2]==0.0&&_mx[3][3]==1.0);
+  return _mx == MatrixT{{ {{1., 0., 0., 0.}},
+                          {{0., 1., 0., 0.}},
+                          {{0., 0., 1., 0.}},
+                          {{0., 0., 0., 1.}} }};
 }
 
 // inverse ( inverse is gamma0 S dagger gamma 0)
 SpinHalfLorentzRotation SpinHalfLorentzRotation::inverse() const {
   return SpinHalfLorentzRotation
-    (conj(_mx[2][2]),conj(_mx[3][2]),conj(_mx[0][2]),conj(_mx[1][2]),
+    {conj(_mx[2][2]),conj(_mx[3][2]),conj(_mx[0][2]),conj(_mx[1][2]),
      conj(_mx[2][3]),conj(_mx[3][3]),conj(_mx[0][3]),conj(_mx[1][3]),
      conj(_mx[2][0]),conj(_mx[3][0]),conj(_mx[0][0]),conj(_mx[1][0]),
-     conj(_mx[2][1]),conj(_mx[3][1]),conj(_mx[0][1]),conj(_mx[1][1]));
+     conj(_mx[2][1]),conj(_mx[3][1]),conj(_mx[0][1]),conj(_mx[1][1])};
 }
 
 // specify the components of a lorentz boost
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::
 setBoost (double bx, double by, double bz, double) {
   // work out beta and chi
-  static double eps=1e-10;
+  const double eps=1e-10;
   double beta(sqrt(bx*bx+by*by+bz*bz));
   double chi(atanh(beta)), chc(cosh(0.5*chi)), shc(0.5);
   if ( beta > eps ) 
     shc=sinh(0.5*chi)/beta;
   Complex ii(0.,1.),nxminy(bx-ii*by),nxplny(bx+ii*by);
   _mx[0][0]= chc-shc*bz;
   _mx[0][1]=-shc*nxminy;
   _mx[0][2]= 0.        ;
   _mx[0][3]= 0.        ;
 
   _mx[1][0]=-shc*nxplny;
   _mx[1][1]= chc+shc*bz;
   _mx[1][2]= 0.        ;
   _mx[1][3]= 0.        ;
 
   _mx[2][0]= 0.        ;
   _mx[2][1]= 0.        ;
   _mx[2][2]= chc+shc*bz;
   _mx[2][3]=+shc*nxminy;
   
   _mx[3][0]= 0.        ;
   _mx[3][1]= 0.        ;
   _mx[3][2]=+shc*nxplny;
   _mx[3][3]= chc-shc*bz;
   return *this;
 }
 
 // specify a boost vector
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::setBoost (const Boost & b,double) {
   // work out beta and chi
-  static double eps=1e-10;
+  const double eps=1e-10;
   double bx(b.x()),by(b.y()),bz(b.z()),beta(b.mag()),chi(atanh(beta)),
     chc(cosh(0.5*chi)),shc(0.5);
   if(beta>eps){shc=sinh(0.5*chi)/beta;}
   Complex ii(0.,1.),nxminy(bx-ii*by),nxplny(bx+ii*by);
   _mx[0][0]= chc-shc*bz;_mx[0][1]=-shc*nxminy;_mx[0][2]= 0.        ;_mx[0][3]= 0.        ;
   _mx[1][0]=-shc*nxplny;_mx[1][1]= chc+shc*bz;_mx[1][2]= 0.        ;_mx[1][3]= 0.        ;
   _mx[2][0]= 0.        ;_mx[2][1]= 0.        ;_mx[2][2]= chc+shc*bz;_mx[2][3]=+shc*nxminy;
   _mx[3][0]= 0.        ;_mx[3][1]= 0.        ;_mx[3][2]=+shc*nxplny;_mx[3][3]= chc-shc*bz;
   return *this;
 }
 
 // lorentz boost in x direction
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::setBoostX (double & bx) {
   // work out beta and chi
   double chi(atanh(bx)),shc(sinh(0.5*chi)),chc(cosh(0.5*chi));
   _mx[0][0]= chc;_mx[0][1]=-shc;_mx[0][2]= 0. ;_mx[0][3]= 0. ;
   _mx[1][0]=-shc;_mx[1][1]= chc;_mx[1][2]= 0. ;_mx[1][3]= 0. ;
   _mx[2][0]= 0  ;_mx[2][1]= 0. ;_mx[2][2]= chc;_mx[2][3]=+shc;
   _mx[3][0]= 0  ;_mx[3][1]= 0. ;_mx[3][2]=+shc;_mx[3][3]= chc;
   return *this;
 }
 
 // lorentz boost in y direction
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::setBoostY (double & by)
 {
   // work out beta and chi
   double chi(atanh(by)),chc(cosh(0.5*chi));
   Complex shc(0.,sinh(0.5*chi));
   _mx[0][0]= chc;_mx[0][1]= shc;_mx[0][2]= 0. ;_mx[0][3]= 0. ;
   _mx[1][0]=-shc;_mx[1][1]= chc;_mx[1][2]= 0. ;_mx[1][3]= 0  ;
   _mx[2][0]= 0. ;_mx[2][1]= 0. ;_mx[2][2]= chc;_mx[2][3]=-shc;
   _mx[3][0]= 0. ;_mx[3][1]= 0. ;_mx[3][2]=+shc;_mx[3][3]= chc;
   return *this;
 }
 
 // lorentz boost in z direction 
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::setBoostZ (double & bz)
 {
   // work out beta and chi
   double chi(atanh(bz)),shc(sinh(0.5*chi)),chc(cosh(0.5*chi));
   _mx[0][0]= chc-shc;_mx[0][1]= 0.     ;_mx[0][2]= 0.        ;_mx[0][3]= 0.     ;
   _mx[1][0]= 0.     ;_mx[1][1]= chc+shc;_mx[1][2]= 0.        ;_mx[1][3]= 0.     ;
   _mx[2][0]= 0.     ;_mx[2][1]= 0.     ;_mx[2][2]= chc+shc   ;_mx[2][3]= 0.     ;
   _mx[3][0]= 0.     ;_mx[3][1]= 0.     ;_mx[3][2]=+0.        ;_mx[3][3]= chc-shc;
   return *this;
 }
 
 // Pure boost along the x-axis; equivalent to LT = BoostX(beta) * LT
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::boostX(double bx)
 {
   double chi(atanh(bx)),shc(sinh(0.5*chi)),chc(cosh(0.5*chi));
-  Complex temp[4][4];
-  unsigned int ix,iy;
-  for(ix=0;ix<4;++ix)
+  MatrixT temp;
+  for(size_t ix=0;ix<4;++ix)
     {
       temp[0][ix]= chc*_mx[0][ix]-shc*_mx[1][ix];
       temp[1][ix]= chc*_mx[1][ix]-shc*_mx[0][ix];
       temp[2][ix]= chc*_mx[2][ix]+shc*_mx[3][ix];
       temp[3][ix]= chc*_mx[3][ix]+shc*_mx[2][ix];
     }
-  for(ix=0;ix<4;++ix)
-    for(iy=0;iy<4;++iy)
-      _mx[ix][iy]=temp[ix][iy];
+  _mx = temp;
   return *this;
 }
 // Pure boost along the y-axis; equivalent to LT = BoostY(beta) * LT
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::boostY(double by)
 {
   double chi(atanh(by)),chc(cosh(0.5*chi));
   Complex shc(0.,sinh(0.5*chi));
-  Complex temp[4][4];
-  unsigned int ix,iy;
-  for(ix=0;ix<4;++ix)
+  MatrixT temp;
+  for(size_t ix=0;ix<4;++ix)
     {
       temp[0][ix]= chc*_mx[0][ix]+shc*_mx[1][ix];
       temp[1][ix]= chc*_mx[1][ix]-shc*_mx[0][ix];
       temp[2][ix]= chc*_mx[2][ix]-shc*_mx[3][ix];
       temp[3][ix]= chc*_mx[3][ix]+shc*_mx[2][ix];
     }
-  for(ix=0;ix<4;++ix)
-    for(iy=0;iy<4;++iy)
-      _mx[ix][iy]=temp[ix][iy];
+  _mx = temp;
   return *this;
 }
 
 // Pure boost along the z-axis; equivalent to LT = BoostX(beta) * LT
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::boostZ(double bz)
 {
   double chi(atanh(bz)),shc(sinh(0.5*chi)),chc(cosh(0.5*chi));
-  Complex temp[4][4];
-  unsigned int ix,iy;
-  for(ix=0;ix<4;++ix)
+  MatrixT temp;
+  for(size_t ix=0;ix<4;++ix)
     {
       temp[0][ix]=(chc-shc)*_mx[0][ix];
       temp[1][ix]=(chc+shc)*_mx[1][ix];
       temp[2][ix]=(chc+shc)*_mx[2][ix];
       temp[3][ix]=( chc-shc)*_mx[3][ix];
     }
-  for(ix=0;ix<4;++ix)
-    for(iy=0;iy<4;++iy)
-      _mx[ix][iy]=temp[ix][iy];
+  _mx = temp;
   return *this;
 }
 
 // General boost equivalent to LT = Boost(bx,by,bz) * LT
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::boost(double bx, double by, double bz, double gamma) {
   // calculation of gamma and beta
   double b2(bx*bx+by*by+bz*bz);
   if(gamma<1.) gamma = 1./sqrt(1.-b2);
   double beta = sqrt(b2);
   // work out beta and chi
-  static double eps=1e-8;
+  const double eps=1e-8;
   double chc = sqrt(0.5*(1+gamma));
   double shc = beta>eps ? sqrt(0.5*(gamma-1))/beta : 0.5+b2*(0.1875+0.12109375*b2);
-  Complex ii(0.,1.),nxminy(bx-ii*by),nxplny(bx+ii*by),temp[4][4];
-  unsigned int ix,iy;
-  for(ix=0;ix<4;++ix) {
+  Complex ii(0.,1.),nxminy(bx-ii*by),nxplny(bx+ii*by);
+  MatrixT temp;
+  for(size_t ix=0;ix<4;++ix) {
     temp[0][ix]= (chc-shc*bz)*_mx[0][ix]-shc*nxminy  *_mx[1][ix];
     temp[1][ix]=-shc*nxplny  *_mx[0][ix]+(chc+shc*bz)*_mx[1][ix];
     temp[2][ix]= (chc+shc*bz)*_mx[2][ix]+shc*nxminy  *_mx[3][ix];
     temp[3][ix]= shc*nxplny  *_mx[2][ix]+(chc-shc*bz)*_mx[3][ix];
   }
-  for(ix=0;ix<4;++ix)
-    for(iy=0;iy<4;++iy)
-      _mx[ix][iy]=temp[ix][iy];
+  _mx = temp;
   return *this;
 }
 
 // General boost equivalent to LT = Boost(bv) * LT
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::boost(const Boost & b, double gamma) {
   // calculation of gamma and beta
   double b2(b.mag2());
   if(gamma<1.) gamma = 1./sqrt(1.-b2);
   double beta = sqrt(b2);
   // work out chi
-  static double eps=1e-8;
+  const double eps=1e-8;
   double bx(b.x()),by(b.y()),bz(b.z());
   double chc = sqrt(0.5*(1+gamma));
   double shc = beta>eps ? sqrt(0.5*(gamma-1))/beta : 0.5+b2*(0.1875+0.12109375*b2);
-  Complex ii(0.,1.),nxminy(bx-ii*by),nxplny(bx+ii*by),temp[4][4];
-  unsigned int ix,iy;
-  for(ix=0;ix<4;++ix) {
+  Complex ii(0.,1.),nxminy(bx-ii*by),nxplny(bx+ii*by);
+  MatrixT temp;
+  for(size_t ix=0;ix<4;++ix) {
     temp[0][ix]= (chc-shc*bz)*_mx[0][ix]-shc*nxminy  *_mx[1][ix];
     temp[1][ix]=-shc*nxplny  *_mx[0][ix]+(chc+shc*bz)*_mx[1][ix];
     temp[2][ix]= (chc+shc*bz)*_mx[2][ix]+shc*nxminy  *_mx[3][ix];
     temp[3][ix]= shc*nxplny  *_mx[2][ix]+(chc-shc*bz)*_mx[3][ix];
   }
-  for(ix=0;ix<4;++ix)
-    for(iy=0;iy<4;++iy)
-      _mx[ix][iy]=temp[ix][iy];
+  _mx = temp;
   return *this;
 }
 
 std::ostream & SpinHalfLorentzRotation::print( std::ostream & os ) const {
   os << "\n   [ ( " <<
     std::setw(14) << std::setprecision(6) << s1s1() << "   " <<
     std::setw(14) << std::setprecision(6) << s1s2() << "   " <<
     std::setw(14) << std::setprecision(6) << s1s3() << "   " <<
     std::setw(14) << std::setprecision(6) << s1s4() << ")\n"
      << "     ( " <<
     std::setw(14) << std::setprecision(6) << s2s1() << "   " <<
     std::setw(14) << std::setprecision(6) << s2s2() << "   " <<
     std::setw(14) << std::setprecision(6) << s2s3() << "   " <<
     std::setw(14) << std::setprecision(6) << s2s4() << ")\n"
      << "     ( " <<
     std::setw(14) << std::setprecision(6) << s3s1() << "   " <<
     std::setw(14) << std::setprecision(6) << s3s2() << "   " <<
     std::setw(14) << std::setprecision(6) << s3s3() << "   " <<
     std::setw(14) << std::setprecision(6) << s3s4() << ")\n"
      << "     ( " <<
     std::setw(14) << std::setprecision(6) << s4s1() << "   " <<
     std::setw(14) << std::setprecision(6) << s4s2() << "   " <<
     std::setw(14) << std::setprecision(6) << s4s3() << "   " <<
     std::setw(14) << std::setprecision(6) << s4s4() << ") ]\n";
   return os;
 }
 
 
 // general rotation
 SpinHalfLorentzRotation & 
 SpinHalfLorentzRotation::setRotate(double phi, const Axis & axis) {
   double cp(cos(0.5*phi));
   // get the normalised components of the vector
   double amag(axis.mag()),ax(axis.x()/amag),ay(axis.y()/amag),az(axis.z()/amag);
   Complex ii(0.,1.),nxminy(ax-ii*ay),nxplny(ax+ii*ay),isp(0.,sin(0.5*phi));
   // rotatation matrix is the same in both conventions
   _mx[0][0]= cp-isp*az ;_mx[0][1]=-isp*nxminy;_mx[0][2]= 0.        ;_mx[0][3]= 0.        ;
   _mx[1][0]=-isp*nxplny;_mx[1][1]= cp+isp*az ;_mx[1][2]= 0.        ;_mx[1][3]= 0.        ;
   _mx[2][0]= 0.        ;_mx[2][1]= 0.        ;_mx[2][2]= cp-isp*az ;_mx[2][3]=-isp*nxminy;
   _mx[3][0]= 0.        ;_mx[3][1]= 0.        ;_mx[3][2]=-isp*nxplny;_mx[3][3]= cp+isp*az ;
   return *this;
 }
 
 // rotation about x
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::setRotateX(double& phi) {
   double cp(cos(0.5*phi));
   Complex isp(0.,sin(0.5*phi));
   // rotatation matrix is the same in both conventions
   _mx[0][0]= cp ;_mx[0][1]=-isp;_mx[0][2]= 0. ;_mx[0][3]= 0. ;
   _mx[1][0]=-isp;_mx[1][1]= cp ;_mx[1][2]= 0. ;_mx[1][3]= 0. ;
   _mx[2][0]= 0. ;_mx[2][1]= 0. ;_mx[2][2]= cp ;_mx[2][3]=-isp;
   _mx[3][0]= 0. ;_mx[3][1]= 0. ;_mx[3][2]=-isp;_mx[3][3]= cp ;
   return *this;
 }
 
 // rotation about y
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::setRotateY(double& phi) {
   double cp(cos(0.5*phi)),sp(sin(0.5*phi));
   // rotatation matrix is the same in both conventions
   _mx[0][0]= cp;_mx[0][1]=-sp;_mx[0][2]= 0.;_mx[0][3]= 0.;
   _mx[1][0]= sp;_mx[1][1]= cp;_mx[1][2]= 0.;_mx[1][3]= 0.;
   _mx[2][0]= 0.;_mx[2][1]= 0.;_mx[2][2]= cp;_mx[2][3]=-sp;
   _mx[3][0]= 0.;_mx[3][1]= 0.;_mx[3][2]= sp;_mx[3][3]= cp;
   return *this;
 }
 
 // rotation about z
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::setRotateZ(double& phi)
 {
   double cp(cos(0.5*phi));
   Complex isp(0.,sin(0.5*phi));
   // rotatation matrix is the same in both conventions
   _mx[0][0]= cp-isp ;_mx[0][1]= 0.    ;_mx[0][2]= 0.    ;_mx[0][3]= 0.    ;
   _mx[1][0]= 0.     ;_mx[1][1]= cp+isp;_mx[1][2]= 0.    ;_mx[1][3]= 0.    ;
   _mx[2][0]= 0.     ;_mx[2][1]= 0.    ;_mx[2][2]= cp-isp;_mx[2][3]= 0.    ;
   _mx[3][0]= 0.     ;_mx[3][1]= 0.    ;_mx[3][2]= 0.    ;_mx[3][3]= cp+isp;
   return *this;
 }
 
 
 // product
 SpinHalfLorentzRotation 
 SpinHalfLorentzRotation::operator * (const SpinHalfLorentzRotation & lt) const {
-  Complex output[4][4];
-  unsigned int ix,iy,iz;
-  for(ix=0;ix<4;++ix) {
-    for(iy=0;iy<4;++iy) {
-      output[ix][iy]=0.;
-      for(iz=0;iz<4;++iz){output[ix][iy]+=_mx[ix][iz]*lt._mx[iz][iy];}
-    }
-  }
-  return SpinHalfLorentzRotation(output[0][0],output[0][1],output[0][2],output[0][3],
-				 output[1][0],output[1][1],output[1][2],output[1][3],
-				 output[2][0],output[2][1],output[2][2],output[2][3],
-				 output[3][0],output[3][1],output[3][2],output[3][3]);
+  SpinHalfLorentzRotation temp{_mx};
+  temp *= lt;
+  return temp;
 }
 
 // multiply and assign
- SpinHalfLorentzRotation & 
+SpinHalfLorentzRotation & 
 SpinHalfLorentzRotation::operator *= (const SpinHalfLorentzRotation & lt) {
-  Complex output[4][4];
-  unsigned int ix,iy,iz;
-  for(ix=0;ix<4;++ix) {
-      for(iy=0;iy<4;++iy) {
-	output[ix][iy]=0.;
-	for(iz=0;iz<4;++iz){output[ix][iy]+=_mx[ix][iz]*lt._mx[iz][iy];}
-      }
-  }
-  for(ix=0;ix<4;++ix){for(iy=0;iy<4;++iy){_mx[ix][iy]=output[ix][iy];}}
+  MatrixT temp;
+  for(size_t ix=0;ix<4;++ix)
+    for(size_t iy=0;iy<4;++iy)
+      for(size_t iz=0;iz<4;++iz)
+        temp[ix][iy] += _mx[ix][iz] * lt._mx[iz][iy];
+  _mx = temp;
   return *this;
 }
 
 // transform method
- SpinHalfLorentzRotation & 
+SpinHalfLorentzRotation & 
 SpinHalfLorentzRotation::transform(const SpinHalfLorentzRotation & lt) {
-  Complex output[4][4];
-  unsigned int ix,iy,iz;
-  for(ix=0;ix<4;++ix) {
-    for(iy=0;iy<4;++iy) {
-      output[ix][iy]=0.;
-      for(iz=0;iz<4;++iz){output[ix][iy]+=lt._mx[ix][iz]*_mx[iz][iy];}
-    }
-  }
-  for(ix=0;ix<4;++ix){for(iy=0;iy<4;++iy){_mx[ix][iy]=output[ix][iy];}}
+  SpinHalfLorentzRotation temp{lt._mx};
+  temp *= (*this);
+  _mx = temp._mx;
   return *this;
 }
 
 // Rotation around the x-axis; equivalent to LT = RotationX(delta) * LT
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::rotateX(double phi) {
   double cp(cos(0.5*phi));
-  Complex isp(0.,sin(0.5*phi)),temp[4][4];
-  unsigned int ix,iy;
-  for(ix=0;ix<4;++ix) {
+  Complex isp(0.,sin(0.5*phi));
+  MatrixT temp;
+  for(size_t ix=0;ix<4;++ix) {
     temp[0][ix]=  cp*_mx[0][ix]-isp*_mx[1][ix];
     temp[1][ix]=-isp*_mx[0][ix]+ cp*_mx[1][ix];
     temp[2][ix]=  cp*_mx[2][ix]-isp*_mx[3][ix];
     temp[3][ix]=-isp*_mx[2][ix]+ cp*_mx[3][ix];
   }
-  for(ix=0;ix<4;++ix){for(iy=0;iy<4;++iy){_mx[ix][iy]=temp[ix][iy];}}
+  _mx = temp;
   return *this;
 }
 
 // Rotation around the y-axis; equivalent to LT = RotationY(delta) * LT
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::rotateY(double phi) {
   double cp(cos(0.5*phi)),sp(sin(0.5*phi));
-  Complex temp[4][4];
-  unsigned int ix,iy;
-  for(ix=0;ix<4;++ix) {
+  MatrixT temp;
+  for(size_t ix=0;ix<4;++ix) {
     temp[0][ix]= cp*_mx[0][ix]-sp*_mx[1][ix];
     temp[1][ix]= sp*_mx[0][ix]+cp*_mx[1][ix];
     temp[2][ix]= cp*_mx[2][ix]-sp*_mx[3][ix];
     temp[3][ix]= sp*_mx[2][ix]+cp*_mx[3][ix];
   }
-  for(ix=0;ix<4;++ix){for(iy=0;iy<4;++iy){_mx[ix][iy]=temp[ix][iy];}}
+  _mx = temp;
   return *this;
 }
 
 // Rotation around the z-axis; equivalent to LT = RotationZ(delta) * LT
 SpinHalfLorentzRotation & SpinHalfLorentzRotation::rotateZ(double phi) {
   double cp(cos(0.5*phi));
-  Complex isp(0.,sin(0.5*phi)),temp[4][4];
-  unsigned int ix,iy;
-  for(ix=0;ix<4;++ix) {
+  Complex isp(0.,sin(0.5*phi));
+  MatrixT temp;
+  for(size_t ix=0;ix<4;++ix) {
     temp[0][ix]= (cp-isp)*_mx[0][ix];
     temp[1][ix]= (cp+isp)*_mx[1][ix];
     temp[2][ix]= (cp-isp)*_mx[2][ix];
     temp[3][ix]= (cp+isp)*_mx[3][ix];
   }
-  for(ix=0;ix<4;++ix){for(iy=0;iy<4;++iy){_mx[ix][iy]=temp[ix][iy];}}
+  _mx = temp;
   return *this;
 }
 
 // Rotation around specified vector - LT = Rotation(delta,axis)*LT
 SpinHalfLorentzRotation & 
 SpinHalfLorentzRotation::rotate(double phi, const Axis & axis) {
   double cp(cos(0.5*phi)),amag(axis.mag()),
     ax(axis.x()/amag),ay(axis.y()/amag),az(axis.z()/amag);
-  Complex ii(0.,1.),nxminy(ax-ii*ay),nxplny(ax+ii*ay),isp(0.,sin(0.5*phi)),temp[4][4];
-  unsigned int ix,iy;
-  for(ix=0;ix<4;++ix) {
+  Complex ii(0.,1.),nxminy(ax-ii*ay),nxplny(ax+ii*ay),isp(0.,sin(0.5*phi));
+  MatrixT temp;
+  for(size_t ix=0;ix<4;++ix) {
     temp[0][ix]= (cp-isp*az)*_mx[0][ix]-isp*nxminy *_mx[1][ix];
     temp[1][ix]=-isp*nxplny *_mx[0][ix]+(cp+isp*az)*_mx[1][ix];
     temp[2][ix]= (cp-isp*az)*_mx[2][ix]-isp*nxminy *_mx[3][ix];
     temp[3][ix]=-isp*nxplny *_mx[2][ix]+(cp+isp*az)*_mx[3][ix];
   }
-  for(ix=0;ix<4;++ix){for(iy=0;iy<4;++iy){_mx[ix][iy]=temp[ix][iy];}}
+  _mx = temp;
   return *this;
 }
diff --git a/Vectors/SpinHalfLorentzRotation.h b/Vectors/SpinHalfLorentzRotation.h
--- a/Vectors/SpinHalfLorentzRotation.h
+++ b/Vectors/SpinHalfLorentzRotation.h
@@ -1,350 +1,353 @@
 // -*- C++ -*-
 //
 // SpinHalfLorentzRotation.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_SpinHalfLorentzRotation_H
 #define THEPEG_SpinHalfLorentzRotation_H
 //
 // This is the declaration of the SpinHalfLorentzRotation class.
 //
 #include "ThePEG/Helicity/HelicityDefinitions.h"
 #include "ThreeVector.h"
 
 namespace ThePEG {
 
 /**
  * The SpinHalfLorentzRotation class is designed to offer the same
  * features as the HepLorentzRotation class of CLHEP but for the spin-\f$\frac12\f$
  * Lorentz transformation. This is then combined into the general LorentzRotation
  * class of ThePEG to provide the Lorentz transformation for any object as the
  * transformations for higher spin objects can be built from the spin-\f$\frac12\f$
  * and spin-1 transformations.
  *
  * The boost matrix is calculated using the default Dirac matrix representation.
  * Any conversion to other Dirac matrix representations must be handled when the
  * transformation is used.
  */
 class SpinHalfLorentzRotation {
 
   /**
    * The external inverseOf needs to be a friend
    */
   friend SpinHalfLorentzRotation inverseOf ( const SpinHalfLorentzRotation & lt );
 
 public:
 
   /** @name Constructors and destructor. */
   //@{
 
   /**
    * Default constructor. Gives a unit matrix.
    */
   SpinHalfLorentzRotation();
 
   /**
    * Constructor giving the components of a Lorentz boost.
    * @param bx The x-component of the boost
    * @param by The y-component of the boost
    * @param bz The z-component of the boost
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinHalfLorentzRotation (double bx, double by, double bz, double gamma=-1.);
 
   /**
    * Constructor giving the vector for a Lorentz boost.
    * @param b The boost vector
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinHalfLorentzRotation (const Boost & b,double gamma=-1.);
   //@}
 
   /**
    * Returns true if the Identity matrix.
    */
   bool isIdentity() const;
 
   /**
    * Return the inverse.
    */
   SpinHalfLorentzRotation inverse() const;
 
   /**
    * Inverts the SpinHalfLorentzRotation matrix.
    */
   SpinHalfLorentzRotation & invert() { return *this = inverse(); }
 
   /**
    *  output operator
    */
   std::ostream & print( std::ostream & os ) const;
 
   /** @name Set methods for speical cases of simple rotations and boosts */
   //@{
 
   /**
    * Specify the components of a Lorentz Boost
    * @param bx The x-component of the boost
    * @param by The y-component of the boost
    * @param bz The z-component of the boost
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinHalfLorentzRotation & setBoost (double bx, double by, double bz,double gamma=-1.);
 
   /**
    * Specify a Lorentz Boost as a vector
    * @param b The boost vector
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinHalfLorentzRotation & setBoost (const Boost & b,double gamma=-1.);
 
   /**
    * Specify a boost by the given factor along the x-axis
    * @param boost The Lorentz boost  
    */
   SpinHalfLorentzRotation & setBoostX (double & boost);
 
   /**
    * Specify a boost by the given factor along the y-axis
    * @param boost The Lorentz boost  
    */
   SpinHalfLorentzRotation & setBoostY (double & boost);
 
   /**
    * Specify a boost by the given factor along the z-axis
    * @param boost The Lorentz boost  
    */
   SpinHalfLorentzRotation & setBoostZ (double & boost);
 
   /**
    * Specify a rotation about a general axis by the angle given.
    * @param delta The angle
    * @param axis The axis
    */
   SpinHalfLorentzRotation & setRotate(double delta, const Axis & axis);
 
   /**
    * Specify a rotation by the given angle about the x-axis
    * @param angle The rotation angle 
    */
   SpinHalfLorentzRotation & setRotateX (double & angle);
 
   /**
    * Specify a rotation by the given angle about the y-axis
    * @param angle The rotation angle 
    */
   SpinHalfLorentzRotation & setRotateY (double & angle);
 
   /**
    * Specify a rotation by the given angle about the z-axis
    * @param angle The rotation angle 
    */
   SpinHalfLorentzRotation & setRotateZ (double & angle);
   
   //@}
 
 
   /** @name Access methods for the components */
   //@{
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s1s1() const { return _mx[0][0]; }
 
   /**
    *   The \f$(1,2)\f$ component
    */
   Complex s1s2() const { return _mx[0][1]; }
 
   /**
    *   The \f$(1,3)\f$ component
    */
   Complex s1s3() const { return _mx[0][2]; }
 
   /**
    *   The \f$(1,4)\f$ component
    */
   Complex s1s4() const { return _mx[0][3]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s2s1() const { return _mx[1][0]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s2s2() const { return _mx[1][1]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s2s3() const { return _mx[1][2]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s2s4() const { return _mx[1][3]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s3s1() const { return _mx[2][0]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s3s2() const { return _mx[2][1]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s3s3() const { return _mx[2][2]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s3s4() const { return _mx[2][3]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s4s1() const { return _mx[3][0]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s4s2() const { return _mx[3][1]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s4s3() const { return _mx[3][2]; }
 
   /**
    *   The \f$(1,1)\f$ component
    */
   Complex s4s4() const { return _mx[3][3]; }
 
   /**
    *  Fortran style subscript operator
    */
   Complex operator()(unsigned int i, unsigned int j) const {
     assert(i<=3 && j<=3);
     return _mx[i][j];
   }
   //@}
 
 
   /** @name Transformation and product members */
   //@{
 
   /**
    * Product of two SpinHalfLorentzRotations (this) * lt - matrix multiplication  
    * @param lt The SpinHalfLorentzRotation we are multiplying
    */
   SpinHalfLorentzRotation operator * (const SpinHalfLorentzRotation & lt) const;
 
   /**
    * Multiply by and assign a*=b becomes a= a*b
    */
    SpinHalfLorentzRotation & operator *= (const SpinHalfLorentzRotation & );
 
   /**
    *  Transform  (similar to *= but a.transform(b) becomes a = b*a
    */
    SpinHalfLorentzRotation & transform   (const SpinHalfLorentzRotation & );
 
   /**
    * Rotation around the x-axis; equivalent to LT = RotationX(delta) * LT
    */
   SpinHalfLorentzRotation & rotateX(double delta);
 
   /**
    * Rotation around the y-axis; equivalent to LT = RotationY(delta) * LT
    */
   SpinHalfLorentzRotation & rotateY(double delta);
 
   /**
    * Rotation around the z-axis; equivalent to LT = RotationZ(delta) * LT
    */
   SpinHalfLorentzRotation & rotateZ(double delta);
   
   /**
    *  Rotation around specified vector - LT = Rotation(delta,axis)*LT
    */
   SpinHalfLorentzRotation & rotate(double delta, const Axis & axis);
 
   /**
    * Pure boost along the x-axis; equivalent to LT = BoostX(beta) * LT
    */
   SpinHalfLorentzRotation & boostX(double beta);
 
   /**
    * Pure boost along the y-axis; equivalent to LT = BoostX(beta) * LT
    */
   SpinHalfLorentzRotation & boostY(double beta);
 
   /**
    * Pure boost along the z-axis; equivalent to LT = BoostX(beta) * LT
    */
   SpinHalfLorentzRotation & boostZ(double beta);
 
   /**
    * General boost equivalent to LT = Boost(bx,by,bz) * LT
    * @param bx The x-component of the boost
    * @param by The y-component of the boost
    * @param bz The z-component of the boost
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinHalfLorentzRotation & boost(double bx, double by, double bz, double gamma=-1.);
 
   /**
    * General boost equivalent to LT = Boost(bv) * LT
    * @param bv The boost vector
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinHalfLorentzRotation & boost(const Boost & bv, double gamma=-1.);
   //@}
 
 protected:
 
   /**
    *  Protected constructor giving all the members, no check it is a valid
    *  transformation
    */
   SpinHalfLorentzRotation(Complex s1s1,Complex s1s2,Complex s1s3,Complex s1s4,
 				 Complex s2s1,Complex s2s2,Complex s2s3,Complex s2s4,
 				 Complex s3s1,Complex s3s2,Complex s3s3,Complex s3s4,
 				 Complex s4s1,Complex s4s2,Complex s4s3,Complex s4s4);
 
 private:
 
+  using MatrixT = array<array<Complex,4>,4>;
+
+  SpinHalfLorentzRotation(const MatrixT & m) : _mx(m) {}
 
   /**
    * The members of the transformation matrix.
    */
-  vector<vector<Complex> > _mx;
+  MatrixT _mx;
 };
 
 /**
  *  Global method to get the inverse
  */
 inline SpinHalfLorentzRotation inverseOf ( const SpinHalfLorentzRotation & lt ) {
   return lt.inverse();
 }
 
 /**
  *  output operator
  */
 inline std::ostream & operator<< ( std::ostream & os,
 				   const  SpinHalfLorentzRotation& lt ) {
   return lt.print(os);
 }
 
 }
 
 #endif /* THEPEG_SpinHalfLorentzRotation_H */
diff --git a/Vectors/SpinOneLorentzRotation.cc b/Vectors/SpinOneLorentzRotation.cc
--- a/Vectors/SpinOneLorentzRotation.cc
+++ b/Vectors/SpinOneLorentzRotation.cc
@@ -1,188 +1,188 @@
 // -*- C++ -*-
 //
 // SpinOneLorentzRotation.cc 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.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the LorentzRotation class.
 //
 
 #include "SpinOneLorentzRotation.h"
 
 using namespace ThePEG;
 SpinOneLorentzRotation::
 SpinOneLorentzRotation(double xx, double xy, double xz, double xt,
 		       double yx, double yy, double yz, double yt,
 		       double zx, double zy, double zz, double zt,
 		       double tx, double ty, double tz, double tt)
-  : matrix_(16) {
+{
   xx_() = xx; xy_() = xy; xz_() = xz; xt_() = xt;
   yx_() = yx; yy_() = yy; yz_() = yz; yt_() = yt;
   zx_() = zx; zy_() = zy; zz_() = zz; zt_() = zt;
   tx_() = tx; ty_() = ty; tz_() = tz; tt_() = tt;
 }
 
 bool SpinOneLorentzRotation::isIdentity() const {
   return 
        1.0 == xx() && 0.0 == xy() && 0.0 == xz() && 0.0 == xt()
     && 0.0 == yx() && 1.0 == yy() && 0.0 == yz() && 0.0 == yt()
     && 0.0 == zx() && 0.0 == zy() && 1.0 == zz() && 0.0 == zt()
     && 0.0 == tx() && 0.0 == ty() && 0.0 == tz() && 1.0 == tt();
 }
 
 SpinOneLorentzRotation SpinOneLorentzRotation::inverse() const {
   return SpinOneLorentzRotation( xx(), yx(), zx(),-tx(),
 				 xy(), yy(), zy(),-ty(), 
 				 xz(), yz(), zz(),-tz(),
 				-xt(),-yt(),-zt(), tt());
 }
 
 // output operator
 std::ostream & SpinOneLorentzRotation::print( std::ostream &  os) const {  
   os << "\n   [ ( " <<
     std::setw(14) << std::setprecision(6) << xx() << "   " <<
     std::setw(14) << std::setprecision(6) << xy() << "   " <<
     std::setw(14) << std::setprecision(6) << xz() << "   " <<
     std::setw(14) << std::setprecision(6) << xt() << ")\n"
      << "     ( " <<
     std::setw(14) << std::setprecision(6) << yx() << "   " <<
     std::setw(14) << std::setprecision(6) << yy() << "   " <<
     std::setw(14) << std::setprecision(6) << yz() << "   " <<
     std::setw(14) << std::setprecision(6) << yt() << ")\n"
      << "     ( " <<
     std::setw(14) << std::setprecision(6) << zx() << "   " <<
     std::setw(14) << std::setprecision(6) << zy() << "   " <<
     std::setw(14) << std::setprecision(6) << zz() << "   " <<
     std::setw(14) << std::setprecision(6) << zt() << ")\n"
      << "     ( " <<
     std::setw(14) << std::setprecision(6) << tx() << "   " <<
     std::setw(14) << std::setprecision(6) << ty() << "   " <<
     std::setw(14) << std::setprecision(6) << tz() << "   " <<
     std::setw(14) << std::setprecision(6) << tt() << ") ]\n";
   return os;
 }
 
 SpinOneLorentzRotation & SpinOneLorentzRotation::
 setBoost (double bx, double by, double bz, double gamma) {
   double beta2 = sqr(bx) + sqr(by) + sqr(bz);
   if (beta2 >= 1.0) {
     throw Exception() << "Invalid boost (" << bx << ',' << by << ',' << bz 
 		      << ") in SpinOneLorentzRotatio::setBoost" << Exception::eventerror;
   }
   if(gamma<0.) gamma = 1.0 / sqrt((1.-bz)*(1.+bz)-sqr(bx)-sqr(by));
   double bgamma = sqr(gamma) / (1.0 + gamma);
   xx_() = 1.0 + bgamma * bx * bx;
   yy_() = 1.0 + bgamma * by * by;
   zz_() = 1.0 + bgamma * bz * bz;
   xy_() = yx_() = bgamma * bx * by;
   xz_() = zx_() = bgamma * bx * bz;
   yz_() = zy_() = bgamma * by * bz;
   xt_() = tx_() = gamma * bx;
   yt_() = ty_() = gamma * by;
   zt_() = tz_() = gamma * bz;
   tt_() = gamma;
   return *this;
 }
 
 SpinOneLorentzRotation & 
 SpinOneLorentzRotation::setRotate(double delta, const Axis & axis) {
 
   double sinDelta = sin(delta), cosDelta = cos(delta);
   double oneMinusCosDelta = 1.0 - cosDelta;
 
   Axis u = unitVector(axis);
 
   double uX = u.x();
   double uY = u.y();
   double uZ = u.z();
 
   double rxx = oneMinusCosDelta * uX * uX  +  cosDelta;
   double rxy = oneMinusCosDelta * uX * uY  -  sinDelta * uZ;
   double rxz = oneMinusCosDelta * uX * uZ  +  sinDelta * uY;
 
   double ryx = oneMinusCosDelta * uY * uX  +  sinDelta * uZ;
   double ryy = oneMinusCosDelta * uY * uY  +  cosDelta;
   double ryz = oneMinusCosDelta * uY * uZ  -  sinDelta * uX;
 
   double rzx = oneMinusCosDelta * uZ * uX  -  sinDelta * uY;
   double rzy = oneMinusCosDelta * uZ * uY  +  sinDelta * uX;
   double rzz = oneMinusCosDelta * uZ * uZ  +  cosDelta;
 
   xx_() = rxx; xy_() = rxy; xz_() = rxz; xt_() = 0.0;
   yx_() = ryx; yy_() = ryy; yz_() = ryz; yt_() = 0.0;
   zx_() = rzx; zy_() = rzy; zz_() = rzz; zt_() = 0.0;
   tx_() = 0.0; ty_() = 0.0; tz_() = 0.0; tt_() = 1.0;
 
   return  *this;
 }
 
 SpinOneLorentzRotation & 
 SpinOneLorentzRotation::setRotateX(double delta) {
   double sinDelta = sin(delta), cosDelta = cos(delta);
 
   double ryy = cosDelta, ryz = -sinDelta;
   double rzy = sinDelta, rzz =  cosDelta;
 
   xx_() = 1.0; xy_() = 0.0; xz_() = 0.0; xt_() = 0.0;
   yx_() = 0.0; yy_() = ryy; yz_() = ryz; yt_() = 0.0;
   zx_() = 0.0; zy_() = rzy; zz_() = rzz; zt_() = 0.0;
   tx_() = 0.0; ty_() = 0.0; tz_() = 0.0; tt_() = 1.0;
   return  *this;
 }
 
 SpinOneLorentzRotation & 
 SpinOneLorentzRotation::setRotateY(double delta) {
   double sinDelta = sin(delta), cosDelta = cos(delta);
 
   double rxx =  cosDelta, rxz =  sinDelta;
   double rzx = -sinDelta, rzz =  cosDelta;
 
   xx_() = rxx; xy_() = 0.0; xz_() = rxz; xt_() = 0.0;
   yx_() = 0.0; yy_() = 1.0; yz_() = 0.0; yt_() = 0.0;
   zx_() = rzx; zy_() = 0.0; zz_() = rzz; zt_() = 0.0;
   tx_() = 0.0; ty_() = 0.0; tz_() = 0.0; tt_() = 1.0;
   return  *this;
 }
 
 SpinOneLorentzRotation & 
 SpinOneLorentzRotation::setRotateZ(double delta) {
   double sinDelta = sin(delta), cosDelta = cos(delta);
 
   double rxx = cosDelta, rxy = -sinDelta;
   double ryx = sinDelta, ryy =  cosDelta;
 
   xx_() = rxx; xy_() = rxy; xz_() = 0.0; xt_() = 0.0;
   yx_() = ryx; yy_() = ryy; yz_() = 0.0; yt_() = 0.0;
   zx_() = 0.0; zy_() = 0.0; zz_() = 1.0; zt_() = 0.0;
   tx_() = 0.0; ty_() = 0.0; tz_() = 0.0; tt_() = 1.0;
   return  *this;
 }
 
 SpinOneLorentzRotation 
 SpinOneLorentzRotation::operator*(const SpinOneLorentzRotation & b) const {
   return SpinOneLorentzRotation
     (xx()*b.xx() + xy()*b.yx() + xz()*b.zx() + xt()*b.tx(),
      xx()*b.xy() + xy()*b.yy() + xz()*b.zy() + xt()*b.ty(),
      xx()*b.xz() + xy()*b.yz() + xz()*b.zz() + xt()*b.tz(),
      xx()*b.xt() + xy()*b.yt() + xz()*b.zt() + xt()*b.tt(),
      
      yx()*b.xx() + yy()*b.yx() + yz()*b.zx() + yt()*b.tx(),
      yx()*b.xy() + yy()*b.yy() + yz()*b.zy() + yt()*b.ty(),
      yx()*b.xz() + yy()*b.yz() + yz()*b.zz() + yt()*b.tz(),
      yx()*b.xt() + yy()*b.yt() + yz()*b.zt() + yt()*b.tt(),
      
      zx()*b.xx() + zy()*b.yx() + zz()*b.zx() + zt()*b.tx(),
      zx()*b.xy() + zy()*b.yy() + zz()*b.zy() + zt()*b.ty(),
      zx()*b.xz() + zy()*b.yz() + zz()*b.zz() + zt()*b.tz(),
      zx()*b.xt() + zy()*b.yt() + zz()*b.zt() + zt()*b.tt(),
 
      tx()*b.xx() + ty()*b.yx() + tz()*b.zx() + tt()*b.tx(),
      tx()*b.xy() + ty()*b.yy() + tz()*b.zy() + tt()*b.ty(),
      tx()*b.xz() + ty()*b.yz() + tz()*b.zz() + tt()*b.tz(),		
      tx()*b.xt() + ty()*b.yt() + tz()*b.zt() + tt()*b.tt());
 }
 
diff --git a/Vectors/SpinOneLorentzRotation.h b/Vectors/SpinOneLorentzRotation.h
--- a/Vectors/SpinOneLorentzRotation.h
+++ b/Vectors/SpinOneLorentzRotation.h
@@ -1,394 +1,394 @@
 // -*- C++ -*-
 //
 // SpinOneLorentzRotation.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_SpinOneLorentzRotation_H
 #define ThePEG_SpinOneLorentzRotation_H
 
 #include "ThePEG/Helicity/HelicityDefinitions.h"
 #include "ThePEG/Helicity/LorentzTensor.fh"
 #include "ThePEG/Helicity/LorentzRSSpinor.fh"
 #include "ThePEG/Helicity/LorentzRSSpinorBar.fh"
 #include "ThreeVector.h"
 #include <vector>
 
 namespace ThePEG {
 
 /**
  * The SpinOneLorentzRotation class is ... */
 
 class SpinOneLorentzRotation {
 public:
 
   /** @name Constructors and destructor. */
   //@{
 
   /**
    * Default constructor. Gives a unit matrix.
    */
-  SpinOneLorentzRotation() : matrix_(16) {
+  SpinOneLorentzRotation() {
     xx_() = yy_() = zz_() = tt_() = 1.0;
   }
   
   /**
    * Constructor giving the components of a Lorentz boost.
    * @param bx The x-component of the boost
    * @param by The y-component of the boost
    * @param bz The z-component of the boost
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinOneLorentzRotation (double bx, double by, double bz, double gamma=-1.) 
-    : matrix_(16) {
+  {
     setBoost(bx,by,bz,gamma);
   }
 
   /**
    * Constructor giving the vector for a Lorentz boost.
    * @param b The boost vector
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   explicit SpinOneLorentzRotation (const Boost & b, double gamma=-1.)
-    : matrix_(16) {
+  {
     setBoost(b.x(), b.y(), b.z(),gamma);
   }
   //@}
 
   /**
    * Returns true if the Identity matrix.
    */
   bool isIdentity() const;
 
   /**
    * Return the inverse.
    */
   SpinOneLorentzRotation inverse() const;
 
   /**
    * Inverts the SpinOneLorentzRotation matrix.
    */
   SpinOneLorentzRotation & invert() { return *this = inverse(); }
 
   /**
    *  output operator
    */
   std::ostream & print( std::ostream & os ) const;
 
   /** @name Set methods for speical cases of simple rotations and boosts */
   //@{
 
   /**
    * Specify the components of a Lorentz Boost
    * @param bx The x-component of the boost
    * @param by The y-component of the boost
    * @param bz The z-component of the boost
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinOneLorentzRotation & setBoost (double bx, double by, double bz, double gamma=-1.);
 
   /**
    * Specify a Lorentz Boost as a vector
    * @param b The boost vector
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinOneLorentzRotation & setBoost (const Boost & b, double gamma=-1.) {
     return setBoost(b.x(), b.y(), b.z(),gamma); 
   }
 
   /**
    * Specify a rotation about a general axis by the angle given.
    * @param delta The angle
    * @param axis The axis
    */
   SpinOneLorentzRotation & setRotate(double delta, const Axis & axis);
 
   /**
    * Specify a rotation by the given angle about the x-axis
    * @param angle The rotation angle 
    */
   SpinOneLorentzRotation & setRotateX (double angle);
 
   /**
    * Specify a rotation by the given angle about the y-axis
    * @param angle The rotation angle 
    */
   SpinOneLorentzRotation & setRotateY (double angle);
 
   /**
    * Specify a rotation by the given angle about the z-axis
    * @param angle The rotation angle 
    */
   SpinOneLorentzRotation & setRotateZ (double angle);
   
   //@}
 
   /** @name Access methods for the components of the spin-1 rotation */
   //@{
 
   /**
    *   The xx component
    */
   double xx() const { return matrix_[ 0]; }
 
   /**
    *   The xy component
    */
   double xy() const { return matrix_[ 1]; }
 
   /**
    *   The xz component
    */
   double xz() const { return matrix_[ 2]; }
 
   /**
    *   The xt component
    */
   double xt() const { return matrix_[ 3]; }
 
   /**
    *   The yx component
    */
   double yx() const { return matrix_[ 4]; }
 
   /**
    *   The yy component
    */
   double yy() const { return matrix_[ 5]; }
 
   /**
    *   The yz component
    */
   double yz() const { return matrix_[ 6]; }
 
   /**
    *   The yt component
    */
   double yt() const { return matrix_[ 7]; }
 
   /**
    *   The zx component
    */
   double zx() const { return matrix_[ 8]; }
 
   /**
    *   The zy component
    */
   double zy() const { return matrix_[ 9]; }
 
   /**
    *   The zz component
    */
   double zz() const { return matrix_[10]; }
 
   /**
    *   The zt component
    */
   double zt() const { return matrix_[11]; }
 
   /**
    *   The tx component
    */
   double tx() const { return matrix_[12]; }
 
   /**
    *   The ty component
    */
   double ty() const { return matrix_[13]; }
 
   /**
    *   The tz component
    */
   double tz() const { return matrix_[14]; }
 
   /**
    *   The tt component
    */
   double tt() const { return matrix_[15]; }
   //@}
 
   /** @name Transformation and product members */
   //@{
 
   /**
    * Product with a LorentzVector simply returns the rotated vector.
    */
   template <typename Value>
   LorentzVector<Value>
   operator*(const LorentzVector<Value> & v) const {
     return LorentzVector<Value>
       (xx()*v.x() + xy()*v.y() + xz()*v.z() + xt()*v.t(),
        yx()*v.x() + yy()*v.y() + yz()*v.z() + yt()*v.t(),
        zx()*v.x() + zy()*v.y() + zz()*v.z() + zt()*v.t(),
        tx()*v.x() + ty()*v.y() + tz()*v.z() + tt()*v.t());
   }
 
   /**
    * Product with a Lorentz5Vector simply returns the rotated vector.
    */
   template <typename Value>
   Lorentz5Vector<Value>
   operator*(const Lorentz5Vector<Value> & v) const {
     return Lorentz5Vector<Value>
       (xx()*v.x() + xy()*v.y() + xz()*v.z() + xt()*v.t(),
        yx()*v.x() + yy()*v.y() + yz()*v.z() + yt()*v.t(),
        zx()*v.x() + zy()*v.y() + zz()*v.z() + zt()*v.t(),
        tx()*v.x() + ty()*v.y() + tz()*v.z() + tt()*v.t());
   }
 
   /**
    * Product of two LorentzRotations (this) * lt - matrix multiplication  
    * @param lt The LorentzRotation we are multiplying
    */
   SpinOneLorentzRotation operator * (const SpinOneLorentzRotation & lt) const;
 
   /**
    * Multiply by and assign a*=b becomes a= a*b
    */
   SpinOneLorentzRotation & operator *= (const SpinOneLorentzRotation & lt) {
     return *this = *this * lt;
   }
 
   /**
    *  Transform  (similar to *= but a.transform(b) becomes a = b*a
    */
   SpinOneLorentzRotation & transform   (const SpinOneLorentzRotation & lt) {
     return *this = lt * *this;
   }
 
   /**
    * Rotation around the x-axis; equivalent to LT = RotationX(delta) * LT
    */
   SpinOneLorentzRotation & rotateX(double delta) {
     SpinOneLorentzRotation tmp;
     tmp.setRotateX(delta);
     return *this = tmp * *this;
   }
 
   /**
    * Rotation around the y-axis; equivalent to LT = RotationY(delta) * LT
    */
   SpinOneLorentzRotation & rotateY(double delta) {
     SpinOneLorentzRotation tmp;
     tmp.setRotateY(delta);
     return *this = tmp * *this;
   }
 
   /**
    * Rotation around the z-axis; equivalent to LT = RotationZ(delta) * LT
    */
   SpinOneLorentzRotation & rotateZ(double delta) {
     SpinOneLorentzRotation tmp;
     tmp.setRotateZ(delta);
     return *this = tmp * *this;
   }
   
   /**
    *  Rotation around specified vector - LT = Rotation(delta,axis)*LT
    */
   SpinOneLorentzRotation & rotate(double delta, const Axis & axis) {
     SpinOneLorentzRotation tmp;
     tmp.setRotate(delta, axis);
     return *this = tmp * *this;
   }
 
   /**
    * Pure boost along the x-axis; equivalent to LT = BoostX(beta) * LT
    */
   SpinOneLorentzRotation & boostX(double beta) {
     return *this = SpinOneLorentzRotation(beta,0,0) * *this;
   }
 
   /**
    * Pure boost along the y-axis; equivalent to LT = BoostX(beta) * LT
    */
   SpinOneLorentzRotation & boostY(double beta) {
     return *this = SpinOneLorentzRotation(0,beta,0) * *this;
   }
 
   /**
    * Pure boost along the z-axis; equivalent to LT = BoostX(beta) * LT
    */
   SpinOneLorentzRotation & boostZ(double beta) {
     return *this = SpinOneLorentzRotation(0,0,beta) * *this;
   }
 
   /**
    *  boost equivalent to LT = Boost(bx,by,bz) * LT
    * @param bx The x-component of the boost
    * @param by The y-component of the boost
    * @param bz The z-component of the boost
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinOneLorentzRotation & boost(double bx, double by, double bz,
 				 double gamma=-1.) {
     return *this = SpinOneLorentzRotation(bx,by,bz,gamma) * *this;
   }
 
   /**
    *  boost equivalent to LT = Boost(bv) * LT
    * @param b The boost vector
    * @param gamma The \f$\gamma\f$ factor (optional)
    */
   SpinOneLorentzRotation & boost(const Boost & b, double gamma=-1.) {
     return *this = SpinOneLorentzRotation(b.x(),b.y(),b.z(),gamma) * *this;
   }
   //@}
 
 private:
 
   template<typename Value> friend class Helicity::LorentzTensor;
   template<typename Value> friend class Helicity::LorentzRSSpinor;
   template<typename Value> friend class Helicity::LorentzRSSpinorBar;
 
   /// Matrix components, order: \f$(xx, xy, \ldots, tz, tt)\f$.
-  vector<double> matrix_;
+  array<double,16> matrix_ = {};
 
   /// Constructor from doubles.
   SpinOneLorentzRotation (double xx, double xy, double xz, double xt,
 			  double yx, double yy, double yz, double yt,
 			  double zx, double zy, double zz, double zt,
 			  double tx, double ty, double tz, double tt);
 
   /// Component access by index: x=0, t=3.
   double operator()(unsigned int i, unsigned int j) const {
     return matrix_[4*i + j];
   }
 
   /// @name Component access.
   //@{
   double & xx_() { return matrix_[ 0]; }
   double & xy_() { return matrix_[ 1]; }
   double & xz_() { return matrix_[ 2]; }
   double & xt_() { return matrix_[ 3]; }
 
   double & yx_() { return matrix_[ 4]; }
   double & yy_() { return matrix_[ 5]; }
   double & yz_() { return matrix_[ 6]; }
   double & yt_() { return matrix_[ 7]; }
 
   double & zx_() { return matrix_[ 8]; }
   double & zy_() { return matrix_[ 9]; }
   double & zz_() { return matrix_[10]; }
   double & zt_() { return matrix_[11]; }
 
   double & tx_() { return matrix_[12]; }
   double & ty_() { return matrix_[13]; }
   double & tz_() { return matrix_[14]; }
   double & tt_() { return matrix_[15]; }
   //@}
 };
 
 /**
  *  output operator
  */
 inline std::ostream & operator<< ( std::ostream & os,
 				   const  SpinOneLorentzRotation& lt ) {
   return lt.print(os);
 }
 
 }
 
 #endif /* ThePEG_SpinOneLorentzRotation_H */