diff --git a/EventRecord/Particle.cc b/EventRecord/Particle.cc
--- a/EventRecord/Particle.cc
+++ b/EventRecord/Particle.cc
@@ -1,521 +1,523 @@
 // -*- C++ -*-
 //
 // Particle.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 member functions of
 // the Particle class.
 //
 #include "Particle.h"
 #include "ThePEG/EventRecord/Step.h"
 #include "ThePEG/EventRecord/Event.h"
 #include "ThePEG/EventRecord/ColourLine.h"
 #include "ThePEG/Utilities/Rebinder.h"
 #include "ThePEG/Config/algorithm.h"
 #include "ThePEG/EventRecord/ParticleTraits.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/PDT/DecayMode.h"
 #include <iostream>
 #include <iomanip>
 #include <cctype>
 
 #ifdef ThePEG_TEMPLATES_IN_CC_FILE
 #include "Particle.tcc"
 #endif
 
 using namespace ThePEG;
 
 Particle::ParticleRep::ParticleRep(const ParticleRep & p)
   : theParents(p.theParents), theChildren(p.theChildren),
     thePrevious(p.thePrevious), theNext(p.theNext),
     theBirthStep(p.theBirthStep), theVertex(p.theVertex),
     theLifeLength(p.theLifeLength), theScale(p.theScale),
     theVetoScale(p.theVetoScale), theNumber(p.theNumber), 
     theExtraInfo(p.theExtraInfo.size()) {
   if ( p.theColourInfo )
     theColourInfo = dynamic_ptr_cast<CBPtr>(p.theColourInfo->clone());
   if ( p.theSpinInfo )
     theSpinInfo = dynamic_ptr_cast<SpinPtr>(p.theSpinInfo->clone());
   for ( int i = 0, N = p.theExtraInfo.size(); i < N; ++i )
     theExtraInfo[i] = p.theExtraInfo[i]->clone();
 }
 
 Particle::Particle(const Particle & p)
   : Base(p), theData(p.theData), theMomentum(p.theMomentum), theRep(p.theRep) {
   if ( p.theRep ) {
     theRep = new ParticleRep(*p.theRep);
     theRep->theParents.clear();
   }
 }
 
 Particle::~Particle() {
   if ( theRep ) {
     if ( colourLine() ) colourLine()->removeColoured(this);
     if ( antiColourLine() ) antiColourLine()->removeAntiColoured(this);
     delete theRep;
   }
   theRep = 0;
   theData = cEventPDPtr();
 }
 
 void Particle::initFull() {
   if ( theRep ) return;
   theRep = new ParticleRep;
 
   Energy width = data().generateWidth(mass());
 
   if ( width > ZERO ) {
     Time lifetime = data().generateLifeTime(mass(), width);
     theRep->theLifeLength.setTau(lifetime);
     theRep->theLifeLength.
       setVect((momentum().vect()*(lifetime /
 				max(mass(), Constants::epsilon*GeV))));
     theRep->theLifeLength.rescaleEnergy();
   }
 }
 
 PPtr Particle::clone() const {
   return ptr_new<PPtr>(*this);
 }
 
 void Particle::rebind(const EventTranslationMap & trans) {
   for ( ParticleVector::iterator pit = rep().theChildren.begin();
 	pit != rep().theChildren.end(); ++pit ) *pit = trans.translate(*pit);
   for ( tParticleVector::iterator pit = rep().theParents.begin();
 	pit != rep().theParents.end(); ++pit ) *pit = trans.translate(*pit);
   rep().thePrevious = trans.translate(rep().thePrevious);
   rep().theNext = trans.translate(rep().theNext);
   if ( hasColourInfo() ) colourInfo()->rebind(trans);
   if ( spinInfo() ) spinInfo()->rebind(trans);
   rep().theBirthStep = trans.translate(rep().theBirthStep);
   for ( EIVector::const_iterator ie = rep().theExtraInfo.begin();
 	ie != rep().theExtraInfo.end(); ++ie ) (**ie).rebind(trans);
 }
 
 tParticleSet Particle::siblings() const {
   tParticleSet theSiblings;
   for ( tParticleVector::const_iterator pit = parents().begin();
 	pit != parents().end(); ++pit )
     theSiblings.insert((*pit)->children().begin(), (*pit)->children().end());
   theSiblings.erase(const_cast<Particle *>(this));
   return theSiblings;
 }
 
 void Particle::colourNeighbour(tPPtr p, bool anti) {
   tColinePtr line = colourLine(!anti);
   if ( !line ) line = ColourLine::create(this, !anti);
   line->addColoured(p, anti);  
 }
 
 void Particle::outgoingColour(tPPtr p, bool anti) {
   tColinePtr line = colourLine(anti);
   if ( !line ) line = ColourLine::create(this, anti);
   line->addColoured(p, anti);  
 }
 
 tPPtr Particle::incomingColour(bool anti) const {
   if ( !hasColourInfo() ) return tPPtr();
   tColinePtr line = colourLine(anti);
   if ( !line ) return tPPtr();
   for ( int i = 0, N = parents().size(); i < N; ++i )
     if ( parents()[i]->hasColourLine(line, anti) ) return parents()[i];
   return tPPtr();
 }
 
 tPPtr Particle::outgoingColour(bool anti) const {
   if ( !hasColourInfo() ) return tPPtr();
   tColinePtr line = colourLine(anti);
   if ( !line ) return tPPtr();
   for ( int i = 0, N = children().size(); i < N; ++i )
     if ( children()[i]->hasColourLine(line, anti) ) return children()[i];
   return tPPtr();
 }
 
 LorentzPoint Particle::labVertex() const {
   LorentzPoint r(rep().theBirthStep && rep().theBirthStep->collision()?
 		 vertex() + rep().theBirthStep->collision()->vertex():
 		 vertex());
   return r;
 }
 
 void Particle::setLabVertex(const LorentzPoint & p) {
   rep().theVertex = ( rep().theBirthStep && rep().theBirthStep->collision()?
 		      p - rep().theBirthStep->collision()->vertex() : p );
 }
 
 void Particle::transform(const LorentzRotation & r) {
   if ( hasRep() && spinInfo() ) spinInfo()->transform(momentum(), r);
   theMomentum.transform(r);
   if ( !hasRep() ) return;
   rep().theVertex.transform(r);
   rep().theLifeLength.transform(r);
 }
 
 void Particle::deepTransform(const LorentzRotation & r) {
   transform(r);
   if ( !theRep ) return;
   for ( int i = 0, N = children().size(); i < N; ++i )
     rep().theChildren[i]->deepTransform(r);
   if ( rep().theNext ) rep().theNext->deepTransform(r);
 }
 
 void Particle::rotateX(double a) {
   LorentzRotation r;
   r.rotateX(a);
   transform(r);
 }
 
 void Particle::deepRotateX(double a) {
   LorentzRotation r;
   r.rotateX(a);
   deepTransform(r);
 }
 
 void Particle::rotateY(double a) {
   LorentzRotation r;
   r.rotateY(a);
   transform(r);
 }
 
 void Particle::deepRotateY(double a) {
   LorentzRotation r;
   r.rotateY(a);
   deepTransform(r);
 }
 
 void Particle::rotateZ(double a) {
   LorentzRotation r;
   r.rotateZ(a);
   transform(r);
 }
 
 void Particle::deepRotateZ(double a) {
   LorentzRotation r;
   r.rotateZ(a);
   deepTransform(r);
 }
 
 void Particle::rotate(double a, const Axis & axis) {
   LorentzRotation r;
   r.rotate(a, axis);
   transform(r);
 }
 
 void Particle::deepRotate(double a, const Axis & axis) {
   LorentzRotation r;
   r.rotate(a, axis);
   deepTransform(r);
 }
 
 string Particle::outputFormat =
 "%n3%s10 %i7 %p[,]0 %c(,) %^^0%vv0 %>>0%<>0 %l{,}0\n"
 "                            %x10.3%y10.3%z10.3%e10.3%m10.3\n";
 
 int getNumber(string::const_iterator & pos, int def) {
   if ( !isdigit(*pos) ) return def;
   def = *pos++ - '0';
   while ( isdigit(*pos) ) def = 10*def + *pos++ - '0';
   return def;
 }
 
 void writePrecision(ostream & os, string::const_iterator & pos,
 			      int defw, int defp, double x) {
   defw = getNumber(pos, defw);
   if ( *pos == '.' ) defp = getNumber(++pos, defp);
   int oldp = os.precision();
   os << setprecision(defp) << setw(defw) << x << setprecision(oldp);
 }
 
 void writeStringAdjusted(ostream & os, bool left, int w, string str) {
   while ( !left && w-- > int(str.size()) ) os << ' ';
   os << str;
   while ( left && w-- > int(str.size()) ) os << ' ';
 }
 
 template <typename Container>
 void writeParticleRanges(ostream & os, const Container & co, char sep, int w) {
   set<int> cnum;
   for ( typename Container::const_iterator it = co.begin();
 	it != co.end(); ++it) cnum.insert((**it).number());
 
   bool elipsis = false;
   int last = -10;
   for ( set<int>::iterator it = cnum.begin(); it != cnum.end(); ++it) {
     int n = *it;
     int next = 0;
     set<int>::iterator itn = it;
     if ( ++itn != cnum.end() ) next = *itn;
     bool writeit = true;
     bool writesep = false;
     if ( elipsis && ( n != last + 1 || n != next - 1 ) )
       elipsis = false;
     else if ( !elipsis && n == last + 1 && n == next -1 ) {
       os << "..";
       elipsis = true;
       writeit = false;
     }
     else if ( elipsis && n == last + 1 && n == next -1 )
       writeit = false;
     else if ( it != cnum.begin() )
       writesep = true;
     if ( writeit ) {
       if ( writesep ) os << sep;
       os << setw(w) << n;
     }
     last = n;
   }
 }
 
 ostream & ThePEG::operator<<(ostream & os, const Particle & p) {
   return p.print(os, p.birthStep());
 }
 
 ostream & Particle::print(ostream & os, tcStepPtr step) const {
   if ( !step ) step = birthStep();
   tCollPtr coll = step? step->collision(): tCollPtr();
   tEventPtr event = coll? coll->event(): tEventPtr();
   string::const_iterator pos = Particle::outputFormat.begin();
   ios::fmtflags saveflags = os.setf(ios::fixed, ios::floatfield);
   while ( pos != Particle::outputFormat.end() ) {
     if ( *pos == '%' && ++pos != Particle::outputFormat.end() ) {
       bool left = false;
       if ( *pos == '-' ) {
 	left = true;
 	os.setf(ios::left, ios::adjustfield);
 	++pos;
       } else {
 	os.setf(ios::right, ios::adjustfield);
       }
       char mark;
       char open;
       char close;
       char sep;
       int w;
       string str;
       string fill;
       if ( pos == Particle::outputFormat.end() ) break;
       bool fullColour = false;
       switch ( *pos ) {
       case 'n':
 	os << setw(getNumber(++pos, 3)) << number();
 	break;
       case 'i':
 	os << setw(getNumber(++pos, 8)) << id();
 	break;
       case 's':
 	writeStringAdjusted(os, left, getNumber(++pos, 8), PDGName());
 	break;
       case 'x':
 	writePrecision(os, ++pos, 10, 3, momentum().x()/GeV);
 	break;
       case 'y':
 	writePrecision(os, ++pos, 10, 3, momentum().y()/GeV);
 	break;
       case 'z':
 	writePrecision(os, ++pos, 10, 3, momentum().z()/GeV);
 	break;
       case 'e':
 	writePrecision(os, ++pos, 10, 3, momentum().e()/GeV);
 	break;
       case 'm':
 	writePrecision(os, ++pos, 10, 3, momentum().mass()/GeV);
 	break;
       case 'P':
 	fullColour = true;
+	[[fallthrough]];
       case 'p':
 	open = *++pos;
 	sep = *++pos;
 	close = *++pos;
 	w = getNumber(++pos, 0);
 	if ( parents().empty() ) break;
 	if ( open ) os << open;
 	writeParticleRanges(os, parents(), sep, w);
 	if ( fullColour && hasColourInfo() &&
 	     ( incomingColour() || incomingAntiColour() ) ) {
 	  if ( close ) os << open;
 	  if ( incomingColour() )
 	    os << "+" << incomingColour()->number();
 	  if ( incomingAntiColour() )
 	    os << "-" << incomingAntiColour()->number();
 	  if ( close ) os << close;
 	}
 	if ( close ) os << close;
 	break;
       case 'l':
 	open = *++pos;
 	sep = *++pos;
 	close = *++pos;
 	w = getNumber(++pos, 0);
 	if ( hasColourInfo() &&
 	     ( colourLine() || antiColourLine() ) && event) {
 	  if ( open ) os << open;
 	  vector<tcColinePtr> clines = colourInfo()->colourLines();
 	  for ( int i = 0, N = clines.size(); i < N; ++i ) {
 	    if ( i > 0 && sep )  os << sep;
 	    clines[i]->write(os, event, false);
 	  }
 	  vector<tcColinePtr> aclines = colourInfo()->antiColourLines();
 	  for ( int i = 0, N = aclines.size(); i < N; ++i ) {
 	    if ( ( i > 0 || clines.size() ) && sep )  os << sep;
 	    aclines[i]->write(os, event, true);
 	  }
 	  if ( close ) os << close;
 	}
 	break;
       case 'C':
 	fullColour = true;
+	[[fallthrough]];
       case 'c':
 	open = *++pos;
 	sep = *++pos;
 	close = *++pos;
 	w = getNumber(++pos, 0);
 	if ( children().empty() ) break;
 	if ( open ) os << open;
 	writeParticleRanges(os, children(), sep, w);
 	if ( fullColour && hasColourInfo() &&
 	     ( outgoingColour() || outgoingAntiColour() ) ) {
 	  if ( close ) os << open;
 	  if ( outgoingColour() )
 	    os << "+" << outgoingColour()->number();
 	  if ( outgoingAntiColour() )
 	    os << "-" << outgoingAntiColour()->number();
 	  if ( close ) os << close;
 	}
 	if ( close ) os << close;
 	break;
       case '>':
 	mark = *++pos;
 	w = getNumber(++pos, 0);
 	if ( hasColourInfo() && step && step->colourNeighbour(this) ) {
 	  os << setw(w-1) << step->colourNeighbour(this)->number() << mark;
 	}
 	break;
       case '<':
 	mark = *++pos;
 	w = getNumber(++pos, 0);
 	if ( hasColourInfo() && step && step->antiColourNeighbour(this) ) {
 	  int n = step->antiColourNeighbour(this)->number();
 	  ostringstream oss;
 	  oss << mark << n;
 	  writeStringAdjusted(os, left, w, oss.str());
 	}
 	break;
       case 'v':
 	mark = *++pos;
 	w = getNumber(++pos, 0);
 	if ( next() ) {
 	  if ( left && mark ) os << mark;
 	  os << setw(w) << next()->number();
 	  if ( !left && mark ) os << mark;
 	}
 	break;
       case '^':
 	mark = *++pos;
 	w = getNumber(++pos, 0);
 	if ( previous() ) {
 	  if ( left && mark ) os << mark;
 	  os << setw(w) << previous()->number();
 	  if ( !left && mark ) os << mark;
 	}
 	break;
       case 'd':
 	switch ( *++pos ) {
 	case 'x':
 	  writePrecision(os, ++pos, 10, 3, lifeLength().x()/mm);
 	  break;
 	case 'y':
 	  writePrecision(os, ++pos, 10, 3, lifeLength().y()/mm);
 	  break;
 	case 'z':
 	  writePrecision(os, ++pos, 10, 3, lifeLength().z()/mm);
 	  break;
 	case 't':
 	  writePrecision(os, ++pos, 10, 3, lifeLength().e()/mm);
 	  break;
 	case 'T':
 	  writePrecision(os, ++pos, 10, 3, lifeLength().tau()/mm);
 	  break;
 	}
 	break;
       case 'V':
 	switch ( *++pos ) {
 	case 'x':
 	  writePrecision(os, ++pos, 10, 3, vertex().x()/mm);
 	  break;
 	case 'y':
 	  writePrecision(os, ++pos, 10, 3, vertex().y()/mm);
 	  break;
 	case 'z':
 	  writePrecision(os, ++pos, 10, 3, vertex().z()/mm);
 	  break;
 	case 't':
 	  writePrecision(os, ++pos, 10, 3, vertex().e()/mm);
 	  break;
 	}
       case 'L':
 	switch ( *++pos ) {
 	case 'x':
 	  writePrecision(os, ++pos, 10, 3, labVertex().x()/mm);
 	  break;
 	case 'y':
 	  writePrecision(os, ++pos, 10, 3, labVertex().y()/mm);
 	  break;
 	case 'z':
 	  writePrecision(os, ++pos, 10, 3, labVertex().z()/mm);
 	  break;
 	case 't':
 	  writePrecision(os, ++pos, 10, 3, labVertex().e()/mm);
 	  break;
 	}
 	break;
       default:
 	os << *pos++;
       }
     } else {
       if ( pos != Particle::outputFormat.end() ) os << *pos++;
     }
   }
   os.flags(saveflags);
   return os;
 }
 
 void Particle::debugme() const {
   cerr << *this;
   EventRecordBase::debugme();
 }
 
 void Particle::persistentOutput(PersistentOStream & os) const {
   EventConfig::putParticleData(os, theData);
   os << ounit(theMomentum, GeV) << bool( theRep != 0 );
   if ( !theRep ) return;
   os << rep().theParents << rep().theChildren
      << rep().thePrevious << rep().theNext << rep().theBirthStep
      << ounit(rep().theVertex, mm) << ounit(rep().theLifeLength, mm)
      << ounit(rep().theScale, GeV2) << ounit(rep().theVetoScale, GeV2) 
      << rep().theNumber << rep().theDecayMode
      << rep().theColourInfo << rep().theSpinInfo << rep().theExtraInfo;
 }
 
 void Particle::persistentInput(PersistentIStream & is, int) {
   bool readRep;
   EventConfig::getParticleData(is, theData);
   is >> iunit(theMomentum, GeV) >> readRep;
   if ( !readRep ) return;
   if ( !hasRep() ) theRep = new ParticleRep;
 
   is >> rep().theParents >> rep().theChildren
      >> rep().thePrevious >> rep().theNext >> rep().theBirthStep
      >> iunit(rep().theVertex, mm) >> iunit(rep().theLifeLength, mm)
      >> iunit(rep().theScale, GeV2) >> iunit(rep().theVetoScale, GeV2) 
      >> rep().theNumber >> rep().theDecayMode
      >> rep().theColourInfo >> rep().theSpinInfo >> rep().theExtraInfo;
 }
 
 ClassDescription<Particle> Particle::initParticle;
 
 void Particle::Init() {}
diff --git a/MatrixElement/ME2to2Base.cc b/MatrixElement/ME2to2Base.cc
--- a/MatrixElement/ME2to2Base.cc
+++ b/MatrixElement/ME2to2Base.cc
@@ -1,202 +1,202 @@
 // -*- C++ -*-
 //
 // ME2to2Base.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 ME2to2Base class.
 //
 
 #include "ME2to2Base.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Utilities/SimplePhaseSpace.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Handlers/StandardXComb.h"
 #include "ThePEG/Cuts/Cuts.h"
 
 using namespace ThePEG;
 
 ME2to2Base::~ME2to2Base() {}
 
 Energy2 ME2to2Base::scale() const {
   switch ( scaleChoice() ) {
   case 1:
     return -tHat()*uHat()/(tHat() + uHat());
   default:
     return tHat()*uHat()/sHat();
   }
 }
 
 void ME2to2Base::setKinematics() {
   MEBase::setKinematics();
   theLastTHat = (meMomenta()[0] - meMomenta()[2]).m2();
   theLastUHat = (meMomenta()[1] - meMomenta()[2]).m2();
   theLastPhi = meMomenta()[2].phi();
 }
 
 bool ME2to2Base::generateKinematics(const double * r) {
   // generate the masses of the particles
   for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
     meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->generateMass());
   }
 
   double ctmin = -1.0;
   double ctmax = 1.0;
 
   Energy q = ZERO;
   try {
     q = SimplePhaseSpace::
       getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass());
