diff --git a/Decay/DecayIntegrator.cc b/Decay/DecayIntegrator.cc
--- a/Decay/DecayIntegrator.cc
+++ b/Decay/DecayIntegrator.cc
@@ -1,292 +1,292 @@
 // -*- C++ -*-
 //
 // DecayIntegrator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2019 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.
 //
 
 #include "DecayIntegrator.h"
 #include "PhaseSpaceMode.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/EventRecord/Particle.h"
 #include "ThePEG/Repository/UseRandom.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "Herwig/PDT/WidthCalculatorBase.h"
 
 using namespace Herwig;
 
 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,unsigned int> interfaceIteration
     ("Iteration",
      "Number of iterations for the initialization of the phase space",
      &DecayIntegrator::nIter_, 10, 0, 100,
      false, false, true);  
   
   static Parameter<DecayIntegrator,unsigned int> interfacePoints
     ("Points",
      "number of phase space points to generate in the initialisation.",
      &DecayIntegrator::nPoint_, 10000, 1, 1000000000,
      false, false, true);
   
   static Parameter<DecayIntegrator,unsigned 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);
 }
 
 double DecayIntegrator::oneLoopVirtualME(unsigned int ,
 					  const Particle &, 
 					  const ParticleVector &) {
   throw Exception()
     << "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 Exception()
     << "DecayIntegrator::realEmmisionME() called. This should"
     << " have been overidden in an inheriting class if it is used"
     << Exception::runerror;
 }
 
 ParticleVector DecayIntegrator::decay(const Particle & parent,
-				       const tPDVector & children) const {
+				      const tPDVector & children) const {
   // return empty vector if products heavier than parent
   Energy mout(ZERO);
   for(tPDPtr pd : children)  mout += pd->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_]->generateDecay(parent,this,generateInter_,cc);
 }
 
 void DecayIntegrator::doinitrun() {
   HwDecayerBase::doinitrun();
   if ( initialize() && Debug::level > 1 ) 
     CurrentGenerator::current().log() << "Start of the initialisation for " 
 				      << name() << "\n";
   for(unsigned int ix=0;ix<modes_.size();++ix) {
     if(!modes_[ix]) continue;
     modes_[ix]->initrun();
     iMode_=ix;
     modes_[ix]->initializePhaseSpace(initialize(),this);
   }
 }
 
 // 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";
   }
 }
 
 // set the code for the partial width
 void DecayIntegrator::setPartialWidth(const DecayMode & dm, int imode) {
   int ifound = findMode(dm);
   if(ifound>=0) modes_[ifound]->setPartialWidth(imode);
 }
 
 WidthCalculatorBasePtr 
 DecayIntegrator::threeBodyMEIntegrator(const DecayMode &) const {
   return WidthCalculatorBasePtr();
 }
 
 int DecayIntegrator::findMode(const DecayMode & dm) {
   int imode(-1);
   vector<int> extid;
   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;
     }
     tcPDPtr in = modes_[ix]->incoming().first;
     tcPDPtr cc = modes_[ix]->incoming().first->CC();
     tmax=1;if(!cc){++tmax;}
     for(iz=0;iz<tmax;++iz) {
       extid.clear();
       // check the parent
       if(dm.parent()!=in && dm.parent()!=cc) continue;
       if(dm.parent()->id()==in->id()&&iz==0) {
 	for(iy=0,N=modes_[ix]->numberOfParticles();iy<N;++iy) {
 	  extid.push_back(modes_[ix]->outgoing()[iy]->id());
 	}
       }
       else if(dm.parent()->id()==in->id()&&iz==1) {
 	for(iy=0,N=modes_[ix]->numberOfParticles();iy<N;++iy) {
 	  tcPDPtr cc2=modes_[ix]->outgoing()[iy]->CC();
 	  extid.push_back( cc2 ? cc2->id() : modes_[ix]->outgoing()[iy]->id());
 	}
       }
       else if(cc&&dm.parent()->id()==cc->id()) {
 	for(iy=0,N=modes_[ix]->numberOfParticles();iy<N;++iy) {
 	  tcPDPtr cc2 = modes_[ix]->outgoing()[iy]->CC();
 	  extid.push_back( cc2 ? cc2->id() : modes_[ix]->outgoing()[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=0;
 	  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()) {
 	  imode=ix;
 	  found=true;
 	}
       }
     }
     ++ix;
   }
   while(!found&&ix<modes_.size());
   return imode;
 }
 
 // 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 Exception() 
     << "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 Exception() 
     << "Calling the virtual DecayIntegrator::threeBodydGammads()" 
     <<"method. This must be overwritten in the classes "
     << "inheriting from DecayIntegrator where it is needed"
     << Exception::runerror;
 }
 
 // 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]->generateDecay(inpart,this,inter,cc);
 }
 
 void  DecayIntegrator::addMode(PhaseSpaceModePtr mode) const {
   modes_.push_back(mode);
   if(mode) mode->init();
 }
 
 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;
 }
   
 // 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);
   }
 }
 
 Energy DecayIntegrator::initializePhaseSpaceMode(unsigned int imode,bool init, bool onShell) const{
   tcPhaseSpaceModePtr cmodeptr=mode(imode);
   tPhaseSpaceModePtr modeptr = const_ptr_cast<tPhaseSpaceModePtr>(cmodeptr);
   modeptr->init();
   return modeptr->initializePhaseSpace(init,this,onShell);
 }
