diff --git a/MatrixElement/ColourLines.cc b/MatrixElement/ColourLines.cc
--- a/MatrixElement/ColourLines.cc
+++ b/MatrixElement/ColourLines.cc
@@ -1,143 +1,142 @@
 // -*- C++ -*-
 //
 // ColourLines.cc is a part of ThePEG - Toolkit for HEP Event Generation
 // Copyright (C) 1999-2011 Leif Lonnblad
 //
 // ThePEG is licenced under version 2 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 ColourLines class.
 //
 
 #include "ColourLines.h"
 #include "ColourLines.xh"
 #include "ThePEG/EventRecord/ColourLine.h"
 #include "ThePEG/EventRecord/MultiColour.h"
 #include "ThePEG/EventRecord/Particle.h"
 #include "ThePEG/Utilities/StringUtils.h"
 
 using namespace ThePEG;
 
 ColourLines::ColourLines(string s) {
   reset(s); 
 } 
 
 void ColourLines::reset(string s) { 
   theLines.clear(); 
   while ( true ) {
     string line = StringUtils::car(s, ",");
     line = StringUtils::stripws(line);
     Line l;
     while (line!="") {
       string loc = StringUtils::car(line);
       string first  = StringUtils::car(loc,":");
       string second = StringUtils::cdr(loc,":");
       if(second!="") {
 	int i;
 	istringstream is(first);
 	is >> i;
 	int j;
 	istringstream is2(second);
 	is2 >> j;
 	l.push_back(make_pair(i,j));
       }
       else {
 	int i;
 	istringstream is(first);
 	is >> i;
 	l.push_back(make_pair(i,0));
       }
       line = StringUtils::cdr(line);
     };
     if ( l.empty() ) return;
     theLines.push_back(l);
     s = StringUtils::cdr(s, ",");
   }
 }
 
 void ColourLines::connect(const tPVector & partons) const {
   VertexVector sinks;
   VertexVector sources;
   long np = partons.size();
 
   // Create each line and connect the specified partons to them. Save
   // all lines coming from a source or ending in a sink.
   for ( LineVector::size_type il = 0; il < theLines.size(); ++il ) {
     const Line & line = theLines[il];
     ColinePtr cline = new_ptr(ColourLine());
     for ( Line::size_type i = 0; i < line.size(); ++i ) {
       if ( line[i].first > np ) {
 	// this is a colour source.
 	int is = line[i].first - np;
 	sources.resize(is);
 	sources[is - 1].push_back(cline);
       } else if ( -line[i].first > np ) {
 	// this is a colour sink.
 	int is = -line[i].first - np;
-	sources.resize(is);
-	sources[is - 1].push_back(cline);
+	sinks.resize(is);
+	sinks[is - 1].push_back(cline);
       } else if ( line[i].first > 0 ) {
 	// This is a coloured particle.
 	if ( !partons[line[i].first - 1]->hasColour() )
 	  throw ColourGeometryException(partons, line);
 	if(line[i].second==0) {
 	  cline->addColoured(partons[line[i].first - 1]);
 	}
 	else {
 	  Ptr<MultiColour>::pointer colour = 
 	    dynamic_ptr_cast<Ptr<MultiColour>::pointer>(partons[line[i].first - 1]->colourInfo());
 	  assert(colour);
 	  colour->colourLine(cline,line[i].second);
 	}
       } else {
 	if ( !partons[-line[i].first - 1]->hasAntiColour() )
 	  throw ColourGeometryException(partons, line);
 	if(line[i].second==0) {
 	  cline->addAntiColoured(partons[-line[i].first - 1]);
 	}
 	else {
 	  Ptr<MultiColour>::pointer colour = 
 	    dynamic_ptr_cast<Ptr<MultiColour>::pointer>(partons[-line[i].first - 1]->colourInfo());
 	  assert(colour);
 	  colour->antiColourLine(cline,line[i].second);
 	}
       }
     }
   }
-
   // Now connect up all lines steming from sources.
   for ( VertexVector::size_type i = 0; i < sources.size(); ++i ) {
     if ( sources[i].empty() ) continue;
     if ( sources[i].size() != 3 ) throw ColourGeometryException(partons,
                                                                 vector<pair<int,int> >() );
     sources[i][0]->setSourceNeighbours(sources[i][1], sources[i][2]);
   }
 
   // Now connect up all lines ending in sinks.
   for ( VertexVector::size_type i = 0; i < sinks.size(); ++i ) {
     if ( sinks[i].empty() ) continue;
     if ( sinks[i].size() != 3 ) throw ColourGeometryException(partons,
                                                               vector<pair<int,int> >());
     sinks[i][0]->setSinkNeighbours(sinks[i][1], sinks[i][2]);
   }
 }
 
 ColourGeometryException::
 ColourGeometryException(const tPVector & p, const vector<pair<int,int> > & c) {
   if ( c.empty() )
     theMessage << "The number of colour lines steming from one colour source "
 	       << "or ending in one colour sink was not equal to three.\n";
   else {
     theMessage << "Cannot connect the following partons:\n";
     for ( unsigned i = 0; i < p.size(); ++i )
       theMessage << " " << p[i]->PDGName();
     theMessage << "\n to the following colour line:\n";
     for ( unsigned i = 0; i < c.size(); ++i ) theMessage << " (" << c[i].first << ","
 							 << c[i].second << ") ";
     theMessage << endl;
   }
   severity(maybeabort);
 }
 