-  } catch ( ImpossibleKinematics ) {
+  } catch ( ImpossibleKinematics & e ) {
     return false;
   }
 
   Energy e = sqrt(sHat())/2.0;
 		    
   Energy2 m22 = meMomenta()[2].mass2();
   Energy2 m32 = meMomenta()[3].mass2();
   Energy2 e0e2 = 2.0*e*sqrt(sqr(q) + m22);
   Energy2 e1e2 = 2.0*e*sqrt(sqr(q) + m22);
   Energy2 e0e3 = 2.0*e*sqrt(sqr(q) + m32);
   Energy2 e1e3 = 2.0*e*sqrt(sqr(q) + m32);
   Energy2 pq = 2.0*e*q;
 
   Energy2 thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[2]);
   if ( thmin > ZERO ) ctmax = min(ctmax, (e0e2 - m22 - thmin)/pq);
 
   thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[2]);
   if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m22 - e1e2)/pq);
 
   thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[3]);
   if ( thmin > ZERO ) ctmax = min(ctmax, (e1e3 - m32 - thmin)/pq);
 
   thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[3]);
   if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m32 - e0e3)/pq);
 
   Energy ptmin = max(lastCuts().minKT(mePartonData()[2]),
    		     lastCuts().minKT(mePartonData()[3]));
   if ( ptmin > ZERO ) {
     double ctm = 1.0 - sqr(ptmin/q);
     if ( ctm <= 0.0 ) return false;
     ctmin = max(ctmin, -sqrt(ctm));
     ctmax = min(ctmax, sqrt(ctm));
   }
 
   double ymin2 = lastCuts().minYStar(mePartonData()[2]);
   double ymax2 = lastCuts().maxYStar(mePartonData()[2]);
   double ymin3 = lastCuts().minYStar(mePartonData()[3]);
   double ymax3 = lastCuts().maxYStar(mePartonData()[3]);
   double ytot = lastCuts().Y() + lastCuts().currentYHat();
   if ( ymin2 + ytot > -0.9*Constants::MaxRapidity )
     ctmin = max(ctmin, sqrt(sqr(q) +  m22)*tanh(ymin2)/q);
   if ( ymax2 + ytot < 0.9*Constants::MaxRapidity )
     ctmax = min(ctmax, sqrt(sqr(q) +  m22)*tanh(ymax2)/q);
   if ( ymin3 + ytot > -0.9*Constants::MaxRapidity )
     ctmax = min(ctmax, sqrt(sqr(q) +  m32)*tanh(-ymin3)/q);
   if ( ymax3 + ytot < 0.9*Constants::MaxRapidity )
     ctmin = max(ctmin, sqrt(sqr(q) +  m32)*tanh(-ymax3)/q);
 
   if ( ctmin >= ctmax ) return false;
     
   double cth = getCosTheta(ctmin, ctmax, r);
 
   Energy pt = q*sqrt(1.0-sqr(cth));
   theLastPhi = rnd(2.0*Constants::pi);
   meMomenta()[2].setVect(Momentum3( pt*sin(phi()),  pt*cos(phi()),  q*cth));
 
   meMomenta()[3].setVect(Momentum3(-pt*sin(phi()), -pt*cos(phi()), -q*cth));
 
   meMomenta()[2].rescaleEnergy();
   meMomenta()[3].rescaleEnergy();
 
   vector<LorentzMomentum> out(2);
   out[0] = meMomenta()[2];
   out[1] = meMomenta()[3];
   tcPDVector tout(2);
   tout[0] = mePartonData()[2];
   tout[1] = mePartonData()[3];
   if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) )
     return false;
 
   theLastTHat = pq*cth + m22 - e0e2;
   theLastUHat = m22 + m32 - sHat() - theLastTHat;
   jacobian((pq/sHat())*Constants::pi*jacobian());
   return true;
 }
 
 double ME2to2Base::getCosTheta(double ctmin, double ctmax, const double * r) {
   double cth = 0.0;
   static const double eps = 1.0e-6;
   if ( 1.0 + ctmin <= eps && 1.0 - ctmax <= eps ) {
     jacobian(ctmax - ctmin);
     cth = ctmin + (*r)*jacobian();
   } else if (  1.0 + ctmin <= eps ) {
     cth = 1.0 - (1.0 - ctmax)*pow((1.0 - ctmin)/(1.0 - ctmax), *r);
     jacobian(log((1.0 - ctmin)/(1.0 - ctmax))*(1.0 - cth));
   } else if (  1.0 - ctmax <= eps ) {
     cth = -1.0 + (1.0 + ctmin)*pow((1.0 + ctmax)/(1.0 + ctmin), *r);
     jacobian(log((1.0 + ctmax)/(1.0 + ctmin))*(1.0 + cth));
   } else {
     double zmin = 0.5*(1.0 - ctmax);
     double zmax = 0.5*(1.0 - ctmin);
     double A1 = -ctmin/(zmax*(1.0-zmax));
     double A0 = -ctmax/(zmin*(1.0-zmin));
     double A = *r*(A1 - A0) + A0;
     double z = A < 2.0? 2.0/(sqrt(sqr(A) + 4.0) + 2 - A):
       0.5*(A - 2.0 + sqrt(sqr(A) + 4.0))/A;
     cth = 1.0 - 2.0*z;
     jacobian(2.0*(A1 - A0)*sqr(z)*sqr(1.0 - z)/(sqr(z) + sqr(1.0 - z)));
   }
   return cth;
 }
 
 CrossSection ME2to2Base::dSigHatDR() const {
   return me2()*jacobian()/(16.0*sqr(Constants::pi)*sHat())*sqr(hbarc);
 }
 
 void ME2to2Base::persistentOutput(PersistentOStream & os) const {
   os << theScaleChoice << ounit(theLastTHat, GeV2) << ounit(theLastUHat, GeV2)
      << theLastPhi;
 }
 
 void ME2to2Base::persistentInput(PersistentIStream & is, int) {
   is >> theScaleChoice >> iunit(theLastTHat, GeV2) >> iunit(theLastUHat, GeV2)
      >> theLastPhi;
 }
 
 AbstractClassDescription<ME2to2Base> ME2to2Base::initME2to2Base;
 // Definition of the static class description member.
 
 Switch<ME2to2Base,int> & ME2to2Base::interfaceScaleChoice() {
   static Switch<ME2to2Base,int> dummy
     ("ScaleChoice",
      "Different options for calculating the scale of the generated "
      "hard sub-process.",
      &ME2to2Base::theScaleChoice, 0, false, false);
   return dummy;
 }
 
 void ME2to2Base::Init() {
 
   static ClassDocumentation<ME2to2Base> documentation
     ("The ThePEG::ME2to2Base class may be used as a base class "
      "for all \\f$2\\rightarrow 2\\f$ matrix elements.");
 
   static SwitchOption interfaceScaleChoice0
     (interfaceScaleChoice(),
      "that.uhat/shat", "\\f$\\hat{t}\\hat{u}/\\hat{s}\\f$", 0);
 
   static SwitchOption interfaceScaleChoice1
     (interfaceScaleChoice(),
      "that.uhat/(that+uhat)",
      "\\f$-\\hat{t}\\hat{u}/(\\hat{t}+\\hat{u})\\f$", 1);
 
 }
 
diff --git a/PDF/PartonExtractor.cc b/PDF/PartonExtractor.cc
--- a/PDF/PartonExtractor.cc
+++ b/PDF/PartonExtractor.cc
@@ -1,653 +1,653 @@
 // -*- C++ -*-
 //
 // PartonExtractor.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 PartonExtractor class.
 //
 
 #include "PartonExtractor.h"
 #include "ThePEG/Handlers/XComb.h"
 #include "ThePEG/Cuts/Cuts.h"
 #include "ThePEG/PDF/NoPDF.h"
 #include "ThePEG/PDF/RemnantHandler.h"
 #include "ThePEG/PDF/BeamParticleData.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/EventRecord/Particle.h"
 #include "ThePEG/EventRecord/Step.h"
 #include "ThePEG/EventRecord/SubProcess.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/RefVector.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/SimplePhaseSpace.h"
 #include "ThePEG/Utilities/UtilityBase.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/PDT/EnumParticles.h"
 
 using namespace ThePEG;
 
 PartonExtractor::PartonExtractor()
   : theMaxTries(100), flatSHatY(false) {}
 
 PartonExtractor::~PartonExtractor() {}
 
 IBPtr PartonExtractor::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr PartonExtractor::fullclone() const {
   return new_ptr(*this);
 }
 
 PartonPairVec PartonExtractor::
 getPartons(Energy maxEnergy, const cPDPair & incoming,
 	   const Cuts & kc) const {
   PartonPairVec result;
   PartonVector first;
   PDFCuts cuts1(kc, true, maxEnergy);
   PBPtr p1 =
     new_ptr(PartonBin(PDPtr(), PBPtr(), incoming.first, PDFPtr(), cuts1));
   addPartons(p1, cuts1,  theFirstPDF, first);
   PartonVector second;
   PDFCuts cuts2(kc, false, maxEnergy);
   PBPtr p2 =
     new_ptr(PartonBin(PDPtr(), PBPtr(), incoming.second, PDFPtr(), cuts2));
   addPartons(p2, cuts2, theSecondPDF, second);
   for ( PartonVector::iterator it1 = first.begin();
 	it1 != first.end(); ++it1 )
     for ( PartonVector::iterator it2 = second.begin();
 	it2 != second.end(); ++it2 )
       result.push_back(PBPair(*it1, *it2));
 
   // We add the original parton bins as well to avoid them being
   // deleted.
   result.push_back(PBPair(p1, p2));
   return result;
 }
 
 void PartonExtractor::
 addPartons(tPBPtr incoming, const PDFCuts & cuts, tcPDFPtr pdf,
 	   PartonVector & pbins) const {
   if(!pdf) pdf = getPDF(incoming->parton());
   if ( dynamic_ptr_cast<Ptr<NoPDF>::tcp>(pdf) ||
        incoming->parton() == incoming->particle() ) {
     pbins.push_back(incoming);
     return;
   }
   cPDVector partons = pdf->partons(incoming->parton());
   for ( int i = 0, N = partons.size(); i < N; ++i ) {
     PBPtr pb =
       new_ptr(PartonBin(incoming->parton(), incoming, partons[i], pdf, cuts));
     incoming->addOutgoing(pb);
     addPartons(pb, cuts, PDFPtr(), pbins);
   }
 
 }
 
 tcPDFPtr PartonExtractor::getPDF(tcPDPtr particle) const {
   for ( vector<PDFPtr>::const_iterator it = theSpecialDensities.begin();
 	it != theSpecialDensities.end(); ++it )
     if ( (**it).canHandle(particle) ) return *it;
   Ptr<BeamParticleData>::tcp p = 
     dynamic_ptr_cast<Ptr<BeamParticleData>::tcp>(particle);
   if ( !p || !p->pdf() ) return noPDF();
   return p->pdf();
 }
 
 void PartonExtractor::select(tXCombPtr newXComb) {
   theLastXComb = newXComb;
 }
 
 tPBIPtr PartonExtractor::partonBinInstance(tcPPtr p) const {
   PartonBinInstanceMap::const_iterator it = partonBinInstances().find(p);
   return it == partonBinInstances().end()? PBIPtr(): it->second;
 }
 
 void PartonExtractor::
 colourConnect(tPPtr particle, tPPtr parton, const tPVector & remnants) const {
 
   // Sorry cannot handle coloured resolved particles.
   if ( particle->coloured() ) throw RemColException(*this);
 
   // First connect the loose colour line from the extacted parton.
   if ( parton->hasColour() )
     findConnect(parton->colourLine(), parton, true,
 		remnants.rbegin(), remnants.rend());
 
   // First connect the loose anti-colour line from the extacted parton.
   if ( parton->hasAntiColour() )
     findConnect(parton->antiColourLine(), parton, false,
 		remnants.begin(), remnants.end());
 
   // Go through the rest of the remnants and create new colour lines
   // if needed. Go through it forwards and backwards to catch possible
   // inconsistencies.
   for ( tPVector::const_iterator it = remnants.begin();
 	it != remnants.end(); ++it ) {
     if ( (**it).hasAntiColour() && !(**it).antiColourLine() )
       findConnect(ColourLine::create(*it, true), *it, false,
 		it + 1, remnants.end());
   }
   for ( tPVector::const_reverse_iterator it = remnants.rbegin();
 	it != remnants.rend(); ++it ) {
     if ( (**it).hasColour() && !(**it).colourLine() )
       findConnect(ColourLine::create(*it), *it, true, it + 1, remnants.rend());
   } 
 }
 
 Energy2 PartonExtractor::newScale() {
   return lastScale();
 }
 
 pair<int,int> PartonExtractor::nDims(const PBPair & pbins) {
   // if photon from a lepton or proton generate scale
   bool genscale[2]={false,false};
   for(unsigned int ix=0;ix<2;++ix) {
     PBPtr bin = ix==0 ? pbins.first : pbins.second;
     if (!bin || !bin->particle() || !bin->parton()) continue;
     if(bin->pdf()->partons(bin->particle()).size()==1 &&
        bin->particle()->id()!=bin->parton()->id())
       genscale[ix]=true;
   }
   return make_pair(pbins.first ->nDim(genscale[0]),
 		   pbins.second->nDim(genscale[1]));
 }
 
 void PartonExtractor::prepare(const PBIPair & pbins) {
   partonBinInstances().clear();
   pbins.first->prepare();
   pbins.second->prepare();
 }
 
 void PartonExtractor::updatePartonBinInstances(const PBIPair & pbins) {
   partonBinInstances().clear();
   tPBIPtr current = pbins.first;
   while ( current->incoming() ) {
     partonBinInstances()[current->parton()] = current;
     current = current->incoming();
   }
   current = pbins.second;
   while ( current->incoming() ) {
     partonBinInstances()[current->parton()] = current;
     current = current->incoming();
   }
 }
 
 bool PartonExtractor::
 generateL(const PBIPair & pbins, const double * r1, const double * r2) {
   Direction<0> dir(true);
   generateL(*pbins.first, r1);
   dir.reverse();
   generateL(*pbins.second, r2);
   if ( !flatSHatY || pbins.first->hasPoleIn1() || pbins.second->hasPoleIn1() )
     return true;
 
   Energy2 shmax = lastCuts().sHatMax();
   Energy2 shmin = lastCuts().sHatMin();
   Energy2 sh = shmin*pow(shmax/shmin, *r1);
   double ymax = lastCuts().yHatMax();
   double ymin = lastCuts().yHatMin();
   double km = log(shmax/shmin);
   ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/sh)));
   ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/sh)));
 
   double y = ymin + (*r2)*(ymax - ymin);
   double l1 = 0.5*log(lastS()/sh) - y;
   double l2 = 0.5*log(lastS()/sh) + y;
 
   pbins.first->li(l1 - pbins.first->l() + pbins.first->li());
   pbins.first->l(l1);
   pbins.first->jacobian(km*(ymax - ymin));
   pbins.second->li(l2 - pbins.second->l() + pbins.second->li());
   pbins.second->l(l2);
   pbins.second->jacobian(1.0);
   return ( pbins.first->li() >= 0.0 && pbins.second->li() >= 0.0 );
 }
 
 Energy2 PartonExtractor::
 generateSHat(Energy2, const PBIPair & pbins,
 	     const double * r1, const double * r2,
 	     bool haveMEPartons) {
   Direction<0> dir(true);
   if(pbins.first->bin()->pdfDim()<=1) pbins.first->scale(-lastScale());
   if ( !generate(*pbins.first, r1, lastSHat(),
 		 pbins.first->getFirst()->parton()->momentum(),
 		 haveMEPartons) )
     return -1.0*GeV2;
   dir.reverse();
   if(pbins.second->bin()->pdfDim()<=1) pbins.second->scale(-lastScale());
   if ( !generate(*pbins.second, r2, lastSHat(),
 		 pbins.second->getFirst()->parton()->momentum(),
 		 haveMEPartons) )
     return -1.0*GeV2;
   
   return (pbins.first->parton()->momentum() +
 	  pbins.second->parton()->momentum()).m2();
 }
 
 void PartonExtractor::
 generateL(PartonBinInstance & pb, const double * r) {
   if ( !pb.incoming() ) return;
 
   pb.parton(pb.partonData()->produceParticle(Lorentz5Momentum()));
   generateL(*pb.incoming(), r + pb.bin()->pdfDim() + pb.bin()->remDim());
   pb.particle(pb.incoming()->parton());
 
   if ( pb.li() >= 0 ) return;
 
   double jac = 1.0;
   if ( pb.bin()->pdfDim() )
     pb.li(pb.pdf()->flattenL(pb.particleData(), pb.partonData(),
 			     pb.bin()->cuts(), *r++, jac));
   pb.scale(-1.0*GeV2);
   if ( pb.bin()->pdfDim() > 1 )
     pb.scale(pb.pdf()->flattenScale(pb.particleData(), pb.partonData(),
 				    pb.bin()->cuts(), pb.li(), *r++, jac)
 	     *pb.bin()->cuts().scaleMaxL(pb.li()));
   pb.jacobian(jac);
   pb.l(pb.incoming()->l() + pb.li());
 }
 
 bool PartonExtractor::
 generate(PartonBinInstance & pb, const double * r,
 	 Energy2 shat, const Lorentz5Momentum & first,
 	 bool haveMEPartons) {
   if ( !pb.incoming() ) return true;
   if ( !generate(*pb.incoming(), r + pb.bin()->pdfDim() + pb.bin()->remDim(),
 		 shat/pb.xi(), first) )
     return false;
   pb.remnantWeight(1.0);
   pb.parton()->setMomentum
     (pb.remnantHandler()->generate(pb, r + pb.bin()->pdfDim(), pb.scale(), shat,
 				   pb.particle()->momentum(),haveMEPartons));
   if ( pb.remnantWeight() <= 0.0 ) return false;
   partonBinInstances()[pb.parton()] = &pb;
   return true;
 }
 
 void PartonExtractor::
 constructRemnants(const PBIPair & pbins, tSubProPtr sub, tStepPtr step) const {
   partonBinInstances().clear();
   LorentzMomentum k1 = pbins.first->parton()->momentum();
   LorentzMomentum k2 = pbins.second->parton()->momentum();
   LorentzMomentum Ph = k1 + k2;
   LorentzMomentum Phold = Ph;
   LorentzRotation Rh = Utilities::getBoostToCM(make_pair(k1, k2));
 
   bool pickside = rndbool();
   if ( pickside && pbins.first->incoming() ) {
     Direction<0> dir(true);
     constructRemnants(*pbins.first, Ph, k2);
     construct(*pbins.first, step, false);
   }
   if ( pbins.second->incoming() ) {
     Direction<0> dir(false);
     constructRemnants(*pbins.second, Ph, pbins.first->parton()->momentum());
     construct(*pbins.second, step, false);
   }
   if ( (!pickside) && pbins.first->incoming() ) {
     Direction<0> dir(true);
     constructRemnants(*pbins.first, Ph, pbins.second->parton()->momentum());
     construct(*pbins.first, step, false);
   }
   //  LorentzRotation rot = Utilities::transformToMomentum(Phold, Ph);
   k1 = pbins.first->parton()->momentum();
   k2 = pbins.second->parton()->momentum();
   LorentzRotation rot = Utilities::getBoostFromCM(make_pair(k1, k2))*Rh;
   Utilities::transform(sub->outgoing(), rot);
   Utilities::transform(sub->intermediates(), rot);
   Ph = k1 + k2;
   if ( abs(Ph.m2() - Phold.m2())/Phold.m2() > 0.000001 )
     cerr << Ph.m2()/GeV2 << " was (" << Phold.m2()/GeV2 << ")" << endl;
 }
 
 void PartonExtractor::
 constructRemnants(PartonBinInstance & pb, LorentzMomentum & Ph,
 		  const LorentzMomentum & k) const {
   LorentzMomentum P = pb.particle()->momentum();
   DVector r = UseRandom::rndvec(pb.bin()->remDim());
   if ( r.empty() ) r.push_back(0.0);
   pb.parton()->setMomentum(pb.remnantHandler()->
 			   generate(pb, &r[0], pb.scale(), Ph.m2(), P));
   if ( pb.remnantWeight() <= 0.0 ) throw Veto();
   pb.remnantHandler()->boostRemnants(pb);
   LorentzMomentum Pr = Utilities::sumMomentum(pb.remnants());
   transformRemnants(Ph, Pr, k, pb.particle()->momentum());
   pb.parton()->setMomentum(pb.particle()->momentum() - Pr);
   try {
     Utilities::setMomentum(pb.remnants().begin(),
 			   pb.remnants().end(),
 			   static_cast<const LorentzMomentum &>(Pr));
   }
   catch ( ThePEG::Exception & e) {
     throw e;
   }
   catch ( ThePEG::Veto ) {
     throw;
   }
   catch ( std::exception & e ) {
     throw Exception() << "Caught non-ThePEG exception " << e.what() << "in "
 		      << "PartonExtractor::constructRemnants"
 		      << Exception::eventerror;
   }
   partonBinInstances()[pb.parton()] = &pb;
   if ( !pb.incoming()->incoming() ) return;
 
   // We get here if we need to construct remnants recursively.
   LorentzMomentum Phnew = Ph + Pr;
   constructRemnants(*pb.incoming(), Phnew, k);
   LorentzRotation rot = Utilities::transformToMomentum(Ph + Pr, Phnew);
   Utilities::transform(pb.remnants(), rot);
   Ph.transform(rot);
 }
 
 LorentzRotation PartonExtractor::
 boostRemnants(PBIPair & bins, LorentzMomentum k1, LorentzMomentum k2,
 	 bool side1, bool side2) const {
   if ( !side1 && !side2 ) return LorentzRotation();
 
   LorentzMomentum P1 =
     bins.first? LorentzMomentum(bins.first->parton()->momentum()): k1;
   LorentzMomentum Pr1;
   if ( side1 ) {
     P1 = bins.first->particle()->momentum();
     Pr1 = Utilities::sumMomentum(bins.first->remnants());
   }
   LorentzMomentum P2 =
     bins.second? LorentzMomentum(bins.second->parton()->momentum()): k2;
   LorentzMomentum Pr2;
   if ( side2 ) {
     P2 = bins.second->particle()->momentum();
     Pr2 = Utilities::sumMomentum(bins.second->remnants());
   }
 
   LorentzRotation Rh = Utilities::getBoostToCM(make_pair(k1, k2));
   LorentzMomentum Ph = k1 + k2;
   //  LorentzMomentum Phold = Ph;
 
   bool otherside = rndbool();
   if ( otherside && side2 ){
     Direction<0> dir(false);
     transformRemnants(Ph, Pr2, k1, P2);
     k2 = P2 - Pr2;
   }
   if ( side1 ){
     Direction<0> dir(true);
     transformRemnants(Ph, Pr1, k2, P1);
     k1 = P1 - Pr1;
   }
   if ( side2 && !otherside ) {
     Direction<0> dir(false);
     transformRemnants(Ph, Pr2, k1, P2);
     k2 = P2 - Pr2;
   }
 
   if ( bins.first ) {
     if ( bins.first->remnants().size() == 1 )
       bins.first->remnants()[0]->setMomentum(Pr1);
     else
       Utilities::setMomentum(bins.first->remnants().begin(),
 			     bins.first->remnants().end(),
 			     static_cast<const LorentzMomentum &>(Pr1));
     bins.first->parton()->setMomentum(k1);
   }
 
   if ( bins.second ) {
     if ( bins.second->remnants().size() == 1 )
       bins.second->remnants()[0]->setMomentum(Pr2);
     else
       Utilities::setMomentum(bins.second->remnants().begin(),
 			     bins.second->remnants().end(),
 			     static_cast<const LorentzMomentum &>(Pr2));
     bins.second->parton()->setMomentum(k2);
   }
 
   Rh.transform(Utilities::getBoostFromCM(make_pair(k1, k2)));
 
   // LorentzMomentum phh = Rh*Phold;
 
   return Rh;
 
   //  return Utilities::transformToMomentum(Phold, Ph);
 
 }
 
 void PartonExtractor::
 transformRemnants(LorentzMomentum & Ph, LorentzMomentum & Pr,
 		  const LorentzMomentum & k, const LorentzMomentum & P) const {
   // don't do this for very soft remnants, as
   // we may run into numerical troubles; threshold
   // needs to become a parameter at some point
   if ( Pr.vect().mag2()/k.vect().mag2() < 1e-10 &&
        sqr(Pr.e()/k.e()) < 1e-10 )
     return;
   TransverseMomentum pt = Pr;
   try {
     if ( Direction<0>::pos() )
       SimplePhaseSpace::CMS(Pr, Ph, (P + k).m2(), 1.0, 0.0);
     else
       SimplePhaseSpace::CMS(Ph, Pr, (k + P).m2(), 1.0, 0.0);
     LorentzRotation rpt;
     if ( sqr(Pr.z()) > ZERO ) rpt.rotateY(asin(pt.pt()/Pr.z()));
     rpt.rotateZ(pt.phi());
     rpt = Direction<0>::pos()?
       Utilities::getBoostFromCM(make_pair(P, k))*rpt:
       Utilities::getBoostFromCM(make_pair(k, P))*rpt;
     Ph.transform(rpt);
     Pr.transform(rpt);
-  } catch ( ImpossibleKinematics ) {}
+  } catch ( ImpossibleKinematics & e ) {}
 }
 
 
 double PartonExtractor::fullFn(const PBIPair & pbins, Energy2 scale,
 			       pair<bool,bool> noLastPDF) {
   if(pbins.first->bin()->pdfDim()<=1) pbins.first->scale(scale);
   if(pbins.second->bin()->pdfDim()<=1) pbins.second->scale(scale);
   return fullFn(*pbins.first,noLastPDF.first)*fullFn(*pbins.second,noLastPDF.second);
 }
 
 double PartonExtractor::fullFn(const PartonBinInstance & pb,
 			       bool noLastPDF) {
   if ( !pb.incoming() ) return 1.0;
   if (noLastPDF)
     return 
       fullFn(*pb.incoming(),false) * pb.jacobian() * 
       pb.remnantWeight() * exp(-pb.li());
   return fullFn(*pb.incoming(),false) * pb.jacobian() * pb.remnantWeight() *
     pb.pdf()->xfl(pb.particleData(), pb.partonData(), pb.scale(),
 		  pb.li(), pb.incoming()->scale());
 }
 
 void PartonExtractor::
 construct(const PBIPair & pbins, tStepPtr step) const {
   // if a long chain we need to break some mother/child relationships
   if(pbins.first->incoming()) {
     if(pbins.first->incoming()->incoming()) {
       if(!pbins.first->parton()->parents().empty()) {
 	tParticleVector parents=pbins.first->parton()->parents();
 	tPPtr parton = pbins.first->parton();
 	for(unsigned int ix=0;ix<parents.size();++ix) parents[ix]->abandonChild(parton);
       }
     }
   }
   if(pbins.second->incoming()) {
     if(pbins.second->incoming()->incoming()) {
       if(!pbins.second->parton()->parents().empty()) {
 	tParticleVector parents=pbins.second->parton()->parents();
 	tPPtr parton = pbins.second->parton();
 	for(unsigned int ix=0;ix<parents.size();++ix) parents[ix]->abandonChild(parton);
       }
     }
   }
   Direction<0> dir(true);
   construct(*pbins.first, step);
   dir.reverse();
   construct(*pbins.second, step);
 }
 
 void PartonExtractor::
 construct(PartonBinInstance & pb, tStepPtr step, bool boost) const {
   if ( !pb.incoming() ) return;
   if ( boost ) pb.remnantHandler()->boostRemnants(pb);
   if ( pb.incoming()->incoming() ) {
     step->insertIntermediate(pb.particle(),pb.incoming()->particle(),pb.parton());
   }
   tPVector rem(pb.remnants().begin(), pb.remnants().end());
   if ( !step->addDecayProduct(pb.particle(), rem.begin(), rem.end(), false) )
     {}
   colourConnect(pb.particle(), pb.parton(), rem);
   construct(*pb.incoming(), step);
 }
 
 PBIPair PartonExtractor::newRemnants(tPPair oldp, tPPair newp, tStepPtr step) {
   PBIPair pb;
   Direction<0> dir(true);
   pb.first = newRemnants(partonBinInstance(oldp.first),
 			 newp.first, newp.second->momentum());
   dir.reverse();
   pb.second = newRemnants(partonBinInstance(oldp.second),
 			  newp.second, newp.first->momentum());
   addNewRemnants(partonBinInstance(oldp.first), pb.first, step);
   addNewRemnants(partonBinInstance(oldp.second), pb.second, step);
   return pb;
 }
 
 PBIPtr PartonExtractor::
 newRemnants(tPBIPtr oldpb, tPPtr newp, const LorentzMomentum & k) {
   if ( ! oldpb || !oldpb->incoming() ) return oldpb;
   Energy2 shat = (k + newp->momentum()).m2();
 
   // Loop over all possible PartonBin sisters to find the one
   // corresponding to the newly extracted parton.
   const PartonBin::PBVector & sisters = oldpb->incoming()->bin()->outgoing();
   for ( int i = 0, N = sisters.size(); i < N; ++i )
     if ( sisters[i]->parton() == newp->dataPtr() ) {
       // Setup necessary info in new PartonBinInstance object.
       PBIPtr newpb = new_ptr(PartonBinInstance(sisters[i], oldpb->incoming()));
       newpb->particle(oldpb->particle());
       newpb->parton(newp);
       newpb->li(log(oldpb->particle()->momentum().dirPlus()/
 		    newp->momentum().dirPlus()));
       newpb->l(oldpb->l() - oldpb->li() + newpb->li());
       Energy2 sc = -newp->scale();
       newpb->scale(newp->scale());
       if ( oldpb->incoming()->incoming() )
 	sc = -newpb->particle()->momentum().m2();
       // Now we can construct the new remnants.
       newpb->remnantWeight(1.0);
       if ( !newpb->remnantHandler()->
 	   recreateRemnants(*newpb, oldpb->parton(), newp, newpb->li(),
 			    sc, shat, newpb->particle()->momentum()) )
 	throw Veto();
       if ( newpb->remnantWeight() <= 0.0 ) throw Veto();
       return newpb;
     }
   throw Veto();
 }
 
 void PartonExtractor::
 addNewRemnants(tPBIPtr oldpb, tPBIPtr newpb, tStepPtr step) {
   if ( oldpb == newpb ) return;
   if ( oldpb->parton() != newpb->parton() ) {
     step->removeDecayProduct(newpb->particle(), oldpb->parton());
     if ( !step->addDecayProduct(newpb->particle(), newpb->parton()) )
       throw Veto();
   }
   tPVector rem(newpb->remnants().begin(), newpb->remnants().end());
   colourConnect(newpb->particle(), newpb->parton(), rem);
   partonBinInstances()[newpb->parton()] = newpb;
   if ( !step->addDecayProduct(oldpb->remnants().begin(),
 			      oldpb->remnants().end(),
  			      rem.begin(), rem.end()) )
     throw Veto();
 }
 
 void PartonExtractor::persistentOutput(PersistentOStream & os) const {
   os << theLastXComb << theSpecialDensities << theNoPDF << theMaxTries
      << flatSHatY << theFirstPDF << theSecondPDF;
 }
 
 void PartonExtractor::persistentInput(PersistentIStream & is, int) {
   is >> theLastXComb >> theSpecialDensities >> theNoPDF >> theMaxTries
      >> flatSHatY >> theFirstPDF >> theSecondPDF;
 }
 
 ClassDescription<PartonExtractor> PartonExtractor::initPartonExtractor;
 
 void PartonExtractor::Init() {
 
   static ClassDocumentation<PartonExtractor> documentation
     ("There is no documentation for the ThePEG::PartonExtractor class");
 
   static RefVector<PartonExtractor,PDFBase> interfaceSpecialDensities
     ("SpecialDensities",
      "A list of parton density objects to be used for incoming particles "
      "overriding possible densities given for particles of the "
      "BeamParticleData class.",
      &PartonExtractor::theSpecialDensities, 0, false, false, true, false);
 
   static Reference<PartonExtractor,PDFBase> interfaceNoPDF
     ("NoPDF",
      "A fixed reference to a NoPDF object to be used for particles without "
      "substructure.",
      &PartonExtractor::theNoPDF, true, true, true, false);
 
   static Parameter<PartonExtractor,int> interfaceMaxTries
     ("MaxTries",
      "The maximum number of attempts allowed when trying to generate "
      "remnants.",
      &PartonExtractor::theMaxTries, 100, 1, 1000, false, false, true);
 
   static Switch<PartonExtractor,bool> interfaceFlatSHatY
     ("FlatSHatY",
      "The possibility to override the l-generation in the PDFs and generate "
      "a flat distribution in \\f$\\log(\\hat{s})\\f$ and \\f$y\\f$. This only "
      "applies if the parton densities do not have poles in \\f$x=1\\f$.",
      &PartonExtractor::flatSHatY, false, false, false);
 
   static SwitchOption interfaceFlatSHatY0
     (interfaceFlatSHatY,
      "Off", "Use the l-generation defined by the PDFs", false);
 
   static SwitchOption interfaceFlatSHatY1
     (interfaceFlatSHatY,
      "On", "Generate flat rapidity and \\f$\\log(\\hat{s})\\f$", true);
 
   static SwitchOption interfaceFlatSHatNo
     (interfaceFlatSHatY,
      "No", "Use the l-generation defined by the PDFs", false);
 
   static SwitchOption interfaceFlatSHatYes
     (interfaceFlatSHatY,
      "Yes", "Generate flat rapidity and \\f$\\log(\\hat{s})\\f$", true);
 
   static Reference<PartonExtractor,PDFBase> interfaceFirstPDF
     ("FirstPDF",
      "PDF to override the default PDF for the first beam particle",
      &PartonExtractor::theFirstPDF, false, false, true, true, false);
 
   static Reference<PartonExtractor,PDFBase> interfaceSecondPDF
     ("SecondPDF",
      "PDF to override the default PDF for the second beam particle",
      &PartonExtractor::theSecondPDF, false, false, true, true, false);
 
 }
 
 RemColException::RemColException(const PartonExtractor & pe) {
   theMessage << "Parton extractor '" << pe.name() << "' failed to connect "
 	     << "the colours of the outgoing partons and the remnants.";
   severity(maybeabort);
 }
   
 void PartonExtractor::dofinish() {
   partonBinInstances().clear();
   HandlerBase::dofinish();
 }
