diff --git a/MatrixElement/DiagramBase.cc b/MatrixElement/DiagramBase.cc
--- a/MatrixElement/DiagramBase.cc
+++ b/MatrixElement/DiagramBase.cc
@@ -1,73 +1,73 @@
 // -*- C++ -*-
 //
 // DiagramBase.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 DiagramBase class.
 //
 
 #include "DiagramBase.h"
 #include "ThePEG/Utilities/ClassDescription.h"
 #include "ThePEG/Utilities/DescriptionList.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace ThePEG;
 
 DiagramBase::~DiagramBase() {}
 
 struct ParticleOrdering {
-  bool operator()(tcPDPtr p1, tcPDPtr p2) {
+  bool operator()(tcPDPtr p1, tcPDPtr p2) const {
     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() );
   }
 };
 
 string DiagramBase::getTag() const {
   string tag;
   if ( !done() ) throw DiagramBaseSetupException(*this);
   for ( int i = 0; i < nIncoming(); ++i ) {
     if ( i ) tag += ",";
     tag += partons()[i]->PDGName();
   }
   tag += "->";
 
   multiset<tcPDPtr,ParticleOrdering> out;
   for ( int i = nIncoming(), N = partons().size(); i < N; ++i )
     out.insert(partons()[i]);
   for ( multiset<tcPDPtr,ParticleOrdering>::iterator i = out.begin();
 	i != out.end(); ++i ) {
     if ( i != out.begin() ) tag += ",";
     tag += (**i).PDGName();
   }
   return tag;
 }
 
 void DiagramBase::persistentOutput(PersistentOStream & os) const {
   os << theNIncoming << thePartons << theId;
 }
 
 void DiagramBase::persistentInput(PersistentIStream & is, int) {
   is >> theNIncoming >> thePartons >> theId;
 }
 
 AbstractClassDescription<DiagramBase> DiagramBase::initDiagramBase;
 // Definition of the static class description member.
 
 void DiagramBase::Init() {}
 
 DiagramBaseSetupException::DiagramBaseSetupException(const DiagramBase & db) {
   const ClassDescriptionBase * cd = DescriptionList::find(typeid(db));
   if ( !cd )
     theMessage << "Tried to use an unknown sub class of DiagramBase.";
   else
     theMessage << "The '" << cd->name()
 	       << "' sub class did not setup the DiagramBase class correctly.";
   severity(abortnow);
 }
