diff --git a/MatrixElement/General/MEff2ff.cc b/MatrixElement/General/MEff2ff.cc
--- a/MatrixElement/General/MEff2ff.cc
+++ b/MatrixElement/General/MEff2ff.cc
@@ -1,615 +1,615 @@
 // -*- C++ -*-
 //
 // MEff2ff.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2017 The Herwig Collaboration
 //
 // Herwig 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 MEff2ff class.
 //
 
 #include "MEff2ff.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 using ThePEG::Helicity::VectorWaveFunction;
 using ThePEG::Helicity::ScalarWaveFunction;
 using ThePEG::Helicity::TensorWaveFunction;
 using ThePEG::Helicity::incoming;
 using ThePEG::Helicity::outgoing;
 
 
 void MEff2ff::doinit() {
   GeneralHardME::doinit();
   scalar_.resize(numberOfDiags());
   vector_.resize(numberOfDiags());
   tensor_.resize(numberOfDiags());
   four_  .resize(numberOfDiags());
   initializeMatrixElements(PDT::Spin1Half, PDT::Spin1Half, 
 			   PDT::Spin1Half, PDT::Spin1Half);
   for(size_t ix = 0;ix < numberOfDiags(); ++ix) {
     const HPDiagram & current = getProcessInfo()[ix];
     tcPDPtr offshell = current.intermediate;
     if(!offshell) {
       four_[ix] = dynamic_ptr_cast<AbstractFFFFVertexPtr>
 	(current.vertices.first);
     }
     else if(offshell->iSpin() == PDT::Spin0) {
       AbstractFFSVertexPtr vert1 = dynamic_ptr_cast<AbstractFFSVertexPtr>
 	(current.vertices.first);
       AbstractFFSVertexPtr vert2 = dynamic_ptr_cast<AbstractFFSVertexPtr>
 	(current.vertices.second);
       scalar_[ix] = make_pair(vert1, vert2);
     }
     else if(offshell->iSpin() == PDT::Spin1) {
       AbstractFFVVertexPtr vert1 = dynamic_ptr_cast<AbstractFFVVertexPtr>
 	(current.vertices.first);
       AbstractFFVVertexPtr vert2 = dynamic_ptr_cast<AbstractFFVVertexPtr>
 	(current.vertices.second);
       vector_[ix] = make_pair(vert1, vert2);
     }
     else if(offshell->iSpin() == PDT::Spin2) {
       AbstractFFTVertexPtr vert1 = dynamic_ptr_cast<AbstractFFTVertexPtr>
 	(current.vertices.first);
       AbstractFFTVertexPtr vert2 = dynamic_ptr_cast<AbstractFFTVertexPtr>
 	(current.vertices.second);
       tensor_[ix] = make_pair(vert1, vert2);
     }
   }
 }
 
 double MEff2ff::me2() const {
   for(unsigned int ix=0;ix<4;++ix) {
     spin_[ix].clear();
     sbar_[ix].clear();
     for(unsigned int ih=0;ih<2;++ih) {
       spin_[ix].push_back(SpinorWaveFunction   (rescaledMomenta()[ix],
 						mePartonData()[ix],
 						ih, ix<2 ? incoming : outgoing));
       sbar_[ix].push_back(SpinorBarWaveFunction(rescaledMomenta()[ix],
 						mePartonData()[ix],
 						ih, ix<2 ? incoming : outgoing));
     }
   }
   double full_me(0.);
   if( mePartonData()[0]->id() > 0 && mePartonData()[1]->id() < 0) {
     ffb2ffbHeME (full_me,true);
   }
   else if( mePartonData()[0]->id() > 0 && mePartonData()[1]->id() > 0 )
     ff2ffHeME(full_me,true);
   else if( mePartonData()[0]->id() < 0 && mePartonData()[1]->id() < 0 )
     fbfb2fbfbHeME(full_me,true);
   else
     assert(false);
 
 #ifndef NDEBUG
   if( debugME() ) debug(full_me);
 #endif
 
   return full_me;
 }
 
 ProductionMatrixElement 
 MEff2ff::ffb2ffbHeME(double & me2, bool first) const {
   const Energy2 q2(scale());
   // weights for the selection of the diagram
   vector<double> me(numberOfDiags(), 0.);
   // weights for the selection of the colour flow
   vector<double> flow(numberOfFlows(),0.);
   // flow over the helicities and diagrams
   for(unsigned int ifhel1 = 0; ifhel1 < 2; ++ifhel1) {
     for(unsigned int ifhel2 = 0; ifhel2 < 2; ++ifhel2) {
       for(unsigned int ofhel1 = 0; ofhel1 < 2; ++ofhel1) {
 	for(unsigned int ofhel2 = 0; ofhel2 < 2; ++ofhel2) {
 	  vector<Complex> flows(numberOfFlows(),0.);
 	  for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
 	    Complex diag(0.);
 	    const HPDiagram & current = getProcessInfo()[ix];
 	    tcPDPtr offshell = current.intermediate;
 	    if(current.channelType == HPDiagram::tChannel) {
 	      if(offshell->iSpin() == PDT::Spin0) {
 		if(current.ordered.second) {
 		  ScalarWaveFunction interS = scalar_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[3][ofhel2],sbar_[1][ifhel2]);
 		  diag = -scalar_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1],sbar_[2][ofhel1],interS);
 		}
 		else {
 		  ScalarWaveFunction interS = scalar_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[2][ofhel1],sbar_[1][ifhel2]);
 		  diag = scalar_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1],sbar_[3][ofhel2],interS);
 		}
 	      }
 	      else if(offshell->iSpin() == PDT::Spin1) {
 		if(current.ordered.second) {
 		  VectorWaveFunction interV = vector_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[3][ofhel2],sbar_[1][ifhel2]);
 		  diag = -vector_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1], sbar_[2][ofhel1], interV);
 		}
 		else {
 		  VectorWaveFunction interV = vector_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[2][ofhel1],sbar_[1][ifhel2]);
 		  diag = vector_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1], sbar_[3][ofhel2], interV);
  		}
 	      }
 	      else if(offshell->iSpin() == PDT::Spin2) {
 		if(current.ordered.second) {
 		  TensorWaveFunction interT = tensor_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[3][ofhel2],sbar_[1][ifhel2]);
 		  diag = -tensor_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1], sbar_[2][ofhel1], interT);
 		}
 		else {
 		  TensorWaveFunction interT = tensor_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[2][ofhel1],sbar_[1][ifhel2]);
 		  diag = tensor_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1], sbar_[3][ofhel2], interT);
  		}
 	      }
 	    }
 	    else if(current.channelType == HPDiagram::sChannel) {
 	      if(offshell->iSpin() == PDT::Spin0) {
 		ScalarWaveFunction interS = scalar_[ix].second->
 		  evaluate(q2, 1, offshell,spin_[3][ofhel2],sbar_[2][ofhel1]);
 		diag = scalar_[ix].first->
 		  evaluate(q2, spin_[0][ifhel1], sbar_[1][ifhel2], interS);
 	      }
 	      else if(offshell->iSpin() == PDT::Spin1) {
 		VectorWaveFunction interV = vector_[ix].second->
 		  evaluate(q2, 1, offshell,spin_[3][ofhel2],sbar_[2][ofhel1]);
 		diag = vector_[ix].first->
 		  evaluate(q2, spin_[0][ifhel1], sbar_[1][ifhel2], interV);
 	      }
 	      else if(offshell->iSpin() == PDT::Spin2) {
 		TensorWaveFunction interT = tensor_[ix].second->
 		  evaluate(q2, 1, offshell,spin_[3][ofhel2],sbar_[2][ofhel1]);
 		diag = tensor_[ix].first->
 		  evaluate(q2, spin_[0][ifhel1], sbar_[1][ifhel2], interT);
 	      }
 	    }
 	    else {
 	      diag = four_[ix]->evaluate(q2,spin_[0][ifhel1], sbar_[1][ifhel2],
 					 spin_[3][ofhel2],sbar_[2][ofhel1]);
 	    }
 	    me[ix] += norm(diag);
 	    diagramME()[ix](ifhel1, ifhel2, ofhel1, ofhel2) = diag;
 	    //Compute flows
 	    for(size_t iy = 0; iy < current.colourFlow.size(); ++iy) {
 	      assert(current.colourFlow[iy].first<flows.size());
 	      flows[current.colourFlow[iy].first] += 
 		current.colourFlow[iy].second * diag;
 	    }
 	  }
 	  // MEs for the different colour flows
 	  for(unsigned int iy = 0; iy < numberOfFlows(); ++iy) 
 	    flowME()[iy](ifhel1, ifhel2, ofhel1, ofhel2) = flows[iy];
 	  //Now add flows to me2 with appropriate colour factors
 	  for(size_t ii = 0; ii < numberOfFlows(); ++ii) {
 	    for(size_t ij = 0; ij < numberOfFlows(); ++ij) {
 	      me2 += getColourFactors()[ii][ij]*(flows[ii]*conj(flows[ij])).real();
 	    }
 	  }
 	  // contribution to the colour flow
 	  for(unsigned int ii = 0; ii < numberOfFlows(); ++ii) {
 	    flow[ii] += getColourFactors()[ii][ii]*norm(flows[ii]);
 	  }
 	}
       }
     }
   }
   // if not computing the cross section return the selected colour flow
   if(!first) return flowME()[colourFlow()];
   me2 = selectColourFlow(flow,me,me2);
   return flowME()[colourFlow()];
 }
 
 ProductionMatrixElement 
 MEff2ff:: ff2ffHeME(double & me2, bool first) const {
   const Energy2 q2(scale());
   // weights for the selection of the diagram
   vector<double> me(numberOfDiags(), 0.);
   // weights for the selection of the colour flow
   vector<double> flow(numberOfFlows(),0.);
   // flow over the helicities and diagrams
   for(unsigned int ifhel1 = 0; ifhel1 < 2; ++ifhel1) {
     for(unsigned int ifhel2 = 0; ifhel2 < 2; ++ifhel2) {
       for(unsigned int ofhel1 = 0; ofhel1 < 2; ++ofhel1) {
 	for(unsigned int ofhel2 = 0; ofhel2 < 2; ++ofhel2) {
 	  vector<Complex> flows(numberOfFlows(),0.);
 	  for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
 	    Complex diag(0.);
 	    const HPDiagram & current = getProcessInfo()[ix];
 	    tcPDPtr offshell = current.intermediate;
+	    if(offshell->CC()) offshell=offshell->CC();
 	    if(current.channelType == HPDiagram::tChannel) {
 	      if(offshell->iSpin() == PDT::Spin0) {
 		if(current.ordered.second) {
 		  ScalarWaveFunction interS = scalar_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[1][ifhel2],sbar_[3][ofhel2]);
 		  diag = scalar_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1], sbar_[2][ofhel1], interS);
 		}
 		else {
 		  ScalarWaveFunction interS = scalar_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[1][ifhel2],sbar_[2][ofhel1]);
 		  diag = -scalar_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1], sbar_[3][ofhel2], interS);
 		}
 	      }
 	      else if(offshell->iSpin() == PDT::Spin1) {
 		if(current.ordered.second) {
 		  VectorWaveFunction interV = vector_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[1][ifhel2],sbar_[3][ofhel2]);
 		  diag = vector_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1], sbar_[2][ofhel1], interV);
 		}
 		else {
 		  VectorWaveFunction interV = vector_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[1][ifhel2],sbar_[2][ofhel1]);
 		  diag = -vector_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1], sbar_[3][ofhel2], interV);
 		}
 	      }
 	      else if(offshell->iSpin() == PDT::Spin2) {
 		if(current.ordered.second) {
 		  TensorWaveFunction interT = tensor_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[1][ifhel2],sbar_[3][ofhel2]);
 		  diag = tensor_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1], sbar_[2][ofhel1], interT);
 		}
 		else {
 		  TensorWaveFunction interT = tensor_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[1][ifhel2],sbar_[2][ofhel1]);
 		  diag = -tensor_[ix].first->
 		    evaluate(q2, spin_[0][ifhel1], sbar_[3][ofhel2], interT);
 		}
 	      }
 	    }
 	    else if(current.channelType == HPDiagram::sChannel) {
-	      if(offshell->CC()) offshell=offshell->CC();
 	      if(offshell->iSpin() == PDT::Spin0) {
 		ScalarWaveFunction interS = scalar_[ix].second->
 		  evaluate(q2, 1, offshell,spin_[3][ofhel2],sbar_[2][ofhel1]);
 		diag = scalar_[ix].first->
 		  evaluate(q2, spin_[0][ifhel1], sbar_[1][ifhel2], interS);
 	      }
 	      else if(offshell->iSpin() == PDT::Spin1) {
 		VectorWaveFunction interV = vector_[ix].second->
 		  evaluate(q2, 1, offshell,spin_[3][ofhel2],sbar_[2][ofhel1]);
 		diag = vector_[ix].first->
 		  evaluate(q2, spin_[0][ifhel1], sbar_[1][ifhel2], interV);
 	      }
 	      else if(offshell->iSpin() == PDT::Spin2) {
 		TensorWaveFunction interT = tensor_[ix].second->
 		  evaluate(q2, 1, offshell,spin_[3][ofhel2],sbar_[2][ofhel1]);
 		diag = tensor_[ix].first->
 		  evaluate(q2, spin_[0][ifhel1], sbar_[1][ifhel2], interT);
 	      }
 	    }
 	    else if(current.channelType == HPDiagram::fourPoint) {
 	      diag= four_[ix]->evaluate(q2,spin_[0][ifhel1],sbar_[2][ofhel1],
 					spin_[1][ifhel2],sbar_[3][ofhel2]);
 	    }
 	    me[ix] += norm(diag);
 	    diagramME()[ix](ifhel1, ifhel2, ofhel1, ofhel2) = diag;
 	    //Compute flows
 	    for(size_t iy = 0; iy < current.colourFlow.size(); ++iy) {
 	      assert(current.colourFlow[iy].first<flows.size());
 	      flows[current.colourFlow[iy].first] += 
 		current.colourFlow[iy].second * diag;
 	    }
 	  }
 	  // MEs for the different colour flows
 	  for(unsigned int iy = 0; iy < numberOfFlows(); ++iy) 
 	    flowME()[iy](ifhel1, ifhel2, ofhel1, ofhel2) = flows[iy];
 	  //Now add flows to me2 with appropriate colour factors
 	  for(size_t ii = 0; ii < numberOfFlows(); ++ii)
 	    for(size_t ij = 0; ij < numberOfFlows(); ++ij)
 	      me2 += getColourFactors()[ii][ij]*(flows[ii]*conj(flows[ij])).real();
 	  // contribution to the colour flow
 	  for(unsigned int ii = 0; ii < numberOfFlows(); ++ii) {
 	    flow[ii] += getColourFactors()[ii][ii]*norm(flows[ii]);
 	  }
 	}
       }
     }
   }
   // if not computing the cross section return the selected colour flow
   if(!first) return flowME()[colourFlow()];
   me2 = selectColourFlow(flow,me,me2);
   return flowME()[colourFlow()];
 }
 
 ProductionMatrixElement
 MEff2ff::fbfb2fbfbHeME(double & me2, bool first) const {
   const Energy2 q2(scale());
   // weights for the selection of the diagram
   vector<double> me(numberOfDiags(), 0.);
   // weights for the selection of the colour flow
   vector<double> flow(numberOfFlows(),0.);
   // flow over the helicities and diagrams
   for(unsigned int ifhel1 = 0; ifhel1 < 2; ++ifhel1) {
     for(unsigned int ifhel2 = 0; ifhel2 < 2; ++ifhel2) {
       for(unsigned int ofhel1 = 0; ofhel1 < 2; ++ofhel1) {
 	for(unsigned int ofhel2 = 0; ofhel2 < 2; ++ofhel2) {
 	  vector<Complex> flows(numberOfFlows(),0.);
 	  for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
 	    Complex diag(0.);
 	    const HPDiagram & current = getProcessInfo()[ix];
 	    tcPDPtr offshell = current.intermediate;
+	    if(offshell->CC()) offshell=offshell->CC();
 	    if(current.channelType == HPDiagram::tChannel) {
 	      if(offshell->iSpin() == PDT::Spin0) {
 		if(current.ordered.second) {
 		  ScalarWaveFunction interS = scalar_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[3][ofhel2],sbar_[1][ifhel2]);
 		  diag = scalar_[ix].first->
 		    evaluate(q2, spin_[2][ofhel1], sbar_[0][ifhel1], interS);
 		}
 		else {
 		  ScalarWaveFunction interS = scalar_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[2][ofhel1],sbar_[1][ifhel2]);
 		  diag = -scalar_[ix].first->
 		    evaluate(q2, spin_[3][ofhel2], sbar_[0][ifhel1], interS);
 		}
 	      }
 	      else if(offshell->iSpin() == PDT::Spin1) {
 		if(current.ordered.second) {
 		  VectorWaveFunction interV = vector_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[3][ofhel2],sbar_[1][ifhel2]);
 		  diag = vector_[ix].first->
 		    evaluate(q2, spin_[2][ofhel1], sbar_[0][ifhel1], interV);
 		}
 		else {
 		  VectorWaveFunction interV = vector_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[2][ofhel1],sbar_[1][ifhel2]);
 		  diag = -vector_[ix].first->
 		    evaluate(q2, spin_[3][ofhel2], sbar_[0][ifhel1], interV);
 		}
 	      }
 	      else if(offshell->iSpin() == PDT::Spin2) {
 		if(current.ordered.second) {
 		  TensorWaveFunction interT = tensor_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[3][ofhel2],sbar_[1][ifhel2]);
 		  diag = tensor_[ix].first->
 		    evaluate(q2, spin_[2][ofhel1], sbar_[0][ifhel1], interT);
 		}
 		else {
 		  TensorWaveFunction interT = tensor_[ix].second->
 		    evaluate(q2, 3, offshell,spin_[2][ofhel1],sbar_[1][ifhel2]);
 		  diag = -tensor_[ix].first->
 		    evaluate(q2, spin_[3][ofhel2], sbar_[0][ifhel1], interT);
 		}
 	      }
 	    }
 	    else if(current.channelType == HPDiagram::sChannel) {
-	      if(offshell->CC()) offshell=offshell->CC();
 	      if(offshell->iSpin() == PDT::Spin0) {
 		ScalarWaveFunction interS = scalar_[ix].second->
 		  evaluate(q2, 1, offshell,spin_[3][ofhel2],sbar_[2][ofhel1]);
 		diag = scalar_[ix].first->
 		  evaluate(q2, spin_[0][ifhel1], sbar_[1][ifhel2], interS);
 	      }
 	      else if(offshell->iSpin() == PDT::Spin1) {
 		VectorWaveFunction interV = vector_[ix].second->
 		  evaluate(q2, 1, offshell,spin_[3][ofhel2],sbar_[2][ofhel1]);
 		diag = vector_[ix].first->
 		  evaluate(q2, spin_[0][ifhel1], sbar_[1][ifhel2], interV);
 	      }
 	      else if(offshell->iSpin() == PDT::Spin2) {
 		TensorWaveFunction interT = tensor_[ix].second->
 		  evaluate(q2, 1, offshell,spin_[3][ofhel2],sbar_[2][ofhel1]);
 		diag = tensor_[ix].first->
 		  evaluate(q2, spin_[0][ifhel1], sbar_[1][ifhel2], interT);
 	      }
 	    }
 	    else if(current.channelType == HPDiagram::fourPoint) {
 	      diag= four_[ix]->evaluate(q2,spin_[2][ofhel1],sbar_[0][ifhel1],
 					spin_[3][ofhel2],sbar_[1][ifhel2]);
 	    }
 	    me[ix] += norm(diag);
 	    diagramME()[ix](ifhel1, ifhel2, ofhel1, ofhel2) = diag;
 	    //Compute flows
 	    for(size_t iy = 0; iy < current.colourFlow.size(); ++iy) {
 	      assert(current.colourFlow[iy].first<flows.size());
 	      flows[current.colourFlow[iy].first] += 
 		current.colourFlow[iy].second * diag;
 	    }
 	  }
 	  // MEs for the different colour flows
 	  for(unsigned int iy = 0; iy < numberOfFlows(); ++iy) 
 	    flowME()[iy](ifhel1, ifhel2, ofhel1, ofhel2) = flows[iy];
 	  //Now add flows to me2 with appropriate colour factors
 	  for(size_t ii = 0; ii < numberOfFlows(); ++ii)
 	    for(size_t ij = 0; ij < numberOfFlows(); ++ij)
 	      me2 += getColourFactors()[ii][ij]*(flows[ii]*conj(flows[ij])).real();
 	  // contribution to the colour flow
 	  for(unsigned int ii = 0; ii < numberOfFlows(); ++ii) {
 	    flow[ii] += getColourFactors()[ii][ii]*norm(flows[ii]);
 	  }
 	}
       }
     }
   }
   // if not computing the cross section return the selected colour flow
   if(!first) return flowME()[colourFlow()];
   me2 = selectColourFlow(flow,me,me2);
   return flowME()[colourFlow()];
 }
 
 void MEff2ff::constructVertex(tSubProPtr subp) {
   // Hard process external particles
   ParticleVector hardpro = hardParticles(subp);
   //Need to use rescale momenta to calculate matrix element
   setRescaledMomenta(hardpro);
   for(unsigned int ix=0;ix<4;++ix) {
     spin_[ix].clear();
     sbar_[ix].clear();
     SpinorWaveFunction   (spin_[ix],hardpro[ix],
 			  ix<2 ? incoming : outgoing,ix>1);
     SpinorBarWaveFunction(sbar_[ix],hardpro[ix],
 			  ix<2 ? incoming : outgoing,ix>1);
   }
   double dummy(0.);
   //pick which process we are doing
   if( hardpro[0]->id() > 0) {
     //ffbar->ffbar
     if( hardpro[1]->id() < 0 ) {
       ProductionMatrixElement prodME = ffb2ffbHeME(dummy,false);
       createVertex(prodME,hardpro);
     }
     //ff2ff
     else {
       ProductionMatrixElement prodME = ff2ffHeME(dummy,false);
       createVertex(prodME,hardpro);
     }
   } 
   //fbarfbar->fbarfbar
   else {
     ProductionMatrixElement prodME = fbfb2fbfbHeME(dummy,false);
     createVertex(prodME,hardpro);
   }
   
 #ifndef NDEBUG
   if( debugME() ) debug(dummy);
 #endif
 
 }
 
 void MEff2ff::persistentOutput(PersistentOStream & os) const {
   os << scalar_ << vector_ << tensor_ << four_;
 }
 
 void MEff2ff::persistentInput(PersistentIStream & is, int) {
   is >> scalar_ >> vector_ >> tensor_ >> four_;
   initializeMatrixElements(PDT::Spin1Half, PDT::Spin1Half, 
 			   PDT::Spin1Half, PDT::Spin1Half);
 }
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<MEff2ff,GeneralHardME>
 describeHerwigMEff2ff("Herwig::MEff2ff", "Herwig.so");
 
 void MEff2ff::Init() {
 
   static ClassDocumentation<MEff2ff> documentation
     ("This is the implementation of the matrix element for fermion-"
      "antifermion -> fermion-antifermion.");
 
 }
 
 void MEff2ff::debug(double me2) const {
   if( !generator()->logfile().is_open() ) return;
   long id1 = mePartonData()[0]->id();
   long id2 = mePartonData()[1]->id();
   long id3 = mePartonData()[2]->id();
   long id4 = mePartonData()[3]->id();
   long aid1 = abs(mePartonData()[0]->id());
   long aid2 = abs(mePartonData()[1]->id());
   long aid3 = abs(mePartonData()[2]->id());
   long aid4 = abs(mePartonData()[3]->id());
   if( (aid1 != 1 && aid1 != 2) || (aid2 != 1 && aid2 != 2) ) return;
   double analytic(0.);
   if( id3 == id4 && id3 == 1000021 ) {
     tcSMPtr sm = generator()->standardModel();
     double gs4 = sqr( 4.*Constants::pi*sm->alphaS(scale()) );
     int Nc = sm->Nc();
     double Cf = (sqr(Nc) - 1.)/2./Nc;
     Energy2 mgo2 = meMomenta()[3].m2();
     long squark = (aid1 == 1) ? 1000001 : 1000002;
     Energy2 muL2 = sqr(getParticleData(squark)->mass());
     Energy2 deltaL = muL2 - mgo2;
     Energy2 muR2 = sqr(getParticleData(squark + 1000000)->mass());
     Energy2 deltaR = muR2 - mgo2;
     Energy2 s(sHat());
     Energy2 m3s = meMomenta()[2].m2();
     Energy2 m4s = meMomenta()[3].m2();
     Energy4 spt2 = uHat()*tHat() - m3s*m4s;
     Energy2 t3(tHat() - m3s), u4(uHat() - m4s);
     
     double Cl = 2.*spt2*( (u4*u4 - deltaL*deltaL) + (t3*t3 - deltaL*deltaL)
 			  - (s*s/Nc/Nc) )/s/s/(u4 - deltaL)/(t3 - deltaL);
     Cl += deltaL*deltaL*( (1./sqr(t3 - deltaL)) + (1./sqr(u4 - deltaL))
 			  - ( sqr( (1./(t3 - deltaL)) - 
 				   (1./(u4 - deltaL)) )/Nc/Nc ) );
     
     double Cr = 2.*spt2*( (u4*u4 - deltaR*deltaR) + (t3*t3 - deltaR*deltaR)
 			  - (s*s/Nc/Nc) )/s/s/(u4 - deltaR)/(t3 - deltaR);
     Cr += deltaR*deltaR*( (1./sqr(t3 - deltaR)) + (1./sqr(u4 - deltaR))
 			  - ( sqr( (1./(t3 - deltaR)) 
 				   - (1./(u4 - deltaR)) )/Nc/Nc ) );
     analytic = gs4*Cf*(Cl + Cr)/4.;
   }
   else if( (aid3 == 5100001 || aid3 == 5100002 ||
 	    aid3 == 6100001 || aid3 == 6100002) &&
 	   (aid4 == 5100001 || aid4 == 5100002 ||
 	    aid4 == 6100001 || aid4 == 6100002) ) {
     tcSMPtr sm = generator()->standardModel();
     double gs4 = sqr( 4.*Constants::pi*sm->alphaS(scale()) );
     Energy2 s(sHat());
     Energy2 mf2 = meMomenta()[2].m2();
     Energy2 t3(tHat() - mf2), u4(uHat() - mf2);
     Energy4 s2(sqr(s)), t3s(sqr(t3)), u4s(sqr(u4));
     
     bool iflav = (aid2 - aid1 == 0);
     int alpha(aid3/1000000), beta(aid4/1000000);
     bool oflav = ((aid3 - aid1) % 10  == 0);
     if( alpha != beta ) {
       if( ( id1 > 0 && id2 > 0) ||
 	  ( id1 < 0 && id2 < 0) ) {
 	if( iflav )
 	  analytic = gs4*( mf2*(2.*s2*s/t3s/u4s - 4.*s/t3/u4) 
 			   + 2.*sqr(s2)/t3s/u4s - 8.*s2/t3/u4 + 5. )/9.;
 	else
 	  analytic = gs4*( -2.*mf2*(1./t3 + u4/t3s) + 0.5 + 2.*u4s/t3s)/9.;
       }
       else
 	analytic = gs4*( 2.*mf2*(1./t3 + u4/t3s) + 5./2. + 4.*u4/t3 
 			 + 2.*u4s/t3s)/9.;
     }
     else {
       if( ( id1 > 0 && id2 > 0) ||
 	  ( id1 < 0 && id2 < 0) ) {
 	if( iflav ) {
 	  analytic = gs4*( mf2*(6.*t3/u4s + 6.*u4/t3s - s/t3/u4) 
 			   + 2.*(3.*t3s/u4s + 3.*u4s/t3s 
 				 + 4.*s2/t3/u4 - 5.) )/27.;
 	}
 	else
 	  analytic = 2.*gs4*( -mf2*s/t3s + 0.25 + s2/t3s )/9.;
       }
       else {
 	if( iflav ) {
 	  if( oflav )	  
 	    analytic = gs4*( 2.*mf2*(4./s + s/t3s - 1./t3) + 23./6.+ 2.*s2/t3s
 			     + 8.*s/3./t3 + 6.*t3/s + 8.*t3s/s2 )/9.;
 	  else
 	    analytic = 4.*gs4*( 2.*mf2/s + (t3s + u4s)/s2)/9.;
 	}
 	else 
 	  analytic = gs4*(4.*mf2*s/t3s + 5. + 4.*s2/t3s + 8.*s/t3 )/18.;
       }
     }
     if( id3 == id4 ) analytic /= 2.;
 
   }
   else return;
   double diff = abs(analytic - me2);
   if( diff  > 1e-4 ) {
     generator()->log() 
       << mePartonData()[0]->PDGName() << ","
       << mePartonData()[1]->PDGName() << "->"
       << mePartonData()[2]->PDGName() << ","
       << mePartonData()[3]->PDGName() << "   difference: " 
       << setprecision(10) << diff  << "  ratio: " << analytic/me2 
       << '\n';
   }
 }