diff --git a/PDT/DecayMode.cc b/PDT/DecayMode.cc
--- a/PDT/DecayMode.cc
+++ b/PDT/DecayMode.cc
@@ -1,691 +1,692 @@
 // -*- C++ -*-
 //
 // DecayMode.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 DecayMode class.
 //
 
 #include "DecayMode.h"
 #include "DecayMode.xh"
 #include "ThePEG/Repository/Repository.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Utilities/Rebinder.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/EnumIO.h"
 
 using namespace ThePEG;
 
 DecayMode::DecayMode()
   : theBrat(0.0), isOn(false) {}
 
 DecayMode::DecayMode(tPDPtr newParticle, double newBrat, bool newOn)
   : theBrat(newBrat), isOn(newOn), theParent(newParticle) {}
 
 DecayMode::DecayMode(const DecayMode & dm)
   : Interfaced(dm), theTag(dm.theTag), theBrat(dm.theBrat), isOn(dm.isOn),
     theParent(dm.theParent), theProducts(dm.theProducts),
     theOrderedProducts(dm.theOrderedProducts),
     theCascadeProducts(dm.theCascadeProducts), theMatchers(dm.theMatchers),
     theWildMatcher(dm.theWildMatcher), theExcluded(dm.theExcluded),
     theOverlap(dm.theOverlap), theDecayer(dm.theDecayer),
     theAntiPartner(dm.theAntiPartner), theLinks(dm.theLinks) {}
 
 DecayMode::~DecayMode() {}
 
 DMPtr DecayMode::Create(tPDPtr newParent, double newBrat, bool newOn) {
   DMPtr dm = new_ptr(DecayMode(newParent, newBrat, newOn));
   Repository::Register(dm, newParent->fullName() + "/NEWMODE");
   if ( !newParent->CC() ) return dm;
   DMPtr adm = new_ptr(DecayMode(newParent->CC(), newBrat, newOn));
   Repository::Register(adm, newParent->CC()->fullName() + "/NEWMODE");
   dm->theAntiPartner = adm;
   adm->theAntiPartner = dm;
   return dm;
 }
 
 void DecayMode::readSetup(istream & is) {
   string decnam;
   is >> theBrat >> ienum(isOn) >> decnam;
   if ( decnam.empty() ) return;
   BaseRepository::DirectoryAppend(decnam);
   setDecayer(BaseRepository::GetObject<tDecayerPtr>(decnam));	       
 }
 
 IBPtr DecayMode::clone() const {
   return dmclone();
 }
 
 DMPtr DecayMode::dmclone() const {
   return new_ptr(*this);
 }
 
 IBPtr DecayMode::fullclone() const {
   DMPtr dm = dmclone();
   Repository::Register(dm);
   if ( !CC() ) return dm;
   DMPtr adm = CC()->dmclone();
   Repository::Register(adm);
   dm->theAntiPartner = adm;
   adm->theAntiPartner = dm;
   return dm;
 }
 
 void DecayMode::doupdate() {
   Interfaced::doupdate();
   bool redo = touched();
   UpdateChecker::check(decayer(), redo);
   if ( !redo ) return;
   if ( theBrat > 0.0 && isOn && !decayer()->accept(*this) )
     throw DecModNoAccept(tag(), decayer()->name());
 }
 
 void DecayMode::rebind(const TranslationMap & trans) {
   try {
     theParent = trans.alwaysTranslate(theParent);
     ParticleMSet newProducts;
     trans.alwaysTranslate(inserter(newProducts),
 			  products().begin(), products().end());
     products().swap(newProducts);
     tPDVector newOrdered;
     trans.alwaysTranslate(inserter(newOrdered),
 			  orderedProducts().begin(), orderedProducts().end());
     theOrderedProducts.swap(newOrdered);
     ModeMSet newCasc;
     trans.alwaysTranslate(inserter(newCasc),
 			  cascadeProducts().begin(), cascadeProducts().end());
     cascadeProducts().swap(newCasc);
     MatcherMSet newMatchers;
     trans.alwaysTranslate(inserter(newMatchers),
 			  productMatchers().begin(), productMatchers().end());
     productMatchers().swap(newMatchers);
     wildProductMatcher() = trans.alwaysTranslate(wildProductMatcher());
     ParticleMSet newExclude;
     trans.alwaysTranslate(inserter(newExclude),
 			  excluded().begin(), excluded().end());
     excluded().swap(newExclude);
     theAntiPartner = trans.alwaysTranslate(CC());
 
     for ( int i = 0, N = theLinks.size(); i < N; ++i ) {
       theLinks[i].first = trans.alwaysTranslate(theLinks[i].first);
       theLinks[i].second = trans.alwaysTranslate(theLinks[i].second);
     }
   }
   catch (IBPtr ip) {
     throw DecModRebind(name(), ip->name());
   }
   catch (...) {
     throw DecModRebind(name(), "<unknown>");
   }
 }
 
 IVector DecayMode::getReferences() {
   IVector ret;
   ret.push_back(theParent);
   ret.insert(ret.end(), products().begin(), products().end());
   ret.insert(ret.end(), cascadeProducts().begin(), cascadeProducts().end());
   ret.insert(ret.end(), productMatchers().begin(), productMatchers().end());
   if ( wildProductMatcher() ) ret.push_back(wildProductMatcher());
   ret.insert(ret.end(), excluded().begin(), excluded().end());
   if ( CC() ) ret.push_back(CC());
   return ret;
 }
 
 bool DecayMode::addOverlap(tcDMPtr d) {
   bool inc = includes(*d);
   if ( !inc ) return false;
   if ( find(theOverlap.begin(), theOverlap.end(), d) != theOverlap.end() )
     return true;
   theOverlap.push_back(d);
   return true;
 }
 
 void DecayMode::resetOverlap() {
   theOverlap.clear();
 }
 
 bool DecayMode::
 compareId(const ParticleMSet & s1, const ParticleMSet & si2) const {
   if ( generator() ) return s1 == si2;
   ParticleMSet s2 = si2;
   for ( ParticleMSet::const_iterator p1 = s1.begin(); p1 != s1.end(); ++p1 ) {
     ParticleMSet::const_iterator p2 = s2.begin();
     while ( p2 != s2.end() && (**p2).id() != (**p1).id() ) ++p2;
     if ( p2 == s2.end() ) return false;
     s2.erase(p2);
   }
   return s2.empty();
 }
 
 ParticleMSet::const_iterator
 DecayMode::findId(const ParticleMSet & s, const ParticleData & p) const {
   for ( ParticleMSet::const_iterator pit = s.begin(); pit != s.end(); ++pit )
     if ( (**pit).id() == p.id() ) return pit;
   return s.end();
 }
 
 struct IdCmp {
   bool operator()(tcPDPtr p1, tcPDPtr p2) { return p1->id() == p2->id(); }
 };
 
 bool DecayMode::includes(const DecayMode & d) const {
   // Fast check for ordinary decay modes.
   if (  cascadeProducts().empty() && productMatchers().empty() &&
 	excluded().empty() &&  !wildProductMatcher() &&
 	d.cascadeProducts().empty() &&  d.productMatchers().empty() &&
 	d.excluded().empty() &&  !d.wildProductMatcher() ) {
     if ( links().size() != d.links().size() ) return false;
     if ( !compareId(products(), d.products()) ) return false;
     LinkVector dlinks = d.links();
     for ( int i = 0, N = links().size(); i < N; ++i ) {
       for ( int j = 0, M = dlinks.size(); j < M; ++j ) {
 	if ( ( links()[i].first->id() == dlinks[j].first->id() &&
 	       links()[i].second->id() == dlinks[j].second->id() ) ||
 	     ( links()[i].first->id() == dlinks[j].second->id() &&
 	       links()[i].second->id() == dlinks[j].first->id() ) ) {
 	  dlinks.erase(dlinks.begin() + j);
 	  break;
 	}
       }
       return false;
     }
     return dlinks.empty();
   }
 
   // First check that none of the excluded products in this are
   // present in the other.
   ParticleMSet::const_iterator pit;
   for ( pit = excluded().begin(); pit != excluded().end(); ++pit ) {
     if ( findId(d.products(), **pit ) != d.products().end() ) return false;
     for ( ModeMSet::const_iterator mit = d.cascadeProducts().begin();
 	  mit != d.cascadeProducts().end(); ++mit )
       if ( (**pit).id() == (**mit).parent()->id() ) return false;
   }
 
   // Check that all cascade decays in this overlaps with one in the
   // other. Save the ones that are left
   ModeMSet cascleft = d.cascadeProducts();
   for ( ModeMSet::iterator mit = cascadeProducts().begin();
 	mit != cascadeProducts().end(); ++mit ) {
     ModeMSet::iterator mit2 = cascleft.begin();
     while ( mit2 != cascleft.end() && !(**mit).includes(**mit2) ) ++mit2;
     if ( mit2 == cascleft.end() ) return false;
   }
 
 
   // Check that all cascade product parents in the other matches
   // something in this. Otherwise expand the cascade product.
   ParticleMSet partleft = d.products();
   MatcherMSet matchleft = d.productMatchers();
   ParticleMSet excludeleft = d.excluded();
   MatcherMSet wildleft;
   if ( d.wildProductMatcher() ) wildleft.insert(wildProductMatcher());
   ParticleMSet part = products();
   MatcherMSet match = productMatchers();
   while ( cascleft.size() ) {
     ModeMSet::iterator cdmit = cascleft.begin();
     cDMPtr cdm = *cdmit;
     ParticleMSet::iterator pit = findId(part, *(cdm->parent()));
     if ( pit != part.end() ) {
       cascleft.erase(cdmit);
       part.erase(pit);
     } else {
       MatcherMSet::iterator mit = match.begin();
       while ( mit != match.end() && !(**mit).matches(*(cdm->parent())) ) ++mit;
       if ( mit != match.end() ) {
 	cascleft.erase(cdmit);
 	match.erase(mit);
       } else {
 	if ( wildProductMatcher() &&
 	     wildProductMatcher()->matches(*(cdm->parent())) ) {
 	  cascleft.erase(cdmit);
 	} else {
 	  cascleft.erase(cdmit);
 	  partleft.insert(cdm->products().begin(), cdm->products().end());
 	  matchleft.insert(cdm->productMatchers().begin(),
 			   cdm->productMatchers().end());
 	  if ( cdm->wildProductMatcher() )
 	    wildleft.insert(cdm->wildProductMatcher());
 	  excludeleft.insert(cdm->excluded().begin(), cdm->excluded().end());
 	  cascleft.insert(cdm->cascadeProducts().begin(),
 			  cdm->cascadeProducts().end());
 	}
       }
     }
   }
 
   // Check that all excluded left in the other are absent in this.
   if ( find_first_of(excludeleft.begin(), excludeleft.end(), part.begin(),
 		     part.end(), IdCmp()) != excludeleft.end() ) return false;
 
   // Now all particles and matches left in this must match something
   // in the other.
   pit = part.begin();
   while ( pit != part.end() ) {
     ParticleMSet::iterator pit2 = findId(partleft, **pit++);
     if ( pit2 == partleft.end() ) return false;
     partleft.erase(pit2);
   }
   MatcherMSet::const_iterator pmit = match.begin();
   while ( pmit != match.end() ) {
     ParticleMSet::iterator pit2 = partleft.begin();
     while ( pit2 != partleft.end() && ! (**pmit).matches(**pit2) ) ++pit2;
     if ( pit2 != partleft.end() ) {
       partleft.erase(pit2);
     } else {
       MatcherMSet::iterator pmit2 = matchleft.begin();
       while ( pmit2 != matchleft.end() && ! (**pmit).matches(**pmit2) ) ++pmit2;
       if ( pmit2 != matchleft.end() ) {
 	matchleft.erase(pmit2);
       } else
 	return false;
     }
   }
 
   // Now all particles and matchers left in the other must be matched
   // by the wild match in this.
   if ( wildProductMatcher() ) {
     pit = partleft.begin();
     while ( pit != partleft.end() )
       if ( !(wildProductMatcher()->matches(**pit++)) ) return false;
     pmit = matchleft.begin();
     while (pmit != matchleft.end() )
       if ( !(wildProductMatcher()->matches(**pmit++)) ) return false;
     pmit = wildleft.begin();
     while (pmit != wildleft.end() )
       if ( !(wildProductMatcher()->matches(**pmit++)) ) return false;
   } else
     return partleft.empty() && matchleft.empty() && wildleft.empty();
   return true;
 }
 
 DMPtr DecayMode::clone(tPDPtr pd) const {
   DMPtr dm = dmclone();
   dm->theParent = pd;
   Repository::Register(dm, pd->fullName() + "/" + dm->name());
   if ( !theDecayer || !theDecayer->accept(*dm) ) dm->isOn = false;
   if ( pd->CC() ) {
     DMPtr adm = CC()? CC()->dmclone(): dmclone();
     adm->theParent = pd->CC();
     Repository::Register(adm, pd->CC()->fullName() + "/" + adm->name());
     dm->theAntiPartner = adm;
     adm->theAntiPartner = dm;
     if ( !adm->theDecayer->accept(*adm) ) adm->isOn = false;
   } else
     dm->theAntiPartner = DMPtr();
   return dm;
 }
 
 void DecayMode::synchronize() {
   if ( !CC() ) return;
   theBrat = CC()->theBrat;
   isOn = CC()->isOn;
   theDecayer = CC()->theDecayer;
 }
 
 struct ParticleOrdering {
   bool operator()(tcPDPtr p1, tcPDPtr p2) {
     return abs(p1->id()) > abs(p2->id()) ||
       ( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
       ( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
   }
 };
 
 struct ModeOrdering {
   bool operator()(tcDMPtr d1, tcDMPtr d2) {
     ParticleOrdering ord;
     return ord(d1->parent(), d2->parent()) ||
       ( !ord(d2->parent(), d1->parent()) &&
 	( d1->tag() < d2->tag() ||
 	  ( d1->tag() == d2->tag() && d1->fullName() < d2->fullName() ) ) );
   }
 };
 
 struct MatcherOrdering {
   bool operator()(tcPMPtr m1, tcPMPtr m2) {
     return m1->name() < m2->name() ||
       ( m1->name() == m2->name() && m1->fullName() < m2->fullName() );
   }
 };
 
 PVector DecayMode::produceProducts() const {
   PVector ret;
   for ( int i = 0, N = orderedProducts().size(); i < N; ++i )
     ret.push_back(orderedProducts()[i]->produceParticle());
   return ret;
 }
 
 string DecayMode::makeTag() const {
   string ret;
   typedef multiset<tcPDPtr,ParticleOrdering> OrderedParticles;
   typedef multiset<tcPMPtr,MatcherOrdering> OrderedMatchers;
   typedef multiset<tcDMPtr,ModeOrdering> OrderedModes;
   LinkVector dlinks = links();
 
   ret = theParent->PDGName() + "->";
   if ( dlinks.empty() ) {
     OrderedParticles prod(products().begin(), products().end());
     for (OrderedParticles::iterator pit = prod.begin();
 	 pit != prod.end(); ++pit )
       ret += (**pit).PDGName() + ",";
   } else {
     unsigned int dl = 0;
     for ( int i = 0, N = orderedProducts().size(); i < N; ++i ) {
       if ( dl < dlinks.size() && orderedProducts()[i] == dlinks[dl].first ) {
 	ret +=  orderedProducts()[i]->PDGName() + "=";
 	++dl;
       } else 
 	ret +=  orderedProducts()[i]->PDGName() + ",";
     }
   }
 
   OrderedModes casc(cascadeProducts().begin(), cascadeProducts().end());
   for ( OrderedModes::iterator dmit = casc.begin();dmit != casc.end(); ++dmit )
     ret += "[" + (**dmit).tag() + "],";
   OrderedMatchers match(productMatchers().begin(), productMatchers().end());
   for ( OrderedMatchers::iterator mit = match.begin();
 	mit != match.end(); ++mit ) ret += "?" +(**mit).name() + ",";
   if ( theWildMatcher )
     ret += "*" + theWildMatcher->name() + ",";
   OrderedParticles ex(excluded().begin(), excluded().end());
   for ( OrderedParticles::iterator pit = ex.begin(); pit != ex.end(); ++pit )
     ret += "!" + (**pit).PDGName() + ",";
   ret[ret.size()-1] = ';';
   return ret;
 }
 
 void DecayMode::brat(double newBrat) {
   theBrat = newBrat;
   if ( theBrat <= 0.0 ) switchOff();
   if ( CC() && parent()->synchronized() ) CC()->theBrat = newBrat;
 }
 
 double DecayMode::brat() const {
   return isOn? theDecayer->brat(*this, *theParent, theBrat): 0.0;
 }
 
 double DecayMode::brat(const Particle & p) const {
   return isOn && p.dataPtr() == parent()?
     theDecayer->brat(*this, p, theBrat): 0.0;
 }
 
 void DecayMode::switchOn() {
   isOn = true;
   if ( CC() && parent()->synchronized() ) CC()->isOn = true;
 }
 
 void DecayMode::switchOff() {
   isOn = false;
   if ( CC() && parent()->synchronized() ) CC()->isOn = false;
 }
 
 void DecayMode::addProduct(tPDPtr pd) {
   products().insert(pd);
   theOrderedProducts.push_back(pd);
   if ( CC() ) {
     CC()->products().insert(pd->CC()? pd->CC(): pd);
     CC()->theOrderedProducts.push_back(pd->CC()? pd->CC(): pd);
   }
   resetTag();
 }
 
 void DecayMode::addLink(tPDPtr a, tPDPtr b) {
   theLinks.push_back(make_pair(a, b));
   if ( CC() ) CC()->theLinks.push_back(make_pair(a->CC()? a->CC(): a,
 						 b->CC()? b->CC(): b));
   resetTag();
 }
 
 void DecayMode::addCascadeProduct(tDMPtr dm) {
   cascadeProducts().insert(dm);
   if ( CC() ) CC()->cascadeProducts().insert(dm->CC()? dm->CC(): dm);
   resetTag();
 }
 
 void DecayMode::addProductMatcher( tPMPtr pm) {
   productMatchers().insert(pm);
   if ( CC() ) CC()->productMatchers().insert(pm->CC()? pm->CC(): pm);
   resetTag();
 }
 
 void DecayMode::setWildMatcher(tPMPtr pm) {
   wildProductMatcher() = pm;
   if ( CC() ) CC()->wildProductMatcher() = pm->CC()? pm->CC(): pm;
   resetTag();
 }
 
 void DecayMode::addExcluded(tPDPtr pd) {
   excluded().insert(pd);
   if ( CC() ) CC()->excluded().insert(pd->CC()? pd->CC(): pd);
   resetTag();
 }
 
 void DecayMode::decayer(tDecayerPtr dec) {
   if ( !dec || !dec->accept(*this) )
     throw DecModSetupNoAccept(tag(), dec->name());
   if ( CC() && parent()->synchronized() ) {
     if ( !dec->accept(*CC()) )
       throw DecModSetupNoAccept(CC()->tag(), dec->name());
     CC()->theDecayer = dec;
   }
   theDecayer = dec;
 }
 
 DMPtr DecayMode::constructDecayMode(string & tag, vector<DMPtr> * save) {
   DMPtr rdm;
   DMPtr adm;
   int level = 0;
   string::size_type end = 0;
   while ( end < tag.size() && ( tag[end] != ']' || level ) ) {
     switch ( tag[end++] ) {
     case '[':
       ++level;
       break;
     case ']':
       --level;
       break;
     }
   }
   string::size_type next = tag.find("->");
   if ( next == string::npos ) return rdm;
   if ( tag.find(';') == string::npos ) return rdm;
   tPDPtr pd = Repository::findParticle(tag.substr(0,next));
   if ( !pd ) return rdm;
 
   rdm = ptr_new<DMPtr>();
   rdm->parent(pd);
   if ( pd->CC() ) {
     adm = ptr_new<DMPtr>();
     adm->parent(pd->CC());
     rdm->theAntiPartner = adm;
     adm->theAntiPartner = rdm;
   }
   bool error = false;
   tag = tag.substr(next+2);
   tPDPtr lastprod;
   bool dolink = false;
   do {
     switch ( tag[0] ) {
     case '[':
       {
 	tag = tag.substr(1);
 	DMPtr cdm = constructDecayMode(tag, save);
 	if ( save ) save->push_back(cdm);
 	if ( cdm ) rdm->addCascadeProduct(cdm);
 	else error = true;
       } break;
     case '=':
       dolink = true;
+      [[fallthrough]]; 
     case ',':
     case ']':
       tag = tag.substr(1);
       break;
     case '?':
       {
 	next = min(tag.find(','), tag.find(';'));
 	tPMPtr pm = Repository::findMatcher(tag.substr(1,next-1));
 	if ( pm ) rdm->addProductMatcher(pm);
 	else error = true;
 	tag = tag.substr(next);
       } break;
     case '!':
       {
 	next = min(tag.find(','), tag.find(';'));
 	tPDPtr pd = Repository::findParticle(tag.substr(1,next-1));
 	if ( pd ) rdm->addExcluded(pd);
 	else error = true;
 	tag = tag.substr(next);
       } break;
     case '*':
       {
 	next = min(tag.find(','), tag.find(';'));
 	tPMPtr pm = Repository::findMatcher(tag.substr(1,next-1));
 	if ( pm ) rdm->setWildMatcher(pm);
 	else error = true;
 	tag = tag.substr(next);
       } break;
     default:
       {
 	next = min(tag.find('='), min(tag.find(','), tag.find(';')));
 	tPDPtr pdp = Repository::findParticle(tag.substr(0,next));
 	if ( pdp ) rdm->addProduct(pdp);
 	else error = true;
 	tag = tag.substr(next);
 	if ( dolink && lastprod ) {
 	  rdm->addLink(lastprod, pdp);
 	  dolink = false;
 	}
 	lastprod = pdp;
       } break;
     }
   } while ( tag[0] != ';' && tag.size() );
   if ( tag[0] != ';' || error ) {
     return DMPtr();
   }
 
   tag = tag.substr(1);
   
   for ( DecaySet::const_iterator dit = pd->decayModes().begin();
 	dit != pd->decayModes().end(); ++dit )
     if ( (**dit).tag() == rdm->tag() ) return *dit;
 
   if ( save ) {
     save->push_back(rdm);
     save->push_back(adm);
   } else {
     pd->addDecayMode(rdm);
     Repository::Register(rdm, pd->fullName() + "/" + rdm->tag());
     if ( adm )
       Repository::Register(adm, pd->CC()->fullName() + "/" + adm->tag());
   }
   return rdm;
 }
 
 void DecayMode::persistentOutput(PersistentOStream & os) const {
   multiset<tcPDPtr,ParticleOrdering> prod(products().begin(), products().end());
   multiset<tcDMPtr,ModeOrdering>
     casc(cascadeProducts().begin(), cascadeProducts().end());
   multiset<tcPMPtr,MatcherOrdering>
     match(productMatchers().begin(), productMatchers().end());
   multiset<tcPDPtr,ParticleOrdering>  ex(excluded().begin(), excluded().end());
   multiset<tcDMPtr,ModeOrdering> ovlap(overlap().begin(), overlap().end());
 
   os << theTag << theBrat << isOn << theParent << prod << theOrderedProducts
      << casc << match << theWildMatcher << ex << ovlap << theDecayer
      << theAntiPartner << theLinks;
 }
 
 void DecayMode::persistentInput(PersistentIStream & is, int) {
   is >> theTag >> theBrat >> isOn >> theParent >> theProducts
      >> theOrderedProducts >> theCascadeProducts >> theMatchers
      >> theWildMatcher >> theExcluded >> theOverlap >> theDecayer
      >> theAntiPartner >> theLinks;
 }
 
 ClassDescription<DecayMode> DecayMode::initDecayMode;
 
 void DecayMode::setOn(long i) {
   isOn = i;
 }
 
 long DecayMode::getOn() const {
   return isOn;
 }
 
 void DecayMode::setDecayer(DecayerPtr dp) {
   decayer(dp);
 }
 
 void DecayMode::Init() {
 
   static ClassDocumentation<DecayMode> documentation
     ("Represents a specific decay channel of a particle.");
 
   static Parameter<DecayMode,double> interfaceBrat
     ("BranchingRatio",
      "The branching fraction for this decay mode. Note that if the sum of "
      "branching ratios for one particle is always renormalized to 1. Also, "
      "the decaying particle may change this branching ratio if it has a "
      "ThePEG::WidthGenerator object assigned to it. ",
      &DecayMode::theBrat, 0.0, 0.0, 1.0, false, false, true);
   interfaceBrat.setHasDefault(false);
 
   static Switch<DecayMode> interfaceOn
     ("OnOff",
      "Indicates if the decay mode is switched on or off.",
      0, 0, false, false, &DecayMode::setOn, &DecayMode::getOn);
   static SwitchOption interfaceOnYes
     (interfaceOn, "On", "The decay channel is switched on.", 1);
   static SwitchOption interfaceOnNo
     (interfaceOn, "Off", "The decay channel is switched off.", 0);
   interfaceOn.setHasDefault(false);
 
   static Switch<DecayMode> interfaceActive
     ("Active",
      "Indicates if the decay mode is switched on or off.",
      0, 0, false, false, &DecayMode::setOn, &DecayMode::getOn);
   static SwitchOption interfaceActiveYes
     (interfaceActive, "Yes", "The decay channel is switched on.", 1);
   static SwitchOption interfaceActiveNo
     (interfaceActive, "No", "The decay channel is switched off.", 0);
   interfaceActive.setHasDefault(false);
 
   static Reference<DecayMode,Decayer> interfaceDecayer
     ("Decayer",
      "The ThePEG::Decayer object responsible for performing this decay.",
      &DecayMode::theDecayer, false, false, true, false,
      &DecayMode::setDecayer);
 
   interfaceBrat.rank(10);
   interfaceDecayer.rank(9);
   interfaceOn.rank(8);
   interfaceActive.rank(8);
 
 }
 
 DecModNoAccept::DecModNoAccept(string tag, string dec) {
   theMessage << "The Decayer '" << dec << "' is not capable to "
 	     << "perform the decay in the DecayMode '" << tag << "'.";
   severity(warning);
 }
 
 DecModSetupNoAccept::DecModSetupNoAccept(string tag, string dec) {
   theMessage << "The Decayer '" << dec << "' is not capable to "
 	     << "perform the decay in the DecayMode '" << tag << "'.";
   severity(warning);
 }
 
 DecModRebind::DecModRebind(string tag, string obj) {
   theMessage << "'Rebind' of DecayMode '" << tag << "' failed because "
 	     << "the object '" << obj << "' refered to lacked a translation.";
   severity(abortnow);
 }
diff --git a/PDT/ParticleData.cc b/PDT/ParticleData.cc
--- a/PDT/ParticleData.cc
+++ b/PDT/ParticleData.cc
@@ -1,1038 +1,1038 @@
 // -*- C++ -*-
 //
 // ParticleData.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 ParticleData class.
 //
 
 #include "ParticleData.h"
 #include "ParticleData.xh"
 #include "ThePEG/PDT/DecayMode.h"
 #include "ThePEG/Utilities/HoldFlag.h"
 #include "ThePEG/Utilities/Rebinder.h"
 #include "ThePEG/Utilities/StringUtils.h"
 #include "ThePEG/EventRecord/Particle.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/RefVector.h"
 #include "ThePEG/Interface/Command.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Repository/Repository.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Config/algorithm.h"
 #include "ThePEG/Utilities/Exception.h"
 #include "ThePEG/Utilities/EnumIO.h"
 #include "ThePEG/Repository/UseRandom.h"
 
 namespace ThePEG {
 
 ParticleData::ParticleData()
   : theId(0), thePDGName(""), theMass(-1.0*GeV), theWidth(-1.0*GeV),
     theHardProcessMass(-1.0*GeV), hardProcessMassSet(false),
     theHardProcessWidth(-1.0*GeV), hardProcessWidthSet(false),
     theWidthUpCut(-1.0*GeV), theWidthLoCut(-1.0*GeV), theCTau(-1.0*mm),
     theCharge(PDT::ChargeUnknown),
     theSpin(PDT::SpinUnknown), theColour(PDT::ColourUnknown), isStable(true),
     theVariableRatio(false), syncAnti(false), theDefMass(-1.0*GeV),
     theDefWidth(-1.0*GeV), theDefCut(-1.0*GeV), theDefCTau(-1.0*mm),
     theDefCharge(PDT::ChargeUnknown), theDefSpin(PDT::SpinUnknown),
     theDefColour(PDT::ColourUnknown) {}
 
 ParticleData::
 ParticleData(PID newId, const string & newPDGName)
   : theId(newId), thePDGName(newPDGName), theMass(-1.0*GeV), theWidth(-1.0*GeV),
     theHardProcessMass(-1.0*GeV), hardProcessMassSet(false),
     theHardProcessWidth(-1.0*GeV), hardProcessWidthSet(false),
     theWidthUpCut(-1.0*GeV), theWidthLoCut(-1.0*GeV), theCTau(-1.0*mm),
     theCharge(PDT::ChargeUnknown),
     theSpin(PDT::SpinUnknown), theColour(PDT::ColourUnknown), isStable(true),
     theVariableRatio(false), syncAnti(false), theDefMass(-1.0*GeV),
     theDefWidth(-1.0*GeV), theDefCut(-1.0*GeV), theDefCTau(-1.0*mm),
     theDefCharge(PDT::ChargeUnknown), theDefSpin(PDT::SpinUnknown),
     theDefColour(PDT::ColourUnknown) {}
 
 ParticleData::~ParticleData() {}
 
 PDPtr ParticleData::Create(PID newId, const string & newPDGName) {
   return new_ptr(ParticleData(newId, newPDGName));
 }
 
 PDPair ParticleData::
 Create(PID newId, const string & newPDGName, const string & newAntiPDGName) {
   PDPair pap;
   pap.first = new_ptr(ParticleData(newId, newPDGName));
   pap.second = new_ptr(ParticleData(-newId, newAntiPDGName));
   antiSetup(pap);
   return pap;
 }
 
 void ParticleData::readSetup(istream & is) {
   long id;
   is >> id >> thePDGName >> iunit(theDefMass, GeV) >> iunit(theDefWidth, GeV)
      >> iunit(theDefCut, GeV) >> iunit(theDefCTau, mm) >> ienum(theDefCharge)
      >> ienum(theDefColour) >> ienum(theDefSpin) >> ienum(isStable);
   theId = id;
   theMass = theDefMass;
   theWidth = theDefWidth;
   theWidthUpCut = theDefCut;
   theWidthLoCut = theDefCut;
   theCTau = theDefCTau;
   theCharge = theDefCharge;
   theColour = theDefColour;
   theSpin = theDefSpin;
   if ( PDGName() == "-" ) thePDGName = name();
   return;
 }
   
 void ParticleData::antiSetup(const PDPair & pap) {
   pap.first->theAntiPartner = pap.second;
   pap.second->theAntiPartner = pap.first;
   pap.first->syncAnti = pap.second->syncAnti = true;
 }
 
 PDPtr ParticleData::pdclone() const {
   return new_ptr(*this);
 }
 
 IBPtr ParticleData::clone() const {
   return pdclone();
 }
   
 IBPtr ParticleData::fullclone() const {
   PDPtr pd = pdclone();
   Repository::Register(pd);
   pd->theDecaySelector.clear();
   pd->theDecayModes.clear();
   pd->isStable = true;
   PDPtr apd;
   if ( CC() ) {
     apd = CC()->pdclone();
     Repository::Register(apd);
     apd->theDecaySelector.clear();
     apd->theDecayModes.clear();
     apd->isStable = true;
     pd->theAntiPartner = apd;
     apd->theAntiPartner = pd;
     pd->syncAnti = syncAnti;
     apd->syncAnti = CC()->syncAnti;
   }
   HoldFlag<> dosync(pd->syncAnti, true);
   for ( DecaySet::const_iterator it = theDecayModes.begin();
 	it != theDecayModes.end(); ++it )
     pd->addDecayMode(*it);
   return pd;
 }
 
 Energy ParticleData::mass(Energy mi) {
   theMass = mi;
   if ( synchronized() && CC() ) CC()->theMass = theMass;
   return theMass;
 }
 
 Energy ParticleData::width(Energy wi) {
   theWidth = wi;
   if ( synchronized() && CC() ) CC()->theWidth = theWidth;
   return theWidth;
 }
 
 Energy ParticleData::widthUpCut(Energy wci) {
   theWidthUpCut = wci;
   if ( synchronized() && CC() ) CC()->theWidthUpCut = theWidthUpCut;
   return theWidthUpCut;
 }
 
 Energy ParticleData::widthLoCut(Energy wci) {
   theWidthLoCut = wci;
   if ( synchronized() && CC() ) CC()->theWidthLoCut = theWidthLoCut;
   return theWidthLoCut;
 }
 
 Length ParticleData::cTau(Length ti) {
   theCTau = ti;
   if ( synchronized() && CC() ) CC()->theCTau = theCTau;
   return theCTau;
 }
 
 PDT::Charge ParticleData::iCharge(PDT::Charge ci) {
   theCharge = ci;
   if ( synchronized() && CC() ) CC()->theCharge = PDT::Charge(-ci);
   return theCharge;
 }
 
 PDT::Spin ParticleData::iSpin(PDT::Spin si) {
   theSpin = si;
   if ( synchronized() && CC() ) CC()->theSpin = si;
   return si;
 }
 
 PDT::Colour ParticleData::iColour(PDT::Colour ci) {
   theColour = ci;
   if ( synchronized() && CC() ) CC()->theColour = PDT::Colour(-ci);
   return theColour;
 }
 
 void ParticleData::stable(bool s) {
   isStable = s;
   if ( synchronized() && CC() ) CC()->isStable = s;
 }
 
 void ParticleData::synchronized(bool h) {
   syncAnti = h;
   if ( CC() ) CC()->syncAnti = h;
 }
 
 void ParticleData::variableRatio(bool varRatio) {
   theVariableRatio=varRatio;
 }
 
 void ParticleData::addDecayMode(tDMPtr dm) {
   if ( member(theDecayModes, dm) ) return;
   cPDPtr parent = dm->parent();
   if ( !parent ) parent = this;
   if ( parent != this ) {
     dm = dm->clone(this);
   }
   theDecayModes.insert(dm);
   theDecaySelector.insert(dm->brat(), dm);
   if ( CC() ) {
     if ( !synchronized() ) dm->CC()->switchOff();
     CC()->theDecayModes.insert(dm->CC());
     CC()->theDecaySelector.insert(dm->CC()->brat(), dm->CC());
   }
 }
 
 void ParticleData::removeDecayMode(tDMPtr dm) {
   theDecayModes.erase(theDecayModes.find(dm));
   if(theDecayModes.empty()) isStable = true;
   theDecaySelector.erase(dm);
   if ( !CC() ) return;
   CC()->theDecayModes.erase(dm->CC());
   if(CC()->theDecayModes.empty()) CC()->isStable = true;
   CC()->theDecaySelector.erase(dm->CC());
 }
 
 void ParticleData::synchronize() {
   if ( !CC() ) return;
   isStable = CC()->isStable;
   theMass = CC()->theMass;
   theHardProcessMass = CC()->theHardProcessMass;
   hardProcessMassSet = CC()->hardProcessMassSet;
   theWidth = CC()->theWidth;
   theHardProcessWidth = CC()->theHardProcessWidth;
   hardProcessWidthSet = CC()->hardProcessWidthSet;
   theWidthUpCut = CC()->theWidthUpCut;
   theWidthLoCut = CC()->theWidthLoCut;
   theCTau = CC()->theCTau;
   theCharge = PDT::Charge(-CC()->theCharge);
   theSpin = CC()->theSpin;
   theColour = PDT::antiColour(CC()->theColour);
   theMassGenerator = CC()->theMassGenerator;
   theWidthGenerator = CC()->theWidthGenerator;
   syncAnti = CC()->syncAnti;
   theDecaySelector.clear();
   for ( DecaySet::iterator it = theDecayModes.begin();
 	it != theDecayModes.end(); ++it ) {
     (*it)->synchronize();
     theDecaySelector.insert((*it)->brat(), *it);
   }
 }
 
 void ParticleData::doupdate() {
   Interfaced::doupdate();
   bool redo = touched();
   for_each(theDecayModes, UpdateChecker(redo));
   UpdateChecker::check(theMassGenerator, redo);
   UpdateChecker::check(theWidthGenerator, redo);
   if ( !redo ) return;
 
   theDecaySelector.clear();
   for ( DecaySet::const_iterator dit = theDecayModes.begin();
 	dit != theDecayModes.end(); ++dit ) {
     tDMPtr dm = *dit;
     dm->resetOverlap();
     for ( DecaySet::const_iterator dit2 = theDecayModes.begin();
 	  dit2 != theDecayModes.end(); ++dit2 )
       if ( dit2 != dit ) dm->addOverlap(dm);
     if ( dm->brat() > 0.0 ) theDecaySelector.insert(dm->brat(), dm);
   }
   if ( theMassGenerator && !theMassGenerator->accept(*this) )
     throw UpdateException();
   if ( theWidthGenerator &&
        !theWidthGenerator->accept(*this) )
     throw UpdateException();
   if ( theWidthGenerator ) theDecaySelector = theWidthGenerator->rate(*this);
   touch();
 }
 
 tDMPtr ParticleData::selectMode(Particle & p) const {
   if ( &(p.data()) != this ) return tDMPtr();
   try {
     if ( !theWidthGenerator || !theVariableRatio )
       return theDecaySelector.select(UseRandom::current());
     DecaySelector local;
     if ( theWidthGenerator )
       local = theWidthGenerator->rate(p);
     else
       for ( DecaySet::const_iterator mit = theDecayModes.begin();
 	    mit != theDecayModes.end(); ++mit  )
 	local.insert((*mit)->brat(p), *mit);
     return local.select(UseRandom::current());
   }
-  catch (range_error) {
+  catch (range_error & e) {
     return tDMPtr();
   }
 }
 
 void ParticleData::rebind(const TranslationMap & trans) {
   if ( CC() ) theAntiPartner = trans.translate(theAntiPartner);
   DecaySet newModes;
   DecaySelector newSelector;
   for ( DecaySet::iterator it = theDecayModes.begin();
 	it != theDecayModes.end(); ++it ) {
     DMPtr dm;
     dm = trans.translate(*it);
     if ( !dm ) throw RebindException();
     newModes.insert(dm);
     newSelector.insert(dm->brat(), dm);
   }
   theDecayModes.swap(newModes);
   theDecaySelector.swap(newSelector);
 }
 
 IVector ParticleData::getReferences() {
   IVector refs = Interfaced::getReferences();
   if ( CC() ) refs.push_back(CC());
   refs.insert(refs.end(), theDecayModes.begin(), theDecayModes.end());
   return refs;
 }
 
 
 
 void ParticleData::massGenerator(tMassGenPtr mg) {
   if ( mg && !mg->accept(*this) ) return;
   if ( mg && synchronized() && CC() && !mg->accept(*CC()) ) return;
   theMassGenerator = mg;
   if ( synchronized() && CC() ) CC()->theMassGenerator = mg;
 }
 
 void ParticleData::widthGenerator(tWidthGeneratorPtr newGen) {
   if ( newGen && !newGen->accept(*this) ) return;
   if ( newGen && synchronized() && CC() && !newGen->accept(*CC()) ) return;
   theWidthGenerator = newGen;
   if ( synchronized() && CC() ) CC()->theWidthGenerator = newGen;
 }
 
 Energy ParticleData::generateMass() const {
   return massGenerator()? massGenerator()->mass(*this): mass();
 }
 
 Energy ParticleData::generateWidth(Energy m) const {
   return widthGenerator()? widthGenerator()->width(*this, m): width();
 }
 
 Length ParticleData::generateLifeTime(Energy m, Energy w) const {
   return widthGenerator() ? 
     widthGenerator()->lifeTime(*this, m, w) :
     UseRandom::rndExp(cTau());
 }
 
 PPtr ParticleData::produceParticle(const Lorentz5Momentum & pp) const {
   PPtr p = new_ptr(Particle(this));
   p->set5Momentum(pp);
   return p;
 }
 
 PPtr ParticleData::produceParticle(const LorentzMomentum & pp) const {
   PPtr p(produceParticle(Lorentz5Momentum(pp)));
   return p;
 }
 
 PPtr ParticleData::produceParticle(const LorentzMomentum & pp, Energy m) const {
   PPtr p(produceParticle(Lorentz5Momentum(pp, m)));
   return p;
 }
 
 PPtr ParticleData::produceParticle(Energy m, const Momentum3 & pp) const {
   PPtr p(produceParticle(Lorentz5Momentum(m, pp)));
   return p;
 }
 
 PPtr ParticleData::produceParticle(const Momentum3 & pp) const {
   PPtr p(produceParticle(Lorentz5Momentum(generateMass(), pp)));
   return p;
 }
 
 PPtr ParticleData::
 produceParticle(Energy plus, Energy minus, Energy px, Energy py) const {
   PPtr p(produceParticle(LorentzMomentum(px, py, 0.5*(plus-minus),
 					 0.5*(plus+minus))));
   return p;
 }
 
 void ParticleData::setMass(Energy mi) {
   mass(mi);
 }
 
 void ParticleData::setHardProcessMass(Energy mi) {
   theHardProcessMass = mi;
   hardProcessMassSet = true;
   ParticleData * apd = CC().operator->();
   if ( synchronized() && apd ) {
     apd->theHardProcessMass = theHardProcessMass;
     apd->hardProcessMassSet = true;
   }
 }
 
 void ParticleData::setHardProcessWidth(Energy mi) {
   theHardProcessWidth = mi;
   hardProcessWidthSet = true;
   ParticleData * apd = CC().operator->();
   if ( synchronized() && apd ) {
     apd->theHardProcessWidth = theHardProcessWidth;
     apd->hardProcessWidthSet = true;
   }
 }
 
 string ParticleData::doUnsetHardProcessMass(string) {
   hardProcessMassSet = false;  
   theHardProcessMass = -1.*GeV;
   return "";
 }
 
 string ParticleData::doAdjustNominalMass(string) {
   if ( hardProcessMassSet )
     setMass(theHardProcessMass);
   return "";
 }
   
 string ParticleData::doUnsetHardProcessWidth(string) {
   hardProcessWidthSet = false;  
   theHardProcessWidth = -1.*GeV;
   return "";
 }
 
 Energy ParticleData::defMass() const {
   return theDefMass;
 }
 
 void ParticleData::setWidth(Energy wi) {
   width(wi);
 }
 
 Energy ParticleData::getWidth() const {
   return width();
 }
 
 Energy ParticleData::defWidth() const {
   return theDefWidth;
 }
 
 void ParticleData::setCut(Energy ci) {
   widthCut(ci);
 }
 
 Energy ParticleData::getCut() const {
   return (theWidthUpCut >= ZERO && theWidthLoCut >= ZERO)?
     max(theWidthUpCut, theWidthLoCut): min(theWidthUpCut, theWidthLoCut);
 }
 
 Energy ParticleData::defCut() const {
   return theDefCut;
 }
 
 void ParticleData::setUpCut(Energy ci) {
   widthUpCut(ci);
 }
 
 Energy ParticleData::getUpCut() const {
   return theWidthUpCut;
 }
 
 void ParticleData::setLoCut(Energy ci) {
   widthLoCut(ci);
 }
 
 Energy ParticleData::getLoCut() const {
   return theWidthLoCut;
 }
 
 void ParticleData::setCTau(Length ti) {
   cTau(ti);
 }
 
 Length ParticleData::getCTau() const {
   return cTau();
 }
 
 Length ParticleData::defCTau() const {
   return theDefCTau;
 }
 
 void ParticleData::setStable(long is) {
   stable(is);
 }
 
 long ParticleData::getStable() const {
   return stable();
 }
 
 void ParticleData::setSync(long is) {
   synchronized(is);
 }
 
 long ParticleData::getSync() const {
   return synchronized();
 }
 
 void ParticleData::setVariableRatio(long is) {
   variableRatio(is);
 }
 
 long ParticleData::getVariableRatio() const {
   return variableRatio();
 }
 
 string ParticleData::doSync(string) {
   synchronize();
   return "";
 }
 
 void ParticleData::setMassGenerator(MassGenPtr gi) {
   massGenerator(gi);
 }
 
 void ParticleData::setWidthGenerator(WidthGeneratorPtr wg) {
   widthGenerator(wg);
 }
 
 void ParticleData::setColour(long c) {
   theColour = PDT::Colour(c);
 }
 
 long ParticleData::getColour() const {
   return theColour;
 }
 
 long ParticleData::defColour() const {
   return theDefColour;
 }
 
 void ParticleData::setCharge(int c) {
   theCharge = PDT::Charge(c);
 }
 
 string ParticleData::ssetCharge(string arg) {
   istringstream is(arg);
   long i;
   if ( is >> i ) {
     theCharge = PDT::Charge(i);
     return "New charge is " + arg;
   }
   if ( arg == "unknown" )
     theCharge = PDT::ChargeUnknown;
   else if ( arg == "charged" )
     theCharge = PDT::Charged;
   else if ( arg == "positive" )
     theCharge = PDT::Positive;
   else if ( arg == "negative" )
     theCharge = PDT::Negative;
   else throw ParticleChargeCommand(*this, arg);
   return  "New charge is " + arg;
 }
 
 int ParticleData::getCharge() const {
   return theCharge;
 }
 
 int ParticleData::defCharge() const {
   return theDefCharge;
 }
 
 void ParticleData::setSpin(int s) {
   theSpin = PDT::Spin(s);
 }
 
 int ParticleData::getSpin() const {
   return theSpin;
 }
 
 int ParticleData::defSpin() const {
   return theDefSpin;
 }
 
 ClassDescription<ParticleData> ParticleData::initParticleData;
 
 struct ParticleOrdering {
   bool operator()(tcPDPtr p1, tcPDPtr p2) {
     return abs(p1->id()) > abs(p2->id()) ||
       ( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
       ( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
   }
 };
 
 struct ModeOrdering {
   bool operator()(const tcDMPtr & d1, const tcDMPtr & d2) {
     ParticleOrdering ord;
     return ord(d1->parent(), d2->parent()) ||
       ( !ord(d2->parent(), d1->parent()) &&
 	( d1->tag() < d2->tag() ||
 	  ( d1->tag() == d2->tag() && d1->fullName() < d2->fullName() ) ) );
   }
 };
 
 void ParticleData::persistentOutput(PersistentOStream & os) const {
   multiset<tcDMPtr,ModeOrdering>
     modes(theDecayModes.begin(), theDecayModes.end());
   os << long(theId) << thePDGName << ounit(theMass, GeV) << ounit(theWidth, GeV)
      << ounit(theHardProcessMass,GeV) << hardProcessMassSet
      << ounit(theHardProcessWidth,GeV) << hardProcessWidthSet
      << ounit(theWidthUpCut, GeV) << ounit(theWidthLoCut, GeV)
      << ounit(theCTau, mm) << oenum(theCharge) << oenum(theSpin)
      << oenum(theColour);
   os << theMassGenerator << isStable << modes << theDecaySelector
      << theWidthGenerator << theVariableRatio << theAntiPartner << syncAnti
      << ounit(theDefMass, GeV) << ounit(theDefWidth, GeV)
      << ounit(theDefCut, GeV) << ounit(theDefCTau, mm) << oenum(theDefColour)
      << oenum(theDefCharge) << oenum(theDefSpin);
 }
 
 void ParticleData::persistentInput(PersistentIStream & is, int) {
   long id;
   is >> id >> thePDGName >> iunit(theMass, GeV) >> iunit(theWidth, GeV)
      >> iunit(theHardProcessMass,GeV) >> hardProcessMassSet
      >> iunit(theHardProcessWidth,GeV) >> hardProcessWidthSet
      >> iunit(theWidthUpCut, GeV) >> iunit(theWidthLoCut, GeV)
      >> iunit(theCTau, mm) >> ienum(theCharge) >> ienum(theSpin)
      >> ienum(theColour) >> theMassGenerator >> isStable
      >> theDecayModes >> theDecaySelector >> theWidthGenerator >> theVariableRatio
      >> theAntiPartner >> syncAnti >> iunit(theDefMass, GeV)
      >> iunit(theDefWidth, GeV) >> iunit(theDefCut, GeV)
      >> iunit(theDefCTau, mm) >> ienum(theDefColour) >> ienum(theDefCharge)
      >> ienum(theDefSpin);
   theId = id;
 }
 
 void ParticleData::Init() {
 
   static ClassDocumentation<ParticleData> documentation
     ("There is no documentation for the ThePEG::ParticleData class");
 
   static Parameter<ParticleData,Energy> interfaceMass
     ("NominalMass",
      "The nominal mass in GeV of the particle. The actual mass "
      "of a particle instance is generated depending on the "
      "nominal mass and the width and is generated by the "
      "<interface>Mass_generator</interface> object associated with the "
      "particle.",
      &ParticleData::theMass, GeV, ZERO, ZERO, Constants::MaxEnergy,
      false, false, Interface::lowerlim,
      &ParticleData::setMass, 0, 0, 0, &ParticleData::defMass);
 
   static Parameter<ParticleData,Energy> interfaceHardProcessMass
     ("HardProcessMass",
      "The mass in GeV of the particle to be used in calculating hard process cross sections.",
      &ParticleData::theHardProcessMass, GeV, ZERO, ZERO, Constants::MaxEnergy,
      false, false, Interface::lowerlim,
      &ParticleData::setHardProcessMass, 0, 0, 0, 0);
 
   static Parameter<ParticleData,Energy> interfaceDefMass
     ("DefaultMass",
      "The default nominal mass in GeV of the particle. The actual mass "
      "of a particle instance is generated depending on the "
      "nominal mass and the width and is generated by the "
      "<interface>Mass_generator</interface> object associated with the "
      "particle.",
      &ParticleData::theDefMass, GeV, ZERO, ZERO, Constants::MaxEnergy,
      false, true, Interface::lowerlim);
   interfaceDefMass.setHasDefault(false);
 
   static Parameter<ParticleData,Energy> interfaceWidth
     ("Width",
      "The width of the particle in GeV.",
      0, GeV, ZERO, ZERO, Constants::MaxEnergy,
      false, false, Interface::lowerlim,
      &ParticleData::setWidth, &ParticleData::getWidth,
      0, 0, &ParticleData::defWidth);
 
   static Parameter<ParticleData,Energy> interfaceHardProcessWidth
     ("HardProcessWidth",
      "The width in GeV of the particle to be used in calculating hard process cross sections.",
      &ParticleData::theHardProcessWidth, GeV, ZERO, ZERO, Constants::MaxEnergy,
      false, false, Interface::lowerlim,
      &ParticleData::setHardProcessWidth, 0, 0, 0, 0);
 
   static Parameter<ParticleData,Energy> interfaceDefWidth
     ("DefaultWidth",
      "The default width of the particle in GeV.",
      &ParticleData::theDefWidth, GeV, ZERO, ZERO, Constants::MaxEnergy,
      false, true, Interface::lowerlim);
   interfaceDefWidth.setHasDefault(false);
 
   static Parameter<ParticleData,Energy> interfaceWidthUpCut
     ("WidthUpCut",
      "The upper hard cutoff in GeV in generated mass, which is the maximum "
      "allowed upwards deviation from the nominal mass. A negative value "
      "corresponds to no limit.",
      0, GeV, ZERO, -1.0*GeV, Constants::MaxEnergy,
      false, false, Interface::lowerlim,
      &ParticleData::setUpCut, &ParticleData::getUpCut,
      0, 0, &ParticleData::defCut);
   
   static Parameter<ParticleData,Energy> interfaceWidthLoCut
     ("WidthLoCut",
      "The lower hard cutoff in GeV in generated mass, which is the maximum "
      "allowed downwards deviation from the nominal mass. A negative value "
      "corresponds to no limit.",
      0, GeV, ZERO, -1.0*GeV, Constants::MaxEnergy,
      false, false, Interface::lowerlim,
      &ParticleData::setLoCut, &ParticleData::getLoCut,
      0, 0, &ParticleData::defCut);
   
   static Parameter<ParticleData,Energy> interfaceWidthCut
     ("WidthCut",
      "The hard cutoff in GeV in generated mass, which is the maximum "
      "allowed deviation from the nominal mass. Sets both the upper and lower "
      "cut. (The displayed value is the maximium of the upper and lower cut.) "
      "A negative value corresponds to no limit.",
      0, GeV, ZERO, -1.0*GeV, Constants::MaxEnergy,
      false, false, Interface::lowerlim,
      &ParticleData::setCut, &ParticleData::getCut,
      0, 0, &ParticleData::defCut);
   interfaceWidthCut.setHasDefault(false);
   
   static Parameter<ParticleData,Energy> interfaceDefWidthCut
     ("DefaultWidthCut",
      "The default hard cutoff in GeV in generated mass, which is the maximum "
      "allowed deviation from the nominal mass. For the actual cutoff, the "
      "upper and lower cut can be set separately.",
      &ParticleData::theDefCut, GeV, ZERO, ZERO, Constants::MaxEnergy,
      false, true, Interface::lowerlim);
   interfaceDefWidthCut.setHasDefault(false);
   
   static Parameter<ParticleData,Length> interfaceCTau
     ("LifeTime",
      "c times the average lifetime of the particle measuerd in mm."
      "The actual lifetime of a particle instance is generated "
      "from this number by the <interface>Mass_generator</interface> "
      "object associated with the particle.",
      0, mm, ZERO, ZERO, Constants::MaxLength,
      false, false, Interface::lowerlim,
      &ParticleData::setCTau, &ParticleData::getCTau,
      0, 0, &ParticleData::defCTau);
   interfaceCTau.setHasDefault(false);
 
   static Parameter<ParticleData,Length> interfaceDefCTau
     ("DefaultLifeTime",
      "c times the default average lifetime of the particle measuerd in mm."
      "The actual lifetime of a particle instance is generated "
      "from this number by the <interface>Mass_generator</interface> "
      "object associated with the particle.",
      &ParticleData::theDefCTau, mm, ZERO, ZERO, Constants::MaxLength,
      false, true, Interface::lowerlim);
   interfaceDefCTau.setHasDefault(false);
 
   static Switch<ParticleData> interfaceColour
     ("Colour",
      "The colour quantum number of this particle type.",
      0, -1, false, false, &ParticleData::setColour, &ParticleData::getColour,
      &ParticleData::defColour);
   static SwitchOption interfaceColourUndefined
     (interfaceColour, "Undefined", "The coulur is undefined.", -1);
   static SwitchOption interfaceColourNeutral
     (interfaceColour, "Neutral", "This particle is colour neutral.", 0);
   static SwitchOption interfaceColour3
     (interfaceColour, "Triplet", "This particle is a colour triplet.", 3);
   static SwitchOption interfaceColour3bar
     (interfaceColour, "AntiTriplet",
      "This particle is a colour anti-triplet.", -3);
   static SwitchOption interfaceColour6
     (interfaceColour, "Sextet", "This particle is a colour sextet.", 6);
   static SwitchOption interfaceColour6bar
     (interfaceColour, "AntiSextet",
      "This particle is a colour anti-sextet.", -6);
   static SwitchOption interfaceColour8
     (interfaceColour, "Octet", "This particle is a colour octet.", 8);
   
   static Switch<ParticleData,PDT::Colour> interfaceDefColour
     ("DefaultColour",
      "The default colour quantum number of this particle type.",
      &ParticleData::theDefColour, PDT::Colour(-1), false, true);
   static SwitchOption interfaceDefColourUndefined
     (interfaceDefColour, "Undefined", "The coulur is undefined.", -1);
   static SwitchOption interfaceDefColourNeutral
     (interfaceDefColour, "Neutral", "This particle is colour neutral.", 0);
   static SwitchOption interfaceDefColour3
     (interfaceDefColour, "Triplet", "This particle is a colour triplet.", 3);
   static SwitchOption interfaceDefColour3bar
     (interfaceDefColour, "AntiTriplet",
      "This particle is a colour anti-triplet.", -3);
   static SwitchOption interfaceDefColour6
     (interfaceDefColour, "Sextet", "This particle is a colour sextet.", 6);
   static SwitchOption interfaceDefColour6bar
     (interfaceDefColour, "AntiSextet",
      "This particle is a colour anti-sextet.", -6);
   static SwitchOption interfaceDefColour8
     (interfaceDefColour, "Octet", "This particle is a colour octet.", 8);
   interfaceDefColour.setHasDefault(false);  
 
   static Parameter<ParticleData, int> interfaceCharge
     ("Charge",
      "The charge of this particle in units of e/3. "
      "See also the command interface <interface>SetCharge</interface>.",
      0, 0, -24, 24, false, false, true,
      &ParticleData::setCharge, &ParticleData::getCharge, 0, 0,
      &ParticleData::defCharge);
 
   static Parameter<ParticleData, PDT::Charge> interfaceDefCharge
     ("DefaultCharge",
      "The default charge of this particle in units of e/3. "
      "See also the command interface <interface>SetCharge</interface>.",
      &ParticleData::theDefCharge, PDT::Charge(0), PDT::Charge(-24),
      PDT::Charge(24), false, true, true);
   interfaceDefCharge.setHasDefault(false);
 
   static Command<ParticleData> interfaceSetCharge
     ("SetCharge",
      "Set the charge of this particle. The argument should be given as an "
      "interger giving three times the unit charge, or 'unknown', "
      "'charged', 'positive' or 'negative'", &ParticleData::ssetCharge);
 
   static Parameter<ParticleData, int> interfaceSpin
     ("Spin",
      "The spin quantim number of this particle on the form 2j+1.",
      0, 0, 0, 9, false, false, true,
      &ParticleData::setSpin, &ParticleData::getSpin, 0, 0,
      &ParticleData::defSpin);
 
   static Parameter<ParticleData, PDT::Spin> interfaceDefSpin
     ("DefaultSpin",
      "The default spin quantim number of this particle on the form 2j+1.",
      &ParticleData::theDefSpin, PDT::Spin(0), PDT::Spin(0), PDT::Spin(9),
      false, true, true);
   interfaceDefSpin.setHasDefault(false);
 
   static Switch<ParticleData> interfaceStable
     ("Stable",
      "Indicates if the particle is stable or not.",
      0, 0, false, false,
      &ParticleData::setStable, &ParticleData::getStable, 0);
   static SwitchOption interfaceStableYes
     (interfaceStable,
      "Stable",
      "This particle is stable",
      1);
   static SwitchOption interfaceStableNo
     (interfaceStable,
      "Unstable",
      "This particle is not stable",
      0);
   interfaceStable.setHasDefault(false);
 
   static Switch<ParticleData> interfaceVariableRatio
     ("VariableRatio",
      "Indicates if the branching ratios of the particle are allowed"
      " to vary for given Particle instances depending on the mass of the instance.",
      0, 0, false, false,
      &ParticleData::setVariableRatio, &ParticleData::getVariableRatio, 0);
   static SwitchOption interfaceVariableRatioYes
     (interfaceVariableRatio,
      "Yes",
      "The branching ratio varies.",
      1);
   static SwitchOption interfaceVariableRatioNo
     (interfaceVariableRatio,
      "No",
      "The branching ratio does not vary.",
      0);
 
   static Switch<ParticleData> interfaceSync
     ("Synchronized",
      "Indicates if the changes to this particle is propagated to "
      "its anti-partner or not. Note that setting this switch does not "
      "actually synchronize the properties with the anti-partner, "
      "it only assures that following changes are propagated. "
      "To sync the particle with its anti-particle, use the "
      "<interface>Synchronize</interface> command.",
      0, 1, false, false,
      &ParticleData::setSync, &ParticleData::getSync, 0);
   static SwitchOption interfaceSyncYes
     (interfaceSync,
      "Synchronized",
      "Changes to this particle will propagate to its "
      "anti-partner",
      1);
   static SwitchOption interfaceSyncNo
     (interfaceSync,
      "Not_synchronized",
      "Changes to this particle will propagate to its "
      "anti-partner",
      0);
   interfaceSync.setHasDefault(false);
 
   static Command<ParticleData> interfaceSynchronize
     ("Synchronize",
      "Synchronizes this particle so that all its properties "
      "correspond to those of its anti-partner",
      &ParticleData::doSync, false);
 
   static Reference<ParticleData,MassGenerator> interfaceMassGenerator
     ("Mass_generator",
      "An object derived from the ThePEG::MassGenerator"
      "class, which is able to generate a mass for a given "
      "particle instance",
      &ParticleData::theMassGenerator, false, false, true, true,
      &ParticleData::setMassGenerator, 0, 0);
 
   static Reference<ParticleData,WidthGenerator> interfaceWidthGenerator
     ("Width_generator",
      "An object derived from the ThePEG::WidthGenerator class, "
      "which is able to calculate the full and partial widths for"
      "this particle type and for a given instance of this "
      "particle type.",
      &ParticleData::theWidthGenerator, false, false, true, true,
      &ParticleData::setWidthGenerator, 0, 0);
 
   static RefVector<ParticleData,DecayMode> interfaceDecayModes
     ("DecayModes",
      "The list of decay modes defined for this particle type.",
      0, -1, false, false, false, false, 0,
      &ParticleData::insDecayModes, &ParticleData::delDecayModes,
      &ParticleData::getDecayModes);
 
 
   static Command<ParticleData> interfaceSelectDecayModes
     ("SelectDecayModes",
      "Only the decay modes which are given as (white-space separated) "
      "decay tags will be switched on, all others will be switched off. "
      "If no argument or 'none' is given, all decay modes are switched off. "
      "If the argument is 'all', all decay modes are switched on.",
      &ParticleData::doSelectDecayModes, false);
 
 
   static Command<ParticleData> interfacePrintDecayModes
     ("PrintDecayModes",
      "Print all decay modes of this particle.",
      &ParticleData::doPrintDecayModes, true);
 
 
 
   static Command<ParticleData> interfaceUnsetHardProcessMass
     ("UnsetHardProcessMass",
      "Unset a previously set hard process mass.",
      &ParticleData::doUnsetHardProcessMass, false);
 
   static Command<ParticleData> interfaceAdjustNominalMass
     ("AdjustNominalMass",
      "Unset a previously set hard process mass.",
      &ParticleData::doAdjustNominalMass, false);
 
   static Command<ParticleData> interfaceUnsetHardProcessWidth
     ("UnsetHardProcessWidth",
      "Unset a previously set hard process width.",
      &ParticleData::doUnsetHardProcessWidth, false);
 
   interfaceStable.rank(14);
   interfaceDecayModes.rank(13);
   interfaceMass.rank(12);
   interfaceWidth.rank(11);
   interfaceWidthCut.rank(10);
   interfaceCTau.rank(9);
   interfaceMassGenerator.rank(8);
   interfaceWidthGenerator.rank(7);
   interfaceWidthUpCut.rank(-0.1);
   interfaceWidthLoCut.rank(-0.1);
 
 
 }
 
 string ParticleData::doPrintDecayModes(string) {
   multimap<double,tDMPtr, std::greater<double> > sorted;
   for ( DecaySet::iterator it = decayModes().begin();
 	it != decayModes().end(); ++it )
     sorted.insert(make_pair((**it).brat(), *it));
   ostringstream os;
   for ( multimap<double,tDMPtr, 
 	  std::greater<double> >::iterator it = sorted.begin();
 	it != sorted.end(); ++it )
     os << it->second->tag()
        << (it->second->on()? " ": " (off) ")
        << it->first << endl;
   return os.str();
 }
 
 
 string ParticleData::doSelectDecayModes(string args) {
   DecaySet on;
   while ( !args.empty() ) {
     string arg = StringUtils::car(args);
     if ( arg == "all" ) {
       on = decayModes();
       break;
     }
     if ( arg == "none" ) {
       on.clear();
       break;
     }
     string name = arg;
     args = StringUtils::cdr(args);
     if ( arg.empty() ) continue;
     if ( arg[0] != '/' ) arg = fullName() + "/" + arg;
     DMPtr dm = Repository::GetPtr<DMPtr>(arg);
     if ( !dm ) return "Error: No decay mode with tag '" + name + "' exists.";
     on.insert(dm);
   }
   for ( DecaySet::iterator it = decayModes().begin();
 	it != decayModes().end(); ++it ) {
     if ( on.find(*it) != on.end() ) {
       (**it).switchOn();
       on.erase(*it);
     } else {
       (**it).switchOff();
     }
   }
   if ( !on.empty() )
     return "Error: decay mode '" + (**on.begin()).tag() + "'was not available.";
   return "";
 }
 
 void ParticleData::insDecayModes(DMPtr dm, int) {
   addDecayMode(dm);
 }
 
 void ParticleData::delDecayModes(int i) {
   vector<DMPtr> mv = getDecayModes();
   if ( i >= 0 && static_cast<unsigned int>(i) < mv.size() ) removeDecayMode(mv[i]);
 }
 
 vector<DMPtr> ParticleData::getDecayModes() const {
   return vector<DMPtr>(theDecayModes.begin(), theDecayModes.end());
 }
 
 ParticleChargeCommand::
 ParticleChargeCommand(const ParticleData & pd, string arg) {
   theMessage << "Cannot set the charge of particle '" << pd.name()
 	     << "' to '" << arg << "'.";
   severity(warning);
 }
 
 void ParticleData::doinit() {
   Interfaced::doinit();
   if( theMassGenerator )  theMassGenerator->init();
   if( theWidthGenerator ) theWidthGenerator->init();
 }
 
 void ParticleData::doinitrun() {
   Interfaced::doinitrun();
   if( theMassGenerator )  theMassGenerator->initrun();
   if( theWidthGenerator ) theWidthGenerator->initrun();
 }
 
 }
diff --git a/PDT/QuarksToHadronsDecayer.cc b/PDT/QuarksToHadronsDecayer.cc
--- a/PDT/QuarksToHadronsDecayer.cc
+++ b/PDT/QuarksToHadronsDecayer.cc
@@ -1,235 +1,235 @@
 // -*- C++ -*-
 //
 // QuarksToHadronsDecayer.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 QuarksToHadronsDecayer class.
 //
 
 #include "QuarksToHadronsDecayer.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/PDT/DecayMode.h"
 #include "ThePEG/PDT/StandardMatchers.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Utilities/SimplePhaseSpace.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace ThePEG;
 
 QuarksToHadronsDecayer::~QuarksToHadronsDecayer() {}
 
 IBPtr QuarksToHadronsDecayer::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr QuarksToHadronsDecayer::fullclone() const {
   return new_ptr(*this);
 }
 
 bool QuarksToHadronsDecayer::accept(const DecayMode & dm) const {
   int col = 0;
   int acol = 0;
   if ( !dm.productMatchers().empty() ) {
     for ( MatcherMSet::const_iterator it = dm.productMatchers().begin();
 	  it != dm.productMatchers().end(); ++it ) {
       if ( typeid(**it) == typeid(MatchLightQuark) ) ++col;
       else if ( typeid(**it) == typeid(MatchLightAntiQuark) ) ++acol;
       else return false;
     }
     if ( col != 1 || col != acol ) return false;
   }
   if ( dm.orderedProducts().size() + col + acol < 2 ||
        !dm.cascadeProducts().empty() || dm.wildProductMatcher() ) return false;
   for ( int i = 0, N = dm.orderedProducts().size(); i < N; ++i ) {
     if ( DiquarkMatcher::Check(*dm.orderedProducts()[i]) ) {
       if ( i + 1 != N ) return false;
       if ( dm.orderedProducts()[i]->id() < 0 ) ++col;
       else ++acol;
     }
     if ( QuarkMatcher::Check(*dm.orderedProducts()[i]) ) {
       if ( dm.orderedProducts()[i]->id() > 0 ) ++col;
       else ++acol;
     }
   }
   if ( acol != col || col < 1 || col > 2 ) return false;
   return true;
 }
 
 PVector QuarksToHadronsDecayer::decay(const DecayMode & dm,
 				      const Particle & parent) const {
   PVector children;
   tcPDVector quarks;
   if ( !dm.productMatchers().empty() ) {
     tcPDPtr pd = getParticleData(flavourGenerator()->selectQuark());
     quarks.push_back(pd);
     quarks.push_back(pd->CC());
   }
   Energy summq = ZERO;
   Energy summp = ZERO;
   tPDVector prods = dm.orderedProducts();
   for ( int i = 0, N = prods.size(); i < N; ++i )
     if ( QuarkMatcher::Check(*prods[i]) || DiquarkMatcher::Check(*prods[i])) {
       quarks.push_back(prods[i]);
       summq += quarks.back()->mass();
     } else {
       children.push_back(prods[i]->produceParticle());
       summp += children.back()->mass();
     }
 
   Energy summh = ZERO;
   PVector hadrons;
 
   if ( !quarks.empty() ) do {
 
     hadrons = getHadrons(getN(parent.mass(), summq, quarks.size()), quarks);
 
     summh = ZERO;
     for ( int i = 0, N = hadrons.size(); i < N; ++i )
       summh += hadrons[i]->mass();
 
   } while ( hadrons.empty() || summp + summh >= parent.mass() );
 
   children.insert(children.end(), hadrons.begin(), hadrons.end());
 
   distribute(parent, children);
 
   finalBoost(parent, children);
   setScales(parent, children);
 
   return children;
 
 }
 
 int QuarksToHadronsDecayer::getN(Energy m0, Energy summq, int Nq) const {
   int Nh = fixedN();
   if ( Nh >= 2 ) return Nh;
 
   double c = c1()*log((m0 - summq)/c2()) + c3();
   if ( c < 0.0 ) return minN();
   while ( true ) {
     using namespace Constants;
     Nh = int(0.5 + double(Nq)/4.0 + c +
 	     sqrt(-2.0*c*log(max(1.0e-10, rnd())))*sin(2.0*pi*rnd()));
     if ( Nh >= minN() ) return Nh;
   }
 }
 
 PVector QuarksToHadronsDecayer::
 getHadrons(int Nh, tcPDVector quarks) const {
   PVector hadrons;
   Nh -= quarks.size()/2;
   while ( Nh-- > 0 ) {
     int i = irnd(quarks.size() - 1);
     tcPDPair hq = flavourGenerator()->alwaysGenerateHadron(quarks[i]);
     hadrons.push_back(hq.first->produceParticle());
     quarks[i] = hq.second;
   }
   if ( DiquarkMatcher::Check(*quarks[0]) && DiquarkMatcher::Check(*quarks[1]) )
     return PVector();
   tcPDPtr h = flavourGenerator()->alwaysGetHadron(quarks[0], quarks[1]);
   hadrons.push_back(h->produceParticle());
   if ( quarks.size() <= 2 ) return hadrons;
   if ( DiquarkMatcher::Check(*quarks[2]) && DiquarkMatcher::Check(*quarks[3]) )
     return PVector();
   h = flavourGenerator()->alwaysGetHadron(quarks[2], quarks[3]);
   hadrons.push_back(h->produceParticle());
   return hadrons;
 }
 
 void QuarksToHadronsDecayer::
 distribute(const Particle & parent, PVector & children) const {
   do {
     try {
       SimplePhaseSpace::CMSn(children, parent.mass());
     }
-    catch ( ImpossibleKinematics ) {
+    catch ( ImpossibleKinematics & e) {
       children.clear();
       return;
     }
   } while ( reweight(parent, children) < rnd() );
 }
 
 double QuarksToHadronsDecayer::
 reweight(const Particle &, const PVector &) const {
   return 1.0;
 }
 
 void QuarksToHadronsDecayer::persistentOutput(PersistentOStream & os) const {
   os << theFixedN << theMinN << theC1 << ounit(theC2,GeV) << theC3 << theFlavourGenerator;
 }
 
 void QuarksToHadronsDecayer::persistentInput(PersistentIStream & is, int) {
   is >> theFixedN >> theMinN >> theC1 >> iunit(theC2,GeV) >> theC3 >> theFlavourGenerator;
 }
 
 ClassDescription<QuarksToHadronsDecayer> QuarksToHadronsDecayer::initQuarksToHadronsDecayer;
 // Definition of the static class description member.
 
 void QuarksToHadronsDecayer::Init() {
 
   static ClassDocumentation<QuarksToHadronsDecayer> documentation
     ("This class decays particles to nq (2 or 4) quarks which then are "
      "decayes to hadrons according to phase space. The number of final "
      "hadrons can either be given by a fixed number or as a Gaussian "
      "multiplicity distribution centered around c+nq/4+c3 and a width "
      "sqrt(c), where c = c1 log((m - summ)/c2), m is the mass of the "
      "decaying particle, summ the sum of the quark masses and ci real "
      "parameters.");
 
   static Parameter<QuarksToHadronsDecayer,int> interfaceFixedN
     ("FixedN",
      "The fixed number of hadrons to be produced. If less than 2, the "
      "number is instead given by a gaussian multiplicity distribution.",
      &QuarksToHadronsDecayer::theFixedN, 0, 0, 10,
      true, false, true);
 
   static Parameter<QuarksToHadronsDecayer,int> interfaceMinN
     ("MinN",
      "The minimum hadrons to be produced.",
      &QuarksToHadronsDecayer::theMinN, 2, 2, 10,
      true, false, true);
 
   static Parameter<QuarksToHadronsDecayer,double> interfaceC1
     ("C1",
      "The c1 parameter of the gaussian multiplicity distribution centered "
      "around c1 log((m - summ)/c2) +c3.",
      &QuarksToHadronsDecayer::theC1, 4.5, 0.0, 10.0,
      true, false, true);
 
   static Parameter<QuarksToHadronsDecayer,Energy> interfaceC2
     ("C2",
      "The c2 parameter of the gaussian multiplicity distribution centered "
      "around c1 log((m - summ)/c2) +c3.",
      &QuarksToHadronsDecayer::theC2, GeV, 0.7*GeV, ZERO, 10.0*GeV,
      true, false, true);
 
   static Parameter<QuarksToHadronsDecayer,double> interfaceC3
     ("C3",
      "The c3 parameter of the gaussian multiplicity distribution centered "
      "around c1 log((m - summ)/c2) +c3.",
      &QuarksToHadronsDecayer::theC3, 0.0, 0.0, 10.0,
      true, false, true);
 
   static Reference<QuarksToHadronsDecayer,FlavourGenerator>
     interfaceFlavourGenerator
     ("FlavourGenerator",
      "The object in charge of generating hadrons spieces from given quark "
      "flavours.",
      &QuarksToHadronsDecayer::theFlavourGenerator,
      true, false, true, false, true);
 
   interfaceFixedN.rank(10);
   interfaceMinN.rank(9);
   interfaceFlavourGenerator.rank(8);
   interfaceMinN.setHasDefault(false);;
 
 }
 
diff --git a/Repository/EventGenerator.cc b/Repository/EventGenerator.cc
--- a/Repository/EventGenerator.cc
+++ b/Repository/EventGenerator.cc
@@ -1,1400 +1,1401 @@
 // -*- C++ -*-
 //
 // EventGenerator.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 EventGenerator class.
 //
 
 #include "EventGenerator.h"
 #include "EventGenerator.xh"
 #include "ThePEG/Handlers/EventHandler.h"
 #include "Repository.h"
 #include "ThePEG/Utilities/HoldFlag.h"
 #include "ThePEG/Utilities/Debug.h"
 #include "ThePEG/Utilities/DebugItem.h"
 #include "ThePEG/Interface/Interfaced.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/RefVector.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/Command.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/PDT/ParticleData.h"
 #include "ThePEG/PDT/MatcherBase.h"
 #include "ThePEG/PDT/DecayMode.h"
 #include "ThePEG/StandardModel/StandardModelBase.h"
 #include "ThePEG/Repository/Strategy.h"
 #include "ThePEG/Repository/CurrentGenerator.h"
 #include "ThePEG/Handlers/AnalysisHandler.h"
 #include "ThePEG/Analysis/FactoryBase.h"
 #include "ThePEG/Handlers/EventManipulator.h"
 #include "ThePEG/Handlers/LuminosityFunction.h"
 #include "ThePEG/MatrixElement/MEBase.h"
 #include "ThePEG/EventRecord/Event.h"
 #include "ThePEG/Handlers/SubProcessHandler.h"
 #include "ThePEG/Handlers/CascadeHandler.h"
 #include "ThePEG/Handlers/HadronizationHandler.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Config/algorithm.h"
 #include "ThePEG/Utilities/DynamicLoader.h"
 #include <cstdlib>
 #include "ThePEG/Repository/Main.h"
 #include <csignal>
 
 #ifdef ThePEG_TEMPLATES_IN_CC_FILE
 #include "EventGenerator.tcc"
 #endif
 
 using namespace ThePEG;
 
 namespace {
   volatile sig_atomic_t THEPEG_SIGNAL_STATE = 0;
 }
 
 // signal handler function
 // very restricted in what it is allowed do
 // without causing undefined behaviour
 extern "C" {
   void thepegSignalHandler(int id) {
     THEPEG_SIGNAL_STATE=id;
     signal(id,SIG_DFL);
   }
 }
 
 void EventGenerator::checkSignalState() {
   if (THEPEG_SIGNAL_STATE) {
     log() << "Caught signal " << THEPEG_SIGNAL_STATE << ". Exiting ..." << std::endl;
     finalize();
     exit(0);
   }
 }
 
 EventGenerator::EventGenerator()
   : thePath("."), theNumberOfEvents(1000), theQuickSize(7000),
     preinitializing(false), ieve(0), weightSum(0.0),
     theDebugLevel(0), logNonDefault(-1), printEvent(0), dumpPeriod(0),
     keepAllDumps(false),
     debugEvent(0), maxWarnings(10), maxErrors(10), theCurrentRandom(0),
     theCurrentGenerator(0), useStdout(false), theIntermediateOutput(false) {}
 
 EventGenerator::EventGenerator(const EventGenerator & eg)
   : Interfaced(eg), theDefaultObjects(eg.theDefaultObjects),
     theLocalParticles(eg.theLocalParticles),
     theStandardModel(eg.theStandardModel),
     theStrategy(eg.theStrategy), theRandom(eg.theRandom),
     theEventHandler(eg.theEventHandler),
     theAnalysisHandlers(eg.theAnalysisHandlers),
     theHistogramFactory(eg.theHistogramFactory),
     theEventManipulator(eg.theEventManipulator),
     thePath(eg.thePath), theRunName(eg.theRunName),
     theNumberOfEvents(eg.theNumberOfEvents), theObjects(eg.theObjects),
     theObjectMap(eg.theObjectMap),
     theParticles(eg.theParticles), theQuickParticles(eg.theQuickParticles),
     theQuickSize(eg.theQuickSize), preinitializing(false),
     theMatchers(eg.theMatchers),
     usedObjects(eg.usedObjects), ieve(eg.ieve), weightSum(eg.weightSum),
     theDebugLevel(eg.theDebugLevel), logNonDefault(eg.logNonDefault),
     printEvent(eg.printEvent), dumpPeriod(eg.dumpPeriod),
     keepAllDumps(eg.keepAllDumps),
     debugEvent(eg.debugEvent),
     maxWarnings(eg.maxWarnings), maxErrors(eg.maxErrors), theCurrentRandom(0),
     theCurrentGenerator(0),
     theCurrentEventHandler(eg.theCurrentEventHandler),
     theCurrentStepHandler(eg.theCurrentStepHandler),
     useStdout(eg.useStdout),
     theIntermediateOutput(eg.theIntermediateOutput) {}
 
 EventGenerator::~EventGenerator() {
   if ( theCurrentRandom ) delete theCurrentRandom;
   if ( theCurrentGenerator ) delete theCurrentGenerator;
 }
 
 IBPtr EventGenerator::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr EventGenerator::fullclone() const {
   return new_ptr(*this);
 }
 
 tcEventPtr EventGenerator::currentEvent() const {
   return eventHandler()->currentEvent();
 }
 
 CrossSection EventGenerator::histogramScale() const {
   return eventHandler()->histogramScale();
 }
 
 CrossSection EventGenerator::integratedXSec() const {
   return eventHandler()->integratedXSec();
 }
 
 CrossSection EventGenerator::integratedXSecErr() const {
   return eventHandler()->integratedXSecErr();
 }
 
 void EventGenerator::setSeed(long seed) {
   random().setSeed(seed);
   ostringstream s;
   s << seed;
   const InterfaceBase * ifb = BaseRepository::FindInterface(theRandom, "Seed");
   ifb->exec(*theRandom, "set", s.str());
 }
 
 void
 EventGenerator::setup(string newRunName,
 		      ObjectSet & newObjects,
 		      ParticleMap & newParticles,
 		      MatcherSet & newMatchers) {
   HoldFlag<int> debug(Debug::level, Debug::isset? Debug::level: theDebugLevel);
   theRunName = newRunName;
   theObjects.swap(newObjects);
   theParticles.swap(newParticles);
   theMatchers.swap(newMatchers);
   theObjectMap.clear();
   for ( ObjectSet::const_iterator it = objects().begin();
 	it != objects().end(); ++it ) theObjectMap[(**it).fullName()] = *it;
   UseRandom currentRandom(theRandom);
   CurrentGenerator currentGenerator(this);
 
   // Force update of all objects and then reset.
   touch();
   for_each(theObjects, mem_fun(&InterfacedBase::touch));
   update();
   for_each(theObjects, mem_fun(&InterfacedBase::update));
   clear();
   BaseRepository::clearAll(theObjects);
 
   init();
 
 }
 
 IBPtr EventGenerator::getPointer(string name) const {
   ObjectMap::const_iterator it = objectMap().find(name);
   if ( it == objectMap().end() ) return IBPtr();
   else return it->second;
 }
 
 void EventGenerator::openOutputFiles() {
   if ( !useStdout ) {
     logfile().open((filename() + ".log").c_str());
     theOutFileName = filename() + ".out";
     outfile().open(theOutFileName.c_str());
     outfile().close();
     theOutStream.str("");
   }
   out() << Repository::banner() << endl;
   log() << Repository::banner() << endl;
 }
 
 void EventGenerator::closeOutputFiles() {
   flushOutputFile();
   if ( !useStdout ) logfile().close();
 }
 
 void EventGenerator::flushOutputFile() {
   if ( !useStdout ) {
     outfile().open(theOutFileName.c_str(), ios::out|ios::app);
     outfile() << theOutStream.str();
     outfile().close();
   } else
     BaseRepository::cout() << theOutStream.str();
   theOutStream.str("");
 }
 
 void EventGenerator::doinit() {
 
   HoldFlag<int> debug(Debug::level, Debug::isset? Debug::level: theDebugLevel);
 
   // First initialize base class and random number generator.
   Interfaced::doinit();
   random().init();
 
   // Make random generator and this available in standard static
   // classes.
   UseRandom useRandom(theRandom);
   CurrentGenerator currentGenerator(this);
 
   // First initialize all objects which have requested this by
   // implementing a InterfacedBase::preInitialize() function which
   // returns true.
   while ( true ) {
     HoldFlag<bool> hold(preinitializing, true);
     ObjectSet preinits;
     for ( ObjectSet::iterator it = objects().begin();
 	  it != objects().end(); ++it )
       if ( (**it).preInitialize() &&
 	   (**it).state() == InterfacedBase::uninitialized )
 	preinits.insert(*it);
     if ( preinits.empty() ) break;
     for_each(preinits, mem_fun(&InterfacedBase::init));
   }
 
   // Initialize the quick access to particles.
   theQuickParticles.clear();
   theQuickParticles.resize(2*theQuickSize);
   for ( ParticleMap::const_iterator pit = theParticles.begin();
 	pit != theParticles.end(); ++pit )
     if ( abs(pit->second->id()) < theQuickSize )
       theQuickParticles[pit->second->id()+theQuickSize] = pit->second;
 
   // Then call the init method for all objects. Start with the
   // standard model and the strategy.
   standardModel()->init();
   if ( strategy() ) strategy()->init();
   eventHandler()->init();
 
   // initialize particles first
   for(ParticleMap::const_iterator pit = particles().begin();
       pit != particles().end(); ++pit) pit->second->init();
   
   for_each(objects(), mem_fun(&InterfacedBase::init));
 
   // Then initialize the Event Handler calculating initial cross
   // sections and stuff.
   eventHandler()->initialize();
 }
 
 void EventGenerator::doinitrun() {
 
   HoldFlag<int> debug(Debug::level, Debug::isset? Debug::level: theDebugLevel);
 
   signal(SIGHUP, thepegSignalHandler);
   signal(SIGINT, thepegSignalHandler);
   signal(SIGTERM,thepegSignalHandler);
 
   currentEventHandler(eventHandler());
 
   Interfaced::doinitrun();
   random().initrun();
 
   // Then call the init method for all objects. Start with the
   // standard model and the strategy.
   standardModel()->initrun();
   if ( strategy() ) { 
     strategy()->initrun();
     if ( ! strategy()->versionstring().empty() ) {
       out() << ">> " << strategy()->versionstring() << '\n' << endl;
       log() << ">> " << strategy()->versionstring() << '\n' << endl;
     }
   }
   // initialize particles first
   for(ParticleMap::const_iterator pit = particles().begin();
       pit != particles().end(); ++pit) {
     pit->second->initrun();
   }
   eventHandler()->initrun();
 
   
   for_each(objects(), mem_fun(&InterfacedBase::initrun));
 
   if ( logNonDefault > 0 || ( ThePEG_DEBUG_LEVEL && logNonDefault == 0 ) ) {
     vector< pair<IBPtr, const InterfaceBase *> > changed =
       Repository::getNonDefaultInterfaces(objects());
     if ( changed.size() ) {
       log() << string(78, '=') << endl
 	    << "The following interfaces have non-default values (default):"
 	    << endl << string(78, '-') << endl;
       for ( int i = 0, N = changed.size(); i < N; ++i ) {
 	log() << changed[i].first->fullName() << ":"
 	      << changed[i].second->name() << " = "
 	      << changed[i].second->exec(*changed[i].first, "notdef", "")
 	      << endl;
       }
       log() << string(78,'=') << endl;
     }
   }
 
   weightSum = 0.0;
 
 }
 
 PDPtr EventGenerator::getParticleData(PID id) const {
   long newId = id;
   if ( abs(newId) < theQuickSize && theQuickParticles.size() )
     return theQuickParticles[newId+theQuickSize];
   ParticleMap::const_iterator it = theParticles.find(newId);
   if ( it == theParticles.end() ) return PDPtr();
   return it->second;
 }
 
 PPtr EventGenerator::getParticle(PID newId) const {
   tcPDPtr pd = getParticleData(newId);
   if ( !pd ) return PPtr();
   return pd->produceParticle();
 }
 
 void EventGenerator::finalize() {
   UseRandom currentRandom(theRandom);
   CurrentGenerator currentGenerator(this);
   finish();
   finally();
 }
 
 void EventGenerator::dofinish() {
 
   HoldFlag<int> debug(Debug::level, Debug::isset? Debug::level: theDebugLevel);
 
   // first write out statistics from the event handler.
   eventHandler()->statistics(out());
 
   // Call the finish method for all other objects.
   for_each(objects(), mem_fun(&InterfacedBase::finish));
 
   if ( theExceptions.empty() ) {
     log() << "No exceptions reported in this run.\n";
   } else {
 
     log() << "\nThe following exception classes were reported in this run:\n";
 
     for ( ExceptionMap::iterator it = theExceptions.begin();
 	  it != theExceptions.end(); ++it ) {
       string severity;
       switch ( it->first.second ) {
       case Exception::info       : severity="info"; break;
       case Exception::warning    : severity="warning"; break;
       case Exception::setuperror : severity="setuperror"; break;
       case Exception::eventerror : severity="eventerror"; break;
       case Exception::runerror   : severity="runerror"; break;
       case Exception::maybeabort : severity="maybeabort"; break;
       case Exception::abortnow   : severity="abortnow"; break;
       default                    : severity="unknown";
       }
       log() << it->first.first << ' ' << severity
 	    << " (" << it->second << " times)\n";
     }
   }
 
   theExceptions.clear();
 
   const string & msg = theMiscStream.str();
   if ( ! msg.empty() ) {
     log() << endl 
 	  << "Miscellaneous output from modules to the standard output:\n\n"
 	  << msg;
     theMiscStream.str("");
   }
 
   flushOutputFile();
 
 }
 
 void EventGenerator::finally() {
 
   generateReferences();
 
   closeOutputFiles();
 
   if ( theCurrentRandom ) delete theCurrentRandom;
   if ( theCurrentGenerator ) delete theCurrentGenerator;
   theCurrentRandom = 0;
   theCurrentGenerator = 0;
 
 }
 
 void EventGenerator::initialize(bool initOnly) {
   UseRandom currentRandom(theRandom);
   CurrentGenerator currentGenerator(this);
   doInitialize(initOnly);
 }
 
 bool EventGenerator::loadMain(string file) {
   initialize();
   UseRandom currentRandom(theRandom);
   CurrentGenerator currentGenerator(this);
   Main::eventGenerator(this);
   bool ok = DynamicLoader::load(file);
   finish();
   finally();
   return ok;
 }
 
 
 void EventGenerator::go(long next, long maxevent, bool tics) {
   UseRandom currentRandom(theRandom);
   CurrentGenerator currentGenerator(this);
   doGo(next, maxevent, tics);
 }
 
 EventPtr EventGenerator::shoot() {
   static DebugItem debugfpu("ThePEG::FPU", 1);
   if ( debugfpu ) Debug::unmaskFpuErrors();
   UseRandom currentRandom(theRandom);
   CurrentGenerator currentGenerator(this);
   checkSignalState();
   EventPtr event = doShoot();
   if ( event ) weightSum += event->weight();
   DebugItem::tic();
   return event;
 }
 
 EventPtr EventGenerator::doShoot() {
   EventPtr event;
   if ( N() >= 0 && ++ieve > N() ) return event;
   HoldFlag<int> debug(Debug::level, Debug::isset? Debug::level: theDebugLevel);
   do { 
     int state = 0;
     int loop = 1;
     eventHandler()->clearEvent();
     try {
       do {
 	// Generate a full event or part of an event
 	if ( eventHandler()->empty() ) event = eventHandler()->generateEvent();
 	else event = eventHandler()->continueEvent();
 
 	if ( eventHandler()->empty() ) loop = -loop;
 	
 	// Analyze the possibly uncomplete event
 	for ( AnalysisVector::iterator it = analysisHandlers().begin();
 	      it != analysisHandlers().end(); ++it )
 	  (**it).analyze(event, ieve, loop, state);
 	
 	// Manipulate the current event, possibly deleting some steps
 	// and telling the event handler to redo them.
 	if ( manipulator() )
 	  state = manipulator()->manipulate(eventHandler(), event);
 	
 	// If the event was not completed, continue generation and continue.
 	loop = abs(loop) + 1;
       } while ( !eventHandler()->empty() );
     }
     catch (Exception & ex) {
       if ( logException(ex, eventHandler()->currentEvent()) ) throw;
     }
     catch (...) {
       event = eventHandler()->currentEvent();
       if ( event )
 	log() << *event;
       else
 	log() << "An exception occurred before any event object was created!";
       log() << endl;
       dump();
       throw;
     }
     if ( ThePEG_DEBUG_LEVEL ) {
       if ( ( ThePEG_DEBUG_LEVEL == Debug::printEveryEvent ||
 	     ieve < printEvent ) && event ) log() << *event;
       if ( debugEvent > 0 && ieve + 1 >= debugEvent )
 	Debug::level = Debug::full;
     }
   } while ( !event );
 
   // If scheduled, dump a clean state between events
   if ( ThePEG_DEBUG_LEVEL && dumpPeriod > 0 && ieve%dumpPeriod == 0 ) {
     eventHandler()->clearEvent();
     eventHandler()->clean();
     dump();
   }
 
   return event;
 }
 
 EventPtr EventGenerator::doGenerateEvent(tEventPtr e) {
   if ( N() >= 0 && ++ieve > N() ) return EventPtr();
   EventPtr event = e;
   try {
     event = eventHandler()->generateEvent(e);
   }
   catch (Exception & ex) {
     if ( logException(ex, eventHandler()->currentEvent()) ) throw;
   }
   catch (...) {
     event = eventHandler()->currentEvent();
     if ( !event ) event = e;
     log() << *event << endl;
     dump();
     throw;
   }
   return event;
 }
 
 EventPtr EventGenerator::doGenerateEvent(tStepPtr s) {
   if ( N() >= 0 && ++ieve > N() ) return EventPtr();
   EventPtr event;
   try {
     event = eventHandler()->generateEvent(s);
   }
   catch (Exception & ex) {
     if ( logException(ex, eventHandler()->currentEvent()) ) throw;
   }
   catch (...) {
     event = eventHandler()->currentEvent();
     if ( event ) log() << *event << endl;
     dump();
     throw;
   }
   return event;
 }
 
 EventPtr EventGenerator::generateEvent(Event & e) {
   UseRandom currentRandom(theRandom);
   CurrentGenerator currentGenerator(this);
   EventPtr event = doGenerateEvent(tEventPtr(&e));
   if ( event ) weightSum += event->weight();
   return event;
 }
 
 EventPtr EventGenerator::generateEvent(Step & s) {
   UseRandom currentRandom(theRandom);
   CurrentGenerator currentGenerator(this);
   EventPtr event = doGenerateEvent(tStepPtr(&s));
   if ( event ) weightSum += event->weight();
   return event;
 }
 
 Energy EventGenerator::maximumCMEnergy() const {
   tcEHPtr eh = eventHandler();
   return eh->lumiFnPtr()? eh->lumiFn().maximumCMEnergy(): ZERO;
 }
 
 void EventGenerator::doInitialize(bool initOnly) {
 
   if ( !initOnly )
     openOutputFiles();
 
   init();
 
   if ( !initOnly )
     initrun();
 
   if ( !ThePEG_DEBUG_LEVEL ) Exception::noabort = true;
 
 }
 
 void EventGenerator::doGo(long next, long maxevent, bool tics) {
 
   if ( maxevent >= 0 ) N(maxevent);
 
   if ( next >= 0 ) {
     if ( tics ) 
       cerr << "event> " << setw(9) << "init\r" << flush;
     initialize();
     ieve = next-1;
   } else {
     openOutputFiles();
   }
 
   if ( tics ) tic();
   try {
     while ( shoot() ) {
       if ( tics ) tic();
     }
   }
   catch ( ... ) {
     finish();
     throw;
   }
 
   finish();
 
   finally();
 
 }
 
 void EventGenerator::tic(long currev, long totev) const {
   if ( !currev ) currev = ieve;
   if ( !totev ) totev = N();
   long i = currev;
   long n = totev;
   bool skip = currev%(max(totev/100, 1L));
   if ( i > n/2 ) i = n-i;
   while ( skip && i >= 10 && !(i%10) ) i /= 10;
   if ( i == 1 || i == 2 || i == 5 ) skip = false;
   if (!theIntermediateOutput) { //default
     if ( skip ) return;
     cerr << "event> " << setw(8) << currev << " " << setw(8) << totev << "\r";
   }
   else if (theIntermediateOutput) {
     if ( skip && currev%10000!=0) return;
     cerr <<  "event> " << setw(9) << right << currev << "/" << totev 
          << "; xs = " << integratedXSec()/picobarn << " pb +- " 
          << integratedXSecErr()/picobarn << " pb" << endl;
   }
   cerr.flush();
   if ( currev == totev ) cerr << endl;
 }
   
 
 void EventGenerator::dump() const {
   if ( dumpPeriod > -1 ) {
     string dumpfile;
     if ( keepAllDumps ) {
       ostringstream number;
       number << ieve;
       dumpfile = filename() + "-" + number.str() + ".dump";
     }
     else
       dumpfile = filename() + ".dump";
     PersistentOStream file(dumpfile, globalLibraries());
     file << tcEGPtr(this);
   }
 }
 
 void EventGenerator::use(const Interfaced & i) {
   IBPtr ip = getPtr(i);
   if ( ip ) usedObjects.insert(ip);
 }
 
 void EventGenerator::generateReferences() {
   typedef map<string,string> StringMap;
   StringMap references;
 
   // First get all model descriptions and model references from the
   // used objects. Put them in a map indexed by the description to
   // avoid duplicates.
   for ( ObjectSet::iterator it = usedObjects.begin();
 	it != usedObjects.end(); ++it ) {
     if ( *it == strategy() ) continue;
     string desc = Repository::getModelDescription(*it);
     if ( desc.empty() ) continue;
     if ( dynamic_ptr_cast<cEHPtr>(*it) ) desc = "A " + desc;
     else if ( dynamic_ptr_cast<cSMPtr>(*it) ) desc = "B " + desc;
     else if ( dynamic_ptr_cast<cMEPtr>(*it) ) desc = "C " + desc;
     else if ( dynamic_ptr_cast<cCascHdlPtr>(*it) ) desc = "D " + desc;
     else if ( dynamic_ptr_cast<cHadrHdlPtr>(*it) ) desc = "E " + desc;
     else if ( dynamic_ptr_cast<cStepHdlPtr>(*it) ) desc = "F " + desc;
     else if ( dynamic_ptr_cast<cDecayerPtr>(*it) ) desc = "Y " + desc;
     else if ( dynamic_ptr_cast<cAnalysisHdlPtr>(*it) ) desc = "Z " + desc;
     else if ( dynamic_ptr_cast<Ptr<HandlerBase>::const_pointer>(*it) )
       desc = "G " + desc;
     else  desc = "H " + desc;
     references[desc] = Repository::getModelReferences(*it);
   }
 
   // Now get the main strategy description which should put first and
   // remove it from the map.
   string stratdesc;
   string stratref;
   if ( strategy() ) {
     stratdesc = Repository::getModelDescription(strategy());
     stratref = Repository::getModelReferences(strategy());
     references.erase(stratdesc);
   }
 
   // Open the file and write out an appendix header
   if ( !useStdout )
     reffile().open((filename() + ".tex").c_str());
   ref() << "\\documentclass{article}\n"
 	<< "\\usepackage{graphics}\n"
 	<< "\\begin{document}\n"
 	<< "\\appendix\n"
 	<< "\\section[xxx]{\\textsc{ThePEG} version " << Repository::version()
 	<< " \\cite{ThePEG} Run Information}\n"
 	<< "Run name: \\textbf{" << runName()
 	<< "}:\\\\\n";
   if ( !stratdesc.empty() )
     ref() << "This run was generated using " << stratdesc 
 	  << " and the following models:\n";
   else
     ref() << "The following models were used:\n";
 
   ref() << "\\begin{itemize}\n";
 
   // Write out all descriptions.
   for ( StringMap::iterator it = references.begin();
 	it != references.end(); ++it )
     ref() << "\\item " << it->first.substr(2) << endl;
 
   // Write out thebibliography header and all references.
   ref() << "\\end{itemize}\n\n"
 	    << "\\begin{thebibliography}{99}\n"
 	    << "\\bibitem{ThePEG} L.~L\\\"onnblad, "
 	    << "Comput.~Phys.~Commun.\\ {\\bf 118} (1999) 213.\n";
   if ( !stratref.empty() ) ref() << stratref << '\n';
   for ( StringMap::iterator it = references.begin();
 	it != references.end(); ++it )
     ref() << it->second << '\n';
   ref() << "\\end{thebibliography}\n"
 	    << "\\end{document}" << endl;
   if ( !useStdout )
     reffile().close();
 }
 
 void EventGenerator::strategy(StrategyPtr s) {
   theStrategy = s;
 }
 
 int EventGenerator::count(const Exception & ex) {
   return ++theExceptions[make_pair(StringUtils::typeName(typeid(ex)),
 				   ex.severity())];
 }
 
 void EventGenerator::printException(const Exception & ex) {
   switch ( ex.severity() ) {
   case Exception::info:
     log() << "* An information";
         break;
   case Exception::warning:
     log() << "* A warning";
     break;
   case Exception::setuperror:
     log() << "** A setup";
     break;
   case Exception::eventerror:
     log() << "** An event";
     break;
   case Exception::runerror:
     log() << "*** An run";
     break;
   case Exception::maybeabort:
   case Exception::abortnow:
     log() << "**** A serious";
     break;
   default:
     log() << "**** An unknown";
     break;
   }
   if ( ieve > 0 )
     log() << " exception of type " << StringUtils::typeName(typeid(ex))
 	  << " occurred while generating event number "
 	  << ieve << ": \n" << ex.message() << endl;
   else
     log() << " exception occurred in the initialization of "
 	  << name() << ": \n" << ex.message() << endl;
   if ( ex.severity() == Exception::eventerror )
     log() << "The event will be discarded." << endl;
 }
 
 void EventGenerator::logWarning(const Exception & ex) {
   if ( ex.severity() != Exception::info &&
        ex.severity() != Exception::warning ) throw ex;
   ex.handle();  
   int c = count(ex);
   if ( c > maxWarnings ) return;
   printException(ex);
   if ( c == maxWarnings )
     log() << "No more warnings of this kind will be reported." << endl;
 }
  
 bool EventGenerator::
 logException(const Exception & ex, tcEventPtr event) {
   bool noEvent = !event;
   ex.handle();
   int c = count(ex);
   if ( c <= maxWarnings ) {
     printException(ex);
     if ( c == maxWarnings )
       log() << "No more warnings of this kind will be reported." << endl;
   }
   if ( ex.severity() == Exception::info ||
        ex.severity() == Exception::warning ) {
     ex.handle();
     return false;
   }
   if ( ex.severity() == Exception::eventerror ) {
     if ( c < maxErrors || maxErrors <= 0 ) {
       ex.handle();
       if ( ThePEG_DEBUG_LEVEL > 0 && !noEvent ) log() << *event;
       return false;
     }
     if ( c > maxErrors ) printException(ex);
     log() << "Too many (" << c << ") exceptions of this kind has occurred. "
       "Execution will be stopped.\n";
   } else {
     log() << "This exception is too serious. Execution will be stopped.\n";
   }
   if ( !noEvent ) log() << *event;
   else log()
     << "An exception occurred before any event object was created!\n";
   dump();
   return true;
 }
 
 struct MatcherOrdering {
   bool operator()(tcPMPtr m1, tcPMPtr m2) {
     return m1->name() < m2->name() ||
       ( m1->name() == m2->name() && m1->fullName() < m2->fullName() );
   }
 };
 
 struct ObjectOrdering {
   bool operator()(tcIBPtr i1, tcIBPtr i2) {
     return i1->fullName() < i2->fullName();
   }
 };
 
 void EventGenerator::persistentOutput(PersistentOStream & os) const {
   set<tcPMPtr,MatcherOrdering> match(theMatchers.begin(), theMatchers.end());
   set<tcIBPtr,ObjectOrdering> usedset(usedObjects.begin(), usedObjects.end());
   os << theDefaultObjects << theLocalParticles << theStandardModel
      << theStrategy << theRandom << theEventHandler << theAnalysisHandlers
      << theHistogramFactory << theEventManipulator << thePath << theRunName
      << theNumberOfEvents << theObjectMap << theParticles
      << theQuickParticles << theQuickSize << match << usedset
      << ieve << weightSum << theDebugLevel << logNonDefault << printEvent
      << dumpPeriod << keepAllDumps << debugEvent
      << maxWarnings << maxErrors << theCurrentEventHandler
      << theCurrentStepHandler << useStdout << theIntermediateOutput << theMiscStream.str()
      << Repository::listReadDirs();
 }
 
 void EventGenerator::persistentInput(PersistentIStream & is, int) {
   string dummy;
   vector<string> readdirs;
   theGlobalLibraries = is.globalLibraries();
   is >> theDefaultObjects >> theLocalParticles >> theStandardModel
      >> theStrategy >> theRandom >> theEventHandler >> theAnalysisHandlers
      >> theHistogramFactory >> theEventManipulator >> thePath >> theRunName
      >> theNumberOfEvents >> theObjectMap >> theParticles
      >> theQuickParticles >> theQuickSize >> theMatchers >> usedObjects
      >> ieve >> weightSum >> theDebugLevel >> logNonDefault >> printEvent
      >> dumpPeriod >> keepAllDumps >> debugEvent
      >> maxWarnings >> maxErrors >> theCurrentEventHandler
      >> theCurrentStepHandler >> useStdout >> theIntermediateOutput >> dummy
      >> readdirs;
   theMiscStream.str(dummy);
   theMiscStream.seekp(0, std::ios::end);
   theObjects.clear();
   for ( ObjectMap::iterator it = theObjectMap.begin();
 	it != theObjectMap.end(); ++it ) theObjects.insert(it->second);
   Repository::appendReadDir(readdirs);
 }
 
 void EventGenerator::setLocalParticles(PDPtr pd, int) {
   localParticles()[pd->id()] = pd;
 }
   
 void EventGenerator::insLocalParticles(PDPtr pd, int) {
   localParticles()[pd->id()] = pd;
 }
   
 void EventGenerator::delLocalParticles(int place) {
   ParticleMap::iterator it = localParticles().begin();
   while ( place-- && it != localParticles().end() ) ++it;
   if ( it != localParticles().end() ) localParticles().erase(it);
 }
 
 vector<PDPtr> EventGenerator::getLocalParticles() const {
   vector<PDPtr> ret;
   for ( ParticleMap::const_iterator it = localParticles().begin();
  	it != localParticles().end(); ++it ) ret.push_back(it->second);
   return ret;
 }
 
 void EventGenerator::setPath(string newPath) {
   if ( std::system(("mkdir -p " + newPath).c_str()) ) throw EGNoPath(newPath);
   if ( std::system(("touch " + newPath + "/.ThePEG").c_str()) )
     throw EGNoPath(newPath);
   if ( std::system(("rm -f " + newPath + "/.ThePEG").c_str()) )
     throw EGNoPath(newPath);
   thePath = newPath;
 }
 
 string EventGenerator::defPath() const {
   char * env = std::getenv("ThePEG_RUN_DIR");
   if ( env ) return string(env);
   return string(".");
 }
 
 ostream & EventGenerator::out() {
   return theOutStream;
 }
 
 ostream & EventGenerator::log() {
   return logfile().is_open()? logfile(): BaseRepository::cout();
 }
 
 ostream & EventGenerator::ref() {
   return reffile().is_open()? reffile(): BaseRepository::cout();
 }
 
 string EventGenerator::doSaveRun(string runname) {
   runname = StringUtils::car(runname);
   if ( runname.empty() ) runname = theRunName;
   if ( runname.empty() ) runname = name();
   EGPtr eg = Repository::makeRun(this, runname);
   string file =  eg->filename() + ".run";
   PersistentOStream os(file);
   os << eg;
   if ( !os ) return "Error: Save failed! (I/O error)";
   return "";
 }
 
 string EventGenerator::doMakeRun(string runname) {
   runname = StringUtils::car(runname);
   if ( runname.empty() ) runname = theRunName;
   if ( runname.empty() ) runname = name();
   Repository::makeRun(this, runname);
   return "";
 }
 
 bool EventGenerator::preinitRegister(IPtr obj, string fullname) {
   if ( !preinitializing ) throw InitException()
     << "Tried to register a new object in the initialization of an "
     << "EventGenerator outside of the pre-initialization face. "
     << "The preinitRegister() can only be called from a doinit() function "
     << "in an object for which preInitialize() returns true.";
   if ( objectMap().find(fullname) != objectMap().end() ) return false;
   obj->name(fullname);
   objectMap()[fullname] = obj;
   objects().insert(obj);
   obj->theGenerator = this;
   PDPtr pd = dynamic_ptr_cast<PDPtr>(obj);
   if ( pd ) theParticles[pd->id()] = pd;
   PMPtr pm = dynamic_ptr_cast<PMPtr>(obj);
   if ( pm ) theMatchers.insert(pm);
   return true;
 }
 
 IPtr EventGenerator::
 preinitCreate(string classname, string fullname, string libraries) {
   if ( !preinitializing ) throw InitException()
     << "Tried to create a new object in the initialization of an "
     << "EventGenerator outside of the pre-initialization face. "
     << "The preinitCreate() can only be called from a doinit() function "
     << "in an object for which preInitialize() returns true.";
   if ( objectMap().find(fullname) != objectMap().end() ) return IPtr();
   const ClassDescriptionBase * db = DescriptionList::find(classname);
   while ( !db && libraries.length() ) {
     string library = StringUtils::car(libraries);
     libraries = StringUtils::cdr(libraries);
     DynamicLoader::load(library);
     db = DescriptionList::find(classname);
   }
   if ( !db ) return IPtr();
   IPtr obj = dynamic_ptr_cast<IPtr>(db->create());
   if ( !obj ) return IPtr();
   if ( !preinitRegister(obj, fullname) ) return IPtr();
   return obj;
 }
 
 string EventGenerator::
 preinitInterface(IPtr obj, string ifcname, string cmd, string value) {
   if ( !preinitializing ) throw InitException()
     << "Tried to manipulate an external object in the initialization of an "
     << "EventGenerator outside of the pre-initialization face. "
     << "The preinitSet() can only be called from a doinit() function "
     << "in an object for which preInitialize() returns true.";
   if ( !obj ) return "Error: No object found.";
   const InterfaceBase * ifc = Repository::FindInterface(obj, ifcname);
   if ( !ifc ) return "Error: No such interface found.";
   try {
     return ifc->exec(*obj, cmd, value);
   }
   catch ( const InterfaceException & ex) {
     ex.handle();
     return "Error: " + ex.message();
   }
 }
 
 string EventGenerator::
 preinitInterface(IPtr obj, string ifcname, int index,
 		 string cmd, string value) {
   ostringstream os;
   os << index;
   return preinitInterface(obj, ifcname, cmd, os.str() + " " + value);
 }
 
 string EventGenerator::
 preinitInterface(string fullname, string ifcname, string cmd, string value) {
   return preinitInterface(getObject<Interfaced>(fullname), ifcname, cmd, value);
 }
 
 string EventGenerator::
 preinitInterface(string fullname, string ifcname, int index,
 		 string cmd, string value) {
   return preinitInterface(getObject<Interfaced>(fullname), ifcname, index,
 			  cmd, value);
 }
 
 tDMPtr EventGenerator::findDecayMode(string tag) const {
   for ( ObjectSet::const_iterator it = objects().begin();
 	it != objects().end(); ++it ) {
     tDMPtr dm = dynamic_ptr_cast<tDMPtr>(*it);
     if ( dm && dm->tag() == tag ) return dm;
   }
   return tDMPtr();
 }
 
 tDMPtr EventGenerator::preinitCreateDecayMode(string tag) {
   return constructDecayMode(tag);
 }
 
 DMPtr EventGenerator::constructDecayMode(string & tag) {
   DMPtr rdm;
   DMPtr adm;
   int level = 0;
   string::size_type end = 0;
   while ( end < tag.size() && ( tag[end] != ']' || level ) ) {
     switch ( tag[end++] ) {
     case '[':
       ++level;
       break;
     case ']':
       --level;
       break;
     }
   }
   rdm = findDecayMode(tag.substr(0,end));
   if ( rdm ) return rdm;
 
   string::size_type next = tag.find("->");
   if ( next == string::npos ) return rdm;
   if ( tag.find(';') == string::npos ) return rdm;
   tPDPtr pd = getObject<ParticleData>(tag.substr(0,next));
   if ( !pd ) pd = findParticle(tag.substr(0,next));
   if ( !pd ) return rdm;
   rdm = ptr_new<DMPtr>();
   rdm->parent(pd);
   if ( pd->CC() ) {
     adm = ptr_new<DMPtr>();
     adm->parent(pd->CC());
     rdm->theAntiPartner = adm;
     adm->theAntiPartner = rdm;
   }
   bool error = false;
   tag = tag.substr(next+2);
   tPDPtr lastprod;
   bool dolink = false;
   do {
     switch ( tag[0] ) {
     case '[':
       {
 	tag = tag.substr(1);
 	tDMPtr cdm = constructDecayMode(tag);
 	if ( cdm ) rdm->addCascadeProduct(cdm);
 	else error = true;
       } break;
     case '=':
       dolink = true;
+      [[fallthrough]];
     case ',':
     case ']':
       tag = tag.substr(1);
       break;
     case '?':
       {
 	next = min(tag.find(','), tag.find(';'));
 	tPMPtr pm = findMatcher(tag.substr(1,next-1));
 	if ( pm ) rdm->addProductMatcher(pm);
 	else error = true;
 	tag = tag.substr(next);
       } break;
     case '!':
       {
 	next = min(tag.find(','), tag.find(';'));
 	tPDPtr pd = findParticle(tag.substr(1,next-1));
 	if ( pd ) rdm->addExcluded(pd);
 	else error = true;
 	tag = tag.substr(next);
       } break;
     case '*':
       {
 	next = min(tag.find(','), tag.find(';'));
 	tPMPtr pm = findMatcher(tag.substr(1,next-1));
 	if ( pm ) rdm->setWildMatcher(pm);
 	else error = true;
 	tag = tag.substr(next);
       } break;
     default:
       {
 	next = min(tag.find('='), min(tag.find(','), tag.find(';')));
 	tPDPtr pdp = findParticle(tag.substr(0,next));
 	if ( pdp ) rdm->addProduct(pdp);
 	else error = true;
 	tag = tag.substr(next);
 	if ( dolink && lastprod ) {
 	  rdm->addLink(lastprod, pdp);
 	  dolink = false;
 	}
 	lastprod = pdp;
       } break;
     }
   } while ( tag[0] != ';' && tag.size() );
   if ( tag[0] != ';' || error ) {
     return DMPtr();
   }
 
   tag = tag.substr(1);
   
   DMPtr ndm = findDecayMode(rdm->tag());
   if ( ndm ) return ndm;
   pd->addDecayMode(rdm);
   if ( !preinitRegister(rdm, pd->fullName() + "/" + rdm->tag()) )
     return DMPtr();
   if ( adm ) {
     preinitRegister(adm, pd->CC()->fullName() + "/" + adm->tag());
     rdm->CC(adm);
     adm->CC(rdm);
   }
 
   return rdm;
 }
 
 tPDPtr EventGenerator::findParticle(string pdgname) const {
   for ( ParticleMap::const_iterator it = particles().begin();
 	it != particles().end(); ++it )
     if ( it->second->PDGName() == pdgname ) return it->second;
   return tPDPtr();
 }
 
 tPMPtr EventGenerator::findMatcher(string name) const {
   for ( MatcherSet::const_iterator it = matchers().begin();
 	it != matchers().end(); ++it )
     if ( (**it).name() == name ) return *it;
   return tPMPtr();
 }
 
 ClassDescription<EventGenerator> EventGenerator::initEventGenerator;
 
 void EventGenerator::Init() {
   
   static ClassDocumentation<EventGenerator> documentation
     ("This is the main class used to administer an event generation run. "
      "The actual generation of each event is handled by the assigned "
      "<interface>EventHandler</interface> object. When the event generator"
      "is properly set up it can be initialized with the command "
      "<interface>MakeRun</interface> and/or saved to a file with the command "
      "<interface>SaveRun</interface>. If saved to a file, the event generator "
      "can be read into another program to produce events. The file can also "
      "be read into the <tt>runThePEG</tt> program where a number of events "
      "determined by the parameter <interface>NumberOfEvents</interface> is "
      "generated with each event analysed by the list of assigned "
      "<interface>AnalysisHandlers</interface>.");
 
   static Reference<EventGenerator,StandardModelBase> interfaceStandardModel
     ("StandardModelParameters",
      "The ThePEG::StandardModelBase object to be used to access standard "
      "model parameters in this run.",
      &EventGenerator::theStandardModel, false, false, true, false);
 
   static Reference<EventGenerator,EventHandler> interfaceEventHandler
     ("EventHandler",
      "The ThePEG::EventHandler object to be used to generate the "
      "individual events in this run.",
      &EventGenerator::theEventHandler, false, false, true, false);
 
   static RefVector<EventGenerator,AnalysisHandler> interfaceAnalysisHandlers
     ("AnalysisHandlers",
      "ThePEG::AnalysisHandler objects to be used to analyze the produced "
      "events in this run.",
      &EventGenerator::theAnalysisHandlers, 0, true, false, true, false);
 
   static Reference<EventGenerator,FactoryBase> interfaceHistogramFactory
     ("HistogramFactory",
      "An associated factory object for handling histograms to be used by "
      "<interface>AnalysisHandlers</interface>.",
      &EventGenerator::theHistogramFactory, true, false, true, true, true);
 
   static Reference<EventGenerator,EventManipulator> interfaceEventManip
     ("EventManipulator",
      "An ThePEG::EventManipulator called each time the generation of an "
      "event is stopped. The ThePEG::EventManipulator object is able to "
      "manipulate the generated event, as opposed to an "
      "ThePEG::AnalysisHandler which may only look at the event.",
      &EventGenerator::theEventManipulator, true, false, true, true);
 
   static RefVector<EventGenerator,ParticleData> interfaceLocalParticles
     ("LocalParticles",
      "Special versions of ThePEG::ParticleData objects to be used "
      "in this run. Note that to delete an object, its number in the list "
      "should be given, rather than its id number.",
      0, 0, false, false, true, false,
      &EventGenerator::setLocalParticles, &EventGenerator::insLocalParticles,
      &EventGenerator::delLocalParticles, &EventGenerator::getLocalParticles);
 
   static RefVector<EventGenerator,Interfaced> interfaceDefaultObjects
     ("DefaultObjects",
      "A vector of pointers to default objects. In a ThePEG::Reference or "
      "ThePEG::RefVector interface with the defaultIfNull() flag set, if a "
      "null pointer is encountered this vector is gone through until an "
      "acceptable object is found in which case the null pointer is replaced "
      "by a pointer to this object.",
      &EventGenerator::theDefaultObjects, 0, true, false, true, false, false);
 
   static Reference<EventGenerator,Strategy> interfaceStrategy
     ("Strategy",
      "An ThePEG::Strategy with additional ThePEG::ParticleData objects to "
      "be used in this run.",
      &EventGenerator::theStrategy, false, false, true, true);
 
   static Reference<EventGenerator,RandomGenerator> interfaceRandomGenerator
     ("RandomNumberGenerator",
      "An ThePEG::RandomGenerator object which should typically interaface to "
      "a CLHEP Random object. This will be the default random number generator "
      "for the run, but individual objects may use their own random generator "
      "if they wish.",
      &EventGenerator::theRandom, true, false, true, false);
 
   static Parameter<EventGenerator,string> interfacePath
     ("Path",
      "The directory where the output files are put.",
      &EventGenerator::thePath, ".", true, false,
       &EventGenerator::setPath, 0, &EventGenerator::defPath);
   interfacePath.directoryType();
 
   static Parameter<EventGenerator,string> interfaceRunName
     ("RunName",
      "The name of this run. This name will be used in the output filenames. "
      "The files wil be placed in the directory specified by the "
      "<interface>Path</interface> parameter"
      "If empty the name of the event generator will be used instead.",
      &EventGenerator::theRunName, "", true, false,
      0, 0, &EventGenerator::name);
 
   static Parameter<EventGenerator,long> interfaceNumberOfEvents
     ("NumberOfEvents",
      "The number of events to be generated in this run. If less than zero, "
      "the number of events is unlimited",
      &EventGenerator::theNumberOfEvents, 1000, -1, Constants::MaxInt,
      true, false, Interface::lowerlim);
 
   static Parameter<EventGenerator,int> interfaceDebugLevel
     ("DebugLevel",
      "The level of debug information sent out to the log file in the run. "
      "Level 0 only gives a limited ammount of warnings and error messages. "
      "Level 1 will print the first few events. "
      "Level 5 will print every event. "
      "Level 9 will print every step in every event.",
      &EventGenerator::theDebugLevel, 0, 0, 9, true, false, true);
 
   static Parameter<EventGenerator,int> interfacePrintEvent
     ("PrintEvent",
      "If the debug level is above zero, print the first 'PrintEvent' events.",
      &EventGenerator::printEvent, 0, 0, 1000, true, false, Interface::lowerlim);
 
   static Parameter<EventGenerator,long> interfaceDumpPeriod
     ("DumpPeriod",
      "If the debug level is above zero, dump the full state of the run every "
      "'DumpPeriod' events. Set it to -1 to disable dumping even in the case of errors.",
      &EventGenerator::dumpPeriod, 0, -1, Constants::MaxInt,
      true, false, Interface::lowerlim);
 
   static Switch<EventGenerator,bool> interfaceKeepAllDumps
     ("KeepAllDumps",
      "Whether all dump files should be kept, labelled by event number.",
      &EventGenerator::keepAllDumps, false, true, false);
   static SwitchOption interfaceKeepAllDumpsYes
     (interfaceKeepAllDumps,
      "Yes",
      "Keep all dump files, labelled by event number.",
      true);
   static SwitchOption interfaceKeepAllDumpsNo
     (interfaceKeepAllDumps,
      "No",
      "Keep only the latest dump file.",
      false);
 
   static Parameter<EventGenerator,long> interfaceDebugEvent
     ("DebugEvent",
      "If the debug level is above zero, step up to the highest debug level "
      "befor event number 'DebugEvent'.",
      &EventGenerator::debugEvent, 0, 0, Constants::MaxInt,
      true, false, Interface::lowerlim);
 
   static Parameter<EventGenerator,int> interfaceMaxWarnings
     ("MaxWarnings",
      "The maximum number of warnings of each type which will be printed.",
      &EventGenerator::maxWarnings,
      10, 1, 100, true, false, Interface::lowerlim);
 
   static Parameter<EventGenerator,int> interfaceMaxErrors
     ("MaxErrors",
      "The maximum number of errors of each type which will be tolerated. "
      "If more errors are reported, the run will be aborted.",
      &EventGenerator::maxErrors,
      10, -1, 100000, true, false, Interface::lowerlim);
 
   static Parameter<EventGenerator,long> interfaceQuickSize
     ("QuickSize",
      "The max absolute id number of particle data objects which are accessed "
      "quickly through a vector indexed by the id number.",
      &EventGenerator::theQuickSize,
      7000, 0, 50000, true, false, Interface::lowerlim);
 
 
   static Command<EventGenerator> interfaceSaveRun
     ("SaveRun",
      "Isolate, initialize and save this event generator to a file, from which "
      "it can be read in and run in another program. If an agument is given "
      "this is used as the run name, otherwise the run name is taken from the "
      "<interface>RunName</interface> parameter.",
      &EventGenerator::doSaveRun, true);
 
 
   static Command<EventGenerator> interfaceMakeRun
     ("MakeRun",
      "Isolate and initialize this event generator and give it a run name. "
      "If no argument is given, the run name is taken from the "
      "<interface>RunName</interface> parameter.",
      &EventGenerator::doMakeRun, true);
 
   interfaceEventHandler.rank(11.0);
   interfaceSaveRun.rank(10.0);
   interfaceMakeRun.rank(9.0);
   interfaceRunName.rank(8.0);
   interfaceNumberOfEvents.rank(7.0);
   interfaceAnalysisHandlers.rank(6.0);
 
 
   static Switch<EventGenerator,bool> interfaceUseStdout
     ("UseStdout",
      "Redirect the logging and output to stdout instead of files.",
      &EventGenerator::useStdout, false, true, false);
   static SwitchOption interfaceUseStdoutYes
     (interfaceUseStdout,
      "Yes",
      "Use stdout instead of log files.",
      true);
   static SwitchOption interfaceUseStdoutNo
     (interfaceUseStdout,
      "No",
      "Use log files.",
      false);
   
 
   static Switch<EventGenerator,int> interfaceLogNonDefault
     ("LogNonDefault",
      "Controls the printout of important interfaces which has been changed from their default values.",
      &EventGenerator::logNonDefault, -1, true, false);
   static SwitchOption interfaceLogNonDefaultYes
     (interfaceLogNonDefault,
      "Yes",
      "Always print changed interfaces.",
      1);
   static SwitchOption interfaceLogNonDefaultOnDebug
     (interfaceLogNonDefault,
      "OnDebug",
      "Only print changed interfaces if debugging is turned on.",
      0);
   static SwitchOption interfaceLogNonDefaultNo
     (interfaceLogNonDefault,
      "No",
      "Don't print changed interfaces.",
      -1);
   interfaceLogNonDefault.setHasDefault(false);
 
   static Switch<EventGenerator,bool> interfaceIntermediateOutput
     ("IntermediateOutput",
      "Modified event number count with the number of events processed so far, "
      "which updates at least every 10000 events, together with the corresponding "
      "intermediate estimate for the cross section plus the integration error.",
      &EventGenerator::theIntermediateOutput, false, true, false);
   static SwitchOption interfaceIntermediateOutputYes
     (interfaceIntermediateOutput,
      "Yes",
      "Show the modified event number count with the number of events processed so far, "
      "plus further information on the intermediate cross section estimate.",
      true);
   static SwitchOption interfaceIntermediateOutputNo
     (interfaceIntermediateOutput,
      "No",
      "Show the usual event number count with the number of events processed so far, "
      "but no further information on the intermediate cross section estimate.",
      false);
 
 }
 
 EGNoPath::EGNoPath(string path) {
   theMessage << "Cannot set the directory path for output files to '" << path
 	     << "' because the directory did not exist and could not be "
 	     << "created.";
   severity(warning);
 }