diff --git a/MatrixElement/Tree2toNDiagram.cc b/MatrixElement/Tree2toNDiagram.cc
--- a/MatrixElement/Tree2toNDiagram.cc
+++ b/MatrixElement/Tree2toNDiagram.cc
@@ -1,462 +1,454 @@
 // -*- C++ -*-
 //
 // Tree2toNDiagram.cc is a part of ThePEG - Toolkit for HEP Event Generation
 // Copyright (C) 1999-2011 Leif Lonnblad
 //
 // ThePEG is licenced under version 2 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 Tree2toNDiagram class.
 //
 
 #include "Tree2toNDiagram.h"
 #include "ThePEG/EventRecord/SubProcess.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Utilities/Rebinder.h"
 #include "ThePEG/Utilities/UtilityBase.h"
 #include "ThePEG/Handlers/StandardXComb.h"
 #include "ThePEG/PDT/EnumParticles.h"
 
 using namespace ThePEG;
 
 Tree2toNDiagram::~Tree2toNDiagram() {}
 
 Tree2toNDiagram & Tree2toNDiagram::add(tcPDPtr pd) {
   if ( thePartons.size() < theNSpace ) addSpacelike(pd);
   else addTimelike(pd, nextOrig);
   return *this;
 }
 
 void Tree2toNDiagram::addTimelike(tcPDPtr pd, size_type orig) {
   if ( allPartons().size() < theNSpace ||
        orig >= allPartons().size())
     throw Tree2toNDiagramError();
   thePartons.push_back(pd);
   theParents.push_back(orig);
 }
 
 tPVector Tree2toNDiagram::
 construct(SubProPtr sp, const StandardXComb & xc, const ColourLines & cl) const {
   tPVector out;
   vector<Lorentz5Momentum> pout(xc.meMomenta().begin() + 2,
 				xc.meMomenta().end());
   if ( xc.needsReshuffling() )
     xc.reshuffle(pout);
   tPPair in = xc.lastPartons();
   if ( xc.mirror() ) swap(in.first, in.second);
 
   tPVector ret;
   if ( in.first->dataPtr() != allPartons()[0] ||
        in.second->dataPtr() != allPartons()[nSpace() - 1] )
     throw Tree2toNDiagramError();
 
   PVector slike;
   slike.push_back(in.first);
   for ( int i = 1; i < nSpace() - 1; ++i )
     slike.push_back(allPartons()[i]->produceParticle());
   slike.push_back(in.second);
   ret = tPVector(slike.begin(), slike.end());
   for ( size_type i = 1; i < slike.size() - 1; ++i ) {
     slike[i-1]->addChild(slike[i]);
     sp->addIntermediate(slike[xc.mirror()? i: slike.size() - 1 - i], false);
   }
   int io = pout.size();
   PVector tlike(allPartons().size() - nSpace());
   ParticleSet done;
   for ( int i = allPartons().size() - 1; i >=  nSpace(); --i ) {
     int it = i - nSpace();
     pair<int,int> ch = children(i);
     bool iso = ch.first < 0;
     if ( iso ) {
       tlike[it] = allPartons()[i]->produceParticle(pout[--io]);
       done.insert(tlike[it]);
     } else {
       Lorentz5Momentum p = tlike[ch.first - nSpace()]->momentum() +
 	tlike[ch.second - nSpace()]->momentum();
       tlike[it] = allPartons()[i]->produceParticle(p);
     }
     if ( parent(i) < nSpace() ) {
       slike[parent(i)]->addChild(tlike[it]);
       if ( parent(i) == nSpace() - 2 )
 	slike[parent(i) + 1]->addChild(tlike[it]);
     }
     if ( !iso ) {
       tlike[it]->addChild(tlike[ch.first - nSpace()]);
       tlike[it]->addChild(tlike[ch.second - nSpace()]);
     }
     if ( iso ) out.push_back(tlike[it]);
     else sp->addIntermediate(tlike[it], false);
   }
   ret.insert(ret.end(), tlike.begin(), tlike.end());
   for ( int i = 0, N = out.size(); i < N; ++i )
     sp->addOutgoing(out[xc.mirror()? i: out.size() - i - 1], false);
   for ( PVector::size_type i = 0; i < slike.size() - 2; ++i ) {
     pair<int,int> ch = children(i);
     slike[ch.first]->set5Momentum(slike[i]->momentum() -
 				  tlike[ch.second - nSpace()]->momentum());
   }
 
   cl.connect(ret);
 
   return out;
 }
 
 tcPDVector Tree2toNDiagram::outgoing() const {
   tcPDVector pdv;
   for ( size_type i = nSpace(); i < allPartons().size(); ++i )
     if ( children(i).first < 0 ) pdv.push_back(allPartons()[i]);
   return pdv;
 }
 
 tcPDVector Tree2toNDiagram::external() const {
   tcPDVector pdv;
   pdv.push_back(allPartons()[0]);
   pdv.push_back(allPartons()[nSpace() - 1]);
   for ( size_type i = nSpace(); i < allPartons().size(); ++i )
     if ( children(i).first < 0 ) pdv.push_back(allPartons()[i]);
   return pdv;
 }
 
 tcPDPair Tree2toNDiagram::incoming() const {
   return tcPDPair(allPartons()[0], allPartons()[nSpace() - 1]);
 }
 
 pair<int,int> Tree2toNDiagram::children(int ii) const {
   pair<int,int> ret = make_pair(-1, -1);
   for ( size_type i = 0; i < theParents.size(); ++i ) {
     if ( parent(i) == ii ) {
       if ( ret.first < 0 ) ret.first = i;
       else if ( ret.second < 0 ) ret.second = i;
       else throw Tree2toNDiagramError();
     }
   }
   return ret;
 }
 
 void Tree2toNDiagram::check() {
   vector< pair<int,int> > children(allPartons().size(), make_pair(-1, -1));
   theNOutgoing = 0;
   for ( size_type i = nSpace(); i < allPartons().size(); ++i ) {
     if ( children[parent(i)].first < 0 ) children[parent(i)].first = i;
     else if ( children[parent(i)].second < 0 ) children[parent(i)].second = i;
     else throw Tree2toNDiagramError();
   }
   for ( size_type i = nSpace(); i < allPartons().size(); ++i ) {
     if ( children[i].first < 0 && children[i].second < 0 ) ++theNOutgoing;
     else if ( children[i].first < 0 || children[i].second < 0 )
       throw Tree2toNDiagramError();
   }
   cPDVector parts(2);
   parts[0] = incoming().first;
   parts[1] = incoming().second;
   tcPDVector out(outgoing());
   parts.insert(parts.end(), out.begin(), out.end());
   partons(2, parts, nextOrig + 1);
 }
 
 bool Tree2toNDiagram::isSame (tcDiagPtr diag) const {
   Ptr<Tree2toNDiagram>::tcptr cmp = 
     dynamic_ptr_cast<Ptr<Tree2toNDiagram>::tcptr>( diag );
   if ( !cmp )
     return false;
   if ( cmp->nSpace() != nSpace() )
     return false;
   return equals(cmp) && external() == cmp->external();
 }
 
 bool Tree2toNDiagram::isSame (tcDiagPtr diag, map<int,int>& remap) const {
   Ptr<Tree2toNDiagram>::tcptr cmp = 
     dynamic_ptr_cast<Ptr<Tree2toNDiagram>::tcptr>( diag );
   if ( !cmp )
     return false;
   if ( cmp->nSpace() != nSpace() )
     return false;
   remap.clear();
   remap[0] = 0;
   return equals(cmp,remap);
 }
 
 bool Tree2toNDiagram::equals(Ptr<Tree2toNDiagram>::tcptr diag, 
 			     int start, int startCmp) const {
 
   // one leg ended externally while the other still has children left
   if ( start < 0 || startCmp < 0 )
     return false;
 
-  // match, if both are external legs
-  if ( start < 0 && startCmp < 0 )
-    return true;
-
   // no match, if the legs are not the same
   if ( allPartons()[start] != diag->allPartons()[startCmp] )
     return false;
 
   pair<int,int> ch = children(start);
   pair<int,int> chCmp = diag->children(startCmp);
 
   // start and startCmp are matching external legs
   if ( ch.first < 0 && chCmp.first < 0 ) {
     return true;
   }
 
   // check the first combination of children
   bool match =
     equals(diag,ch.first,chCmp.first) &&
     equals(diag,ch.second,chCmp.second);
 
   // also try swapped legs on same vertex for time-like legs
   if ( !match && start > nSpace() - 1 )
     match = 
       equals(diag,ch.first,chCmp.second) &&
       equals(diag,ch.second,chCmp.first);
 
   return match;
     
 }
 
 bool Tree2toNDiagram::equals(Ptr<Tree2toNDiagram>::tcptr diag, 
 			     map<int,int>& remap,
 			     int start, int startCmp) const {
 
   // one leg ended externally while the other still has children left
   if ( start < 0 || startCmp < 0 )
     return false;
 
-  // match, if both are external legs
-  if ( start < 0 && startCmp < 0 )
-    return true;
-
   // no match, if the legs are not the same
   if ( allPartons()[start] != diag->allPartons()[startCmp] )
     return false;
 
   pair<int,int> ch = children(start);
   pair<int,int> chCmp = diag->children(startCmp);
 
   // start and startCmp are matching external legs, which require remapping of
   // external labels
   if ( ch.first < 0 && chCmp.first < 0 ) {
     remap[externalId(start)] = diag->externalId(startCmp);
     return true;
   }
 
   // check the first combination of children
   bool match =
     equals(diag,remap,ch.first,chCmp.first) &&
     equals(diag,remap,ch.second,chCmp.second);
 
   // also try swapped legs on same vertex for time-like legs
   if ( !match && start > nSpace() - 1 )
     match = 
       equals(diag,remap,ch.first,chCmp.second) &&
       equals(diag,remap,ch.second,chCmp.first);
 
   return match;
 
 }
 
 int Tree2toNDiagram::externalId(int id) const {
   if ( id < 0 )
     return -1;
   if ( id == 0 )
     return 0;
   if ( id == nSpace() - 1 )
     return 1;
   int k = 1;
   for ( size_type i = nSpace(); i < allPartons().size(); ++i ) {
     if ( children(i).first < 0 ) ++k;
     if ( i == size_type(id) )
       break;
   }
   return k;
 }
 
 int Tree2toNDiagram::diagramId(int id) const {
   if ( id < 0 )
     return -1;
   if ( id == 0 ) return 0;
   if ( id == 1 ) return nSpace() - 1;
   int k = 1;
   size_type i = nSpace();
   for ( ; i < allPartons().size(); ++i ) {
     if ( children(i).first < 0 ) ++k;
     if ( k == id )
       break;
   }
   return i;
 }
 
 int Tree2toNDiagram::mergeEmission(int emitter, int id, map<int,int>& remap) {
 
   if ( id < 2 )
     return -1;
 
   if ( remap.find(emitter) != remap.end() ) {
     remap.erase(emitter);
   }
   if ( remap.find(id) != remap.end() ) {
     remap.erase(id);
   }
 
   for ( map<int,int>::iterator rm = remap.begin();
 	rm != remap.end(); ++rm ) {
     if ( rm->first == 0 || rm->first == 1 ) {
       rm->second = rm->first;
     } else {
       rm->second = diagramId(rm->first);
     }
   }
 
   // translate to diagram id
   int did = diagramId(id);
   int demitter = diagramId(emitter);
 
   if ( children(did) != make_pair(-1,-1) )
     return -1;
 
   // now get the parent
   int p = parent(did);
 
   int npos = -1;
   if ( p == 0 || p == nSpace() - 2 ) {
     npos = ( p == 0 ? 0 : 1 );
   } else if ( p >= nSpace() ) {
     if ( id > emitter )
       npos = emitter;
     else
       npos = emitter - 1;
   }
 
   pair<int,int> remove;
 
   size_type theNSpaceBackup = theNSpace;
   int theNOutgoingBackup = theNOutgoing;
   int nextOrigBackup = nextOrig;
   cPDVector thePartonsBackup = thePartons;
   vector<int> theParentsBackup = theParents;
 
   int deltaFlow = 0;
   if ( npos == 1 ) {
     if ( thePartons[did]->CC() )
       deltaFlow -= ( thePartons[did]->id() < 0 ? -1 : 1 );
     if ( thePartons[nSpace()-1]->CC() )
       deltaFlow += ( thePartons[nSpace()-1]->id() < 0 ? -1 : 1 );
   }
 
   // emitted from spacelike
   if ( p == 0 || p == nSpace() - 2 ) {
     if ( p == 0 && p != demitter )
       return -1;
     if ( p == nSpace() - 2 && demitter != nSpace()-1 )
       return -1;
     if ( p == 0 )
       remove = make_pair(p,did);
     else
       remove = make_pair(nSpace()-1,did);
     --theNSpace;
     --theNOutgoing;
   } else if ( p >= nSpace() ) {
     remove = children(p);
     if ( remove.first != demitter )
       swap(remove.first,remove.second);
     if ( remove != make_pair(demitter,did) )
       return -1;
     --theNOutgoing;
   } else {
     return -1;
   }
 
   if ( remove.first > remove.second )
     swap(remove.first,remove.second);
 
   for ( map<int,int>::iterator rm = remap.begin();
 	rm != remap.end(); ++rm ) {
     if ( rm->first > 1 ) {
       if ( rm->second > remove.first &&
 	   rm->second < remove.second )
 	rm->second -= 1;
       else if ( rm->second > remove.second )
 	rm->second -= 2;
     }
   }
 
   for ( unsigned int k = remove.first + 1; k < theParents.size(); ++k ) {
     if ( theParents[k] >= remove.first && 
 	 theParents[k] < remove.second &&
 	 theParents[k] >= 0 )
       theParents[k] -= 1;
     else if ( theParents[k] > remove.second && theParents[k] > 0 )
       theParents[k] -= 2;
   }
   thePartons.erase(thePartons.begin() + remove.second);
   theParents.erase(theParents.begin() + remove.second);
   thePartons.erase(thePartons.begin() + remove.first);
   theParents.erase(theParents.begin() + remove.first);
 
   if ( npos > 1 )
     if ( npos != externalId(p) ) {
       pair<int,int> swapDiagIds(p,diagramId(npos));
       swap(thePartons[swapDiagIds.first],thePartons[swapDiagIds.second]);
       swap(theParents[swapDiagIds.first],theParents[swapDiagIds.second]);
       for ( map<int,int>::iterator rm = remap.begin();
 	    rm != remap.end(); ++rm ) {
 	if ( rm->first > 1 ) {
 	  if ( rm->second == swapDiagIds.first ) {
 	    rm->second = swapDiagIds.second;
 	  } else if ( rm->second == swapDiagIds.second ) {
 	    rm->second = swapDiagIds.first;
 	  }
 	}
       }
     }
 
   for ( map<int,int>::iterator rm = remap.begin();
 	rm != remap.end(); ++rm ) {
     if ( rm->first > 1 ) {
       rm->second = externalId(rm->second);
     }
   }
 
   if ( npos == 1 ) {
     if ( thePartons[nSpace()-1]->CC() )
       deltaFlow -= ( thePartons[nSpace()-1]->id() < 0 ? -1 : 1 );
 
     if ( deltaFlow != 0 )
       thePartons[nSpace()-1] = thePartons[nSpace()-1]->CC();
 
   }
 
   try {
     check();
   } catch (Tree2toNDiagramError&) {
     theNSpace = theNSpaceBackup;
     theNOutgoing = theNOutgoingBackup;
     nextOrig = nextOrigBackup;
     thePartons = thePartonsBackup;
     theParents = theParentsBackup;
     return -1;
   }
 
   return npos;
 
 }
  
 ClassDescription<Tree2toNDiagram> Tree2toNDiagram::initTree2toNDiagram;
 
 void Tree2toNDiagram::persistentInput(PersistentIStream & is, int) {
   is >> theNSpace >> theNOutgoing >> thePartons >> theParents >> nextOrig;
 }
 
 void Tree2toNDiagram::persistentOutput(PersistentOStream & os) const {
   os << theNSpace << theNOutgoing << thePartons << theParents << nextOrig;
 }
 
 Tree2toNDiagramError::Tree2toNDiagramError() {
   theMessage << "An error occurred while setting up a diagram of class "
 	     << "'Tree2toNDiagram'.";
   severity(abortnow);
 }