diff --git a/Models/General/HardProcessConstructor.cc b/Models/General/HardProcessConstructor.cc
--- a/Models/General/HardProcessConstructor.cc
+++ b/Models/General/HardProcessConstructor.cc
@@ -1,956 +1,956 @@
 // -*- C++ -*-
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the HardProcessConstructor class.
 //
 
 #include "HardProcessConstructor.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 void HardProcessConstructor::persistentOutput(PersistentOStream & os) const {
   os << debug_ << subProcess_ << model_;
 }
 
 void HardProcessConstructor::persistentInput(PersistentIStream & is, int) {
   is >> debug_ >> subProcess_ >> model_;
 }
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeAbstractClass<HardProcessConstructor,Interfaced>
 describeHerwigHardProcessConstructor("Herwig::HardProcessConstructor", "Herwig.so");
 
 void HardProcessConstructor::Init() {
 
   static ClassDocumentation<HardProcessConstructor> documentation
     ("Base class for implementation of the automatic generation of hard processes");
 
   static Switch<HardProcessConstructor,bool> interfaceDebugME
     ("DebugME",
      "Print comparison with analytical ME",
      &HardProcessConstructor::debug_, false, false, false);
   static SwitchOption interfaceDebugMEYes
     (interfaceDebugME,
      "Yes",
      "Print the debug information",
      true);
   static SwitchOption interfaceDebugMENo
     (interfaceDebugME,
      "No",
      "Do not print the debug information",
      false);
 
 }
 
 void HardProcessConstructor::doinit() {
   Interfaced::doinit();
   EGPtr eg = generator();
   model_ = dynamic_ptr_cast<HwSMPtr>(eg->standardModel());
   if(!model_)
     throw InitException() << "HardProcessConstructor:: doinit() - "
 			  << "The model pointer is null!"
 			  << Exception::abortnow;
   if(!eg->eventHandler()) {
     throw
       InitException() << "HardProcessConstructor:: doinit() - "
 		      << "The eventHandler pointer was null therefore "
 		      << "could not get SubProcessHandler pointer " 
 		      << Exception::abortnow;
   }
   string subProcessName = 
     eg->preinitInterface(eg->eventHandler(), "SubProcessHandlers", "get","");
   subProcess_ = eg->getObject<SubProcessHandler>(subProcessName);
   if(!subProcess_) {
     ostringstream s;
     s << "HardProcessConstructor:: doinit() - "
       << "There was an error getting the SubProcessHandler "
       << "from the current event handler. ";
     generator()->logWarning( Exception(s.str(), Exception::warning) );
   }
 }
 
 GeneralHardME::ColourStructure HardProcessConstructor::
 colourFlow(const tcPDVector & extpart) const {
   PDT::Colour ina = extpart[0]->iColour();
   PDT::Colour inb = extpart[1]->iColour();
   PDT::Colour outa = extpart[2]->iColour();
   PDT::Colour outb = extpart[3]->iColour();
 
   // incoming colour neutral
   if(ina == PDT::Colour0 && inb == PDT::Colour0) {
     if( outa == PDT::Colour0 && outb == PDT::Colour0 ) {
       return GeneralHardME::Colour11to11;
     }
     else if( outa == PDT::Colour3 && outb == PDT::Colour3bar ) {
       return GeneralHardME::Colour11to33bar;
     } 
     else if( outa == PDT::Colour8 && outb == PDT::Colour8 ) {
       return GeneralHardME::Colour11to88;
     } 
     else
       assert(false);
   }
   // incoming 3 3
   else if(ina == PDT::Colour3 && inb == PDT::Colour3 ) {
     if( outa == PDT::Colour3 && outb == PDT::Colour3 ) {
       return GeneralHardME::Colour33to33;
     }
     else if( outa == PDT::Colour6 && outb == PDT::Colour0 ) {
       return GeneralHardME::Colour33to61;
     }
     else if( outa == PDT::Colour0 && outb == PDT::Colour6 ) {
       return GeneralHardME::Colour33to16;
     }
     else if ( outa == PDT::Colour0 && outb == PDT::Colour3bar) {
       return GeneralHardME::Colour33to13bar;
     }
     else if ( outb == PDT::Colour0 && outa == PDT::Colour3bar) {
       return GeneralHardME::Colour33to3bar1;
     }
     else if ( outa == PDT::Colour8 && outb == PDT::Colour3bar) {
       return GeneralHardME::Colour33to83bar;
     }
     else if ( outb == PDT::Colour8 && outa == PDT::Colour3bar) {
       return GeneralHardME::Colour33to3bar8;
     }
     else
       assert(false);
   }
   // incoming 3bar 3bar
   else if(ina == PDT::Colour3bar && inb == PDT::Colour3bar ) {
     if( outa == PDT::Colour3bar && outb == PDT::Colour3bar ) {
       return GeneralHardME::Colour3bar3barto3bar3bar;
     }
     else if( outa == PDT::Colour6bar && outb == PDT::Colour0) {
       return GeneralHardME::Colour3bar3barto6bar1;
     }
     else if ( outa == PDT::Colour0 && outb == PDT::Colour6bar ) {
       return GeneralHardME::Colour3bar3barto16bar;
     }
     else if ( outa == PDT::Colour0 && outb == PDT::Colour3) {
       return GeneralHardME::Colour3bar3barto13;
     }
     else if ( outb == PDT::Colour0 && outa == PDT::Colour3) {
       return GeneralHardME::Colour3bar3barto31;
     }
     else if ( outa == PDT::Colour8 && outb == PDT::Colour3) {
       return GeneralHardME::Colour3bar3barto83;
     }
     else if ( outb == PDT::Colour8 && outa == PDT::Colour3) {
       return GeneralHardME::Colour3bar3barto38;
     }
     else
       assert(false);
   }
   // incoming 3 3bar
   else if(ina == PDT::Colour3 && inb == PDT::Colour3bar ) {
     if( outa == PDT::Colour0 && outb == PDT::Colour0 ) {
       return GeneralHardME::Colour33barto11;
     }
     else if( outa == PDT::Colour3 && outb == PDT::Colour3bar ) {
       return GeneralHardME::Colour33barto33bar;
     }
     else if( outa == PDT::Colour8 && outb == PDT::Colour8 ) {
       return GeneralHardME::Colour33barto88;
     }
     else if( outa == PDT::Colour8 && outb == PDT::Colour0 ) {
       return GeneralHardME::Colour33barto81;
     }
     else if( outa == PDT::Colour0 && outb == PDT::Colour8 ) {
       return GeneralHardME::Colour33barto18;
     }
     else if( outa == PDT::Colour6 && outb == PDT::Colour6bar) {
       return GeneralHardME::Colour33barto66bar;
     }
     else if( outa == PDT::Colour6bar && outb == PDT::Colour6) {
       return GeneralHardME::Colour33barto6bar6;
     }
     else
       assert(false);
   }
   // incoming 88
   else if(ina == PDT::Colour8 && inb == PDT::Colour8 ) {
     if( outa == PDT::Colour0 && outb == PDT::Colour0 ) {
       return GeneralHardME::Colour88to11;
     }
     else if( outa == PDT::Colour3 && outb == PDT::Colour3bar ) {
       return GeneralHardME::Colour88to33bar;
     }
     else if( outa == PDT::Colour8 && outb == PDT::Colour8 ) {
       return GeneralHardME::Colour88to88;
     }
     else if( outa == PDT::Colour8 && outb == PDT::Colour0 ) {
       return GeneralHardME::Colour88to81;
     }
     else if( outa == PDT::Colour0 && outb == PDT::Colour8 ) {
       return GeneralHardME::Colour88to18;
     }
     else if( outa == PDT::Colour6 && outb == PDT::Colour6bar ) {
       return GeneralHardME::Colour88to66bar;
     }    
     else
       assert(false);
   }
   // incoming 38
   else if(ina == PDT::Colour3 && inb == PDT::Colour8 ) {
     if(outa == PDT::Colour3 && outb == PDT::Colour0) {
       return GeneralHardME::Colour38to31;
     }
     else if(outa == PDT::Colour0 && outb == PDT::Colour3) {
       return GeneralHardME::Colour38to13;
     }
     else if(outa == PDT::Colour3 && outb == PDT::Colour8) {
       return GeneralHardME::Colour38to38;
     }
     else if(outa == PDT::Colour8 && outb == PDT::Colour3) {
       return GeneralHardME::Colour38to83;
     }
     else if(outa == PDT::Colour3bar && outb == PDT::Colour6){
       return GeneralHardME::Colour38to3bar6;
     }
     else if(outa == PDT::Colour6 && outb == PDT::Colour3bar) {
       return GeneralHardME::Colour38to63bar;
     }
     else if(outa == PDT::Colour3bar && outb == PDT::Colour3bar) {
       return GeneralHardME::Colour38to3bar3bar;
     }
     else
       assert(false);
   }
   // incoming 3bar8
   else if(ina == PDT::Colour3bar && inb == PDT::Colour8 ) {   
     if(outa == PDT::Colour3bar && outb == PDT::Colour0 ) {
       return GeneralHardME::Colour3bar8to3bar1;
     }
     else if(outa == PDT::Colour0 && outb == PDT::Colour3bar) {
       return GeneralHardME::Colour3bar8to13bar;
     }
     else if(outa == PDT::Colour3bar && outb == PDT::Colour8 ) {
       return GeneralHardME::Colour3bar8to3bar8;
     }
     else if(outa == PDT::Colour8 && outb == PDT::Colour3bar) {
       return GeneralHardME::Colour3bar8to83bar;
     }
     else if(outa == PDT::Colour3 && outb == PDT::Colour3) {
       return GeneralHardME::Colour3bar8to33;
     }
     else
       assert(false);
   }
   // unknown colour flow
   else 
     assert(false);
   return GeneralHardME::UNDEFINED;
 }
 
 
 void HardProcessConstructor::fixFSOrder(HPDiagram & diag) {
   tcPDPtr psa = getParticleData(diag.incoming.first);
   tcPDPtr psb = getParticleData(diag.incoming.second);
   tcPDPtr psc = getParticleData(diag.outgoing.first);
   tcPDPtr psd = getParticleData(diag.outgoing.second);
 
   //fix a spin order
   if( psc->iSpin() < psd->iSpin() ) {
     swap(diag.outgoing.first, diag.outgoing.second);
     if(diag.channelType == HPDiagram::tChannel) {
       diag.ordered.second = !diag.ordered.second;
     }
     return;
   }
   
   if( psc->iSpin() == psd->iSpin() && 
       psc->id() < 0 && psd->id() > 0 ) {
     swap(diag.outgoing.first, diag.outgoing.second);
     if(diag.channelType == HPDiagram::tChannel) {
       diag.ordered.second = !diag.ordered.second;
     }
     return;
   }
 }
 
 void HardProcessConstructor::assignToCF(HPDiagram & diag) {
   if(diag.channelType == HPDiagram::tChannel) {
     if(diag.ordered.second) tChannelCF(diag);
     else                    uChannelCF(diag);
   }
   else if(diag.channelType == HPDiagram::sChannel) {
     sChannelCF(diag);
   }
   else if (diag.channelType == HPDiagram::fourPoint) {
     fourPointCF(diag);
   }
   else 
     assert(false);
 }
 
 void HardProcessConstructor::tChannelCF(HPDiagram & diag) {
   tcPDPtr ia = getParticleData(diag.incoming.first );
   tcPDPtr ib = getParticleData(diag.incoming.second);
   tcPDPtr oa = getParticleData(diag.outgoing.first );
   tcPDPtr ob = getParticleData(diag.outgoing.second);
   PDT::Colour ina  = ia->iColour();
   PDT::Colour inb  = ib->iColour();
   PDT::Colour outa = oa->iColour();
   PDT::Colour outb = ob->iColour();
   vector<CFPair> cfv(1, make_pair(0, 1.));
   if(diag.intermediate->iColour() == PDT::Colour0) {
     if(ina==PDT::Colour0) {
       cfv[0] = make_pair(0, 1);
     }
     else if(ina==PDT::Colour3 || ina==PDT::Colour3bar) {
       if( inb == PDT::Colour0 ) {
 	cfv[0] = make_pair(0, 1);
       }
       else if(inb==PDT::Colour3 || outb==PDT::Colour3bar) {
 	cfv[0] = make_pair(2, 1);
       }
       else if(inb==PDT::Colour8) {
 	cfv[0] = make_pair(2, 1);
       }
     }
     else if(ina==PDT::Colour8) {
       if( inb == PDT::Colour0 ) {
 	cfv[0] = make_pair(0, 1);
       }
       else if(inb==PDT::Colour3 || outb==PDT::Colour3bar) {
 	cfv[0] = make_pair(2, 1);
       }
       else if(inb==PDT::Colour8) {
 	cfv[0] = make_pair(7, -1);
       }
     }
   }
   else if(diag.intermediate->iColour() == PDT::Colour8) {
     if(ina==PDT::Colour8&&outa==PDT::Colour8&&
        inb==PDT::Colour8&&outb==PDT::Colour8) {
       cfv[0]=make_pair(2,  2.);
       cfv.push_back(make_pair(3, -2.));
       cfv.push_back(make_pair(1, -2.));
       cfv.push_back(make_pair(4,  2.));
     }
     else if(ina==PDT::Colour8&&outa==PDT::Colour0&&
 	    inb==PDT::Colour8&&outb==PDT::Colour8&&
 	    (oa->iSpin()==PDT::Spin0||oa->iSpin()==PDT::Spin1Half||
 	     oa->iSpin()==PDT::Spin3Half)) {
       cfv[0] = make_pair(0,-1);
     }
     else if(ina==PDT::Colour8&&outa==PDT::Colour8&&
 	    inb==PDT::Colour8&&outb==PDT::Colour0&&
 	    (ob->iSpin()==PDT::Spin0||ob->iSpin()==PDT::Spin1Half||
 	     ob->iSpin()==PDT::Spin3Half)) {
       cfv[0] = make_pair(0,-1);
     }
   } 
   else if(diag.intermediate->iColour() == PDT::Colour3 ||
 	  diag.intermediate->iColour() == PDT::Colour3bar) {
     if(outa == PDT::Colour0 || outb == PDT::Colour0) {
       if( outa == PDT::Colour6    || outb == PDT::Colour6   ||
 	  outa == PDT::Colour6bar || outb == PDT::Colour6bar) {
 	cfv[0] = make_pair(0,0.5);
 	cfv.push_back(make_pair(1,0.5));
       }
       else if ((ina==PDT::Colour3 && inb == PDT::Colour3 &&
 		(outa == PDT::Colour3bar || outb == PDT::Colour3bar))||
 	       (ina==PDT::Colour3bar && inb == PDT::Colour3bar &&
 		(outa == PDT::Colour3 || outb == PDT::Colour3 ))) {
 	cfv[0] = make_pair(0,1.);
       }
       else {
 	cfv[0] = make_pair(0,1.);
       }
     }
     else if(outa==PDT::Colour6 && outb==PDT::Colour3bar) {
       cfv[0] = make_pair(4,1.);
       cfv.push_back(make_pair(5,1.));
     }
     else if(outa==PDT::Colour6 && outb==PDT::Colour6bar) {
       cfv[0] = make_pair(4, 1.);
       for(unsigned int ix=5;ix<8;++ix)
 	cfv.push_back(make_pair(ix,1.));
     }
     else if(outa==PDT::Colour6 || outa ==PDT::Colour6bar ||
 	    outb==PDT::Colour6 || outb ==PDT::Colour6bar ) {
       assert(false);
     }
     else if(ina==PDT::Colour3    && inb==PDT::Colour3    ) {
       if((outa==PDT::Colour0 && outb==PDT::Colour3bar)||
 	 (outb==PDT::Colour0 && outa==PDT::Colour3bar))
 	cfv[0] = make_pair(0,1.);
       else if((outa==PDT::Colour8 && outb==PDT::Colour3bar)||
 	      (outb==PDT::Colour8 && outa==PDT::Colour3bar)) {
-	double sign = diag.intermediate->iSpin()==PDT::Spin0 ? -1. : 1.;
-	cfv[0] = make_pair(1,sign);
+	cfv[0] = make_pair(1,1.);
       }
     }
     else if(ina==PDT::Colour3bar && inb==PDT::Colour3bar ) {
       if((outa==PDT::Colour0 && outb==PDT::Colour3)||
 	 (outb==PDT::Colour0 && outa==PDT::Colour3))
 	cfv[0] = make_pair(0,1.);
       else if((outa==PDT::Colour8 && outb==PDT::Colour3)||
 	      (outb==PDT::Colour8 && outa==PDT::Colour3)) {
-	cfv[0] = make_pair(1,1.);
+	double sign = diag.intermediate->iSpin()==PDT::Spin0 ? -1. : 1.;
+	cfv[0] = make_pair(1,sign);
       }
     }
     else if((ina==PDT::Colour3    && inb==PDT::Colour8) ||
 	    (ina==PDT::Colour3bar && inb==PDT::Colour8) ||
 	    (inb==PDT::Colour3    && ina==PDT::Colour8) ||
 	    (inb==PDT::Colour3bar && ina==PDT::Colour8) ) {
       if((outa==PDT::Colour3    && outb==PDT::Colour3    ) ||
 	 (outa==PDT::Colour3bar && outb==PDT::Colour3bar)) {
 	cfv[0] = make_pair(1,1.);
       }
     }
   }
   else if(diag.intermediate->iColour() == PDT::Colour6 ||
 	  diag.intermediate->iColour() == PDT::Colour6bar) {
     if(ina==PDT::Colour8 && inb==PDT::Colour8) {
       cfv[0] = make_pair(0, 1.);
       for(unsigned int ix=1;ix<4;++ix)
 	cfv.push_back(make_pair(ix,1.));
       for(unsigned int ix=4;ix<8;++ix)
 	cfv.push_back(make_pair(ix,1.));
     }
     else if(outa==PDT::Colour3bar && outb==PDT::Colour6) {
       cfv[0] = make_pair(0,1.);
       for(unsigned int ix=1;ix<4;++ix)
 	cfv.push_back(make_pair(ix,1.));
     }
     else if(outa==PDT::Colour6 && outb==PDT::Colour3bar) {
       cfv[0] = make_pair(4,1.);
       cfv.push_back(make_pair(5,1.));
     }
   }
   diag.colourFlow = cfv;
 }
  
 void HardProcessConstructor::uChannelCF(HPDiagram & diag) {
   tcPDPtr ia = getParticleData(diag.incoming.first );
   tcPDPtr ib = getParticleData(diag.incoming.second);
   tcPDPtr oa = getParticleData(diag.outgoing.first );
   tcPDPtr ob = getParticleData(diag.outgoing.second);
   PDT::Colour ina  = ia->iColour();
   PDT::Colour inb  = ib->iColour();
   PDT::Colour outa = oa->iColour();
   PDT::Colour outb = ob->iColour();
   PDT::Colour offshell = diag.intermediate->iColour();
   vector<CFPair> cfv(1, make_pair(1, 1.));
   if(offshell == PDT::Colour8) {
     if(outa == PDT::Colour0 &&
        outb == PDT::Colour0) {
       cfv[0].first = 0;
     }
     else if( outa != outb ) {
       if(outa == PDT::Colour0 || 
 	 outb == PDT::Colour0) {
 	cfv[0].first = 0;
       }
       else if(ina  == PDT::Colour3 && inb  == PDT::Colour8 &&
 	      outb == PDT::Colour3 && outa == PDT::Colour8) {
 	tPDPtr off = diag.intermediate;
 	if(off->CC()) off=off->CC();
 	if(off->iSpin()!=PDT::Spin1Half ||
 	   diag.vertices.second->allowed(off->id(),diag.outgoing.first,diag.incoming.second)) {
 	  cfv[0].first = 0;
 	  cfv.push_back(make_pair(1, -1.));
 	}
 	else {
 	  cfv[0].first = 1;
 	  cfv.push_back(make_pair(0, -1.));
 	}
       }
       else if(ina  == PDT::Colour3bar && inb  == PDT::Colour8 &&
 	      outb == PDT::Colour3bar && outa == PDT::Colour8) {
 	tPDPtr off = diag.intermediate;
 	if(off->CC()) off=off->CC();
 	if(off->iSpin()!=PDT::Spin1Half ||
 	   diag.vertices.second->allowed(diag.outgoing.first,off->id(),diag.incoming.second)) {
 	  cfv[0].first = 0;
 	  cfv.push_back(make_pair(1, -1.));
 	}
 	else {
 	  cfv[0].first = 1;
 	  cfv.push_back(make_pair(0, -1.));
 	}
       }
       else {
 	cfv[0].first = 0;
 	cfv.push_back(make_pair(1, -1.));
       }
     }
     else if(outa==PDT::Colour8&&ina==PDT::Colour8) {
       cfv[0]=make_pair(4, 2.);
       cfv.push_back(make_pair(5, -2.));
       cfv.push_back(make_pair(0, -2.));
       cfv.push_back(make_pair(2,  2.));
     }
   }
   else if(offshell == PDT::Colour3 || offshell == PDT::Colour3bar) {
     if( outa == PDT::Colour0 || outb == PDT::Colour0 ) {
       if( outa == PDT::Colour6    || outb == PDT::Colour6   ||
 	  outa == PDT::Colour6bar || outb == PDT::Colour6bar) {
 	cfv[0] = make_pair(0,0.5);
 	cfv.push_back(make_pair(1,0.5));
       }
       else if ((ina==PDT::Colour3 && inb == PDT::Colour3 &&
 		(outa == PDT::Colour3bar || outb == PDT::Colour3bar))||
 	       (ina==PDT::Colour3bar && inb == PDT::Colour3bar &&
 		(outa == PDT::Colour3 || outb == PDT::Colour3 ))) {
-	double sign = diag.intermediate->iSpin()==PDT::Spin0 ? -1 : 1.;
-	cfv[0] = make_pair(0,sign);
+	cfv[0] = make_pair(0,1.);
       }
       else {
 	cfv[0] = make_pair(0,1.);
       }
     }
     else if(outa==PDT::Colour3bar && outb==PDT::Colour6) {
       cfv[0] = make_pair(4,1.);
       cfv.push_back(make_pair(5,1.));
     }
     else if(outa==PDT::Colour6 && outb==PDT::Colour3bar) {
       cfv[0] = make_pair(0,1.);
       for(int ix=0; ix<4;++ix)
 	cfv.push_back(make_pair(ix,1.));
     }
     else if(outa==PDT::Colour6bar && outb==PDT::Colour6) {
       cfv[0] = make_pair(4,1.);
       for(int ix=5; ix<8;++ix)
 	cfv.push_back(make_pair(ix,1.));
     }
     else if(ina==PDT::Colour0 && inb==PDT::Colour0) {
       cfv[0] = make_pair(0,1.);
     }
     else if(ina==PDT::Colour3    && inb==PDT::Colour3    ) {
       if((outa==PDT::Colour0 && outb==PDT::Colour3bar)||
 	 (outb==PDT::Colour0 && outa==PDT::Colour3bar))
 	cfv[0] = make_pair(0,1.);
       else if((outa==PDT::Colour8 && outb==PDT::Colour3bar)||
 	      (outb==PDT::Colour8 && outa==PDT::Colour3bar)) {
 	double sign = diag.intermediate->iSpin()==PDT::Spin0 ? -1. : 1.;
 	cfv[0] = make_pair(2,sign);
       }
     }
     else if(ina==PDT::Colour3bar && inb==PDT::Colour3bar ) {
       if((outa==PDT::Colour0 && outb==PDT::Colour3)||
 	 (outb==PDT::Colour0 && outa==PDT::Colour3))
 	cfv[0] = make_pair(0,1.);
       else if((outa==PDT::Colour8 && outb==PDT::Colour3)||
-	      (outb==PDT::Colour8 && outa==PDT::Colour3))
+	      (outb==PDT::Colour8 && outa==PDT::Colour3)) {
 	cfv[0] = make_pair(2,1.);
+      }
     }
     else if(((ina==PDT::Colour3    && inb==PDT::Colour8) ||
 	     (ina==PDT::Colour3bar && inb==PDT::Colour8) ||
 	     (inb==PDT::Colour3    && ina==PDT::Colour8) ||
 	     (inb==PDT::Colour3bar && ina==PDT::Colour8)) &&
 	    ((outa==PDT::Colour3    && outb==PDT::Colour3    ) ||
 	     (outa==PDT::Colour3bar && outb==PDT::Colour3bar))) {
       cfv[0] = make_pair(2, 1.);
     }
     else if(( ina==PDT::Colour3    &&  inb==PDT::Colour3bar && 
 	      outa==PDT::Colour3    && outb==PDT::Colour3bar)) {
       cfv[0] = make_pair(2, 1.);
       cfv.push_back(make_pair(3,-1.));
     }
   }
   else if( offshell == PDT::Colour0 ) {
     if(ina==PDT::Colour0) {
       cfv[0] = make_pair(0, 1);
     }
     else if(ina==PDT::Colour3 || ina==PDT::Colour3bar) {
       if( inb == PDT::Colour0 ) {
 	cfv[0] = make_pair(0, 1);
       }
       else if(inb==PDT::Colour3 || inb==PDT::Colour3bar) {
 	cfv[0] = make_pair(3, 1);
       }
       else if(inb==PDT::Colour8) {
 	cfv[0] = make_pair(2, 1);
       }
     }
     else if(ina==PDT::Colour8) {
       if( inb == PDT::Colour0 ) {
 	cfv[0] = make_pair(0, 1);
       }
       else if(inb==PDT::Colour3 || outb==PDT::Colour3bar) {
 	cfv[0] = make_pair(2, 1);
       }
       else if(inb==PDT::Colour8) {
 	cfv[0] = make_pair(8, -1);
       }
     }
   }
   else if(diag.intermediate->iColour() == PDT::Colour6 ||
 	  diag.intermediate->iColour() == PDT::Colour6bar) {
     if(ina==PDT::Colour8 && inb==PDT::Colour8) {
       cfv[0] = make_pair(0, 1.);
       for(unsigned int ix=1;ix<4;++ix)
 	cfv.push_back(make_pair(ix,1.));
       for(unsigned int ix=8;ix<12;++ix)
 	cfv.push_back(make_pair(ix,1.));
     }
     else if(outa==PDT::Colour3bar && outa==PDT::Colour6) {
       cfv[0] = make_pair(4, 1.);
       cfv.push_back(make_pair(5,1.));
     }
     else if(outa==PDT::Colour6 && outa==PDT::Colour3bar) {
       cfv[0] = make_pair(0, 1.);
       for(unsigned int ix=1;ix<4;++ix)
 	cfv.push_back(make_pair(ix,1.));
     }
   }
   diag.colourFlow = cfv;
 }
 
 void HardProcessConstructor::sChannelCF(HPDiagram & diag) {
   tcPDPtr pa = getParticleData(diag.incoming.first);
   tcPDPtr pb = getParticleData(diag.incoming.second);
   PDT::Colour ina = pa->iColour();
   PDT::Colour inb = pb->iColour();
   PDT::Colour offshell = diag.intermediate->iColour();
   tcPDPtr pc = getParticleData(diag.outgoing.first);
   tcPDPtr pd = getParticleData(diag.outgoing.second);
   PDT::Colour outa = pc->iColour();
   PDT::Colour outb = pd->iColour();
   vector<CFPair> cfv(1);
   if(offshell == PDT::Colour8) {
     if(ina  == PDT::Colour0 || inb  == PDT::Colour0 || 
        outa == PDT::Colour0 || outb == PDT::Colour0) {
       cfv[0] = make_pair(0, 1);
     }
     else {
       bool incol   = ina  == PDT::Colour8 && inb  == PDT::Colour8;
       bool outcol  = outa == PDT::Colour8 && outb == PDT::Colour8;
       bool intrip  = ina  == PDT::Colour3 && inb  == PDT::Colour3bar;
       bool outtrip = outa == PDT::Colour3 && outb == PDT::Colour3bar;
       bool outsex  = outa == PDT::Colour6 && outb == PDT::Colour6bar;
       bool outsexb = outa == PDT::Colour6bar && outb == PDT::Colour6;
       if(incol || outcol) {
 	// Require an additional minus sign for a scalar/fermion
 	// 33bar final state due to the way the vertex rules are defined.
 	int prefact(1);
 	if( ((pc->iSpin() == PDT::Spin1Half && pd->iSpin() == PDT::Spin1Half) ||
 	     (pc->iSpin() == PDT::Spin0     && pd->iSpin() == PDT::Spin0    )) &&
 	    (outa        == PDT::Colour3   && outb        == PDT::Colour3bar) )
 	  prefact = -1;
 	if(incol && outcol) {
 	  cfv[0] = make_pair(0, -2.);
 	  cfv.push_back(make_pair(1,  2.));
 	  cfv.push_back(make_pair(3,  2.));
 	  cfv.push_back(make_pair(5,  -2.));
 	}
 	else if(incol && outsex) {
 	  cfv[0].first = 4;
 	  cfv[0].second =  prefact;
 	  for(unsigned int ix=1;ix<4;++ix)
 	    cfv.push_back(make_pair(4+ix, prefact));
 	  for(unsigned int ix=0;ix<4;++ix)
 	    cfv.push_back(make_pair(8+ix,-prefact));
 	}
 	else {
 	  cfv[0].first = 0;
 	  cfv[0].second = -prefact;
 	  cfv.push_back(make_pair(1, prefact));
 	}
       }
       else if( (  intrip && !outtrip ) || 
 	       ( !intrip &&  outtrip ) ) {
 	if(!outsex)
 	  cfv[0] = make_pair(0, 1);
 	else {
 	  cfv[0] = make_pair(0, 1.);
 	  for(unsigned int ix=0;ix<3;++ix)
 	    cfv.push_back(make_pair(ix+1, 1.));
 	}
       }
       else if((intrip && outsex) || (intrip && outsexb)) {
 	cfv[0] = make_pair(0,1.);
 	for(int ix=1; ix<4; ++ix)
 	  cfv.push_back(make_pair(ix,1.));
       }
       else
 	cfv[0] = make_pair(1, 1);
     }
   }
   else if(offshell == PDT::Colour0) {
     if( ina == PDT::Colour0 ) {
       cfv[0] = make_pair(0, 1);
     }
     else if(ina==PDT::Colour3 || ina==PDT::Colour3bar) {
       if( outa == PDT::Colour0 ) {
 	cfv[0] = make_pair(0, 1);
       }
       else if(outa==PDT::Colour3 || outa==PDT::Colour3bar) {
 	cfv[0] = make_pair(3, 1);
       }
       else if(outa==PDT::Colour8) {
 	cfv[0] = make_pair(2, 1);
       }
       else if(outa==PDT::Colour6 || outa==PDT::Colour6bar) {
 	cfv[0] = make_pair(8, 1.);
 	cfv.push_back(make_pair(9,1.));
       }
       else
 	assert(false);
     }
     else if(ina==PDT::Colour8) {
       if( outa == PDT::Colour0 ) {
 	cfv[0] = make_pair(0, 1);
       }
       else if(outa==PDT::Colour3 || outb==PDT::Colour3bar) {
 	cfv[0] = make_pair(2, 1);
       }
       else if(outa==PDT::Colour8) {
 	cfv[0] = make_pair(6, 1);
       }
     }
   }
   else if(offshell == PDT::Colour3 || offshell == PDT::Colour3bar) {
     if(outa == PDT::Colour6    || outa == PDT::Colour6bar || 
        outb == PDT::Colour6bar || outb == PDT::Colour6) {
       cfv[0] = make_pair(6, 1.);
       cfv.push_back(make_pair(7,1.));
     }
     else if((ina  == PDT::Colour3    && inb  == PDT::Colour3) ||
 	    (ina  == PDT::Colour3bar && inb  == PDT::Colour3bar)) {
       if((outa == PDT::Colour3    && outb == PDT::Colour3   ) ||
 	 (outa == PDT::Colour3bar && outb == PDT::Colour3bar)) {
 	cfv[0]      = make_pair(2, 1.);
 	cfv.push_back(make_pair(3,-1.));
       }
       else
 	cfv[0] = make_pair(0,1.);
     }
     else if(((ina==PDT::Colour3    && inb==PDT::Colour8) ||
 	     (ina==PDT::Colour3bar && inb==PDT::Colour8) ||
 	     (inb==PDT::Colour3    && ina==PDT::Colour8) ||
 	     (inb==PDT::Colour3bar && ina==PDT::Colour8) ) &&
 	    ((outa==PDT::Colour3    && outb==PDT::Colour3    ) ||
 	     (outa==PDT::Colour3bar && outb==PDT::Colour3bar))) {
       cfv[0] = make_pair(0,1.);
     }
     else {
       if(outa == PDT::Colour0 || outb == PDT::Colour0)
 	cfv[0] = make_pair(0, 1);
       else
 	cfv[0] = make_pair(1, 1);
     }
   }
   else if( offshell == PDT::Colour6 || offshell == PDT::Colour6bar) {
     if((ina  == PDT::Colour3    && inb  == PDT::Colour3    &&
 	outa == PDT::Colour3    && outb == PDT::Colour3   ) ||
        (ina  == PDT::Colour3bar && inb  == PDT::Colour3bar &&
 	outa == PDT::Colour3bar && outb == PDT::Colour3bar)) {
       cfv[0]      = make_pair(2,0.5);
       cfv.push_back(make_pair(3,0.5));
     }
     else if((ina  == PDT::Colour3    && inb  == PDT::Colour3    &&
 	     ((outa == PDT::Colour6    && outb == PDT::Colour0)||
 	      (outb == PDT::Colour6    && outa == PDT::Colour0))) ||
 	    (ina  == PDT::Colour3bar && inb  == PDT::Colour3bar &&
 	     ((outa == PDT::Colour6bar && outb == PDT::Colour0)||
 	      (outb == PDT::Colour6bar && outa == PDT::Colour0)))) {
       cfv[0]      = make_pair(0,0.5);
       cfv.push_back(make_pair(1,0.5));
     }
     else
       assert(false);
   }
   else {
     if(outa == PDT::Colour0 || outb == PDT::Colour0)
       cfv[0] = make_pair(0, 1);
     else
       cfv[0] = make_pair(1, 1);
   }  
   diag.colourFlow = cfv; 
 }
 
 void HardProcessConstructor::fourPointCF(HPDiagram & diag) {
   using namespace ThePEG::Helicity;
   // count the colours
   unsigned int noct(0),ntri(0),nsng(0),nsex(0),nf(0);
   vector<tcPDPtr> particles;
   for(unsigned int ix=0;ix<4;++ix) {
     particles.push_back(getParticleData(diag.ids[ix]));
     PDT::Colour col = particles.back()->iColour();
     if(col==PDT::Colour0)                            ++nsng;
     else if(col==PDT::Colour3||col==PDT::Colour3bar) ++ntri;
     else if(col==PDT::Colour8)                       ++noct;
     else if(col==PDT::Colour6||col==PDT::Colour6bar) ++nsex;
     if(particles.back()->iSpin()==2) nf+=1;
   }
   if(nsng==4 || (ntri==2&&nsng==2) || 
      (noct==3            && nsng==1) ||
      (ntri==2 && noct==1 && nsng==1) ||
      (noct == 2 && nsng == 2) ) {
     vector<CFPair> cfv(1,make_pair(0,1));
     diag.colourFlow = cfv;
   }
   else if(noct==4) {
     // flows for SSVV, VVVV is handled in me class
     vector<CFPair> cfv(6);
     cfv[0] = make_pair(0, -2.);
     cfv[1] = make_pair(1, -2.);
     cfv[2] = make_pair(2, +4.);
     cfv[3] = make_pair(3, -2.);
     cfv[4] = make_pair(4, +4.);
     cfv[5] = make_pair(5, -2.);
     diag.colourFlow = cfv;
   }
   else if(ntri==2&&noct==2) {
     vector<CFPair> cfv(2);
     cfv[0] = make_pair(0, 1);
     cfv[1] = make_pair(1, 1);
     if(nf==2) cfv[1].second = -1.;
     diag.colourFlow = cfv;
   }
   else if(nsex==2&&noct==2) {
     vector<CFPair> cfv;
     for(unsigned int ix=0;ix<4;++ix)
       cfv.push_back(make_pair(ix  ,2.));
     for(unsigned int ix=0;ix<8;++ix)
       cfv.push_back(make_pair(4+ix,1.));
     diag.colourFlow = cfv;
   }
   else if(ntri==4) {
     // get the order from the vertex
     vector<long> temp;
     for(unsigned int ix=0;ix<4;++ix) {
       temp = diag.vertices.first->search(ix,diag.outgoing.first);
       if(!temp.empty()) break;
     }
     // compute the mapping
     vector<long> ids;
     ids.push_back( particles[0]->CC() ? -diag.incoming.first  : diag.incoming.first );
     ids.push_back( particles[1]->CC() ? -diag.incoming.second : diag.incoming.second);
     ids.push_back(  diag.outgoing.first );
     ids.push_back(  diag.outgoing.second);
     vector<unsigned int> order = {0,1,2,3};
     vector<bool> matched(4,false);
     for(unsigned int ix=0;ix<temp.size();++ix) {
       for(unsigned int iy=0;iy<ids.size();++iy) {
 	if(matched[iy]) continue;
 	if(temp[ix]==ids[iy]) {
 	  matched[iy] = true;
 	  order[ix]=iy;
 	  break;
 	}
       }
     }
     // 3 3 -> 3 3
     if((particles[0]->iColour()==PDT::Colour3 &&
 	particles[1]->iColour()==PDT::Colour3) ||
        (particles[0]->iColour()==PDT::Colour3bar &&
 	particles[1]->iColour()==PDT::Colour3bar) ) {
       if(diag.vertices.first->colourStructure()==ColourStructure::SU3I12I34) {
 	if( (order[0]==0 && order[1]==2) || (order[2]==0 && order[3]==2) ||
 	    (order[0]==2 && order[1]==0) || (order[2]==2 && order[3]==0))
       	  diag.colourFlow = vector<CFPair>(1,make_pair(2,1.));
       	else
       	  diag.colourFlow = vector<CFPair>(1,make_pair(3,1.));
       }
       else if(diag.vertices.first->colourStructure()==ColourStructure::SU3I14I23) {
       	if( (order[0]==0 && order[3]==2) || (order[1]==0 && order[2]==2) ||
       	    (order[0]==2 && order[3]==0) || (order[1]==2 && order[2]==0))
       	  diag.colourFlow = vector<CFPair>(1,make_pair(2,1.));
       	else
       	  diag.colourFlow = vector<CFPair>(1,make_pair(3,1.));
       }
       else if(diag.vertices.first->colourStructure()==ColourStructure::SU3T21T43) {
 	if( (order[1]==0 && order[0]==2) || (order[3]==0 && order[2]==2) ||
 	    (order[1]==2 && order[0]==0) || (order[3]==2 && order[2]==0))
 	  diag.colourFlow = vector<CFPair>(1,make_pair(0,1.));
 	else
 	  diag.colourFlow = vector<CFPair>(1,make_pair(1,1.));
       }
       else if(diag.vertices.first->colourStructure()==ColourStructure::SU3T23T41) {
 	if( (order[1]==0 && order[2]==2) || (order[3]==0 && order[0]==2) ||
 	    (order[1]==2 && order[2]==0) || (order[3]==2 && order[0]==0))
 	  diag.colourFlow = vector<CFPair>(1,make_pair(0,1.));
 	else
 	  diag.colourFlow = vector<CFPair>(1,make_pair(1,1.));
       }
       else
 	assert(false);
     }
     else if((particles[0]->iColour()==PDT::Colour3 &&
 	     particles[1]->iColour()==PDT::Colour3bar) |
 	    (particles[0]->iColour()==PDT::Colour3bar &&
 	     particles[1]->iColour()==PDT::Colour3)) {
       if(diag.vertices.first->colourStructure()==ColourStructure::SU3I12I34) {
       	if( (order[0]==0 && order[1]==1) || (order[2]==0 && order[3]==0) ||
       	    (order[0]==1 && order[1]==0) || (order[2]==1 && order[3]==1))
       	  diag.colourFlow = vector<CFPair>(1,make_pair(3,1.));
       	else
       	  diag.colourFlow = vector<CFPair>(1,make_pair(2,1.));
       }
       else if(diag.vertices.first->colourStructure()==ColourStructure::SU3I14I23) {
       	if( (order[0]==0 && order[3]==1) || (order[0]==2 && order[3]==3) ||
       	    (order[0]==1 && order[3]==0) || (order[0]==3 && order[3]==2))
       	  diag.colourFlow = vector<CFPair>(1,make_pair(3,1.));
       	else
       	  diag.colourFlow = vector<CFPair>(1,make_pair(2,1.));
       }
       else if(diag.vertices.first->colourStructure()==ColourStructure::SU3T21T43) {
        	if( (order[1]==0 && order[0]==1) || (order[3]==0 && order[2]==1) ||
       	    (order[1]==1 && order[0]==0) || (order[3]==1 && order[2]==0))
        	  diag.colourFlow = vector<CFPair>(1,make_pair(1,1.));
       	else
       	  diag.colourFlow = vector<CFPair>(1,make_pair(0,1.));
       }
       else if(diag.vertices.first->colourStructure()==ColourStructure::SU3T23T41) {
 	if( (order[1]==0 && order[2]==1) || (order[1]==1 && order[2]==0) ||
 	    (order[1]==3 && order[2]==2) || (order[1]==2 && order[2]==3))
 	  diag.colourFlow = vector<CFPair>(1,make_pair(1,1.));
 	else
 	  diag.colourFlow = vector<CFPair>(1,make_pair(2,1.));
       }
       else
 	assert(false);
     }
     else {
       assert(false);
     }
   }
   else {
     assert(false);
   }
 }
 
 namespace {
   // Helper functor for find_if in duplicate function.
   class SameDiagramAs {
   public:
     SameDiagramAs(const HPDiagram & diag) : a(diag) {}
     bool operator()(const HPDiagram & b) const {
       return a == b;
     }
   private:
     HPDiagram a;
   };
 }
 
 bool HardProcessConstructor::duplicate(const HPDiagram & diag, 
 				       const HPDVector & group) const {
   //find if a duplicate diagram exists
   HPDVector::const_iterator it = 
     find_if(group.begin(), group.end(), SameDiagramAs(diag));
   return it != group.end();
 } 
 
 bool HardProcessConstructor::checkOrder(const HPDiagram & diag) const {
   for(map<string,pair<unsigned int,int> >::const_iterator it=model_->couplings().begin();
       it!=model_->couplings().end();++it) {
     int order=0;
     if(diag.vertices.first ) order += diag.vertices.first ->orderInCoupling(it->second.first);
     if(diag.vertices.second&&diag.vertices.first->getNpoint()==3)
       order += diag.vertices.second->orderInCoupling(it->second.first);
     if(order>it->second.second) return false;
   }
   return true;
 }