diff --git a/Decay/DecayIntegrator.cc b/Decay/DecayIntegrator.cc
--- a/Decay/DecayIntegrator.cc
+++ b/Decay/DecayIntegrator.cc
@@ -1,381 +1,382 @@
 // -*- C++ -*-
 //
 // DecayIntegrator.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 DecayIntegrator class.
 //
 // Author: Peter Richardson
 // 
 
 #include "DecayIntegrator.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/PDT/DecayMode.h"
 #include "ThePEG/PDT/EnumParticles.h"
 #include "ThePEG/EventRecord/Particle.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Interface/ParVector.h"
 #include "ThePEG/Interface/RefVector.h"
 #include "ThePEG/Utilities/Debug.h"
 #include "Herwig/Utilities/Kinematics.h"
 #include "Herwig/PDT/GenericMassGenerator.h"
 #include "DecayPhaseSpaceMode.h"
 #include "Herwig/PDT/WidthCalculatorBase.h"
 #include "ThePEG/Interface/Reference.h"
 
 using namespace Herwig; 
 
 ParticleVector DecayIntegrator::decay(const Particle & parent,
 				      const tPDVector & children) const {
   // return empty vector if products heavier than parent
   Energy mout(ZERO);
   for(tPDVector::const_iterator it=children.begin();
       it!=children.end();++it) mout+=(**it).massMin();
   if(mout>parent.mass()) return ParticleVector();
   // generate the decay
   bool cc;
   _imode = modeNumber(cc,parent.dataPtr(),children);
+  if(numberModes()==0) return ParticleVector();
   return _modes[_imode]->generate(_generateinter,cc,parent);
 }
   
 void DecayIntegrator::persistentOutput(PersistentOStream & os) const {
   os << _modes << _niter << _npoint << _ntry << _photongen << _generateinter << ounit(_eps,GeV);
 }
   
 void DecayIntegrator::persistentInput(PersistentIStream & is, int) {
   is >> _modes >> _niter >> _npoint >> _ntry >> _photongen >> _generateinter >> iunit(_eps,GeV);
 }
   
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeAbstractClass<DecayIntegrator,HwDecayerBase>
 describeHerwigDecayIntegrator("Herwig::DecayIntegrator", "Herwig.so");
   
 void DecayIntegrator::Init() {
   
   static ClassDocumentation<DecayIntegrator> documentation
     ("The DecayIntegrator class is a base decayer class "
      "including a multi-channel integrator.");
   
   static Parameter<DecayIntegrator,int> interfaceIteration
     ("Iteration",
      "Number of iterations for the initialization of the phase space",
      &DecayIntegrator::_niter, 10, 0, 100,
      false, false, true);  
   
   static Parameter<DecayIntegrator,int> interfacePoints
     ("Points",
      "number of phase space points to generate in the initialisation.",
      &DecayIntegrator::_npoint, 10000, 1, 1000000000,
      false, false, true);
   
   static Parameter<DecayIntegrator,int> interfaceNtry
     ("Ntry",
      "Number of attempts to generate the decay",
      &DecayIntegrator::_ntry, 500, 0, 100000,
      false, false, true);
 
   static Reference<DecayIntegrator,DecayRadiationGenerator> interfacePhotonGenerator
     ("PhotonGenerator",
      "Object responsible for generating photons in the decay.",
      &DecayIntegrator::_photongen, false, false, true, true, false);
  
   static Switch<DecayIntegrator,bool> interfaceGenerateIntermediates
     ("GenerateIntermediates",
      "Whether or not to include intermediate particles in the output",
      &DecayIntegrator::_generateinter, false, false, false);
   static SwitchOption interfaceGenerateIntermediatesNoIntermediates
     (interfaceGenerateIntermediates,
      "No",
      "Don't include the intermediates",
      false);
   static SwitchOption interfaceGenerateIntermediatesIncludeIntermediates
     (interfaceGenerateIntermediates,
      "Yes",
      "include the intermediates",
      true);
 }
 
 // output info on the integrator
 ostream & Herwig::operator<<(ostream & os, const DecayIntegrator & decay) {
   os << "The integrator has " << decay._modes.size() << " modes"  << endl;
   for(unsigned int ix=0;ix<decay._modes.size();++ix) {
     os << "Information on mode " << ix << endl;
     os << *(decay._modes[ix]);
   }
   return os;
 }
 
 // generate the momenta for the decay
 ParticleVector DecayIntegrator::generate(bool inter,bool cc,
 					 const unsigned int & imode,
 					 const Particle & inpart) const {
   _imode=imode;
   return _modes[imode]->generate(inter,cc,inpart);
 }  
 
 // initialization for a run
 void DecayIntegrator::doinitrun() {
   HwDecayerBase::doinitrun();
   if ( initialize() && Debug::level > 1 ) 
     CurrentGenerator::current().log() << "Start of the initialisation for " 
 				      << this->name() << "\n";
   for(unsigned int ix=0;ix<_modes.size();++ix) {
     if(!_modes[ix]) continue;
     _modes[ix]->initrun();
     _imode=ix;
     _modes[ix]->initializePhaseSpace(initialize());
   }
   // CurrentGenerator::current().log() << *this << "\n";
 }
 
 // add a new mode
 void DecayIntegrator::addMode(DecayPhaseSpaceModePtr in,double maxwgt,
 				     const vector<double> inwgt) const {
   _modes.push_back(in);
   if(in) {
     in->setMaxWeight(maxwgt);
     in->setWeights(inwgt);
     in->setIntegrate(_niter,_npoint,_ntry);
     in->init();
   }
 }
 
 // reset the properities of all intermediates
 void DecayIntegrator::resetIntermediate(tcPDPtr part, Energy mass, Energy width) {
   if(!part) return;
   for(unsigned int ix=0,N=_modes.size();ix<N;++ix) {
     _modes[ix]->resetIntermediate(part,mass,width);
   }
 } 
 
 WidthCalculatorBasePtr 
 DecayIntegrator::threeBodyMEIntegrator(const DecayMode &) const {
   return WidthCalculatorBasePtr();
 }
 
 // set the code for the partial width
 void DecayIntegrator::setPartialWidth(const DecayMode & dm, int imode) {
   vector<int> extid;
   tcPDPtr cc,cc2;
   int nfound(0),ifound,nmax(1),id;
   unsigned int ix(0),iy,N,iz,tmax,nmatched;
   if(dm.parent()->CC()) nmax=2;
   if(_modes.size()==0) return;
   do {
     if(!_modes[ix]) {
       ++ix;
       continue;
     }
     cc = _modes[ix]->externalParticles(0)->CC();
     tmax = cc ? 1 : 2;
     for(iz=0;iz<tmax;++iz) {
       ifound=-1;
       extid.clear();
       // check the parent
       if(dm.parent()->id()==_modes[ix]->externalParticles(0)->id()&&iz==0) {
 	for(iy=0,N=_modes[ix]->numberofParticles();iy<N;++iy) {
 	  extid.push_back(_modes[ix]->externalParticles(iy)->id());
 	}
       }
       else if(dm.parent()->id()==_modes[ix]->externalParticles(0)->id()&&iz==1) {
 	for(iy=0,N=_modes[ix]->numberofParticles();iy<N;++iy) {
 	  cc2=_modes[ix]->externalParticles(iy)->CC();
 	  extid.push_back( cc2 ? cc2->id() : _modes[ix]->externalParticles(iy)->id());
 	}
       }
       else if(cc&&dm.parent()->id()==cc->id()) {
 	for(iy=0,N=_modes[ix]->numberofParticles();iy<N;++iy) {
 	  cc2 = _modes[ix]->externalParticles(iy)->CC();
 	  extid.push_back( cc2 ? cc2->id() : _modes[ix]->externalParticles(iy)->id());
 	}
       }
       // if the parents match
       if(!extid.empty()) {
 	vector<bool> matched(extid.size(),false);
 	bool done;
 	nmatched=0;
 	ParticleMSet::const_iterator pit = dm.products().begin();
 	do {
 	  id=(**pit).id();
 	  done=false;
 	  iy=1;
 	  do {
 	    if(id==extid[iy]&&!matched[iy]) {
 	      matched[iy]=true;
 	      ++nmatched;
 	      done=true;
 	    }
 	    ++iy;
 	  }
 	  while(iy<extid.size()&&!done);
 	  ++pit;
 	}
 	while(pit!=dm.products().end());
 	if(nmatched==extid.size()-1) {
 	  ifound=ix;
 	  ++nfound;
 	}
       }
       if(ifound>=0) _modes[ifound]->setPartialWidth(imode);
     }
     ++ix;
   }
   while(nfound<nmax&&ix<_modes.size());
 }
 
 int DecayIntegrator::findMode(const DecayMode & dm) {
   int imode(-1);
   vector<int> extid;
   tcPDPtr cc,cc2;
   bool found(false);
   int id;
   unsigned int ix(0),iy,N,iz,tmax,nmatched;
   if(_modes.size()==0) return -1;
   do {
     if(!_modes[ix]) {
       ++ix;
       continue;
     }
     cc = _modes[ix]->externalParticles(0)->CC();
     tmax=1;if(!cc){++tmax;}
     for(iz=0;iz<tmax;++iz) {
       extid.clear();
       // check the parent
       if(dm.parent()->id()==_modes[ix]->externalParticles(0)->id()&&iz==0) {
 	for(iy=0,N=_modes[ix]->numberofParticles();iy<N;++iy) {
 	  extid.push_back(_modes[ix]->externalParticles(iy)->id());
 	}
       }
       else if(dm.parent()->id()==_modes[ix]->externalParticles(0)->id()&&iz==1) {
 	for(iy=0,N=_modes[ix]->numberofParticles();iy<N;++iy) {
 	  cc2=_modes[ix]->externalParticles(iy)->CC();
 	  extid.push_back( cc2 ? cc2->id() : _modes[ix]->externalParticles(iy)->id());
 	}
       }
       else if(cc&&dm.parent()->id()==cc->id()) {
 	for(iy=0,N=_modes[ix]->numberofParticles();iy<N;++iy) {
 	  cc = _modes[ix]->externalParticles(iy)->CC();
 	  extid.push_back( cc ? cc->id() : _modes[ix]->externalParticles(iy)->id());
 	}
       }
       // if the parents match
       if(!extid.empty()) {
 	vector<bool> matched(extid.size(),false);
 	bool done;
 	nmatched=0;
 	ParticleMSet::const_iterator pit = dm.products().begin();
 	do {
 	  id=(**pit).id();
 	  done=false;
 	  iy=1;
 	  do {
 	    if(id==extid[iy]&&!matched[iy]) {
 	      matched[iy]=true;
 	      ++nmatched;
 	      done=true;
 	    }
 	    ++iy;
 	  }
 	  while(iy<extid.size()&&!done);
 	  ++pit;
 	}
 	while(pit!=dm.products().end());
 	if(nmatched==extid.size()-1) {
 	  imode=ix;
 	  found=true;
 	}
       }
     }
     ++ix;
   }
   while(!found&&ix<_modes.size());
   return imode;
 }
 
 // output the information for the database
 void DecayIntegrator::dataBaseOutput(ofstream & output,bool header) const {
   // header for MySQL
   if(header) output << "update decayers set parameters=\"";
   output << "newdef " << name() << ":Iteration " << _niter << "\n";
   output << "newdef " << name() << ":Ntry " << _ntry << "\n";
   output << "newdef " << name() << ":Points " << _npoint << "\n";
   //if(_photongen){;}
   output << "newdef " << name() << ":GenerateIntermediates " << _generateinter << " \n";
   // footer for MySQL
   if(header) {
     output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";\n";
   }
 }
   
 // pointer to a mode
 tDecayPhaseSpaceModePtr DecayIntegrator::mode(unsigned int ix) {
   return _modes[ix];
 }
 
 tcDecayPhaseSpaceModePtr DecayIntegrator::mode(unsigned int ix) const {
   return _modes[ix];
 }
  
 Energy DecayIntegrator::initializePhaseSpaceMode(unsigned int imode,bool init, bool onShell) const{
   tcDecayPhaseSpaceModePtr cmodeptr=mode(imode);
   tDecayPhaseSpaceModePtr modeptr = const_ptr_cast<tDecayPhaseSpaceModePtr>(cmodeptr);
   modeptr->init();
   return modeptr->initializePhaseSpace(init,onShell);
 }
 
 // the matrix element to be integrated for the me
 double DecayIntegrator::threeBodyMatrixElement(const int,const Energy2,
 					       const Energy2,
 					       const Energy2,const Energy2,
 					       const Energy, const Energy, 
 					       const Energy) const {
   throw DecayIntegratorError() 
     << "Calling the virtual DecayIntegrator::threeBodyMatrixElement"
     << "method. This must be overwritten in the classes "
     << "inheriting from DecayIntegrator where it is needed"
     << Exception::runerror;
 }
 
 // the differential three body decay rate with one integral performed
 InvEnergy DecayIntegrator::threeBodydGammads(const int, const Energy2,
 					     const Energy2,
 					     const Energy, const Energy, 
 					     const Energy) const {
   throw DecayIntegratorError() 
     << "Calling the virtual DecayIntegrator::threeBodydGammads()" 
     <<"method. This must be overwritten in the classes "
     << "inheriting from DecayIntegrator where it is needed"
     << Exception::runerror;
 }
 
 double DecayIntegrator::oneLoopVirtualME(unsigned int ,
 					 const Particle &, 
 					 const ParticleVector &) {
   throw DecayIntegratorError()
     << "DecayIntegrator::oneLoopVirtualME() called. This should"
     << " have been overidden in an inheriting class if it is used"
     << Exception::runerror;
 }
 
 InvEnergy2 DecayIntegrator::realEmissionME(unsigned int,
 					   const Particle &, 
 					   ParticleVector &,
 					   unsigned int,
 					   double, double, 
 					   const LorentzRotation &,
 					   const LorentzRotation &) {
   throw DecayIntegratorError()
     << "DecayIntegrator::realEmmisionME() called. This should"
     << " have been overidden in an inheriting class if it is used"
     << Exception::runerror;
 }