diff --git a/PDT/DecayMode.cc b/PDT/DecayMode.cc
--- a/PDT/DecayMode.cc
+++ b/PDT/DecayMode.cc
@@ -1,692 +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 operator()(tcPDPtr p1, tcPDPtr p2) const { 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) {
+  bool operator()(tcPDPtr p1, tcPDPtr p2) const {
     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) {
+  bool operator()(tcDMPtr d1, tcDMPtr d2) const {
     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) {
+  bool operator()(tcPMPtr m1, tcPMPtr m2) const {
     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/MatcherBase.cc b/PDT/MatcherBase.cc
--- a/PDT/MatcherBase.cc
+++ b/PDT/MatcherBase.cc
@@ -1,201 +1,201 @@
 // -*- C++ -*-
 //
 // MatcherBase.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 MatcherBase class.
 //
 
 #include "MatcherBase.h"
 #include "DecayMode.h"
 #include "ThePEG/Utilities/Rebinder.h"
 #include "ThePEG/Repository/Repository.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Utilities/EnumIO.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 
 using namespace ThePEG;
 
 MatcherBase::MatcherBase()
   : theMaxMass(ZERO), theMinMass(ZERO), commonMass(-1.0*GeV),
     commonWidth(-1.0*GeV), commonCTau(-1.0*mm),
     commonCharge(PDT::ChargeUndefined), commonSpin(PDT::SpinUndefined),
     commonColour(PDT::ColourUndefined), commonStable(-1) {}
 
 MatcherBase::MatcherBase(const MatcherBase & m)
   : Interfaced(m), theMaxMass(m.theMaxMass), theMinMass(m.theMinMass),
     commonMass(m.commonMass), commonWidth(m.commonWidth),
     commonCTau(m.commonCTau), commonCharge(m.commonCharge),
     commonSpin(m.commonSpin), commonColour(m.commonColour),
     commonStable(m.commonStable), theAntiPartner(m.theAntiPartner) {}
 
 MatcherBase::~MatcherBase() {}
 
 void MatcherBase::doupdate() {
   Interfaced::doupdate();
   tPDSet oldParticles;
   tPMSet oldMatchers;
   Energy oldMass = commonMass;
   Energy oldWidth = commonWidth;
   Length oldCTau = commonCTau;
   PDT::Spin  oldSpin = commonSpin;
   PDT::Charge oldCharge = commonCharge;
   PDT::Colour oldColour = commonColour;
   int oldStable = commonStable;
   matchingParticles.swap(oldParticles);
   matchingMatchers.swap(oldMatchers);
   if ( generator() ) {
     for ( ParticleMap::const_iterator it = generator()->particles().begin();
 	  it != generator()->particles().end(); ++it ) addPIfMatch(it->second);
     addMIfMatchFrom(generator()->matchers());
   } else {
     addPIfMatchFrom(Repository::allParticles());
     addMIfMatchFrom(Repository::allMatchers());
   }
   if ( matchingParticles != oldParticles || matchingMatchers != oldMatchers ||
     oldMass != commonMass || oldWidth != commonWidth || oldCTau != commonCTau ||
     oldSpin != commonSpin || oldCharge != commonCharge ||
        oldColour != commonColour || oldStable != commonStable ) touch();
 }
 
 void MatcherBase::clear() {
   matchingParticles.clear();
   matchingMatchers.clear();
   theMaxMass = ZERO;
   theMinMass = ZERO;
   commonMass = -1.0*GeV;
   commonWidth = -1.0*GeV;
   commonCTau = -1.0*mm;
   commonSpin = PDT::SpinUndefined;
   commonCharge = PDT::ChargeUndefined;
   commonColour = PDT::ColourUndefined;
   commonStable = -1;
 }
 
 void MatcherBase::addMIfMatch(tPMPtr pm) {
   if ( member(matchingMatchers, pm) ) return;
   pm->update();
   tPDSet::const_iterator i = pm->matchingParticles.begin();
   while ( i != pm->matchingParticles.end() )
     if ( !member(matchingParticles, *i++) ) return;
   matchingMatchers.insert(pm);
 }
 
 
 void MatcherBase::addPIfMatch(tPDPtr pd) {
   if ( !pd || !check(*pd) || member(matchingParticles, pd) ) return;
   if ( matchingParticles.empty() ) {
     commonMass = pd->mass();
     theMinMass = pd->mass();
     theMaxMass = pd->mass();
     commonWidth = pd->width();
     commonCTau = pd->cTau();
     commonCharge = pd->iCharge();
     commonSpin = pd->iSpin();
     commonColour = pd->iColour();
     commonStable = pd->stable();
   } else {
     if ( commonMass != pd->mass() ) commonMass = -1.0*GeV;
     theMinMass = min(theMinMass, pd->mass());
     theMaxMass = min(theMaxMass, pd->mass());
     if ( commonWidth != pd->width() ) commonWidth = -1.0*GeV;
     if ( commonCTau != pd->cTau() ) commonCTau = -1.0*mm;
     if ( commonCharge != pd->iCharge() ) {
       switch ( commonCharge ) {
       case PDT::ChargeUndefined:
 	break;
       case PDT::Positive:
 	if ( PDT::negative(pd->iCharge()) ) commonCharge = PDT::Charged;
 	else if ( !PDT::positive(pd->iCharge()) )
 	  commonCharge = PDT::ChargeUndefined;
 	break;
       case PDT::Negative:
 	if ( PDT::positive(pd->iCharge()) ) commonCharge = PDT::Charged;
 	else if ( !PDT::negative(pd->iCharge()) )
 	  commonCharge = PDT::ChargeUndefined;
 	break;
       case PDT::Charged:
 	if ( !PDT::charged(pd->iCharge()) )
 	  commonCharge = PDT::ChargeUndefined;
 	break;
       default:
 	if ( PDT::positive(commonCharge) ) {
 	  if ( PDT::positive(pd->iCharge()) )
 	    commonCharge = PDT::Positive;
 	  else if ( PDT::negative(pd->iCharge()) )
 	    commonCharge = PDT::Charged;
 	  else commonCharge = PDT::ChargeUndefined;
 	}
 	else if ( PDT::negative(commonCharge) ) {
 	  if ( PDT::negative(pd->iCharge()) )
 	    commonCharge = PDT::Negative;
 	  else if ( PDT::positive(pd->iCharge()) )
 	    commonCharge = PDT::Charged;
 	  else commonCharge = PDT::ChargeUndefined;
 	}
 	else commonCharge = PDT::ChargeUndefined;
       }
     }
     if ( commonSpin != pd->iSpin() ) commonSpin = PDT::SpinUndefined;
     if ( commonColour != pd->iColour() ) {
       if ( PDT::coloured(commonColour) && PDT::coloured(pd->iColour()) )
 	commonColour = PDT::Coloured;
       else
 	commonColour = PDT::ColourUndefined;
       if ( commonStable != pd->stable() ) commonStable = -1;
     }
   }
   matchingParticles.insert(pd);
 }
 
 struct ParticleOrdering {
-  bool operator()(tcPDPtr p1, tcPDPtr p2) {
+  bool operator()(tcPDPtr p1, tcPDPtr p2) const {
     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 MatcherOrdering {
-  bool operator()(tcPMPtr m1, tcPMPtr m2) {
+  bool operator()(tcPMPtr m1, tcPMPtr m2) const {
     return m1->name() < m2->name() ||
       ( m1->name() == m2->name() && m1->fullName() < m2->fullName() );
   }
 };
 
 void MatcherBase::persistentOutput(PersistentOStream & os ) const {
   multiset<tcPDPtr,ParticleOrdering>
     parts(particles().begin(), particles().end());
   multiset<tcPMPtr,MatcherOrdering>
     match(matchers().begin(), matchers().end());
 
   os << parts << match << ounit(theMaxMass, GeV) << ounit(theMinMass, GeV)
      << ounit(commonMass, GeV) << ounit(commonWidth, GeV)
      << ounit(commonCTau, mm) << oenum(commonCharge) << oenum(commonSpin)
      << oenum(commonColour) << commonStable << theAntiPartner;
 }
 
 void MatcherBase::persistentInput(PersistentIStream & is, int) {
   is >> matchingParticles >> matchingMatchers >> iunit(theMaxMass, GeV)
      >> iunit(theMinMass, GeV) >> iunit(commonMass, GeV)
      >> iunit(commonWidth, GeV) >> iunit(commonCTau, mm) >> ienum(commonCharge)
      >> ienum(commonSpin) >> ienum(commonColour) >> commonStable
      >> theAntiPartner;
 }
 
 AbstractClassDescription<MatcherBase> MatcherBase::initMatcherBase;
 
 void MatcherBase::Init() {
 
   static ClassDocumentation<MatcherBase> documentation
     ("This is the base class for objects representing groups of particle "
      "types.");
 
 }
 
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 & 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) {
+  bool operator()(tcPDPtr p1, tcPDPtr p2) const {
     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) {
+  bool operator()(const tcDMPtr & d1, const tcDMPtr & d2) const  {
     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/Repository/EventGenerator.cc b/Repository/EventGenerator.cc
--- a/Repository/EventGenerator.cc
+++ b/Repository/EventGenerator.cc
@@ -1,1401 +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) {
+  bool operator()(tcPMPtr m1, tcPMPtr m2) const {
     return m1->name() < m2->name() ||
       ( m1->name() == m2->name() && m1->fullName() < m2->fullName() );
   }
 };
 
 struct ObjectOrdering {
-  bool operator()(tcIBPtr i1, tcIBPtr i2) {
+  bool operator()(tcIBPtr i1, tcIBPtr i2) const {
     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);
 }
 
diff --git a/Repository/Repository.cc b/Repository/Repository.cc
--- a/Repository/Repository.cc
+++ b/Repository/Repository.cc
@@ -1,1150 +1,1150 @@
 // -*- C++ -*-
 //
 // Repository.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 Repository class.
 //
 
 // macro is passed in from -D compile flag
 #ifndef THEPEG_PKGLIBDIR
 #error Makefile.am needs to define THEPEG_PKGLIBDIR
 #endif
 
 #include "Repository.h"
 #include "ThePEG/Utilities/Rebinder.h"
 #include "ThePEG/Handlers/EventHandler.h"
 #include "ThePEG/PDT/DecayMode.h"
 #include "ThePEG/Repository/Strategy.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Utilities/Debug.h"
 #include "ThePEG/Config/algorithm.h"
 #include "ThePEG/Utilities/DynamicLoader.h"
 #include "ThePEG/Utilities/StringUtils.h"
 
 #include <iterator>
 #include <chrono>
 
 #include <config.h>
 
 // readline options taken from
 // http://autoconf-archive.cryp.to/vl_lib_readline.html 
 // Copyright © 2008 Ville Laurikari <vl@iki.fi> 
 // Copying and distribution of this file, with or without
 // modification, are permitted in any medium without royalty provided
 // the copyright notice and this notice are preserved.
 
 #ifdef HAVE_LIBREADLINE
 #  if defined(HAVE_READLINE_READLINE_H)
 #    include <readline/readline.h>
 #  elif defined(HAVE_READLINE_H)
 #    include <readline.h>
 #  else
 extern "C" char *readline (const char *);
 #  endif
 #endif
 
 #ifdef HAVE_READLINE_HISTORY
 #  if defined(HAVE_READLINE_HISTORY_H)
 #    include <readline/history.h>
 #  elif defined(HAVE_HISTORY_H)
 #    include <history.h>
 #  else
 extern "C" void add_history (const char *);
 #  endif
 #endif
 
 using namespace ThePEG;
 
 ParticleMap & Repository::defaultParticles() {
   static ParticleMap theMap;
   return theMap;
 }
 
 
 ParticleDataSet & Repository::particles() {
   static ParticleDataSet theSet;
   return theSet;
 }
 
 MatcherSet & Repository::matchers() {
   static MatcherSet theSet;
   return theSet;
 }
 
 Repository::GeneratorMap & Repository::generators() {
   static GeneratorMap theMap;;
   return theMap;
 }
 
 string & Repository::currentFileName() {
   static string theCurrentFileName;
   return theCurrentFileName;
 }
 
 int & Repository::exitOnError() {
   static int exitonerror = 0;
   return exitonerror;
 }
 
 void Repository::cleanup() {
   generators().clear();
 }
 
 void Repository::Register(IBPtr ip) {
   BaseRepository::Register(ip);
   registerParticle(dynamic_ptr_cast<PDPtr>(ip));
   registerMatcher(dynamic_ptr_cast<PMPtr>(ip));
 }
 
 void Repository::Register(IBPtr ip, string newName) {
   DirectoryAppend(newName);
   BaseRepository::Register(ip, newName);
   registerParticle(dynamic_ptr_cast<PDPtr>(ip));
   registerMatcher(dynamic_ptr_cast<PMPtr>(ip));
 }
 
 void Repository::registerParticle(tPDPtr pd) {
   if ( !pd ) return;
   if ( !member(particles(), pd) ) {
     particles().insert(pd);
     CreateDirectory(pd->fullName());
   }
   if ( pd->id() == 0 ) return;
   if ( !member(defaultParticles(), pd->id()) )
     defaultParticles()[pd->id()] = pd;
   for ( MatcherSet::iterator it = matchers().begin();
 	it != matchers().end(); ++it) (*it)->addPIfMatch(pd);
 }
 
 void Repository::registerMatcher(tPMPtr pm) {
   if ( !pm || member(matchers(), pm) ) return;
   pm->addPIfMatchFrom(particles());
   for ( MatcherSet::iterator it = matchers().begin();
 	it != matchers().end(); ++it) {
     (*it)->addMIfMatch(pm);
     pm->addMIfMatch(*it);
   }
   matchers().insert(pm);
 }
 
 tPDPtr Repository::findParticle(string name) {
   tPDPtr pd;
   string path = name;
   DirectoryAppend(path);
   pd = dynamic_ptr_cast<tPDPtr>(GetPointer(path));
   if ( pd ) return pd;
   for ( ParticleMap::iterator pit = defaultParticles().begin();
 	pit != defaultParticles().end(); ++pit )
     if ( pit->second->PDGName() == name ) return pit->second;
   for ( ParticleDataSet::iterator pit = particles().begin();
 	pit != particles().end(); ++pit )
     if ( (**pit).PDGName() == name ) return *pit;
   return pd;
 }
 
 tPMPtr Repository::findMatcher(string name) {
   for ( MatcherSet::iterator mit = matchers().begin();
 	mit != matchers().end(); ++mit )
     if ( name == (**mit).name() ) return *mit;
   return tPMPtr();
 }
 
 void Repository::saveRun(string EGname, string name, string filename) {
   EGPtr eg = BaseRepository::GetObject<EGPtr>(EGname);
   EGPtr run = makeRun(eg, name);
   PersistentOStream os(filename, globalLibraries());
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "Saving event generator '" << name << "'... " << flush;
   os << run;
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "done" << endl;
 }
 
 EGPtr Repository::makeRun(tEGPtr eg, string name) {
 
   // Clone all objects relevant for the EventGenerator. This is
   // the EventGenerator itself, all particles and all particle
   // matchers. 'localObject' is the set of all object refered to by
   // the generator particles and matcher and in the end these are
   // cloned as well.
 
   // Clone all Particle matchers
 
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "Making event generator '" << name << "':" << endl
 	   << "Updating all objects... " << flush;
 
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "done\nCloning matchers and particles... " << flush;
 
   MatcherSet localMatchers;
   ObjectSet localObjects;
   ObjectSet clonedObjects;
   TranslationMap trans;
 
   for ( MatcherSet::iterator mit = matchers().begin();
 	mit != matchers().end(); ++mit ) {
     PMPtr pm = clone(**mit);
     pm->clear();
     trans[*mit] = pm;
     localMatchers.insert(pm);
     clonedObjects.insert(pm);
     localObjects.insert(*mit);
     addReferences(*mit, localObjects);
   }
 
 
   // Clone the particles. But only the ones which should be
   // used. First select the localParticles of the EventGenerator, then
   // add particles from the strategy of the EventGenerator which have
   // not already been selected. Finally add particles from the global
   // default if no default directories has been specified in the
   // strategy which have not already been selected.
   PDVector allParticles;
 
   for ( ParticleMap::const_iterator pit = eg->localParticles().begin();
  	pit != eg->localParticles().end(); ++pit )
     allParticles.push_back(pit->second);
   if ( eg->strategy() ) {
     tcStrategyPtr strat = eg->strategy();
     for ( ParticleMap::const_iterator pit = strat->particles().begin();
  	  pit != strat->particles().end(); ++pit )
       allParticles.push_back(pit->second);
 
     vector<string> pdirs;
     if ( eg->strategy()->localParticlesDir().length() )
       pdirs.push_back(eg->strategy()->localParticlesDir());
     pdirs.insert(pdirs.end(), eg->strategy()->defaultParticlesDirs().begin(),
 		 eg->strategy()->defaultParticlesDirs().end());
 
     for ( int i = 0, N = pdirs.size(); i < N; ++i ) {
       string dir = pdirs[i];
       for ( ParticleDataSet::iterator pit = particles().begin();
 	    pit != particles().end(); ++pit )
 	if ( (**pit).fullName().substr(0, dir.length()) == dir )
 	  allParticles.push_back(*pit);
     }
   }
 
   if ( !eg->strategy() || eg->strategy()->defaultParticlesDirs().empty() )
     for ( ParticleMap::iterator pit = defaultParticles().begin();
 	  pit != defaultParticles().end(); ++pit )
       allParticles.push_back(pit->second);
 
   for ( ParticleDataSet::iterator pit = particles().begin();
 	pit != particles().end(); ++pit )
     allParticles.push_back(*pit);
 
   ParticleMap localParticles;
 
   for ( PDVector::iterator pit = allParticles.begin();
 	pit != allParticles.end(); ++pit ) {
     ParticleMap::iterator it = localParticles.find((**pit).id());
     if ( it == localParticles.end() ) {
       PDPtr pd = clone(**pit);
       trans[*pit] = pd;
       localParticles[pd->id()] = pd;
       clonedObjects.insert(pd);
       localObjects.insert(*pit);
       addReferences(*pit, localObjects);
     } else {
       trans[*pit] = it->second;
     }
   }
 
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "done\nCloning other objects... " << flush;
 
   // Clone the OldEventGenerator object to be used:
   localObjects.insert(eg);
   addReferences(eg, localObjects);
   EGPtr egrun = clone(*eg);
   clonedObjects.insert(egrun);
   trans[eg] = egrun;
 
   for ( ObjectSet::iterator it = localObjects.begin();
 	it != localObjects.end(); ++it ) {
     if ( member(trans.map(), *it) ) continue;
     IBPtr ip = clone(**it);
     trans[*it] = ip;
     clonedObjects.insert(ip);
   }
 
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "done\nRebind references... " << flush;
 
   IVector defaults;
 
   trans.translate(inserter(defaults), eg->defaultObjects().begin(),
 		  eg->defaultObjects().end());
   if ( eg->strategy() )
     trans.translate(inserter(defaults),
 		    eg->strategy()->defaultObjects().begin(),
 		    eg->strategy()->defaultObjects().end());
 
   for ( ObjectSet::iterator it = clonedObjects.begin();
 	it != clonedObjects.end(); ++it ) {
     dynamic_cast<Interfaced &>(**it).theGenerator = egrun;
     rebind(**it, trans, defaults);
   }
 
   // Now, dependencies may have changed, so we do a final round of
   // updates.
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "done\nUpdating cloned objects... " << flush;
 
 
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "done\nInitializing... " << flush;
 
   clonedObjects.erase(egrun);
   egrun->setup(name, clonedObjects, localParticles, localMatchers);
 
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "done" << endl;
 
   generators()[name] = egrun;
 
   return egrun;
 
 }
 
 PDPtr Repository::defaultParticle(PID id) {
   ParticleMap::iterator pit = defaultParticles().find(id);
   return pit == defaultParticles().end()? PDPtr(): pit->second;
 }
 
 void Repository::defaultParticle(tPDPtr pdp) {
   if ( pdp ) defaultParticles()[pdp->id()] = pdp;
 }
 
 struct ParticleOrdering {
-  bool operator()(tcPDPtr p1, tcPDPtr p2) {
+  bool operator()(tcPDPtr p1, tcPDPtr p2) const {
     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 MatcherOrdering {
-  bool operator()(tcPMPtr m1, tcPMPtr m2) {
+  bool operator()(tcPMPtr m1, tcPMPtr m2) const {
     return m1->name() < m2->name() ||
       ( m1->name() == m2->name() && m1->fullName() < m2->fullName() );
   }
 };
 
 struct InterfaceOrdering {
-  bool operator()(tcIBPtr i1, tcIBPtr i2) {
+  bool operator()(tcIBPtr i1, tcIBPtr i2) const {
     return i1->fullName() < i2->fullName();
   }
 };
 
 void Repository::save(string filename) {
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "saving '" << filename << "'... " << flush;
   PersistentOStream os(filename, globalLibraries());
   set<tcPDPtr,ParticleOrdering>
     part(particles().begin(), particles().end());
   set<tcPMPtr,MatcherOrdering>  match(matchers().begin(), matchers().end());
 
   os << objects().size();
   for ( ObjectMap::iterator it = objects().begin();
 	it != objects().end(); ++it ) os << it->second;
   os << defaultParticles() << part << match << generators()
      << directories() << directoryStack() << globalLibraries() << readDirs();
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "(" << objects().size() << " objects in " << directories().size()
 	   << " directories) done" << endl;
 }
 
 string Repository::load(string filename) {
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "loading '" << filename << "'... " << flush;
   currentFileName() = filename;
   PersistentIStream * is = new PersistentIStream(filename);
   if ( !*is ) {
     delete is;
     // macro is passed in from -D compile flag
     string fullpath = string(THEPEG_PKGLIBDIR) + '/' + filename;
     is = new PersistentIStream(fullpath);
     if ( !*is ) {
       delete is;
       return "Error: Could not find repository '" + filename + "'.";
     }
   }
   *is >> allObjects() >> defaultParticles()
       >> particles() >> matchers() >> generators()
       >> directories() >> directoryStack() >> globalLibraries() >> readDirs();
   delete is;
   objects().clear();
   for ( ObjectSet::iterator it = allObjects().begin();
 	it != allObjects().end(); ++it )
     objects()[(**it).fullName()] = *it;
 
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "(" << objects().size() << " objects in " << directories().size()
 	   << " directories) done\nUpdating... " << flush;
   BaseRepository::resetAll(allObjects());
   BaseRepository::update();
 
   if ( ThePEG_DEBUG_ITEM(3) )
     clog() << "done" << endl;
   return "";
 }
 
 void Repository::stats(ostream & os) {
   os << "number of objects:        " << setw(6) << objects().size() << endl;
   os << "number of objects (all):  " << setw(6) << allObjects().size() << endl;
   os << "number of particles:        " << setw(6) << particles().size() << endl;
   os << "number of matchers:         " << setw(6) << matchers().size() << endl;
 }
 
 string Repository::read(string filename, ostream & os) {
   ifstream is;
   string file = filename;
   if ( file[0] == '/' ) {
     if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= trying to open " << file << " =)" << endl;
     is.open(file.c_str());
   }
   else {
     vector<string> dirs(readDirs().rbegin(), readDirs().rend());
     dirs.push_back(currentReadDirStack().top()); 
     if ( ThePEG_DEBUG_LEVEL > 1 ) {
       os << "(= search path order =)\n(== ";
       std::copy(dirs.rbegin(), dirs.rend(), std::ostream_iterator<string>(os, " ==)\n(== "));
       os << ")" << endl;
     }
     while ( dirs.size() ) {
       string dir = dirs.back();
       if ( dir != "" && dir[dir.length() -1] != '/' ) dir += '/';
       file = dir + filename;
       is.clear();
       if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= trying to open " << file << " =)" << endl;
       is.open(file.c_str());
       if ( is ) break;
       if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= no, try next search path =)" << endl;
       dirs.pop_back();
     }
   }
   if ( !is ) {
     return "Error: Could not find input file '" + filename + "'";
   }
   if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= yes =)" << endl;
   const string dir = StringUtils::dirname(file);
   if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= pushing <" << dir << "> to stack =)" << endl;
   currentReadDirStack().push(dir);
   try {
     Repository::read(is, os);
     if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= popping <" << currentReadDirStack().top() << "> from stack =)" << endl;
     currentReadDirStack().pop();
   }
   catch ( ... ) {
     if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= popping <" << currentReadDirStack().top() << "> from stack =)" << endl;
     currentReadDirStack().pop();
     throw;
   }
   return "";
 }
 
 string Repository::
 modifyEventGenerator(EventGenerator & eg, string filename, 
 		     ostream & os, bool initOnly) {
   ObjectSet objs = eg.objects();
   objs.insert(&eg);
   for ( ObjectSet::iterator it = objs.begin(); it != objs.end(); ++it ) {
     string name = (**it).fullName();
     if ( name.rfind('/') != string::npos )
       CreateDirectory(name.substr(0, name.rfind('/') + 1));
     objects()[name] = *it;
     allObjects().insert(*it);
   }
   
   string msg = read(filename, os);
 
   if ( !msg.empty() )
     return msg;
  
   for_each(objs, mem_fun(&InterfacedBase::reset)); 
   eg.initialize(initOnly);
 
   if ( !generators().empty() )
     msg += "Warning: new generators were initialized while modifying "
       + eg.fullName() + ".\n";
 
   return msg;
 }
 
 void Repository::resetEventGenerator(EventGenerator & eg) {
 
   ObjectSet objs = eg.objects();
   objs.insert(&eg);
   for ( ObjectSet::iterator it = objs.begin(); it != objs.end(); ++it ) {
     string name = (**it).fullName();
     if ( name.rfind('/') != string::npos )
       CreateDirectory(name.substr(0, name.rfind('/') + 1));
     objects()[name] = *it;
     allObjects().insert(*it);
   }
   
   for_each(objs, mem_fun(&InterfacedBase::reset)); 
   eg.initialize(true);
 
 }
 
 void Repository::execAndCheckReply(string line, ostream & os) {
   string reply = exec(line, os);
   if ( reply.size() ) 
     os << reply;
   if ( reply.size() && reply[reply.size()-1] != '\n' ) 
     os << endl;
   if ( exitOnError() && reply.size() >= 7 
        && reply.substr(0, 7) == "Error: " )
     exit(exitOnError());
 }
 
 void Repository::read(istream & is, ostream & os, string prompt) {
 #ifdef HAVE_LIBREADLINE
   if ( &is == &std::cin ) {
     char * line_read = 0;
     do {
       if ( line_read ) {
 	free(line_read);
 	line_read = 0;
       }
       
       line_read = readline(prompt.c_str());
       
       if ( line_read && *line_read ) {
 	string line = line_read;
 	while ( !line.empty() && line[line.size() - 1] == '\\' ) {
 	  line[line.size() - 1] = ' ';
 	  char * cont_read = readline("... ");
 	  if ( cont_read ) {
 	    line += cont_read;
 	    free(cont_read);
 	  }
 	}
 	if ( prompt.empty() && ThePEG_DEBUG_LEVEL > 0 )
 	  os << "(" << line << ")" << endl;
 #ifdef HAVE_READLINE_HISTORY
 	add_history(line.c_str());
 #endif // HAVE_READLINE_HISTORY
 	execAndCheckReply(line, os);
       }
     }
     while ( line_read );
   }
   else {
 #endif // HAVE_LIBREADLINE
     string line;
     if ( prompt.size() ) os << prompt;
     while ( getline(is, line) ) {
       while ( !line.empty() && line[line.size() - 1] == '\\' ) {
 	line[line.size() - 1] = ' ';
 	string cont;
 	if ( prompt.size() ) os << "... ";
 	getline(is, cont);
 	line += cont;
       }
       if ( prompt.empty() && ThePEG_DEBUG_LEVEL > 0 )
 	os << "(" << line << ")" << endl;
       execAndCheckReply(line, os);
       if ( prompt.size() ) os << prompt;
     }
 #ifdef HAVE_LIBREADLINE
   }
 #endif
   if ( prompt.size() ) os << endl;
 }
 
 
 string Repository::copyParticle(tPDPtr p, string newname) {
   DirectoryAppend(newname);
   
   string newdir = newname.substr(0, newname.rfind('/')+1);
   newname =newname.substr(newname.rfind('/')+1);
   if ( newname.empty() ) newname = p->name();
   if ( GetPointer(newdir + newname) )
     return "Error: Cannot create particle " + newdir + newname +
       ". Object already exists.";
   if ( p->CC() && GetPointer(newdir + p->CC()->name()) )
     return "Error: Cannot create anti-particle " + newdir + newname +
       ". Object already exists.";
   PDPtr pd = p->pdclone();
   Register(pd, newdir + newname);
   pd->theDecaySelector.clear();
   pd->theDecayModes.clear();
   pd->isStable = true;
   if ( p->CC() ) {
     PDPtr apd = p->CC()->pdclone();
     Register(apd, newdir + apd->name());
     apd->theDecaySelector.clear();
     apd->theDecayModes.clear();
     apd->isStable = true;
     pd->theAntiPartner = apd;
     apd->theAntiPartner = pd;
     pd->syncAnti = p->syncAnti;
     apd->syncAnti = p->CC()->syncAnti;
   }
   HoldFlag<> dosync(pd->syncAnti, true);
   for ( DecaySet::const_iterator it = p->theDecayModes.begin();
 	it != p->theDecayModes.end(); ++it )
     pd->addDecayMode(*it);
   return "";
 }
 
 void Repository::remove(tIBPtr ip) {
   ObjectMap::iterator it = objects().find(ip->fullName());
   if ( it == objects().end() || ip != it->second ) return;
   objects().erase(it);
   allObjects().erase(ip);
   if ( dynamic_ptr_cast<tPDPtr>(ip) ) {
     particles().erase(dynamic_ptr_cast<tPDPtr>(ip));
     defaultParticles().erase(dynamic_ptr_cast<tPDPtr>(ip)->id());
   }
   if ( dynamic_ptr_cast<tPMPtr>(ip) )
     matchers().erase(dynamic_ptr_cast<tPMPtr>(ip));
 }
 
 string Repository::remove(const ObjectSet & rmset) {
   ObjectSet refset;
   for ( ObjectMap::const_iterator i = objects().begin();
 	i != objects().end(); ++i ) {
     if ( member(rmset, i->second) ) continue;
     IVector ov = DirectReferences(i->second);
     for ( int j = 0, M = ov.size(); j < M; ++j )
       if ( member(rmset, ov[j]) ) {
 	refset.insert(i->second);
 	break;
       }
   }
   if ( refset.empty() ) {
     for ( ObjectSet::iterator oi = rmset.begin(); oi != rmset.end(); ++oi )
       remove(*oi);
     return "";
   }
   string ret = "Error: cannot remove the objects because the following "
     "objects refers to some of them:\n";
   for ( ObjectSet::iterator oi = refset.begin(); oi != refset.end(); ++oi )
     ret += (**oi).fullName() + "\n";
   return ret;
 }
    
 string Repository::exec(string command, ostream & os) {
   string cpcmd = command;
   try {
     string verb = StringUtils::car(command);
     command = StringUtils::cdr(command);
     if ( verb == "help" ) {
       help(command, os);
       return "";
     }
     if ( verb == "rm" ) {
       ObjectSet rmset;
       while ( !command.empty() ) {
 	string name = StringUtils::car(command);
 	DirectoryAppend(name);
 	IBPtr obj = GetPointer(name);
 	if ( !obj ) return "Error: Could not find object named " + name;
 	rmset.insert(obj);
 	command = StringUtils::cdr(command);
       }
       return remove(rmset);
     }
     if ( verb == "rmdir" || verb == "rrmdir" ) {
       string dir = StringUtils::car(command);
       DirectoryAppend(dir);
       if ( dir[dir.size() - 1] != '/' ) dir += '/';
       if ( !member(directories(), dir) )
 	return verb == "rmdir"? "Error: No such directory.": "";
       IVector ov = SearchDirectory(dir);
       if ( ov.size() && verb == "rmdir" )
 	return "Error: Cannot remove a non-empty directory. "
 	  "(Use rrmdir do remove all object and subdirectories.)";
       ObjectSet rmset(ov.begin(), ov.end());
       string ret = remove(rmset);
       if ( !ret.empty() ) return ret;
       StringVector dirs(directories().begin(), directories().end());
       for ( int i = 0, N = dirs.size(); i < N; ++ i )
 	if ( dirs[i].substr(0, dir.size()) == dir )
 	  directories().erase(dirs[i]);
       for ( int i = 0, N = directoryStack().size(); i < N; ++i )
 	if ( directoryStack()[i].substr(0, dir.size()) == dir )
 	  directoryStack()[i] = '/';
       return "";
     }
     if ( verb == "cp" ) {
       string name = StringUtils::car(command);
       DirectoryAppend(name);
       tPDPtr p = dynamic_ptr_cast<tPDPtr>(GetPointer(name));
       if ( p ) return copyParticle(p, StringUtils::cdr(command));
       return BaseRepository::exec(cpcmd, os);
     }
     if ( verb == "setup" ) {
       string name = StringUtils::car(command);
       DirectoryAppend(name);
       IBPtr obj = GetPointer(name);
       if ( !obj ) return "Error: Could not find object named " + name;
       istringstream is(StringUtils::cdr(command));
       readSetup(obj, is);
       // A particle may have been registered before but under the wrong id().
       PDPtr pd = dynamic_ptr_cast<PDPtr>(obj);
       if(pd) {
 	registerParticle(pd);
 	checkDuplicatePDGName(pd);
       }
       return "";
     }
     if ( verb == "decaymode" ) {
       string tag = StringUtils::car(command);
       DMPtr dm = DecayMode::constructDecayMode(tag);
       if ( !dm ) return "Error: Could not create decay mode from the tag " +
 		   StringUtils::car(command);
       istringstream is(StringUtils::cdr(command));
       readSetup(dm, is);
       if ( !dm->CC() ) return "";
 
       if ( dm->CC()->parent()->synchronized() ) {
 	dm->CC()->synchronize();
 	return "";
       }
 
       if ( !dm->CC()->decayer() )
 	return FindInterface(dm, "Decayer")->
 	  exec(*dm->CC(), "set", dm->decayer()->fullName());
       return "";
     }
     if ( verb == "makeanti" ) {
       string name = StringUtils::car(command);
       DirectoryAppend(name);
       tPDPtr p = dynamic_ptr_cast<tPDPtr>(GetPointer(name));
       if ( !p ) return "Error: No particle named " + name;
       name = StringUtils::car(StringUtils::cdr(command));
       DirectoryAppend(name);
       tPDPtr ap = dynamic_ptr_cast<tPDPtr>(GetPointer(name));
       if ( !ap ) return "Error: No particle named " + name;
       ParticleData::antiSetup(PDPair(p, ap));
       return "";
     }
     if ( verb == "read" ) {
       // remember directory we're in
       string cwd = directoryStack().back();
       string filename = StringUtils::car(command);
       string msg = read(filename, os);
       // Return to the original directory, so that
       // calling 'read' in an input file will not change the 
       // repository directory you're in
       ChangeDirectory(cwd);
       return msg;
     }
     if ( verb == "load" ) {
       return load(StringUtils::car(command));
     }      
     if ( verb == "save" ) {
       save(StringUtils::car(command));
       return "";
     }
     if ( verb == "lsruns" ) {
       string ret;
       for ( GeneratorMap::iterator ieg = generators().begin();
 	    ieg != generators().end(); ++ieg ) ret += ieg->first + "\n";
       return ret;
     }
     if ( verb == "makerun" ) {
       string runname = StringUtils::car(command);
       string generator = StringUtils::car(StringUtils::cdr(command));
       DirectoryAppend(generator);
       EGPtr eg = BaseRepository::GetObject<EGPtr>(generator);
       makeRun(eg, runname);
       return "";
     }
     if ( verb == "rmrun" ) {
       string runname = StringUtils::car(command);
       generators().erase(runname);
       return "";
     }
     if ( verb == "saverun" || verb == "saverunfile" || verb == "run" ) {
       string runname = StringUtils::car(command);
       string generator = StringUtils::car(StringUtils::cdr(command));
       DirectoryAppend(generator);
       GeneratorMap::iterator ieg = generators().find(runname);
       EGPtr eg;
       if ( ieg == generators().end() ) {
 	eg = BaseRepository::GetObject<EGPtr>(generator);
 	eg = makeRun(eg, runname);
       } else
 	eg = ieg->second;
       if ( !eg )
 	return "Error: Could not create/find run named'" + runname + "'.";
       if ( verb == "run" ) 
 	eg->go();
       else if ( verb == "saverunfile" ) {
 	string file = generator;
 	PersistentOStream os(file, globalLibraries());
 	os << eg;
 	if ( !os ) return "Save failed! (I/O error)";
       } else {
 	string file = eg->filename() + ".run";
 	PersistentOStream os(file, globalLibraries());
 	os << eg;
 	if ( !os ) return "Save failed! (I/O error)";
       }
       return "";
     }
     if ( verb == "removerun" ) {
       string runname = StringUtils::car(command);
       GeneratorMap::iterator ieg = generators().find(runname);
       if ( ieg != generators().end() ) {
 	generators().erase(ieg);
 	return "";
       } else
 	return "Error: No run named '" + runname + "' available.";
     }
     if ( verb == "create" ) {
       string className = StringUtils::car(command);
       command = StringUtils::cdr(command);
       string name = StringUtils::car(command);
       const ClassDescriptionBase * db = DescriptionList::find(className);
       command = StringUtils::cdr(command);
       while ( !db && command.length() ) {
 	string library = StringUtils::car(command);
 	command = StringUtils::cdr(command);
 	DynamicLoader::load(library);
 	db = DescriptionList::find(className);
       }
       if ( !db ) {
 	string msg = "Error: " + className + ": No such class found.";
 	if ( !DynamicLoader::lastErrorMessage.empty() )
 	  msg += "\nerror message from dynamic loader:\n" +
 	    DynamicLoader::lastErrorMessage;
 	return msg;
       }
       IBPtr obj = dynamic_ptr_cast<IBPtr>(db->create());
       if ( !obj ) return "Error: Could not create object of this class class.";
       if ( name.empty() ) return "Error: No name specified.";
       Register(obj, name);
       return "";
     }
     if ( verb == "defaultparticle" ) {
       while ( !command.empty() ) {
 	string name = StringUtils::car(command);
 	DirectoryAppend(name);
 	tPDPtr p = dynamic_ptr_cast<tPDPtr>(GetPointer(name));
 	if ( !p ) return "Error: No particle named " + name;
 	defaultParticle(p);
 	command = StringUtils::cdr(command);
       }
       return "";
     }
     if ( verb == "EXITONERROR" ) {
       exitOnError() = 1;
       return "";
     }
   }
   catch (const Exception & e) {
     e.handle();
     return "Error: " + e.message();
   }
 
   return BaseRepository::exec(cpcmd, os);
 }
 
 void Repository::help(string cmd, ostream & os) {
  
   cmd = StringUtils::car(cmd);
 
   if ( cmd == "cd" )
     os << "Usage: cd <directory>" << endl
        << "Set the current directory to <directory>." << endl;
   else if ( cmd == "mkdir" )
     os << "Usage: mkdir <path-name>" << endl
        << "Create a new directory called with the given path name." << endl;
   else if ( cmd == "rmdir" )
     os << "Usage: rmdir <directory>" << endl
        << "Remove an empty directory." << endl;
   else if ( cmd == "rrmdir" )
     os << "Usage: rrmdir <directory>" << endl
        << "Remove a directory and everything that is in it recursively." << endl
        << "Will only succeed if no other objects refers to the ones to "
        << "be deleted." << endl;
   else if ( cmd == "cp" )
     os << "Usage: cp <object> <path-name>" << endl
        << "Copy the given object to a new object with the given name." << endl;
   else if ( cmd == "setup" )
     os << "Usage: setup <object> <arguments> ..." << endl
        << "Tell a given object to read information given by the arguments."
        << endl;
   else if ( cmd == "decaymode" )
     os << "Usage: decaymode <tag> <branching fraction> <on|off> <decayer-object>"
        << endl
        << "Construct a decay mode from the given decay tag. The resulting "
        << "object will be inserted in the directory with the same path as "
        << "the decaying particle object. The given brancing fraction will "
        << "be set as well as the given decayer object. If the mode should "
        << "be switched on by default 1(on) should be specified (otherwise "
        << "0(off))." << endl;
   else if ( cmd == "makeanti" )
     os << "Usage: makeanti <particle-object> <particle-object>" << endl
        << "Indicate that the two given particle objects are eachothers "
        << "anti-partnets." << endl;
   else if ( cmd == "read" )
     os << "Usage: read <file-name>" << endl
        << "Read more commands from the given file. The file name can be "
        << "given relative to the current directory in the shell, or "
        << "relative to standard directories, or as an absolute path." << endl;
   else if ( cmd == "load" )
     os << "Usage: load <repository-file-name>" << endl
        << "Discard everything in the reopsitory and read in a completely "
        << "new repository from the given file." << endl;
   else if ( cmd == "save" )
     os << "Usage: save <file-name>" << endl
        << "Save the complete repository to the given file." << endl;
   else if ( cmd == "lsruns" )
     os << "Usage: lsruns" << endl
        << "List the run names of all initialized event generators." << endl;
   else if ( cmd == "makerun" )
     os << "Usage: makerun <run-name> <event-generator-object>" << endl
        << "Initialize the given event generator and assign a run name." << endl;
   else if ( cmd == "rmrun" )
     os << "Usage: rmrun <run-name>" << endl
        << "Remove the initialized event generator given by the run name."
        << endl;
   else if ( cmd == "saverun" )
     os << "Usage: saverun <run-name> <event-generator-object>" << endl
        << "Initialize the given event generator and assign a run name "
        << "and save it to a file named <run-name>.run" << endl;
   else if ( cmd == "run" )
     os << "Usage: run <run-name>" << endl
        << "Run the initialized event generator given b the run name." << endl;
   else if ( cmd == "create" )
     os << "Usage: create <class-name> <name> {<dynamic-library>}" << endl
        << "Create an object of the given class and assign the given name. "
        << "Optionally supply a dynamically loaded library where the class "
        << "is included." << endl;
   else if ( cmd == "pushd" )
     os << "Usage: pushd <directory>" << endl
        << "Set the current directory to <directory>, but keep the previous "
        << "working directory on the directory stack." << endl;
   else if ( cmd == "popd" )
     os << "Usage: popd" << endl
        << "Leave the current working directory and set the current "
        << "directory to the previous one on the directory stack." << endl;
   else if ( cmd == "pwd" )
     os << "Usage: pwd" << endl
        << "Print the current working directory." << endl;
   else if ( cmd == "dirs" )
     os << "Usage: dirs" << endl
        << " Print the contents of the directory stack." << endl;
   else if ( cmd == "mv" )
     os << "Usage: mv  <object> <path-name>" << endl
        << "Rename the given object to a new path name." << endl;
   else if ( cmd == "ls" )
     os << "Usage: ls {<directory>}" << endl
        << "List the objects and subdirectories in the current or given "
        << "directory." << endl;
   else if ( cmd == "library" )
     os << "Usage: library <dynamic-library>" << endl
        << "Make new classes available to the repository by dynamically "
        << "linking the given library." << endl;
   else if ( cmd == "globallibrary" )
     os << "Usage: globallibrary <dynamic-library>" << endl
        << "Make new classes available to the repository by dynamically "
        << "linking the given library. If this repository is saved and read "
        << "in again, this library will be linked in from the beginning." << endl;
   else if ( cmd == "rmgloballibrary" )
     os << "Usage: rmgloballibrary <dynamic-library>" << endl
        << "Remove a dynamic library previously added with globallibrary."
        << endl;
   else if ( cmd == "appendpath" )
     os << "Usage: appendpath <unix-directory>" << endl
        << "Add a search path for dynamic libraries to the end of the "
        << "search list." << endl;
   else if ( cmd == "lspaths" )
     os << "Usage: lspaths" << endl
        << "List search paths for dynamic libraries." << endl;
   else if ( cmd == "prependpath" )
     os << "Usage: prependpath <unix-directory>" << endl
        << "Add a search path for dynamic libraries to the beginning of the "
        << "search list." << endl;
   else if ( cmd == "doxygendump" )
     os << "Usage: doxygendump <namespace> <filename>" << endl
        << "Extract doxygen documentation of all loaded classes in the "
        << "given name space and weite it to a file.." << endl;
   else if ( cmd == "mset" || cmd == "minsert" || cmd == "mdo" )
     os << "Usage: " << cmd << " <directory> <class> <interface> <value>" << endl
        << "Recursively find in the given directory all objects of the "
        << "given class and call '" << cmd.substr(1)
        << "' with the given value for the given interface." << endl;
   else if ( cmd == "msetdef" || cmd == "mget" || cmd == "mdef" ||
 	    cmd == "mmin" || cmd == "mmax" || cmd == "merase" )
     os << "Usage: " << cmd << " <directory> <class> <interface>" << endl
        << "Recursively find in the given directory all objects of the given "
        << "class and call '" << cmd.substr(1)
        << "' for the given interface." << endl;
   else if ( cmd == "set" )
     os << "Usage: set <object>:<interface> <value>" << endl
        << "Set the interface for the given object to the given value." << endl;
   else if ( cmd == "setdef" )
     os << "Usage: setdef <object>:<interface>" << endl
        << "Set the interface for the given object to its default value." << endl;
   else if ( cmd == "insert" )
     os << "Usage: insert <object>:<interface> <value>" << endl
        << "Insert a value in the vector interface of the given object." << endl;
   else if ( cmd == "erase" )
     os << "Usage: erase <object>:<interface>" << endl
        << "Erase a value from the vector interface of the given object." << endl;
   else if ( cmd == "do" )
     os << "Usage: do <object>:<command-interface> <arguments>" << endl
        << "Call the command interface of the given object with the "
        << "given arguments." << endl;
   else if ( cmd == "get" )
     os << "Usage: get <object>:<interface>" << endl
        << "Print the value of the interface of the given object." << endl;
   else if ( cmd == "def" )
     os << "Usage: def <object>:<interface>" << endl
        << "Print the default value of the interface of the given object."
        << endl;
   else if ( cmd == "min" )
     os << "Usage: min <object>:<interface>" << endl
        << "Print the minimum value of the interface of the given object."
        << endl;
   else if ( cmd == "max" )
     os << "Usage: max <object>:<interface>" << endl
        << "Print the maximum value of the interface of the given object."
        << endl;
   else if ( cmd == "describe" )
     os << "Usage: describe <object>{:<interface>}" << endl
        << "Describe the given object or an interface of the object." << endl;
   else if ( cmd == "lsclass" )
     os << "Usage: lsclass" << endl
        << "List all classes available in the repository." << endl;
   else if ( cmd == "all" ) {
     os << "Available commands:"
        << endl
        << "* cd, mkdir, rmdir, rrmdir, pwd, cp, mv, rm, pushd, popd, dirs, ls:\n"
        << "  Manipulate the repository structure. Analogous to unix "
        << "shell commands."
        << endl
        << "* create, setup, decaymode makeanti:\n"
        << "  Create or setup an object."
        << endl
        << "* set, get, insert, erase, do, detdef, def, min, max, describe\n"
        << "  mset, minsert, mdo, msetdef, mdef, mmin, mmax, merase:\n"
        << "  Manipulate interfaces to objects."
        << endl
        << "* makerun, saverun, run, lsruns, rmrun:\n"
        << "  Create and handle initialized event genrators which can be run."
        << endl
        << "* read, load, library globallibrary, rmgloballibrary,\n"
        << "  appendpath, prependpath, lspaths, doxygendump:\n"
        << "  Handle files external files and libraries."
        << endl;
     os << "Do 'help syntax' for help on syntax." << endl
        << "Do 'help <command>' for help on a particular command." << endl;
   }
   else if ( cmd == "syntax" )
     os << "* <directory> = '/' | <name> | <directory>/<name>" << endl
        << "  <object> = <name> | <directory>/<name> | <object>:<ref-interface>\n"
        << "  Analogous to a unix file structure, an object can be "
        << "specified with an\n  absolute path or a path relative to "
        << "the current directory." << endl
        << "* <interface> = <interface-name>|<interface-name>[<index>]" << endl
        << "  An interface can be a parameter (floating point, integer or "
        << "string),\n  a switch (integer, possibly named), a reference to "
        << "another object in the\n  repository or a command which takes "
        << "an arbitrary string as argument.\n  There are also vector interfaces "
        << "of parameters and references for which\n  an index must be supplied."
        << endl;
   else {
     if ( !cmd.empty() ) os << "No command '" << cmd << "' found." << endl;
     os << "Common commands:" << endl
        << "* cd, mkdir, rmdir, pwd, cp, mv, rm:\n"
        << "  Manipulate the repository structure. Analogous to unix "
        << "shell commands." << endl
        << "* create, setup:\n"
        << " Create an object." << endl
        << "set, get, insert, erase, do:\n"
        << " Manipulate interfaces to objects." << endl
        << "* makerun, saverun, run, lsruns:\n"
        << " Create and handle initialized event genrators which can be run."
        << endl;
     os << "Do 'help all' for a complete list of commands." << endl
        << "Do 'help syntax' for help on syntax." << endl
        << "Do 'help <command>' for help on a particular command." << endl;
   }
 
 }
 
 Repository::Repository() {
   ++ninstances;
 }
 
 Repository::~Repository() {
   --ninstances;
   if ( ninstances <= 0 ) {
     generators().clear();
   }
 }
 
 void Repository::checkDuplicatePDGName(PDPtr pd) {
   string name = pd->PDGName();
   for ( ParticleMap::iterator pit = defaultParticles().begin();
 	pit != defaultParticles().end(); ++pit ) {
     if( pit->second == pd) continue;
     if ( pit->second->PDGName() == name ) {
       std::cerr << "Using duplicate PDGName " << pd->PDGName()
 		<< " for a new particle.\n This can cause problems and is not "
 		<< "recommended.\n If this second particle is a new particle "
 		<< "in a BSM Model we recommend you change the name of the particle.\n";
     }
   }
   for ( ParticleDataSet::iterator pit = particles().begin();
 	pit != particles().end(); ++pit ) {
     if( *pit == pd) continue;
     if ( (**pit).PDGName() == name ) {
       std::cerr << "Using duplicate PDGName " << pd->PDGName()
 		<< " for a new particle.\n This can cause problems and is not "
 		<< "recommended.\n If this second particle is a new particle "
 		<< "in a BSM Model we recommend you change the name of the particle.\n";
     }
   }
 }
 
 int Repository::ninstances = 0;
 
 namespace {
 static string version_ =
 #include "versionstamp.inc"
 "";
 }
 
 string Repository::version() {
   return ::version_;
 }
 
 string Repository::banner() {
   const auto now    = std::chrono::system_clock::now();
   const auto now_c  = std::chrono::system_clock::to_time_t(now);
   string time = ">>>> " ;
   time += StringUtils::stripws(string(std::ctime(&now_c))) + ' ';
   time += string(max(0,74 - int(time.size())), ' ');
   time += "<<<<";
 
   string line = ">>>> Toolkit for HEP Event Generation - "
     + Repository::version() + ' ';
   line += string(max(0,78 - int(line.size())), '<');
 
   string block = string(78, '>') + '\n' 
                  + line + '\n' 
                  + time + '\n'
                  + string(78, '<') + '\n';
   return block;
 }