diff --git a/Decay/Perturbative/SMWDecayer.cc b/Decay/Perturbative/SMWDecayer.cc
--- a/Decay/Perturbative/SMWDecayer.cc
+++ b/Decay/Perturbative/SMWDecayer.cc
@@ -1,740 +1,740 @@
 // -*- C++ -*-
 //
 // SMWDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2019 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 SMWDecayer class.
 //
 
 #include "SMWDecayer.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/ParVector.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/PDT/DecayMode.h"
 #include "Herwig/Decay/DecayVertex.h"
 #include "ThePEG/Helicity/VectorSpinInfo.h"
 #include "ThePEG/Helicity/FermionSpinInfo.h"
 #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
 #include "Herwig/Models/StandardModel/StandardModel.h"
 #include "Herwig/Shower/RealEmissionProcess.h"
 #include "Herwig/Decay/GeneralDecayMatrixElement.h"
 #include <numeric>
 
 using namespace Herwig;
 using namespace ThePEG::Helicity;
 
 const double SMWDecayer::EPS_=0.00000001;
 
 SMWDecayer::SMWDecayer()
   : quarkWeight_(6,0.), leptonWeight_(3,0.), CF_(4./3.),
     NLO_(false) {
   quarkWeight_[0]  = 1.01596;
   quarkWeight_[1]  = 0.0537308;
   quarkWeight_[2]  = 0.0538085;
   quarkWeight_[3]  = 1.01377;
   quarkWeight_[4]  = 1.45763e-05;
   quarkWeight_[5]  = 0.0018143;
   leptonWeight_[0] = 0.356594;
   leptonWeight_[1] = 0.356593;
   leptonWeight_[2] = 0.356333;
   // intermediates
   generateIntermediates(false);
 }
 
 void SMWDecayer::doinit() {
   PerturbativeDecayer::doinit();
   // get the vertices from the Standard Model object
   tcHwSMPtr hwsm=dynamic_ptr_cast<tcHwSMPtr>(standardModel());
   if(!hwsm) throw InitException() << "Must have Herwig StandardModel object in"
 				  << "SMWDecayer::doinit()"
 				  << Exception::runerror;
   FFWVertex_ = hwsm->vertexFFW();
   FFGVertex_ = hwsm->vertexFFG();
   WWWVertex_ = hwsm->vertexWWW();
   FFPVertex_ = hwsm->vertexFFP();
   // make sure they are initialized
   FFGVertex_->init();
   FFWVertex_->init();
   WWWVertex_->init();
   FFPVertex_->init();
   // now set up the decay modes
   // W modes
   tPDPtr Wp = getParticleData(ParticleID::Wplus);
   // loop for the quarks
   unsigned int iz=0;
   for(int ix=1;ix<6;ix+=2) {
     for(int iy=2;iy<6;iy+=2) {
       // check that the combination of particles is allowed
       if(!FFWVertex_->allowed(-ix,iy,ParticleID::Wminus))
 	throw InitException() << "SMWDecayer::doinit() the W vertex" 
 			      << "cannot handle all the quark modes" 
 			      << Exception::abortnow;
       tPDVector out = {getParticleData(-ix),
 		       getParticleData( iy)};
       PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(Wp,out,quarkWeight_[iz]));
       addMode(mode);
       ++iz;
     }
   }
   // loop for the leptons
   for(int ix=11;ix<17;ix+=2) {
     tPDVector out = {getParticleData(-ix  ),
 		     getParticleData( ix+1)};
     PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(Wp,out,leptonWeight_[(ix-11)/2]));
     addMode(mode);
   }
   gluon_ = getParticleData(ParticleID::g);
 }
 
 int SMWDecayer::modeNumber(bool & cc,tcPDPtr parent, 
 			    const tPDVector & children) const {
   int imode(-1);
   if(children.size()!=2) return imode;
   int id0=parent->id();
   tPDVector::const_iterator pit = children.begin();
   int id1=(**pit).id();
   ++pit;
   int id2=(**pit).id();
   if(abs(id0)!=ParticleID::Wplus) return imode;
   int idd(0),idu(0);
   if(abs(id1)%2==1&&abs(id2)%2==0) {
     idd=abs(id1);
     idu=abs(id2);
   }
   else if(abs(id1)%2==0&&abs(id2)%2==1) {
     idd=abs(id2);
     idu=abs(id1);
   }
   if(idd==0&&idu==0) {
     return imode;
   }
   else if(idd<=5) {
     imode=idd+idu/2-2;
   }
   else {
     imode=(idd-1)/2+1;
   }
   cc= (id0==ParticleID::Wminus);
   return imode;
 }
 
 ParticleVector SMWDecayer::decay(const Particle & parent,
 				 const tPDVector & children) const {
   // generate the decay
   bool cc;
   unsigned int imode = modeNumber(cc,parent.dataPtr(),children);
-  ParticleVector output = generate(false,false,imode,parent);
+  ParticleVector output = generate(false,cc,imode,parent);
   if(output[0]->hasColour())      output[0]->antiColourNeighbour(output[1]);
   else if(output[1]->hasColour()) output[1]->antiColourNeighbour(output[0]);
   return output;
 }
 
 void SMWDecayer::persistentOutput(PersistentOStream & os) const {
   os << FFWVertex_ << quarkWeight_ << leptonWeight_
      << FFGVertex_ << gluon_ << NLO_
      << WWWVertex_ << FFPVertex_;  
 }
 
 void SMWDecayer::persistentInput(PersistentIStream & is, int) {
   is >> FFWVertex_ >> quarkWeight_ >> leptonWeight_
      >> FFGVertex_ >> gluon_ >> NLO_
      >> WWWVertex_ >> FFPVertex_;
 }
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<SMWDecayer,PerturbativeDecayer>
 describeHerwigSMWDecayer("Herwig::SMWDecayer", "HwPerturbativeDecay.so");
 
 void SMWDecayer::Init() {
 
   static ClassDocumentation<SMWDecayer> documentation
     ("The SMWDecayer class is the implementation of the decay"
      " of the W boson to the Standard Model fermions.");
 
   static ParVector<SMWDecayer,double> interfaceWquarkMax
     ("QuarkMax",
      "The maximum weight for the decay of the W to quarks",
      &SMWDecayer::quarkWeight_,
      0, 0, 0, -10000, 10000, false, false, true);
 
   static ParVector<SMWDecayer,double> interfaceWleptonMax
     ("LeptonMax",
      "The maximum weight for the decay of the W to leptons",
      &SMWDecayer::leptonWeight_,
      0, 0, 0, -10000, 10000, false, false, true);
 
   static Switch<SMWDecayer,bool> interfaceNLO
     ("NLO",
      "Whether to return the LO or NLO result",
      &SMWDecayer::NLO_, false, false, false);
   static SwitchOption interfaceNLOLO
     (interfaceNLO,
      "No",
      "Leading-order result",
      false);
   static SwitchOption interfaceNLONLO
     (interfaceNLO,
      "Yes",
      "NLO result",
      true);
 
 }
 void SMWDecayer::
 constructSpinInfo(const Particle & part, ParticleVector decay) const {
   int iferm(1),ianti(0);
   if(decay[0]->id()>0) swap(iferm,ianti);
   VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast<tPPtr>(&part),
 					incoming,true,false);
   SpinorBarWaveFunction::
     constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
   SpinorWaveFunction::
     constructSpinInfo(wave_   ,decay[ianti],outgoing,true);
 }
 
 // return the matrix element squared
 double SMWDecayer::me2(const int,const Particle & part,
 		       const tPDVector & outgoing,
 		       const vector<Lorentz5Momentum> & momenta,
 		       MEOption meopt) const {
   if(!ME()) 
     ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin1Half,PDT::Spin1Half)));
   int iferm(1),ianti(0);
   if(outgoing[0]->id()>0) swap(iferm,ianti);
   if(meopt==Initialize) {
     VectorWaveFunction::calculateWaveFunctions(vectors_,rho_,
   					       const_ptr_cast<tPPtr>(&part),
   					       incoming,false);
     // fix rho if no correlations
     fixRho(rho_);
   }
   SpinorBarWaveFunction wbar(momenta[iferm],outgoing[iferm],Helicity::outgoing);
   SpinorWaveFunction    w   (momenta[ianti],outgoing[ianti],Helicity::outgoing);
   wavebar_.resize(2);
   wave_   .resize(2);
   for(unsigned int ihel=0;ihel<2;++ihel) {
     wbar.reset(ihel);
     wavebar_[ihel] = wbar;
     w.reset(ihel);
     wave_   [ihel] = w;
   }
   // compute the matrix element
   Energy2 scale(sqr(part.mass()));
   for(unsigned int ifm=0;ifm<2;++ifm) {
     for(unsigned int ia=0;ia<2;++ia) {
       for(unsigned int vhel=0;vhel<3;++vhel) {
   	if(iferm>ianti) (*ME())(vhel,ia,ifm)=
   	  FFWVertex_->evaluate(scale,wave_[ia],wavebar_[ifm],vectors_[vhel]);
   	else            (*ME())(vhel,ifm,ia)=
   	  FFWVertex_->evaluate(scale,wave_[ia],wavebar_[ifm],vectors_[vhel]);
       }
     }
   }
   double output=(ME()->contract(rho_)).real()*UnitRemoval::E2/scale;
   if(abs(outgoing[0]->id())<=6) output*=3.;
   // // leading-order result
   if(!NLO_) return output;
   // check decay products coloured, otherwise return
   if(!outgoing[0]->coloured()) return output;
   // inital masses, couplings  etc
   // W mass
   mW_ = part.mass();
   // strong coupling
   aS_ = SM().alphaS(sqr(mW_));
   // reduced mass
   double mu1  = outgoing[0]->mass()/mW_;
   double mu2  = outgoing[1]->mass()/mW_;
   // scale
   scale_ = sqr(mW_);
   // now for the nlo loop correction
   double virt = CF_*aS_/Constants::pi;
   // now for the real correction
   double realFact=0.;
   for(int iemit=0;iemit<2;++iemit) {
     double phi  = UseRandom::rnd()*Constants::twopi;
     // set the emitter and the spectator
     double muj  = iemit==0 ? mu1 : mu2;
     double muk  = iemit==0 ? mu2 : mu1;
     double muj2 = sqr(muj);
     double muk2 = sqr(muk);
     // calculate y
     double yminus = 0.; 
     double yplus  = 1.-2.*muk*(1.-muk)/(1.-muj2-muk2);
     double y = yminus + UseRandom::rnd()*(yplus-yminus);
     double v = sqrt(sqr(2.*muk2 + (1.-muj2-muk2)*(1.-y))-4.*muk2)
       /(1.-muj2-muk2)/(1.-y);
     double zplus  = (1.+v)*(1.-muj2-muk2)*y/2./(muj2+(1.-muj2-muk2)*y);
     double zminus = (1.-v)*(1.-muj2-muk2)*y/2./(muj2+(1.-muj2-muk2)*y);
     double z = zminus + UseRandom::rnd()*(zplus-zminus);
     double jac = (1.-y)*(yplus-yminus)*(zplus-zminus);
     // calculate x1,x2,x3,xT
     double x2 = 1.-y*(1.-muj2-muk2)-muj2+muk2;
     double x1 = 1.+muj2-muk2-z*(x2-2.*muk2);
     // copy the particle objects over for calculateRealEmission
     tcPDVector outgoing = {part.dataPtr(),outgoing[0],outgoing[1],gluon_};
     vector<Lorentz5Momentum> mom = {part.momentum(),momenta[0],momenta[1]};
     realFact += 0.25*jac*sqr(1.-muj2-muk2)/
       sqrt((1.-sqr(muj-muk))*(1.-sqr(muj+muk)))/Constants::twopi
       *2.*CF_*aS_*calculateRealEmission(x1, x2, outgoing,mom, phi, 
   					muj, muk, iemit, true);
   }
   // the born + virtual + real
   output *= (1. + virt + realFact);
   return output;
 }
 
 void SMWDecayer::doinitrun() {
   PerturbativeDecayer::doinitrun();
   if(initialize()) {
     for(unsigned int ix=0;ix<numberModes();++ix) {
       if(ix<6) quarkWeight_ [ix]=mode(ix)->maxWeight();
       else     leptonWeight_[ix-6]=mode(ix)->maxWeight();
     }
   }
 }
 
 void SMWDecayer::dataBaseOutput(ofstream & output,
 				 bool header) const {
   if(header) output << "update decayers set parameters=\"";
   for(unsigned int ix=0;ix<quarkWeight_.size();++ix) {
     output << "newdef " << name() << ":QuarkMax " << ix << " "
 	   << quarkWeight_[ix] << "\n";
   }
   for(unsigned int ix=0;ix<leptonWeight_.size();++ix) {
     output << "newdef " << name() << ":LeptonMax " << ix << " "
 	   << leptonWeight_[ix] << "\n";
   }
   // parameters for the PerturbativeDecayer base class
   PerturbativeDecayer::dataBaseOutput(output,false);
   if(header) output << "\n\" where BINARY ThePEGName=\"" 
 		    << fullName() << "\";" << endl;
 }
 
 
 void SMWDecayer::
 initializeMECorrection(RealEmissionProcessPtr born, double & initial,
 		       double & final) {
   // get the quark and antiquark
   ParticleVector qq; 
   for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix)
     qq.push_back(born->bornOutgoing()[ix]);
   // ensure quark first
   if(qq[0]->id()<0) swap(qq[0],qq[1]);
   // centre of mass energy
   d_Q_ = (qq[0]->momentum() + qq[1]->momentum()).m();
   // quark mass
   d_m_ = 0.5*(qq[0]->momentum().m()+qq[1]->momentum().m());
   // set the other parameters
   setRho(sqr(d_m_/d_Q_));
   setKtildeSymm();
   // otherwise can do it
   initial=1.;
   final  =1.;
 }
 
 bool SMWDecayer::softMatrixElementVeto(PPtr parent,
 				       PPtr progenitor,
 				       const bool & ,
 				       const Energy & highestpT,
 				       const vector<tcPDPtr> & ids,
 				       const double & d_z,
 				       const Energy & d_qt,
 				       const Energy & ) {
   // check we should be applying the veto
   if(parent->id()!=progenitor->id()||
      ids[0]!=ids[1]||
      ids[2]->id()!=ParticleID::g) return false;
   // calculate pt
   Energy2 d_m2 = parent->momentum().m2();
   Energy2 pPerp2 = sqr(d_z*d_qt) - d_m2;
   if(pPerp2<ZERO) return true;
   Energy pPerp = (1.-d_z)*sqrt(pPerp2);
   // if not hardest so far don't apply veto
   if(pPerp<highestpT) return false;
   // calculate the weight
   double weight = 0.;
   if(parent->id()>0) weight = qWeightX(d_qt, d_z);
   else weight = qbarWeightX(d_qt, d_z);
   // compute veto from weight and return
   return !UseRandom::rndbool(weight);
 }
 
 
 void SMWDecayer::setRho(double r) 
 { 
   d_rho_ = r;
   d_v_ = sqrt(1.-4.*d_rho_);
 }
 
 void SMWDecayer::setKtildeSymm() { 
   d_kt1_ = (1. + sqrt(1. - 4.*d_rho_))/2.;
   setKtilde2();
 }
 
 void SMWDecayer::setKtilde2() { 
    double num = d_rho_ * d_kt1_ + 0.25 * d_v_ *(1.+d_v_)*(1.+d_v_);
    double den = d_kt1_ - d_rho_;
    d_kt2_ = num/den;
 }
 
 double SMWDecayer::getZfromX(double x1, double x2) {
   double uval = u(x2);
   double num = x1 - (2. - x2)*uval;
   double den = sqrt(x2*x2 - 4.*d_rho_);
   return uval + num/den;
 }
 
 double SMWDecayer::getKfromX(double x1, double x2) {
    double zval = getZfromX(x1, x2);
    return (1.-x2)/(zval*(1.-zval));
 }
 
 double SMWDecayer::MEV(double x1, double x2) {
   // Vector part
   double num = (x1+2.*d_rho_)*(x1+2.*d_rho_) + (x2+2.*d_rho_)*(x2+2.*d_rho_) 
     - 8.*d_rho_*(1.+2.*d_rho_);
   double den = (1.+2.*d_rho_)*(1.-x1)*(1.-x2);
   return (num/den - 2.*d_rho_/((1.-x1)*(1.-x1)) 
 	  - 2*d_rho_/((1.-x2)*(1.-x2)))/d_v_;
 }
 
 double SMWDecayer::MEA(double x1, double x2) {
   // Axial part
   double num = (x1+2.*d_rho_)*(x1+2.*d_rho_) + (x2+2.*d_rho_)*(x2+2.*d_rho_) 
     + 2.*d_rho_*((5.-x1-x2)*(5.-x1-x2) - 19.0 + 4*d_rho_);
   double den = d_v_*d_v_*(1.-x1)*(1.-x2);
   return (num/den - 2.*d_rho_/((1.-x1)*(1.-x1)) 
 	  - 2*d_rho_/((1.-x2)*(1.-x2)))/d_v_;
 }
 
 double SMWDecayer::u(double x2) {
   return 0.5*(1. + d_rho_/(1.-x2+d_rho_));
 }
 
 void SMWDecayer::
 getXXbar(double kti, double z, double &x, double &xbar) {
   double w = sqr(d_v_) + kti*(-1. + z)*z*(2. + kti*(-1. + z)*z);
   if (w < 0) {
     x = -1.; 
     xbar = -1;
   } else {
     x = (1. + sqr(d_v_)*(-1. + z) + sqr(kti*(-1. + z))*z*z*z 
 	 + z*sqrt(w)
 	 - kti*(-1. + z)*z*(2. + z*(-2 + sqrt(w))))/
       (1. - kti*(-1. + z)*z + sqrt(w));
     xbar = 1. + kti*(-1. + z)*z;
   }
 }
 
 double SMWDecayer::qWeight(double x, double xbar) {
   double rval; 
   double xg = 2. - xbar - x;
   // always return one in the soft gluon region
   if(xg < EPS_) return 1.0;
   // check it is in the phase space
   if((1.-x)*(1.-xbar)*(1.-xg) < d_rho_*xg*xg) return 0.0;
   double k1 = getKfromX(x, xbar);
   double k2 = getKfromX(xbar, x);
   // Is it in the quark emission zone?
   if(k1 < d_kt1_) {
     rval = MEV(x, xbar)/PS(x, xbar);
     // is it also in the anti-quark emission zone?
     if(k2 < d_kt2_) rval *= 0.5;
     return rval;
   }
   return 1.0;
 }
 
 double SMWDecayer::qbarWeight(double x, double xbar) {
   double rval; 
   double xg = 2. - xbar - x;
   // always return one in the soft gluon region
   if(xg < EPS_) return 1.0;
   // check it is in the phase space
   if((1.-x)*(1.-xbar)*(1.-xg) < d_rho_*xg*xg) return 0.0;
   double k1 = getKfromX(x, xbar);
   double k2 = getKfromX(xbar, x);
   // Is it in the antiquark emission zone?
   if(k2 < d_kt2_) {
     rval = MEV(x, xbar)/PS(xbar, x);
     // is it also in the quark emission zone?
     if(k1 < d_kt1_) rval *= 0.5;
     return rval;
   }
   return 1.0;
 }
 
 double SMWDecayer::qWeightX(Energy qtilde, double z) {
   double x, xb;
   getXXbar(sqr(qtilde/d_Q_), z, x, xb);
   // if exceptionally out of phase space, leave this emission, as there 
   // is no good interpretation for the soft ME correction. 
   if (x < 0 || xb < 0) return 1.0; 
   return qWeight(x, xb); 
 }
 
 double SMWDecayer::qbarWeightX(Energy qtilde, double z) {
   double x, xb;
   getXXbar(sqr(qtilde/d_Q_), z, xb, x);
   // see above in qWeightX. 
   if (x < 0 || xb < 0) return 1.0; 
   return qbarWeight(x, xb); 
 }
 
 double SMWDecayer::PS(double x, double xbar) {
   double u = 0.5*(1. + d_rho_ / (1.-xbar+d_rho_));
   double z = u + (x - (2.-xbar)*u)/sqrt(xbar*xbar - 4.*d_rho_);
   double brack = (1.+z*z)/(1.-z)- 2.*d_rho_/(1-xbar);
   // interesting: the splitting function without the subtraction
   // term. Actually gives a much worse approximation in the collinear
   // limit.  double brack = (1.+z*z)/(1.-z);
   double den = (1.-xbar)*sqrt(xbar*xbar - 4.*d_rho_);
   return brack/den;
 }
 
 double SMWDecayer::matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
 				      const ParticleVector & decay3, MEOption,
 				      ShowerInteraction inter) {
   // extract partons and LO momentas
   vector<tcPDPtr> partons(1,inpart.dataPtr());
   vector<Lorentz5Momentum> lomom(1,inpart.momentum());
   for(unsigned int ix=0;ix<2;++ix) {
     partons.push_back(decay2[ix]->dataPtr());
     lomom.push_back(decay2[ix]->momentum());
   }
   vector<Lorentz5Momentum> realmom(1,inpart.momentum());
   for(unsigned int ix=0;ix<3;++ix) {
     if(ix==2) partons.push_back(decay3[ix]->dataPtr());
     realmom.push_back(decay3[ix]->momentum());
   }
   if(partons[0]->id()<0) {
     swap(partons[1],partons[2]);
     swap(lomom[1],lomom[2]);
     swap(realmom[1],realmom[2]);
   }
   scale_ = sqr(inpart.mass());
   double     lome = loME(partons,lomom);
   InvEnergy2 reme = realME(partons,realmom,inter);
   double ratio = reme/lome*sqr(inpart.mass())*4.*Constants::pi;
   if(inter==ShowerInteraction::QCD) ratio *= CF_;
   return ratio;
 }
 
 double SMWDecayer::meRatio(tcPDVector partons, 
 			   vector<Lorentz5Momentum> momenta,
 			   unsigned int iemitter, bool subtract) const {
   Lorentz5Momentum q = momenta[1]+momenta[2]+momenta[3];
   Energy2 Q2=q.m2();
   Energy2 lambda = sqrt((Q2-sqr(momenta[1].mass()+momenta[2].mass()))*
 			(Q2-sqr(momenta[1].mass()-momenta[2].mass())));
   InvEnergy2 D[2];
   double lome(0.);
   for(unsigned int iemit=0;iemit<2;++iemit) {
     unsigned int ispect = iemit==0 ? 1 : 0;    
     Energy2 pipj = momenta[3      ] * momenta[1+iemit ];
     Energy2 pipk = momenta[3      ] * momenta[1+ispect];
     Energy2 pjpk = momenta[1+iemit] * momenta[1+ispect];
     double y = pipj/(pipj+pipk+pjpk); 
     double z = pipk/(     pipk+pjpk);
     Energy mij = sqrt(2.*pipj+sqr(momenta[1+iemit].mass()));
     Energy2 lamB = sqrt((Q2-sqr(mij+momenta[1+ispect].mass()))*
 			(Q2-sqr(mij-momenta[1+ispect].mass())));
     Energy2 Qpk = q*momenta[1+ispect];
     Lorentz5Momentum pkt = 
       lambda/lamB*(momenta[1+ispect]-Qpk/Q2*q)
       +0.5/Q2*(Q2+sqr(momenta[1+ispect].mass())-sqr(momenta[1+ispect].mass()))*q;
     Lorentz5Momentum pijt = 
       q-pkt;
     double muj = momenta[1+iemit ].mass()/sqrt(Q2);
     double muk = momenta[1+ispect].mass()/sqrt(Q2);
     double vt = sqrt((1.-sqr(muj+muk))*(1.-sqr(muj-muk)))/(1.-sqr(muj)-sqr(muk));
     double v  = sqrt(sqr(2.*sqr(muk)+(1.-sqr(muj)-sqr(muk))*(1.-y))-4.*sqr(muk))
       /(1.-y)/(1.-sqr(muj)-sqr(muk));
     // dipole term
     D[iemit] = 0.5/pipj*(2./(1.-(1.-z)*(1.-y))
 			 -vt/v*(2.-z+sqr(momenta[1+iemit].mass())/pipj));
     // matrix element
     vector<Lorentz5Momentum> lomom(3);
     lomom[0] = momenta[0];
     if(iemit==0) {
       lomom[1] = pijt;
       lomom[2] = pkt ;
     }
     else {
       lomom[2] = pijt;
       lomom[1] = pkt ;
     }
     if(iemit==0) lome  = loME(partons,lomom);
   }
   InvEnergy2 ratio = realME(partons,momenta,ShowerInteraction::QCD)/lome*abs(D[iemitter])
     /(abs(D[0])+abs(D[1]));
   if(subtract)
     return Q2*(ratio-2.*D[iemitter]);
   else
     return Q2*ratio;
 }
 
 double SMWDecayer::loME(const vector<tcPDPtr> & partons, 
 			const vector<Lorentz5Momentum> & momenta) const {
   // compute the spinors
   vector<VectorWaveFunction>    vin;
   vector<SpinorWaveFunction>    aout;
   vector<SpinorBarWaveFunction> fout;
   VectorWaveFunction    win  (momenta[0],partons[0],incoming);
   SpinorBarWaveFunction qkout(momenta[1],partons[1],outgoing);
   SpinorWaveFunction    qbout(momenta[2],partons[2],outgoing);
   for(unsigned int ix=0;ix<2;++ix){
     qkout.reset(ix);
     fout.push_back(qkout);
     qbout.reset(ix);
     aout.push_back(qbout);
   }
   for(unsigned int ix=0;ix<3;++ix){
     win.reset(ix);
     vin.push_back(win);
   }
   // temporary storage of the different diagrams
   // sum over helicities to get the matrix element
   double total(0.);
   for(unsigned int inhel=0;inhel<3;++inhel) {
     for(unsigned int outhel1=0;outhel1<2;++outhel1) {
       for(unsigned int outhel2=0;outhel2<2;++outhel2) {
 	Complex diag1 = FFWVertex_->evaluate(scale_,aout[outhel2],fout[outhel1],vin[inhel]);
 	total += norm(diag1);
       }
     }
   }
   // return the answer
   return total;
 }
  
 InvEnergy2 SMWDecayer::realME(const vector<tcPDPtr> & partons, 
 			      const vector<Lorentz5Momentum> & momenta,
 			      ShowerInteraction inter) const {
   // compute the spinors
   vector<VectorWaveFunction>     vin;
   vector<SpinorWaveFunction>     aout;
   vector<SpinorBarWaveFunction>  fout;
   vector<VectorWaveFunction>     gout;
   VectorWaveFunction    win  (momenta[0],partons[0],incoming);
   SpinorBarWaveFunction qkout(momenta[1],partons[1],outgoing);
   SpinorWaveFunction    qbout(momenta[2],partons[2],outgoing);
   VectorWaveFunction    gluon(momenta[3],partons[3],outgoing);
   for(unsigned int ix=0;ix<2;++ix){
     qkout.reset(ix);
     fout.push_back(qkout);
     qbout.reset(ix);
     aout.push_back(qbout);
     gluon.reset(2*ix);
     gout.push_back(gluon);
   }
   for(unsigned int ix=0;ix<3;++ix){
     win.reset(ix);
     vin.push_back(win);
   }
   vector<Complex> diag(3,0.);
 
   double total(0.);
 
   AbstractFFVVertexPtr vertex = inter==ShowerInteraction::QCD ? FFGVertex_ : FFPVertex_;
   
   for(unsigned int inhel1=0;inhel1<3;++inhel1) {
     for(unsigned int outhel1=0;outhel1<2;++outhel1) {
       for(unsigned int outhel2=0;outhel2<2;++outhel2) {
 	for(unsigned int outhel3=0;outhel3<2;++outhel3) {
 	  SpinorBarWaveFunction off1 =
 	    vertex->evaluate(scale_,3,partons[1]->CC(),fout[outhel1],gout[outhel3]);
 	  diag[0] = FFWVertex_->evaluate(scale_,aout[outhel2],off1,vin[inhel1]);
 	  
 	  SpinorWaveFunction off2 = 
 	    vertex->evaluate(scale_,3,partons[2]->CC(),aout[outhel2],gout[outhel3]);
 	  diag[1] = FFWVertex_->evaluate(scale_,off2,fout[outhel1],vin[inhel1]);
 
 	  if(inter==ShowerInteraction::QED) {
 	    VectorWaveFunction off3 =
 	      WWWVertex_->evaluate(scale_,3,partons[0],vin[inhel1],gout[outhel3]);
 	    diag[2] = FFWVertex_->evaluate(scale_,aout[outhel2],fout[outhel1],off3);
 	  }
 	  
 	  // sum of diagrams
 	  Complex sum = std::accumulate(diag.begin(),diag.end(),Complex(0.));
 	  // me2
 	  total += norm(sum);
 	}
       }
     }
   }
   // divide out the coupling
   total /= norm(vertex->norm());
   // double g = sqrt(2.)*abs(FFWVertex_->norm());
   // double xg = 2.*momenta[3].t()/momenta[0].mass();
   // double xe,mue2;
   // if(abs(partons[1]->id())==ParticleID::eminus) {
   //   xe = 2.*momenta[1].t()/momenta[0].mass();
   //   mue2 = sqr(momenta[1].mass()/momenta[0].mass());
   // }
   // else {
   //   xe = 2.*momenta[2].t()/momenta[0].mass();
   //   mue2 = sqr(momenta[2].mass()/momenta[0].mass());
   // }
   // double cg = -4. * g * g * (-pow(mue2, 3.) / 2. + (xg * xg / 4. + (xe / 2. + 1.) * xg + 5. / 2. * xe - 2.) * mue2 * mue2
   // 			     + (pow(xg, 3.) / 4. + (xe / 4. - 5. / 4.) * xg * xg + (-7. / 2. * xe + 3.) * xg - 3. * xe * xe
   // 				+ 11. / 2. * xe - 7. / 2.) * mue2 + (xg * xg / 2. + (xe - 2.) * xg + xe * xe - 2. * xe + 2.) * (-1. + xg + xe)) * (xe - mue2 - 1.) *
   //   pow(xg, -2.) * pow(-1. + xg + xe - mue2, -2.);
   
   // cerr << "real " << cg/total << "\n";
   // return the total
   return total*UnitRemoval::InvE2;
 }
 
 double SMWDecayer::calculateRealEmission(double x1, double x2, 
 					 tcPDVector partons,
 					 vector<Lorentz5Momentum> pin,
 					 double phi, double muj,
 					 double muk, int iemit, 
 					 bool subtract) const {
   // calculate x3
   double x3 = 2.-x1-x2;
   double xT = sqrt(max(0.,sqr(x3)-0.25*sqr(sqr(x2)+sqr(x3)-sqr(x1)-4.*sqr(muk)+4.*sqr(muj))
 		       /(sqr(x2)-4.*sqr(muk))));
   // calculate the momenta
   Energy M = mW_;
   Lorentz5Momentum pspect(ZERO,ZERO,-0.5*M*sqrt(max(sqr(x2)-4.*sqr(muk),0.)),
 			  0.5*M*x2,M*muk); 
   Lorentz5Momentum pemit (-0.5*M*xT*cos(phi),-0.5*M*xT*sin(phi),
 			  0.5*M*sqrt(max(sqr(x1)-sqr(xT)-4.*sqr(muj),0.)),
 			  0.5*M*x1,M*muj);
   Lorentz5Momentum pgluon(0.5*M*xT*cos(phi), 0.5*M*xT*sin(phi),
 			  0.5*M*sqrt(max(sqr(x3)-sqr(xT),0.)),0.5*M*x3,ZERO);
   if(abs(pspect.z()+pemit.z()-pgluon.z())/M<1e-6) 
     pgluon.setZ(-pgluon.z());
   else if(abs(pspect.z()-pemit.z()+pgluon.z())/M<1e-6) 
     pemit .setZ(- pemit.z());
   // boost and rotate momenta
   LorentzRotation eventFrame( ( pin[1] +
 				pin[2] ).findBoostToCM() );
   unsigned int ispect = iemit==0 ? 2 : 1;
   Lorentz5Momentum spectator = eventFrame*pin[ispect];
   eventFrame.rotateZ( -spectator.phi()    );
   eventFrame.rotateY( -spectator.theta()  );
   eventFrame.invert();
   vector<Lorentz5Momentum> momenta(3);
   momenta[0]   = pin[0];
   momenta[ispect ] = eventFrame*pspect;
   momenta[iemit+1] = eventFrame*pemit ;
   momenta.push_back(eventFrame*pgluon);
   // calculate the weight
   double realwgt(0.);
   if(1.-x1>1e-5 && 1.-x2>1e-5) 
     realwgt = meRatio(partons,momenta,iemit,subtract);
   return realwgt;
 }