diff --git a/Decay/General/GeneralThreeBodyDecayer.cc b/Decay/General/GeneralThreeBodyDecayer.cc
--- a/Decay/General/GeneralThreeBodyDecayer.cc
+++ b/Decay/General/GeneralThreeBodyDecayer.cc
@@ -1,1053 +1,1054 @@
 // -*- C++ -*-
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the GeneralThreeBodyDecayer class.
 //
 
 #include "GeneralThreeBodyDecayer.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "Herwig/Decay/DecayPhaseSpaceMode.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
 
 using namespace Herwig;
 
 void GeneralThreeBodyDecayer::persistentOutput(PersistentOStream & os) const {
   os << incoming_ << outgoing_ << diagrams_ << diagmap_ << colour_ << colourLargeNC_
      << nflow_ << widthOpt_ << refTag_ << refTagCC_ << intOpt_ << relerr_;
 }
 
 void GeneralThreeBodyDecayer::persistentInput(PersistentIStream & is, int) {
   is >> incoming_ >> outgoing_ >> diagrams_ >> diagmap_ >> colour_ >> colourLargeNC_
      >> nflow_ >> widthOpt_ >> refTag_ >> refTagCC_ >> intOpt_ >> relerr_;
 }
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeAbstractClass<GeneralThreeBodyDecayer,DecayIntegrator>
 describeHerwigGeneralThreeBodyDecayer("Herwig::GeneralThreeBodyDecayer", "Herwig.so");
 
 void GeneralThreeBodyDecayer::Init() {
 
   static ClassDocumentation<GeneralThreeBodyDecayer> documentation
     ("The GeneralThreeBodyDecayer class is the base class for the implementation of"
      " all three body decays based on spin structures in Herwig.");
 
   static Switch<GeneralThreeBodyDecayer,unsigned int> interfaceWidthOption
     ("WidthOption",
      "Option for the treatment of the widths of the intermediates",
      &GeneralThreeBodyDecayer::widthOpt_, 1, false, false);
   static SwitchOption interfaceWidthOptionFixed
     (interfaceWidthOption,
      "Fixed",
      "Use fixed widths",
      1);
   static SwitchOption interfaceWidthOptionRunning
     (interfaceWidthOption,
      "Running",
      "Use running widths",
      2);
   static SwitchOption interfaceWidthOptionZero
     (interfaceWidthOption,
      "Zero",
      "Set the widths to zero",
      3);
 
   static Switch<GeneralThreeBodyDecayer,unsigned int> interfacePartialWidthIntegration
     ("PartialWidthIntegration",
      "Switch to control the partial width integration",
      &GeneralThreeBodyDecayer::intOpt_, 0, false, false);
   static SwitchOption interfacePartialWidthIntegrationAllPoles
     (interfacePartialWidthIntegration,
      "AllPoles",
      "Include all potential poles",
      0);
   static SwitchOption interfacePartialWidthIntegrationShallowestPole
     (interfacePartialWidthIntegration,
      "ShallowestPole",
      "Only include the shallowest pole",
      1);
 
   static Parameter<GeneralThreeBodyDecayer,double> interfaceRelativeError
     ("RelativeError",
      "The relative error for the GQ integration of the partial width",
      &GeneralThreeBodyDecayer::relerr_, 1e-2, 1e-10, 1.,
      false, false, Interface::limited);
 
 }
 
 ParticleVector GeneralThreeBodyDecayer::decay(const Particle & parent,
 					      const tPDVector & children) const {
   // return empty vector if products heavier than parent
   Energy mout(ZERO);
   for(tPDVector::const_iterator it=children.begin();
       it!=children.end();++it) mout+=(**it).massMin();
   if(mout>parent.mass()) return ParticleVector();
   // generate the decay
   bool cc;
   int imode=modeNumber(cc,parent.dataPtr(),children);
   // generate the kinematics
   ParticleVector decay=generate(generateIntermediates(),cc,imode,parent);
   // make the colour connections
   colourConnections(parent, decay);
   // return the answer
   return decay;
 }
 
 int  GeneralThreeBodyDecayer::
 modeNumber(bool & cc, tcPDPtr in, const tPDVector & outin) const {
   assert( !refTag_.empty() && !refTagCC_.empty() );
   // check number of outgoing particles
   if( outin.size() != 3 || abs(in->id()) != abs(incoming_->id()) ) return -1;
   OrderedParticles testmode(outin.begin(), outin.end());
   OrderedParticles::const_iterator dit = testmode.begin();
   string testtag(in->name() + "->");
   for( unsigned int i = 1; dit != testmode.end(); ++dit, ++i) {
     testtag += (**dit).name();
     if( i != 3 ) testtag += string(",");
   }
   if( testtag == refTag_ ) {
     cc = false;
     return 0;
   }
   else if ( testtag == refTagCC_ ) {
     cc = true;
     return 0;
   }
   else return -1;
 }
 
 bool GeneralThreeBodyDecayer::setDecayInfo(PDPtr incoming,
 					   vector<PDPtr> outgoing,
 					   const vector<TBDiagram> & process,
 					   double symfac) {
   // set the member variables from the info supplied
   incoming_        = incoming;
   outgoing_        = outgoing;
   diagrams_        = process;
   assert( outgoing_.size() == 3 );
   // Construct reference tags for testing in modeNumber function
   OrderedParticles refmode(outgoing_.begin(), outgoing_.end());
   OrderedParticles::const_iterator dit = refmode.begin();
   refTag_ = incoming_->name() + "->";
   for( unsigned int i = 1; dit != refmode.end(); ++dit, ++i) {
     refTag_ += (**dit).name();
     if( i != 3 )  refTag_ += string(",");
   }
   //CC-mode
   refmode.clear();
   refTagCC_ = incoming_->CC() ? incoming_->CC()->name() : 
     incoming_->name();
   refTagCC_ += "->";
   for( unsigned int i = 0;  i < 3; ++i ) {
     if( outgoing_[i]->CC() ) refmode.insert( outgoing_[i]->CC() );
     else refmode.insert( outgoing_[i] );
   }
   dit = refmode.begin();
   for( unsigned int i = 1; dit != refmode.end(); ++dit , ++i) {
     refTagCC_ += (**dit).name();
     if( i != 3 ) refTagCC_ += string(",");
   }
   // check if intermeidates or only four point diagrams
   bool intermediates(false);
   for(auto diagram : diagrams_) {
     if(diagram.intermediate) {
       intermediates=true;
       break;
     }
   }
   if(!intermediates) {
     incoming_= PDPtr();
     outgoing_.clear();
     generator()->log() << "Only four body diagrams for decay "
 		       << refTag_  << " in GeneralThreeBodyDecayer::"
 		       << "setDecayInfo(), omitting decay\n";
     return false;
   }
   // set the colour factors and return the answer
   if(setColourFactors(symfac)) return true;
   incoming_= PDPtr();
   outgoing_.clear();
   return false;
 }
 
 void GeneralThreeBodyDecayer::doinit() {
   DecayIntegrator::doinit();
   if(outgoing_.empty()) return;
   setupDiagrams(false);
 }
 
 void GeneralThreeBodyDecayer::doinitrun() {
   setupDiagrams(true);
   DecayIntegrator::doinitrun();
 }
 
 double GeneralThreeBodyDecayer::
 threeBodyMatrixElement(const int imode,  const Energy2 q2,
 		       const Energy2 s3, const Energy2 s2, 
 		       const Energy2 s1, const Energy  m1, 
 		       const Energy  m2, const Energy  m3) const {
   // calculate the momenta of the outgoing particles
   Energy m0=sqrt(q2);
   // energies
   Energy eout[3] = {0.5*(q2+sqr(m1)-s1)/m0,
 		    0.5*(q2+sqr(m2)-s2)/m0,
 		    0.5*(q2+sqr(m3)-s3)/m0};
   // magnitudes of the momenta
   Energy pout[3] = {sqrt(sqr(eout[0])-sqr(m1)),
 		    sqrt(sqr(eout[1])-sqr(m2)),
 		    sqrt(sqr(eout[2])-sqr(m3))};
   double cos2 = 0.5*(sqr(pout[0])+sqr(pout[1])-sqr(pout[2]))/pout[0]/pout[1];
   double cos3 = 0.5*(sqr(pout[0])-sqr(pout[1])+sqr(pout[2]))/pout[0]/pout[2];
   double sin2 = sqrt(1.-sqr(cos2)), sin3 = sqrt(1.-sqr(cos3));
   Lorentz5Momentum out[3]=
     {Lorentz5Momentum(      ZERO   , ZERO ,  pout[0]      , eout[0] , m1),
      Lorentz5Momentum(  pout[1]*sin2 , ZERO , -pout[1]*cos2 , eout[1] , m2),
      Lorentz5Momentum( -pout[2]*sin3 , ZERO , -pout[2]*cos3 , eout[2] , m3)};
   // create the incoming
   PPtr inpart=mode(imode)->externalParticles(0)->
     produceParticle(Lorentz5Momentum(sqrt(q2)));
   // and outgoing particles
   ParticleVector decay;
   for(unsigned int ix=1;ix<4;++ix)
     decay.push_back(mode(imode)->externalParticles(ix)->produceParticle(out[ix-1]));
   // return the matrix element
   return me2(-1,*inpart,decay,Initialize);
 }
 
 double GeneralThreeBodyDecayer::brat(const DecayMode &, const Particle & p,
 				     double oldbrat) const {
   ParticleVector children = p.children();
   if( children.size() != 3 || !p.data().widthGenerator() ) 
     return oldbrat;
   // partial width for this mode
   Energy scale = p.mass();
   Energy pwidth = 
     partialWidth( make_pair(p.dataPtr(), scale),
 		  make_pair(children[0]->dataPtr(), children[0]->mass()),
 		  make_pair(children[1]->dataPtr(), children[1]->mass()),
 		  make_pair(children[2]->dataPtr(), children[2]->mass()) );
   Energy width = p.data().widthGenerator()->width(p.data(), scale);
   return pwidth/width;
 }
 
 Energy GeneralThreeBodyDecayer::partialWidth(PMPair inpart, PMPair outa, 
 					     PMPair outb, PMPair outc) const {
   if(inpart.second<outa.second+outb.second+outc.second) return ZERO;
   // create the object to calculate the width if it doesn't all ready exist
   if(!widthCalc_) {
     string tag = incoming_->name() + "->";
     tag += outgoing_[0]->name() + "," + outgoing_[1]->name() + ","
       + outgoing_[2]->name() + ";";
     DMPtr dm = generator()->findDecayMode(tag);
     widthCalc_ = threeBodyMEIntegrator(*dm);
   }
   return widthCalc_->partialWidth(sqr(inpart.second));
 }
 
 void GeneralThreeBodyDecayer::
 colourConnections(const Particle & parent,
 		  const ParticleVector & out) const {
   // first extract the outgoing particles and intermediate
   PPtr inter;
   ParticleVector outgoing;
   if(!generateIntermediates()) {
     outgoing=out;
   }
   else {
     // find the diagram
     unsigned int idiag = diagramMap()[mode(imode())->selectedChannel()];
     PPtr child;
     for(unsigned int ix=0;ix<out.size();++ix) {
       if(out[ix]->children().empty()) child = out[ix];
       else                            inter = out[ix];
     }
     outgoing.resize(3);
     switch(diagrams_[idiag].channelType) {
     case TBDiagram::channel23:
       outgoing[0] = child;
       outgoing[1] = inter->children()[0];
       outgoing[2] = inter->children()[1];
       break;
     case TBDiagram::channel13:
       outgoing[0] = inter->children()[0];
       outgoing[1] = child;
       outgoing[2] = inter->children()[1];
       break;
     case TBDiagram::channel12:
       outgoing[0] = inter->children()[0];
       outgoing[1] = inter->children()[1];
       outgoing[2] = child;
       break;
     default:
       throw Exception() << "unknown diagram type in GeneralThreeBodyDecayer::"
 			<< "colourConnections()" << Exception::runerror;
     }
   }
   // extract colour of the incoming and outgoing particles
   PDT::Colour inColour(parent.data().iColour());
   vector<PDT::Colour> outColour;
   vector<int> singlet,octet,triplet,antitriplet;
   for(unsigned int ix=0;ix<outgoing.size();++ix) {
     outColour.push_back(outgoing[ix]->data().iColour());
     switch(outColour.back()) {
     case PDT::Colour0   :     
       singlet.push_back(ix);
       break;
     case PDT::Colour3   :     
       triplet.push_back(ix);
       break;
     case PDT::Colour3bar: 
       antitriplet.push_back(ix);
       break;
     case PDT::Colour8   :     
       octet.push_back(ix);
       break;
     default:
       throw Exception() << "Unknown colour for particle in GeneralThreeBodyDecayer::"
 			<< "colourConnections()" << Exception::runerror;
     }
   }
   // colour neutral decaying particle
   if     ( inColour == PDT::Colour0) {
     // options are all neutral or triplet/antitriplet+ neutral
     if(singlet.size()==3) return;
     else if(singlet.size()==1&&triplet.size()==1&&antitriplet.size()==1) {
       outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
       // add intermediate if needed
       if(inter&&inter->coloured()) {
 	if(inter->dataPtr()->iColour()==PDT::Colour3)
 	  outgoing[triplet[0]]->colourLine()->addColoured(inter);
 	else if(inter->dataPtr()->iColour()==PDT::Colour3bar)
 	  outgoing[triplet[0]]->colourLine()->addAntiColoured(inter);
       }
     }
     else if(octet.size()==1&&triplet.size()==1&&antitriplet.size()==1) {
       outgoing[    triplet[0]]->antiColourNeighbour(outgoing[octet[0]]);
       outgoing[antitriplet[0]]->    colourNeighbour(outgoing[octet[0]]);
       if(inter&&inter->coloured()) {
 	if(inter->dataPtr()->iColour()==PDT::Colour3)
 	  outgoing[antitriplet[0]]->antiColourLine()->addColoured(inter);
 	else if(inter->dataPtr()->iColour()==PDT::Colour3bar)
 	  outgoing[    triplet[0]]->    colourLine()->addAntiColoured(inter);
 	else if(inter->dataPtr()->iColour()==PDT::Colour8) {
 	  outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter);
 	  outgoing[    triplet[0]]->    colourLine()->addColoured(inter);
 	}
       }
     }
     else if(triplet.size()==3) {
       tColinePtr col[3] = {ColourLine::create(outgoing[0]),
 			   ColourLine::create(outgoing[1]),
 			   ColourLine::create(outgoing[2])};
       col[0]->setSourceNeighbours(col[1],col[2]);
     }
     else if(antitriplet.size()==3) {
       tColinePtr col[3] = {ColourLine::create(outgoing[0],true),
 			   ColourLine::create(outgoing[1],true),
 			   ColourLine::create(outgoing[2],true)};
       col[0]->setSinkNeighbours(col[1],col[2]);
     }
     else {
       string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
 	+ out[1]->PDGName() + " " + out[2]->PDGName();
       throw Exception() 
 	<< "Unknown colour structure in GeneralThreeBodyDecayer::"
 	<< "colourConnections() for singlet decaying particle "
 	<< mode << Exception::runerror;
     } 
   }
   // colour triplet decaying particle
   else if( inColour == PDT::Colour3) {
     if(singlet.size()==2&&triplet.size()==1) {
       outgoing[triplet[0]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
       if(inter&&inter->coloured()) 
 	outgoing[triplet[0]]->colourLine()->addColoured(inter);
     }
     else if(antitriplet.size()==1&&triplet.size()==2) {
       if(colourFlow()==0) {
 	outgoing[triplet[0]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
 	outgoing[antitriplet[0]]->colourNeighbour(outgoing[triplet[1]]);
 	if(inter&&inter->coloured()) {
 	  switch (inter->dataPtr()->iColour()) {
 	  case PDT::Colour8:
 	    inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
 	    outgoing[triplet[1]]->colourLine()->addAntiColoured(inter);
 	    break;
 	  default:
 	    string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
 	      + out[1]->PDGName() + " " + out[2]->PDGName();
 	    throw Exception() << "Unknown colour for intermediate in "
 			      << "GeneralThreeBodyDecayer::"
 			      << "colourConnections() for "
 			      << "decaying colour triplet " 
 			      << mode << Exception::runerror;
 	  }
 	}
       }
       else {
 	outgoing[triplet[1]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
 	outgoing[antitriplet[0]]->colourNeighbour(outgoing[triplet[0]]);
 	if(inter&&inter->coloured()) {
 	  switch (inter->dataPtr()->iColour()) {
 	  case PDT::Colour8:
 	    inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
 	    outgoing[triplet[0]]->colourLine()->addAntiColoured(inter);
 	    break;
 	  default:
 	    string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
 	      + out[1]->PDGName() + " " + out[2]->PDGName();
 	    throw Exception() << "Unknown colour for intermediate in "
 			      << "GeneralThreeBodyDecayer::"
 			      << "colourConnections() for "
 			      << "decaying colour triplet " 
 			      << mode << Exception::runerror;
 	  }
 	}
       }
     }
     else if (singlet.size()==1&&triplet.size()==1&&octet.size()==1) {
       if(inter) {
 	if(inter->children()[0]->dataPtr()->iColour()==PDT::Colour8 ||
 	   inter->children()[1]->dataPtr()->iColour()==PDT::Colour8) {
 	  inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
 	  outgoing[octet[0]]->incomingColour(inter);
 	  outgoing[octet[0]]->colourNeighbour(outgoing[triplet[0]]);
 	}
 	else {
 	  outgoing[octet[0]]->incomingColour(inter);
 	  outgoing[octet[0]]->colourNeighbour(inter);
 	  outgoing[triplet[0]]->incomingColour(inter);
 	}
       }
       else {
 	outgoing[octet[0]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
 	outgoing[octet[0]]->colourNeighbour(outgoing[triplet[0]]);
       }
     }
     else if (singlet.size()==1&&antitriplet.size()==2) {
       tColinePtr col[2] = {ColourLine::create(outgoing[antitriplet[0]],true),
 			   ColourLine::create(outgoing[antitriplet[1]],true)};
       parent.colourLine()->setSinkNeighbours(col[0],col[1]);
     }
     else {
       string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
 	+ out[1]->PDGName() + " " + out[2]->PDGName();
       throw Exception() 
 	<< "Unknown colour structure in GeneralThreeBodyDecayer::"
 	<< "colourConnections() for triplet decaying particle " 
 	<< mode << Exception::runerror;
     }
   }
   else if( inColour == PDT::Colour3bar) {
     if(singlet.size()==2&&antitriplet.size()==1) {
       outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
     }
     else if(antitriplet.size()==2&&triplet.size()==1) {
       if(colourFlow()==0) {
 	outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
 	outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[1]]);
 	if(inter&&inter->coloured()) {
 	  switch (inter->dataPtr()->iColour()) {
 	  case PDT::Colour8:
 	    inter->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
 	    outgoing[antitriplet[1]]->antiColourLine()->addAntiColoured(inter);
 	    break;
 	  default:
 	    string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
 	      + out[1]->PDGName() + " " + out[2]->PDGName();
 	    throw Exception() << "Unknown colour for intermediate in"
 			      << " GeneralThreeBodyDecayer::"
 			      << "colourConnections() for "
 			      << "decaying colour antitriplet " 
 			      << mode << Exception::runerror;
 	  }
 	}
       }
       else {
 	outgoing[antitriplet[1]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
 	outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
 	if(inter&&inter->coloured()) {
 	  switch (inter->dataPtr()->iColour()) {
 	  case PDT::Colour8:
 	    inter->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
 	    outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter);
 	    break;
 	  default:
 	    string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
 	      + out[1]->PDGName() + " " + out[2]->PDGName();
 	    throw Exception() << "Unknown colour for intermediate in "
 			      << "GeneralThreeBodyDecayer::"
 			      << "colourConnections() for "
 			      << "decaying colour antitriplet " 
 			      << mode << Exception::runerror;
 	  }
 	}
       }
     }
     else if (singlet.size()==1&&antitriplet.size()==1&&octet.size()==1) {
       if(inter) {
 	if(inter->children()[0]->dataPtr()->iColour()==PDT::Colour8 ||
 	   inter->children()[1]->dataPtr()->iColour()==PDT::Colour8) {
 	  inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
 	  outgoing[octet[0]]->incomingAntiColour(inter);
 	  outgoing[octet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
 	}
 	else {
 	  outgoing[octet[0]]->incomingAntiColour(inter);
 	  outgoing[octet[0]]->antiColourNeighbour(inter);
 	  outgoing[antitriplet[0]]->incomingAntiColour(inter);
 	}
       }
       else {
 	outgoing[octet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
 	outgoing[octet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
       }
     }
     else if (singlet.size()==1&&triplet.size()==2) {
       tColinePtr col[2] = {ColourLine::create(outgoing[triplet[0]]),
 			   ColourLine::create(outgoing[triplet[1]])};
       parent.antiColourLine()->setSourceNeighbours(col[0],col[1]);
     }
     else {
       string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
 	+ out[1]->PDGName() + " " + out[2]->PDGName();
       throw Exception() 
 	<< "Unknown colour structure in GeneralThreeBodyDecayer::"
 	<< "colourConnections() for anti-triplet decaying particle" 
 	<< mode << Exception::runerror;
     }
   }
   else if( inColour == PDT::Colour8) {
     if(triplet.size()==1&&antitriplet.size()==1&&singlet.size()==1) {
       outgoing[    triplet[0]]->incomingColour    (const_ptr_cast<tPPtr>(&parent));
       outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
       if(inter&&inter->coloured()) {
 	switch (inter->dataPtr()->iColour()) {
 	case PDT::Colour3:
 	  outgoing[triplet[0]]->colourLine()->addColoured(inter);
 	  break;
 	case PDT::Colour3bar:
 	  outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter);
 	  break;
 	default:
 	  string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
 	    + out[1]->PDGName() + " " + out[2]->PDGName();
 	  throw Exception() << "Unknown colour for intermediate"
 			    << " in GeneralThreeBodyDecayer::"
 			    << "colourConnections() for "
 			    << "decaying colour octet " 
 			    << mode << Exception::runerror;
 	}
       }
     }
     else if(triplet.size()==3) {
       tColinePtr col[2];
       if(colourFlow()==0) {
 	outgoing[0]->incomingColour    (const_ptr_cast<tPPtr>(&parent));
 	col[0] = ColourLine::create(outgoing[1]);
 	col[1] = ColourLine::create(outgoing[2]);
       }
       else if(colourFlow()==1) {
 	outgoing[1]->incomingColour    (const_ptr_cast<tPPtr>(&parent));
 	col[0] = ColourLine::create(outgoing[0]);
 	col[1] = ColourLine::create(outgoing[2]);
       }
       else if(colourFlow()==2) {
 	outgoing[2]->incomingColour    (const_ptr_cast<tPPtr>(&parent));
 	col[0] = ColourLine::create(outgoing[0]);
 	col[1] = ColourLine::create(outgoing[1]);
       }
       else
 	assert(false);
       parent.antiColourLine()->setSourceNeighbours(col[0],col[1]);
     }
     else if(antitriplet.size()==3) {
       tColinePtr col[2];
       if(colourFlow()==0) {
 	outgoing[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
 	col[0] = ColourLine::create(outgoing[1],true);
 	col[1] = ColourLine::create(outgoing[2],true);
       }
       else if(colourFlow()==1) {
 	outgoing[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
 	col[0] = ColourLine::create(outgoing[0],true);
 	col[1] = ColourLine::create(outgoing[2],true);
       }
       else if(colourFlow()==2) {
 	outgoing[2]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
 	col[0] = ColourLine::create(outgoing[0],true);
 	col[1] = ColourLine::create(outgoing[1],true);
       }
       else
 	assert(false);
       parent.colourLine()->setSinkNeighbours(col[0],col[1]);
     }
     else {
       string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
 	+ out[1]->PDGName() + " " + out[2]->PDGName();
       throw Exception() 
 	<< "Unknown colour structure in GeneralThreeBodyDecayer::"
 	<< "colourConnections() for octet decaying particle" 
 	<< mode << Exception::runerror;
     }
   }
 }
 
 void GeneralThreeBodyDecayer::
 constructIntegratorChannels(vector<int> & intype, vector<Energy> & inmass,
 			    vector<Energy> & inwidth, vector<double> & inpow,
 			    vector<double> & inweights) const {
   // check if any intermediate photons
   bool hasPhoton=false;
   for(unsigned int iy=0;iy<diagmap_.size();++iy) {
     unsigned int ix=diagmap_[iy];
     if(getProcessInfo()[ix].intermediate->id()==ParticleID::gamma)
       hasPhoton = true;
   }
   // loop over channels
   Energy min = incoming()->mass();
   int nchannel(0);
   pair<int,Energy> imin[4]={make_pair(-1,-1.*GeV),make_pair(-1,-1.*GeV),
 			    make_pair(-1,-1.*GeV),make_pair(-1,-1.*GeV)};
   Energy absmin = -1e20*GeV;
   int minType   = -1;
   for(unsigned int iy=0;iy<diagmap_.size();++iy) {
     unsigned int ix=diagmap_[iy];
     if(getProcessInfo()[ix].channelType==TBDiagram::fourPoint) continue;
     Energy dm1(min-getProcessInfo()[ix].intermediate->mass());
     Energy dm2(getProcessInfo()[ix].intermediate->mass());
     int itype(0);
     if     (getProcessInfo()[ix].channelType==TBDiagram::channel23) {
       dm1 -= outgoing()[0]->mass();
       dm2 -= outgoing()[1]->mass()+outgoing()[2]->mass();
       itype = 3;
     }
     else if(getProcessInfo()[ix].channelType==TBDiagram::channel13) {
       dm1 -= outgoing()[1]->mass();
       dm2 -= outgoing()[0]->mass()+outgoing()[2]->mass();
       itype = 2;
     }
     else if(getProcessInfo()[ix].channelType==TBDiagram::channel12) {
       dm1 -= outgoing()[2]->mass();
       dm2 -= outgoing()[0]->mass()+outgoing()[1]->mass();
       itype = 1;
     }
     if((dm1<ZERO||dm2<ZERO)&&!hasPhoton) {
       if (imin[itype].first < 0  ||
 	  (dm1<ZERO && imin[itype].second < dm1)  ) {
 	imin[itype] = make_pair(ix,dm1);
 	if(dm1<ZERO&&absmin<dm1) {
 	  absmin = dm1;
 	  minType = itype;
 	}
       }
       continue;
     }
     if(getProcessInfo()[ix].intermediate->id()!=ParticleID::gamma) {
       intype.push_back(itype);
       inpow.push_back(0.);
       inmass.push_back(getProcessInfo()[ix].intermediate->mass());
       inwidth.push_back(widthOption() ==3 ? ZERO : getProcessInfo()[ix].intermediate->width());
       ++nchannel;
     }
     else if(getProcessInfo()[ix].intermediate->id()==ParticleID::gamma) {
       intype.push_back(itype);
       inpow.push_back(-2.);
       inmass.push_back(-1.*GeV);
       inwidth.push_back(-1.*GeV);
       ++nchannel;
     }
   }
   // physical poles, use them and return
   if(nchannel>0) {
     inweights = vector<double>(nchannel,1./double(nchannel));
     return;
   }
   // use shallowest pole
   else if(intOpt_==1&&minType>0&&getProcessInfo()[imin[minType].first].intermediate->id()!=ParticleID::gamma) {
     intype.push_back(minType);
     inpow.push_back(0.);
     inmass.push_back(getProcessInfo()[imin[minType].first].intermediate->mass());
     inwidth.push_back(widthOption() ==3 ? ZERO : getProcessInfo()[imin[minType].first].intermediate->width());
     inweights = vector<double>(1,1.);
     return;
   }
   for(unsigned int ix=1;ix<4;++ix) {
     if(imin[ix].first>=0) {
       intype.push_back(ix);
       if(getProcessInfo()[imin[ix].first].intermediate->id()!=ParticleID::gamma) {
 	inpow.push_back(0.);
 	inmass.push_back(getProcessInfo()[imin[ix].first].intermediate->mass());
 	inwidth.push_back(widthOption() ==3 ? ZERO : getProcessInfo()[imin[ix].first].intermediate->width());
       }
       else {
 	inpow.push_back(-2.);
 	inmass.push_back(-1.*GeV);
 	inwidth.push_back(-1.*GeV);
       }
       ++nchannel;
     }
   }
   inweights = vector<double>(nchannel,1./double(nchannel));
 }
 
 bool GeneralThreeBodyDecayer::setColourFactors(double symfac) {
   string name = incoming_->PDGName() + "->";
   vector<int> sng,trip,atrip,oct;
   unsigned int iloc(0);
 
   for(vector<PDPtr>::const_iterator it = outgoing_.begin();
       it != outgoing_.end();++it) {
     name += (**it).PDGName() + " ";
     if     ((**it).iColour() == PDT::Colour0    ) sng.push_back(iloc) ;
     else if((**it).iColour() == PDT::Colour3    ) trip.push_back(iloc) ;
     else if((**it).iColour() == PDT::Colour3bar ) atrip.push_back(iloc);
     else if((**it).iColour() == PDT::Colour8    ) oct.push_back(iloc) ;
     ++iloc;
   }
   // colour neutral decaying particle
   if     ( incoming_->iColour() == PDT::Colour0) {
     // options are all neutral or triplet/antitriplet+ neutral
     if(sng.size()==3) {
       nflow_ = 1;
       colour_        = vector<DVector>(1,DVector(1,1.));
       colourLargeNC_ = vector<DVector>(1,DVector(1,1.));
     }
     else if(sng.size()==1&&trip.size()==1&&atrip.size()==1) {
       nflow_ = 1;
       colour_         = vector<DVector>(1,DVector(1,3.));
       colourLargeNC_  = vector<DVector>(1,DVector(1,3.));
     }
     else if(trip.size()==1&&atrip.size()==1&&oct.size()==1) {
       nflow_ = 1;
       colour_         = vector<DVector>(1,DVector(1,4.));
       colourLargeNC_  = vector<DVector>(1,DVector(1,4.));
     }
     else if( trip.size() == 3 || atrip.size() == 3 ) {
       nflow_ = 1;
       colour_         = vector<DVector>(1,DVector(1,6.));
       colourLargeNC_  = vector<DVector>(1,DVector(1,6.));
       for(unsigned int ix=0;ix<diagrams_.size();++ix) {
 	double sign = diagrams_[ix].channelType == TBDiagram::channel13 ? -1. : 1.;
 	diagrams_[ix].       colourFlow = vector<CFPair>(1,make_pair(1,sign));
 	diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(1,sign));
       }
     }
     else {
       generator()->log() << "Unknown colour flow structure for "
 			 << "colour neutral decay " 
 			 << name  << " in GeneralThreeBodyDecayer::"
 			 << "setColourFactors(), omitting decay\n";
       return false;
     }
   }
   // colour triplet decaying particle
   else if( incoming_->iColour() == PDT::Colour3) {
     if(sng.size()==2&&trip.size()==1) {
       nflow_ = 1;
       colour_        = vector<DVector>(1,DVector(1,1.));
       colourLargeNC_ = vector<DVector>(1,DVector(1,1.));
     }
     else if(trip.size()==2&&atrip.size()==1) {
       nflow_ = 2;
       colour_.clear();
       colour_.resize(2,DVector(2,0.));
       colour_[0][0] = 3.; colour_[0][1] = 1.;
       colour_[1][0] = 1.; colour_[1][1] = 3.;
       colourLargeNC_.clear();
       colourLargeNC_.resize(2,DVector(2,0.));
       colourLargeNC_[0][0] = 3.; colourLargeNC_[1][1] = 3.;
       // sort out the contribution of the different diagrams to the colour
       // flows
       for(unsigned int ix=0;ix<diagrams_.size();++ix) {
 	// colour singlet intermediate
 	if(diagrams_[ix].intermediate->iColour()==PDT::Colour0) {
 	  if(diagrams_[ix].channelType==trip[0]) {
 	    diagrams_[ix].       colourFlow = vector<CFPair>(1,make_pair(1,1.));
 	    diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
 	  }
 	  else {
 	    diagrams_[ix].colourFlow        = vector<CFPair>(1,make_pair(2,1.));
 	    diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
 	  }
 	}
 	// colour octet intermediate
 	else if(diagrams_[ix].intermediate->iColour()==PDT::Colour8) {
 	  if(diagrams_[ix].channelType==trip[0]) {
 	    vector<CFPair> flow(1,make_pair(2, 0.5  ));
 	    diagrams_[ix].largeNcColourFlow = flow;
 	    flow.push_back(       make_pair(1,-1./6.));
 	    diagrams_[ix].colourFlow=flow;
 	  }
 	  else {
 	    vector<CFPair> flow(1,make_pair(1, 0.5  ));
 	    diagrams_[ix].largeNcColourFlow = flow;
 	    flow.push_back(       make_pair(2,-1./6.));
 	    diagrams_[ix].colourFlow=flow;
 	  }
 	}
 	else {
 	  generator()->log() << "Unknown colour for the intermediate in "
 			     << "triplet -> triplet triplet antitriplet in "
 			     << "GeneralThreeBodyDecayer::setColourFactors()"
 			     << " for " << name << " omitting decay\n";
 	  return false;
 	}
       }
     }
     else if(trip.size()==1&&oct.size()==1&&sng.size()==1) {
       nflow_ = 1;
       colour_        = vector<DVector>(1,DVector(1,4./3.));
       colourLargeNC_ = vector<DVector>(1,DVector(1,4./3.));
     }
     else if(sng.size()==1&&atrip.size()==2) {
       nflow_ = 1;
       colour_        = vector<DVector>(1,DVector(1,2.));
       colourLargeNC_ = vector<DVector>(1,DVector(1,2.));
     }
     else {
       generator()->log() << "Unknown colour structure for "
 			 << "triplet decay in "
 			 << "GeneralThreeBodyDecayer::setColourFactors()"
 			 << " for " << name << " omitting decay\n";
       return false;
     }
   }
   // colour antitriplet decaying particle
   else if( incoming_->iColour() == PDT::Colour3bar) {
     if(sng.size()==2&&atrip.size()==1) {
       nflow_ = 1;
       colour_        = vector<DVector>(1,DVector(1,1.));
       colourLargeNC_ = vector<DVector>(1,DVector(1,1.));
     }
     else if(atrip.size()==2&&trip.size()==1) {
       nflow_ = 2;
       colour_.clear();
       colour_.resize(2,DVector(2,0.));
       colour_[0][0] = 3.; colour_[0][1] = 1.;
       colour_[1][0] = 1.; colour_[1][1] = 3.;
       colourLargeNC_.clear();
       colourLargeNC_.resize(2,DVector(2,0.));
       colourLargeNC_[0][0] = 3.; colourLargeNC_[1][1] = 3.;
       // sort out the contribution of the different diagrams to the colour
       // flows
       for(unsigned int ix=0;ix<diagrams_.size();++ix) {
 	// colour singlet intermediate
 	if(diagrams_[ix].intermediate->iColour()==PDT::Colour0) {
 	  if(diagrams_[ix].channelType==atrip[0]) {
 	    diagrams_[ix].       colourFlow = vector<CFPair>(1,make_pair(1,1.));
 	    diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
 	  }
 	  else {
 	    diagrams_[ix].colourFlow        = vector<CFPair>(1,make_pair(2,1.));
 	    diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
 	  }
 	}
 	// colour octet intermediate
 	else if(diagrams_[ix].intermediate->iColour()==PDT::Colour8) {
 	  if(diagrams_[ix].channelType==atrip[0]) {
 	    vector<CFPair> flow(1,make_pair(2, 0.5  ));
 	    diagrams_[ix].largeNcColourFlow = flow;
 	    flow.push_back(       make_pair(1,-1./6.));
 	    diagrams_[ix].colourFlow=flow;
 	  }
 	  else {
 	    vector<CFPair> flow(1,make_pair(1, 0.5  ));
 	    diagrams_[ix].largeNcColourFlow = flow;
 	    flow.push_back(       make_pair(2,-1./6.));
 	    diagrams_[ix].colourFlow=flow;
 	  }
 	}
 	else {
 	  generator()->log() << "Unknown colour for the intermediate in "
 			     << "antitriplet -> antitriplet antitriplet triplet in "
 			     << "GeneralThreeBodyDecayer::setColourFactors()"
 			     << " for " << name << " omitting decay\n";
 	  return false;
 	}
       }
     }
     else if(atrip.size()==1&&oct.size()==1&&sng.size()==1) {
       nflow_ = 1;
       colour_        = vector<DVector>(1,DVector(1,4./3.));
       colourLargeNC_ = vector<DVector>(1,DVector(1,4./3.));
     }
     else if(sng.size()==1&&trip.size()==2) {
       nflow_ = 1;
       colour_        = vector<DVector>(1,DVector(1,2.));
       colourLargeNC_ = vector<DVector>(1,DVector(1,2.));
     }
     else {
       generator()->log() << "Unknown colour antitriplet decay in "
 			 << "GeneralThreeBodyDecayer::setColourFactors()"
 			 << " for " << name << " omitting decay\n";
       return false;
     }
   }
   // colour octet particle
   else if( incoming_->iColour() == PDT::Colour8) {
     // triplet antitriplet
     if(trip.size() == 1 && atrip.size() == 1 && sng.size() == 1) {
       nflow_ = 1;
       colour_        = vector<DVector>(1,DVector(1,0.5));
       colourLargeNC_ = vector<DVector>(1,DVector(1,0.5));
     }
     // three (anti)triplets
     else if(trip.size()==3||atrip.size()==3) {
       nflow_ = 3;
       colour_        = vector<DVector>(3,DVector(3,0.));
       colourLargeNC_ = vector<DVector>(3,DVector(3,0.));
       colour_[0][0] = 1.; colour_[1][1] = 1.; colour_[2][2] = 1.;
       colour_[0][1] = -0.5; colour_[1][0] = -0.5;
       colour_[0][2] = -0.5; colour_[2][0] = -0.5;
       colour_[1][2] = -0.5; colour_[2][1] = -0.5;
       colourLargeNC_ = vector<DVector>(3,DVector(3,0.));
       colourLargeNC_[0][0] = 1.; colourLargeNC_[1][1] = 1.; colourLargeNC_[2][2] = 1.;
       // sett the factors for the diagrams
       for(unsigned int ix=0;ix<diagrams_.size();++ix) {
 	tPDPtr inter = diagrams_[ix].intermediate;
 	if(inter->CC()) inter = inter->CC();
 	unsigned int io[2]={1,2};
 	double sign = diagrams_[ix].channelType == TBDiagram::channel13 ? -1. : 1.;
 	for(unsigned int iy=0;iy<3;++iy) {
 	  if     (iy==1) io[0]=0;
 	  else if(iy==2) io[1]=1;
 	  tPDVector decaylist = diagrams_[ix].vertices.second->search(iy, inter);
 	  if(decaylist.empty()) continue;
 	  bool found=false;
 	  for(unsigned int iz=0;iz<decaylist.size();iz+=3) {	    
 	    if(decaylist[iz+io[0]]->id()==diagrams_[ix].outgoingPair.first &&
 	       decaylist[iz+io[1]]->id()==diagrams_[ix].outgoingPair.second) {
 	      sign *= 1.;
 	      found = true;
 	    }
 	    else if(decaylist[iz+io[0]]->id()==diagrams_[ix].outgoingPair.second &&
 		    decaylist[iz+io[1]]->id()==diagrams_[ix].outgoingPair.first ) {
 	      sign *= -1.;
 	      found = true;
 	    }
 	  }
 	  if(found) {
 	    if(iy==1) sign *=-1.;
 	    break;
 	  }
 	}
 	diagrams_[ix].       colourFlow = vector<CFPair>(1,make_pair(diagrams_[ix].channelType+1,sign));
 	diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(diagrams_[ix].channelType+1,sign));
       }
     }
     // unknown
     else {
       generator()->log() << "Unknown colour octet decay in "
 			 << "GeneralThreeBodyDecayer::setColourFactors()"
 			 << " for " << name << " omitting decay\n";
       return false;
     }
   }
   else if (incoming_->iColour() == PDT::Colour6 ) {
     generator()->log() << "Unknown colour sextet decay in "
 		       << "GeneralThreeBodyDecayer::setColourFactors()"
 		       << " for " << name << " omitting decay\n";
     return false;
   }
   else if (incoming_->iColour() == PDT::Colour6bar ) {
     generator()->log() << "Unknown colour anti-sextet decay in "
 		       << "GeneralThreeBodyDecayer::setColourFactors()"
 		       << " for " << name << " omitting decay\n";
     return false;
   }
 
   assert(nflow_ != 999);
 
   for(unsigned int ix=0;ix<nflow_;++ix) {
     for(unsigned int iy=0;iy<nflow_;++iy) {
       colour_       [ix][iy] /= symfac;
       colourLargeNC_[ix][iy] /= symfac;
     }
   }
   if( Debug::level > 1 ) {
     generator()->log() << "Mode: " << name << " has colour factors\n";
     for(unsigned int ix=0;ix<nflow_;++ix) {
       for(unsigned int iy=0;iy<nflow_;++iy) {
 	generator()->log() << colour_[ix][iy] << " ";
       }
       generator()->log() << "\n";
     }
     for(unsigned int ix=0;ix<diagrams_.size();++ix) {
       generator()->log() << "colour flow for diagram : " << ix;
       for(unsigned int iy=0;iy<diagrams_[ix].colourFlow.size();++iy)
 	generator()->log() << "(" << diagrams_[ix].colourFlow[iy].first  << "," 
 			   << diagrams_[ix].colourFlow[iy].second << "); ";
       generator()->log() << "\n";
     }
   }
   return true;
 }
 
 void GeneralThreeBodyDecayer::setupDiagrams(bool kinCheck) {
   clearModes();
   // create the phase space integrator
   tPDVector extpart(1,incoming_);
   extpart.insert(extpart.end(),outgoing_.begin(),outgoing_.end());
   // create the integration channels for the decay
   DecayPhaseSpaceModePtr mode(new_ptr(DecayPhaseSpaceMode(extpart,this,true)));
   DecayPhaseSpaceChannelPtr newchannel;
   // create the phase-space channels for the integration
   unsigned int nmode(0);
   unsigned int idiag(0);
   for(vector<TBDiagram>::iterator it = diagrams_.begin();it!=diagrams_.end();++it) {
     if(it->channelType==TBDiagram::fourPoint||
        it->channelType==TBDiagram::UNDEFINED) {
       idiag+=1;
       continue;
     }
     // create the new channel
     newchannel=new_ptr(DecayPhaseSpaceChannel(mode));
     int jac = 0;
     double power = 0.0;
     if ( it->intermediate->mass() == ZERO ||
   	 it->intermediate->width() == ZERO ) {
       jac = 1;
       power = -2.0;
     }
     if(it->channelType==TBDiagram::channel23) {
       newchannel->addIntermediate(extpart[0],0,0.0,-1,1);
       newchannel->addIntermediate(it->intermediate,jac,power, 2,3);
     }
     else if(it->channelType==TBDiagram::channel13) {
       newchannel->addIntermediate(extpart[0],0,0.0,-1,2);
       newchannel->addIntermediate(it->intermediate,jac,power, 1,3);
     }
     else if(it->channelType==TBDiagram::channel12) {
       newchannel->addIntermediate(extpart[0],0,0.0,-1,3);
       newchannel->addIntermediate(it->intermediate,jac,power, 1,2);
     }
     newchannel->init();
     if(kinCheck&&!newchannel->checkKinematics()) {
       generator()->log() << "Erasing diagram "
 			 << *it
 			 << "from three body decay as zero width propagator can be on-shell,\n"
 			 << "hopefully this diagram is zero in your model, but you should check this\n";
       it = diagrams_.erase(it);
       if(it == diagrams_.end()) break;
       continue;
     }
     diagmap_.push_back(idiag);
     mode->addChannel(newchannel);
     ++nmode;
     ++idiag;
   }
   if(nmode==0) {
     string mode = extpart[0]->PDGName() + "->";
     for(unsigned int ix=1;ix<extpart.size();++ix) mode += extpart[ix]->PDGName() + " ";
-    throw Exception() << "No decay channels in GeneralThreeBodyDecayer::"
-  		      << "doinit() for " << mode << "\n" << Exception::runerror;
+    generator()->log() << "No decay channels in GeneralThreeBodyDecayer::"
+		       << "setupDiagrams() for " << mode << "\n";
+    return;
   }
-  // // add the mode
+  // add the mode
   vector<double> wgt(nmode,1./double(nmode));
   addMode(mode,1.,wgt);
 }